Context Navigation

source: trunk/GSASIIstrMain.py @ 5028

Last change on this file since 5028 was 5025, checked in by vondreele, 4 years ago
Add shell plots to 3D reflection plots. Fix escape character issues in CopyRietveld2Origin Revise LeBail? operation - now works in sequential refinements. Remove "Fit LeBail?" from Calculate menu. OnLeBail? does 10 cycles instead of 3 in the "zero cycle" step. Initial FOSQ set to Scale*Phase fr. instead of 1.0 - improves initial LeBaill? fit. Put newLeBail flag in Controls - as more global. Allow LeBail? refinements in sequential refinements. Remove a wx.Yield() - obsolete Disable deletion of selected parameters in a LeBail? refinement & let the LS handle any singularities that may ensue Derivatives no longer controlled by LeBail? flag.
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 57.8 KB

Line
1	# -- coding: utf-8 --
2	'''
3	GSASIIstrMain: main structure routine
4	---------------------------------------
5
6	'''
7	########### SVN repository information ###################
8	# $Date: 2021-09-14 16:55:59 +0000 (Tue, 14 Sep 2021) $
9	# $Author: vondreele $
10	# $Revision: 5025 $
11	# $URL: trunk/GSASIIstrMain.py $
12	# $Id: GSASIIstrMain.py 5025 2021-09-14 16:55:59Z vondreele $
13	########### SVN repository information ###################
14	from __future__ import division, print_function
15	import platform
16	import sys
17	import os.path as ospath
18	import time
19	import math
20	import copy
21	if '2' in platform.python_version_tuple()[0]:
22	import cPickle
23	else:
24	import pickle as cPickle
25	import numpy as np
26	import numpy.linalg as nl
27	import scipy.optimize as so
28	import GSASIIpath
29	GSASIIpath.SetBinaryPath()
30	GSASIIpath.SetVersionNumber("$Revision: 5025 $")
31	import GSASIIlattice as G2lat
32	import GSASIIspc as G2spc
33	import GSASIImapvars as G2mv
34	import GSASIImath as G2mth
35	import GSASIIstrIO as G2stIO
36	import GSASIIstrMath as G2stMth
37	import GSASIIobj as G2obj
38	import GSASIIfiles as G2fil
39	import GSASIIElem as G2elem
40	import GSASIIscriptable as G2sc
41	import atmdata
42
43	sind = lambda x: np.sin(x*np.pi/180.)
44	cosd = lambda x: np.cos(x*np.pi/180.)
45	tand = lambda x: np.tan(x*np.pi/180.)
46	asind = lambda x: 180.*np.arcsin(x)/np.pi
47	acosd = lambda x: 180.*np.arccos(x)/np.pi
48	atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
49
50	ateln2 = 8.0*math.log(2.0)
51	DEBUG = True
52
53	def ReportProblems(result,Rvals,varyList):
54	'''Create a message based results from the refinement
55	'''
56	#report on SVD 0's and highly correlated variables
57	msg = ''
58	# process singular variables; all vars go to console, first 10 to
59	# dialog window
60	psing = result[2].get('psing',[])
61	if len(psing):
62	if msg: msg += '\n'
63	m = 'Error: {} Parameter(s) dropped:'.format(len(psing))
64	msg += m
65	G2fil.G2Print(m, mode='warn')
66	m = ''
67	for i,val in enumerate(psing):
68	if i == 0:
69	msg += '\n{}'.format(varyList[val])
70	m = ' {}'.format(varyList[val])
71	else:
72	if len(m) > 70:
73	G2fil.G2Print(m, mode='warn')
74	m = ' '
75	else:
76	m += ', '
77	m += '{}'.format(varyList[val])
78	if i == 10:
79	msg += ', {}... see console for full list'.format(varyList[val])
80	elif i > 10:
81	pass
82	else:
83	msg += ', {}'.format(varyList[val])
84	if m: G2fil.G2Print(m, mode='warn')
85	SVD0 = result[2].get('SVD0',0)
86	if SVD0 == 1:
87	msg += 'Warning: Soft (SVD) singularity in the Hessian'
88	elif SVD0 > 0:
89	msg += 'Warning: {} soft (SVD) Hessian singularities'.format(SVD0)
90	SVDsing = result[2].get('SVDsing',[])
91	if len(SVDsing):
92	if msg: msg += '\n'
93	m = 'SVD problem(s) likely from:'
94	msg += m
95	G2fil.G2Print(m, mode='warn')
96	m = ''
97	for i,val in enumerate(SVDsing):
98	if i == 0:
99	msg += '\n{}'.format(varyList[val])
100	m = ' {}'.format(varyList[val])
101	else:
102	if len(m) > 70:
103	G2fil.G2Print(m, mode='warn')
104	m = ' '
105	else:
106	m += ', '
107	m += '{}'.format(varyList[val])
108	if i == 10:
109	msg += ', {}... see console for full list'.format(varyList[val])
110	elif i > 10:
111	pass
112	else:
113	msg += ', {}'.format(varyList[val])
114	if m: G2fil.G2Print(m, mode='warn')
115	#report on highly correlated variables
116	Hcorr = result[2].get('Hcorr',[])
117	if len(Hcorr) > 0:
118	if msg: msg += '\n'
119	m = 'Note highly correlated parameters:'
120	G2fil.G2Print(m, mode='warn')
121	msg += m
122	elif SVD0 > 0:
123	if msg: msg += '\n'
124	m = 'Check covariance matrix for parameter correlation'
125	G2fil.G2Print(m, mode='warn')
126	msg += m
127	for i,(v1,v2,corr) in enumerate(Hcorr):
128	if corr > .95:
129	stars = '**'
130	else:
131	stars = ' '
132	m = ' {} {} and {} (@{:.2f}%)'.format(
133	stars,varyList[v1],varyList[v2],100.*corr)
134	G2fil.G2Print(m, mode='warn')
135	if i == 5:
136	msg += '\n' + m
137	msg += '\n ... check console for more'
138	elif i < 5:
139	msg += '\n' + m
140	if msg:
141	if 'msg' not in Rvals: Rvals['msg'] = ''
142	Rvals['msg'] += msg
143
144	def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
145	calcControls,pawleyLookup,ifSeq,printFile,dlg,refPlotUpdate=None):
146	'''Core optimization routines, shared between SeqRefine and Refine
147
148	:returns: 5-tuple of ifOk (bool), Rvals (dict), result, covMatrix, sig
149	'''
150	#patch (added Oct 2020) convert variable names for parm limits to G2VarObj
151	G2sc.patchControls(Controls)
152	# end patch
153	# print 'current',varyList
154	# for item in parmDict: print item,parmDict[item] ######### show dict just before refinement
155	ifPrint = True
156	if ifSeq:
157	ifPrint = False
158	Rvals = {}
159	chisq0 = None
160	Lastshft = None
161	while True:
162	G2mv.Map2Dict(parmDict,varyList)
163	begin = time.time()
164	values = np.array(G2stMth.Dict2Values(parmDict, varyList))
165	if np.any(np.isnan(values)):
166	raise G2obj.G2Exception('ERROR - nan found in LS parameters - use Calculate/View LS parms to locate')
167	# test code to compute GOF and save for external repeat
168	#args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)
169	#print '* before fit chi2',np.sum(G2stMth.errRefine(values,args)*2)
170	#fl = open('beforeFit.cpickle','wb')
171	#cPickle.dump(values,fl,1)
172	#cPickle.dump(args[:-1],fl,1)
173	#fl.close()
174	Ftol = Controls['min dM/M']
175	Xtol = Controls['SVDtol']
176	Factor = Controls['shift factor']
177	if 'Jacobian' in Controls['deriv type']:
178	result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
179	ftol=Ftol,col_deriv=True,factor=Factor,
180	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
181	ncyc = int(result[2]['nfev']/2)
182	if refPlotUpdate is not None: refPlotUpdate(Histograms) # update plot after completion
183	elif 'analytic Hessian' in Controls['deriv type']:
184	Lamda = Controls.get('Marquardt',-3)
185	maxCyc = Controls['max cyc']
186	result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
187	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
188	refPlotUpdate=refPlotUpdate)
189	ncyc = result[2]['num cyc']+1
190	Rvals['lamMax'] = result[2]['lamMax']
191	if 'Ouch#4' in result[2]:
192	Rvals['Aborted'] = True
193	if 'msg' in result[2]:
194	Rvals['msg'] = result[2]['msg']
195	Controls['Marquardt'] = -3 #reset to default
196	if 'chisq0' in result[2] and chisq0 is None:
197	chisq0 = result[2]['chisq0']
198	elif 'Hessian SVD' in Controls['deriv type']:
199	maxCyc = Controls['max cyc']
200	result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
201	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
202	refPlotUpdate=refPlotUpdate)
203	if result[1] is None:
204	IfOK = False
205	covMatrix = []
206	sig = len(varyList)*[None,]
207	break
208	ncyc = result[2]['num cyc']+1
209	if 'chisq0' in result[2] and chisq0 is None:
210	chisq0 = result[2]['chisq0']
211	else: #'numeric'
212	result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
213	args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
214	ncyc = 1
215	if len(varyList):
216	ncyc = int(result[2]['nfev']/len(varyList))
217	if refPlotUpdate is not None: refPlotUpdate(Histograms) # update plot
218	# table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
219	# for item in table: print item,table[item] #useful debug - are things shifting?
220	runtime = time.time()-begin
221	Rvals['SVD0'] = result[2].get('SVD0',0)
222	Rvals['converged'] = result[2].get('Converged')
223	Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
224	Rvals['chisq'] = np.sum(result[2]['fvec']**2)
225	G2stMth.Values2Dict(parmDict, varyList, result[0])
226	G2mv.Dict2Map(parmDict,varyList)
227	Rvals['Nobs'] = Histograms['Nobs']
228	Rvals['Nvars'] = len(varyList)
229	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to %
230	Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
231	printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'% \
232	(result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
233	if ncyc:
234	printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
235	printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
236	sig = len(varyList)*[None,]
237	if 'None' in str(type(result[1])) and ifSeq: #this bails out of a sequential refinement on singular matrix
238	IfOK = False
239	covMatrix = []
240	G2fil.G2Print ('Warning: ** Refinement failed - singular matrix **')
241	if 'Hessian' in Controls['deriv type']:
242	num = len(varyList)-1
243	# BHT -- I am not sure if this works correctly:
244	for i,val in enumerate(np.flipud(result[2]['psing'])):
245	if val:
246	G2fil.G2Print('Bad parameter: '+varyList[num-i],mode='warn')
247	else:
248	Ipvt = result[2]['ipvt']
249	for i,ipvt in enumerate(Ipvt):
250	if not np.sum(result[2]['fjac'],axis=1)[i]:
251	G2fil.G2Print('Bad parameter: '+varyList[ipvt-1],mode='warn')
252	break
253	IfOK = True
254	if not len(varyList):
255	covMatrix = []
256	break
257	try:
258	covMatrix = result[1]Rvals['GOF']*2
259	sig = np.sqrt(np.diag(covMatrix))
260	Lastshft = result[0]-values #NOT last shift since values is starting set before current refinement
261	Rvals['Max shft/sig'] = np.max(np.nan_to_num(Lastshft/sig))
262	if np.any(np.isnan(sig)) or not sig.shape:
263	G2fil.G2Print ('* Least squares aborted - some invalid esds possible *',mode='error')
264	# report on refinement issues. Result in Rvals['msg']
265	ReportProblems(result,Rvals,varyList)
266	break #refinement succeeded - finish up!
267	except TypeError:
268	# if we get here, no result[1] (covar matrix) was returned or other calc error: refinement failed
269	IfOK = False
270	if 'Hessian' in Controls['deriv type']:
271	SVD0 = result[2].get('SVD0')
272	if SVD0 == -1:
273	G2fil.G2Print ('** Refinement failed - singular matrix **',mode='error')
274	elif SVD0 == -2:
275	G2fil.G2Print ('** Refinement failed - other problem **',mode='error')
276	elif SVD0 > 0:
277	G2fil.G2Print ('** Refinement failed with {} SVD singularities **'.format(SVD0),mode='error')
278	else:
279	G2fil.G2Print ('** Refinement failed **',mode='error')
280	if result[1] is None:
281	IfOK = False
282	covMatrix = []
283	sig = len(varyList)*[None,]
284	# report on highly correlated variables
285	ReportProblems(result,Rvals,varyList)
286	# process singular variables
287	if dlg: break # refining interactively
288	else:
289	G2fil.G2Print ('** Refinement failed - singular matrix **',mode='error')
290	Ipvt = result[2]['ipvt']
291	for i,ipvt in enumerate(Ipvt):
292	if not np.sum(result[2]['fjac'],axis=1)[i]:
293	G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
294	del(varyList[ipvt-1])
295	break
296	if IfOK:
297	if CheckLeBail(Phases): # only needed for LeBail extraction
298	G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],
299	parmDict,[],calcControls,pawleyLookup,dlg)
300	G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
301	if chisq0 is not None:
302	Rvals['GOF0'] = np.sqrt(chisq0/(Histograms['Nobs']-len(varyList)))
303	return IfOK,Rvals,result,covMatrix,sig,Lastshft
304
305	def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None):
306	'''Global refinement -- refines to minimize against all histograms.
307	This can be called in one of three ways, from :meth:`GSASIIdataGUI.GSASII.OnRefine` in an
308	interactive refinement, where dlg will be a wx.ProgressDialog, or non-interactively from
309	:meth:`GSASIIscriptable.G2Project.refine` or from :func:`main`, where dlg will be None.
310	'''
311	import GSASIImpsubs as G2mp
312	G2mp.InitMP()
313	import pytexture as ptx
314	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
315
316	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
317	G2stIO.ShowBanner(printFile)
318	varyList = []
319	parmDict = {}
320	G2mv.InitVars()
321	Controls = G2stIO.GetControls(GPXfile)
322	G2stIO.ShowControls(Controls,printFile)
323	calcControls = {}
324	calcControls.update(Controls)
325	constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
326	restraintDict = G2stIO.GetRestraints(GPXfile)
327	Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
328	if not Phases:
329	G2fil.G2Print (' * ERROR - you have no phases to refine! *')
330	G2fil.G2Print (' * Refine aborted *')
331	return False,{'msg':'No phases'}
332	if not Histograms:
333	G2fil.G2Print (' * ERROR - you have no data to refine with! *')
334	G2fil.G2Print (' * Refine aborted *')
335	return False,{'msg':'No data'}
336	rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
337	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
338	rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
339	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
340	maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
341	calcControls['atomIndx'] = atomIndx
342	calcControls['Natoms'] = Natoms
343	calcControls['FFtables'] = FFtables
344	calcControls['BLtables'] = BLtables
345	calcControls['MFtables'] = MFtables
346	calcControls['maxSSwave'] = maxSSwave
347	hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,Controls=calcControls,pFile=printFile)
348	TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
349	constrDict += TwConstr
350	fixedList += TwFixed
351	calcControls.update(controlDict)
352	histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
353	calcControls.update(controlDict)
354	varyList = rbVary+phaseVary+hapVary+histVary
355	parmDict.update(rbDict)
356	parmDict.update(phaseDict)
357	parmDict.update(hapDict)
358	parmDict.update(histDict)
359	G2stIO.GetFprime(calcControls,Histograms)
360	# do constraint processing
361	varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
362	msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
363	if msg:
364	return False,{'msg':'Unable to interpret multiplier(s): '+msg}
365	try:
366	G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict)
367	#print(G2mv.VarRemapShow(varyList))
368	#print('DependentVars',G2mv.GetDependentVars())
369	#print('IndependentVars',G2mv.GetIndependentVars())
370	except G2mv.ConstraintException:
371	G2fil.G2Print (' * ERROR - your constraints are internally inconsistent *')
372	#errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
373	#print 'Errors',errmsg
374	#if warnmsg: print 'Warnings',warnmsg
375	return False,{'msg':' Constraint error'}
376	# print G2mv.VarRemapShow(varyList)
377
378	# remove frozen vars from refinement
379	if 'parmFrozen' not in Controls:
380	Controls['parmFrozen'] = {}
381	if 'FrozenList' not in Controls['parmFrozen']:
382	Controls['parmFrozen']['FrozenList'] = []
383	parmFrozenList = Controls['parmFrozen']['FrozenList']
384	frozenList = [i for i in varyList if i in parmFrozenList]
385	if len(frozenList) != 0:
386	varyList = [i for i in varyList if i not in parmFrozenList]
387	G2fil.G2Print(
388	'Frozen refined variables (due to exceeding limits)\n\t:{}'
389	.format(frozenList))
390
391	ifSeq = False
392	printFile.write('\n Refinement results:\n')
393	printFile.write(135*'-'+'\n')
394	Rvals = {}
395	try:
396	covData = {}
397	IfOK,Rvals,result,covMatrix,sig,Lastshft = RefineCore(Controls,Histograms,Phases,restraintDict,
398	rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
399	refPlotUpdate=refPlotUpdate)
400	if IfOK:
401	sigDict = dict(zip(varyList,sig))
402	newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
403	newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
404	covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
405	'varyListStart':varyListStart,'Lastshft':Lastshft,
406	'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
407	'newCellDict':newCellDict,'freshCOV':True}
408	# add the uncertainties into the esd dictionary (sigDict)
409	sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
410	# check for variables outside their allowed range, reset and freeze them
411	frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
412	G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
413	G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
414	G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
415	G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
416	G2stIO.PrintISOmodes(printFile,Phases,parmDict,sigDict)
417	G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,pFile=printFile)
418	G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls,pFile=printFile)
419	if len(frozen):
420	if 'msg' in Rvals:
421	Rvals['msg'] += '\n'
422	else:
423	Rvals['msg'] = ''
424	msg = ('Warning: {} variable(s) refined outside limits and were frozen ({} total frozen)'
425	.format(len(frozen),len(parmFrozenList))
426	)
427	G2fil.G2Print(msg)
428	Rvals['msg'] += msg
429	elif len(parmFrozenList):
430	if 'msg' in Rvals:
431	Rvals['msg'] += '\n'
432	else:
433	Rvals['msg'] = ''
434	msg = ('Note: a total of {} variable(s) are frozen due to refining outside limits'
435	.format(len(parmFrozenList))
436	)
437	G2fil.G2Print('Note: ',msg)
438	Rvals['msg'] += msg
439	G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,makeBack)
440	printFile.close()
441	G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
442	G2fil.G2Print (' *** Refinement successful ***')
443	else:
444	G2fil.G2Print ('****ERROR - Refinement failed',mode='error')
445	if 'msg' in Rvals:
446	G2fil.G2Print ('Note refinement problem:',mode='warn')
447	G2fil.G2Print (Rvals['msg'],mode='warn')
448	raise G2obj.G2Exception('** ERROR: Refinement failed **')
449	except G2obj.G2RefineCancel as Msg:
450	printFile.close()
451	G2fil.G2Print (' *** Refinement stopped ***')
452	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
453	if 'msg' in Rvals:
454	Rvals['msg'] += '\n'
455	Rvals['msg'] += Msg.msg
456	if not dlg:
457	G2fil.G2Print ('Note refinement problem:',mode='warn')
458	G2fil.G2Print (Rvals['msg'],mode='warn')
459	else:
460	Rvals['msg'] = Msg.msg
461	return False,Rvals
462	# except G2obj.G2Exception as Msg: # cell metric error, others?
463	except Exception as Msg: # cell metric error, others?
464	if GSASIIpath.GetConfigValue('debug'):
465	import traceback
466	print(traceback.format_exc())
467	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
468	printFile.close()
469	G2fil.G2Print (' *** Refinement error ***')
470	if 'msg' in Rvals:
471	Rvals['msg'] += '\n\n'
472	Rvals['msg'] += Msg.msg
473	if not dlg:
474	G2fil.G2Print ('Note refinement problem:',mode='warn')
475	G2fil.G2Print (Rvals['msg'],mode='warn')
476	else:
477	Rvals['msg'] = Msg.msg
478	return False,Rvals
479
480	#for testing purposes, create a file for testderiv
481	if GSASIIpath.GetConfigValue('debug'): # and IfOK:
482	#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
483	fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
484	cPickle.dump(result[0],fl,1)
485	cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
486	cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
487	cPickle.dump(parmDict,fl,1)
488	cPickle.dump(varyList,fl,1)
489	cPickle.dump(calcControls,fl,1)
490	cPickle.dump(pawleyLookup,fl,1)
491	fl.close()
492	if dlg:
493	return True,Rvals
494	elif 'msg' in Rvals:
495	G2fil.G2Print ('Reported from refinement:',mode='warn')
496	G2fil.G2Print (Rvals['msg'],mode='warn')
497
498	def CheckLeBail(Phases):
499	'''Check if there is a LeBail extraction in any histogram
500
501	:returns: True if there is at least one LeBail flag turned on, False otherwise
502	'''
503	for key in Phases:
504	phase = Phases[key]
505	for h in phase['Histograms']:
506	#phase['Histograms'][h]
507	if not phase['Histograms'][h]['Use']: continue
508	try:
509	if phase['Histograms'][h]['LeBail']:
510	return True
511	except KeyError: #HKLF & old gpx files
512	pass
513	return False
514
515	def DoLeBail(GPXfile,dlg=None,cycles=10,refPlotUpdate=None):
516	'''Fit LeBail intensities without changes to any other refined parameters.
517	This is a stripped-down version of :func:`Refine` that does not perform
518	any refinement cycles
519	'''
520	import GSASIImpsubs as G2mp
521	G2mp.InitMP()
522	import pytexture as ptx
523	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
524
525	#varyList = []
526	parmDict = {}
527	Controls = G2stIO.GetControls(GPXfile)
528	calcControls = {}
529	calcControls.update(Controls)
530	constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
531	restraintDict = {}
532	Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
533	if not Phases:
534	G2fil.G2Print (' * ERROR - you have no phases to refine! *')
535	return False,{'msg':'No phases'}
536	if not Histograms:
537	G2fil.G2Print (' * ERROR - you have no data to refine with! *')
538	return False,{'msg':'No data'}
539	if not CheckLeBail(Phases):
540	msg = 'Warning: There are no histograms with LeBail extraction enabled'
541	G2fil.G2Print ('* '+msg+' *')
542	return False,{'msg': msg}
543	rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
544	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
545	rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,Print=False)
546	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
547	maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,Print=False)
548	calcControls['atomIndx'] = atomIndx
549	calcControls['Natoms'] = Natoms
550	calcControls['FFtables'] = FFtables
551	calcControls['BLtables'] = BLtables
552	calcControls['MFtables'] = MFtables
553	calcControls['maxSSwave'] = maxSSwave
554	hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,Controls=calcControls,Print=False)
555	calcControls.update(controlDict)
556	histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,Print=False)
557	calcControls.update(controlDict)
558	parmDict.update(rbDict)
559	parmDict.update(phaseDict)
560	parmDict.update(hapDict)
561	parmDict.update(histDict)
562	G2stIO.GetFprime(calcControls,Histograms)
563	try:
564	for i in range(cycles):
565	M = G2stMth.errRefine([],[Histograms,Phases,restraintDict,rigidbodyDict],parmDict,[],calcControls,pawleyLookup,dlg)
566	G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
567	if refPlotUpdate is not None: refPlotUpdate(Histograms,i)
568	Rvals = {}
569	Rvals['chisq'] = np.sum(M**2)
570	Rvals['Nobs'] = Histograms['Nobs']
571	Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100. #to %
572	Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs'])) # no variables
573
574	covData = {'variables':0,'varyList':[],'sig':[],'Rvals':Rvals,
575	'varyListStart':[],
576	'covMatrix':np.zeros([0,0]),'title':GPXfile,
577	#'newAtomDict':newAtomDict,'newCellDict':newCellDict,
578	'freshCOV':True}
579
580	G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,[],True)
581	G2fil.G2Print (' *** LeBail fit completed ***')
582	return True,Rvals
583	except Exception as Msg:
584	G2fil.G2Print (' *** LeBail fit error ***')
585	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
586	if GSASIIpath.GetConfigValue('debug'):
587	import traceback
588	print(traceback.format_exc())
589	return False,{'msg':Msg.msg}
590
591	def phaseCheck(phaseVary,Phases,histogram):
592	'''
593	Removes unused parameters from phase varylist if phase not in histogram
594	for seq refinement removes vars in "Fix FXU" and "FixedSeqVars" here
595	'''
596	NewVary = []
597	for phase in Phases:
598	if histogram not in Phases[phase]['Histograms']: continue
599	if Phases[phase]['Histograms'][histogram]['Use']:
600	pId = Phases[phase]['pId']
601	newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
602	FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
603	if 'F' in FixVals:
604	newVary = [item for item in newVary if not 'Afrac' in item]
605	if 'X' in FixVals:
606	newVary = [item for item in newVary if not 'dA' in item]
607	if 'U' in FixVals:
608	newVary = [item for item in newVary if not 'AU' in item]
609	if 'M' in FixVals:
610	newVary = [item for item in newVary if not 'AM' in item]
611	removeVars = Phases[phase]['Histograms'][histogram].get('FixedSeqVars',[])
612	newVary = [item for item in newVary if item not in removeVars]
613	NewVary += newVary
614	return NewVary
615
616	def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
617	'''Perform a sequential refinement -- cycles through all selected histgrams,
618	one at a time
619	'''
620	import GSASIImpsubs as G2mp
621	G2mp.InitMP()
622	import pytexture as ptx
623	ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics
624	msgs = {}
625	printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
626	G2fil.G2Print ('Starting Sequential Refinement')
627	G2stIO.ShowBanner(printFile)
628	Controls = G2stIO.GetControls(GPXfile)
629	preFrozenCount = 0
630	for h in Controls['parmFrozen']:
631	if h == 'FrozenList': continue
632	preFrozenCount += len(Controls['parmFrozen'][h])
633	G2stIO.ShowControls(Controls,printFile,SeqRef=True,
634	preFrozenCount=preFrozenCount)
635	restraintDict = G2stIO.GetRestraints(GPXfile)
636	Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
637	if not Phases:
638	G2fil.G2Print (' * ERROR - you have no phases to refine! *')
639	G2fil.G2Print (' * Refine aborted *')
640	return False,'No phases'
641	if not Histograms:
642	G2fil.G2Print (' * ERROR - you have no data to refine with! *')
643	G2fil.G2Print (' * Refine aborted *')
644	return False,'No data'
645	rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
646	rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
647	rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
648	G2mv.InitVars()
649	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
650	maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,
651	Print=False,pFile=printFile,seqRef=True)
652	for item in phaseVary:
653	if '::A0' in item:
654	G2fil.G2Print ('** WARNING - lattice parameters should not be refined in a sequential refinement **')
655	G2fil.G2Print ('** instead use the Dij parameters for each powder histogram **')
656	return False,'Lattice parameter refinement error - see console message'
657	if '::C(' in item:
658	G2fil.G2Print ('** WARNING - phase texture parameters should not be refined in a sequential refinement **')
659	G2fil.G2Print ('** instead use the C(L,N) parameters for each powder histogram **')
660	return False,'Phase texture refinement error - see console message'
661	if 'Seq Data' in Controls:
662	histNames = Controls['Seq Data']
663	else: # patch from before Controls['Seq Data'] was implemented?
664	histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
665	if Controls.get('Reverse Seq'):
666	histNames.reverse()
667	SeqResult = G2stIO.GetSeqResult(GPXfile)
668	# SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
669	Histo = {}
670	NewparmDict = {}
671	G2stIO.SetupSeqSavePhases(GPXfile)
672	msgs['steepestNum'] = 0
673	msgs['maxshift/sigma'] = []
674	for ihst,histogram in enumerate(histNames):
675	if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
676	G2fil.G2Print('\nRefining with '+str(histogram))
677	G2mv.InitVars()
678	(Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
679	FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
680	Phases,restraintDict,rbIds,
681	Print=False,pFile=printFile,seqRef=True)
682	ifPrint = False
683	if dlg:
684	dlg.SetTitle('Residual for histogram '+str(ihst))
685	calcControls = {}
686	calcControls['atomIndx'] = atomIndx
687	calcControls['Natoms'] = Natoms
688	calcControls['FFtables'] = FFtables
689	calcControls['BLtables'] = BLtables
690	calcControls['MFtables'] = MFtables
691	calcControls['maxSSwave'] = maxSSwave
692	if histogram not in Histograms:
693	G2fil.G2Print("Error: not found!")
694	raise G2obj.G2Exception("refining with invalid histogram {}".
695	format(histogram))
696	hId = Histograms[histogram]['hId']
697	redphaseVary = phaseCheck(phaseVary,Phases,histogram)
698	Histo = {histogram:Histograms[histogram],}
699	hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Controls=calcControls,Print=False)
700	calcControls.update(controlDict)
701	histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
702	calcControls.update(controlDict)
703	varyList = rbVary+redphaseVary+hapVary+histVary
704	# if not ihst:
705	# save the initial vary list, but without histogram numbers on parameters
706	saveVaryList = varyList[:]
707	for i,item in enumerate(saveVaryList):
708	items = item.split(':')
709	if items[1]:
710	items[1] = ''
711	item = ':'.join(items)
712	saveVaryList[i] = item
713	if not ihst:
714	SeqResult['varyList'] = saveVaryList
715	else:
716	SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
717	parmDict = {}
718	parmDict.update(rbDict)
719	parmDict.update(phaseDict)
720	parmDict.update(hapDict)
721	parmDict.update(histDict)
722	if Controls['Copy2Next']: # update with parms from last histogram
723	#parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
724	for parm in NewparmDict:
725	if parm in parmDict:
726	parmDict[parm] = NewparmDict[parm]
727	elif histogram in SeqResult: # update phase from last seq ref
728	NewparmDict = SeqResult[histogram].get('parmDict',{})
729	for parm in NewparmDict:
730	if '::' in parm and parm in parmDict:
731	parmDict[parm] = NewparmDict[parm]
732
733	G2stIO.GetFprime(calcControls,Histo)
734	# do constraint processing
735	#reload(G2mv) # debug
736	constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
737	varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
738	msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
739	if msg:
740	return False,'Unable to interpret multiplier(s): '+msg
741	try:
742	groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
743	# if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
744	# G2mv.VarRemapShow(varyList,True))
745	# print('DependentVars',G2mv.GetDependentVars())
746	# print('IndependentVars',G2mv.GetIndependentVars())
747	constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
748	except G2mv.ConstraintException:
749	G2fil.G2Print (' * ERROR - your constraints are internally inconsistent *')
750	#errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
751	#print 'Errors',errmsg
752	#if warnmsg: print 'Warnings',warnmsg
753	return False,' Constraint error'
754	#print G2mv.VarRemapShow(varyList)
755	if not ihst:
756	# first histogram to refine against
757	firstVaryList = []
758	for item in varyList:
759	items = item.split(':')
760	if items[1]:
761	items[1] = ''
762	item = ':'.join(items)
763	firstVaryList.append(item)
764	newVaryList = firstVaryList
765	else:
766	newVaryList = []
767	for item in varyList:
768	items = item.split(':')
769	if items[1]:
770	items[1] = ''
771	item = ':'.join(items)
772	newVaryList.append(item)
773	if newVaryList != firstVaryList and Controls['Copy2Next']:
774	# variable lists are expected to match between sequential refinements when Copy2Next is on
775	#print '**** ERROR - variable list for this histogram does not match previous'
776	#print ' Copy of variables is not possible'
777	#print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
778	combined = list(set(firstVaryList+newVaryList))
779	c = [var for var in combined if var not in newVaryList]
780	p = [var for var in combined if var not in firstVaryList]
781	G2fil.G2Print('* Variables change *')
782	for typ,vars in [('Removed',c),('Added',p)]:
783	line = ' '+typ+': '
784	if vars:
785	for var in vars:
786	if len(line) > 70:
787	G2fil.G2Print(line)
788	line = ' '
789	line += var + ', '
790	else:
791	line += 'none, '
792	G2fil.G2Print(line[:-2])
793	firstVaryList = newVaryList
794
795	ifSeq = True
796	printFile.write('\n Refinement results for histogram id {}: {}\n'
797	.format(hId,histogram))
798	printFile.write(135*'-'+'\n')
799	# remove frozen vars
800	if 'parmFrozen' not in Controls:
801	Controls['parmFrozen'] = {}
802	if histogram not in Controls['parmFrozen']:
803	Controls['parmFrozen'][histogram] = []
804	parmFrozenList = Controls['parmFrozen'][histogram]
805	frozenList = [i for i in varyList if i in parmFrozenList]
806	if len(frozenList) != 0:
807	varyList = [i for i in varyList if i not in parmFrozenList]
808	s = ''
809	for a in frozenList:
810	if s:
811	s+= ', '
812	s += a
813	printFile.write(
814	' The following refined variables have previously been frozen due to exceeding limits:\n\t{}\n'
815	.format(s))
816	try:
817	IfOK,Rvals,result,covMatrix,sig,Lastshft = RefineCore(Controls,Histo,Phases,restraintDict,
818	rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
819	refPlotUpdate=refPlotUpdate)
820	try:
821	shft = '%.4f'% Rvals['Max shft/sig']
822	except:
823	shft = '?'
824	G2fil.G2Print (' wR = %7.2f%%, chi2 = %12.6g, reduced chi2 = %6.2f, last delta chi = %.4f, last shft/sig = %s'%(
825	Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'],shft))
826	if Rvals.get('lamMax',0) >= 10.:
827	msgs['steepestNum'] += 1
828	if Rvals.get('Max shft/sig'):
829	msgs['maxshift/sigma'].append(Rvals['Max shft/sig'])
830	# add the uncertainties into the esd dictionary (sigDict)
831	if not IfOK:
832	G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
833	break
834	sigDict = dict(zip(varyList,sig))
835	# the uncertainties for dependent constrained parms into the esd dict
836	sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
837	# check for variables outside their allowed range, reset and freeze them
838	frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
839	msg = None
840	if len(frozen) > 0:
841	msg = ('Hist {}: {} variables were outside limits and were frozen (now {} frozen total)'
842	.format(ihst,len(frozen),len(parmFrozenList)))
843	G2fil.G2Print(msg)
844	msg = (' {} variables were outside limits and were frozen (now {} frozen total)'
845	.format(len(frozen),len(parmFrozenList)))
846	for p in frozen:
847	if p not in varyList:
848	print('Frozen Warning: {} not in varyList. This should not happen!'.format(p))
849	continue
850	i = varyList.index(p)
851	result[0][i] = parmDict[p]
852	sig[i] = -0.1
853	# a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
854	depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
855	newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
856	newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
857	histRefData = {
858	'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
859	'varyListStart':varyListStart,'Lastshft':Lastshft,
860	'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
861	'newCellDict':newCellDict,'depParmDict':depParmDict,
862	'constraintInfo':constraintInfo,
863	'parmDict':parmDict,
864	}
865	SeqResult[histogram] = histRefData
866	G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
867	# G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
868	G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
869	G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile,seq=True)
870	G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,Controls['parmFrozen'])
871	if msg:
872	printFile.write(msg+'\n')
873	NewparmDict = {}
874	# make dict of varied parameters in current histogram, renamed to
875	# next histogram, for use in next refinement.
876	if Controls['Copy2Next'] and ihst < len(histNames)-1:
877	hId = Histo[histogram]['hId'] # current histogram
878	nexthId = Histograms[histNames[ihst+1]]['hId']
879	for parm in set(list(varyList)+list(varyListStart)):
880	items = parm.split(':')
881	if len(items) < 3:
882	continue
883	if str(hId) in items[1]:
884	items[1] = str(nexthId)
885	newparm = ':'.join(items)
886	NewparmDict[newparm] = parmDict[parm]
887	else:
888	if items[2].startswith('dA'): parm = parm.replace(':dA',':A')
889	NewparmDict[parm] = parmDict[parm]
890
891	except G2obj.G2RefineCancel as Msg:
892	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
893	printFile.close()
894	G2fil.G2Print (' *** Refinement stopped ***')
895	return False,Msg.msg
896	except G2obj.G2Exception as Msg: # cell metric error, others?
897	if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
898	printFile.close()
899	G2fil.G2Print (' *** Refinement error ***')
900	return False,Msg.msg
901	if GSASIIpath.GetConfigValue('Show_timing'):
902	t2 = time.time()
903	G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
904	t1 = t2
905	SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
906	try:
907	G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
908	except Exception as msg:
909	print('Error reading Sequential results\n',str(msg))
910	if GSASIIpath.GetConfigValue('debug'):
911	import traceback
912	print(traceback.format_exc())
913	postFrozenCount = 0
914	for h in Controls['parmFrozen']:
915	if h == 'FrozenList': continue
916	postFrozenCount += len(Controls['parmFrozen'][h])
917	if postFrozenCount:
918	msgs['Frozen'] = 'Ending refinement with {} Frozen variables ({} added now)\n'.format(postFrozenCount,postFrozenCount-preFrozenCount)
919	printFile.write('\n'+msgs['Frozen'])
920	printFile.close()
921	G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
922	G2fil.G2Print (' *** Sequential refinement successful ***')
923	return True,msgs
924
925	def dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList):
926	'''Find variables in the parameters dict that are outside the ranges
927	(in parmMinDict and parmMaxDict) and set them to the limits values.
928	Add any such variables into the list of frozen variable
929	(parmFrozenList). Returns a list of newly frozen variables, if any.
930	'''
931	parmMinDict = Controls.get('parmMinDict',{})
932	parmMaxDict = Controls.get('parmMaxDict',{})
933	freeze = []
934	if parmMinDict or parmMaxDict:
935	for name in varyList:
936	if name not in parmDict: continue
937	n,val = G2obj.prmLookup(name,parmMinDict)
938	if n is not None:
939	if parmDict[name] < parmMinDict[n]:
940	parmDict[name] = parmMinDict[n]
941	sigDict[name] = 0.0
942	freeze.append(name)
943	continue
944	n,val = G2obj.prmLookup(name,parmMaxDict)
945	if n is not None:
946	if parmDict[name] > parmMaxDict[n]:
947	parmDict[name] = parmMaxDict[n]
948	sigDict[name] = 0.0
949	freeze.append(name)
950	continue
951	for v in freeze:
952	if v not in parmFrozenList:
953	parmFrozenList.append(v)
954	return freeze
955
956	def RetDistAngle(DisAglCtls,DisAglData,dlg=None):
957	'''Compute and return distances and angles
958
959	:param dict DisAglCtls: contains distance/angle radii usually defined using
960	:func:`GSASIIctrlGUI.DisAglDialog`
961	:param dict DisAglData: contains phase data:
962	Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
963	for distance/angle origins and atoms to be used as targets.
964	Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
965
966	:returns: AtomLabels,DistArray,AngArray where:
967
968	AtomLabels is a dict of atom labels, keys are the atom number
969
970	DistArray is a dict keyed by the origin atom number where the value is a list
971	of distance entries. The value for each distance is a list containing:
972
973	0) the target atom number (int);
974	1) the unit cell offsets added to x,y & z (tuple of int values)
975	2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
976	3) an interatomic distance in A (float)
977	4) an uncertainty on the distance in A or 0.0 (float)
978
979	AngArray is a dict keyed by the origin (central) atom number where
980	the value is a list of
981	angle entries. The value for each angle entry consists of three values:
982
983	0) a distance item reference for one neighbor (int)
984	1) a distance item reference for a second neighbor (int)
985	2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
986
987	The AngArray distance reference items refer directly to the index of the items in the
988	DistArray item for the list of distances for the central atom.
989	'''
990	import numpy.ma as ma
991
992	SGData = DisAglData['SGData']
993	Cell = DisAglData['Cell']
994	Amat,Bmat = G2lat.cell2AB(Cell[:6])
995	covData = {}
996	if len(DisAglData.get('covData',{})):
997	covData = DisAglData['covData']
998	covMatrix = covData['covMatrix']
999	varyList = covData['varyList']
1000	pfx = str(DisAglData['pId'])+'::'
1001
1002	Factor = DisAglCtls['Factors']
1003	Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
1004	indices = (-2,-1,0,1,2)
1005	Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
1006	origAtoms = DisAglData['OrigAtoms']
1007	targAtoms = DisAglData['TargAtoms']
1008	AtomLabels = {}
1009	for Oatom in origAtoms:
1010	AtomLabels[Oatom[0]] = Oatom[1]
1011	for Oatom in targAtoms:
1012	AtomLabels[Oatom[0]] = Oatom[1]
1013	DistArray = {}
1014	AngArray = {}
1015	for iO,Oatom in enumerate(origAtoms):
1016	DistArray[Oatom[0]] = []
1017	AngArray[Oatom[0]] = []
1018	OxyzNames = ''
1019	IndBlist = []
1020	Dist = []
1021	Vect = []
1022	VectA = []
1023	angles = []
1024	for Tatom in targAtoms:
1025	Xvcov = []
1026	TxyzNames = ''
1027	if len(DisAglData.get('covData',{})):
1028	OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
1029	TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
1030	Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
1031	BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
1032	AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
1033	for [Txyz,Top,Tunit,Spn] in G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False):
1034	Dx = (Txyz-np.array(Oatom[3:6]))+Units
1035	dx = np.inner(Amat,Dx)
1036	dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
1037	IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
1038	if np.any(IndB):
1039	for indb in IndB:
1040	for i in range(len(indb)):
1041	if str(dx.T[indb][i]) not in IndBlist:
1042	IndBlist.append(str(dx.T[indb][i]))
1043	unit = Units[indb][i]
1044	tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
1045	sig = 0.0
1046	if len(Xvcov):
1047	pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
1048	sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
1049	Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
1050	if (Dist[-1][-2]-AsumR) <= 0.:
1051	Vect.append(dx.T[indb][i]/Dist[-1][-2])
1052	VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
1053	else:
1054	Vect.append([0.,0.,0.])
1055	VectA.append([])
1056	if dlg is not None:
1057	dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
1058	for D in Dist:
1059	DistArray[Oatom[0]].append(D[1:])
1060	Vect = np.array(Vect)
1061	angles = np.zeros((len(Vect),len(Vect)))
1062	angsig = np.zeros((len(Vect),len(Vect)))
1063	for i,veca in enumerate(Vect):
1064	if np.any(veca):
1065	for j,vecb in enumerate(Vect):
1066	if np.any(vecb):
1067	angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
1068	if i <= j: continue
1069	AngArray[Oatom[0]].append((i,j,
1070	G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
1071	return AtomLabels,DistArray,AngArray
1072
1073	def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
1074	'''Print distances and angles
1075
1076	:param dict DisAglCtls: contains distance/angle radii usually defined using
1077	:func:`GSASIIctrlGUI.DisAglDialog`
1078	:param dict DisAglData: contains phase data:
1079	Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
1080	for distance/angle origins and atoms to be used as targets.
1081	Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
1082	:param file out: file object for output. Defaults to sys.stdout.
1083	'''
1084	def MyPrint(s):
1085	out.write(s+'\n')
1086	# print(s,file=out) # use in Python 3
1087
1088	def ShowBanner(name):
1089	MyPrint(80'')
1090	MyPrint(' Interatomic Distances and Angles for phase '+name)
1091	MyPrint((80'')+'\n')
1092
1093	ShowBanner(DisAglCtls['Name'])
1094	SGData = DisAglData['SGData']
1095	SGtext,SGtable = G2spc.SGPrint(SGData)
1096	for line in SGtext: MyPrint(line)
1097	if len(SGtable) > 1:
1098	for i,item in enumerate(SGtable[::2]):
1099	if 2*i+1 == len(SGtable):
1100	line = ' %s'%(item.ljust(30))
1101	else:
1102	line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
1103	MyPrint(line)
1104	else:
1105	MyPrint(' ( 1) %s'%(SGtable[0])) #triclinic case
1106	Cell = DisAglData['Cell']
1107
1108	Amat,Bmat = G2lat.cell2AB(Cell[:6])
1109	covData = {}
1110	if len(DisAglData.get('covData',{})):
1111	covData = DisAglData['covData']
1112	pfx = str(DisAglData['pId'])+'::'
1113	A = G2lat.cell2A(Cell[:6])
1114	cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
1115	names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
1116	valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
1117	line = '\n Unit cell:'
1118	for name,vals in zip(names,valEsd):
1119	line += name+vals
1120	MyPrint(line)
1121	else:
1122	MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
1123	' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
1124	('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
1125
1126	AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
1127	origAtoms = DisAglData['OrigAtoms']
1128	for Oatom in origAtoms:
1129	i = Oatom[0]
1130	Dist = DistArray[i]
1131	nDist = len(Dist)
1132	angles = np.zeros((nDist,nDist))
1133	angsig = np.zeros((nDist,nDist))
1134	for k,j,tup in AngArray[i]:
1135	angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
1136	line = ''
1137	for i,x in enumerate(Oatom[3:6]):
1138	line += ('%12.5f'%x).rstrip('0')
1139	MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
1140	MyPrint(80'')
1141	line = ''
1142	for dist in Dist[:-1]:
1143	line += '%12s'%(AtomLabels[dist[0]].center(12))
1144	MyPrint(' To cell +(sym. op.) dist. '+line.rstrip())
1145	BVS = {}
1146	BVdat = {}
1147	Otyp = G2elem.FixValence(Oatom[2]).split('+')[0].split('-')[0]
1148	BVox = [BV for BV in atmdata.BVSoxid[Otyp] if '+' in BV]
1149	if len(BVox):
1150	BVS = {BV:0.0 for BV in BVox}
1151	BVdat = {BV:dict(zip(['O','F','Cl'],atmdata.BVScoeff[BV])) for BV in BVox}
1152	pvline = 'Bond Valence sums for: '
1153	for i,dist in enumerate(Dist):
1154	line = ''
1155	for j,angle in enumerate(angles[i][0:i]):
1156	sig = angsig[i][j]
1157	if angle:
1158	if sig:
1159	line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
1160	else:
1161	val = '%.3f'%(angle)
1162	line += '%12s'%(val.center(12))
1163	else:
1164	line += 12*' '
1165	if dist[4]: #sig exists!
1166	val = G2mth.ValEsd(dist[3],dist[4])
1167	else:
1168	val = '%8.4f'%(dist[3])
1169	if len(BVox):
1170	Tatm = G2elem.FixValence(DisAglData['TargAtoms'][dist[0]][2]).split('-')[0]
1171	if Tatm in ['O','F','Cl']:
1172	for BV in BVox:
1173	BVS[BV] += np.exp((BVdat[BV][Tatm]-dist[3])/0.37)
1174	tunit = '[%2d%2d%2d]'% dist[1]
1175	MyPrint((' %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
1176	if len(BVox):
1177	MyPrint(80'')
1178	for BV in BVox:
1179	pvline += ' %s: %.2f '%(BV,BVS[BV])
1180	MyPrint(pvline)
1181
1182	def DisAglTor(DATData):
1183	'Needs a doc string'
1184	SGData = DATData['SGData']
1185	Cell = DATData['Cell']
1186
1187	Amat,Bmat = G2lat.cell2AB(Cell[:6])
1188	covData = {}
1189	pfx = ''
1190	if 'covData' in DATData:
1191	covData = DATData['covData']
1192	pfx = str(DATData['pId'])+'::'
1193	Datoms = []
1194	Oatoms = []
1195	for i,atom in enumerate(DATData['Datoms']):
1196	symop = atom[-1].split('+')
1197	if len(symop) == 1:
1198	symop.append('0,0,0')
1199	symop[0] = int(symop[0])
1200	symop[1] = eval(symop[1])
1201	atom.append(symop)
1202	Datoms.append(atom)
1203	oatom = DATData['Oatoms'][i]
1204	names = ['','','']
1205	if pfx:
1206	names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
1207	oatom += [names,]
1208	Oatoms.append(oatom)
1209	atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
1210	if DATData['Natoms'] == 4: #torsion
1211	Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1212	G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
1213	G2fil.G2Print (' NB: Atom sequence determined by selection order')
1214	return # done with torsion
1215	elif DATData['Natoms'] == 3: #angle
1216	Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1217	G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
1218	G2fil.G2Print (' NB: Atom sequence determined by selection order')
1219	else: #2 atoms - distance
1220	Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1221	G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
1222
1223	def BestPlane(PlaneData):
1224	'Needs a doc string'
1225
1226	def ShowBanner(name):
1227	G2fil.G2Print (80'')
1228	G2fil.G2Print (' Best plane result for phase '+name)
1229	G2fil.G2Print (80'','\n')
1230
1231	ShowBanner(PlaneData['Name'])
1232
1233	Cell = PlaneData['Cell']
1234	Amat,Bmat = G2lat.cell2AB(Cell[:6])
1235	Atoms = PlaneData['Atoms']
1236	sumXYZ = np.zeros(3)
1237	XYZ = []
1238	Natoms = len(Atoms)
1239	for atom in Atoms:
1240	xyz = np.array(atom[3:6])
1241	XYZ.append(xyz)
1242	sumXYZ += xyz
1243	sumXYZ /= Natoms
1244	XYZ = np.array(XYZ)-sumXYZ
1245	XYZ = np.inner(Amat,XYZ).T
1246	Zmat = np.zeros((3,3))
1247	for i,xyz in enumerate(XYZ):
1248	Zmat += np.outer(xyz.T,xyz)
1249	G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
1250	Evec,Emat = nl.eig(Zmat)
1251	Evec = np.sqrt(Evec)/(Natoms-3)
1252	Order = np.argsort(Evec)
1253	XYZ = np.inner(XYZ,Emat.T).T
1254	XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
1255	G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
1256	G2fil.G2Print (' Name X Y Z')
1257	for i,xyz in enumerate(XYZ):
1258	G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
1259	G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
1260
1261	def main():
1262	'Called to run a refinement when this module is executed '
1263	starttime = time.time()
1264	arg = sys.argv
1265	if len(arg) > 1:
1266	GPXfile = arg[1]
1267	if not ospath.exists(GPXfile):
1268	G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
1269	exit()
1270	else:
1271	G2fil.G2Print ('ERROR - missing filename')
1272	exit()
1273	# TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
1274	Refine(GPXfile,None)
1275	G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
1276
1277	if __name__ == '__main__':
1278	GSASIIpath.InvokeDebugOpts()
1279	main()

Note: See TracBrowser for help on using the repository browser.

Download in other formats: