source: trunk/GSASIIstrMain.py @ 1633

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

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