source: trunk/GSASIIstrMain.py @ 1141

Last change on this file since 1141 was 1141, checked in by toby, 9 years ago

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