source: trunk/GSASIIstrMain.py @ 2447

Last change on this file since 2447 was 2447, checked in by toby, 7 years ago

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