source: trunk/GSASIIstrMain.py @ 4451

Last change on this file since 4451 was 4451, checked in by vondreele, 3 years ago

Add Histogram distances & angles to Phase/Atom/Compute? menu
fix getAtomSelections for macromolecules & add getAtomPtrs
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 40.8 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2020-06-01 17:40:11 +0000 (Mon, 01 Jun 2020) $
9# $Author: vondreele $
10# $Revision: 4451 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4451 2020-06-01 17:40:11Z vondreele $
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: 4451 $")
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    ifSeq = False
256    printFile.write('\n Refinement results:\n')
257    printFile.write(135*'-'+'\n')
258    try:
259        covData = {}
260        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
261            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
262            refPlotUpdate=refPlotUpdate)
263        if IfOK:
264            sigDict = dict(zip(varyList,sig))
265            newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
266            newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
267            covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
268                       'varyListStart':varyListStart,
269                       'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
270                       'newCellDict':newCellDict,'freshCOV':True}
271            # add the uncertainties into the esd dictionary (sigDict)
272            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
273            G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
274            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
275            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
276            G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
277            G2stIO.PrintISOmodes(printFile,Phases,parmDict,sigDict)
278            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
279            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)
280            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,makeBack)
281            printFile.close()
282            G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
283            G2fil.G2Print (' ***** Refinement successful *****')
284        else:
285            G2fil.G2Print ('****ERROR - Refinement failed')
286            raise G2obj.G2Exception('****ERROR - Refinement failed')
287    except G2obj.G2RefineCancel as Msg:
288        printFile.close()
289        G2fil.G2Print (' ***** Refinement stopped *****')
290        return False,Msg.msg
291    except G2obj.G2Exception as Msg:  # cell metric error, others?
292        printFile.close()
293        G2fil.G2Print (' ***** Refinement error *****')
294        return False,Msg.msg
295
296#for testing purposes, create a file for testderiv
297    if GSASIIpath.GetConfigValue('debug'):   # and IfOK:
298#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
299        fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
300        cPickle.dump(result[0],fl,1)
301        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
302        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
303        cPickle.dump(parmDict,fl,1)
304        cPickle.dump(varyList,fl,1)
305        cPickle.dump(calcControls,fl,1)
306        cPickle.dump(pawleyLookup,fl,1)
307        fl.close()
308    if dlg:
309        return True,Rvals
310
311def phaseCheck(phaseVary,Phases,histogram):
312    '''
313    Removes unused parameters from phase varylist if phase not in histogram
314    #TODO - implement "Fix FXU" for seq refinement here - done?
315    '''
316    NewVary = []
317    for phase in Phases:
318        if histogram not in Phases[phase]['Histograms']: continue
319        if Phases[phase]['Histograms'][histogram]['Use']:
320            pId = Phases[phase]['pId']
321            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
322            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
323            if 'F' in FixVals:
324                newVary = [item for item in newVary if not 'Afrac' in item]
325            if 'X' in FixVals:
326                newVary = [item for item in newVary if not 'dA' in item]
327            if 'U' in FixVals:
328                newVary = [item for item in newVary if not 'AU' in item]
329            if 'M' in FixVals:
330                newVary = [item for item in newVary if not 'AM' in item]
331            NewVary += newVary
332    return NewVary
333
334def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
335    '''Perform a sequential refinement -- cycles through all selected histgrams,
336    one at a time
337    '''
338    import GSASIImpsubs as G2mp
339    G2mp.InitMP()
340    import pytexture as ptx
341    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
342
343    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
344    G2fil.G2Print ('Starting Sequential Refinement')
345    G2stIO.ShowBanner(printFile)
346    Controls = G2stIO.GetControls(GPXfile)
347    G2stIO.ShowControls(Controls,printFile,SeqRef=True)
348    restraintDict = G2stIO.GetRestraints(GPXfile)
349    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
350    if not Phases:
351        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
352        G2fil.G2Print (' *** Refine aborted ***')
353        return False,'No phases'
354    if not Histograms:
355        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
356        G2fil.G2Print (' *** Refine aborted ***')
357        return False,'No data'
358    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
359    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
360    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
361    G2mv.InitVars()
362    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
363         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,
364                                    Print=False,pFile=printFile,seqRef=True)
365    for item in phaseVary:
366        if '::A0' in item:
367            G2fil.G2Print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
368            G2fil.G2Print ('****           instead use the Dij parameters for each powder histogram            ****')
369            return False,'Lattice parameter refinement error - see console message'
370        if '::C(' in item:
371            G2fil.G2Print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
372            G2fil.G2Print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
373            return False,'Phase texture refinement error - see console message'
374    if 'Seq Data' in Controls:
375        histNames = Controls['Seq Data']
376    else: # patch from before Controls['Seq Data'] was implemented?
377        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
378    if Controls.get('Reverse Seq'):
379        histNames.reverse()
380    SeqResult = G2stIO.GetSeqResult(GPXfile)
381#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
382    Histo = {}
383    NewparmDict = {}
384    G2stIO.SetupSeqSavePhases(GPXfile)
385    for ihst,histogram in enumerate(histNames):
386        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
387        G2fil.G2Print('\nRefining with '+str(histogram))
388        G2mv.InitVars()
389        (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
390             FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
391                 Phases,restraintDict,rbIds,
392                 Print=False,pFile=printFile,seqRef=True)
393        ifPrint = False
394        if dlg:
395            dlg.SetTitle('Residual for histogram '+str(ihst))
396        calcControls = {}
397        calcControls['atomIndx'] = atomIndx
398        calcControls['Natoms'] = Natoms
399        calcControls['FFtables'] = FFtables
400        calcControls['BLtables'] = BLtables
401        calcControls['MFtables'] = MFtables
402        calcControls['maxSSwave'] = maxSSwave
403        if histogram not in Histograms:
404            G2fil.G2Print("Error: not found!")
405            continue
406    #TODO - implement "Fix FXU" for seq refinement here - done?
407        hId = Histograms[histogram]['hId']
408        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
409        Histo = {histogram:Histograms[histogram],}
410        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
411        calcControls.update(controlDict)
412        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
413        calcControls.update(controlDict)
414        varyList = rbVary+redphaseVary+hapVary+histVary
415#        if not ihst:
416            # save the initial vary list, but without histogram numbers on parameters
417        saveVaryList = varyList[:]
418        for i,item in enumerate(saveVaryList):
419            items = item.split(':')
420            if items[1]:
421                items[1] = ''
422            item = ':'.join(items)
423            saveVaryList[i] = item
424        if not ihst:
425            SeqResult['varyList'] = saveVaryList
426        else:
427            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
428        parmDict = {}
429        parmDict.update(rbDict)
430        parmDict.update(phaseDict)
431        parmDict.update(hapDict)
432        parmDict.update(histDict)
433        if Controls['Copy2Next']:   # update with parms from last histogram
434            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
435            for parm in NewparmDict:
436                if parm in parmDict:
437                    parmDict[parm] = NewparmDict[parm]
438        elif histogram in SeqResult:  # update phase from last seq ref
439            NewparmDict = SeqResult[histogram].get('parmDict',{})
440            for parm in NewparmDict:
441                if '::' in parm and parm in parmDict:
442                    parmDict[parm] = NewparmDict[parm]
443           
444        G2stIO.GetFprime(calcControls,Histo)
445        # do constraint processing
446        #reload(G2mv) # debug
447        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
448        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
449        msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
450        if msg:
451            return False,'Unable to interpret multiplier(s): '+msg
452        try:
453            groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
454#            if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
455#                G2mv.VarRemapShow(varyList,True))
456#            print('DependentVars',G2mv.GetDependentVars())
457#            print('IndependentVars',G2mv.GetIndependentVars())
458            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
459        except G2mv.ConstraintException:
460            G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
461            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
462            #print 'Errors',errmsg
463            #if warnmsg: print 'Warnings',warnmsg
464            return False,' Constraint error'
465        #print G2mv.VarRemapShow(varyList)
466        if not ihst:
467            # first histogram to refine against
468            firstVaryList = []
469            for item in varyList:
470                items = item.split(':')
471                if items[1]:
472                    items[1] = ''
473                item = ':'.join(items)
474                firstVaryList.append(item)
475            newVaryList = firstVaryList
476        else:
477            newVaryList = []
478            for item in varyList:
479                items = item.split(':')
480                if items[1]:
481                    items[1] = ''
482                item = ':'.join(items)
483                newVaryList.append(item)
484        if newVaryList != firstVaryList and Controls['Copy2Next']:
485            # variable lists are expected to match between sequential refinements when Copy2Next is on
486            #print '**** ERROR - variable list for this histogram does not match previous'
487            #print '     Copy of variables is not possible'
488            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
489            combined = list(set(firstVaryList+newVaryList))
490            c = [var for var in combined if var not in newVaryList]
491            p = [var for var in combined if var not in firstVaryList]
492            G2fil.G2Print('*** Variables change ***')
493            for typ,vars in [('Removed',c),('Added',p)]:
494                line = '  '+typ+': '
495                if vars:
496                    for var in vars:
497                        if len(line) > 70:
498                            G2fil.G2Print(line)
499                            line = '    '
500                        line += var + ', '
501                else:
502                        line += 'none, '
503                G2fil.G2Print(line[:-2])
504            firstVaryList = newVaryList
505
506        ifSeq = True
507        printFile.write('\n Refinement results for histogram: %s\n'%histogram)
508        printFile.write(135*'-'+'\n')
509        try:
510            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
511                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
512                refPlotUpdate=refPlotUpdate)
513            G2fil.G2Print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f, last shft/sig = %.4f'%(
514                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'],Rvals.get('Max shft/sig',np.nan)))
515            # add the uncertainties into the esd dictionary (sigDict)
516            if not IfOK:
517                G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
518                break
519            sigDict = dict(zip(varyList,sig))
520            # the uncertainties for dependent constrained parms into the esd dict
521            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
522
523            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
524            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
525            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
526            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
527            histRefData = {
528                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
529                'varyListStart':varyListStart,
530                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
531                'newCellDict':newCellDict,'depParmDict':depParmDict,
532                'constraintInfo':constraintInfo,
533                'parmDict':parmDict}
534            SeqResult[histogram] = histRefData
535            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
536#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
537            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
538            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
539            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData)
540            NewparmDict = {}
541            # make dict of varied parameters in current histogram, renamed to
542            # next histogram, for use in next refinement.
543            if Controls['Copy2Next'] and ihst < len(histNames)-1:
544                hId = Histo[histogram]['hId'] # current histogram
545                nexthId = Histograms[histNames[ihst+1]]['hId']
546                for parm in set(list(varyList)+list(varyListStart)):
547                    items = parm.split(':')
548                    if len(items) < 3: 
549                        continue
550                    if str(hId) in items[1]:
551                        items[1] = str(nexthId)
552                        newparm = ':'.join(items)
553                        NewparmDict[newparm] = parmDict[parm]
554                    else:
555                        if items[2].startswith('dA'): parm = parm.replace(':dA',':A') 
556                        NewparmDict[parm] = parmDict[parm]
557                   
558        except G2obj.G2RefineCancel as Msg:
559            printFile.close()
560            G2fil.G2Print (' ***** Refinement stopped *****')
561            return False,Msg.msg
562        except G2obj.G2Exception as Msg:  # cell metric error, others?
563            printFile.close()
564            G2fil.G2Print (' ***** Refinement error *****')
565            return False,Msg.msg
566        if GSASIIpath.GetConfigValue('Show_timing'):
567            t2 = time.time()
568            G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
569            t1 = t2
570    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
571    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
572    printFile.close()
573    G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
574    G2fil.G2Print (' ***** Sequential refinement successful *****')
575    return True,'Success'
576
577def RetDistAngle(DisAglCtls,DisAglData,dlg=None):
578    '''Compute and return distances and angles
579
580    :param dict DisAglCtls: contains distance/angle radii usually defined using
581       :func:`GSASIIctrlGUI.DisAglDialog`
582    :param dict DisAglData: contains phase data:
583       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
584       for distance/angle origins and atoms to be used as targets.
585       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
586
587    :returns: AtomLabels,DistArray,AngArray where:
588
589      **AtomLabels** is a dict of atom labels, keys are the atom number
590
591      **DistArray** is a dict keyed by the origin atom number where the value is a list
592      of distance entries. The value for each distance is a list containing:
593
594        0) the target atom number (int);
595        1) the unit cell offsets added to x,y & z (tuple of int values)
596        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
597        3) an interatomic distance in A (float)
598        4) an uncertainty on the distance in A or 0.0 (float)
599
600      **AngArray** is a dict keyed by the origin (central) atom number where
601      the value is a list of
602      angle entries. The value for each angle entry consists of three values:
603
604        0) a distance item reference for one neighbor (int)
605        1) a distance item reference for a second neighbor (int)
606        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
607
608      The AngArray distance reference items refer directly to the index of the items in the
609      DistArray item for the list of distances for the central atom.
610    '''
611    import numpy.ma as ma
612
613    SGData = DisAglData['SGData']
614    Cell = DisAglData['Cell']
615
616    Amat,Bmat = G2lat.cell2AB(Cell[:6])
617    covData = {}
618    if len(DisAglData.get('covData',{})):
619        covData = DisAglData['covData']
620        covMatrix = covData['covMatrix']
621        varyList = covData['varyList']
622        pfx = str(DisAglData['pId'])+'::'
623
624    Factor = DisAglCtls['Factors']
625    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
626    indices = (-2,-1,0,1,2)
627    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
628    origAtoms = DisAglData['OrigAtoms']
629    targAtoms = DisAglData['TargAtoms']
630    AtomLabels = {}
631    for Oatom in origAtoms:
632        AtomLabels[Oatom[0]] = Oatom[1]
633    for Oatom in targAtoms:
634        AtomLabels[Oatom[0]] = Oatom[1]
635    DistArray = {}
636    AngArray = {}
637    for iO,Oatom in enumerate(origAtoms):
638        DistArray[Oatom[0]] = []
639        AngArray[Oatom[0]] = []
640        OxyzNames = ''
641        IndBlist = []
642        Dist = []
643        Vect = []
644        VectA = []
645        angles = []
646        for Tatom in targAtoms:
647            Xvcov = []
648            TxyzNames = ''
649            if len(DisAglData.get('covData',{})):
650                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
651                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
652                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
653            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
654            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
655            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
656            for [Txyz,Top,Tunit,Spn] in result:
657                Dx = (Txyz-np.array(Oatom[3:6]))+Units
658                dx = np.inner(Amat,Dx)
659                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
660                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
661                if np.any(IndB):
662                    for indb in IndB:
663                        for i in range(len(indb)):
664                            if str(dx.T[indb][i]) not in IndBlist:
665                                IndBlist.append(str(dx.T[indb][i]))
666                                unit = Units[indb][i]
667                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
668                                sig = 0.0
669                                if len(Xvcov):
670                                    pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
671                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
672                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
673                                if (Dist[-1][-2]-AsumR) <= 0.:
674                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
675                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
676                                else:
677                                    Vect.append([0.,0.,0.])
678                                    VectA.append([])
679            if dlg is not None:
680                dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
681        for D in Dist:
682            DistArray[Oatom[0]].append(D[1:])
683        Vect = np.array(Vect)
684        angles = np.zeros((len(Vect),len(Vect)))
685        angsig = np.zeros((len(Vect),len(Vect)))
686        for i,veca in enumerate(Vect):
687            if np.any(veca):
688                for j,vecb in enumerate(Vect):
689                    if np.any(vecb):
690                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
691                        if i <= j: continue
692                        AngArray[Oatom[0]].append((i,j,
693                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
694    return AtomLabels,DistArray,AngArray
695
696def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
697    '''Print distances and angles
698
699    :param dict DisAglCtls: contains distance/angle radii usually defined using
700       :func:`GSASIIctrlGUI.DisAglDialog`
701    :param dict DisAglData: contains phase data:
702       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
703       for distance/angle origins and atoms to be used as targets.
704       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
705    :param file out: file object for output. Defaults to sys.stdout.
706    '''
707    def MyPrint(s):
708        out.write(s+'\n')
709        # print(s,file=out) # use in Python 3
710
711    def ShowBanner(name):
712        MyPrint(80*'*')
713        MyPrint('   Interatomic Distances and Angles for phase '+name)
714        MyPrint((80*'*')+'\n')
715
716    ShowBanner(DisAglCtls['Name'])
717    SGData = DisAglData['SGData']
718    SGtext,SGtable = G2spc.SGPrint(SGData)
719    for line in SGtext: MyPrint(line)
720    if len(SGtable) > 1:
721        for i,item in enumerate(SGtable[::2]):
722            if 2*i+1 == len(SGtable):
723                line = ' %s'%(item.ljust(30))
724            else:
725                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
726            MyPrint(line)
727    else:
728        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
729    Cell = DisAglData['Cell']
730
731    Amat,Bmat = G2lat.cell2AB(Cell[:6])
732    covData = {}
733    if len(DisAglData.get('covData',{})):
734        covData = DisAglData['covData']
735        pfx = str(DisAglData['pId'])+'::'
736        A = G2lat.cell2A(Cell[:6])
737        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
738        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
739        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
740        line = '\n Unit cell:'
741        for name,vals in zip(names,valEsd):
742            line += name+vals
743        MyPrint(line)
744    else:
745        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
746            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
747            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
748
749    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
750    origAtoms = DisAglData['OrigAtoms']
751    for Oatom in origAtoms:
752        i = Oatom[0]
753        Dist = DistArray[i]
754        nDist = len(Dist)
755        angles = np.zeros((nDist,nDist))
756        angsig = np.zeros((nDist,nDist))
757        for k,j,tup in AngArray[i]:
758            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
759        line = ''
760        for i,x in enumerate(Oatom[3:6]):
761            line += ('%12.5f'%x).rstrip('0')
762        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
763        MyPrint(80*'*')
764        line = ''
765        for dist in Dist[:-1]:
766            line += '%12s'%(AtomLabels[dist[0]].center(12))
767        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
768        for i,dist in enumerate(Dist):
769            line = ''
770            for j,angle in enumerate(angles[i][0:i]):
771                sig = angsig[i][j]
772                if angle:
773                    if sig:
774                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
775                    else:
776                        val = '%.3f'%(angle)
777                        line += '%12s'%(val.center(12))
778                else:
779                    line += 12*' '
780            if dist[4]:            #sig exists!
781                val = G2mth.ValEsd(dist[3],dist[4])
782            else:
783                val = '%8.4f'%(dist[3])
784            tunit = '[%2d%2d%2d]'% dist[1]
785            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
786
787def DisAglTor(DATData):
788    'Needs a doc string'
789    SGData = DATData['SGData']
790    Cell = DATData['Cell']
791
792    Amat,Bmat = G2lat.cell2AB(Cell[:6])
793    covData = {}
794    pfx = ''
795    if 'covData' in DATData:
796        covData = DATData['covData']
797        pfx = str(DATData['pId'])+'::'
798    Datoms = []
799    Oatoms = []
800    for i,atom in enumerate(DATData['Datoms']):
801        symop = atom[-1].split('+')
802        if len(symop) == 1:
803            symop.append('0,0,0')
804        symop[0] = int(symop[0])
805        symop[1] = eval(symop[1])
806        atom.append(symop)
807        Datoms.append(atom)
808        oatom = DATData['Oatoms'][i]
809        names = ['','','']
810        if pfx:
811            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
812        oatom += [names,]
813        Oatoms.append(oatom)
814    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
815    if DATData['Natoms'] == 4:  #torsion
816        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
817        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
818        G2fil.G2Print (' NB: Atom sequence determined by selection order')
819        return      # done with torsion
820    elif DATData['Natoms'] == 3:  #angle
821        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
822        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
823        G2fil.G2Print (' NB: Atom sequence determined by selection order')
824    else:   #2 atoms - distance
825        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
826        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
827
828def BestPlane(PlaneData):
829    'Needs a doc string'
830
831    def ShowBanner(name):
832        G2fil.G2Print (80*'*')
833        G2fil.G2Print ('   Best plane result for phase '+name)
834        G2fil.G2Print (80*'*','\n')
835
836    ShowBanner(PlaneData['Name'])
837
838    Cell = PlaneData['Cell']
839    Amat,Bmat = G2lat.cell2AB(Cell[:6])
840    Atoms = PlaneData['Atoms']
841    sumXYZ = np.zeros(3)
842    XYZ = []
843    Natoms = len(Atoms)
844    for atom in Atoms:
845        xyz = np.array(atom[3:6])
846        XYZ.append(xyz)
847        sumXYZ += xyz
848    sumXYZ /= Natoms
849    XYZ = np.array(XYZ)-sumXYZ
850    XYZ = np.inner(Amat,XYZ).T
851    Zmat = np.zeros((3,3))
852    for i,xyz in enumerate(XYZ):
853        Zmat += np.outer(xyz.T,xyz)
854    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
855    Evec,Emat = nl.eig(Zmat)
856    Evec = np.sqrt(Evec)/(Natoms-3)
857    Order = np.argsort(Evec)
858    XYZ = np.inner(XYZ,Emat.T).T
859    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
860    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
861    G2fil.G2Print (' Name         X         Y         Z')
862    for i,xyz in enumerate(XYZ):
863        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
864    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
865
866def main():
867    'Called to run a refinement when this module is executed '
868    starttime = time.time()
869    arg = sys.argv
870    if len(arg) > 1:
871        GPXfile = arg[1]
872        if not ospath.exists(GPXfile):
873            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
874            exit()
875    else:
876        G2fil.G2Print ('ERROR - missing filename')
877        exit()
878    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
879    Refine(GPXfile,None)
880    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
881
882if __name__ == '__main__':
883    GSASIIpath.InvokeDebugOpts()
884    main()
Note: See TracBrowser for help on using the repository browser.