source: trunk/GSASIIstrMain.py @ 2546

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

cleanup of all GSASII main routines - remove unused variables, dead code, etc.

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