source: trunk/GSASIIstrMain.py @ 1515

Last change on this file since 1515 was 1493, checked in by vondreele, 11 years ago

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