source: branch/2frame/GSASIIstrMain.py @ 2899

Last change on this file since 2899 was 2899, checked in by toby, 6 years ago

partial reorg
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 34.7 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2017-07-03 21:12:45 +0000 (Mon, 03 Jul 2017) $
9# $Author: toby $
10# $Revision: 2899 $
11# $URL: branch/2frame/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 2899 2017-07-03 21:12:45Z toby $
13########### SVN repository information ###################
14import sys
15import os.path as ospath
16import time
17import math
18import copy
19import numpy as np
20import numpy.linalg as nl
21import scipy.optimize as so
22import GSASIIpath
23GSASIIpath.SetVersionNumber("$Revision: 2899 $")
24import GSASIIlattice as G2lat
25import GSASIIspc as G2spc
26import GSASIImapvars as G2mv
27import GSASIImath as G2mth
28import GSASIIstrIO as G2stIO
29import GSASIIstrMath as G2stMth
30import GSASIIobj as G2obj
31
32sind = lambda x: np.sin(x*np.pi/180.)
33cosd = lambda x: np.cos(x*np.pi/180.)
34tand = lambda x: np.tan(x*np.pi/180.)
35asind = lambda x: 180.*np.arcsin(x)/np.pi
36acosd = lambda x: 180.*np.arccos(x)/np.pi
37atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
38   
39ateln2 = 8.0*math.log(2.0)
40DEBUG = True
41
42def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
43    calcControls,pawleyLookup,ifPrint,printFile,dlg):
44    'Core optimization routines, shared between SeqRefine and Refine'
45#    print 'current',varyList
46#    for item in parmDict: print item,parmDict[item] ######### show dict just before refinement
47    G2mv.Map2Dict(parmDict,varyList)
48    Rvals = {}
49    while True:
50        begin = time.time()
51        values =  np.array(G2stMth.Dict2Values(parmDict, varyList))
52        # test code to compute GOF and save for external repeat
53        #args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)
54        #print '*** before fit chi**2',np.sum(G2stMth.errRefine(values,*args)**2)           
55        #fl = open('beforeFit.cpickle','wb')
56        #import cPickle
57        #cPickle.dump(values,fl,1)
58        #cPickle.dump(args[:-1],fl,1)
59        #fl.close()
60        Ftol = Controls['min dM/M']
61        Xtol = Controls['SVDtol']
62        Factor = Controls['shift factor']
63        if 'Jacobian' in Controls['deriv type']:           
64            result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
65                ftol=Ftol,col_deriv=True,factor=Factor,
66                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
67            ncyc = int(result[2]['nfev']/2)
68        elif 'analytic Hessian' in Controls['deriv type']:
69            Lamda = Controls.get('Marquardt',-3)
70            maxCyc = Controls['max cyc']
71            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
72                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
73            ncyc = result[2]['num cyc']+1
74            Rvals['lamMax'] = result[2]['lamMax']
75            Controls['Marquardt'] = -3  #reset to default
76        elif 'Hessian SVD' in Controls['deriv type']:
77            maxCyc = Controls['max cyc']
78            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
79                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
80            ncyc = result[2]['num cyc']+1
81        else:           #'numeric'
82            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
83                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
84            ncyc = 1
85            if len(varyList):
86                ncyc = int(result[2]['nfev']/len(varyList))
87#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
88#        for item in table: print item,table[item]               #useful debug - are things shifting?
89        runtime = time.time()-begin
90        Rvals['SVD0'] = result[2].get('SVD0',0)
91        Rvals['converged'] = result[2].get('Converged')
92        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
93        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
94        G2stMth.Values2Dict(parmDict, varyList, result[0])
95        G2mv.Dict2Map(parmDict,varyList)
96        Rvals['Nobs'] = Histograms['Nobs']
97        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
98        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
99        print >>printFile,' Number of function calls:',result[2]['nfev'],   \
100            ' No. of observations: ',Histograms['Nobs'],' No. of parameters: ',len(varyList),   \
101            ' User rejected: ',Histograms['Nrej'],' Sp. gp. extinct: ',Histograms['Next']
102        print >>printFile,' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc)
103        print >>printFile,' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF'])
104        IfOK = True
105        try:
106            covMatrix = result[1]*Rvals['GOF']**2
107            sig = np.sqrt(np.diag(covMatrix))
108            if np.any(np.isnan(sig)) or not sig.shape:
109                print '*** Least squares aborted - some invalid esds possible ***'
110#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
111#            for item in table: print item,table[item]               #useful debug - are things shifting?
112            break                   #refinement succeeded - finish up!
113        except TypeError:          #result[1] is None on singular matrix or LinAlgError
114            IfOK = False
115            if not len(varyList):
116                covMatrix = []
117                sig = []
118                break
119            print '**** Refinement failed - singular matrix ****'
120            if 'Hessian' in Controls['deriv type']:
121                num = len(varyList)-1
122                for i,val in enumerate(np.flipud(result[2]['psing'])):
123                    if val:
124                        print 'Removing parameter: ',varyList[num-i]
125                        del(varyList[num-i])                   
126            else:
127                Ipvt = result[2]['ipvt']
128                for i,ipvt in enumerate(Ipvt):
129                    if not np.sum(result[2]['fjac'],axis=1)[i]:
130                        print 'Removing parameter: ',varyList[ipvt-1]
131                        del(varyList[ipvt-1])
132                        break
133    G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
134    return IfOK,Rvals,result,covMatrix,sig
135
136def Refine(GPXfile,dlg=None):
137    'Global refinement -- refines to minimize against all histograms'
138    import pytexture as ptx
139    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
140   
141    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
142    G2stIO.ShowBanner(printFile)
143    varyList = []
144    parmDict = {}
145    G2mv.InitVars()   
146    Controls = G2stIO.GetControls(GPXfile)
147    G2stIO.ShowControls(Controls,printFile)
148    calcControls = {}
149    calcControls.update(Controls)           
150    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
151    restraintDict = G2stIO.GetRestraints(GPXfile)
152    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
153    if not Phases:
154        print ' *** ERROR - you have no phases to refine! ***'
155        print ' *** Refine aborted ***'
156        return False,'No phases'
157    if not Histograms:
158        print ' *** ERROR - you have no data to refine with! ***'
159        print ' *** Refine aborted ***'
160        return False,'No data'
161    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
162    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
163    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
164    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
165        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
166    calcControls['atomIndx'] = atomIndx
167    calcControls['Natoms'] = Natoms
168    calcControls['FFtables'] = FFtables
169    calcControls['BLtables'] = BLtables
170    calcControls['MFtables'] = MFtables
171    calcControls['maxSSwave'] = maxSSwave
172    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
173    TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
174    constrDict += TwConstr
175    fixedList += TwFixed
176    calcControls.update(controlDict)
177    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
178    calcControls.update(controlDict)
179    varyList = rbVary+phaseVary+hapVary+histVary
180    parmDict.update(rbDict)
181    parmDict.update(phaseDict)
182    parmDict.update(hapDict)
183    parmDict.update(histDict)
184    G2stIO.GetFprime(calcControls,Histograms)
185    # do constraint processing
186    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
187    try:
188        groups,parmlist = G2mv.GroupConstraints(constrDict)
189        G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict)
190        #print G2mv.VarRemapShow(varyList)
191        #print 'DependentVars',G2mv.GetDependentVars()
192        #print 'IndependentVars',G2mv.GetIndependentVars()
193    except:
194        print ' *** ERROR - your constraints are internally inconsistent ***'
195        #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
196        #print 'Errors',errmsg
197        #if warnmsg: print 'Warnings',warnmsg
198        return False,' Constraint error'
199#    print G2mv.VarRemapShow(varyList)
200   
201    ifPrint = True
202    print >>printFile,'\n Refinement results:'
203    print >>printFile,135*'-'
204    try:
205        covData = {}
206        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
207            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
208        sigDict = dict(zip(varyList,sig))
209        newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
210        newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
211        covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
212                   'varyListStart':varyListStart,
213                   'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
214                   'newCellDict':newCellDict,'freshCOV':True}
215        # add the uncertainties into the esd dictionary (sigDict)
216        sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
217        G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
218        G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
219        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
220        G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
221        G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls['FFtables'],pFile=printFile)
222        G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls['FFtables'],pFile=printFile)
223        G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData)
224        printFile.close()
225        print ' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
226        print ' ***** Refinement successful *****'
227    except G2obj.G2Exception,Msg:
228        printFile.close()
229        return False,Msg.msg
230   
231#for testing purposes!!!
232    if DEBUG:
233#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
234        import cPickle
235        fl = open('testDeriv.dat','wb')
236        cPickle.dump(result[0],fl,1)
237        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
238        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
239        cPickle.dump(parmDict,fl,1)
240        cPickle.dump(varyList,fl,1)
241        cPickle.dump(calcControls,fl,1)
242        cPickle.dump(pawleyLookup,fl,1)
243        fl.close()
244    if dlg:
245        return True,Rvals
246       
247def phaseCheck(phaseVary,Phases,histogram):
248    '''
249    Removes unused parameters from phase varylist if phase not in histogram
250    '''
251    pIds = []
252    for phase in Phases:
253        if histogram not in Phases[phase]['Histograms']: continue
254        if Phases[phase]['Histograms'][histogram]['Use']:
255            pIds.append(str(Phases[phase]['pId']))
256    return [item for item in phaseVary if item.split(':')[0] in pIds]
257
258def SeqRefine(GPXfile,dlg,PlotFunction=None,G2frame=None):
259    '''Perform a sequential refinement -- cycles through all selected histgrams,
260    one at a time
261    '''
262    import pytexture as ptx
263    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
264   
265    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
266    print 'Starting Sequential Refinement'
267    G2stIO.ShowBanner(printFile)
268    Controls = G2stIO.GetControls(GPXfile)
269    G2stIO.ShowControls(Controls,printFile,SeqRef=True)           
270    restraintDict = G2stIO.GetRestraints(GPXfile)
271    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
272    if not Phases:
273        print ' *** ERROR - you have no phases to refine! ***'
274        print ' *** Refine aborted ***'
275        return False,'No phases'
276    if not Histograms:
277        print ' *** ERROR - you have no data to refine with! ***'
278        print ' *** Refine aborted ***'
279        return False,'No data'
280    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
281    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
282    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
283    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
284        G2stIO.GetPhaseData(Phases,restraintDict,rbIds,False,printFile,seqRef=True)
285    for item in phaseVary:
286        if '::A0' in item:
287            print '**** WARNING - lattice parameters should not be refined in a sequential refinement ****'
288            print '****           instead use the Dij parameters for each powder histogram            ****'
289            return False,'Lattice parameter refinement error - see console message'
290        if '::C(' in item:
291            print '**** WARNING - phase texture parameters should not be refined in a sequential refinement ****'
292            print '****           instead use the C(L,N) parameters for each powder histogram               ****'
293            return False,'Phase texture refinement error - see console message'
294    if 'Seq Data' in Controls:
295        histNames = Controls['Seq Data']
296    else:
297        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
298    if Controls.get('Reverse Seq'):
299        histNames.reverse()
300    SeqResult = G2stIO.GetSeqResult(GPXfile)
301#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
302    makeBack = True
303    Histo = {}
304    NewparmDict = {}
305    for ihst,histogram in enumerate(histNames):
306        print('\nRefining with '+str(histogram))
307        ifPrint = False
308        if dlg:
309            dlg.SetTitle('Residual for histogram '+str(ihst))
310        calcControls = {}
311        calcControls['atomIndx'] = atomIndx
312        calcControls['Natoms'] = Natoms
313        calcControls['FFtables'] = FFtables
314        calcControls['BLtables'] = BLtables
315        calcControls['MFtables'] = MFtables
316        calcControls['maxSSwave'] = maxSSwave
317        if histogram not in Histograms:
318            print("Error: not found!")
319            continue
320        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
321        Histo = {histogram:Histograms[histogram],}
322        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
323        calcControls.update(controlDict)
324        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
325        calcControls.update(controlDict)
326        varyList = rbVary+redphaseVary+hapVary+histVary
327#        if not ihst:
328            # save the initial vary list, but without histogram numbers on parameters
329        saveVaryList = varyList[:]
330        for i,item in enumerate(saveVaryList):
331            items = item.split(':')
332            if items[1]:
333                items[1] = ''
334            item = ':'.join(items)
335            saveVaryList[i] = item
336        if not ihst:
337            SeqResult['varyList'] = saveVaryList
338        else:
339            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
340        parmDict = {}
341        parmDict.update(phaseDict)
342        parmDict.update(hapDict)
343        parmDict.update(histDict)
344        if Controls['Copy2Next']:
345            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
346            for parm in NewparmDict:
347                if parm in parmDict:
348                    parmDict[parm] = NewparmDict[parm]
349        G2stIO.GetFprime(calcControls,Histo)
350        # do constraint processing
351        #reload(G2mv) # debug
352        G2mv.InitVars()   
353        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
354        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
355        try:
356            groups,parmlist = G2mv.GroupConstraints(constrDict)
357            G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,fixedList,parmDict,SeqHist=ihst)
358            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
359        except:
360            print ' *** ERROR - your constraints are internally inconsistent ***'
361            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
362            #print 'Errors',errmsg
363            #if warnmsg: print 'Warnings',warnmsg
364            return False,' Constraint error'
365        #print G2mv.VarRemapShow(varyList)
366        if not ihst:
367            # first histogram to refine against
368            firstVaryList = []
369            for item in varyList:
370                items = item.split(':')
371                if items[1]:
372                    items[1] = ''
373                item = ':'.join(items)
374                firstVaryList.append(item)
375            newVaryList = firstVaryList
376        else:
377            newVaryList = []
378            for item in varyList:
379                items = item.split(':')
380                if items[1]:
381                    items[1] = ''
382                item = ':'.join(items)
383                newVaryList.append(item)
384        if newVaryList != firstVaryList and Controls['Copy2Next']:
385            # variable lists are expected to match between sequential refinements when Copy2Next is on
386            #print '**** ERROR - variable list for this histogram does not match previous'
387            #print '     Copy of variables is not possible'
388            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
389            combined = list(set(firstVaryList+newVaryList))
390            c = [var for var in combined if var not in newVaryList]
391            p = [var for var in combined if var not in firstVaryList]
392            print('*** Variables change ***')
393            for typ,vars in [('Removed',c),('Added',p)]:
394                line = '  '+typ+': '
395                if vars:
396                    for var in vars:
397                        if len(line) > 70:
398                            print(line)
399                            line = '    '
400                        line += var + ', '
401                else:
402                        line += 'none, '
403                print(line[:-2])
404            firstVaryList = newVaryList
405       
406        ifPrint = False
407        print >>printFile,'\n Refinement results for histogram: v'+histogram
408        print >>printFile,135*'-'
409        try:
410            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
411                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifPrint,printFile,dlg)
412            if PlotFunction:
413                PlotFunction(G2frame,Histo[histogram]['Data'],histogram)
414   
415            print '  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f'%(
416                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'])
417            # add the uncertainties into the esd dictionary (sigDict)
418            sigDict = dict(zip(varyList,sig))
419            # the uncertainties for dependent constrained parms into the esd dict
420            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
421   
422            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
423            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
424            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
425            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
426            histRefData = {
427                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
428                'varyListStart':varyListStart,
429                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
430                'newCellDict':newCellDict,'depParmDict':depParmDict,
431                'constraintInfo':constraintInfo,
432                'parmDict':parmDict}
433            SeqResult[histogram] = histRefData
434            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
435    #        G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
436            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
437            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile)
438            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,makeBack)
439            makeBack = False
440            NewparmDict = {}
441            # make dict of varied parameters in current histogram, renamed to
442            # next histogram, for use in next refinement.
443            if Controls['Copy2Next'] and ihst < len(histNames)-1:
444                hId = Histo[histogram]['hId'] # current histogram
445                nexthId = Histograms[histNames[ihst+1]]['hId']
446                for parm in set(list(varyList)+list(varyListStart)):
447                    items = parm.split(':')
448                    if len(items) < 3: continue
449                    if str(hId) in items[1]:
450                        items[1] = str(nexthId)
451                        newparm = ':'.join(items)
452                        NewparmDict[newparm] = parmDict[parm]
453        except G2obj.G2Exception,Msg:
454            printFile.close()
455            print ' ***** Refinement aborted *****'
456            return False,Msg.msg
457    if Controls.get('Reverse Seq'):
458        histNames.reverse()
459    SeqResult['histNames'] = histNames
460    G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
461    printFile.close()
462    print ' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst'
463    print ' ***** Sequential refinement successful *****'
464    return True,'Success'
465
466def RetDistAngle(DisAglCtls,DisAglData):
467    '''Compute and return distances and angles
468
469    :param dict DisAglCtls: contains distance/angle radii usually defined using
470       :func:`GSASIIctrls.DisAglDialog`
471    :param dict DisAglData: contains phase data:
472       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
473       for distance/angle origins and atoms to be used as targets.
474       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
475
476    :returns: AtomLabels,DistArray,AngArray where:
477
478      **AtomLabels** is a dict of atom labels, keys are the atom number
479
480      **DistArray** is a dict keyed by the origin atom number where the value is a list
481      of distance entries. The value for each distance is a list containing:
482     
483        0) the target atom number (int);
484        1) the unit cell offsets added to x,y & z (tuple of int values)
485        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
486        3) an interatomic distance in A (float)
487        4) an uncertainty on the distance in A or 0.0 (float)
488
489      **AngArray** is a dict keyed by the origin (central) atom number where
490      the value is a list of
491      angle entries. The value for each angle entry consists of three values:
492
493        0) a distance item reference for one neighbor (int)
494        1) a distance item reference for a second neighbor (int)
495        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
496
497      The AngArray distance reference items refer directly to the index of the items in the
498      DistArray item for the list of distances for the central atom.
499    '''
500    import numpy.ma as ma
501   
502    SGData = DisAglData['SGData']
503    Cell = DisAglData['Cell']
504   
505    Amat,Bmat = G2lat.cell2AB(Cell[:6])
506    covData = {}
507    if 'covData' in DisAglData:   
508        covData = DisAglData['covData']
509        covMatrix = covData['covMatrix']
510        varyList = covData['varyList']
511        pfx = str(DisAglData['pId'])+'::'
512       
513    Factor = DisAglCtls['Factors']
514    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
515    indices = (-2,-1,0,1,2)
516    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
517    origAtoms = DisAglData['OrigAtoms']
518    targAtoms = DisAglData['TargAtoms']
519    AtomLabels = {}
520    for Oatom in origAtoms:
521        AtomLabels[Oatom[0]] = Oatom[1]
522    for Oatom in targAtoms:
523        AtomLabels[Oatom[0]] = Oatom[1]
524    DistArray = {}
525    AngArray = {}
526    for Oatom in origAtoms:
527        DistArray[Oatom[0]] = []
528        AngArray[Oatom[0]] = []
529        OxyzNames = ''
530        IndBlist = []
531        Dist = []
532        Vect = []
533        VectA = []
534        angles = []
535        for Tatom in targAtoms:
536            Xvcov = []
537            TxyzNames = ''
538            if 'covData' in DisAglData:
539                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
540                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
541                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
542            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
543            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
544            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
545            for Txyz,Top,Tunit in result:
546                Dx = (Txyz-np.array(Oatom[3:6]))+Units
547                dx = np.inner(Amat,Dx)
548                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
549                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
550                if np.any(IndB):
551                    for indb in IndB:
552                        for i in range(len(indb)):
553                            if str(dx.T[indb][i]) not in IndBlist:
554                                IndBlist.append(str(dx.T[indb][i]))
555                                unit = Units[indb][i]
556                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
557                                pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
558                                sig = 0.0
559                                if len(Xvcov):
560                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
561                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
562                                if (Dist[-1][-2]-AsumR) <= 0.:
563                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
564                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
565                                else:
566                                    Vect.append([0.,0.,0.])
567                                    VectA.append([])
568        for D in Dist:
569            DistArray[Oatom[0]].append(D[1:])
570        Vect = np.array(Vect)
571        angles = np.zeros((len(Vect),len(Vect)))
572        angsig = np.zeros((len(Vect),len(Vect)))
573        for i,veca in enumerate(Vect):
574            if np.any(veca):
575                for j,vecb in enumerate(Vect):
576                    if np.any(vecb):
577                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
578                        if i <= j: continue
579                        AngArray[Oatom[0]].append((i,j,
580                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
581    return AtomLabels,DistArray,AngArray
582
583def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
584    '''Print distances and angles
585
586    :param dict DisAglCtls: contains distance/angle radii usually defined using
587       :func:`GSASIIctrls.DisAglDialog`
588    :param dict DisAglData: contains phase data:
589       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
590       for distance/angle origins and atoms to be used as targets.
591       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
592    :param file out: file object for output. Defaults to sys.stdout.   
593    '''
594    def MyPrint(s):
595        out.write(s+'\n')
596        # print(s,file=out) # use in Python 3
597   
598    def ShowBanner(name):
599        MyPrint(80*'*')
600        MyPrint('   Interatomic Distances and Angles for phase '+name)
601        MyPrint((80*'*')+'\n')
602
603    ShowBanner(DisAglCtls['Name'])
604    SGData = DisAglData['SGData']
605    SGtext,SGtable = G2spc.SGPrint(SGData)
606    for line in SGtext: MyPrint(line)
607    if len(SGtable) > 1:
608        for i,item in enumerate(SGtable[::2]):
609            if 2*i+1 == len(SGtable):
610                line = ' %s'%(item.ljust(30))
611            else:
612                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
613            MyPrint(line)   
614    else:
615        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
616    Cell = DisAglData['Cell']
617   
618    Amat,Bmat = G2lat.cell2AB(Cell[:6])
619    covData = {}
620    if 'covData' in DisAglData:   
621        covData = DisAglData['covData']
622        pfx = str(DisAglData['pId'])+'::'
623        A = G2lat.cell2A(Cell[:6])
624        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
625        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
626        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
627        line = '\n Unit cell:'
628        for name,vals in zip(names,valEsd):
629            line += name+vals 
630        MyPrint(line)
631    else: 
632        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
633            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
634            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
635
636    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
637    origAtoms = DisAglData['OrigAtoms']
638    for Oatom in origAtoms:
639        i = Oatom[0]
640        Dist = DistArray[i]
641        nDist = len(Dist)
642        angles = np.zeros((nDist,nDist))
643        angsig = np.zeros((nDist,nDist))
644        for k,j,tup in AngArray[i]:
645            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
646        line = ''
647        for i,x in enumerate(Oatom[3:6]):
648            line += ('%12.5f'%x).rstrip('0')
649        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
650        MyPrint(80*'*')
651        line = ''
652        for dist in Dist[:-1]:
653            line += '%12s'%(AtomLabels[dist[0]].center(12))
654        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
655        for i,dist in enumerate(Dist):
656            line = ''
657            for j,angle in enumerate(angles[i][0:i]):
658                sig = angsig[i][j]
659                if angle:
660                    if sig:
661                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
662                    else:
663                        val = '%.3f'%(angle)
664                        line += '%12s'%(val.center(12))
665                else:
666                    line += 12*' '
667            if dist[4]:            #sig exists!
668                val = G2mth.ValEsd(dist[3],dist[4])
669            else:
670                val = '%8.4f'%(dist[3])
671            tunit = '[%2d%2d%2d]'% dist[1]
672            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
673
674def DisAglTor(DATData):
675    'Needs a doc string'
676    SGData = DATData['SGData']
677    Cell = DATData['Cell']
678   
679    Amat,Bmat = G2lat.cell2AB(Cell[:6])
680    covData = {}
681    pfx = ''
682    if 'covData' in DATData:   
683        covData = DATData['covData']
684        pfx = str(DATData['pId'])+'::'
685    Datoms = []
686    Oatoms = []
687    for i,atom in enumerate(DATData['Datoms']):
688        symop = atom[-1].split('+')
689        if len(symop) == 1:
690            symop.append('0,0,0')       
691        symop[0] = int(symop[0])
692        symop[1] = eval(symop[1])
693        atom.append(symop)
694        Datoms.append(atom)
695        oatom = DATData['Oatoms'][i]
696        names = ['','','']
697        if pfx:
698            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
699        oatom += [names,]
700        Oatoms.append(oatom)
701    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
702    if DATData['Natoms'] == 4:  #torsion
703        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
704        print ' Torsion angle for '+DATData['Name']+' atom sequence: ',atmSeq,'=',G2mth.ValEsd(Tors,sig)
705        print ' NB: Atom sequence determined by selection order'
706        return      # done with torsion
707    elif DATData['Natoms'] == 3:  #angle
708        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
709        print ' Angle in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Ang,sig)
710        print ' NB: Atom sequence determined by selection order'
711    else:   #2 atoms - distance
712        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
713        print ' Distance in '+DATData['Name']+' for atom sequence: ',atmSeq,'=',G2mth.ValEsd(Dist,sig)
714               
715def BestPlane(PlaneData):
716    'Needs a doc string'
717
718    def ShowBanner(name):
719        print 80*'*'
720        print '   Best plane result for phase '+name
721        print 80*'*','\n'
722
723    ShowBanner(PlaneData['Name'])
724
725    Cell = PlaneData['Cell']   
726    Amat,Bmat = G2lat.cell2AB(Cell[:6])       
727    Atoms = PlaneData['Atoms']
728    sumXYZ = np.zeros(3)
729    XYZ = []
730    Natoms = len(Atoms)
731    for atom in Atoms:
732        xyz = np.array(atom[3:6])
733        XYZ.append(xyz)
734        sumXYZ += xyz
735    sumXYZ /= Natoms
736    XYZ = np.array(XYZ)-sumXYZ
737    XYZ = np.inner(Amat,XYZ).T
738    Zmat = np.zeros((3,3))
739    for i,xyz in enumerate(XYZ):
740        Zmat += np.outer(xyz.T,xyz)
741    print ' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2])
742    Evec,Emat = nl.eig(Zmat)
743    Evec = np.sqrt(Evec)/(Natoms-3)
744    Order = np.argsort(Evec)
745    XYZ = np.inner(XYZ,Emat.T).T
746    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
747    print ' Atoms in Cartesian best plane coordinates:'
748    print ' Name         X         Y         Z'
749    for i,xyz in enumerate(XYZ):
750        print ' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2])
751    print '\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]])   
752           
753def main():
754    'Needs a doc string'
755    starttime = time.time()
756    arg = sys.argv
757    if len(arg) > 1:
758        GPXfile = arg[1]
759        if not ospath.exists(GPXfile):
760            print 'ERROR - ',GPXfile," doesn't exist!"
761            exit()
762        Refine(GPXfile,None)
763    else:
764        print 'ERROR - missing filename'
765        exit()
766    print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
767         
768if __name__ == '__main__':
769    main()
Note: See TracBrowser for help on using the repository browser.