source: trunk/GSASIIstrMain.py @ 1643

Last change on this file since 1643 was 1643, checked in by vondreele, 7 years ago

fix input of TOF hkl data
fix bond angle calculations
fix 2D & 3D hkl plotting

  • 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-02-10 18:24:35 +0000 (Tue, 10 Feb 2015) $
9# $Author: vondreele $
10# $Revision: 1643 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 1643 2015-02-10 18:24:35Z 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: 1643 $")
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    return AtomLabels,DistArray,AngArray
531
532def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
533    '''Print distances and angles
534
535    :param dict DisAglCtls: contains distance/angle radii usually defined using
536       :func:`GSASIIgrid.DisAglDialog`
537    :param dict DisAglData: contains phase data:
538       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
539       for distance/angle origins and atoms to be used as targets.
540       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
541    :param file out: file object for output. Defaults to sys.stdout.   
542    '''
543    import numpy.ma as ma
544    def MyPrint(s):
545        out.write(s+'\n')
546        # print(s,file=out) # use in Python 3
547   
548    def ShowBanner(name):
549        MyPrint(80*'*')
550        MyPrint('   Interatomic Distances and Angles for phase '+name)
551        MyPrint((80*'*')+'\n')
552
553    ShowBanner(DisAglCtls['Name'])
554    SGData = DisAglData['SGData']
555    SGtext,SGtable = G2spc.SGPrint(SGData)
556    for line in SGtext: MyPrint(line)
557    if len(SGtable):
558        for i,item in enumerate(SGtable[::2]):
559            line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
560            MyPrint(line)   
561    else:
562        MyPrint(' ( 1)    %s'%(SGtable[0])) 
563    Cell = DisAglData['Cell']
564   
565    Amat,Bmat = G2lat.cell2AB(Cell[:6])
566    covData = {}
567    if 'covData' in DisAglData:   
568        covData = DisAglData['covData']
569        covMatrix = covData['covMatrix']
570        varyList = covData['varyList']
571        pfx = str(DisAglData['pId'])+'::'
572        A = G2lat.cell2A(Cell[:6])
573        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
574        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
575        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
576        line = '\n Unit cell:'
577        for name,vals in zip(names,valEsd):
578            line += name+vals 
579        MyPrint(line)
580    else: 
581        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
582            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
583            ('%.3f'%Cell[5])+' volume = '+('%.3f'%Cell[6]))
584
585    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
586    origAtoms = DisAglData['OrigAtoms']
587    targAtoms = DisAglData['TargAtoms']
588    for Oatom in origAtoms:
589        i = Oatom[0]
590        Dist = DistArray[i]
591        nDist = len(Dist)
592        angles = np.zeros((nDist,nDist))
593        angsig = np.zeros((nDist,nDist))
594        for k,j,tup in AngArray[i]:
595            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
596        line = ''
597        for i,x in enumerate(Oatom[3:6]):
598            line += ('%12.5f'%x).rstrip('0')
599        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
600        MyPrint(80*'*')
601        line = ''
602        for dist in Dist[:-1]:
603            line += '%12s'%(AtomLabels[dist[0]].center(12))
604        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
605        for i,dist in enumerate(Dist):
606            line = ''
607            for j,angle in enumerate(angles[i][0:i]):
608                sig = angsig[i][j]
609                if angle:
610                    if sig:
611                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
612                    else:
613                        val = '%.3f'%(angle)
614                        line += '%12s'%(val.center(12))
615                else:
616                    line += 12*' '
617            if dist[4]:            #sig exists!
618                val = G2mth.ValEsd(dist[3],dist[4])
619            else:
620                val = '%8.4f'%(dist[3])
621            tunit = '[%2d%2d%2d]'% dist[1]
622            MyPrint((%8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
623
624def DisAglTor(DATData):
625    'Needs a doc string'
626    SGData = DATData['SGData']
627    Cell = DATData['Cell']
628   
629    Amat,Bmat = G2lat.cell2AB(Cell[:6])
630    covData = {}
631    pfx = ''
632    if 'covData' in DATData:   
633        covData = DATData['covData']
634        covMatrix = covData['covMatrix']
635        varyList = covData['varyList']
636        pfx = str(DATData['pId'])+'::'
637    Datoms = []
638    Oatoms = []
639    for i,atom in enumerate(DATData['Datoms']):
640        symop = atom[-1].split('+')
641        if len(symop) == 1:
642            symop.append('0,0,0')       
643        symop[0] = int(symop[0])
644        symop[1] = eval(symop[1])
645        atom.append(symop)
646        Datoms.append(atom)
647        oatom = DATData['Oatoms'][i]
648        names = ['','','']
649        if pfx:
650            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
651        oatom += [names,]
652        Oatoms.append(oatom)
653    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
654    if DATData['Natoms'] == 4:  #torsion
655        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
656        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
657        print ' NB: Atom sequence determined by selection order'
658        return      # done with torsion
659    elif DATData['Natoms'] == 3:  #angle
660        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
661        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
662        print ' NB: Atom sequence determined by selection order'
663    else:   #2 atoms - distance
664        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
665        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
666               
667def BestPlane(PlaneData):
668    'Needs a doc string'
669
670    def ShowBanner(name):
671        print 80*'*'
672        print '   Best plane result for phase '+name
673        print 80*'*','\n'
674
675    ShowBanner(PlaneData['Name'])
676
677    Cell = PlaneData['Cell']   
678    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
679    Atoms = PlaneData['Atoms']
680    sumXYZ = np.zeros(3)
681    XYZ = []
682    Natoms = len(Atoms)
683    for atom in Atoms:
684        xyz = np.array(atom[3:6])
685        XYZ.append(xyz)
686        sumXYZ += xyz
687    sumXYZ /= Natoms
688    XYZ = np.array(XYZ)-sumXYZ
689    XYZ = np.inner(Amat,XYZ).T
690    Zmat = np.zeros((3,3))
691    for i,xyz in enumerate(XYZ):
692        Zmat += np.outer(xyz.T,xyz)
693    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
694    Evec,Emat = nl.eig(Zmat)
695    Evec = np.sqrt(Evec)/(Natoms-3)
696    Order = np.argsort(Evec)
697    XYZ = np.inner(XYZ,Emat.T).T
698    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
699    print ' Atoms in Cartesian best plane coordinates:'
700    print ' Name         X         Y         Z'
701    for i,xyz in enumerate(XYZ):
702        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
703    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
704
705           
706def main():
707    'Needs a doc string'
708    arg = sys.argv
709    if len(arg) > 1:
710        GPXfile = arg[1]
711        if not ospath.exists(GPXfile):
712            print 'ERROR - ',GPXfile," doesn't exist!"
713            exit()
714        GPXpath = ospath.dirname(arg[1])
715        Refine(GPXfile,None)
716    else:
717        print 'ERROR - missing filename'
718        exit()
719    print "Done"
720         
721if __name__ == '__main__':
722    main()
Note: See TracBrowser for help on using the repository browser.