source: trunk/GSASIIstrMain.py @ 1247

Last change on this file since 1247 was 1247, checked in by toby, 8 years ago

allow wildcarded-histogram constraint input; implement for seq refinements

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