source: trunk/GSASIIstrMain.py @ 1800

Last change on this file since 1800 was 1800, checked in by vondreele, 7 years ago

show no. space group extinct reflections in HKLF data - rejected from calcs.
calculate a "residual" for texture fitting & show it in ProgressBar?
correct angle derivatives in texture fitting (wrong sign)
trap empty useList in plots
make sure U23 is printed in lst output

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