source: trunk/GSASIIstrMain.py @ 1282

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

sequential refinement updates

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