source: trunk/GSASIIstrMain.py @ 3023

Last change on this file since 3023 was 3023, checked in by odonnell, 6 years ago

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