source: trunk/GSASIIstrMain.py @ 1299

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

Set conf flags consistently for .py files

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