source: trunk/GSASIIstrMain.py @ 4759

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

reduce printed precision of bon valence sum from 3 to 2.
connect +/- keys to z-step in structure plots; can be repeated. Updates view point & status line.
correct status line for n/p atom selection on structure plots; now correctly follows atom selection.
  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 52.5 KB
Line 
1# -*- coding: utf-8 -*-
2'''
3*GSASIIstrMain: main structure routine*
4---------------------------------------
5
6'''
7########### SVN repository information ###################
8# $Date: 2021-01-07 17:21:15 +0000 (Thu, 07 Jan 2021) $
9# $Author: toby $
10# $Revision: 4735 $
11# $URL: trunk/GSASIIstrMain.py $
12# $Id: GSASIIstrMain.py 4735 2021-01-07 17:21:15Z 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: 4735 $")
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
38import GSASIIfiles as G2fil
39import GSASIIElem as G2elem
40import atmdata
41
42sind = lambda x: np.sin(x*np.pi/180.)
43cosd = lambda x: np.cos(x*np.pi/180.)
44tand = lambda x: np.tan(x*np.pi/180.)
45asind = lambda x: 180.*np.arcsin(x)/np.pi
46acosd = lambda x: 180.*np.arccos(x)/np.pi
47atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
48
49ateln2 = 8.0*math.log(2.0)
50DEBUG = True
51
52def ReportProblems(result,Rvals,varyList):
53    '''Create a message based results from the refinement
54    '''
55    #report on SVD 0's and highly correlated variables
56    msg = ''
57    SVD0 = result[2].get('SVD0')
58    if SVD0 > 0: 
59        msg += 'Warning: There were {} singularities in the Hessian'.format(SVD0)
60    # process singular variables
61    psing = result[2].get('psing',[])
62    if len(psing):
63        if msg: msg += '\n'
64        m = '{} Parameters dropped due to singularities:'.format(len(psing))
65        msg += m
66        G2fil.G2Print(m, mode='warn')
67        m = ''
68        for i,val in enumerate(psing):
69            if i == 0:
70                msg += '\n{}'.format(varyList[val])
71                m = '  {}'.format(varyList[val])
72            else:
73                if len(m) > 70:
74                    G2fil.G2Print(m, mode='warn')
75                    m = '  '
76                else:
77                    m += ', '
78                m += '{}'.format(varyList[val])
79                if i == 10:
80                    msg += ', {}... see console for full list'.format(varyList[val])
81                elif i > 10:
82                    pass
83                else:
84                    msg += ', {}'.format(varyList[val])
85        if m: G2fil.G2Print(m, mode='warn')
86    #report on highly correlated variables
87    Hcorr = result[2].get('Hcorr',[])
88    if len(Hcorr) > 0: 
89        if msg: msg += '\n'
90        m = 'Note highly correlated parameters:'
91        G2fil.G2Print(m, mode='warn')
92        msg += m
93    elif SVD0 > 0:
94        if msg: msg += '\n'
95        m = 'Check covariance matrix for parameter correlation'
96        G2fil.G2Print(m, mode='warn')
97        msg += m
98    for i,(v1,v2,corr) in enumerate(Hcorr):
99        if corr > .95:
100            stars = '**'
101        else:
102            stars = '   '
103        m = ' {} {} and {} (@{:.2f}%)'.format(
104            stars,varyList[v1],varyList[v2],100.*corr)
105        G2fil.G2Print(m, mode='warn')
106        if i == 5:
107            msg += '\n' + m
108            msg += '\n   ... check console for more'
109        elif i < 5:
110            msg += '\n' + m
111    if msg:
112        if 'msg' not in Rvals: Rvals['msg'] = ''
113        Rvals['msg'] += msg
114
115def RefineCore(Controls,Histograms,Phases,restraintDict,rigidbodyDict,parmDict,varyList,
116    calcControls,pawleyLookup,ifSeq,printFile,dlg,refPlotUpdate=None):
117    '''Core optimization routines, shared between SeqRefine and Refine
118
119    :returns: 5-tuple of ifOk (bool), Rvals (dict), result, covMatrix, sig
120    '''
121    #patch (added Oct 2020) convert variable names for parm limits to G2VarObj
122    import GSASIIscriptable as G2sc
123    G2sc.patchControls(Controls)
124    # end patch
125#    print 'current',varyList
126#    for item in parmDict: print item,parmDict[item] ######### show dict just before refinement
127    ifPrint = True
128    if ifSeq:
129        ifPrint = False
130    Rvals = {}
131    chisq0 = None
132    while True:
133        G2mv.Map2Dict(parmDict,varyList)
134        begin = time.time()
135        values =  np.array(G2stMth.Dict2Values(parmDict, varyList))
136        if np.any(np.isnan(values)):
137            raise G2obj.G2Exception('ERROR - nan found in LS parameters - use Calculate/View LS parms to locate')
138        # test code to compute GOF and save for external repeat
139        #args = ([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg)
140        #print '*** before fit chi**2',np.sum(G2stMth.errRefine(values,*args)**2)
141        #fl = open('beforeFit.cpickle','wb')
142        #cPickle.dump(values,fl,1)
143        #cPickle.dump(args[:-1],fl,1)
144        #fl.close()
145        Ftol = Controls['min dM/M']
146        Xtol = Controls['SVDtol']
147        Factor = Controls['shift factor']
148        if 'Jacobian' in Controls['deriv type']:
149            result = so.leastsq(G2stMth.errRefine,values,Dfun=G2stMth.dervRefine,full_output=True,
150                ftol=Ftol,col_deriv=True,factor=Factor,
151                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
152            ncyc = int(result[2]['nfev']/2)
153            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot after completion
154        elif 'analytic Hessian' in Controls['deriv type']:
155            Lamda = Controls.get('Marquardt',-3)
156            maxCyc = Controls['max cyc']
157            result = G2mth.HessianLSQ(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,lamda=Lamda,
158                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
159                refPlotUpdate=refPlotUpdate)
160            ncyc = result[2]['num cyc']+1
161            Rvals['lamMax'] = result[2]['lamMax']
162            if 'msg' in result[2]:
163                Rvals['msg'] = result[2]['msg']
164            Controls['Marquardt'] = -3  #reset to default
165            if 'chisq0' in result[2] and chisq0 is None:
166                chisq0 = result[2]['chisq0']
167        elif 'Hessian SVD' in Controls['deriv type']:
168            maxCyc = Controls['max cyc']
169            result = G2mth.HessianSVD(G2stMth.errRefine,values,Hess=G2stMth.HessRefine,ftol=Ftol,xtol=Xtol,maxcyc=maxCyc,Print=ifPrint,
170                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg),
171                refPlotUpdate=refPlotUpdate)
172            if result[1] is None:
173                IfOK = False
174                covMatrix = []
175                sig = len(varyList)*[None,]
176                break
177            ncyc = result[2]['num cyc']+1
178            if 'chisq0' in result[2] and chisq0 is None:
179                chisq0 = result[2]['chisq0']
180        else:           #'numeric'
181            result = so.leastsq(G2stMth.errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor,
182                args=([Histograms,Phases,restraintDict,rigidbodyDict],parmDict,varyList,calcControls,pawleyLookup,dlg))
183            ncyc = 1
184            if len(varyList):
185                ncyc = int(result[2]['nfev']/len(varyList))
186            if refPlotUpdate is not None: refPlotUpdate(Histograms)   # update plot
187#        table = dict(zip(varyList,zip(values,result[0],(result[0]-values))))
188#        for item in table: print item,table[item]               #useful debug - are things shifting?
189        runtime = time.time()-begin
190        Rvals['SVD0'] = result[2].get('SVD0',0)
191        Rvals['converged'] = result[2].get('Converged')
192        Rvals['DelChi2'] = result[2].get('DelChi2',-1.)
193        Rvals['chisq'] = np.sum(result[2]['fvec']**2)
194        G2stMth.Values2Dict(parmDict, varyList, result[0])
195        G2mv.Dict2Map(parmDict,varyList)
196        Rvals['Nobs'] = Histograms['Nobs']
197        Rvals['Rwp'] = np.sqrt(Rvals['chisq']/Histograms['sumwYo'])*100.      #to %
198        Rvals['GOF'] = np.sqrt(Rvals['chisq']/(Histograms['Nobs']-len(varyList)))
199        printFile.write(' Number of function calls: %d No. of observations: %d No. of parameters: %d User rejected: %d Sp. gp. extinct: %d\n'%  \
200            (result[2]['nfev'],Histograms['Nobs'],len(varyList),Histograms['Nrej'],Histograms['Next']))
201        if ncyc:
202            printFile.write(' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles\n'%(runtime,runtime/ncyc,ncyc))
203        printFile.write(' wR = %7.2f%%, chi**2 = %12.6g, GOF = %6.2f\n'%(Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']))
204        sig = len(varyList)*[None,]
205        if 'None' in str(type(result[1])) and ifSeq:    #this bails out of a sequential refinement on singular matrix
206            IfOK = False
207            covMatrix = []
208            G2fil.G2Print ('Warning: **** Refinement failed - singular matrix ****')
209            if 'Hessian' in Controls['deriv type']:
210                num = len(varyList)-1
211                for i,val in enumerate(np.flipud(result[2]['psing'])):
212                    if val:
213                        G2fil.G2Print('Bad parameter: '+varyList[num-i],mode='warn')
214            else:
215                Ipvt = result[2]['ipvt']
216                for i,ipvt in enumerate(Ipvt):
217                    if not np.sum(result[2]['fjac'],axis=1)[i]:
218                        G2fil.G2Print('Bad parameter: '+varyList[ipvt-1],mode='warn')
219            break
220        IfOK = True
221        try:
222            covMatrix = result[1]*Rvals['GOF']**2
223            sig = np.sqrt(np.diag(covMatrix))
224            Lastshft = result[2].get('Xvec',None)
225            if Lastshft is None:
226                Rvals['Max shft/sig'] = 0.0
227            else:
228                Rvals['Max shft/sig'] = np.max(np.nan_to_num(Lastshft/sig))
229            if np.any(np.isnan(sig)) or not sig.shape:
230                G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***',mode='error')
231#            table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig)))
232#            for item in table: print item,table[item]               #useful debug - are things shifting?
233            # report on refinement issues. Result in Rvals['msg']
234            ReportProblems(result,Rvals,varyList)
235            break                   #refinement succeeded - finish up!
236        except TypeError:
237            # if we get here, no result[1] (covar matrix) was returned or other calc error: refinement failed
238            IfOK = False
239            if not len(varyList):
240                covMatrix = []
241                break
242            if 'Hessian' in Controls['deriv type']:
243                SVD0 = result[2].get('SVD0')
244                if SVD0 == -1:
245                    G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
246                elif SVD0 == -2:
247                    G2fil.G2Print ('**** Refinement failed - other problem ****',mode='error')
248                elif SVD0 > 0:
249                    G2fil.G2Print ('**** Refinement failed with {} SVD singularities ****'.format(SVD0),mode='error')
250                else:
251                    G2fil.G2Print ('**** Refinement failed ****',mode='error')
252                if result[1] is None:
253                    IfOK = False
254                    covMatrix = []
255                    sig = len(varyList)*[None,]
256                # report on highly correlated variables
257                ReportProblems(result,Rvals,varyList)
258                # process singular variables
259                psing = result[2].get('psing',[])
260                if dlg: break # refining interactively
261                # non-interactive refinement
262                for val in sorted(psing,reverse=True):
263                    G2fil.G2Print ('Removing parameter: '+varyList[val])
264                    del(varyList[val])
265                if not psing: break    # removed variable(s), try again
266            else:
267                G2fil.G2Print ('**** Refinement failed - singular matrix ****',mode='error')
268                Ipvt = result[2]['ipvt']
269                for i,ipvt in enumerate(Ipvt):
270                    if not np.sum(result[2]['fjac'],axis=1)[i]:
271                        G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1])
272                        del(varyList[ipvt-1])
273                        break
274    if IfOK:
275        G2stMth.GetFobsSq(Histograms,Phases,parmDict,calcControls)
276    if chisq0 is not None:
277        Rvals['GOF0'] = np.sqrt(chisq0/(Histograms['Nobs']-len(varyList)))
278    return IfOK,Rvals,result,covMatrix,sig
279
280def Refine(GPXfile,dlg=None,makeBack=True,refPlotUpdate=None):
281    '''Global refinement -- refines to minimize against all histograms.
282    This can be called in one of three ways, from :meth:`GSASIIdataGUI.GSASII.OnRefine` in an
283    interactive refinement, where dlg will be a wx.ProgressDialog, or non-interactively from
284    :meth:`GSASIIscriptable.G2Project.refine` or from :func:`main`, where dlg will be None.
285    '''
286    import GSASIImpsubs as G2mp
287    G2mp.InitMP()
288    import pytexture as ptx
289    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
290
291    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
292    G2stIO.ShowBanner(printFile)
293    varyList = []
294    parmDict = {}
295    G2mv.InitVars()
296    Controls = G2stIO.GetControls(GPXfile)
297    G2stIO.ShowControls(Controls,printFile)
298    calcControls = {}
299    calcControls.update(Controls)
300    constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
301    restraintDict = G2stIO.GetRestraints(GPXfile)
302    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
303    if not Phases:
304        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
305        G2fil.G2Print (' *** Refine aborted ***')
306        return False,{'msg':'No phases'}
307    if not Histograms:
308        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
309        G2fil.G2Print (' *** Refine aborted ***')
310        return False,{'msg':'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,
315         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,pFile=printFile)
316    calcControls['atomIndx'] = atomIndx
317    calcControls['Natoms'] = Natoms
318    calcControls['FFtables'] = FFtables
319    calcControls['BLtables'] = BLtables
320    calcControls['MFtables'] = MFtables
321    calcControls['maxSSwave'] = maxSSwave
322    hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histograms,pFile=printFile)
323    TwConstr,TwFixed = G2stIO.makeTwinFrConstr(Phases,Histograms,hapVary)
324    constrDict += TwConstr
325    fixedList += TwFixed
326    calcControls.update(controlDict)
327    histVary,histDict,controlDict = G2stIO.GetHistogramData(Histograms,pFile=printFile)
328    calcControls.update(controlDict)
329    varyList = rbVary+phaseVary+hapVary+histVary
330    parmDict.update(rbDict)
331    parmDict.update(phaseDict)
332    parmDict.update(hapDict)
333    parmDict.update(histDict)
334    G2stIO.GetFprime(calcControls,Histograms)
335    # do constraint processing
336    varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
337    msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
338    if msg:
339        return False,{'msg':'Unable to interpret multiplier(s): '+msg}
340    try:
341        G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict)
342        #print(G2mv.VarRemapShow(varyList))
343        #print('DependentVars',G2mv.GetDependentVars())
344        #print('IndependentVars',G2mv.GetIndependentVars())
345    except G2mv.ConstraintException:
346        G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
347        #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
348        #print 'Errors',errmsg
349        #if warnmsg: print 'Warnings',warnmsg
350        return False,{'msg':' Constraint error'}
351#    print G2mv.VarRemapShow(varyList)
352
353    # remove frozen vars from refinement
354    if 'parmFrozen' not in Controls:
355        Controls['parmFrozen'] = {}
356    if 'FrozenList' not in Controls['parmFrozen']: 
357        Controls['parmFrozen']['FrozenList'] = []
358    parmFrozenList = Controls['parmFrozen']['FrozenList']
359    frozenList = [i for i in varyList if i in parmFrozenList]
360    if len(frozenList) != 0: 
361        varyList = [i for i in varyList if i not in parmFrozenList]
362        G2fil.G2Print(
363            'Frozen refined variables (due to exceeding limits)\n\t:{}'
364            .format(frozenList))
365       
366    ifSeq = False
367    printFile.write('\n Refinement results:\n')
368    printFile.write(135*'-'+'\n')
369    Rvals = {}
370    try:
371        covData = {}
372        IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histograms,Phases,restraintDict,
373            rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
374            refPlotUpdate=refPlotUpdate)
375        if IfOK:
376            sigDict = dict(zip(varyList,sig))
377            newCellDict = G2stMth.GetNewCellParms(parmDict,varyList)
378            newAtomDict = G2stMth.ApplyXYZshifts(parmDict,varyList)
379            covData = {'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
380                       'varyListStart':varyListStart,
381                       'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,
382                       'newCellDict':newCellDict,'freshCOV':True}
383            # add the uncertainties into the esd dictionary (sigDict)
384            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
385            # check for variables outside their allowed range, reset and freeze them
386            frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
387            G2mv.PrintIndependentVars(parmDict,varyList,sigDict,pFile=printFile)
388            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
389            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
390            G2stIO.SetPhaseData(parmDict,sigDict,Phases,rbIds,covData,restraintDict,printFile)
391            G2stIO.PrintISOmodes(printFile,Phases,parmDict,sigDict)
392            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,pFile=printFile)
393            G2stIO.SetHistogramData(parmDict,sigDict,Histograms,calcControls,pFile=printFile)
394            if len(frozen):
395                if 'msg' in Rvals:
396                    Rvals['msg'] += '\n'
397                else:
398                    Rvals['msg'] = ''
399                msg = ('Warning: {} variable(s) refined outside limits and were frozen ({} total frozen)'
400                    .format(len(frozen),len(parmFrozenList))
401                    )
402                G2fil.G2Print(msg)
403                Rvals['msg'] += msg
404            elif len(parmFrozenList):
405                if 'msg' in Rvals:
406                    Rvals['msg'] += '\n'
407                else:
408                    Rvals['msg'] = ''
409                msg = ('Note: a total of {} variable(s) are frozen due to refining outside limits'
410                    .format(len(parmFrozenList))
411                    )
412                G2fil.G2Print('Note: ',msg)
413                Rvals['msg'] += msg
414            G2stIO.SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,rigidbodyDict,covData,parmFrozenList,makeBack)
415            printFile.close()
416            G2fil.G2Print (' Refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
417            G2fil.G2Print (' ***** Refinement successful *****')
418        else:
419            G2fil.G2Print ('****ERROR - Refinement failed',mode='error')
420            if 'msg' in Rvals:
421                G2fil.G2Print ('Note refinement problem:',mode='warn')
422                G2fil.G2Print (Rvals['msg'],mode='warn')
423            raise G2obj.G2Exception('**** ERROR: Refinement failed ****')
424    except G2obj.G2RefineCancel as Msg:
425        printFile.close()
426        G2fil.G2Print (' ***** Refinement stopped *****')
427        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
428        if 'msg' in Rvals:
429            Rvals['msg'] += '\n'
430            Rvals['msg'] += Msg.msg
431            if not dlg:
432                G2fil.G2Print ('Note refinement problem:',mode='warn')
433                G2fil.G2Print (Rvals['msg'],mode='warn')
434        else:
435            Rvals['msg'] = Msg.msg
436        return False,Rvals
437#    except G2obj.G2Exception as Msg:  # cell metric error, others?
438    except Exception as Msg:  # cell metric error, others?
439        if GSASIIpath.GetConfigValue('debug'):
440            import traceback
441            print(traceback.format_exc())       
442        if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
443        printFile.close()
444        G2fil.G2Print (' ***** Refinement error *****')
445        if 'msg' in Rvals:
446            Rvals['msg'] += '\n\n'
447            Rvals['msg'] += Msg.msg
448            if not dlg:
449                G2fil.G2Print ('Note refinement problem:',mode='warn')
450                G2fil.G2Print (Rvals['msg'],mode='warn')
451        else:
452            Rvals['msg'] = Msg.msg
453        return False,Rvals
454
455#for testing purposes, create a file for testderiv
456    if GSASIIpath.GetConfigValue('debug'):   # and IfOK:
457#needs: values,HistoPhases,parmDict,varylist,calcControls,pawleyLookup
458        fl = open(ospath.splitext(GPXfile)[0]+'.testDeriv','wb')
459        cPickle.dump(result[0],fl,1)
460        cPickle.dump([Histograms,Phases,restraintDict,rigidbodyDict],fl,1)
461        cPickle.dump([constrDict,fixedList,G2mv.GetDependentVars()],fl,1)
462        cPickle.dump(parmDict,fl,1)
463        cPickle.dump(varyList,fl,1)
464        cPickle.dump(calcControls,fl,1)
465        cPickle.dump(pawleyLookup,fl,1)
466        fl.close()
467    if dlg:
468        return True,Rvals
469    elif 'msg' in Rvals:
470        G2fil.G2Print ('Reported from refinement:',mode='warn')
471        G2fil.G2Print (Rvals['msg'],mode='warn')
472
473def phaseCheck(phaseVary,Phases,histogram):
474    '''
475    Removes unused parameters from phase varylist if phase not in histogram
476    #TODO - implement "Fix FXU" for seq refinement here - done?
477    '''
478    NewVary = []
479    for phase in Phases:
480        if histogram not in Phases[phase]['Histograms']: continue
481        if Phases[phase]['Histograms'][histogram]['Use']:
482            pId = Phases[phase]['pId']
483            newVary = [item for item in phaseVary if item.split(':')[0] == str(pId)]
484            FixVals = Phases[phase]['Histograms'][histogram].get('Fix FXU',' ')
485            if 'F' in FixVals:
486                newVary = [item for item in newVary if not 'Afrac' in item]
487            if 'X' in FixVals:
488                newVary = [item for item in newVary if not 'dA' in item]
489            if 'U' in FixVals:
490                newVary = [item for item in newVary if not 'AU' in item]
491            if 'M' in FixVals:
492                newVary = [item for item in newVary if not 'AM' in item]
493            NewVary += newVary
494    return NewVary
495
496def SeqRefine(GPXfile,dlg,refPlotUpdate=None):
497    '''Perform a sequential refinement -- cycles through all selected histgrams,
498    one at a time
499    '''
500    import GSASIImpsubs as G2mp
501    G2mp.InitMP()
502    import pytexture as ptx
503    ptx.pyqlmninit()            #initialize fortran arrays for spherical harmonics
504    msgs = {}
505    printFile = open(ospath.splitext(GPXfile)[0]+'.lst','w')
506    G2fil.G2Print ('Starting Sequential Refinement')
507    G2stIO.ShowBanner(printFile)
508    Controls = G2stIO.GetControls(GPXfile)
509    preFrozenCount = 0
510    for h in Controls['parmFrozen']:
511        if h == 'FrozenList': continue
512        preFrozenCount += len(Controls['parmFrozen'][h])   
513    G2stIO.ShowControls(Controls,printFile,SeqRef=True,
514                            preFrozenCount=preFrozenCount)
515    restraintDict = G2stIO.GetRestraints(GPXfile)
516    Histograms,Phases = G2stIO.GetUsedHistogramsAndPhases(GPXfile)
517    if not Phases:
518        G2fil.G2Print (' *** ERROR - you have no phases to refine! ***')
519        G2fil.G2Print (' *** Refine aborted ***')
520        return False,'No phases'
521    if not Histograms:
522        G2fil.G2Print (' *** ERROR - you have no data to refine with! ***')
523        G2fil.G2Print (' *** Refine aborted ***')
524        return False,'No data'
525    rigidbodyDict = G2stIO.GetRigidBodies(GPXfile)
526    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
527    rbVary,rbDict = G2stIO.GetRigidBodyModels(rigidbodyDict,pFile=printFile)
528    G2mv.InitVars()
529    (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,
530         maxSSwave) = G2stIO.GetPhaseData(Phases,restraintDict,rbIds,
531                                    Print=False,pFile=printFile,seqRef=True)
532    for item in phaseVary:
533        if '::A0' in item:
534            G2fil.G2Print ('**** WARNING - lattice parameters should not be refined in a sequential refinement ****')
535            G2fil.G2Print ('****           instead use the Dij parameters for each powder histogram            ****')
536            return False,'Lattice parameter refinement error - see console message'
537        if '::C(' in item:
538            G2fil.G2Print ('**** WARNING - phase texture parameters should not be refined in a sequential refinement ****')
539            G2fil.G2Print ('****           instead use the C(L,N) parameters for each powder histogram               ****')
540            return False,'Phase texture refinement error - see console message'
541    if 'Seq Data' in Controls:
542        histNames = Controls['Seq Data']
543    else: # patch from before Controls['Seq Data'] was implemented?
544        histNames = G2stIO.GetHistogramNames(GPXfile,['PWDR',])
545    if Controls.get('Reverse Seq'):
546        histNames.reverse()
547    SeqResult = G2stIO.GetSeqResult(GPXfile)
548#    SeqResult = {'SeqPseudoVars':{},'SeqParFitEqList':[]}
549    Histo = {}
550    NewparmDict = {}
551    G2stIO.SetupSeqSavePhases(GPXfile)
552    msgs['steepestNum'] = 0
553    msgs['maxshift/sigma'] = []
554    for ihst,histogram in enumerate(histNames):
555        if GSASIIpath.GetConfigValue('Show_timing'): t1 = time.time()
556        G2fil.G2Print('\nRefining with '+str(histogram))
557        G2mv.InitVars()
558        (Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,
559             FFtables,BLtables,MFtables,maxSSwave) = G2stIO.GetPhaseData(
560                 Phases,restraintDict,rbIds,
561                 Print=False,pFile=printFile,seqRef=True)
562        ifPrint = False
563        if dlg:
564            dlg.SetTitle('Residual for histogram '+str(ihst))
565        calcControls = {}
566        calcControls['atomIndx'] = atomIndx
567        calcControls['Natoms'] = Natoms
568        calcControls['FFtables'] = FFtables
569        calcControls['BLtables'] = BLtables
570        calcControls['MFtables'] = MFtables
571        calcControls['maxSSwave'] = maxSSwave
572        if histogram not in Histograms:
573            G2fil.G2Print("Error: not found!")
574            continue
575    #TODO - implement "Fix FXU" for seq refinement here - done?
576        hId = Histograms[histogram]['hId']
577        redphaseVary = phaseCheck(phaseVary,Phases,histogram)
578        Histo = {histogram:Histograms[histogram],}
579        hapVary,hapDict,controlDict = G2stIO.GetHistogramPhaseData(Phases,Histo,Print=False)
580        calcControls.update(controlDict)
581        histVary,histDict,controlDict = G2stIO.GetHistogramData(Histo,False)
582        calcControls.update(controlDict)
583        varyList = rbVary+redphaseVary+hapVary+histVary
584#        if not ihst:
585            # save the initial vary list, but without histogram numbers on parameters
586        saveVaryList = varyList[:]
587        for i,item in enumerate(saveVaryList):
588            items = item.split(':')
589            if items[1]:
590                items[1] = ''
591            item = ':'.join(items)
592            saveVaryList[i] = item
593        if not ihst:
594            SeqResult['varyList'] = saveVaryList
595        else:
596            SeqResult['varyList'] = list(set(SeqResult['varyList']+saveVaryList))
597        parmDict = {}
598        parmDict.update(rbDict)
599        parmDict.update(phaseDict)
600        parmDict.update(hapDict)
601        parmDict.update(histDict)
602        if Controls['Copy2Next']:   # update with parms from last histogram
603            #parmDict.update(NewparmDict) # don't use in case extra entries would cause a problem
604            for parm in NewparmDict:
605                if parm in parmDict:
606                    parmDict[parm] = NewparmDict[parm]
607        elif histogram in SeqResult:  # update phase from last seq ref
608            NewparmDict = SeqResult[histogram].get('parmDict',{})
609            for parm in NewparmDict:
610                if '::' in parm and parm in parmDict:
611                    parmDict[parm] = NewparmDict[parm]
612           
613        G2stIO.GetFprime(calcControls,Histo)
614        # do constraint processing
615        #reload(G2mv) # debug
616        constrDict,fixedList = G2stIO.GetConstraints(GPXfile)
617        varyListStart = tuple(varyList) # save the original varyList before dependent vars are removed
618        msg = G2mv.EvaluateMultipliers(constrDict,parmDict)
619        if msg:
620            return False,'Unable to interpret multiplier(s): '+msg
621        try:
622            groups,parmlist = G2mv.GenerateConstraints(varyList,constrDict,fixedList,parmDict,SeqHist=hId)
623#            if GSASIIpath.GetConfigValue('debug'): print("DBG_"+
624#                G2mv.VarRemapShow(varyList,True))
625#            print('DependentVars',G2mv.GetDependentVars())
626#            print('IndependentVars',G2mv.GetIndependentVars())
627            constraintInfo = (groups,parmlist,constrDict,fixedList,ihst)
628        except G2mv.ConstraintException:
629            G2fil.G2Print (' *** ERROR - your constraints are internally inconsistent ***')
630            #errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
631            #print 'Errors',errmsg
632            #if warnmsg: print 'Warnings',warnmsg
633            return False,' Constraint error'
634        #print G2mv.VarRemapShow(varyList)
635        if not ihst:
636            # first histogram to refine against
637            firstVaryList = []
638            for item in varyList:
639                items = item.split(':')
640                if items[1]:
641                    items[1] = ''
642                item = ':'.join(items)
643                firstVaryList.append(item)
644            newVaryList = firstVaryList
645        else:
646            newVaryList = []
647            for item in varyList:
648                items = item.split(':')
649                if items[1]:
650                    items[1] = ''
651                item = ':'.join(items)
652                newVaryList.append(item)
653        if newVaryList != firstVaryList and Controls['Copy2Next']:
654            # variable lists are expected to match between sequential refinements when Copy2Next is on
655            #print '**** ERROR - variable list for this histogram does not match previous'
656            #print '     Copy of variables is not possible'
657            #print '\ncurrent histogram',histogram,'has',len(newVaryList),'variables'
658            combined = list(set(firstVaryList+newVaryList))
659            c = [var for var in combined if var not in newVaryList]
660            p = [var for var in combined if var not in firstVaryList]
661            G2fil.G2Print('*** Variables change ***')
662            for typ,vars in [('Removed',c),('Added',p)]:
663                line = '  '+typ+': '
664                if vars:
665                    for var in vars:
666                        if len(line) > 70:
667                            G2fil.G2Print(line)
668                            line = '    '
669                        line += var + ', '
670                else:
671                        line += 'none, '
672                G2fil.G2Print(line[:-2])
673            firstVaryList = newVaryList
674
675        ifSeq = True
676        printFile.write('\n Refinement results for histogram id {}: {}\n'
677                            .format(hId,histogram))
678        printFile.write(135*'-'+'\n')
679        # remove frozen vars
680        if 'parmFrozen' not in Controls:
681            Controls['parmFrozen'] = {}
682        if histogram not in Controls['parmFrozen']: 
683            Controls['parmFrozen'][histogram] = []
684        parmFrozenList = Controls['parmFrozen'][histogram]
685        frozenList = [i for i in varyList if i in parmFrozenList]
686        if len(frozenList) != 0: 
687           varyList = [i for i in varyList if i not in parmFrozenList]
688           s = ''
689           for a in frozenList:
690               if s:
691                   s+= ', '
692               s += a
693           printFile.write(
694               ' The following refined variables have previously been frozen due to exceeding limits:\n\t{}\n'
695               .format(s))
696        try:
697            IfOK,Rvals,result,covMatrix,sig = RefineCore(Controls,Histo,Phases,restraintDict,
698                rigidbodyDict,parmDict,varyList,calcControls,pawleyLookup,ifSeq,printFile,dlg,
699                refPlotUpdate=refPlotUpdate)
700            G2fil.G2Print ('  wR = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f, last delta chi = %.4f, last shft/sig = %.4f'%(
701                Rvals['Rwp'],Rvals['chisq'],Rvals['GOF']**2,Rvals['DelChi2'],Rvals.get('Max shft/sig',np.nan)))
702            if Rvals.get('lamMax',0) >= 10.:
703                msgs['steepestNum'] += 1
704            if Rvals.get('Max shft/sig'):
705                msgs['maxshift/sigma'].append(Rvals['Max shft/sig'])
706            # add the uncertainties into the esd dictionary (sigDict)
707            if not IfOK:
708                G2fil.G2Print('***** Sequential refinement failed at histogram '+histogram,mode='warn')
709                break
710            sigDict = dict(zip(varyList,sig))
711            # the uncertainties for dependent constrained parms into the esd dict
712            sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict))
713            # check for variables outside their allowed range, reset and freeze them
714            frozen = dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList)
715            msg = None
716            if len(frozen) > 0:
717               msg = ('Hist {}: {} variables were outside limits and were frozen (now {} frozen total)'
718                   .format(ihst,len(frozen),len(parmFrozenList)))
719               G2fil.G2Print(msg)
720               msg = (' {} variables were outside limits and were frozen (now {} frozen total)'
721                   .format(len(frozen),len(parmFrozenList)))
722               for p in frozen:
723                   if p not in varyList:
724                       print('Frozen Warning: {} not in varyList. This should not happen!'.format(p))
725                       continue
726                   i = varyList.index(p)
727                   result[0][i] = parmDict[p]
728                   sig[i] = -0.1
729            # a dict with values & esds for dependent (constrained) parameters - avoid extraneous holds
730            depParmDict = {i:(parmDict[i],sigDict[i]) for i in varyListStart if i in sigDict and i not in varyList}
731            newCellDict = copy.deepcopy(G2stMth.GetNewCellParms(parmDict,varyList))
732            newAtomDict = copy.deepcopy(G2stMth.ApplyXYZshifts(parmDict,varyList))
733            histRefData = {
734                'variables':result[0],'varyList':varyList,'sig':sig,'Rvals':Rvals,
735                'varyListStart':varyListStart,
736                'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,
737                'newCellDict':newCellDict,'depParmDict':depParmDict,
738                'constraintInfo':constraintInfo,
739                'parmDict':parmDict,
740                }
741            SeqResult[histogram] = histRefData
742            G2stMth.ApplyRBModels(parmDict,Phases,rigidbodyDict,True)
743#            G2stIO.SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,printFile)
744            G2stIO.SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,None,ifPrint,printFile)
745            G2stIO.SetHistogramData(parmDict,sigDict,Histo,None,ifPrint,printFile,seq=True)
746            G2stIO.SaveUpdatedHistogramsAndPhases(GPXfile,Histo,Phases,rigidbodyDict,histRefData,Controls['parmFrozen'])
747            if msg: 
748                printFile.write(msg+'\n')
749            NewparmDict = {}
750            # make dict of varied parameters in current histogram, renamed to
751            # next histogram, for use in next refinement.
752            if Controls['Copy2Next'] and ihst < len(histNames)-1:
753                hId = Histo[histogram]['hId'] # current histogram
754                nexthId = Histograms[histNames[ihst+1]]['hId']
755                for parm in set(list(varyList)+list(varyListStart)):
756                    items = parm.split(':')
757                    if len(items) < 3: 
758                        continue
759                    if str(hId) in items[1]:
760                        items[1] = str(nexthId)
761                        newparm = ':'.join(items)
762                        NewparmDict[newparm] = parmDict[parm]
763                    else:
764                        if items[2].startswith('dA'): parm = parm.replace(':dA',':A') 
765                        NewparmDict[parm] = parmDict[parm]
766                   
767        except G2obj.G2RefineCancel as Msg:
768            if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
769            printFile.close()
770            G2fil.G2Print (' ***** Refinement stopped *****')
771            return False,Msg.msg
772        except G2obj.G2Exception as Msg:  # cell metric error, others?
773            if not hasattr(Msg,'msg'): Msg.msg = str(Msg)
774            printFile.close()
775            G2fil.G2Print (' ***** Refinement error *****')
776            return False,Msg.msg
777        if GSASIIpath.GetConfigValue('Show_timing'):
778            t2 = time.time()
779            G2fil.G2Print("Fit step time {:.2f} sec.".format(t2-t1))
780            t1 = t2
781    SeqResult['histNames'] = [itm for itm in G2stIO.GetHistogramNames(GPXfile,['PWDR',]) if itm in SeqResult.keys()]
782    try:
783        G2stIO.SetSeqResult(GPXfile,Histograms,SeqResult)
784    except Exception as msg:
785        print('Error reading Sequential results\n',str(msg))
786        if GSASIIpath.GetConfigValue('debug'):
787            import traceback
788            print(traceback.format_exc())       
789    postFrozenCount = 0
790    for h in Controls['parmFrozen']:
791        if h == 'FrozenList': continue
792        postFrozenCount += len(Controls['parmFrozen'][h])
793    if postFrozenCount:
794        msgs['Frozen'] = 'Ending refinement with {} Frozen variables ({} added now)\n'.format(postFrozenCount,postFrozenCount-preFrozenCount)
795        printFile.write('\n'+msgs['Frozen'])
796    printFile.close()
797    G2fil.G2Print (' Sequential refinement results are in file: '+ospath.splitext(GPXfile)[0]+'.lst')
798    G2fil.G2Print (' ***** Sequential refinement successful *****')
799    return True,msgs
800
801def dropOOBvars(varyList,parmDict,sigDict,Controls,parmFrozenList):
802    '''Find variables in the parameters dict that are outside the ranges
803    (in parmMinDict and parmMaxDict) and set them to the limits values.
804    Add any such variables into the list of frozen variable
805    (parmFrozenList). Returns a list of newly frozen variables, if any.
806    '''
807    parmMinDict = Controls.get('parmMinDict',{})
808    parmMaxDict = Controls.get('parmMaxDict',{})
809    freeze = []
810    if parmMinDict or parmMaxDict:
811        for name in varyList:
812            if name not in parmDict: continue
813            n,val = G2obj.prmLookup(name,parmMinDict)
814            if n is not None:
815                if parmDict[name] < parmMinDict[n]:
816                    parmDict[name] = parmMinDict[n]
817                    sigDict[name] = 0.0
818                    freeze.append(name)
819                    continue
820            n,val = G2obj.prmLookup(name,parmMaxDict)
821            if n is not None:
822                if parmDict[name] > parmMaxDict[n]:
823                    parmDict[name] = parmMaxDict[n]
824                    sigDict[name] = 0.0
825                    freeze.append(name)
826                    continue
827        for v in freeze:
828            if v not in parmFrozenList:
829                parmFrozenList.append(v)
830    return freeze
831
832def RetDistAngle(DisAglCtls,DisAglData,dlg=None):
833    '''Compute and return distances and angles
834
835    :param dict DisAglCtls: contains distance/angle radii usually defined using
836       :func:`GSASIIctrlGUI.DisAglDialog`
837    :param dict DisAglData: contains phase data:
838       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
839       for distance/angle origins and atoms to be used as targets.
840       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
841
842    :returns: AtomLabels,DistArray,AngArray where:
843
844      **AtomLabels** is a dict of atom labels, keys are the atom number
845
846      **DistArray** is a dict keyed by the origin atom number where the value is a list
847      of distance entries. The value for each distance is a list containing:
848
849        0) the target atom number (int);
850        1) the unit cell offsets added to x,y & z (tuple of int values)
851        2) the symmetry operator number (which may be modified to indicate centering and center of symmetry)
852        3) an interatomic distance in A (float)
853        4) an uncertainty on the distance in A or 0.0 (float)
854
855      **AngArray** is a dict keyed by the origin (central) atom number where
856      the value is a list of
857      angle entries. The value for each angle entry consists of three values:
858
859        0) a distance item reference for one neighbor (int)
860        1) a distance item reference for a second neighbor (int)
861        2) a angle, uncertainty pair; the s.u. may be zero (tuple of two floats)
862
863      The AngArray distance reference items refer directly to the index of the items in the
864      DistArray item for the list of distances for the central atom.
865    '''
866    import numpy.ma as ma
867
868    SGData = DisAglData['SGData']
869    Cell = DisAglData['Cell']
870
871    Amat,Bmat = G2lat.cell2AB(Cell[:6])
872    covData = {}
873    if len(DisAglData.get('covData',{})):
874        covData = DisAglData['covData']
875        covMatrix = covData['covMatrix']
876        varyList = covData['varyList']
877        pfx = str(DisAglData['pId'])+'::'
878
879    Factor = DisAglCtls['Factors']
880    Radii = dict(zip(DisAglCtls['AtomTypes'],zip(DisAglCtls['BondRadii'],DisAglCtls['AngleRadii'])))
881    indices = (-2,-1,0,1,2)
882    Units = np.array([[h,k,l] for h in indices for k in indices for l in indices])
883    origAtoms = DisAglData['OrigAtoms']
884    targAtoms = DisAglData['TargAtoms']
885    AtomLabels = {}
886    for Oatom in origAtoms:
887        AtomLabels[Oatom[0]] = Oatom[1]
888    for Oatom in targAtoms:
889        AtomLabels[Oatom[0]] = Oatom[1]
890    DistArray = {}
891    AngArray = {}
892    for iO,Oatom in enumerate(origAtoms):
893        DistArray[Oatom[0]] = []
894        AngArray[Oatom[0]] = []
895        OxyzNames = ''
896        IndBlist = []
897        Dist = []
898        Vect = []
899        VectA = []
900        angles = []
901        for Tatom in targAtoms:
902            Xvcov = []
903            TxyzNames = ''
904            if len(DisAglData.get('covData',{})):
905                OxyzNames = [pfx+'dAx:%d'%(Oatom[0]),pfx+'dAy:%d'%(Oatom[0]),pfx+'dAz:%d'%(Oatom[0])]
906                TxyzNames = [pfx+'dAx:%d'%(Tatom[0]),pfx+'dAy:%d'%(Tatom[0]),pfx+'dAz:%d'%(Tatom[0])]
907                Xvcov = G2mth.getVCov(OxyzNames+TxyzNames,varyList,covMatrix)
908            result = G2spc.GenAtom(Tatom[3:6],SGData,False,Move=False)
909            BsumR = (Radii[Oatom[2]][0]+Radii[Tatom[2]][0])*Factor[0]
910            AsumR = (Radii[Oatom[2]][1]+Radii[Tatom[2]][1])*Factor[1]
911            for [Txyz,Top,Tunit,Spn] in result:
912                Dx = (Txyz-np.array(Oatom[3:6]))+Units
913                dx = np.inner(Amat,Dx)
914                dist = ma.masked_less(np.sqrt(np.sum(dx**2,axis=0)),0.5)
915                IndB = ma.nonzero(ma.masked_greater(dist-BsumR,0.))
916                if np.any(IndB):
917                    for indb in IndB:
918                        for i in range(len(indb)):
919                            if str(dx.T[indb][i]) not in IndBlist:
920                                IndBlist.append(str(dx.T[indb][i]))
921                                unit = Units[indb][i]
922                                tunit = (unit[0]+Tunit[0],unit[1]+Tunit[1],unit[2]+Tunit[2])
923                                sig = 0.0
924                                if len(Xvcov):
925                                    pdpx = G2mth.getDistDerv(Oatom[3:6],Tatom[3:6],Amat,unit,Top,SGData)
926                                    sig = np.sqrt(np.inner(pdpx,np.inner(pdpx,Xvcov)))
927                                Dist.append([Oatom[0],Tatom[0],tunit,Top,ma.getdata(dist[indb])[i],sig])
928                                if (Dist[-1][-2]-AsumR) <= 0.:
929                                    Vect.append(dx.T[indb][i]/Dist[-1][-2])
930                                    VectA.append([OxyzNames,np.array(Oatom[3:6]),TxyzNames,np.array(Tatom[3:6]),unit,Top])
931                                else:
932                                    Vect.append([0.,0.,0.])
933                                    VectA.append([])
934        if dlg is not None:
935            dlg.Update(iO,newmsg='Atoms done=%d'%(iO))
936        for D in Dist:
937            DistArray[Oatom[0]].append(D[1:])
938        Vect = np.array(Vect)
939        angles = np.zeros((len(Vect),len(Vect)))
940        angsig = np.zeros((len(Vect),len(Vect)))
941        for i,veca in enumerate(Vect):
942            if np.any(veca):
943                for j,vecb in enumerate(Vect):
944                    if np.any(vecb):
945                        angles[i][j],angsig[i][j] = G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)
946                        if i <= j: continue
947                        AngArray[Oatom[0]].append((i,j,
948                            G2mth.getAngSig(VectA[i],VectA[j],Amat,SGData,covData)))
949    return AtomLabels,DistArray,AngArray
950
951def PrintDistAngle(DisAglCtls,DisAglData,out=sys.stdout):
952    '''Print distances and angles
953
954    :param dict DisAglCtls: contains distance/angle radii usually defined using
955       :func:`GSASIIctrlGUI.DisAglDialog`
956    :param dict DisAglData: contains phase data:
957       Items 'OrigAtoms' and 'TargAtoms' contain the atoms to be used
958       for distance/angle origins and atoms to be used as targets.
959       Item 'SGData' has the space group information (see :ref:`Space Group object<SGData_table>`)
960    :param file out: file object for output. Defaults to sys.stdout.
961    '''
962    def MyPrint(s):
963        out.write(s+'\n')
964        # print(s,file=out) # use in Python 3
965
966    def ShowBanner(name):
967        MyPrint(80*'*')
968        MyPrint('   Interatomic Distances and Angles for phase '+name)
969        MyPrint((80*'*')+'\n')
970
971    ShowBanner(DisAglCtls['Name'])
972    SGData = DisAglData['SGData']
973    SGtext,SGtable = G2spc.SGPrint(SGData)
974    for line in SGtext: MyPrint(line)
975    if len(SGtable) > 1:
976        for i,item in enumerate(SGtable[::2]):
977            if 2*i+1 == len(SGtable):
978                line = ' %s'%(item.ljust(30))
979            else:
980                line = ' %s %s'%(item.ljust(30),SGtable[2*i+1].ljust(30))
981            MyPrint(line)
982    else:
983        MyPrint(' ( 1)    %s'%(SGtable[0])) #triclinic case
984    Cell = DisAglData['Cell']
985
986    Amat,Bmat = G2lat.cell2AB(Cell[:6])
987    covData = {}
988    if len(DisAglData.get('covData',{})):
989        covData = DisAglData['covData']
990        pfx = str(DisAglData['pId'])+'::'
991        A = G2lat.cell2A(Cell[:6])
992        cellSig = G2stIO.getCellEsd(pfx,SGData,A,covData)
993        names = [' a = ',' b = ',' c = ',' alpha = ',' beta = ',' gamma = ',' Volume = ']
994        valEsd = [G2mth.ValEsd(Cell[i],cellSig[i],True) for i in range(7)]
995        line = '\n Unit cell:'
996        for name,vals in zip(names,valEsd):
997            line += name+vals
998        MyPrint(line)
999    else:
1000        MyPrint('\n Unit cell: a = '+('%.5f'%Cell[0])+' b = '+('%.5f'%Cell[1])+' c = '+('%.5f'%Cell[2])+
1001            ' alpha = '+('%.3f'%Cell[3])+' beta = '+('%.3f'%Cell[4])+' gamma = '+
1002            ('%.3f'%Cell[5])+' Volume = '+('%.3f'%Cell[6]))
1003
1004    AtomLabels,DistArray,AngArray = RetDistAngle(DisAglCtls,DisAglData)
1005    origAtoms = DisAglData['OrigAtoms']
1006    for Oatom in origAtoms:
1007        i = Oatom[0]
1008        Dist = DistArray[i]
1009        nDist = len(Dist)
1010        angles = np.zeros((nDist,nDist))
1011        angsig = np.zeros((nDist,nDist))
1012        for k,j,tup in AngArray[i]:
1013            angles[k][j],angsig[k][j] = angles[j][k],angsig[j][k] = tup
1014        line = ''
1015        for i,x in enumerate(Oatom[3:6]):
1016            line += ('%12.5f'%x).rstrip('0')
1017        MyPrint('\n Distances & angles for '+Oatom[1]+' at '+line.rstrip())
1018        MyPrint(80*'*')
1019        line = ''
1020        for dist in Dist[:-1]:
1021            line += '%12s'%(AtomLabels[dist[0]].center(12))
1022        MyPrint('  To       cell +(sym. op.)      dist.  '+line.rstrip())
1023        BVS = {}
1024        BVdat = {}
1025        Otyp = G2elem.FixValence(Oatom[2]).split('+')[0].split('-')[0]
1026        BVox = [BV for BV in atmdata.BVSoxid[Otyp] if '+' in BV]
1027        if len(BVox):
1028            BVS = {BV:0.0 for BV in BVox}
1029            BVdat = {BV:dict(zip(['O','F','Cl'],atmdata.BVScoeff[BV])) for BV in BVox}
1030            pvline = 'Bond Valence sums for: '
1031        for i,dist in enumerate(Dist):
1032            line = ''
1033            for j,angle in enumerate(angles[i][0:i]):
1034                sig = angsig[i][j]
1035                if angle:
1036                    if sig:
1037                        line += '%12s'%(G2mth.ValEsd(angle,sig,True).center(12))
1038                    else:
1039                        val = '%.3f'%(angle)
1040                        line += '%12s'%(val.center(12))
1041                else:
1042                    line += 12*' '
1043            if dist[4]:            #sig exists!
1044                val = G2mth.ValEsd(dist[3],dist[4])
1045            else:
1046                val = '%8.4f'%(dist[3])
1047            if len(BVox):
1048                Tatm = G2elem.FixValence(DisAglData['TargAtoms'][dist[0]][2]).split('-')[0]
1049                if Tatm in ['O','F','Cl']:
1050                    for BV in BVox:
1051                        BVS[BV] += np.exp((BVdat[BV][Tatm]-dist[3])/0.37)               
1052            tunit = '[%2d%2d%2d]'% dist[1]
1053            MyPrint(('  %8s%10s+(%4d) %12s'%(AtomLabels[dist[0]].ljust(8),tunit.ljust(10),dist[2],val.center(12)))+line.rstrip())
1054        if len(BVox):
1055            MyPrint(80*'*')
1056            for BV in BVox:
1057                pvline += ' %s: %.2f  '%(BV,BVS[BV])
1058            MyPrint(pvline)
1059
1060def DisAglTor(DATData):
1061    'Needs a doc string'
1062    SGData = DATData['SGData']
1063    Cell = DATData['Cell']
1064
1065    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1066    covData = {}
1067    pfx = ''
1068    if 'covData' in DATData:
1069        covData = DATData['covData']
1070        pfx = str(DATData['pId'])+'::'
1071    Datoms = []
1072    Oatoms = []
1073    for i,atom in enumerate(DATData['Datoms']):
1074        symop = atom[-1].split('+')
1075        if len(symop) == 1:
1076            symop.append('0,0,0')
1077        symop[0] = int(symop[0])
1078        symop[1] = eval(symop[1])
1079        atom.append(symop)
1080        Datoms.append(atom)
1081        oatom = DATData['Oatoms'][i]
1082        names = ['','','']
1083        if pfx:
1084            names = [pfx+'dAx:'+str(oatom[0]),pfx+'dAy:'+str(oatom[0]),pfx+'dAz:'+str(oatom[0])]
1085        oatom += [names,]
1086        Oatoms.append(oatom)
1087    atmSeq = [atom[1]+'('+atom[-2]+')' for atom in Datoms]
1088    if DATData['Natoms'] == 4:  #torsion
1089        Tors,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1090        G2fil.G2Print (' Torsion angle for %s atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Tors,sig)))
1091        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1092        return      # done with torsion
1093    elif DATData['Natoms'] == 3:  #angle
1094        Ang,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1095        G2fil.G2Print (' Angle in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Ang,sig)))
1096        G2fil.G2Print (' NB: Atom sequence determined by selection order')
1097    else:   #2 atoms - distance
1098        Dist,sig = G2mth.GetDATSig(Oatoms,Datoms,Amat,SGData,covData)
1099        G2fil.G2Print (' Distance in %s for atom sequence: %s = %s'%(DATData['Name'],str(atmSeq).replace("'","")[1:-1],G2mth.ValEsd(Dist,sig)))
1100
1101def BestPlane(PlaneData):
1102    'Needs a doc string'
1103
1104    def ShowBanner(name):
1105        G2fil.G2Print (80*'*')
1106        G2fil.G2Print ('   Best plane result for phase '+name)
1107        G2fil.G2Print (80*'*','\n')
1108
1109    ShowBanner(PlaneData['Name'])
1110
1111    Cell = PlaneData['Cell']
1112    Amat,Bmat = G2lat.cell2AB(Cell[:6])
1113    Atoms = PlaneData['Atoms']
1114    sumXYZ = np.zeros(3)
1115    XYZ = []
1116    Natoms = len(Atoms)
1117    for atom in Atoms:
1118        xyz = np.array(atom[3:6])
1119        XYZ.append(xyz)
1120        sumXYZ += xyz
1121    sumXYZ /= Natoms
1122    XYZ = np.array(XYZ)-sumXYZ
1123    XYZ = np.inner(Amat,XYZ).T
1124    Zmat = np.zeros((3,3))
1125    for i,xyz in enumerate(XYZ):
1126        Zmat += np.outer(xyz.T,xyz)
1127    G2fil.G2Print (' Selected atoms centered at %10.5f %10.5f %10.5f'%(sumXYZ[0],sumXYZ[1],sumXYZ[2]))
1128    Evec,Emat = nl.eig(Zmat)
1129    Evec = np.sqrt(Evec)/(Natoms-3)
1130    Order = np.argsort(Evec)
1131    XYZ = np.inner(XYZ,Emat.T).T
1132    XYZ = np.array([XYZ[Order[2]],XYZ[Order[1]],XYZ[Order[0]]]).T
1133    G2fil.G2Print (' Atoms in Cartesian best plane coordinates:')
1134    G2fil.G2Print (' Name         X         Y         Z')
1135    for i,xyz in enumerate(XYZ):
1136        G2fil.G2Print (' %6s%10.3f%10.3f%10.3f'%(Atoms[i][1].ljust(6),xyz[0],xyz[1],xyz[2]))
1137    G2fil.G2Print ('\n Best plane RMS X =%8.3f, Y =%8.3f, Z =%8.3f'%(Evec[Order[2]],Evec[Order[1]],Evec[Order[0]]))
1138
1139def main():
1140    'Called to run a refinement when this module is executed '
1141    starttime = time.time()
1142    arg = sys.argv
1143    if len(arg) > 1:
1144        GPXfile = arg[1]
1145        if not ospath.exists(GPXfile):
1146            G2fil.G2Print ('ERROR - '+GPXfile+" doesn't exist!")
1147            exit()
1148    else:
1149        G2fil.G2Print ('ERROR - missing filename')
1150        exit()
1151    # TODO: figure out if this is a sequential refinement and call SeqRefine(GPXfile,None)
1152    Refine(GPXfile,None)
1153    G2fil.G2Print("Done. Execution time {:.2f} sec.".format(time.time()-starttime))
1154
1155if __name__ == '__main__':
1156    GSASIIpath.InvokeDebugOpts()
1157    main()
Note: See TracBrowser for help on using the repository browser.