source: trunk/GSASIIstrMain.py @ 3100

Last change on this file since 3100 was 3100, checked in by vondreele, 4 years ago

implement fix FXU for sequential refinements

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