source: trunk/GSASIIstrMain.py @ 3371

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

constraint fixes: implement wild-card equivalences, fix use of all for phase selection in atom vars; prevent use of Edit Constr menu items for sym-generated equivalences

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