source: trunk/GSASIIstrMain.py @ 930

Last change on this file since 930 was 930, checked in by vondreele, 10 years ago

another missing link
File size: 27.3 KB
Line 
1# -*- coding: utf-8 -*-
2#GSASIIstrMain - structure computation main routine
3########### SVN repository information ###################
4# $Date: 2013-05-17 12:13:15 -0500 (Fri, 17 May 2013) $
5# $Author: vondreele $
6# $Revision: 920 $
7# $URL: https://subversion.xor.aps.anl.gov/pyGSAS/trunk/GSASIIstrMain.py $
8# $Id: GSASIIstrMain.py 920 2013-05-17 17:13:15Z vondreele $
9########### SVN repository information ###################
10import sys
11import os
12import os.path as ospath
13import time
14import math
15import copy
16import random
17import cPickle
18import numpy as np
19import numpy.ma as ma
20import numpy.linalg as nl
21import scipy.optimize as so
22import GSASIIpath
23GSASIIpath.SetVersionNumber("$Revision: 920 $")
24import GSASIIlattice as G2lat
25import GSASIIspc as G2spc
26import GSASIImapvars as G2mv
27import GSASIImath as G2mth
28import GSASIIstrIO as G2stIO
29import GSASIIstrMath as G2stMth
30
31sind = lambda x: np.sin(x*np.pi/180.)
32cosd = lambda x: np.cos(x*np.pi/180.)
33tand = lambda x: np.tan(x*np.pi/180.)
34asind = lambda x: 180.*np.arcsin(x)/np.pi
35acosd = lambda x: 180.*np.arccos(x)/np.pi
36atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
37   
38ateln2 = 8.0*math.log(2.0)
39DEBUG = False
40
41
42def Refine(GPXfile,dlg):
43    import pytexture as ptx
44    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
45   
46    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
47    G2stIO.ShowBanner(printFile)
48    varyList = []
49    parmDict = {}
50    G2mv.InitVars()   
51    Controls = G2stIO.GetControls(GPXfile)
52    G2stIO.ShowControls(Controls,printFile)
53    calcControls = {}
54    calcControls.update(Controls)           
55    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
56    restraintDict = G2stIO.GetRestraints(GPXfile)
57    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
58    if not Phases:
59        print ' *** ERROR - you have no phases! ***'
60        print ' *** Refine aborted ***'
61        raise Exception
62    if not Histograms:
63        print ' *** ERROR - you have no data to refine with! ***'
64        print ' *** Refine aborted ***'
65        raise Exception       
66    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
67    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
68    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
69    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
70    calcControls['atomIndx'] = atomIndx
71    calcControls['Natoms'] = Natoms
72    calcControls['FFtables'] = FFtables
73    calcControls['BLtables'] = BLtables
74    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
75    calcControls.update(controlDict)
76    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
77    calcControls.update(controlDict)
78    varyList = rbVary+phaseVary+hapVary+histVary
79    parmDict.update(rbDict)
80    parmDict.update(phaseDict)
81    parmDict.update(hapDict)
82    parmDict.update(histDict)
83    G2stIO.GetFprime(calcControls,Histograms)
84    # do constraint processing
85    try:
86        groups,parmlist = G2mv.GroupConstraints(constrDict)
87        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
88    except:
89        print ' *** ERROR - your constraints are internally inconsistent ***'
90        errmsg, warnmsg = G2stIO.CheckConstraints(varylist,constrDict,fixedList)
91        print 'Errors',errmsg
92        if warnmsg: print 'Warnings',warnmsg
93        raise Exception(' *** Refine aborted ***')
94    # # check to see which generated parameters are fully varied
95    # msg = G2mv.SetVaryFlags(varyList)
96    # if msg:
97    #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
98    #     print msg
99    #     raise Exception(' *** Refine aborted ***')
100    #print G2mv.VarRemapShow(varyList)
101    G2mv.Map2Dict(parmDict,varyList)
102    Rvals = {}
103    while True:
104        begin = time.time()
105        values =  np.array(G2stMth.Dict2Values(parmDict, varyList))
106        Ftol = Controls['min dM/M']
107        Factor = Controls['shift factor']
108        maxCyc = Controls['max cyc']
109        if 'Jacobian' in Controls['deriv type']:           
110            result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
111                ftol=Ftol,col_deriv=True,factor=Factor,
112                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
113            ncyc = int(result[2]['nfev']/2)
114        elif 'Hessian' in Controls['deriv type']:
115            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,maxcyc=maxCyc,
116                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
117            ncyc = result[2]['num cyc']+1
118            Rvals['lamMax'] = result[2]['lamMax']
119        else:           #'numeric'
120            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
121                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
122            ncyc = int(result[2]['nfev']/len(varyList))
123#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
124#        for item in table: print item,table[item]               #useful debug - are things shifting?
125        runtime = time.time()-begin
126        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
127        G2stMth.Values2Dict(parmDict, varyList, result[0])
128        G2mv.Dict2Map(parmDict,varyList)
129        Rvals['Nobs'] = Histograms['Nobs']
130        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
131        Rvals['GOF'] = Rvals['chisq']/(Histograms['Nobs']-len(varyList))
132        print >>printFile,'\n Refinement results:'
133        print >>printFile,135*'-'
134        print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
135        print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
136        print >>printFile,' wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
137        try:
138            covMatrix = result[1]*Rvals['GOF']
139            sig = np.sqrt(np.diag(covMatrix))
140            if np.any(np.isnan(sig)):
141                print '*** Least squares aborted - some invalid esds possible ***'
142#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
143#            for item in table: print item,table[item]               #useful debug - are things shifting?
144            break                   #refinement succeeded - finish up!
145        except TypeError:          #result[1] is None on singular matrix
146            print '**** Refinement failed - singular matrix ****'
147            if 'Hessian' in Controls['deriv type']:
148                num = len(varyList)-1
149                for i,val in enumerate(np.flipud(result[2]['psing'])):
150                    if val:
151                        print 'Removing parameter: ',varyList[num-i]
152                        del(varyList[num-i])                   
153            else:
154                Ipvt = result[2]['ipvt']
155                for i,ipvt in enumerate(Ipvt):
156                    if not np.sum(result[2]['fjac'],axis=1)[i]:
157                        print 'Removing parameter: ',varyList[ipvt-1]
158                        del(varyList[ipvt-1])
159                        break
160
161#    print 'dependentParmList: ',G2mv.dependentParmList
162#    print 'arrayList: ',G2mv.arrayList
163#    print 'invarrayList: ',G2mv.invarrayList
164#    print 'indParmList: ',G2mv.indParmList
165#    print 'fixedDict: ',G2mv.fixedDict
166#    print 'test1'
167    G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
168#    print 'test2'
169    sigDict = dict(zip(varyList,sig))
170    newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
171    newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
172    covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
173        'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
174    # add the uncertainties into the esd dictionary (sigDict)
175    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
176    G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
177    G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
178    G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
179    G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
180    G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
181    G2stIO.SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
182    G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData)
183    printFile.close()
184    print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
185    print ' ***** Refinement successful *****'
186   
187#for testing purposes!!!
188    if DEBUG:
189#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
190        import cPickle
191        fl = open('testDeriv.dat','wb')
192        cPickle.dump(result[0],fl,1)
193        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
194        cPickle.dump(parmDict,fl,1)
195        cPickle.dump(varyList,fl,1)
196        cPickle.dump(calcControls,fl,1)
197        cPickle.dump(pawleyLookup,fl,1)
198        fl.close()
199
200    if dlg:
201        return Rvals['Rwp']
202
203def SeqRefine(GPXfile,dlg):
204    import pytexture as ptx
205    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
206   
207    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
208    print ' Sequential Refinement'
209    G2stIO.ShowBanner(printFile)
210    G2mv.InitVars()   
211    Controls = G2stIO.GetControls(GPXfile)
212    G2stIO.ShowControls(Controls,printFile)           
213    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
214    restraintDict = G2stIO.GetRestraints(GPXfile)
215    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
216    if not Phases:
217        print ' *** ERROR - you have no histograms to refine! ***'
218        print ' *** Refine aborted ***'
219        raise Exception
220    if not Histograms:
221        print ' *** ERROR - you have no data to refine with! ***'
222        print ' *** Refine aborted ***'
223        raise Exception
224    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
225    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
226    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
227    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile)
228    for item in phaseVary:
229        if '::A0' in item:
230            print '**** WARNING - lattice parameters should not be refined in a sequential refinement ****'
231            print '****           instead use the Dij parameters for each powder histogram            ****'
232    if 'Seq Data' in Controls:
233        histNames = Controls['Seq Data']
234    else:
235        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
236    if 'Reverse Seq' in Controls:
237        if Controls['Reverse Seq']:
238            histNames.reverse()
239    SeqResult = {'histNames':histNames}
240    makeBack = True
241    for ihst,histogram in enumerate(histNames):
242        ifPrint = False
243        if dlg:
244            dlg.SetTitle('Residual for histogram '+str(ihst))
245        calcControls = {}
246        calcControls['atomIndx'] = atomIndx
247        calcControls['Natoms'] = Natoms
248        calcControls['FFtables'] = FFtables
249        calcControls['BLtables'] = BLtables
250        varyList = []
251        parmDict = {}
252        Histo = {histogram:Histograms[histogram],}
253        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,False)
254        calcControls.update(controlDict)
255        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
256        calcControls.update(controlDict)
257        varyList = rbVary+phaseVary+hapVary+histVary
258        if not ihst:
259            saveVaryList = varyList[:]
260            for i,item in enumerate(saveVaryList):
261                items = item.split(':')
262                if items[1]:
263                    items[1] = ''
264                item = ':'.join(items)
265                saveVaryList[i] = item
266            SeqResult['varyList'] = saveVaryList
267        else:
268            newVaryList = varyList[:]
269            for i,item in enumerate(newVaryList):
270                items = item.split(':')
271                if items[1]:
272                    items[1] = ''
273                item = ':'.join(items)
274                newVaryList[i] = item
275            if newVaryList != SeqResult['varyList']:
276                print newVaryList
277                print SeqResult['varyList']
278                print '**** ERROR - variable list for this histogram does not match previous'
279                raise Exception
280        parmDict.update(phaseDict)
281        parmDict.update(hapDict)
282        parmDict.update(histDict)
283        G2stIO.GetFprime(calcControls,Histo)
284        # do constraint processing
285        try:
286            groups,parmlist = G2mv.GroupConstraints(constrDict)
287            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList)
288        except:
289            print ' *** ERROR - your constraints are internally inconsistent ***'
290            raise Exception(' *** Refine aborted ***')
291        # check to see which generated parameters are fully varied
292        # msg = G2mv.SetVaryFlags(varyList)
293        # if msg:
294        #     print ' *** ERROR - you have not set the refine flags for constraints consistently! ***'
295        #     print msg
296        #     raise Exception(' *** Refine aborted ***')
297        #print G2mv.VarRemapShow(varyList)
298        G2mv.Map2Dict(parmDict,varyList)
299        Rvals = {}
300        while True:
301            begin = time.time()
302            values =  np.array(G2stMth.Dict2Values(parmDict,varyList))
303            Ftol = Controls['min dM/M']
304            Factor = Controls['shift factor']
305            maxCyc = Controls['max cyc']
306
307            if 'Jacobian' in Controls['deriv type']:           
308                result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
309                    ftol=Ftol,col_deriv=True,factor=Factor,
310                    args=([Histo,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
311                ncyc = int(result[2]['nfev']/2)
312            elif 'Hessian' in Controls['deriv type']:
313                result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,maxcyc=maxCyc,
314                    args=([Histo,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
315                ncyc = result[2]['num cyc']+1                           
316            else:           #'numeric'
317                result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
318                    args=([Histo,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
319                ncyc = int(result[2]['nfev']/len(varyList))
320
321            runtime = time.time()-begin
322            Rvals['chisq'] = np.sum(result[2]['fvec']**2)
323            G2stMth.Values2Dict(parmDict, varyList, result[0])
324            G2mv.Dict2Map(parmDict,varyList)
325           
326            Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histo['sumwYo'])*100.      #to %
327            Rvals['GOF'] = Rvals['Rwp']/(Histo['Nobs']-len(varyList))
328            Rvals['Nobs'] = Histo['Nobs']
329            print >>printFile,'\n Refinement results for histogram: v'+histogram
330            print >>printFile,135*'-'
331            print >>printFile,' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList)
332            print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
333            print >>printFile,' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
334            try:
335                covMatrix = result[1]*Rvals['GOF']
336                sig = np.sqrt(np.diag(covMatrix))
337                if np.any(np.isnan(sig)):
338                    print '*** Least squares aborted - some invalid esds possible ***'
339                    ifPrint = True
340                break                   #refinement succeeded - finish up!
341            except TypeError:          #result[1] is None on singular matrix
342                print '**** Refinement failed - singular matrix ****'
343                if 'Hessian' in Controls['deriv type']:
344                    num = len(varyList)-1
345                    for i,val in enumerate(np.flipud(result[2]['psing'])):
346                        if val:
347                            print 'Removing parameter: ',varyList[num-i]
348                            del(varyList[num-i])                   
349                else:
350                    Ipvt = result[2]['ipvt']
351                    for i,ipvt in enumerate(Ipvt):
352                        if not np.sum(result[2]['fjac'],axis=1)[i]:
353                            print 'Removing parameter: ',varyList[ipvt-1]
354                            del(varyList[ipvt-1])
355                            break
356   
357        G2stMth.GetFobsSq(Histo,Phases,parmDict,calcControls)
358        sigDict = dict(zip(varyList,sig))
359        newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
360        newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
361        covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
362            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict}
363        # add the uncertainties into the esd dictionary (sigDict)
364        G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
365#??        SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
366        G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint,printFile)
367        G2stIO.SetHistogramData(parmDict,sigDict,Histo,ifPrint,printFile)
368        SeqResult[histogram] = covData
369        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,covData,makeBack)
370        makeBack = False
371    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
372    printFile.close()
373    print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
374    print ' ***** Sequential refinement successful *****'
375
376def DistAngle(DisAglCtls,DisAglData):
377    import numpy.ma as ma
378   
379    def ShowBanner(name):
380        print 80*'*'
381        print '   Interatomic Distances and Angles for phase '+name
382        print 80*'*','\n'
383
384    ShowBanner(DisAglCtls['Name'])
385    SGData = DisAglData['SGData']
386    SGtext = G2spc.SGPrint(SGData)
387    for line in SGtext: print line
388    Cell = DisAglData['Cell']
389   
390    Amat,Bmat = G2lat.cell2AB(Cell[:6])
391    covData = {}
392    if 'covData' in DisAglData:   
393        covData = DisAglData['covData']
394        covMatrix = covData['covMatrix']
395        varyList = covData['varyList']
396        pfx = str(DisAglData['pId'])+'::'
397        A = G2lat.cell2A(Cell[:6])
398        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
399        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
400        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
401        line = '\n Unit cell:'
402        for name,vals in zip(names,valEsd):
403            line += name+vals 
404        print line
405    else: 
406        print '\n Unit cell: a = ','%.5f'%(Cell[0]),' b = ','%.5f'%(Cell[1]),' c = ','%.5f'%(Cell[2]), \
407            ' alpha = ','%.3f'%(Cell[3]),' beta = ','%.3f'%(Cell[4]),' gamma = ', \
408            '%.3f'%(Cell[5]),' volume = ','%.3f'%(Cell[6])
409    Factor = DisAglCtls['Factors']
410    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
411    indices = (-1,0,1)
412    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
413    origAtoms = DisAglData['OrigAtoms']
414    targAtoms = DisAglData['TargAtoms']
415    for Oatom in origAtoms:
416        OxyzNames = ''
417        IndBlist = []
418        Dist = []
419        Vect = []
420        VectA = []
421        angles = []
422        for Tatom in targAtoms:
423            Xvcov = []
424            TxyzNames = ''
425            if 'covData' in DisAglData:
426                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
427                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
428                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
429            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
430            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
431            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
432            for Txyz,Top,Tunit in result:
433                Dx = (Txyz-np.array(Oatom[3:6]))+Units
434                dx = np.inner(Amat,Dx)
435                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
436                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
437                if np.any(IndB):
438                    for indb in IndB:
439                        for i in range(len(indb)):
440                            if str(dx.T[indb][i]) not in IndBlist:
441                                IndBlist.append(str(dx.T[indb][i]))
442                                unit = Units[indb][i]
443                                tunit = '[%2d%2d%2d]'%(unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
444                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
445                                sig = 0.0
446                                if len(Xvcov):
447                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
448                                Dist.append([Oatom[1],Tatom[1],tunit,Top,ma.getdata(dist[indb])[i],sig])
449                                if (Dist[-1][-2]-AsumR) <= 0.:
450                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
451                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
452                                else:
453                                    Vect.append([0.,0.,0.])
454                                    VectA.append([])
455        Vect = np.array(Vect)
456        angles = np.zeros((len(Vect),len(Vect)))
457        angsig = np.zeros((len(Vect),len(Vect)))
458        for i,veca in enumerate(Vect):
459            if np.any(veca):
460                for j,vecb in enumerate(Vect):
461                    if np.any(vecb):
462                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
463        line = ''
464        for i,x in enumerate(Oatom[3:6]):
465            if len(Xvcov):
466                line += '%12s'%(G2mth.ValEsd(x,np.sqrt(Xvcov[i][i]),True))
467            else:
468                line += '%12.5f'%(x)
469        print '\n Distances & angles for ',Oatom[1],' at ',line
470        print 80*'*'
471        line = ''
472        for dist in Dist[:-1]:
473            line += '%12s'%(dist[1].center(12))
474        print '  To       cell +(sym. op.)      dist.  ',line
475        for i,dist in enumerate(Dist):
476            line = ''
477            for j,angle in enumerate(angles[i][0:i]):
478                sig = angsig[i][j]
479                if angle:
480                    if sig:
481                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
482                    else:
483                        val = '%.3f'%(angle)
484                        line += '%12s'%(val.center(12))
485                else:
486                    line += 12*' '
487            if dist[5]:            #sig exists!
488                val = G2mth.ValEsd(dist[4],dist[5])
489            else:
490                val = '%8.4f'%(dist[4])
491            print '  %8s%10s+(%4d) %12s'%(dist[1].ljust(8),dist[2].ljust(10),dist[3],val.center(12)),line
492
493def DisAglTor(DATData):
494    SGData = DATData['SGData']
495    Cell = DATData['Cell']
496   
497    Amat,Bmat = G2lat.cell2AB(Cell[:6])
498    covData = {}
499    pfx = ''
500    if 'covData' in DATData:   
501        covData = DATData['covData']
502        covMatrix = covData['covMatrix']
503        varyList = covData['varyList']
504        pfx = str(DATData['pId'])+'::'
505    Datoms = []
506    Oatoms = []
507    for i,atom in enumerate(DATData['Datoms']):
508        symop = atom[-1].split('+')
509        if len(symop) == 1:
510            symop.append('0,0,0')       
511        symop[0] = int(symop[0])
512        symop[1] = eval(symop[1])
513        atom.append(symop)
514        Datoms.append(atom)
515        oatom = DATData['Oatoms'][i]
516        names = ['','','']
517        if pfx:
518            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
519        oatom += [names,]
520        Oatoms.append(oatom)
521    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
522    if DATData['Natoms'] == 4:  #torsion
523        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
524        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
525        print ' NB: Atom sequence determined by selection order'
526        return      # done with torsion
527    elif DATData['Natoms'] == 3:  #angle
528        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
529        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
530        print ' NB: Atom sequence determined by selection order'
531    else:   #2 atoms - distance
532        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
533        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
534               
535def BestPlane(PlaneData):
536
537    def ShowBanner(name):
538        print 80*'*'
539        print '   Best plane result for phase '+name
540        print 80*'*','\n'
541
542    ShowBanner(PlaneData['Name'])
543
544    Cell = PlaneData['Cell']   
545    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
546    Atoms = PlaneData['Atoms']
547    sumXYZ = np.zeros(3)
548    XYZ = []
549    Natoms = len(Atoms)
550    for atom in Atoms:
551        xyz = np.array(atom[3:6])
552        XYZ.append(xyz)
553        sumXYZ += xyz
554    sumXYZ /= Natoms
555    XYZ = np.array(XYZ)-sumXYZ
556    XYZ = np.inner(Amat,XYZ).T
557    Zmat = np.zeros((3,3))
558    for i,xyz in enumerate(XYZ):
559        Zmat += np.outer(xyz.T,xyz)
560    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
561    Evec,Emat = nl.eig(Zmat)
562    Evec = np.sqrt(Evec)/(Natoms-3)
563    Order = np.argsort(Evec)
564    XYZ = np.inner(XYZ,Emat.T).T
565    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
566    print ' Atoms in Cartesian best plane coordinates:'
567    print ' Name         X         Y         Z'
568    for i,xyz in enumerate(XYZ):
569        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
570    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
571
572           
573def main():
574    arg = sys.argv
575    if len(arg) > 1:
576        GPXfile = arg[1]
577        if not ospath.exists(GPXfile):
578            print 'ERROR - ',GPXfile," doesn't exist!"
579            exit()
580        GPXpath = ospath.dirname(arg[1])
581        Refine(GPXfile,None)
582    else:
583        print 'ERROR - missing filename'
584        exit()
585    print "Done"
586         
587if __name__ == '__main__':
588    main()
Note: See TracBrowser for help on using the repository browser.