Context Navigation

source: trunk/GSASIIstrMain.py @ 4850

Last change on this file since 4850 was 4850, checked in by vondreele, 3 years ago
fix LeBail? test to skip HKLF histograms in G2strMain fix DBW access for old gpx files in G2dataGUI fix Quaternion access for old gpx files in G2plot add new testSytSym.py routine
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 57.9 KB

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

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

Download in other formats: