source: trunk/GSASIIstrMain.py @ 4876

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

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