source: trunk/GSASIIstrMain.py @ 1058

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

implement save of Pawley covariances for MC/SA
File size: 31.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,FloatingPointError:          #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,'freshCOV':True}
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 RetDistAngle(DisAglCtls,DisAglData):
389    '''Compute and return distances and angles
390
391    :param dict DisAglCtls: contains distance/angle radii usually defined using
392       :func:`GSASIIgrid.DisAglDialog`
393    :param dict DisAglData: contains phase data:
394       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
395       for distance/angle origins and atoms to be used as targets.
396       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
397
398    :returns: AtomLabels,DistArray,AngArray where:
399
400      **AtomLabels** is a dict of atom labels, keys are the atom number
401
402      **DistArray** is a dict keyed by the origin atom number where the value is a list
403      of distance entries. The value for each distance is a list containing:
404     
405        0) the target atom number (int);
406        1) the unit cell offsets added to x,y & z (tuple of int values)
407        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
408        3) an interatomic distance in A (float)
409        4) an uncertainty on the distance in A or 0.0 (float)
410
411      **AngArray** is a dict keyed by the origin (central) atom number where
412      the value is a list of
413      angle entries. The value for each angle entry consists of three values:
414
415        0) a distance item reference for one neighbor (int)
416        1) a distance item reference for a second neighbor (int)
417        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
418
419      The AngArray distance reference items refer directly to the index of the items in the
420      DistArray item for the list of distances for the central atom.
421    '''
422    import numpy.ma as ma
423   
424    SGData = DisAglData['SGData']
425    Cell = DisAglData['Cell']
426   
427    Amat,Bmat = G2lat.cell2AB(Cell[:6])
428    covData = {}
429    if 'covData' in DisAglData:   
430        covData = DisAglData['covData']
431        covMatrix = covData['covMatrix']
432        varyList = covData['varyList']
433        pfx = str(DisAglData['pId'])+'::'
434        A = G2lat.cell2A(Cell[:6])
435        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
436        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
437        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
438
439    Factor = DisAglCtls['Factors']
440    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
441    indices = (-1,0,1)
442    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
443    origAtoms = DisAglData['OrigAtoms']
444    targAtoms = DisAglData['TargAtoms']
445    AtomLabels = {}
446    for Oatom in origAtoms:
447        AtomLabels[Oatom[0]] = Oatom[1]
448    for Oatom in targAtoms:
449        AtomLabels[Oatom[0]] = Oatom[1]
450    DistArray = {}
451    AngArray = {}
452    for Oatom in origAtoms:
453        DistArray[Oatom[0]] = []
454        AngArray[Oatom[0]] = []
455        OxyzNames = ''
456        IndBlist = []
457        Dist = []
458        Vect = []
459        VectA = []
460        angles = []
461        for Tatom in targAtoms:
462            Xvcov = []
463            TxyzNames = ''
464            if 'covData' in DisAglData:
465                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
466                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
467                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
468            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
469            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
470            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
471            for Txyz,Top,Tunit in result:
472                Dx = (Txyz-np.array(Oatom[3:6]))+Units
473                dx = np.inner(Amat,Dx)
474                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
475                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
476                if np.any(IndB):
477                    for indb in IndB:
478                        for i in range(len(indb)):
479                            if str(dx.T[indb][i]) not in IndBlist:
480                                IndBlist.append(str(dx.T[indb][i]))
481                                unit = Units[indb][i]
482                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
483                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
484                                sig = 0.0
485                                if len(Xvcov):
486                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
487                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
488                                if (Dist[-1][-2]-AsumR) <= 0.:
489                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
490                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
491                                else:
492                                    Vect.append([0.,0.,0.])
493                                    VectA.append([])
494        for D in Dist:
495            DistArray[Oatom[0]].append(D[1:])
496        Vect = np.array(Vect)
497        angles = np.zeros((len(Vect),len(Vect)))
498        angsig = np.zeros((len(Vect),len(Vect)))
499        for i,veca in enumerate(Vect):
500            if np.any(veca):
501                for j,vecb in enumerate(Vect):
502                    if np.any(vecb):
503                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
504                        if i <= j: continue
505                        AngArray[Oatom[0]].append((i,j,
506                                                   G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
507                                                   ))
508    return AtomLabels,DistArray,AngArray
509
510def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
511    '''Print distances and angles
512
513    :param dict DisAglCtls: contains distance/angle radii usually defined using
514       :func:`GSASIIgrid.DisAglDialog`
515    :param dict DisAglData: contains phase data:
516       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
517       for distance/angle origins and atoms to be used as targets.
518       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
519    :param file out: file object for output. Defaults to sys.stdout.   
520    '''
521    import numpy.ma as ma
522    def MyPrint(s):
523        out.write(s+'\n')
524        # print(s,file=out) # use in Python 3
525   
526    def ShowBanner(name):
527        MyPrint(80*'*')
528        MyPrint('   Interatomic Distances and Angles for phase '+name)
529        MyPrint((80*'*')+'\n')
530
531    ShowBanner(DisAglCtls['Name'])
532    SGData = DisAglData['SGData']
533    SGtext = G2spc.SGPrint(SGData)
534    for line in SGtext: MyPrint(line)
535    Cell = DisAglData['Cell']
536   
537    Amat,Bmat = G2lat.cell2AB(Cell[:6])
538    covData = {}
539    if 'covData' in DisAglData:   
540        covData = DisAglData['covData']
541        covMatrix = covData['covMatrix']
542        varyList = covData['varyList']
543        pfx = str(DisAglData['pId'])+'::'
544        A = G2lat.cell2A(Cell[:6])
545        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
546        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
547        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
548        line = '\n Unit cell:'
549        for name,vals in zip(names,valEsd):
550            line += name+vals 
551        MyPrint(line)
552    else: 
553        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
554            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
555            ('%.3f'%Cell[5])+' volume = '+('%.3f'%Cell[6]))
556
557    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
558    origAtoms = DisAglData['OrigAtoms']
559    targAtoms = DisAglData['TargAtoms']
560    for Oatom in origAtoms:
561        i = Oatom[0]
562        Dist = DistArray[i]
563        nDist = len(Dist)
564        angles = np.zeros((nDist,nDist))
565        angsig = np.zeros((nDist,nDist))
566        for k,j,tup in AngArray[i]:
567            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
568        line = ''
569        for i,x in enumerate(Oatom[3:6]):
570            line += ('%12.5f'%x).rstrip('0')
571        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
572        MyPrint(80*'*')
573        line = ''
574        for dist in Dist[:-1]:
575            line += '%12s'%(AtomLabels[dist[0]].center(12))
576        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
577        for i,dist in enumerate(Dist):
578            line = ''
579            for j,angle in enumerate(angles[i][0:i]):
580                sig = angsig[i][j]
581                if angle:
582                    if sig:
583                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
584                    else:
585                        val = '%.3f'%(angle)
586                        line += '%12s'%(val.center(12))
587                else:
588                    line += 12*' '
589            if dist[4]:            #sig exists!
590                val = G2mth.ValEsd(dist[3],dist[4])
591            else:
592                val = '%8.4f'%(dist[3])
593            tunit = '[%2d%2d%2d]'% dist[1]
594            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
595
596def DisAglTor(DATData):
597    'Needs a doc string'
598    SGData = DATData['SGData']
599    Cell = DATData['Cell']
600   
601    Amat,Bmat = G2lat.cell2AB(Cell[:6])
602    covData = {}
603    pfx = ''
604    if 'covData' in DATData:   
605        covData = DATData['covData']
606        covMatrix = covData['covMatrix']
607        varyList = covData['varyList']
608        pfx = str(DATData['pId'])+'::'
609    Datoms = []
610    Oatoms = []
611    for i,atom in enumerate(DATData['Datoms']):
612        symop = atom[-1].split('+')
613        if len(symop) == 1:
614            symop.append('0,0,0')       
615        symop[0] = int(symop[0])
616        symop[1] = eval(symop[1])
617        atom.append(symop)
618        Datoms.append(atom)
619        oatom = DATData['Oatoms'][i]
620        names = ['','','']
621        if pfx:
622            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
623        oatom += [names,]
624        Oatoms.append(oatom)
625    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
626    if DATData['Natoms'] == 4:  #torsion
627        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
628        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
629        print ' NB: Atom sequence determined by selection order'
630        return      # done with torsion
631    elif DATData['Natoms'] == 3:  #angle
632        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
633        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
634        print ' NB: Atom sequence determined by selection order'
635    else:   #2 atoms - distance
636        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
637        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
638               
639def BestPlane(PlaneData):
640    'Needs a doc string'
641
642    def ShowBanner(name):
643        print 80*'*'
644        print '   Best plane result for phase '+name
645        print 80*'*','\n'
646
647    ShowBanner(PlaneData['Name'])
648
649    Cell = PlaneData['Cell']   
650    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
651    Atoms = PlaneData['Atoms']
652    sumXYZ = np.zeros(3)
653    XYZ = []
654    Natoms = len(Atoms)
655    for atom in Atoms:
656        xyz = np.array(atom[3:6])
657        XYZ.append(xyz)
658        sumXYZ += xyz
659    sumXYZ /= Natoms
660    XYZ = np.array(XYZ)-sumXYZ
661    XYZ = np.inner(Amat,XYZ).T
662    Zmat = np.zeros((3,3))
663    for i,xyz in enumerate(XYZ):
664        Zmat += np.outer(xyz.T,xyz)
665    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
666    Evec,Emat = nl.eig(Zmat)
667    Evec = np.sqrt(Evec)/(Natoms-3)
668    Order = np.argsort(Evec)
669    XYZ = np.inner(XYZ,Emat.T).T
670    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
671    print ' Atoms in Cartesian best plane coordinates:'
672    print ' Name         X         Y         Z'
673    for i,xyz in enumerate(XYZ):
674        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
675    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
676
677           
678def main():
679    'Needs a doc string'
680    arg = sys.argv
681    if len(arg) > 1:
682        GPXfile = arg[1]
683        if not ospath.exists(GPXfile):
684            print 'ERROR - ',GPXfile," doesn't exist!"
685            exit()
686        GPXpath = ospath.dirname(arg[1])
687        Refine(GPXfile,None)
688    else:
689        print 'ERROR - missing filename'
690        exit()
691    print "Done"
692         
693if __name__ == '__main__':
694    main()
Note: See TracBrowser for help on using the repository browser.