source: trunk/GSASIIstrMain.py @ 4600

Last change on this file since 4600 was 4588, checked in by toby, 5 years ago

use G2VarObj for param limits; add more info to seq. ref. done dialog; show Frozen in show LS parameters
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 46.3 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2020-10-11 17:07:42 +0000 (Sun, 11 Oct 2020) $
9# $Author: toby $
10# $Revision: 4588 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4588 2020-10-11 17:07:42Z 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: 4588 $")
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: %s\n'%histogram)
553        printFile.write(135*'-'+'\n')
554        # remove frozen vars
555        if 'parmFrozen' not in Controls:
556            Controls['parmFrozen'] = {}
557        if histogram not in Controls['parmFrozen']: 
558            Controls['parmFrozen'][histogram] = []
559        parmFrozenList = Controls['parmFrozen'][histogram]
560        frozenList = [i for i in varyList if i in parmFrozenList]
561        if len(frozenList) != 0: 
562           varyList = [i for i in varyList if i not in parmFrozenList]
563           printFile.write(
564               ' The following refined variables have previously been frozen due to exceeding limits\n\t:{}\n'
565               .format(frozenList))
566        try:
567            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
568                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
569                refPlotUpdate=refPlotUpdate)
570            G2fil.G2Print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f, last shft/sig = %.4f'%(
571                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'],Rvals.get('Max shft/sig',np.nan)))
572            if Rvals.get('lamMax',0) >= 10.:
573                msgs['steepestNum'] += 1
574            if Rvals.get('Max shft/sig'):
575                msgs['maxshift/sigma'].append(Rvals['Max shft/sig'])
576            # add the uncertainties into the esd dictionary (sigDict)
577            if not IfOK:
578                G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
579                break
580            sigDict = dict(zip(varyList,sig))
581            # the uncertainties for dependent constrained parms into the esd dict
582            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
583            # check for variables outside their allowed range, reset and freeze them
584            frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
585            if len(frozen) > 0:
586               msg = ('Hist {}: {} variables were outside limits and were frozen (now {} frozen total)'
587                   .format(ihst,len(frozen),len(parmFrozenList)))
588               G2fil.G2Print(msg)
589               printFile.write(msg+'\n')
590               for p in frozen:
591                   if p not in varyList:
592                       print('Frozen Warning: {} not in varyList. This should not happen!'.format(p))
593                       continue
594                   i = varyList.index(p)
595                   result[0][i] = parmDict[p]
596                   sig[i] = -0.1
597            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
598            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
599            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
600            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
601            histRefData = {
602                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
603                'varyListStart':varyListStart,
604                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
605                'newCellDict':newCellDict,'depParmDict':depParmDict,
606                'constraintInfo':constraintInfo,
607                'parmDict':parmDict,
608                }
609            SeqResult[histogram] = histRefData
610            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
611#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
612            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
613            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
614            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,Controls['parmFrozen'])
615            NewparmDict = {}
616            # make dict of varied parameters in current histogram, renamed to
617            # next histogram, for use in next refinement.
618            if Controls['Copy2Next'] and ihst < len(histNames)-1:
619                hId = Histo[histogram]['hId'] # current histogram
620                nexthId = Histograms[histNames[ihst+1]]['hId']
621                for parm in set(list(varyList)+list(varyListStart)):
622                    items = parm.split(':')
623                    if len(items) < 3: 
624                        continue
625                    if str(hId) in items[1]:
626                        items[1] = str(nexthId)
627                        newparm = ':'.join(items)
628                        NewparmDict[newparm] = parmDict[parm]
629                    else:
630                        if items[2].startswith('dA'): parm = parm.replace(':dA',':A') 
631                        NewparmDict[parm] = parmDict[parm]
632                   
633        except G2obj.G2RefineCancel as Msg:
634            printFile.close()
635            G2fil.G2Print (' ***** Refinement stopped *****')
636            return False,Msg.msg
637        except G2obj.G2Exception as Msg:  # cell metric error, others?
638            printFile.close()
639            G2fil.G2Print (' ***** Refinement error *****')
640            return False,Msg.msg
641        if GSASIIpath.GetConfigValue('Show_timing'):
642            t2 = time.time()
643            G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
644            t1 = t2
645    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
646    try:
647        G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
648    except Exception as msg:
649        print('Error reading Sequential results')
650        if GSASIIpath.GetConfigValue('debug'):
651            import traceback
652            print(traceback.format_exc())       
653    postFrozenCount = 0
654    for h in Controls['parmFrozen']:
655        if h == 'FrozenList': continue
656        postFrozenCount += len(Controls['parmFrozen'][h])
657    if postFrozenCount:
658        msgs['Frozen'] = 'Ending refinement with {} Frozen variables ({} added here)\n'.format(postFrozenCount,postFrozenCount-preFrozenCount)
659        printFile.write('\n '+msg)
660    printFile.close()
661    G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
662    G2fil.G2Print (' ***** Sequential refinement successful *****')
663    return True,msgs
664
665def dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList):
666    '''Find variables in the parameters dict that are outside the ranges
667    (in parmMinDict and parmMaxDict) and set them to the limits values.
668    Add any such variables into the list of frozen variable
669    (parmFrozenList). Returns a list of newly frozen variables, if any.
670    '''
671    parmMinDict = Controls.get('parmMinDict',{})
672    parmMaxDict = Controls.get('parmMaxDict',{})
673    freeze = []
674    if parmMinDict or parmMaxDict:
675        for name in varyList:
676            if name not in parmDict: continue
677            n,val = G2obj.prmLookup(name,parmMinDict)
678            if n is not None:
679                if parmDict[name] < parmMinDict[n]:
680                    parmDict[name] = parmMinDict[n]
681                    sigDict[name] = 0.0
682                    freeze.append(name)
683                    continue
684            n,val = G2obj.prmLookup(name,parmMaxDict)
685            if n is not None:
686                if parmDict[name] > parmMaxDict[n]:
687                    parmDict[name] = parmMaxDict[n]
688                    sigDict[name] = 0.0
689                    freeze.append(name)
690                    continue
691        for v in freeze:
692            if v not in parmFrozenList:
693                parmFrozenList.append(v)
694    return freeze
695
696def RetDistAngle(DisAglCtls,DisAglData,dlg=None):
697    '''Compute and return distances and angles
698
699    :param dict DisAglCtls: contains distance/angle radii usually defined using
700       :func:`GSASIIctrlGUI.DisAglDialog`
701    :param dict DisAglData: contains phase data:
702       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
703       for distance/angle origins and atoms to be used as targets.
704       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
705
706    :returns: AtomLabels,DistArray,AngArray where:
707
708      **AtomLabels** is a dict of atom labels, keys are the atom number
709
710      **DistArray** is a dict keyed by the origin atom number where the value is a list
711      of distance entries. The value for each distance is a list containing:
712
713        0) the target atom number (int);
714        1) the unit cell offsets added to x,y & z (tuple of int values)
715        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
716        3) an interatomic distance in A (float)
717        4) an uncertainty on the distance in A or 0.0 (float)
718
719      **AngArray** is a dict keyed by the origin (central) atom number where
720      the value is a list of
721      angle entries. The value for each angle entry consists of three values:
722
723        0) a distance item reference for one neighbor (int)
724        1) a distance item reference for a second neighbor (int)
725        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
726
727      The AngArray distance reference items refer directly to the index of the items in the
728      DistArray item for the list of distances for the central atom.
729    '''
730    import numpy.ma as ma
731
732    SGData = DisAglData['SGData']
733    Cell = DisAglData['Cell']
734
735    Amat,Bmat = G2lat.cell2AB(Cell[:6])
736    covData = {}
737    if len(DisAglData.get('covData',{})):
738        covData = DisAglData['covData']
739        covMatrix = covData['covMatrix']
740        varyList = covData['varyList']
741        pfx = str(DisAglData['pId'])+'::'
742
743    Factor = DisAglCtls['Factors']
744    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
745    indices = (-2,-1,0,1,2)
746    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
747    origAtoms = DisAglData['OrigAtoms']
748    targAtoms = DisAglData['TargAtoms']
749    AtomLabels = {}
750    for Oatom in origAtoms:
751        AtomLabels[Oatom[0]] = Oatom[1]
752    for Oatom in targAtoms:
753        AtomLabels[Oatom[0]] = Oatom[1]
754    DistArray = {}
755    AngArray = {}
756    for iO,Oatom in enumerate(origAtoms):
757        DistArray[Oatom[0]] = []
758        AngArray[Oatom[0]] = []
759        OxyzNames = ''
760        IndBlist = []
761        Dist = []
762        Vect = []
763        VectA = []
764        angles = []
765        for Tatom in targAtoms:
766            Xvcov = []
767            TxyzNames = ''
768            if len(DisAglData.get('covData',{})):
769                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
770                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
771                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
772            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
773            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
774            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
775            for [Txyz,Top,Tunit,Spn] in result:
776                Dx = (Txyz-np.array(Oatom[3:6]))+Units
777                dx = np.inner(Amat,Dx)
778                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
779                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
780                if np.any(IndB):
781                    for indb in IndB:
782                        for i in range(len(indb)):
783                            if str(dx.T[indb][i]) not in IndBlist:
784                                IndBlist.append(str(dx.T[indb][i]))
785                                unit = Units[indb][i]
786                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
787                                sig = 0.0
788                                if len(Xvcov):
789                                    pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
790                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
791                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
792                                if (Dist[-1][-2]-AsumR) <= 0.:
793                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
794                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
795                                else:
796                                    Vect.append([0.,0.,0.])
797                                    VectA.append([])
798        if dlg is not None:
799            dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
800        for D in Dist:
801            DistArray[Oatom[0]].append(D[1:])
802        Vect = np.array(Vect)
803        angles = np.zeros((len(Vect),len(Vect)))
804        angsig = np.zeros((len(Vect),len(Vect)))
805        for i,veca in enumerate(Vect):
806            if np.any(veca):
807                for j,vecb in enumerate(Vect):
808                    if np.any(vecb):
809                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
810                        if i <= j: continue
811                        AngArray[Oatom[0]].append((i,j,
812                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
813    return AtomLabels,DistArray,AngArray
814
815def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
816    '''Print distances and angles
817
818    :param dict DisAglCtls: contains distance/angle radii usually defined using
819       :func:`GSASIIctrlGUI.DisAglDialog`
820    :param dict DisAglData: contains phase data:
821       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
822       for distance/angle origins and atoms to be used as targets.
823       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
824    :param file out: file object for output. Defaults to sys.stdout.
825    '''
826    def MyPrint(s):
827        out.write(s+'\n')
828        # print(s,file=out) # use in Python 3
829
830    def ShowBanner(name):
831        MyPrint(80*'*')
832        MyPrint('   Interatomic Distances and Angles for phase '+name)
833        MyPrint((80*'*')+'\n')
834
835    ShowBanner(DisAglCtls['Name'])
836    SGData = DisAglData['SGData']
837    SGtext,SGtable = G2spc.SGPrint(SGData)
838    for line in SGtext: MyPrint(line)
839    if len(SGtable) > 1:
840        for i,item in enumerate(SGtable[::2]):
841            if 2*i+1 == len(SGtable):
842                line = ' %s'%(item.ljust(30))
843            else:
844                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
845            MyPrint(line)
846    else:
847        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
848    Cell = DisAglData['Cell']
849
850    Amat,Bmat = G2lat.cell2AB(Cell[:6])
851    covData = {}
852    if len(DisAglData.get('covData',{})):
853        covData = DisAglData['covData']
854        pfx = str(DisAglData['pId'])+'::'
855        A = G2lat.cell2A(Cell[:6])
856        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
857        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
858        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
859        line = '\n Unit cell:'
860        for name,vals in zip(names,valEsd):
861            line += name+vals
862        MyPrint(line)
863    else:
864        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
865            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
866            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
867
868    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
869    origAtoms = DisAglData['OrigAtoms']
870    for Oatom in origAtoms:
871        i = Oatom[0]
872        Dist = DistArray[i]
873        nDist = len(Dist)
874        angles = np.zeros((nDist,nDist))
875        angsig = np.zeros((nDist,nDist))
876        for k,j,tup in AngArray[i]:
877            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
878        line = ''
879        for i,x in enumerate(Oatom[3:6]):
880            line += ('%12.5f'%x).rstrip('0')
881        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
882        MyPrint(80*'*')
883        line = ''
884        for dist in Dist[:-1]:
885            line += '%12s'%(AtomLabels[dist[0]].center(12))
886        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
887        for i,dist in enumerate(Dist):
888            line = ''
889            for j,angle in enumerate(angles[i][0:i]):
890                sig = angsig[i][j]
891                if angle:
892                    if sig:
893                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
894                    else:
895                        val = '%.3f'%(angle)
896                        line += '%12s'%(val.center(12))
897                else:
898                    line += 12*' '
899            if dist[4]:            #sig exists!
900                val = G2mth.ValEsd(dist[3],dist[4])
901            else:
902                val = '%8.4f'%(dist[3])
903            tunit = '[%2d%2d%2d]'% dist[1]
904            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
905
906def DisAglTor(DATData):
907    'Needs a doc string'
908    SGData = DATData['SGData']
909    Cell = DATData['Cell']
910
911    Amat,Bmat = G2lat.cell2AB(Cell[:6])
912    covData = {}
913    pfx = ''
914    if 'covData' in DATData:
915        covData = DATData['covData']
916        pfx = str(DATData['pId'])+'::'
917    Datoms = []
918    Oatoms = []
919    for i,atom in enumerate(DATData['Datoms']):
920        symop = atom[-1].split('+')
921        if len(symop) == 1:
922            symop.append('0,0,0')
923        symop[0] = int(symop[0])
924        symop[1] = eval(symop[1])
925        atom.append(symop)
926        Datoms.append(atom)
927        oatom = DATData['Oatoms'][i]
928        names = ['','','']
929        if pfx:
930            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
931        oatom += [names,]
932        Oatoms.append(oatom)
933    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
934    if DATData['Natoms'] == 4:  #torsion
935        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
936        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
937        G2fil.G2Print (' NB: Atom sequence determined by selection order')
938        return      # done with torsion
939    elif DATData['Natoms'] == 3:  #angle
940        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
941        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
942        G2fil.G2Print (' NB: Atom sequence determined by selection order')
943    else:   #2 atoms - distance
944        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
945        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
946
947def BestPlane(PlaneData):
948    'Needs a doc string'
949
950    def ShowBanner(name):
951        G2fil.G2Print (80*'*')
952        G2fil.G2Print ('   Best plane result for phase '+name)
953        G2fil.G2Print (80*'*','\n')
954
955    ShowBanner(PlaneData['Name'])
956
957    Cell = PlaneData['Cell']
958    Amat,Bmat = G2lat.cell2AB(Cell[:6])
959    Atoms = PlaneData['Atoms']
960    sumXYZ = np.zeros(3)
961    XYZ = []
962    Natoms = len(Atoms)
963    for atom in Atoms:
964        xyz = np.array(atom[3:6])
965        XYZ.append(xyz)
966        sumXYZ += xyz
967    sumXYZ /= Natoms
968    XYZ = np.array(XYZ)-sumXYZ
969    XYZ = np.inner(Amat,XYZ).T
970    Zmat = np.zeros((3,3))
971    for i,xyz in enumerate(XYZ):
972        Zmat += np.outer(xyz.T,xyz)
973    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
974    Evec,Emat = nl.eig(Zmat)
975    Evec = np.sqrt(Evec)/(Natoms-3)
976    Order = np.argsort(Evec)
977    XYZ = np.inner(XYZ,Emat.T).T
978    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
979    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
980    G2fil.G2Print (' Name         X         Y         Z')
981    for i,xyz in enumerate(XYZ):
982        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
983    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
984
985def main():
986    'Called to run a refinement when this module is executed '
987    starttime = time.time()
988    arg = sys.argv
989    if len(arg) > 1:
990        GPXfile = arg[1]
991        if not ospath.exists(GPXfile):
992            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
993            exit()
994    else:
995        G2fil.G2Print ('ERROR - missing filename')
996        exit()
997    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
998    Refine(GPXfile,None)
999    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
1000
1001if __name__ == '__main__':
1002    GSASIIpath.InvokeDebugOpts()
1003    main()
Note: See TracBrowser for help on using the repository browser.