source: trunk/GSASIIstrMain.py @ 1813

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

revise GSAS-II exception handling - traps user aborts & I hope the bad metric tensor error

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