source: trunk/GSASIIstrMain.py @ 1337

Last change on this file since 1337 was 1335, checked in by toby, 11 years ago

Redo pseudovar computation to handle constraint-generated vars
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 31.0 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2014-05-07 19:38:30 +0000 (Wed, 07 May 2014) $
9# $Author: vondreele $
10# $Revision: 1335 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 1335 2014-05-07 19:38:30Z 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: 1335 $")
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            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
296        except:
297            print ' *** ERROR - your constraints are internally inconsistent ***'
298            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
299            #print 'Errors',errmsg
300            #if warnmsg: print 'Warnings',warnmsg
301            raise Exception(' *** Refine aborted ***')
302        #print G2mv.VarRemapShow(varyList)
303        if not ihst:
304            # first histogram to refine against
305            firstVaryList = []
306            for item in varyList:
307                items = item.split(':')
308                if items[1]:
309                    items[1] = ''
310                item = ':'.join(items)
311                firstVaryList.append(item)
312            newVaryList = firstVaryList
313        else:
314            newVaryList = []
315            for item in varyList:
316                items = item.split(':')
317                if items[1]:
318                    items[1] = ''
319                item = ':'.join(items)
320                newVaryList.append(item)
321        if newVaryList != firstVaryList:
322            # variable lists are expected to match between sequential refinements
323            print '**** ERROR - variable list for this histogram does not match previous'
324            print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
325            combined = list(set(firstVaryList+newVaryList))
326            c = [var for var in combined if var not in newVaryList]
327            p = [var for var in combined if var not in firstVaryList]
328            line = 'Variables in previous but not in current: '
329            if c:
330                for var in c:
331                    if len(line) > 100:
332                        print line
333                        line = '    '
334                    line += var + ', '
335            else:
336                line += 'none'
337            print line
338            print '\nPrevious refinement has',len(firstVaryList),'variables'
339            line = 'Variables in current but not in previous: '
340            if p:
341                for var in p:
342                    if len(line) > 100:
343                        print line
344                        line = '    '
345                    line += var + ', '
346            else:
347                line += 'none'
348            print line
349            raise Exception
350       
351        ifPrint = False
352        print >>printFile,'\n Refinement results for histogram: v'+histogram
353        print >>printFile,135*'-'
354        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
355            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
356
357        print '  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
358            Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'],Rvals['DelChi2'])
359        # add the uncertainties into the esd dictionary (sigDict)
360        sigDict = dict(zip(varyList,sig))
361        # the uncertainties for dependent constrained parms into the esd dict
362        sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
363
364        # a dict with values & esds for dependent (constrained) parameters
365        depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart
366                       if i not in varyList}
367        newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
368        newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
369        histRefData = {
370            'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
371            'varyListStart':varyListStart,
372            'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
373            'newCellDict':newCellDict,'depParmDict':depParmDict,
374            'constraintInfo':constraintInfo,
375            'parmDict':parmDict}
376        SeqResult[histogram] = histRefData
377        G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
378#        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
379        G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint,printFile)
380        G2stIO.SetHistogramData(parmDict,sigDict,Histo,ifPrint,printFile)
381        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
382        makeBack = False
383        NewparmDict = {}
384        # make dict of varied parameters in current histogram, renamed to
385        # next histogram, for use in next refinement.
386        if Controls['Copy2Next'] and ihst < len(histNames)-1:
387            hId = Histo[histogram]['hId'] # current histogram
388            nexthId = Histograms[histNames[ihst+1]]['hId']
389            for parm in set(list(varyList)+list(varyListStart)):
390                items = parm.split(':')
391                if len(items) < 3: continue
392                if str(hId) in items[1]:
393                    items[1] = str(nexthId)
394                    newparm = ':'.join(items)
395                    NewparmDict[newparm] = parmDict[parm]
396    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
397    printFile.close()
398    print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
399    print ' ***** Sequential refinement successful *****'
400
401def RetDistAngle(DisAglCtls,DisAglData):
402    '''Compute and return distances and angles
403
404    :param dict DisAglCtls: contains distance/angle radii usually defined using
405       :func:`GSASIIgrid.DisAglDialog`
406    :param dict DisAglData: contains phase data:
407       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
408       for distance/angle origins and atoms to be used as targets.
409       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
410
411    :returns: AtomLabels,DistArray,AngArray where:
412
413      **AtomLabels** is a dict of atom labels, keys are the atom number
414
415      **DistArray** is a dict keyed by the origin atom number where the value is a list
416      of distance entries. The value for each distance is a list containing:
417     
418        0) the target atom number (int);
419        1) the unit cell offsets added to x,y & z (tuple of int values)
420        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
421        3) an interatomic distance in A (float)
422        4) an uncertainty on the distance in A or 0.0 (float)
423
424      **AngArray** is a dict keyed by the origin (central) atom number where
425      the value is a list of
426      angle entries. The value for each angle entry consists of three values:
427
428        0) a distance item reference for one neighbor (int)
429        1) a distance item reference for a second neighbor (int)
430        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
431
432      The AngArray distance reference items refer directly to the index of the items in the
433      DistArray item for the list of distances for the central atom.
434    '''
435    import numpy.ma as ma
436   
437    SGData = DisAglData['SGData']
438    Cell = DisAglData['Cell']
439   
440    Amat,Bmat = G2lat.cell2AB(Cell[:6])
441    covData = {}
442    if 'covData' in DisAglData:   
443        covData = DisAglData['covData']
444        covMatrix = covData['covMatrix']
445        varyList = covData['varyList']
446        pfx = str(DisAglData['pId'])+'::'
447        A = G2lat.cell2A(Cell[:6])
448        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
449        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
450        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
451
452    Factor = DisAglCtls['Factors']
453    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
454    indices = (-1,0,1)
455    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
456    origAtoms = DisAglData['OrigAtoms']
457    targAtoms = DisAglData['TargAtoms']
458    AtomLabels = {}
459    for Oatom in origAtoms:
460        AtomLabels[Oatom[0]] = Oatom[1]
461    for Oatom in targAtoms:
462        AtomLabels[Oatom[0]] = Oatom[1]
463    DistArray = {}
464    AngArray = {}
465    for Oatom in origAtoms:
466        DistArray[Oatom[0]] = []
467        AngArray[Oatom[0]] = []
468        OxyzNames = ''
469        IndBlist = []
470        Dist = []
471        Vect = []
472        VectA = []
473        angles = []
474        for Tatom in targAtoms:
475            Xvcov = []
476            TxyzNames = ''
477            if 'covData' in DisAglData:
478                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
479                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
480                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
481            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
482            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
483            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
484            for Txyz,Top,Tunit in result:
485                Dx = (Txyz-np.array(Oatom[3:6]))+Units
486                dx = np.inner(Amat,Dx)
487                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
488                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
489                if np.any(IndB):
490                    for indb in IndB:
491                        for i in range(len(indb)):
492                            if str(dx.T[indb][i]) not in IndBlist:
493                                IndBlist.append(str(dx.T[indb][i]))
494                                unit = Units[indb][i]
495                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
496                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
497                                sig = 0.0
498                                if len(Xvcov):
499                                    sig = np.sqrt(np.inner(pdpx,np.inner(Xvcov,pdpx)))
500                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
501                                if (Dist[-1][-2]-AsumR) <= 0.:
502                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
503                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
504                                else:
505                                    Vect.append([0.,0.,0.])
506                                    VectA.append([])
507        for D in Dist:
508            DistArray[Oatom[0]].append(D[1:])
509        Vect = np.array(Vect)
510        angles = np.zeros((len(Vect),len(Vect)))
511        angsig = np.zeros((len(Vect),len(Vect)))
512        for i,veca in enumerate(Vect):
513            if np.any(veca):
514                for j,vecb in enumerate(Vect):
515                    if np.any(vecb):
516                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
517                        if i <= j: continue
518                        AngArray[Oatom[0]].append((i,j,
519                                                   G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
520                                                   ))
521    return AtomLabels,DistArray,AngArray
522
523def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
524    '''Print distances and angles
525
526    :param dict DisAglCtls: contains distance/angle radii usually defined using
527       :func:`GSASIIgrid.DisAglDialog`
528    :param dict DisAglData: contains phase data:
529       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
530       for distance/angle origins and atoms to be used as targets.
531       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
532    :param file out: file object for output. Defaults to sys.stdout.   
533    '''
534    import numpy.ma as ma
535    def MyPrint(s):
536        out.write(s+'\n')
537        # print(s,file=out) # use in Python 3
538   
539    def ShowBanner(name):
540        MyPrint(80*'*')
541        MyPrint('   Interatomic Distances and Angles for phase '+name)
542        MyPrint((80*'*')+'\n')
543
544    ShowBanner(DisAglCtls['Name'])
545    SGData = DisAglData['SGData']
546    SGtext = G2spc.SGPrint(SGData)
547    for line in SGtext: MyPrint(line)
548    Cell = DisAglData['Cell']
549   
550    Amat,Bmat = G2lat.cell2AB(Cell[:6])
551    covData = {}
552    if 'covData' in DisAglData:   
553        covData = DisAglData['covData']
554        covMatrix = covData['covMatrix']
555        varyList = covData['varyList']
556        pfx = str(DisAglData['pId'])+'::'
557        A = G2lat.cell2A(Cell[:6])
558        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
559        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
560        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
561        line = '\n Unit cell:'
562        for name,vals in zip(names,valEsd):
563            line += name+vals 
564        MyPrint(line)
565    else: 
566        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
567            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
568            ('%.3f'%Cell[5])+' volume = '+('%.3f'%Cell[6]))
569
570    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
571    origAtoms = DisAglData['OrigAtoms']
572    targAtoms = DisAglData['TargAtoms']
573    for Oatom in origAtoms:
574        i = Oatom[0]
575        Dist = DistArray[i]
576        nDist = len(Dist)
577        angles = np.zeros((nDist,nDist))
578        angsig = np.zeros((nDist,nDist))
579        for k,j,tup in AngArray[i]:
580            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
581        line = ''
582        for i,x in enumerate(Oatom[3:6]):
583            line += ('%12.5f'%x).rstrip('0')
584        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
585        MyPrint(80*'*')
586        line = ''
587        for dist in Dist[:-1]:
588            line += '%12s'%(AtomLabels[dist[0]].center(12))
589        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
590        for i,dist in enumerate(Dist):
591            line = ''
592            for j,angle in enumerate(angles[i][0:i]):
593                sig = angsig[i][j]
594                if angle:
595                    if sig:
596                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
597                    else:
598                        val = '%.3f'%(angle)
599                        line += '%12s'%(val.center(12))
600                else:
601                    line += 12*' '
602            if dist[4]:            #sig exists!
603                val = G2mth.ValEsd(dist[3],dist[4])
604            else:
605                val = '%8.4f'%(dist[3])
606            tunit = '[%2d%2d%2d]'% dist[1]
607            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
608
609def DisAglTor(DATData):
610    'Needs a doc string'
611    SGData = DATData['SGData']
612    Cell = DATData['Cell']
613   
614    Amat,Bmat = G2lat.cell2AB(Cell[:6])
615    covData = {}
616    pfx = ''
617    if 'covData' in DATData:   
618        covData = DATData['covData']
619        covMatrix = covData['covMatrix']
620        varyList = covData['varyList']
621        pfx = str(DATData['pId'])+'::'
622    Datoms = []
623    Oatoms = []
624    for i,atom in enumerate(DATData['Datoms']):
625        symop = atom[-1].split('+')
626        if len(symop) == 1:
627            symop.append('0,0,0')       
628        symop[0] = int(symop[0])
629        symop[1] = eval(symop[1])
630        atom.append(symop)
631        Datoms.append(atom)
632        oatom = DATData['Oatoms'][i]
633        names = ['','','']
634        if pfx:
635            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
636        oatom += [names,]
637        Oatoms.append(oatom)
638    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
639    if DATData['Natoms'] == 4:  #torsion
640        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
641        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
642        print ' NB: Atom sequence determined by selection order'
643        return      # done with torsion
644    elif DATData['Natoms'] == 3:  #angle
645        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
646        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
647        print ' NB: Atom sequence determined by selection order'
648    else:   #2 atoms - distance
649        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
650        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
651               
652def BestPlane(PlaneData):
653    'Needs a doc string'
654
655    def ShowBanner(name):
656        print 80*'*'
657        print '   Best plane result for phase '+name
658        print 80*'*','\n'
659
660    ShowBanner(PlaneData['Name'])
661
662    Cell = PlaneData['Cell']   
663    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
664    Atoms = PlaneData['Atoms']
665    sumXYZ = np.zeros(3)
666    XYZ = []
667    Natoms = len(Atoms)
668    for atom in Atoms:
669        xyz = np.array(atom[3:6])
670        XYZ.append(xyz)
671        sumXYZ += xyz
672    sumXYZ /= Natoms
673    XYZ = np.array(XYZ)-sumXYZ
674    XYZ = np.inner(Amat,XYZ).T
675    Zmat = np.zeros((3,3))
676    for i,xyz in enumerate(XYZ):
677        Zmat += np.outer(xyz.T,xyz)
678    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
679    Evec,Emat = nl.eig(Zmat)
680    Evec = np.sqrt(Evec)/(Natoms-3)
681    Order = np.argsort(Evec)
682    XYZ = np.inner(XYZ,Emat.T).T
683    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
684    print ' Atoms in Cartesian best plane coordinates:'
685    print ' Name         X         Y         Z'
686    for i,xyz in enumerate(XYZ):
687        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
688    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
689
690           
691def main():
692    'Needs a doc string'
693    arg = sys.argv
694    if len(arg) > 1:
695        GPXfile = arg[1]
696        if not ospath.exists(GPXfile):
697            print 'ERROR - ',GPXfile," doesn't exist!"
698            exit()
699        GPXpath = ospath.dirname(arg[1])
700        Refine(GPXfile,None)
701    else:
702        print 'ERROR - missing filename'
703        exit()
704    print "Done"
705         
706if __name__ == '__main__':
707    main()
Note: See TracBrowser for help on using the repository browser.