source: trunk/GSASIIstrMain.py @ 2713

Last change on this file since 2713 was 2713, checked in by vondreele, 5 years ago

avoid extraneous holds @ line 405 in GSASIIstrMain.py

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