source: trunk/GSASIIstrMain.py @ 2998

Last change on this file since 2998 was 2926, checked in by toby, 6 years ago

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