Context Navigation

source: trunk/GSASIIstrMain.py @ 4785

Last change on this file since 4785 was 4783, checked in by toby, 4 years ago
wx4.1 updates; new CIF data item; scripting: peak fit details returned; improve multistring gui
Property svn:eol-style set to `native` Property svn:keywords set to `Date Author Revision URL Id`
File size: 53.7 KB

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

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

Download in other formats: