source: trunk/GSASIIstrMain.py @ 2500

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

protect a dlg in DoIndexPeaks?
fix a lost G2frame.LimitsTable? - need to go back to tree instead
add depVarList to testDeriv file
more work on mag derivs.
modify testDeriv to show depVar derivatives

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