source: trunk/GSASIIstrMain.py @ 1484

Last change on this file since 1484 was 1474, checked in by vondreele, 11 years ago

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