source: trunk/GSASIIstrMain.py @ 1835

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

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