source: trunk/GSASIIstrMain.py @ 1248

Last change on this file since 1248 was 1248, checked in by vondreele, 8 years ago

add setscale for SASD data
don't square the 1/cos(2-theta) correction to integrated intensities in ImageIntegrate?
scale SASD error bars by Scale
replace ':' with ';' in BkPk? parameter names
fix bug of missing Residuals in LS Refine

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