source: trunk/GSASIIstrMain.py @ 2557

Last change on this file since 2557 was 2557, checked in by toby, 5 years ago

Fix sequential refinement with not-ready histograms

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