source: trunk/GSASIIstrMain.py @ 2847

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

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