source: trunk/GSASIIstrMain.py @ 4534

Last change on this file since 4534 was 4534, checked in by toby, 3 years ago

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