source: trunk/GSASIIstrMain.py @ 956

Last change on this file since 956 was 956, checked in by toby, 10 years ago

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