source: trunk/GSASIIstrMain.py @ 2175

Last change on this file since 2175 was 2126, checked in by vondreele, 9 years ago

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