source: trunk/GSASIIstrMain.py @ 3378

Last change on this file since 3378 was 3378, checked in by vondreele, 4 years ago

fix seq refine movie of structure
fix hide columns in seq table - saves settings & allows unhide of previously hidden columns
set default hides on Back, Dij & dAx,y,z parameters

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