source: trunk/GSASIIstrMain.py @ 1004

Last change on this file since 1004 was 975, checked in by vondreele, 12 years ago

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