source: trunk/GSASIIstrMain.py @ 1782

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

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