source: trunk/GSASIIstrMain.py @ 4628

Last change on this file since 4628 was 4628, checked in by toby, 2 years ago

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