source: trunk/GSASIIstrMain.py @ 2717

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