source: trunk/GSASIIstrMain.py @ 1240

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

add control for copying parms from one histogram to next in a sequential refinement - needs to account for reverse order & gaps

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