source: trunk/GSASIIstrMain.py @ 2721

Last change on this file since 2721 was 2721, checked in by vondreele, 6 years ago

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