source: trunk/GSASIIstrMain.py @ 1241

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

copy next now allows reverse & skipped histograms

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