source: trunk/GSASIIstrMain.py @ 1819

Last change on this file since 1819 was 1819, checked in by toby, 8 years ago

fix constraints for all test cases
  • 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: 2015-04-29 17:22:48 +0000 (Wed, 29 Apr 2015) $
9# $Author: toby $
10# $Revision: 1819 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 1819 2015-04-29 17:22:48Z toby $
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: 1819 $")
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            maxCyc = Controls['max cyc']
73            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,maxcyc=maxCyc,Print=ifPrint,
74                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
75            ncyc = result[2]['num cyc']+1
76            Rvals['lamMax'] = result[2]['lamMax']
77        else:           #'numeric'
78            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
79                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
80            ncyc = 1
81            if len(varyList):
82                ncyc = int(result[2]['nfev']/len(varyList))
83#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
84#        for item in table: print item,table[item]               #useful debug - are things shifting?
85        runtime = time.time()-begin
86        Rvals['converged'] = result[2].get('Converged')
87        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
88        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
89        G2stMth.Values2Dict(parmDict, varyList, result[0])
90        G2mv.Dict2Map(parmDict,varyList)
91        Rvals['Nobs'] = Histograms['Nobs']
92        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
93        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
94        print >>printFile,' Number of function calls:',result[2]['nfev'],   \
95            ' No. of observations: ',Histograms['Nobs'],' No. of parameters: ',len(varyList),   \
96            ' User rejected: ',Histograms['Nrej'],' Sp. gp. extinct: ',Histograms['Next']
97        print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
98        print >>printFile,' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
99        IfOK = True
100        try:
101            covMatrix = result[1]*Rvals['GOF']**2
102            sig = np.sqrt(np.diag(covMatrix))
103            if np.any(np.isnan(sig)):
104                print '*** Least squares aborted - some invalid esds possible ***'
105#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
106#            for item in table: print item,table[item]               #useful debug - are things shifting?
107            break                   #refinement succeeded - finish up!
108        except TypeError,FloatingPointError:          #result[1] is None on singular matrix
109            IfOK = False
110            if not len(varyList):
111                covMatrix = []
112                sig = []
113                break
114            print '**** Refinement failed - singular matrix ****'
115            if 'Hessian' in Controls['deriv type']:
116                num = len(varyList)-1
117                for i,val in enumerate(np.flipud(result[2]['psing'])):
118                    if val:
119                        print 'Removing parameter: ',varyList[num-i]
120                        del(varyList[num-i])                   
121            else:
122                Ipvt = result[2]['ipvt']
123                for i,ipvt in enumerate(Ipvt):
124                    if not np.sum(result[2]['fjac'],axis=1)[i]:
125                        print 'Removing parameter: ',varyList[ipvt-1]
126                        del(varyList[ipvt-1])
127                        break
128    G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
129    return IfOK,Rvals,result,covMatrix,sig
130
131def Refine(GPXfile,dlg):
132    'Global refinement -- refines to minimize against all histograms'
133    import pytexture as ptx
134    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
135   
136    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
137    G2stIO.ShowBanner(printFile)
138    varyList = []
139    parmDict = {}
140    G2mv.InitVars()   
141    Controls = G2stIO.GetControls(GPXfile)
142    G2stIO.ShowControls(Controls,printFile)
143    calcControls = {}
144    calcControls.update(Controls)           
145    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
146    restraintDict = G2stIO.GetRestraints(GPXfile)
147    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
148    if not Phases:
149        print ' *** ERROR - you have no phases to refine! ***'
150        print ' *** Refine aborted ***'
151        return False,'No phases'
152    if not Histograms:
153        print ' *** ERROR - you have no data to refine with! ***'
154        print ' *** Refine aborted ***'
155        return False,'No data'
156    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
157    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
158    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
159    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
160    calcControls['atomIndx'] = atomIndx
161    calcControls['Natoms'] = Natoms
162    calcControls['FFtables'] = FFtables
163    calcControls['BLtables'] = BLtables
164    calcControls['maxSSwave'] = maxSSwave
165    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
166    calcControls.update(controlDict)
167    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
168    calcControls.update(controlDict)
169    varyList = rbVary+phaseVary+hapVary+histVary
170    parmDict.update(rbDict)
171    parmDict.update(phaseDict)
172    parmDict.update(hapDict)
173    parmDict.update(histDict)
174    G2stIO.GetFprime(calcControls,Histograms)
175    # do constraint processing
176    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
177    try:
178        groups,parmlist = G2mv.GroupConstraints(constrDict)
179        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict)
180        #print G2mv.VarRemapShow(varyList)
181        #print 'DependentVars',G2mv.GetDependentVars()
182        #print 'IndependentVars',G2mv.GetIndependentVars()
183    except:
184        print ' *** ERROR - your constraints are internally inconsistent ***'
185        #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
186        #print 'Errors',errmsg
187        #if warnmsg: print 'Warnings',warnmsg
188        return False,' Constraint error'
189#    print G2mv.VarRemapShow(varyList)
190   
191    ifPrint = True
192    print >>printFile,'\n Refinement results:'
193    print >>printFile,135*'-'
194    try:
195        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
196            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
197        sigDict = dict(zip(varyList,sig))
198        newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
199        newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
200        covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
201                   'varyListStart':varyListStart,
202                   'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
203                   'newCellDict':newCellDict,'freshCOV':True}
204        # add the uncertainties into the esd dictionary (sigDict)
205        sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
206        G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
207        G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
208        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
209        G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
210        G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
211        G2stIO.SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
212        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData)
213        printFile.close()
214        print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
215        print ' ***** Refinement successful *****'
216    except G2obj.G2Exception,Msg:
217        printFile.close()
218        return False,Msg.msg
219   
220#for testing purposes!!!
221    if DEBUG:
222#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
223        import cPickle
224        fl = open('testDeriv.dat','wb')
225        cPickle.dump(result[0],fl,1)
226        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
227        cPickle.dump([G2mv.dependentParmList,G2mv.arrayList,G2mv.invarrayList,
228            G2mv.indParmList,G2mv.invarrayList],fl,1)
229        cPickle.dump(parmDict,fl,1)
230        cPickle.dump(varyList,fl,1)
231        cPickle.dump(calcControls,fl,1)
232        cPickle.dump(pawleyLookup,fl,1)
233        fl.close()
234
235    if dlg:
236        return True,Rvals['Rwp']
237
238def SeqRefine(GPXfile,dlg):
239    '''Perform a sequential refinement -- cycles through all selected histgrams,
240    one at a time
241    '''
242    import pytexture as ptx
243    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
244   
245    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
246    print 'Starting Sequential Refinement'
247    G2stIO.ShowBanner(printFile)
248    Controls = G2stIO.GetControls(GPXfile)
249    G2stIO.ShowControls(Controls,printFile,SeqRef=True)           
250    restraintDict = G2stIO.GetRestraints(GPXfile)
251    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
252    if not Phases:
253        print ' *** ERROR - you have no phases to refine! ***'
254        print ' *** Refine aborted ***'
255        return False,'No phases'
256    if not Histograms:
257        print ' *** ERROR - you have no data to refine with! ***'
258        print ' *** Refine aborted ***'
259        return False,'No data'
260    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
261    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
262    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
263    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile)
264    for item in phaseVary:
265        if '::A0' in item:
266            print '**** WARNING - lattice parameters should not be refined in a sequential refinement ****'
267            print '****           instead use the Dij parameters for each powder histogram            ****'
268            return False,'Lattice parameter refinement error - see console message'
269        if '::C(' in item:
270            print '**** WARNING - phase texture parameters should not be refined in a sequential refinement ****'
271            print '****           instead use the C(L,N) parameters for each powder histogram               ****'
272            return False,'Phase texture refinement error - see console message'
273    if 'Seq Data' in Controls:
274        histNames = Controls['Seq Data']
275    else:
276        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
277    if 'Reverse Seq' in Controls:
278        if Controls['Reverse Seq']:
279            histNames.reverse()
280    SeqResult = {'histNames':histNames}
281    makeBack = True
282    Histo = {}
283    NewparmDict = {}
284    for ihst,histogram in enumerate(histNames):
285        print('Refining with '+str(histogram))
286        ifPrint = False
287        if dlg:
288            dlg.SetTitle('Residual for histogram '+str(ihst))
289        calcControls = {}
290        calcControls['atomIndx'] = atomIndx
291        calcControls['Natoms'] = Natoms
292        calcControls['FFtables'] = FFtables
293        calcControls['BLtables'] = BLtables
294        calcControls['maxSSwave'] = maxSSwave
295        Histo = {histogram:Histograms[histogram],}
296        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
297        calcControls.update(controlDict)
298        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
299        calcControls.update(controlDict)
300        varyList = rbVary+phaseVary+hapVary+histVary
301        if not ihst:
302            # save the initial vary list, but without histogram numbers on parameters
303            saveVaryList = varyList[:]
304            for i,item in enumerate(saveVaryList):
305                items = item.split(':')
306                if items[1]:
307                    items[1] = ''
308                item = ':'.join(items)
309                saveVaryList[i] = item
310            SeqResult['varyList'] = saveVaryList
311        origvaryList = varyList[:]
312        parmDict = {}
313        parmDict.update(phaseDict)
314        parmDict.update(hapDict)
315        parmDict.update(histDict)
316        if Controls['Copy2Next']:
317            parmDict.update(NewparmDict)
318        G2stIO.GetFprime(calcControls,Histo)
319        # do constraint processing
320        #reload(G2mv) # debug
321        G2mv.InitVars()   
322        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
323        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
324        try:
325            groups,parmlist = G2mv.GroupConstraints(constrDict)
326            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst)
327            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
328        except:
329            print ' *** ERROR - your constraints are internally inconsistent ***'
330            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
331            #print 'Errors',errmsg
332            #if warnmsg: print 'Warnings',warnmsg
333            return False,' Constraint error'
334        #print G2mv.VarRemapShow(varyList)
335        if not ihst:
336            # first histogram to refine against
337            firstVaryList = []
338            for item in varyList:
339                items = item.split(':')
340                if items[1]:
341                    items[1] = ''
342                item = ':'.join(items)
343                firstVaryList.append(item)
344            newVaryList = firstVaryList
345        else:
346            newVaryList = []
347            for item in varyList:
348                items = item.split(':')
349                if items[1]:
350                    items[1] = ''
351                item = ':'.join(items)
352                newVaryList.append(item)
353        if newVaryList != firstVaryList and Controls['Copy2Next']:
354            # variable lists are expected to match between sequential refinements when Copy2Next is on
355            print '**** ERROR - variable list for this histogram does not match previous'
356            print '     Copy of variables is not possible'
357            print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
358            combined = list(set(firstVaryList+newVaryList))
359            c = [var for var in combined if var not in newVaryList]
360            p = [var for var in combined if var not in firstVaryList]
361            line = 'Variables in previous but not in current: '
362            if c:
363                for var in c:
364                    if len(line) > 100:
365                        print line
366                        line = '    '
367                    line += var + ', '
368            else:
369                line += 'none'
370            print line
371            print '\nPrevious refinement has',len(firstVaryList),'variables'
372            line = 'Variables in current but not in previous: '
373            if p:
374                for var in p:
375                    if len(line) > 100:
376                        print line
377                        line = '    '
378                    line += var + ', '
379            else:
380                line += 'none'
381            print line
382            return False,line
383       
384        ifPrint = False
385        print >>printFile,'\n Refinement results for histogram: v'+histogram
386        print >>printFile,135*'-'
387        try:
388            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
389                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
390   
391            print '  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
392                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'])
393            # add the uncertainties into the esd dictionary (sigDict)
394            sigDict = dict(zip(varyList,sig))
395            # the uncertainties for dependent constrained parms into the esd dict
396            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
397   
398            # a dict with values & esds for dependent (constrained) parameters
399            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart
400                           if i not in varyList}
401            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
402            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
403            histRefData = {
404                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
405                'varyListStart':varyListStart,
406                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
407                'newCellDict':newCellDict,'depParmDict':depParmDict,
408                'constraintInfo':constraintInfo,
409                'parmDict':parmDict}
410            SeqResult[histogram] = histRefData
411            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
412    #        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
413            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint,printFile)
414            G2stIO.SetHistogramData(parmDict,sigDict,Histo,ifPrint,printFile)
415            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
416            makeBack = False
417            NewparmDict = {}
418            # make dict of varied parameters in current histogram, renamed to
419            # next histogram, for use in next refinement.
420            if Controls['Copy2Next'] and ihst < len(histNames)-1:
421                hId = Histo[histogram]['hId'] # current histogram
422                nexthId = Histograms[histNames[ihst+1]]['hId']
423                for parm in set(list(varyList)+list(varyListStart)):
424                    items = parm.split(':')
425                    if len(items) < 3: continue
426                    if str(hId) in items[1]:
427                        items[1] = str(nexthId)
428                        newparm = ':'.join(items)
429                        NewparmDict[newparm] = parmDict[parm]
430        except G2obj.G2Exception,Msg:
431            printFile.close()
432            print ' ***** Refinement aborted *****'
433            return False,Msg.msg
434    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
435    printFile.close()
436    print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
437    print ' ***** Sequential refinement successful *****'
438    return True,'Success'
439
440def RetDistAngle(DisAglCtls,DisAglData):
441    '''Compute and return distances and angles
442
443    :param dict DisAglCtls: contains distance/angle radii usually defined using
444       :func:`GSASIIgrid.DisAglDialog`
445    :param dict DisAglData: contains phase data:
446       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
447       for distance/angle origins and atoms to be used as targets.
448       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
449
450    :returns: AtomLabels,DistArray,AngArray where:
451
452      **AtomLabels** is a dict of atom labels, keys are the atom number
453
454      **DistArray** is a dict keyed by the origin atom number where the value is a list
455      of distance entries. The value for each distance is a list containing:
456     
457        0) the target atom number (int);
458        1) the unit cell offsets added to x,y & z (tuple of int values)
459        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
460        3) an interatomic distance in A (float)
461        4) an uncertainty on the distance in A or 0.0 (float)
462
463      **AngArray** is a dict keyed by the origin (central) atom number where
464      the value is a list of
465      angle entries. The value for each angle entry consists of three values:
466
467        0) a distance item reference for one neighbor (int)
468        1) a distance item reference for a second neighbor (int)
469        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
470
471      The AngArray distance reference items refer directly to the index of the items in the
472      DistArray item for the list of distances for the central atom.
473    '''
474    import numpy.ma as ma
475   
476    SGData = DisAglData['SGData']
477    Cell = DisAglData['Cell']
478   
479    Amat,Bmat = G2lat.cell2AB(Cell[:6])
480    covData = {}
481    if 'covData' in DisAglData:   
482        covData = DisAglData['covData']
483        covMatrix = covData['covMatrix']
484        varyList = covData['varyList']
485        pfx = str(DisAglData['pId'])+'::'
486        A = G2lat.cell2A(Cell[:6])
487        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
488        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
489        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
490
491    Factor = DisAglCtls['Factors']
492    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
493    indices = (-1,0,1)
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(Xvcov,pdpx)))
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    import numpy.ma as ma
573    def MyPrint(s):
574        out.write(s+'\n')
575        # print(s,file=out) # use in Python 3
576   
577    def ShowBanner(name):
578        MyPrint(80*'*')
579        MyPrint('   Interatomic Distances and Angles for phase '+name)
580        MyPrint((80*'*')+'\n')
581
582    ShowBanner(DisAglCtls['Name'])
583    SGData = DisAglData['SGData']
584    SGtext,SGtable = G2spc.SGPrint(SGData)
585    for line in SGtext: MyPrint(line)
586    if len(SGtable):
587        for i,item in enumerate(SGtable[::2]):
588            line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
589            MyPrint(line)   
590    else:
591        MyPrint(' ( 1)    %s'%(SGtable[0])) 
592    Cell = DisAglData['Cell']
593   
594    Amat,Bmat = G2lat.cell2AB(Cell[:6])
595    covData = {}
596    if 'covData' in DisAglData:   
597        covData = DisAglData['covData']
598        covMatrix = covData['covMatrix']
599        varyList = covData['varyList']
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    targAtoms = DisAglData['TargAtoms']
617    for Oatom in origAtoms:
618        i = Oatom[0]
619        Dist = DistArray[i]
620        nDist = len(Dist)
621        angles = np.zeros((nDist,nDist))
622        angsig = np.zeros((nDist,nDist))
623        for k,j,tup in AngArray[i]:
624            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
625        line = ''
626        for i,x in enumerate(Oatom[3:6]):
627            line += ('%12.5f'%x).rstrip('0')
628        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
629        MyPrint(80*'*')
630        line = ''
631        for dist in Dist[:-1]:
632            line += '%12s'%(AtomLabels[dist[0]].center(12))
633        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
634        for i,dist in enumerate(Dist):
635            line = ''
636            for j,angle in enumerate(angles[i][0:i]):
637                sig = angsig[i][j]
638                if angle:
639                    if sig:
640                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
641                    else:
642                        val = '%.3f'%(angle)
643                        line += '%12s'%(val.center(12))
644                else:
645                    line += 12*' '
646            if dist[4]:            #sig exists!
647                val = G2mth.ValEsd(dist[3],dist[4])
648            else:
649                val = '%8.4f'%(dist[3])
650            tunit = '[%2d%2d%2d]'% dist[1]
651            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
652
653def DisAglTor(DATData):
654    'Needs a doc string'
655    SGData = DATData['SGData']
656    Cell = DATData['Cell']
657   
658    Amat,Bmat = G2lat.cell2AB(Cell[:6])
659    covData = {}
660    pfx = ''
661    if 'covData' in DATData:   
662        covData = DATData['covData']
663        covMatrix = covData['covMatrix']
664        varyList = covData['varyList']
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
734           
735def main():
736    'Needs a doc string'
737    arg = sys.argv
738    if len(arg) > 1:
739        GPXfile = arg[1]
740        if not ospath.exists(GPXfile):
741            print 'ERROR - ',GPXfile," doesn't exist!"
742            exit()
743        GPXpath = ospath.dirname(arg[1])
744        Refine(GPXfile,None)
745    else:
746        print 'ERROR - missing filename'
747        exit()
748    print "Done"
749         
750if __name__ == '__main__':
751    main()
Note: See TracBrowser for help on using the repository browser.