source: trunk/GSASIIstrMain.py @ 4578

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

docs for freezing parameters started + docs cleanup; start scriptable for freezing params; record initial chi2; Show more post refinement info; noted but unfixed bkg GUI bug

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