source: trunk/GSASIIstrMain.py @ 1531

Last change on this file since 1531 was 1531, checked in by vondreele, 9 years ago

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