source: trunk/GSASIIstrMain.py @ 1772

Last change on this file since 1772 was 1772, checked in by vondreele, 8 years ago

fix a copy flags error in DData
add Compute average for columns in Seq results.
fix missing GeneralMass? error
fix covariance motion error - now reports values correctly
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 32.1 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2015-03-31 17:44:58 +0000 (Tue, 31 Mar 2015) $
9# $Author: vondreele $
10# $Revision: 1772 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 1772 2015-03-31 17:44:58Z 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: 1772 $")
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 = False
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'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList)
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, reduced chi**2 = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2)
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,FloatingPointError:          #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! ***'
148        print ' *** Refine aborted ***'
149        raise Exception
150    if not Histograms:
151        print ' *** ERROR - you have no data to refine with! ***'
152        print ' *** Refine aborted ***'
153        raise Exception       
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,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
158    calcControls['atomIndx'] = atomIndx
159    calcControls['Natoms'] = Natoms
160    calcControls['FFtables'] = FFtables
161    calcControls['BLtables'] = BLtables
162    calcControls['maxSSwave'] = maxSSwave
163    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
164    calcControls.update(controlDict)
165    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
166    calcControls.update(controlDict)
167    varyList = rbVary+phaseVary+hapVary+histVary
168    parmDict.update(rbDict)
169    parmDict.update(phaseDict)
170    parmDict.update(hapDict)
171    parmDict.update(histDict)
172    G2stIO.GetFprime(calcControls,Histograms)
173    # do constraint processing
174    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
175    try:
176        groups,parmlist = G2mv.GroupConstraints(constrDict)
177        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict)
178    except:
179        print ' *** ERROR - your constraints are internally inconsistent ***'
180        #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
181        #print 'Errors',errmsg
182        #if warnmsg: print 'Warnings',warnmsg
183        raise Exception(' *** Refine aborted ***')
184#    print G2mv.VarRemapShow(varyList)
185   
186    ifPrint = True
187    print >>printFile,'\n Refinement results:'
188    print >>printFile,135*'-'
189    IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
190        rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
191    sigDict = dict(zip(varyList,sig))
192    newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
193    newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
194    covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
195               'varyListStart':varyListStart,
196               'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
197               'newCellDict':newCellDict,'freshCOV':True}
198    # add the uncertainties into the esd dictionary (sigDict)
199    sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
200    G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
201    G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
202    G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
203    G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
204    G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,pFile=printFile)
205    G2stIO.SetHistogramData(parmDict,sigDict,Histograms,pFile=printFile)
206    G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData)
207    printFile.close()
208    print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
209    print ' ***** Refinement successful *****'
210   
211#for testing purposes!!!
212    if DEBUG:
213#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
214        import cPickle
215        fl = open('testDeriv.dat','wb')
216        cPickle.dump(result[0],fl,1)
217        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
218        cPickle.dump([G2mv.dependentParmList,G2mv.arrayList,G2mv.invarrayList,
219            G2mv.indParmList,G2mv.invarrayList],fl,1)
220        cPickle.dump(parmDict,fl,1)
221        cPickle.dump(varyList,fl,1)
222        cPickle.dump(calcControls,fl,1)
223        cPickle.dump(pawleyLookup,fl,1)
224        fl.close()
225
226    if dlg:
227        return Rvals['Rwp']
228
229def SeqRefine(GPXfile,dlg):
230    '''Perform a sequential refinement -- cycles through all selected histgrams,
231    one at a time
232    '''
233    import pytexture as ptx
234    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
235   
236    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
237    print 'Starting Sequential Refinement'
238    G2stIO.ShowBanner(printFile)
239    Controls = G2stIO.GetControls(GPXfile)
240    G2stIO.ShowControls(Controls,printFile,SeqRef=True)           
241    restraintDict = G2stIO.GetRestraints(GPXfile)
242    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
243    if not Phases:
244        print ' *** ERROR - you have no phases to refine! ***'
245        print ' *** Refine aborted ***'
246        raise Exception
247    if not Histograms:
248        print ' *** ERROR - you have no data to refine with! ***'
249        print ' *** Refine aborted ***'
250        raise Exception
251    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
252    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
253    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
254    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,maxSSwave = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile)
255    for item in phaseVary:
256        if '::A0' in item:
257            print '**** WARNING - lattice parameters should not be refined in a sequential refinement ****'
258            print '****           instead use the Dij parameters for each powder histogram            ****'
259    if 'Seq Data' in Controls:
260        histNames = Controls['Seq Data']
261    else:
262        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
263    if 'Reverse Seq' in Controls:
264        if Controls['Reverse Seq']:
265            histNames.reverse()
266    SeqResult = {'histNames':histNames}
267    makeBack = True
268    Histo = {}
269    NewparmDict = {}
270    for ihst,histogram in enumerate(histNames):
271        print('Refining with '+str(histogram))
272        ifPrint = False
273        if dlg:
274            dlg.SetTitle('Residual for histogram '+str(ihst))
275        calcControls = {}
276        calcControls['atomIndx'] = atomIndx
277        calcControls['Natoms'] = Natoms
278        calcControls['FFtables'] = FFtables
279        calcControls['BLtables'] = BLtables
280        calcControls['maxSSwave'] = maxSSwave
281        Histo = {histogram:Histograms[histogram],}
282        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
283        calcControls.update(controlDict)
284        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
285        calcControls.update(controlDict)
286        varyList = rbVary+phaseVary+hapVary+histVary
287        if not ihst:
288            # save the initial vary list, but without histogram numbers on parameters
289            saveVaryList = varyList[:]
290            for i,item in enumerate(saveVaryList):
291                items = item.split(':')
292                if items[1]:
293                    items[1] = ''
294                item = ':'.join(items)
295                saveVaryList[i] = item
296            SeqResult['varyList'] = saveVaryList
297        origvaryList = varyList[:]
298        parmDict = {}
299        parmDict.update(phaseDict)
300        parmDict.update(hapDict)
301        parmDict.update(histDict)
302        if Controls['Copy2Next']:
303            parmDict.update(NewparmDict)
304        G2stIO.GetFprime(calcControls,Histo)
305        # do constraint processing
306        #reload(G2mv) # debug
307        G2mv.InitVars()   
308        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
309        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
310        try:
311            groups,parmlist = G2mv.GroupConstraints(constrDict)
312            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst)
313            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
314        except:
315            print ' *** ERROR - your constraints are internally inconsistent ***'
316            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
317            #print 'Errors',errmsg
318            #if warnmsg: print 'Warnings',warnmsg
319            raise Exception(' *** Refine aborted ***')
320        #print G2mv.VarRemapShow(varyList)
321        if not ihst:
322            # first histogram to refine against
323            firstVaryList = []
324            for item in varyList:
325                items = item.split(':')
326                if items[1]:
327                    items[1] = ''
328                item = ':'.join(items)
329                firstVaryList.append(item)
330            newVaryList = firstVaryList
331        else:
332            newVaryList = []
333            for item in varyList:
334                items = item.split(':')
335                if items[1]:
336                    items[1] = ''
337                item = ':'.join(items)
338                newVaryList.append(item)
339        if newVaryList != firstVaryList and Controls['Copy2Next']:
340            # variable lists are expected to match between sequential refinements when Copy2Next is on
341            print '**** ERROR - variable list for this histogram does not match previous'
342            print '     Copy of variables is not possible'
343            print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
344            combined = list(set(firstVaryList+newVaryList))
345            c = [var for var in combined if var not in newVaryList]
346            p = [var for var in combined if var not in firstVaryList]
347            line = 'Variables in previous but not in current: '
348            if c:
349                for var in c:
350                    if len(line) > 100:
351                        print line
352                        line = '    '
353                    line += var + ', '
354            else:
355                line += 'none'
356            print line
357            print '\nPrevious refinement has',len(firstVaryList),'variables'
358            line = 'Variables in current but not in previous: '
359            if p:
360                for var in p:
361                    if len(line) > 100:
362                        print line
363                        line = '    '
364                    line += var + ', '
365            else:
366                line += 'none'
367            print line
368            raise Exception
369       
370        ifPrint = False
371        print >>printFile,'\n Refinement results for histogram: v'+histogram
372        print >>printFile,135*'-'
373        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
374            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
375
376        print '  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
377            Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'])
378        # add the uncertainties into the esd dictionary (sigDict)
379        sigDict = dict(zip(varyList,sig))
380        # the uncertainties for dependent constrained parms into the esd dict
381        sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
382
383        # a dict with values & esds for dependent (constrained) parameters
384        depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart
385                       if i not in varyList}
386        newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
387        newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
388        histRefData = {
389            'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
390            'varyListStart':varyListStart,
391            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
392            'newCellDict':newCellDict,'depParmDict':depParmDict,
393            'constraintInfo':constraintInfo,
394            'parmDict':parmDict}
395        SeqResult[histogram] = histRefData
396        G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
397#        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
398        G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint,printFile)
399        G2stIO.SetHistogramData(parmDict,sigDict,Histo,ifPrint,printFile)
400        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
401        makeBack = False
402        NewparmDict = {}
403        # make dict of varied parameters in current histogram, renamed to
404        # next histogram, for use in next refinement.
405        if Controls['Copy2Next'] and ihst < len(histNames)-1:
406            hId = Histo[histogram]['hId'] # current histogram
407            nexthId = Histograms[histNames[ihst+1]]['hId']
408            for parm in set(list(varyList)+list(varyListStart)):
409                items = parm.split(':')
410                if len(items) < 3: continue
411                if str(hId) in items[1]:
412                    items[1] = str(nexthId)
413                    newparm = ':'.join(items)
414                    NewparmDict[newparm] = parmDict[parm]
415    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
416    printFile.close()
417    print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
418    print ' ***** Sequential refinement successful *****'
419
420def RetDistAngle(DisAglCtls,DisAglData):
421    '''Compute and return distances and angles
422
423    :param dict DisAglCtls: contains distance/angle radii usually defined using
424       :func:`GSASIIgrid.DisAglDialog`
425    :param dict DisAglData: contains phase data:
426       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
427       for distance/angle origins and atoms to be used as targets.
428       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
429
430    :returns: AtomLabels,DistArray,AngArray where:
431
432      **AtomLabels** is a dict of atom labels, keys are the atom number
433
434      **DistArray** is a dict keyed by the origin atom number where the value is a list
435      of distance entries. The value for each distance is a list containing:
436     
437        0) the target atom number (int);
438        1) the unit cell offsets added to x,y & z (tuple of int values)
439        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
440        3) an interatomic distance in A (float)
441        4) an uncertainty on the distance in A or 0.0 (float)
442
443      **AngArray** is a dict keyed by the origin (central) atom number where
444      the value is a list of
445      angle entries. The value for each angle entry consists of three values:
446
447        0) a distance item reference for one neighbor (int)
448        1) a distance item reference for a second neighbor (int)
449        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
450
451      The AngArray distance reference items refer directly to the index of the items in the
452      DistArray item for the list of distances for the central atom.
453    '''
454    import numpy.ma as ma
455   
456    SGData = DisAglData['SGData']
457    Cell = DisAglData['Cell']
458   
459    Amat,Bmat = G2lat.cell2AB(Cell[:6])
460    covData = {}
461    if 'covData' in DisAglData:   
462        covData = DisAglData['covData']
463        covMatrix = covData['covMatrix']
464        varyList = covData['varyList']
465        pfx = str(DisAglData['pId'])+'::'
466        A = G2lat.cell2A(Cell[:6])
467        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
468        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
469        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
470
471    Factor = DisAglCtls['Factors']
472    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
473    indices = (-1,0,1)
474    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
475    origAtoms = DisAglData['OrigAtoms']
476    targAtoms = DisAglData['TargAtoms']
477    AtomLabels = {}
478    for Oatom in origAtoms:
479        AtomLabels[Oatom[0]] = Oatom[1]
480    for Oatom in targAtoms:
481        AtomLabels[Oatom[0]] = Oatom[1]
482    DistArray = {}
483    AngArray = {}
484    for Oatom in origAtoms:
485        DistArray[Oatom[0]] = []
486        AngArray[Oatom[0]] = []
487        OxyzNames = ''
488        IndBlist = []
489        Dist = []
490        Vect = []
491        VectA = []
492        angles = []
493        for Tatom in targAtoms:
494            Xvcov = []
495            TxyzNames = ''
496            if 'covData' in DisAglData:
497                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
498                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
499                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
500            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
501            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
502            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
503            for Txyz,Top,Tunit in result:
504                Dx = (Txyz-np.array(Oatom[3:6]))+Units
505                dx = np.inner(Amat,Dx)
506                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
507                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
508                if np.any(IndB):
509                    for indb in IndB:
510                        for i in range(len(indb)):
511                            if str(dx.T[indb][i]) not in IndBlist:
512                                IndBlist.append(str(dx.T[indb][i]))
513                                unit = Units[indb][i]
514                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
515                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
516                                sig = 0.0
517                                if len(Xvcov):
518                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
519                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
520                                if (Dist[-1][-2]-AsumR) <= 0.:
521                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
522                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
523                                else:
524                                    Vect.append([0.,0.,0.])
525                                    VectA.append([])
526        for D in Dist:
527            DistArray[Oatom[0]].append(D[1:])
528        Vect = np.array(Vect)
529        angles = np.zeros((len(Vect),len(Vect)))
530        angsig = np.zeros((len(Vect),len(Vect)))
531        for i,veca in enumerate(Vect):
532            if np.any(veca):
533                for j,vecb in enumerate(Vect):
534                    if np.any(vecb):
535                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
536                        if i <= j: continue
537                        AngArray[Oatom[0]].append((i,j,
538                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
539    return AtomLabels,DistArray,AngArray
540
541def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
542    '''Print distances and angles
543
544    :param dict DisAglCtls: contains distance/angle radii usually defined using
545       :func:`GSASIIgrid.DisAglDialog`
546    :param dict DisAglData: contains phase data:
547       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
548       for distance/angle origins and atoms to be used as targets.
549       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
550    :param file out: file object for output. Defaults to sys.stdout.   
551    '''
552    import numpy.ma as ma
553    def MyPrint(s):
554        out.write(s+'\n')
555        # print(s,file=out) # use in Python 3
556   
557    def ShowBanner(name):
558        MyPrint(80*'*')
559        MyPrint('   Interatomic Distances and Angles for phase '+name)
560        MyPrint((80*'*')+'\n')
561
562    ShowBanner(DisAglCtls['Name'])
563    SGData = DisAglData['SGData']
564    SGtext,SGtable = G2spc.SGPrint(SGData)
565    for line in SGtext: MyPrint(line)
566    if len(SGtable):
567        for i,item in enumerate(SGtable[::2]):
568            line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
569            MyPrint(line)   
570    else:
571        MyPrint(' ( 1)    %s'%(SGtable[0])) 
572    Cell = DisAglData['Cell']
573   
574    Amat,Bmat = G2lat.cell2AB(Cell[:6])
575    covData = {}
576    if 'covData' in DisAglData:   
577        covData = DisAglData['covData']
578        covMatrix = covData['covMatrix']
579        varyList = covData['varyList']
580        pfx = str(DisAglData['pId'])+'::'
581        A = G2lat.cell2A(Cell[:6])
582        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
583        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
584        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
585        line = '\n Unit cell:'
586        for name,vals in zip(names,valEsd):
587            line += name+vals 
588        MyPrint(line)
589    else: 
590        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
591            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
592            ('%.3f'%Cell[5])+' volume = '+('%.3f'%Cell[6]))
593
594    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
595    origAtoms = DisAglData['OrigAtoms']
596    targAtoms = DisAglData['TargAtoms']
597    for Oatom in origAtoms:
598        i = Oatom[0]
599        Dist = DistArray[i]
600        nDist = len(Dist)
601        angles = np.zeros((nDist,nDist))
602        angsig = np.zeros((nDist,nDist))
603        for k,j,tup in AngArray[i]:
604            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
605        line = ''
606        for i,x in enumerate(Oatom[3:6]):
607            line += ('%12.5f'%x).rstrip('0')
608        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
609        MyPrint(80*'*')
610        line = ''
611        for dist in Dist[:-1]:
612            line += '%12s'%(AtomLabels[dist[0]].center(12))
613        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
614        for i,dist in enumerate(Dist):
615            line = ''
616            for j,angle in enumerate(angles[i][0:i]):
617                sig = angsig[i][j]
618                if angle:
619                    if sig:
620                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
621                    else:
622                        val = '%.3f'%(angle)
623                        line += '%12s'%(val.center(12))
624                else:
625                    line += 12*' '
626            if dist[4]:            #sig exists!
627                val = G2mth.ValEsd(dist[3],dist[4])
628            else:
629                val = '%8.4f'%(dist[3])
630            tunit = '[%2d%2d%2d]'% dist[1]
631            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
632
633def DisAglTor(DATData):
634    'Needs a doc string'
635    SGData = DATData['SGData']
636    Cell = DATData['Cell']
637   
638    Amat,Bmat = G2lat.cell2AB(Cell[:6])
639    covData = {}
640    pfx = ''
641    if 'covData' in DATData:   
642        covData = DATData['covData']
643        covMatrix = covData['covMatrix']
644        varyList = covData['varyList']
645        pfx = str(DATData['pId'])+'::'
646    Datoms = []
647    Oatoms = []
648    for i,atom in enumerate(DATData['Datoms']):
649        symop = atom[-1].split('+')
650        if len(symop) == 1:
651            symop.append('0,0,0')       
652        symop[0] = int(symop[0])
653        symop[1] = eval(symop[1])
654        atom.append(symop)
655        Datoms.append(atom)
656        oatom = DATData['Oatoms'][i]
657        names = ['','','']
658        if pfx:
659            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
660        oatom += [names,]
661        Oatoms.append(oatom)
662    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
663    if DATData['Natoms'] == 4:  #torsion
664        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
665        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
666        print ' NB: Atom sequence determined by selection order'
667        return      # done with torsion
668    elif DATData['Natoms'] == 3:  #angle
669        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
670        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
671        print ' NB: Atom sequence determined by selection order'
672    else:   #2 atoms - distance
673        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
674        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
675               
676def BestPlane(PlaneData):
677    'Needs a doc string'
678
679    def ShowBanner(name):
680        print 80*'*'
681        print '   Best plane result for phase '+name
682        print 80*'*','\n'
683
684    ShowBanner(PlaneData['Name'])
685
686    Cell = PlaneData['Cell']   
687    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
688    Atoms = PlaneData['Atoms']
689    sumXYZ = np.zeros(3)
690    XYZ = []
691    Natoms = len(Atoms)
692    for atom in Atoms:
693        xyz = np.array(atom[3:6])
694        XYZ.append(xyz)
695        sumXYZ += xyz
696    sumXYZ /= Natoms
697    XYZ = np.array(XYZ)-sumXYZ
698    XYZ = np.inner(Amat,XYZ).T
699    Zmat = np.zeros((3,3))
700    for i,xyz in enumerate(XYZ):
701        Zmat += np.outer(xyz.T,xyz)
702    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
703    Evec,Emat = nl.eig(Zmat)
704    Evec = np.sqrt(Evec)/(Natoms-3)
705    Order = np.argsort(Evec)
706    XYZ = np.inner(XYZ,Emat.T).T
707    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
708    print ' Atoms in Cartesian best plane coordinates:'
709    print ' Name         X         Y         Z'
710    for i,xyz in enumerate(XYZ):
711        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
712    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
713
714           
715def main():
716    'Needs a doc string'
717    arg = sys.argv
718    if len(arg) > 1:
719        GPXfile = arg[1]
720        if not ospath.exists(GPXfile):
721            print 'ERROR - ',GPXfile," doesn't exist!"
722            exit()
723        GPXpath = ospath.dirname(arg[1])
724        Refine(GPXfile,None)
725    else:
726        print 'ERROR - missing filename'
727        exit()
728    print "Done"
729         
730if __name__ == '__main__':
731    main()
Note: See TracBrowser for help on using the repository browser.