source: trunk/GSASIIstrMain.py @ 3049

Last change on this file since 3049 was 3049, checked in by toby, 4 years ago

reset InitMP after SelectConfigSetting?.OnApplyChanges? to allow number of cores to change on the fly

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