source: trunk/GSASIIstrIO.py @ 4840

Last change on this file since 4840 was 4840, checked in by toby, 8 months ago

LeBail? issues: bug on zero cycles fixed; new LeBail? fit menu item; HessianLSQ: remove extra func eval on 0 cycles; fix noprint options; misc cleanups

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 173.0 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-03-05 19:54:34 +0000 (Fri, 05 Mar 2021) $
4# $Author: toby $
5# $Revision: 4840 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 4840 2021-03-05 19:54:34Z toby $
8########### SVN repository information ###################
9'''
10*GSASIIstrIO: structure I/O routines*
11-------------------------------------
12
13Contains routines for reading from GPX files and printing to the .LST file.
14Used for refinements and in G2scriptable.
15
16Should not contain any wxpython references as this should be able to be used
17in non-GUI settings.
18
19'''
20from __future__ import division, print_function
21import platform
22import os
23import os.path as ospath
24import time
25import math
26import copy
27if '2' in platform.python_version_tuple()[0]:
28    import cPickle
29else:
30    import pickle as cPickle
31import numpy as np
32import numpy.ma as ma
33import GSASIIpath
34GSASIIpath.SetVersionNumber("$Revision: 4840 $")
35import GSASIIElem as G2el
36import GSASIIlattice as G2lat
37import GSASIIspc as G2spc
38import GSASIIobj as G2obj
39import GSASIImapvars as G2mv
40import GSASIImath as G2mth
41import GSASIIfiles as G2fil
42import GSASIIpy3 as G2py3
43
44sind = lambda x: np.sin(x*np.pi/180.)
45cosd = lambda x: np.cos(x*np.pi/180.)
46tand = lambda x: np.tan(x*np.pi/180.)
47asind = lambda x: 180.*np.arcsin(x)/np.pi
48acosd = lambda x: 180.*np.arccos(x)/np.pi
49atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi
50   
51ateln2 = 8.0*math.log(2.0)
52
53#===============================================================================
54# Support for GPX file reading
55#===============================================================================
56def cPickleLoad(fp):
57    if '2' in platform.python_version_tuple()[0]:
58        return cPickle.load(fp)
59    else:
60       return cPickle.load(fp,encoding='latin-1')
61
62gpxIndex = {}; gpxNamelist = []; gpxSize = -1
63'''Global variables used in :func:`IndexGPX` to see if file has changed
64(gpxSize) and to index where to find each 1st-level tree item in the file.
65'''
66
67def GetFullGPX(GPXfile):
68    ''' Returns complete contents of GSASII gpx file.
69    Used in :func:`GSASIIscriptable.LoadDictFromProjFile`.
70
71    :param str GPXfile: full .gpx file name
72    :returns: Project,nameList, where
73
74      * Project (dict) is a representation of gpx file following the GSAS-II
75        tree structure for each item: key = tree name (e.g. 'Controls',
76        'Restraints', etc.), data is dict
77      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
78    '''
79    return IndexGPX(GPXfile,read=True)
80
81def IndexGPX(GPXfile,read=False):
82    '''Create an index to a GPX file, optionally the file into memory.
83    The byte size of the GPX file is saved. If this routine is called
84    again, and if this size does not change, indexing is not repeated
85    since it is assumed the file has not changed (this can be overriden
86    by setting read=True).
87
88    :param str GPXfile: full .gpx file name
89    :returns: Project,nameList if read=, where
90
91      * Project (dict) is a representation of gpx file following the GSAS-II
92        tree structure for each item: key = tree name (e.g. 'Controls',
93        'Restraints', etc.), data is dict
94      * nameList (list) has names of main tree entries & subentries used to reconstruct project file
95    '''
96    global gpxSize
97    if gpxSize == os.path.getsize(GPXfile) and not read:
98        return
99    global gpxIndex
100    gpxIndex = {}
101    global gpxNamelist
102    gpxNamelist = []
103    if GSASIIpath.GetConfigValue('debug'): print("DBG: Indexing GPX file")
104    gpxSize = os.path.getsize(GPXfile)
105    fp = open(GPXfile,'rb')
106    Project = {}
107    try:
108        while True:
109            pos = fp.tell()
110            data = cPickleLoad(fp)
111            datum = data[0]
112            gpxIndex[datum[0]] = pos
113            if read: Project[datum[0]] = {'data':datum[1]}
114            gpxNamelist.append([datum[0],])
115            for datus in data[1:]:
116                if read: Project[datum[0]][datus[0]] = datus[1]
117                gpxNamelist[-1].append(datus[0])
118        # print('project load successful')
119    except EOFError:
120        pass
121    except Exception as msg:
122        G2fil.G2Print('Read Error:',msg)
123        raise Exception("Error reading file "+str(GPXfile)+". This is not a GSAS-II .gpx file")
124    finally:
125        fp.close()
126    if read: return Project,gpxNamelist   
127   
128def GetControls(GPXfile):
129    ''' Returns dictionary of control items found in GSASII gpx file
130
131    :param str GPXfile: full .gpx file name
132    :return: dictionary of control items
133    '''
134    Controls = copy.deepcopy(G2obj.DefaultControls)
135    IndexGPX(GPXfile)
136    pos = gpxIndex.get('Controls')
137    if pos is None:
138        G2fil.G2Print('Warning: Controls not found in gpx file {}'.format(GPXfile))
139        return Controls
140    fp = open(GPXfile,'rb')
141    fp.seek(pos)
142    datum = cPickleLoad(fp)[0]
143    fp.close()
144    Controls.update(datum[1])
145    return Controls
146
147def GetConstraints(GPXfile):
148    '''Read the constraints from the GPX file and interpret them
149
150    called in :func:`ReadCheckConstraints`, :func:`GSASIIstrMain.Refine`
151    and :func:`GSASIIstrMain.SeqRefine`.
152    '''
153    IndexGPX(GPXfile)
154    fl = open(GPXfile,'rb')
155    pos = gpxIndex.get('Constraints')
156    if pos is None:
157        raise Exception("No constraints in GPX file")
158    fl.seek(pos)
159    datum = cPickleLoad(fl)[0]
160    fl.close()
161    constList = []
162    for item in datum[1]:
163        if item.startswith('_'): continue
164        constList += datum[1][item]
165    constDict,fixedList,ignored = ProcessConstraints(constList)
166    if ignored:
167        G2fil.G2Print ('Warning: {} Constraints were rejected. Was a constrained phase, histogram or atom deleted?'.format(ignored))
168    return constDict,fixedList
169   
170def ProcessConstraints(constList):
171    """Interpret the constraints in the constList input into a dictionary, etc.
172    All :class:`GSASIIobj.G2VarObj` objects are mapped to the appropriate
173    phase/hist/atoms based on the object internals (random Ids). If this can't be
174    done (if a phase has been deleted, etc.), the variable is ignored.
175    If the constraint cannot be used due to too many dropped variables,
176    it is counted as ignored.
177   
178    :param list constList: a list of lists where each item in the outer list
179      specifies a constraint of some form, as described in the :mod:`GSASIIobj`
180      :ref:`Constraint definition <Constraint_definitions_table>`.
181
182    :returns:  a tuple of (constDict,fixedList,ignored) where:
183     
184      * constDict (list of dicts) contains the constraint relationships
185      * fixedList (list) contains the fixed values for each type
186        of constraint.
187      * ignored (int) counts the number of invalid constraint items
188        (should always be zero!)
189    """
190    constDict = []
191    fixedList = []
192    ignored = 0
193    namedVarList = []
194    for item in constList:
195        if item[-1] == 'h':
196            # process a hold
197            fixedList.append('0')
198            var = str(item[0][1])
199            if '?' not in var:
200                constDict.append({var:0.0})
201            else:
202                ignored += 1
203        elif item[-1] == 'f':
204            # process a new variable
205            fixedList.append(None)
206            D = {}
207            varyFlag = item[-2]
208            varname = item[-3]
209            for term in item[:-3]:
210                var = str(term[1])
211                if '?' not in var:
212                    D[var] = term[0]
213            if len(D) > 1:
214                # add extra dict terms for input variable name and vary flag
215                if varname is not None:
216                    varname = str(varname) # in case this is a G2VarObj
217                    if varname.startswith(':'):
218                        D['_name'] = varname
219                    else:
220                        D['_name'] = '::nv-' + varname
221                    D['_name'] = G2obj.MakeUniqueLabel(D['_name'],namedVarList)
222                D['_vary'] = varyFlag == True # force to bool
223                constDict.append(D)
224            else:
225                ignored += 1
226            #constFlag[-1] = ['Vary']
227        elif item[-1] == 'c': 
228            # process a contraint relationship
229            D = {}
230            for term in item[:-3]:
231                var = str(term[1])
232                if '?' not in var:
233                    D[var] = term[0]
234            if len(D) >= 1:
235                fixedList.append(str(item[-3]))
236                constDict.append(D)
237            else:
238                ignored += 1
239        elif item[-1] == 'e':
240            # process an equivalence
241            firstmult = None
242            eqlist = []
243            for term in item[:-3]:
244                if term[0] == 0: term[0] = 1.0
245                var = str(term[1])
246                if '?' in var: continue
247                if firstmult is None:
248                    firstmult = term[0]
249                    firstvar = var
250                else:
251                    eqlist.append([var,firstmult/term[0]])
252            if len(eqlist) > 0:
253                G2mv.StoreEquivalence(firstvar,eqlist,False)
254            else:
255                ignored += 1
256        else:
257            ignored += 1
258    return constDict,fixedList,ignored
259
260def ReadCheckConstraints(GPXfile, seqHist=None):
261    '''Load constraints and related info and return any error or warning messages
262    This is done from the GPX file rather than the tree.
263
264    :param dict seqHist: defines a specific histogram to be loaded for a sequential
265       refinement, if None (default) all are loaded.
266    '''
267    # init constraints
268    G2mv.InitVars()   
269    # get variables
270    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
271    if not Phases:
272        return 'Error: No phases or no histograms in phases!',''
273    if not Histograms:
274        return 'Error: no diffraction data',''
275    constrDict,fixedList = GetConstraints(GPXfile) # load user constraints before internally generated ones
276    if seqHist: Histograms = {seqHist:Histograms[seqHist]}  # sequential fit: only need one histogram
277    rigidbodyDict = GetRigidBodies(GPXfile)
278    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
279    rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
280    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave = \
281        GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False) # generates atom symmetry constraints
282    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False,resetRefList=False)
283    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
284    varyList = rbVary+phaseVary+hapVary+histVary
285    msg = G2mv.EvaluateMultipliers(constrDict,phaseDict,hapDict,histDict)
286    if msg:
287        return 'Unable to interpret multiplier(s): '+msg,''
288    errmsg, warnmsg = G2mv.CheckConstraints(varyList,constrDict,fixedList)
289    if errmsg:
290        # print some diagnostic info on the constraints
291        G2fil.G2Print('Error in constraints:\n'+errmsg+
292              '\nRefinement not possible due to conflict in constraints, see below:\n')
293        G2fil.G2Print(G2mv.VarRemapShow(varyList,True),mode='error')
294    return errmsg, warnmsg
295   
296def makeTwinFrConstr(Phases,Histograms,hapVary):
297    TwConstr = []
298    TwFixed = []
299    for Phase in Phases:
300        pId = Phases[Phase]['pId']
301        for Histogram in Phases[Phase]['Histograms']:
302            try:
303                hId = Histograms[Histogram]['hId']
304                phfx = '%d:%d:'%(pId,hId)
305                if phfx+'TwinFr:0' in hapVary:
306                    TwFixed.append('1.0')     #constraint value
307                    nTwin = len(Phases[Phase]['Histograms'][Histogram]['Twins'])
308                    TwConstr.append({phfx+'TwinFr:'+str(i):'1.0' for i in range(nTwin)})
309            except KeyError:    #unused histograms?
310                pass
311    return TwConstr,TwFixed   
312   
313def GetRestraints(GPXfile):
314    '''Read the restraints from the GPX file.
315    Throws an exception if not found in the .GPX file
316    '''
317    IndexGPX(GPXfile)
318    fl = open(GPXfile,'rb')
319    pos = gpxIndex.get('Restraints')
320    if pos is None:
321        raise Exception("No Restraints in GPX file")
322    fl.seek(pos)
323    datum = cPickleLoad(fl)[0]
324    fl.close()
325    return datum[1]
326   
327def GetRigidBodies(GPXfile):
328    '''Read the rigid body models from the GPX file
329    '''
330    IndexGPX(GPXfile)
331    fl = open(GPXfile,'rb')
332    pos = gpxIndex.get('Rigid bodies')
333    if pos is None:
334        raise Exception("No Rigid bodies in GPX file")
335    fl.seek(pos)
336    datum = cPickleLoad(fl)[0]
337    fl.close()
338    return datum[1]
339       
340def GetFprime(controlDict,Histograms):
341    'Needs a doc string'
342    FFtables = controlDict['FFtables']
343    if not FFtables:
344        return
345    histoList = list(Histograms.keys())
346    histoList.sort()
347    for histogram in histoList:
348        if histogram[:4] in ['PWDR','HKLF']:
349            Histogram = Histograms[histogram]
350            hId = Histogram['hId']
351            hfx = ':%d:'%(hId)
352            if 'X' in controlDict[hfx+'histType']:
353                keV = controlDict[hfx+'keV']
354                for El in FFtables:
355                    Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
356                    FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
357                    FFtables[El][hfx+'FP'] = FP
358                    FFtables[El][hfx+'FPP'] = FPP
359                   
360def PrintFprime(FFtables,pfx,pFile):
361    pFile.write('\n Resonant form factors:(ref: D.T. Cromer & D.A. Liberman (1981), Acta Cryst. A37, 267-268.)\n')
362    Elstr = ' Element:'
363    FPstr = " f'     :"
364    FPPstr = ' f"     :'
365    for El in FFtables:
366        Elstr += ' %8s'%(El)
367        FPstr += ' %8.3f'%(FFtables[El][pfx+'FP'])
368        FPPstr += ' %8.3f'%(FFtables[El][pfx+'FPP'])
369    pFile.write(Elstr+'\n')
370    pFile.write(FPstr+'\n')
371    pFile.write(FPPstr+'\n')
372   
373def PrintBlength(BLtables,wave,pFile):
374    pFile.write('\n Resonant neutron scattering lengths:\n')
375    Elstr = ' Element:'
376    FPstr = " b'     :"
377    FPPstr = ' b"     :'
378    for El in BLtables:
379        BP,BPP = G2el.BlenResCW([El,],BLtables,wave)
380        Elstr += ' %8s'%(El)
381        FPstr += ' %8.3f'%(BP)
382        FPPstr += ' %8.3f'%(BPP)
383    pFile.write(Elstr+'\n')
384    pFile.write(FPstr+'\n')
385    pFile.write(FPPstr+'\n')
386
387def PrintISOmodes(pFile,Phases,parmDict,sigDict):
388    '''Prints the values for the ISODISTORT modes into the project's
389    .lst file after a refinement.
390    '''
391    for phase in Phases:
392        if 'ISODISTORT' not in Phases[phase]: continue
393        data = Phases[phase]
394        ISO = data['ISODISTORT']
395        if 'G2VarList' in ISO:
396            deltaList = []
397            for gv,Ilbl in zip(ISO['G2VarList'],ISO['IsoVarList']):
398                dvar = gv.varname()
399                var = dvar.replace('::dA','::A')
400                albl = Ilbl[:Ilbl.rfind('_')]
401                v = Ilbl[Ilbl.rfind('_')+1:]
402                pval = ISO['ParentStructure'][albl][['dx','dy','dz'].index(v)]
403                if var in parmDict:
404                    cval = parmDict[var]
405                else:
406                    print('PrintISOmodes Error: Atom not found',"No value found for parameter "+str(var))
407                    return
408                deltaList.append(cval-pval)
409            modeVals = np.inner(ISO['Var2ModeMatrix'],deltaList)
410           
411        if 'G2OccVarList' in ISO:
412            deltaOccList = []
413            for gv,Ilbl in zip(ISO['G2OccVarList'],ISO['OccVarList']):
414                var = gv.varname()
415                albl = Ilbl[:Ilbl.rfind('_')]
416                pval = ISO['BaseOcc'][albl]
417                if var in parmDict:
418                    cval = parmDict[var]
419                else:
420                    print('PrintISOmodes Error: Atom not found',"No value found for parameter "+str(var))
421                    return
422                deltaOccList.append(cval-pval)
423            modeOccVals = np.inner(ISO['Var2OccMatrix'],deltaOccList)
424           
425        if 'G2VarList' in ISO:
426            pFile.write('\n ISODISTORT Displacive Modes for phase {}\n'.format(data['General'].get('Name','')))
427            l = str(max([len(i) for i in ISO['IsoModeList']])+3)
428            fmt = '  {:'+l+'}{}'
429            for var,val,norm,G2mode in zip(
430                    ISO['IsoModeList'],modeVals,ISO['NormList'],ISO['G2ModeList'] ):
431                try:
432                    value = G2py3.FormatSigFigs(val/norm)
433                    if str(G2mode) in sigDict:
434                        value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm)
435                except TypeError:
436                    value = '?'
437                pFile.write(fmt.format(var,value)+'\n')
438               
439        if 'G2OccVarList' in ISO:
440            pFile.write('\n ISODISTORT Occupancy Modes for phase {}\n'.format(data['General'].get('Name','')))
441            l = str(max([len(i) for i in ISO['OccModeList']])+3)
442            fmt = '  {:'+l+'}{}'
443            for var,val,norm,G2mode in zip(
444                    ISO['OccModeList'],modeOccVals,ISO['OccNormList'],ISO['G2OccModeList'] ):
445                try:
446                    value = G2py3.FormatSigFigs(val/norm)
447                    if str(G2mode) in sigDict:
448                        value = G2mth.ValEsd(val/norm,sigDict[str(G2mode)]/norm)
449                except TypeError:
450                    value = '?'
451                pFile.write(fmt.format(var,value)+'\n')
452   
453def GetPhaseNames(GPXfile):
454    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
455
456    :param str GPXfile: full .gpx file name
457    :return: list of phase names
458    '''
459    IndexGPX(GPXfile)
460    fl = open(GPXfile,'rb')
461    pos = gpxIndex.get('Phases')
462    if pos is None:
463        raise Exception("No Phases in GPX file")
464    fl.seek(pos)
465    data = cPickleLoad(fl)
466    fl.close()
467    return [datus[0] for datus in data[1:]]
468
469def GetAllPhaseData(GPXfile,PhaseName):
470    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
471
472    :param str GPXfile: full .gpx file name
473    :param str PhaseName: phase name
474    :return: phase dictionary or None if PhaseName is not present
475    '''       
476    IndexGPX(GPXfile)
477    fl = open(GPXfile,'rb')
478    pos = gpxIndex.get('Phases')
479    if pos is None:
480        raise Exception("No Phases in GPX file")
481    fl.seek(pos)
482    data = cPickleLoad(fl)
483    fl.close()
484
485    for datus in data[1:]:
486        if datus[0] == PhaseName:
487            return datus[1]
488   
489def GetHistograms(GPXfile,hNames):
490    """ Returns a dictionary of histograms found in GSASII gpx file
491
492    :param str GPXfile: full .gpx file name
493    :param str hNames: list of histogram names
494    :return: dictionary of histograms (types = PWDR & HKLF)
495
496    """
497    IndexGPX(GPXfile)
498    fl = open(GPXfile,'rb')
499    Histograms = {}
500    for hist in hNames:
501        pos = gpxIndex.get(hist)
502        if pos is None:
503            raise Exception("Histogram {} not found in GPX file".format(hist))
504        fl.seek(pos)
505        data = cPickleLoad(fl)
506        datum = data[0]
507        if 'PWDR' in hist[:4]:
508            PWDRdata = {}
509            PWDRdata.update(datum[1][0])        #weight factor
510            PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1]))          #masked powder data arrays/clear previous masks
511            PWDRdata[data[2][0]] = data[2][1]       #Limits & excluded regions (if any)
512            PWDRdata[data[3][0]] = data[3][1]       #Background
513            PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
514            PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
515            try:
516               PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
517            except IndexError:
518                PWDRdata['Reflection Lists'] = {}
519            PWDRdata['Residuals'] = {}
520            Histograms[hist] = PWDRdata
521        elif 'HKLF' in hist[:4]:
522            HKLFdata = {}
523            HKLFdata.update(datum[1][0])        #weight factor
524#patch
525            if 'list' in str(type(datum[1][1])):
526            #if isinstance(datum[1][1],list):
527                RefData = {'RefList':[],'FF':{}}
528                for ref in datum[1][1]:
529                    RefData['RefList'].append(ref[:11]+[ref[13],])
530                RefData['RefList'] = np.array(RefData['RefList'])
531                datum[1][1] = RefData
532#end patch
533            datum[1][1]['FF'] = {}
534            HKLFdata['Data'] = datum[1][1]
535            HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
536            HKLFdata['Reflection Lists'] = None
537            HKLFdata['Residuals'] = {}
538            Histograms[hist] = HKLFdata           
539    fl.close()
540    return Histograms
541   
542def GetHistogramNames(GPXfile,hTypes):
543    """ Returns a list of histogram names found in a GSAS-II .gpx file that
544    match specifed histogram types. Names are returned in the order they
545    appear in the file.
546
547    :param str GPXfile: full .gpx file name
548    :param str hTypes: list of histogram types
549    :return: list of histogram names (types = PWDR & HKLF)
550
551    """
552    IndexGPX(GPXfile)
553    return [n[0] for n in gpxNamelist if n[0][:4] in hTypes]
554   
555def GetUsedHistogramsAndPhases(GPXfile):
556    ''' Returns all histograms that are found in any phase
557    and any phase that uses a histogram. This also
558    assigns numbers to used phases and histograms by the
559    order they appear in the file.
560
561    :param str GPXfile: full .gpx file name
562    :returns: (Histograms,Phases)
563
564     * Histograms = dictionary of histograms as {name:data,...}
565     * Phases = dictionary of phases that use histograms
566
567    '''
568    phaseNames = GetPhaseNames(GPXfile)
569    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
570    allHistograms = GetHistograms(GPXfile,histoList)
571    phaseData = {}
572    for name in phaseNames: 
573        phaseData[name] =  GetAllPhaseData(GPXfile,name)
574    Histograms = {}
575    Phases = {}
576    for phase in phaseData:
577        Phase = phaseData[phase]
578        if Phase['General']['Type'] == 'faulted': continue      #don't use faulted phases!
579        if Phase['Histograms']:
580            for hist in Phase['Histograms']:
581                if 'Use' not in Phase['Histograms'][hist]:      #patch
582                    Phase['Histograms'][hist]['Use'] = True
583                if Phase['Histograms'][hist]['Use'] and phase not in Phases:
584                    pId = phaseNames.index(phase)
585                    Phase['pId'] = pId
586                    Phases[phase] = Phase
587                if hist not in Histograms and Phase['Histograms'][hist]['Use']:
588                    try:
589                        Histograms[hist] = allHistograms[hist]
590                        hId = histoList.index(hist)
591                        Histograms[hist]['hId'] = hId
592                    except KeyError: # would happen if a referenced histogram were
593                        # renamed or deleted
594                        G2fil.G2Print('Warning: For phase "'+phase+
595                              '" unresolved reference to histogram "'+hist+'"')
596    # load the fix background info into the histograms
597    for hist in Histograms:
598        if 'Background' not in Histograms[hist]: continue
599        fixedBkg = Histograms[hist]['Background'][1].get('background PWDR')
600        if fixedBkg:
601            if not fixedBkg[0]: continue
602            # patch: add refinement flag, if needed
603            if len(fixedBkg) == 2: fixedBkg += [False]
604            h = Histograms[hist]['Background'][1]
605            try:
606                Limits = Histograms[hist]['Limits'][1]
607                x = Histograms[hist]['Data'][0]
608                xB = np.searchsorted(x,Limits[0])
609                xF = np.searchsorted(x,Limits[1])+1
610                h['fixback'] = allHistograms[fixedBkg[0]]['Data'][1][xB:xF]
611            except KeyError: # would happen if a referenced histogram were renamed or deleted
612                G2fil.G2Print('Warning: For hist "{}" unresolved background reference to hist "{}"'
613                          .format(hist,fixedBkg[0]))
614    G2obj.IndexAllIds(Histograms=Histograms,Phases=Phases)
615    return Histograms,Phases
616   
617def getBackupName(GPXfile,makeBack):
618    '''
619    Get the name for the backup .gpx file name
620   
621    :param str GPXfile: full .gpx file name
622    :param bool makeBack: if True the name of a new file is returned, if
623      False the name of the last file that exists is returned
624    :returns: the name of a backup file
625   
626    '''
627    GPXpath,GPXname = ospath.split(GPXfile)
628    if GPXpath == '': GPXpath = '.'
629    Name = ospath.splitext(GPXname)[0]
630    files = os.listdir(GPXpath)
631    last = 0
632    for name in files:
633        name = name.split('.')
634        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
635            if makeBack:
636                last = max(last,int(name[1].strip('bak'))+1)
637            else:
638                last = max(last,int(name[1].strip('bak')))
639    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
640    return GPXback   
641       
642def GPXBackup(GPXfile,makeBack=True):
643    '''
644    makes a backup of the specified .gpx file
645   
646    :param str GPXfile: full .gpx file name
647    :param bool makeBack: if True (default), the backup is written to
648      a new file; if False, the last backup is overwritten
649    :returns: the name of the backup file that was written
650    '''
651    import distutils.file_util as dfu
652    GPXback = getBackupName(GPXfile,makeBack)
653    tries = 0
654    while True:
655        try:
656            dfu.copy_file(GPXfile,GPXback)
657            break
658        except:
659            tries += 1
660            if tries > 10:
661                return GPXfile  #failed!
662            time.sleep(1)           #just wait a second!         
663    return GPXback
664
665def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,parmFrozenList,makeBack=True):
666    ''' Updates gpxfile from all histograms that are found in any phase
667    and any phase that used a histogram. Also updates rigid body definitions.
668    This is used for non-sequential fits, but not for sequential fitting.
669
670    :param str GPXfile: full .gpx file name
671    :param dict Histograms: dictionary of histograms as {name:data,...}
672    :param dict Phases: dictionary of phases that use histograms
673    :param dict RigidBodies: dictionary of rigid bodies
674    :param dict CovData: dictionary of refined variables, varyList, & covariance matrix
675    :param list parmFrozenList: list of parameters (as str) that are frozen
676      due to limits; converted to :class:`GSASIIobj.G2VarObj` objects.
677    :param bool makeBack: True if new backup of .gpx file is to be made; else
678      use the last one made
679    '''
680                       
681    import distutils.file_util as dfu
682    GPXback = GPXBackup(GPXfile,makeBack)
683    G2fil.G2Print ('Read from file:'+GPXback)
684    G2fil.G2Print ('Save to file  :'+GPXfile)
685    infile = open(GPXback,'rb')
686    outfile = open(GPXfile,'wb')
687    while True:
688        try:
689            data = cPickleLoad(infile)
690        except EOFError:
691            break
692        datum = data[0]
693#        print 'read: ',datum[0]
694        if datum[0] == 'Phases':
695            for iphase in range(len(data)):
696                if data[iphase][0] in Phases:
697                    phaseName = data[iphase][0]
698                    data[iphase][1].update(Phases[phaseName])
699        elif datum[0] == 'Covariance':
700            data[0][1] = CovData
701        elif datum[0] == 'Rigid bodies':
702            data[0][1] = RigidBodies
703        elif datum[0] == 'Controls':
704            Controls = data[0][1]
705            if 'parmFrozen' not in Controls:
706                Controls['parmFrozen'] = {}
707            Controls['parmFrozen']['FrozenList'] = [
708                    i if type(i) is G2obj.G2VarObj
709                    else G2obj.G2VarObj(i)
710                    for i in parmFrozenList]
711        try:
712            histogram = Histograms[datum[0]]
713#            print 'found ',datum[0]
714            data[0][1][1] = histogram['Data']
715            data[0][1][0].update(histogram['Residuals'])
716            for datus in data[1:]:
717#                print '    read: ',datus[0]
718                if datus[0] in ['Instrument Parameters','Sample Parameters','Reflection Lists']:
719                    datus[1] = histogram[datus[0]]
720                if datus[0] == 'Background': # remove fixed background from file
721                    d1 = {key:histogram['Background'][1][key]
722                              for key in histogram['Background'][1]
723                              if not key.startswith('_fixed')}               
724                    datus[1] = copy.deepcopy(histogram['Background'])
725                    datus[1][1] = d1
726        except KeyError:
727            pass
728        try:                       
729            cPickle.dump(data,outfile,1)
730        except AttributeError:
731            G2fil.G2Print ('ERROR - bad data in least squares result')
732            infile.close()
733            outfile.close()
734            dfu.copy_file(GPXback,GPXfile)
735            G2fil.G2Print ('GPX file save failed - old version retained',mode='error')
736            return
737       
738    infile.close()
739    outfile.close()
740           
741    G2fil.G2Print ('GPX file save successful')
742   
743def GetSeqResult(GPXfile):
744    '''
745    Returns the sequential results table information from a GPX file.
746    Called at the beginning of :meth:`GSASIIstrMain.SeqRefine`
747   
748    :param str GPXfile: full .gpx file name
749    :returns: a dict containing the sequential results table
750    '''
751    IndexGPX(GPXfile)
752    pos = gpxIndex.get('Sequential results')
753    if pos is None:
754        return {}
755    fl = open(GPXfile,'rb')
756    fl.seek(pos)
757    datum = cPickleLoad(fl)[0]
758    fl.close()
759    return datum[1]
760   
761def SetupSeqSavePhases(GPXfile):
762    '''Initialize the files used to save intermediate results from
763    sequential fits.
764    '''
765    IndexGPX(GPXfile)
766    # load initial Phase results from GPX
767    fl = open(GPXfile,'rb')
768    pos = gpxIndex.get('Phases')
769    if pos is None:
770        raise Exception("No Phases in GPX file")
771    fl.seek(pos)
772    data = cPickleLoad(fl)
773    fl.close()
774    # create GPX-like file to store latest Phase info; init with start vals
775    GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
776    fp = open(GPXphase,'wb')
777    cPickle.dump(data,fp,1)
778    fp.close()
779    # create empty file for histogram info
780    GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
781    fp = open(GPXhist,'wb')
782    fp.close()
783
784def SaveUpdatedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,parmFrozen):
785    '''
786    Save phase and histogram information into "pseudo-gpx" files. The phase
787    information is overwritten each time this is called, but histogram information is
788    appended after each sequential step.
789
790    :param str GPXfile: full .gpx file name
791    :param dict Histograms: dictionary of histograms as {name:data,...}
792    :param dict Phases: dictionary of phases that use histograms
793    :param dict RigidBodies: dictionary of rigid bodies
794    :param dict CovData: dictionary of refined variables, varyList, & covariance matrix
795    :param dict parmFrozen: dict with frozen parameters for all phases
796      and histograms (specified as str values)
797    '''
798                       
799    import distutils.file_util as dfu
800   
801    GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
802    fp = open(GPXphase,'rb')
803    data = cPickleLoad(fp) # first block in file should be Phases
804    if data[0][0] != 'Phases':
805        raise Exception('Unexpected block in {} file. How did this happen?'
806                            .format(GPXphase))
807    fp.close()
808    # update previous phase info
809    for datum in data[1:]: 
810        if datum[0] in Phases:
811            datum[1].update(Phases[datum[0]])
812    # save latest Phase/refinement info
813    fp = open(GPXphase,'wb')
814    cPickle.dump(data,fp,1)
815    cPickle.dump([['Covariance',CovData]],fp,1)
816    cPickle.dump([['Rigid bodies',RigidBodies]],fp,1)
817    cPickle.dump([['parmFrozen',parmFrozen]],fp,1)
818    fp.close()
819    # create an entry that looks like a PWDR tree item
820    for key in Histograms:
821        if key.startswith('PWDR '):
822            break
823    else:
824        raise Exception('No PWDR entry in Histogram dict!')
825    histname = key
826    hist = copy.deepcopy(Histograms[key])
827    xfer_dict = {'Index Peak List': [[], []],
828                   'Comments': [],
829                   'Unit Cells List': [],
830                   'Peak List': {'peaks': [], 'sigDict': {}},
831                   }
832    histData = hist['Data']
833    del hist['Data']
834    for key in ('Limits','Background','Instrument Parameters',
835                    'Sample Parameters','Reflection Lists'):
836        xfer_dict[key] = hist[key]
837        if key == 'Background':  # remove fixed background from file
838            xfer_dict['Background'][1] = {k:hist['Background'][1][k]
839                      for k in hist['Background'][1]
840                      if not k.startswith('_fixed')}               
841        del hist[key]
842    # xform into a gpx-type entry
843    data = []
844    data.append([histname,[hist,histData,histname]])       
845    for key in ['Comments','Limits','Background','Instrument Parameters',
846             'Sample Parameters','Peak List','Index Peak List',
847             'Unit Cells List','Reflection Lists']:
848        data.append([key,xfer_dict[key]])
849    # append histogram to histogram info
850    GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
851    fp = open(GPXhist,'ab')
852    cPickle.dump(data,fp,1)
853    fp.close()
854    return
855   
856def SetSeqResult(GPXfile,Histograms,SeqResult):
857    '''
858    Places the sequential results information into a GPX file
859    after a refinement has been completed.
860    Called at the end of :meth:`GSASIIstrMain.SeqRefine`
861
862    :param str GPXfile: full .gpx file name
863    '''
864    GPXback = GPXBackup(GPXfile)
865    G2fil.G2Print ('Read from file:'+GPXback)
866    G2fil.G2Print ('Save to file  :'+GPXfile)
867    GPXphase = os.path.splitext(GPXfile)[0]+'.seqPhase'
868    fp = open(GPXphase,'rb')
869    data = cPickleLoad(fp) # first block in file should be Phases
870    if data[0][0] != 'Phases':
871        raise Exception('Unexpected block in {} file. How did this happen?'.format(GPXphase))
872    Phases = {}
873    for name,vals in data[1:]:
874        Phases[name] = vals       
875    name,CovData = cPickleLoad(fp)[0] # 2nd block in file should be Covariance
876    name,RigidBodies = cPickleLoad(fp)[0] # 3rd block in file should be Rigid Bodies
877    name,parmFrozenDict = cPickleLoad(fp)[0] # 4th block in file should be frozen parameters
878    fp.close()
879    GPXhist = os.path.splitext(GPXfile)[0]+'.seqHist'
880    hist = open(GPXhist,'rb')
881    # build an index to the GPXhist file
882    histIndex = {}
883    while True:
884        loc = hist.tell()
885        try:
886            datum = cPickleLoad(hist)[0]
887        except EOFError:
888            break
889        histIndex[datum[0]] = loc
890
891    infile = open(GPXback,'rb')
892    outfile = open(GPXfile,'wb')
893    while True:
894        try:
895            data = cPickleLoad(infile)
896        except EOFError:
897            break
898        datum = data[0]
899        if datum[0] == 'Sequential results':
900            data[0][1] = SeqResult
901        elif datum[0] == 'Phases':
902            for pdata in data[1:]:
903                if pdata[0] in Phases:
904                    pdata[1].update(Phases[pdata[0]])
905        elif datum[0] == 'Covariance':
906            data[0][1] = CovData
907        elif datum[0] == 'Rigid bodies':
908            data[0][1] = RigidBodies
909        elif datum[0] == 'Controls': # reset the Copy Next flag after a sequential fit
910            Controls = data[0][1]
911            Controls['Copy2Next'] = False
912            for key in parmFrozenDict:
913                Controls['parmFrozen'][key] = [
914                    i if type(i) is G2obj.G2VarObj
915                    else G2obj.G2VarObj(i)
916                    for i in parmFrozenDict[key]]
917        elif datum[0] in histIndex:
918            hist.seek(histIndex[datum[0]])
919            hdata = cPickleLoad(hist)
920            if data[0][0] != hdata[0][0]:
921                G2fil.G2Print('Error! Updating {} with {}'.format(data[0][0],hdata[0][0]))
922            data[0] = hdata[0]
923            xferItems = ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']
924            hItems = {name:j+1 for j,(name,val) in enumerate(hdata[1:]) if name in xferItems}
925            for j,(name,val) in enumerate(data[1:]):
926                if name not in xferItems: continue
927                data[j+1][1] = hdata[hItems[name]][1]
928        cPickle.dump(data,outfile,1)
929    hist.close()
930    infile.close()
931    outfile.close()
932    # clean up tmp files
933    try:
934        os.remove(GPXphase)
935    except:
936        G2fil.G2Print('Warning: unable to delete {}'.format(GPXphase))
937    try:
938        os.remove(GPXhist)
939    except:
940        G2fil.G2Print('Warning: unable to delete {}'.format(GPXhist))
941    G2fil.G2Print ('GPX file merge completed')
942
943#==============================================================================
944# Refinement routines
945#==============================================================================
946def ShowBanner(pFile=None):
947    'Print authorship, copyright and citation notice'
948    pFile.write(80*'*'+'\n')
949    pFile.write('   General Structure Analysis System-II Crystal Structure Refinement\n')
950    pFile.write('              by Robert B. Von Dreele & Brian H. Toby\n')
951    pFile.write('                Argonne National Laboratory(C), 2010\n')
952    pFile.write(' This product includes software developed by the UChicago Argonne, LLC,\n')
953    pFile.write('            as Operator of Argonne National Laboratory.\n')
954    pFile.write('                          Please cite:\n')
955    pFile.write('   B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)\n')
956
957    pFile.write(80*'*'+'\n')
958
959def ShowControls(Controls,pFile=None,SeqRef=False,preFrozenCount=0):
960    'Print controls information'
961    pFile.write(' Least squares controls:\n')
962    pFile.write(' Refinement type: %s\n'%Controls['deriv type'])
963    if 'Hessian' in Controls['deriv type']:
964        pFile.write(' Maximum number of cycles: %d\n'%Controls['max cyc'])
965    else:
966        pFile.write(' Minimum delta-M/M for convergence: %.2g\n'%(Controls['min dM/M']))
967    pFile.write(' Regularize hydrogens (if any): %s\n'%Controls.get('HatomFix',False))
968    pFile.write(' Initial shift factor: %.3f\n'%(Controls['shift factor']))
969    if SeqRef:
970        pFile.write(' Sequential refinement controls:\n')
971        pFile.write(' Copy of histogram results to next: %s\n'%(Controls['Copy2Next']))
972        pFile.write(' Process histograms in reverse order: %s\n'%(Controls['Reverse Seq']))
973    if preFrozenCount:
974        pFile.write('\n Starting refinement with {} Frozen variables\n\n'.format(preFrozenCount))
975   
976def GetPawleyConstr(SGLaue,PawleyRef,im,pawleyVary):
977    'needs a doc string'
978#    if SGLaue in ['-1','2/m','mmm']:
979#        return                      #no Pawley symmetry required constraints
980    eqvDict = {}
981    for i,varyI in enumerate(pawleyVary):
982        eqvDict[varyI] = []
983        refI = int(varyI.split(':')[-1])
984        ih,ik,il = PawleyRef[refI][:3]
985        dspI = PawleyRef[refI][4+im]
986        for varyJ in pawleyVary[i+1:]:
987            refJ = int(varyJ.split(':')[-1])
988            jh,jk,jl = PawleyRef[refJ][:3]
989            dspJ = PawleyRef[refJ][4+im]
990            if SGLaue in ['4/m','4/mmm']:
991                isum = ih**2+ik**2
992                jsum = jh**2+jk**2
993                if abs(il) == abs(jl) and isum == jsum:
994                    eqvDict[varyI].append(varyJ) 
995            elif SGLaue in ['3R','3mR']:
996                isum = ih**2+ik**2+il**2
997                jsum = jh**2+jk**2+jl**2
998                isum2 = ih*ik+ih*il+ik*il
999                jsum2 = jh*jk+jh*jl+jk*jl
1000                if isum == jsum and isum2 == jsum2:
1001                    eqvDict[varyI].append(varyJ) 
1002            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
1003                isum = ih**2+ik**2+ih*ik
1004                jsum = jh**2+jk**2+jh*jk
1005                if abs(il) == abs(jl) and isum == jsum:
1006                    eqvDict[varyI].append(varyJ) 
1007            elif SGLaue in ['m3','m3m']:
1008                isum = ih**2+ik**2+il**2
1009                jsum = jh**2+jk**2+jl**2
1010                if isum == jsum:
1011                    eqvDict[varyI].append(varyJ)
1012            elif abs(dspI-dspJ)/dspI < 1.e-4:
1013                eqvDict[varyI].append(varyJ) 
1014    for item in pawleyVary:
1015        if eqvDict[item]:
1016            for item2 in pawleyVary:
1017                if item2 in eqvDict[item]:
1018                    eqvDict[item2] = []
1019            G2mv.StoreEquivalence(item,eqvDict[item])
1020                   
1021def cellVary(pfx,SGData): 
1022    '''Creates equivalences for a phase based on the Laue class.
1023    Returns a list of A tensor terms that are non-zero.
1024    '''
1025    if SGData['SGLaue'] in ['-1',]:
1026        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1027    elif SGData['SGLaue'] in ['2/m',]:
1028        if SGData['SGUniq'] == 'a':
1029            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
1030        elif SGData['SGUniq'] == 'b':
1031            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
1032        else:
1033            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
1034    elif SGData['SGLaue'] in ['mmm',]:
1035        return [pfx+'A0',pfx+'A1',pfx+'A2']
1036    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1037        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',))
1038        return [pfx+'A0',pfx+'A1',pfx+'A2']
1039    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1040        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A3',))
1041        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
1042    elif SGData['SGLaue'] in ['3R', '3mR']:
1043        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',))
1044        G2mv.StoreEquivalence(pfx+'A3',(pfx+'A4',pfx+'A5',))
1045        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']                       
1046    elif SGData['SGLaue'] in ['m3m','m3']:
1047        G2mv.StoreEquivalence(pfx+'A0',(pfx+'A1',pfx+'A2',))
1048        return [pfx+'A0',pfx+'A1',pfx+'A2']
1049       
1050def modVary(pfx,SSGData):
1051    vary = []
1052    for i,item in enumerate(SSGData['modSymb']):
1053        if item in ['a','b','g']:
1054            vary.append(pfx+'mV%d'%(i))
1055    return vary
1056       
1057################################################################################
1058##### Rigid Body Models and not General.get('doPawley')
1059################################################################################
1060       
1061def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
1062    'Get Rigid body info from tree entry and print it to .LST file'
1063   
1064    def PrintResRBModel(RBModel):
1065        pFile.write('Residue RB name: %s no.atoms: %d, No. times used: %d\n'%
1066            (RBModel['RBname'],len(RBModel['rbTypes']),RBModel['useCount']))
1067        for i in WriteResRBModel(RBModel):
1068            pFile.write(i)
1069       
1070    def PrintVecRBModel(RBModel):
1071        pFile.write('Vector RB name: %s no.atoms: %d No. times used: %d\n'%
1072            (RBModel['RBname'],len(RBModel['rbTypes']),RBModel['useCount']))
1073        for i in WriteVecRBModel(RBModel):
1074            pFile.write(i)
1075        pFile.write('Orientation defined by: atom %s -> atom %s & atom %s -> atom %s\n'%
1076            (RBModel['rbRef'][0],RBModel['rbRef'][1],RBModel['rbRef'][0],RBModel['rbRef'][2]))
1077           
1078    if Print and pFile is None: raise Exception("specify pFile or Print=False")
1079    rbVary = []
1080    rbDict = {}
1081    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
1082    if len(rbIds['Vector']):
1083        for irb,item in enumerate(rbIds['Vector']):
1084            if rigidbodyDict['Vector'][item]['useCount']:
1085                RBmags = rigidbodyDict['Vector'][item]['VectMag']
1086                RBrefs = rigidbodyDict['Vector'][item]['VectRef']
1087                for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
1088                    pid = '::RBV;'+str(i)+':'+str(irb)
1089                    rbDict[pid] = mag
1090                    if ref:
1091                        rbVary.append(pid)
1092                if Print:
1093                    pFile.write('\nVector rigid body model:\n')
1094                    PrintVecRBModel(rigidbodyDict['Vector'][item])
1095    if len(rbIds['Residue']):
1096        for item in rbIds['Residue']:
1097            if rigidbodyDict['Residue'][item]['useCount']:
1098                if Print:
1099                    pFile.write('\nResidue rigid body model:\n')
1100                    PrintResRBModel(rigidbodyDict['Residue'][item])
1101    return rbVary,rbDict
1102   
1103def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
1104    'needs a doc string'
1105   
1106    def PrintRBVectandSig(VectRB,VectSig):
1107        pFile.write('\n Rigid body vector magnitudes for %s:\n'%VectRB['RBname'])
1108        namstr = '  names :'
1109        valstr = '  values:'
1110        sigstr = '  esds  :'
1111        for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
1112            namstr += '%12s'%('Vect '+str(i))
1113            valstr += '%12.4f'%(val)
1114            if sig:
1115                sigstr += '%12.4f'%(sig)
1116            else:
1117                sigstr += 12*' '
1118        pFile.write(namstr+'\n')
1119        pFile.write(valstr+'\n')
1120        pFile.write(sigstr+'\n')       
1121       
1122    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
1123    if not RBIds['Vector']:
1124        return
1125    for irb,item in enumerate(RBIds['Vector']):
1126        if rigidbodyDict['Vector'][item]['useCount']:
1127            VectSig = []
1128            RBmags = rigidbodyDict['Vector'][item]['VectMag']
1129            for i,mag in enumerate(RBmags):
1130                name = '::RBV;'+str(i)+':'+str(irb)
1131                if name in sigDict:
1132                    VectSig.append(sigDict[name])
1133            PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)   
1134       
1135################################################################################
1136##### Phase data
1137################################################################################                   
1138def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None,seqRef=False):
1139    '''Setup the phase information for a structural refinement, used for
1140    regular and sequential refinements, optionally printing information
1141    to the .lst file (if Print is True)
1142    '''
1143           
1144    def PrintFFtable(FFtable):
1145        pFile.write('\n X-ray scattering factors:\n')
1146        pFile.write('   Symbol     fa                                      fb                                      fc\n')
1147        pFile.write(99*'-'+'\n')
1148        for Ename in FFtable:
1149            ffdata = FFtable[Ename]
1150            fa = ffdata['fa']
1151            fb = ffdata['fb']
1152            pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1153                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc']))
1154               
1155    def PrintMFtable(MFtable):
1156        pFile.write('\n <j0> Magnetic scattering factors:\n')
1157        pFile.write('   Symbol     mfa                                    mfb                                     mfc\n')
1158        pFile.write(99*'-'+'\n')
1159        for Ename in MFtable:
1160            mfdata = MFtable[Ename]
1161            fa = mfdata['mfa']
1162            fb = mfdata['mfb']
1163            pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1164                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['mfc']))
1165        pFile.write('\n <j2> Magnetic scattering factors:\n')
1166        pFile.write('   Symbol     nfa                                    nfb                                     nfc\n')
1167        pFile.write(99*'-'+'\n')
1168        for Ename in MFtable:
1169            mfdata = MFtable[Ename]
1170            fa = mfdata['nfa']
1171            fb = mfdata['nfb']
1172            pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1173                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['nfc']))
1174               
1175    def PrintBLtable(BLtable):
1176        pFile.write('\n Neutron scattering factors:\n')
1177        pFile.write('   Symbol   isotope       mass       b       resonant terms\n')
1178        pFile.write(99*'-'+'\n')
1179        for Ename in BLtable:
1180            bldata = BLtable[Ename]
1181            isotope = bldata[0]
1182            mass = bldata[1]['Mass']
1183            if 'BW-LS' in bldata[1]:
1184                bres = bldata[1]['BW-LS']
1185                blen = 0
1186            else:
1187                blen = bldata[1]['SL'][0]
1188                bres = []
1189            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
1190            for item in bres:
1191                line += '%10.5g'%(item)
1192            pFile.write(line+'\n')
1193           
1194    def PrintRBObjects(resRBData,vecRBData):
1195               
1196        def PrintRBThermals():
1197            tlstr = ['11','22','33','12','13','23']
1198            sstr = ['12','13','21','23','31','32','AA','BB']
1199            TLS = RB['ThermalMotion'][1]
1200            TLSvar = RB['ThermalMotion'][2]
1201            if 'T' in RB['ThermalMotion'][0]:
1202                pFile.write('TLS data\n')
1203                text = ''
1204                for i in range(6):
1205                    text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
1206                pFile.write(text+'\n')
1207                if 'L' in RB['ThermalMotion'][0]: 
1208                    text = ''
1209                    for i in range(6,12):
1210                        text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
1211                    pFile.write(text+'\n')
1212                if 'S' in RB['ThermalMotion'][0]:
1213                    text = ''
1214                    for i in range(12,20):
1215                        text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
1216                    pFile.write(text+'\n')
1217            if 'U' in RB['ThermalMotion'][0]:
1218                pFile.write('Uiso data\n')
1219                text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])           
1220                pFile.write(text+'\n')
1221           
1222        if len(resRBData):
1223            for RB in resRBData:
1224                Oxyz = RB['Orig'][0]
1225                Qrijk = RB['Orient'][0]
1226                Angle = 2.0*acosd(Qrijk[0])
1227                pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'%
1228                    (RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1]))
1229                pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'%
1230                    (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1]))
1231                pFile.write('Atom site frac: %10.3f Refine? %s\n'%(RB['AtomFrac'][0],RB['AtomFrac'][1]))
1232                Torsions = RB['Torsions']
1233                if len(Torsions):
1234                    text = 'Torsions: '
1235                    for torsion in Torsions:
1236                        text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
1237                    pFile.write(text+'\n')
1238                PrintRBThermals()
1239        if len(vecRBData):
1240            for RB in vecRBData:
1241                Oxyz = RB['Orig'][0]
1242                Qrijk = RB['Orient'][0]
1243                Angle = 2.0*acosd(Qrijk[0])
1244                pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'%
1245                    (RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1]))
1246                pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'%
1247                    (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1]))
1248                pFile.write('Atom site frac: %10.3f Refine? %s\n'%(RB['AtomFrac'][0],RB['AtomFrac'][1]))
1249                PrintRBThermals()
1250               
1251    def PrintAtoms(General,Atoms):
1252        cx,ct,cs,cia = General['AtomPtrs']
1253        pFile.write('\n Atoms:\n')
1254        line = '   name    type  refine?   x         y         z    '+ \
1255            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
1256        if General['Type'] == 'macromolecular':
1257            line = ' res no residue chain'+line
1258        pFile.write(line+'\n')
1259        if General['Type'] in ['nuclear','magnetic','faulted',]:
1260            pFile.write(135*'-'+'\n')
1261            for i,at in enumerate(Atoms):
1262                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
1263                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
1264                if at[cia] == 'I':
1265                    line += '%8.5f'%(at[cia+1])+48*' '
1266                else:
1267                    line += 8*' '
1268                    for j in range(6):
1269                        line += '%8.5f'%(at[cia+2+j])
1270                pFile.write(line+'\n')
1271        elif General['Type'] == 'macromolecular':
1272            pFile.write(135*'-'+'\n')           
1273            for i,at in enumerate(Atoms):
1274                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
1275                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
1276                if at[cia] == 'I':
1277                    line += '%8.4f'%(at[cia+1])+48*' '
1278                else:
1279                    line += 8*' '
1280                    for j in range(6):
1281                        line += '%8.4f'%(at[cia+2+j])
1282                pFile.write(line+'\n')
1283               
1284    def PrintMoments(General,Atoms):
1285        cx,ct,cs,cia = General['AtomPtrs']
1286        cmx = cx+4
1287        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1288        pFile.write('\n Magnetic moments:\n')
1289        line = '   name    type  refine?  Mx        My        Mz    '
1290        pFile.write(line+'\n')
1291        pFile.write(135*'-'+'\n')
1292        for i,at in enumerate(Atoms):
1293            if AtInfo[at[ct]]:
1294                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cmx])+'%10.5f'%(at[cmx+1])+ \
1295                    '%10.5f'%(at[cmx+2])
1296                pFile.write(line+'\n')
1297       
1298               
1299    def PrintWaves(General,Atoms):
1300        cx,ct,cs,cia = General['AtomPtrs']
1301        pFile.write('\n Modulation waves\n')
1302        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1303            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1304            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1305        pFile.write(135*'-'+'\n')
1306        for i,at in enumerate(Atoms):
1307            AtomSS = at[-1]['SS1']
1308            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1309                Waves = AtomSS[Stype]
1310                if len(Waves):
1311                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
1312                        (at[ct-1],at[cs],Stype,Waves[0]))
1313                else:
1314                    continue
1315                for iw,wave in enumerate(Waves[1:]):                   
1316                    line = ''
1317                    if Waves[0] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
1318                        for item in names[Stype][6:]:
1319                            line += '%8s '%(item)                       
1320                    else:
1321                        if Stype == 'Spos':
1322                            for item in names[Stype][:6]:
1323                                line += '%8s '%(item)
1324                        else:
1325                            for item in names[Stype]:
1326                                line += '%8s '%(item)
1327                    pFile.write(line+'\n')
1328                    line = ''
1329                    for item in wave[0]:
1330                        line += '%8.4f '%(item)
1331                    line += ' Refine? '+str(wave[1])
1332                    pFile.write(line+'\n')
1333       
1334    def PrintTexture(textureData):
1335        topstr = '\n Spherical harmonics texture: Order:' + \
1336            str(textureData['Order'])
1337        if textureData['Order']:
1338            pFile.write('%s Refine? %s\n'%(topstr,textureData['SH Coeff'][0]))
1339        else:
1340            pFile.write(topstr+'\n')
1341            return
1342        names = ['omega','chi','phi']
1343        line = '\n'
1344        for name in names:
1345            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1346        pFile.write(line+'\n')
1347        pFile.write('\n Texture coefficients:\n')
1348        SHcoeff = textureData['SH Coeff'][1]
1349        SHkeys = list(SHcoeff.keys())
1350        nCoeff = len(SHcoeff)
1351        nBlock = nCoeff//10+1
1352        iBeg = 0
1353        iFin = min(iBeg+10,nCoeff)
1354        for block in range(nBlock):
1355            ptlbls = ' names :'
1356            ptstr =  ' values:'
1357            for item in SHkeys[iBeg:iFin]:
1358                ptlbls += '%12s'%(item)
1359                ptstr += '%12.4f'%(SHcoeff[item]) 
1360            pFile.write(ptlbls+'\n')
1361            pFile.write(ptstr+'\n')
1362            iBeg += 10
1363            iFin = min(iBeg+10,nCoeff)
1364       
1365    def MakeRBParms(rbKey,phaseVary,phaseDict):
1366        #### patch 2/24/21 BHT: new param, AtomFrac in RB
1367        if 'AtomFrac' not in RB: raise Exception('out of date RB: edit in RB Models')
1368        # end patch
1369        rbid = str(rbids.index(RB['RBId']))
1370        pfxRB = pfx+'RB'+rbKey+'P'
1371        pstr = ['x','y','z']
1372        ostr = ['a','i','j','k']
1373        for i in range(3):
1374            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1375            phaseDict[name] = RB['Orig'][0][i]
1376            if RB['Orig'][1]:
1377                phaseVary += [name,]
1378        pfxRB = pfx+'RB'+rbKey+'O'
1379        A,V = G2mth.Q2AV(RB['Orig'][0])
1380        fixAxis = [0, np.abs(V).argmax()+1]
1381        for i in range(4):
1382            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1383            phaseDict[name] = RB['Orient'][0][i]
1384            if RB['Orient'][1] == 'AV' and i:
1385                phaseVary += [name,]
1386            elif RB['Orient'][1] == 'A' and not i:
1387                phaseVary += [name,]
1388            elif RB['Orient'][1] == 'V' and i not in fixAxis:
1389                phaseVary += [name,]
1390        name = pfx+'RB'+rbKey+'f:'+str(iRB)+':'+rbid
1391        phaseDict[name] = RB['AtomFrac'][0]
1392        if RB['AtomFrac'][1]:
1393            phaseVary += [name,]
1394               
1395    def MakeRBThermals(rbKey,phaseVary,phaseDict):
1396        rbid = str(rbids.index(RB['RBId']))
1397        tlstr = ['11','22','33','12','13','23']
1398        sstr = ['12','13','21','23','31','32','AA','BB']
1399        if 'T' in RB['ThermalMotion'][0]:
1400            pfxRB = pfx+'RB'+rbKey+'T'
1401            for i in range(6):
1402                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1403                phaseDict[name] = RB['ThermalMotion'][1][i]
1404                if RB['ThermalMotion'][2][i]:
1405                    phaseVary += [name,]
1406        if 'L' in RB['ThermalMotion'][0]:
1407            pfxRB = pfx+'RB'+rbKey+'L'
1408            for i in range(6):
1409                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1410                phaseDict[name] = RB['ThermalMotion'][1][i+6]
1411                if RB['ThermalMotion'][2][i+6]:
1412                    phaseVary += [name,]
1413        if 'S' in RB['ThermalMotion'][0]:
1414            pfxRB = pfx+'RB'+rbKey+'S'
1415            for i in range(8):
1416                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1417                phaseDict[name] = RB['ThermalMotion'][1][i+12]
1418                if RB['ThermalMotion'][2][i+12]:
1419                    phaseVary += [name,]
1420        if 'U' in RB['ThermalMotion'][0]:
1421            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1422            phaseDict[name] = RB['ThermalMotion'][1][0]
1423            if RB['ThermalMotion'][2][0]:
1424                phaseVary += [name,]
1425               
1426    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1427        rbid = str(rbids.index(RB['RBId']))
1428        pfxRB = pfx+'RB'+rbKey+'Tr;'
1429        for i,torsion in enumerate(RB['Torsions']):
1430            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1431            phaseDict[name] = torsion[0]
1432            if torsion[1]:
1433                phaseVary += [name,]
1434                   
1435    if Print and pFile is None: raise Exception("specify pFile or Print=False")
1436    if Print:
1437        pFile.write('\n Phases:\n')
1438    phaseVary = []
1439    phaseDict = {}
1440    pawleyLookup = {}
1441    FFtables = {}                   #scattering factors - xrays
1442    MFtables = {}                   #Mag. form factors
1443    BLtables = {}                   # neutrons
1444    Natoms = {}
1445    maxSSwave = {}
1446    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1447    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1448    atomIndx = {}
1449    for name in PhaseData:
1450        General = PhaseData[name]['General']
1451        pId = PhaseData[name]['pId']
1452        pfx = str(pId)+'::'
1453        FFtable = G2el.GetFFtable(General['AtomTypes'])
1454        BLtable = G2el.GetBLtable(General)
1455        FFtables.update(FFtable)
1456        BLtables.update(BLtable)
1457        phaseDict[pfx+'isMag'] = False
1458        SGData = General['SGData']
1459        SGtext,SGtable = G2spc.SGPrint(SGData)
1460        if General['Type'] == 'magnetic':
1461            MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g'])
1462            MFtables.update(MFtable)
1463            phaseDict[pfx+'isMag'] = True
1464            SpnFlp = SGData['SpnFlp']
1465        Atoms = PhaseData[name]['Atoms']
1466        if Atoms and not General.get('doPawley'):
1467            cx,ct,cs,cia = General['AtomPtrs']
1468            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1469        PawleyRef = PhaseData[name].get('Pawley ref',[])
1470        cell = General['Cell']
1471        A = G2lat.cell2A(cell[1:7])
1472        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1473            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1474        if cell[0]:
1475            phaseVary += cellVary(pfx,SGData)       #also fills in symmetry required constraints
1476        SSGtext = []    #no superstructure
1477        im = 0
1478        if General.get('Modulated',False):
1479            im = 1      #refl offset
1480            Vec,vRef,maxH = General['SuperVec']
1481            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1482            SSGData = General['SSGData']
1483            SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1484            if vRef:
1485                phaseVary += modVary(pfx,SSGData)       
1486        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1487        if resRBData:
1488            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
1489            for iRB,RB in enumerate(resRBData):
1490                MakeRBParms('R',phaseVary,phaseDict)
1491                MakeRBThermals('R',phaseVary,phaseDict)
1492                MakeRBTorsions('R',phaseVary,phaseDict)
1493       
1494        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1495        if vecRBData:
1496            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
1497            for iRB,RB in enumerate(vecRBData):
1498                MakeRBParms('V',phaseVary,phaseDict)
1499                MakeRBThermals('V',phaseVary,phaseDict)
1500                   
1501        Natoms[pfx] = 0
1502        maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1503        if Atoms and not General.get('doPawley'):
1504            cx,ct,cs,cia = General['AtomPtrs']
1505            Natoms[pfx] = len(Atoms)
1506            for i,at in enumerate(Atoms):
1507                atomIndx[at[cia+8]] = [pfx,i]      #lookup table for restraints
1508                phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1509                    pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1510                    pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1511                    pfx+'AI/A:'+str(i):at[cia],})
1512                if at[cia] == 'I':
1513                    phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1514                else:
1515                    phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1516                        pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1517                if General['Type'] == 'magnetic':
1518                    phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]})
1519                if 'F' in at[ct+1]:
1520                    phaseVary.append(pfx+'Afrac:'+str(i))
1521                if 'X' in at[ct+1]:
1522                    try:    #patch for sytsym name changes
1523                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1524                    except KeyError:
1525                        Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1526                        at[cs] = Sytsym
1527                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1528                    names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1529                    equivs = [[],[],[]]
1530                    for j in range(3):
1531                        if xId[j] > 0:                               
1532                            phaseVary.append(names[j])
1533                            equivs[xId[j]-1].append([names[j],xCoef[j]])
1534                    for equiv in equivs:
1535                        if len(equiv) > 1:
1536                            name = equiv[0][0]
1537                            coef = equiv[0][1]
1538                            for eqv in equiv[1:]:
1539                                eqv[1] /= coef
1540                                G2mv.StoreEquivalence(name,(eqv,))
1541                if 'U' in at[ct+1]:
1542                    if at[cia] == 'I':
1543                        phaseVary.append(pfx+'AUiso:'+str(i))
1544                    else:
1545                        try:    #patch for sytsym name changes
1546                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1547                        except KeyError:
1548                            Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1549                            at[cs] = Sytsym
1550                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1551                        names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1552                            pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1553                        equivs = [[],[],[],[],[],[]]
1554                        for j in range(6):
1555                            if uId[j] > 0:                               
1556                                phaseVary.append(names[j])
1557                                equivs[uId[j]-1].append([names[j],uCoef[j]])
1558                        for equiv in equivs:
1559                            if len(equiv) > 1:
1560                                name = equiv[0][0]
1561                                coef = equiv[0][1]
1562                                for eqv in equiv[1:]:
1563                                    eqv[1] /= coef
1564                                    G2mv.StoreEquivalence(name,(eqv,))
1565                if 'M' in at[ct+1]:
1566                    SytSym,Mul,Nop,dupDir = G2spc.SytSym(at[cx:cx+3],SGData)
1567                    mId,mCoef = G2spc.GetCSpqinel(SpnFlp,dupDir)
1568                    names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1569                    equivs = [[],[],[]]
1570                    for j in range(3):
1571                        if mId[j] > 0:
1572                            phaseVary.append(names[j])
1573                            equivs[mId[j]-1].append([names[j],mCoef[j]])
1574                    for equiv in equivs:
1575                        if len(equiv) > 1:
1576                            name = equiv[0][0]
1577                            coef = equiv[0][1]
1578                            for eqv in equiv[1:]:
1579                                eqv[1] /= coef
1580                                G2mv.StoreEquivalence(name,(eqv,))
1581                if General.get('Modulated',False):
1582                    AtomSS = at[-1]['SS1']
1583                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1584                        Waves = AtomSS[Stype]
1585                        if len(Waves):
1586                            waveType = Waves[0]
1587                        else:
1588                            continue
1589                        phaseDict[pfx+Stype[1].upper()+'waveType:'+str(i)] = waveType
1590                        nx = 0
1591                        for iw,wave in enumerate(Waves[1:]):
1592                            if not iw:
1593                                if waveType in ['ZigZag','Block']:
1594                                    nx = 1
1595                                CSI = G2spc.GetSSfxuinel(waveType,Stype,1,at[cx:cx+3],SGData,SSGData)
1596                            else:
1597                                CSI = G2spc.GetSSfxuinel('Fourier',Stype,iw+1-nx,at[cx:cx+3],SGData,SSGData)
1598                            uId,uCoef = CSI[0]
1599                            stiw = str(i)+':'+str(iw)
1600                            if Stype == 'Spos':
1601                                if waveType in ['ZigZag','Block',] and not iw:
1602                                    names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1603                                    equivs = [[],[], [],[],[]]
1604                                else:
1605                                    names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1606                                        pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1607                                    equivs = [[],[],[], [],[],[]]
1608                            elif Stype == 'Sadp':
1609                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1610                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1611                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1612                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1613                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1614                            elif Stype == 'Sfrac':
1615                                equivs = [[],[]]
1616                                if 'Crenel' in waveType and not iw:
1617                                    names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1618                                else:
1619                                    names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1620                            elif Stype == 'Smag':
1621                                equivs = [[],[],[], [],[],[]]
1622                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1623                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1624                            phaseDict.update(dict(zip(names,wave[0])))
1625                            if wave[1]: #what do we do here for multiple terms in modulation constraints?
1626                                for j in range(len(equivs)):
1627                                    if uId[j][0] > 0:                               
1628                                        phaseVary.append(names[j])
1629                                        equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1630                                for equiv in equivs:
1631                                    if len(equiv) > 1:
1632                                        name = equiv[0][0]
1633                                        coef = equiv[0][1]
1634                                        for eqv in equiv[1:]:
1635                                            eqv[1] /= coef
1636                                            G2mv.StoreEquivalence(name,(eqv,))
1637                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1638            textureData = General['SH Texture']
1639            if textureData['Order'] and not seqRef:
1640                phaseDict[pfx+'SHorder'] = textureData['Order']
1641                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1642                for item in ['omega','chi','phi']:
1643                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1644                    if textureData['Sample '+item][0]:
1645                        phaseVary.append(pfx+'SH '+item)
1646                for item in textureData['SH Coeff'][1]:
1647                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1648                    if textureData['SH Coeff'][0]:
1649                        phaseVary.append(pfx+item)
1650               
1651            if Print:
1652                pFile.write('\n Phase name: %s\n'%General['Name'])
1653                pFile.write(135*'='+'\n')
1654                PrintFFtable(FFtable)
1655                PrintBLtable(BLtable)
1656                if General['Type'] == 'magnetic':
1657                    PrintMFtable(MFtable)
1658                pFile.write('\n')
1659                #how do we print magnetic symmetry table? TBD
1660                if len(SSGtext):    #if superstructure
1661                    for line in SSGtext: pFile.write(line+'\n')
1662                    if len(SSGtable):
1663                        for item in SSGtable:
1664                            line = ' %s '%(item)
1665                            pFile.write(line+'\n') 
1666                    else:
1667                        pFile.write(' ( 1)    %s\n'%(SSGtable[0]))
1668                else:
1669                    for line in SGtext: pFile.write(line+'\n')
1670                    if len(SGtable):
1671                        for item in SGtable:
1672                            line = ' %s '%(item)
1673                            pFile.write(line+'\n') 
1674                    else:
1675                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1676                PrintRBObjects(resRBData,vecRBData)
1677                PrintAtoms(General,Atoms)
1678                if General['Type'] == 'magnetic':
1679                    PrintMoments(General,Atoms)
1680                if General.get('Modulated',False):
1681                    PrintWaves(General,Atoms)
1682                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1683                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1684                if len(SSGtext):    #if superstructure
1685                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1686                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1687                if not seqRef:
1688                    PrintTexture(textureData)
1689                if name in RestraintDict:
1690                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1691                        textureData,RestraintDict[name],pFile)
1692                   
1693        elif PawleyRef:
1694            if Print:
1695                pFile.write('\n Phase name: %s\n'%General['Name'])
1696                pFile.write(135*'='+'\n')
1697                pFile.write('\n')
1698                if len(SSGtext):    #if superstructure
1699                    for line in SSGtext: pFile.write(line+'\n')
1700                    if len(SSGtable):
1701                        for item in SSGtable:
1702                            line = ' %s '%(item)
1703                            pFile.write(line+'\n') 
1704                    else:
1705                        pFile.write(' ( 1)    %s\n'%SSGtable[0])
1706                else:
1707                    for line in SGtext: pFile.write(line+'\n')
1708                    if len(SGtable):
1709                        for item in SGtable:
1710                            line = ' %s '%(item)
1711                            pFile.write(line+'\n')
1712                    else:
1713                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1714                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1715                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1716                if len(SSGtext):    #if superstructure
1717                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1718                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1719            pawleyVary = []
1720            for i,refl in enumerate(PawleyRef):
1721                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1722                if im:
1723                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1724                else:
1725                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1726                if refl[5+im]:
1727                    pawleyVary.append(pfx+'PWLref:'+str(i))
1728            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1729            phaseVary += pawleyVary
1730               
1731    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1732   
1733def cellFill(pfx,SGData,parmDict,sigDict): 
1734    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1735    from the parameter and sig dictionaries.
1736
1737    :param str pfx: parameter prefix ("n::", where n is a phase number)
1738    :param dict SGdata: a symmetry object
1739    :param dict parmDict: a dictionary of parameters
1740    :param dict sigDict:  a dictionary of uncertainties on parameters
1741
1742    :returns: A,sigA where each is a list of six terms with the A terms
1743    '''
1744    if SGData['SGLaue'] in ['-1',]:
1745        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1746            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1747    elif SGData['SGLaue'] in ['2/m',]:
1748        if SGData['SGUniq'] == 'a':
1749            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1750                0,0,parmDict[pfx+'A5']]
1751        elif SGData['SGUniq'] == 'b':
1752            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1753                0,parmDict[pfx+'A4'],0]
1754        else:
1755            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1756                parmDict[pfx+'A3'],0,0]
1757    elif SGData['SGLaue'] in ['mmm',]:
1758        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1759    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1760        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1761    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1762        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1763            parmDict[pfx+'A0'],0,0]
1764    elif SGData['SGLaue'] in ['3R', '3mR']:
1765        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1766            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1767    elif SGData['SGLaue'] in ['m3m','m3']:
1768        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1769
1770    try:
1771        if SGData['SGLaue'] in ['-1',]:
1772            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1773                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1774        elif SGData['SGLaue'] in ['2/m',]:
1775            if SGData['SGUniq'] == 'a':
1776                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1777                    0,0,sigDict[pfx+'A5']]
1778            elif SGData['SGUniq'] == 'b':
1779                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1780                    0,sigDict[pfx+'A4'],0]
1781            else:
1782                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1783                    sigDict[pfx+'A3'],0,0]
1784        elif SGData['SGLaue'] in ['mmm',]:
1785            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1786        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1787            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1788        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1789            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1790        elif SGData['SGLaue'] in ['3R', '3mR']:
1791            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1792        elif SGData['SGLaue'] in ['m3m','m3']:
1793            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1794    except KeyError:
1795        sigA = [0,0,0,0,0,0]
1796    return A,sigA
1797       
1798def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1799    'needs a doc string'
1800    if phaseRest:
1801        Amat = G2lat.cell2AB(cell)[0]
1802        cx,ct,cs = AtPtrs[:3]
1803        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1804            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1805            ['ChemComp','Sites'],['Texture','HKLs']]
1806        for name,rest in names:
1807            itemRest = phaseRest[name]
1808            if rest in itemRest and itemRest[rest] and itemRest['Use']:
1809                pFile.write('\n  %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use'])))
1810                if name in ['Bond','Angle','Plane','Chiral']:
1811                    pFile.write('     calc       obs      sig   delt/sig  atoms(symOp): \n')
1812                    for indx,ops,obs,esd in itemRest[rest]:
1813                        try:
1814                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1815                            AtName = ''
1816                            for i,Aname in enumerate(AtNames):
1817                                AtName += Aname
1818                                if ops[i] == '1':
1819                                    AtName += '-'
1820                                else:
1821                                    AtName += '+('+ops[i]+')-'
1822                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1823                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1824                            if name == 'Bond':
1825                                calc = G2mth.getRestDist(XYZ,Amat)
1826                            elif name == 'Angle':
1827                                calc = G2mth.getRestAngle(XYZ,Amat)
1828                            elif name == 'Plane':
1829                                calc = G2mth.getRestPlane(XYZ,Amat)
1830                            elif name == 'Chiral':
1831                                calc = G2mth.getRestChiral(XYZ,Amat)
1832                            pFile.write(' %9.3f %9.3f %8.3f %8.3f   %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1]))
1833                        except KeyError:
1834                            del itemRest[rest]
1835                elif name in ['Torsion','Rama']:
1836                    pFile.write('  atoms(symOp)  calc  obs  sig  delt/sig  torsions: \n')
1837                    coeffDict = itemRest['Coeff']
1838                    for indx,ops,cofName,esd in itemRest[rest]:
1839                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1840                        AtName = ''
1841                        for i,Aname in enumerate(AtNames):
1842                            AtName += Aname+'+('+ops[i]+')-'
1843                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1844                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1845                        if name == 'Torsion':
1846                            tor = G2mth.getRestTorsion(XYZ,Amat)
1847                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1848                            pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1]))
1849                        else:
1850                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1851                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1852                            pFile.write(' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1]))
1853                elif name == 'ChemComp':
1854                    pFile.write('     atoms   mul*frac  factor     prod\n')
1855                    for indx,factors,obs,esd in itemRest[rest]:
1856                        try:
1857                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1858                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1859                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1860                            mulfrac = mul*frac
1861                            calcs = mul*frac*factors
1862                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1863                                pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc))
1864                            pFile.write(' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd))
1865                        except KeyError:
1866                            del itemRest[rest]
1867                elif name == 'Texture' and textureData['Order']:
1868                    Start = False
1869                    SHCoef = textureData['SH Coeff'][1]
1870                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1871                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1872                    pFile.write ('    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd \n')
1873                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1874                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1875                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1876                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1877                        Z = ma.masked_greater(Z,0.0)
1878                        num = ma.count(Z)
1879                        sum = 0
1880                        if num:
1881                            sum = np.sum(Z)
1882                        pFile.write ('   %d %d %d  %d %8.3f %8.3f %8d   %s    %8.3f\n'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2))
1883       
1884def getCellEsd(pfx,SGData,A,covData):
1885    'needs a doc string'
1886    rVsq = G2lat.calc_rVsq(A)
1887    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1888    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1889    varyList = covData['varyList']
1890    covMatrix = covData['covMatrix']
1891    if len(covMatrix):
1892        vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1893        if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']:
1894            vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0]
1895            vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0]
1896            vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2]
1897        elif SGData['SGLaue'] in ['m3','m3m']:
1898            vcov[0:3,0:3] = vcov[0,0]
1899        elif SGData['SGLaue'] in ['4/m', '4/mmm']:
1900            vcov[0:2,0:2] = vcov[0,0]
1901            vcov[1,2] = vcov[2,1] = vcov[0,2]
1902        elif SGData['SGLaue'] in ['3R','3mR']:
1903            vcov[0:3,0:3] = vcov[0,0]
1904    #        vcov[4,4] = vcov[5,5] = vcov[3,3]
1905            vcov[3:6,3:6] = vcov[3,3]
1906            vcov[0:3,3:6] = vcov[0,3]
1907            vcov[3:6,0:3] = vcov[3,0]
1908    else:
1909        vcov = np.eye(6)
1910    delt = 1.e-9
1911    drVdA = np.zeros(6)
1912    for i in range(6):
1913        A[i] += delt
1914        drVdA[i] = G2lat.calc_rVsq(A)
1915        A[i] -= 2*delt
1916        drVdA[i] -= G2lat.calc_rVsq(A)
1917        A[i] += delt
1918    drVdA /= 2.*delt   
1919    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1920    Vol = 1/np.sqrt(rVsq)
1921    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1922   
1923    dcdA = np.zeros((6,6))
1924    for i in range(6):
1925        pdcdA =np.zeros(6)
1926        A[i] += delt
1927        pdcdA += G2lat.A2cell(A)
1928        A[i] -= 2*delt
1929        pdcdA -= G2lat.A2cell(A)
1930        A[i] += delt
1931        dcdA[i] = pdcdA/(2.*delt)
1932   
1933    sigMat = np.inner(dcdA,np.inner(dcdA,vcov))
1934    var = np.diag(sigMat)
1935    CS = np.where(var>0.,np.sqrt(var),0.)
1936    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']:
1937        CS[3:6] = 0.0
1938    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1939   
1940def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1941    '''Called after a refinement to transfer parameters from the parameter dict to
1942    the phase(s) information read from a GPX file. Also prints values to the .lst file
1943    '''
1944   
1945    def PrintAtomsAndSig(General,Atoms,atomsSig):
1946        pFile.write('\n Atoms:\n')
1947        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1948        if General['Type'] == 'macromolecular':
1949            line = ' res no residue chain '+line
1950        cx,ct,cs,cia = General['AtomPtrs']
1951        pFile.write(line+'\n')
1952        pFile.write(135*'-'+'\n')
1953        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1954            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1955        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1956        for atyp in General['AtomTypes']:       #zero composition
1957            General['NoAtoms'][atyp] = 0.
1958        for i,at in enumerate(Atoms):
1959            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1960            if General['Type'] == 'macromolecular':
1961                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1962                valstr = ' values:          '
1963                sigstr = ' sig   :          '
1964            else:
1965                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1966                valstr = ' values:'
1967                sigstr = ' sig   :'
1968            for ind in range(cx,cx+4):
1969                sigind = str(i)+':'+str(ind)
1970                valstr += fmt[ind]%(at[ind])                   
1971                if sigind in atomsSig:
1972                    sigstr += fmt[ind]%(atomsSig[sigind])
1973                else:
1974                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1975            if at[cia] == 'I':
1976                valstr += fmt[cia+1]%(at[cia+1])
1977                if '%d:%d'%(i,cia+1) in atomsSig:
1978                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1979                else:
1980                    sigstr += 8*' '
1981            else:
1982                valstr += 8*' '
1983                sigstr += 8*' '
1984                for ind in range(cia+2,cia+8):
1985                    sigind = str(i)+':'+str(ind)
1986                    valstr += fmt[ind]%(at[ind])
1987                    if sigind in atomsSig:                       
1988                        sigstr += fmt[ind]%(atomsSig[sigind])
1989                    else:
1990                        sigstr += 8*' '
1991            pFile.write(name+'\n')
1992            pFile.write(valstr+'\n')
1993            pFile.write(sigstr+'\n')
1994           
1995    def PrintMomentsAndSig(General,Atoms,atomsSig):
1996        cell = General['Cell'][1:7]
1997        G = G2lat.fillgmat(cell)
1998        ast = np.sqrt(np.diag(G))
1999        GS = G/np.outer(ast,ast)
2000        pFile.write('\n Magnetic Moments:\n')    #add magnitude & angle, etc.? TBD
2001        line = '   name       Mx        My        Mz       |Mag|'
2002        cx,ct,cs,cia = General['AtomPtrs']
2003        cmx = cx+4
2004        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
2005        pFile.write(line+'\n')
2006        pFile.write(135*'-'+'\n')
2007        fmt = {0:'%7s',ct:'%7s',cmx:'%10.3f',cmx+1:'%10.3f',cmx+2:'%10.3f'}
2008        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
2009        for i,at in enumerate(Atoms):
2010            if AtInfo[at[ct]]:
2011                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
2012                valstr = ' values:'
2013                sigstr = ' sig   :'
2014                for ind in range(cmx,cmx+3):
2015                    sigind = str(i)+':'+str(ind)
2016                    valstr += fmt[ind]%(at[ind])                   
2017                    if sigind in atomsSig:
2018                        sigstr += fmt[ind]%(atomsSig[sigind])
2019                    else:
2020                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
2021                mag = np.array(at[cmx:cmx+3])
2022                Mag = np.sqrt(np.inner(mag,np.inner(mag,GS)))
2023                valstr += '%10.3f'%Mag
2024                sigstr += 10*' '
2025                pFile.write(name+'\n')
2026                pFile.write(valstr+'\n')
2027                pFile.write(sigstr+'\n')
2028           
2029    def PrintWavesAndSig(General,Atoms,wavesSig):
2030        cx,ct,cs,cia = General['AtomPtrs']
2031        pFile.write('\n Modulation waves\n')
2032        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
2033            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
2034            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
2035        pFile.write(135*'-'+'\n')
2036        for i,at in enumerate(Atoms):
2037            AtomSS = at[-1]['SS1']
2038            for Stype in ['Sfrac','Spos','Sadp','Smag']:
2039                Waves = AtomSS[Stype]
2040                if len(Waves) > 1:
2041                    waveType = Waves[0]
2042                else:
2043                    continue
2044                if len(Waves):
2045                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
2046                        (at[ct-1],at[cs],Stype,waveType))
2047                    for iw,wave in enumerate(Waves[1:]):
2048                        stiw = ':'+str(i)+':'+str(iw)
2049                        namstr = '  names :'
2050                        valstr = '  values:'
2051                        sigstr = '  esds  :'
2052                        if Stype == 'Spos':
2053                            nt = 6
2054                            ot = 0
2055                            if waveType in ['ZigZag','Block',] and not iw:
2056                                nt = 5
2057                                ot = 6
2058                            for j in range(nt):
2059                                name = names['Spos'][j+ot]
2060                                namstr += '%12s'%(name)
2061                                valstr += '%12.4f'%(wave[0][j])
2062                                if name+stiw in wavesSig:
2063                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2064                                else:
2065                                    sigstr += 12*' '
2066                        elif Stype == 'Sfrac':
2067                            ot = 0
2068                            if 'Crenel' in waveType and not iw:
2069                                ot = 2
2070                            for j in range(2):
2071                                name = names['Sfrac'][j+ot]
2072                                namstr += '%12s'%(names['Sfrac'][j+ot])
2073                                valstr += '%12.4f'%(wave[0][j])
2074                                if name+stiw in wavesSig:
2075                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2076                                else:
2077                                    sigstr += 12*' '
2078                        elif Stype == 'Sadp':
2079                            for j in range(12):
2080                                name = names['Sadp'][j]
2081                                namstr += '%10s'%(names['Sadp'][j])
2082                                valstr += '%10.6f'%(wave[0][j])
2083                                if name+stiw in wavesSig:
2084                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
2085                                else:
2086                                    sigstr += 10*' '
2087                        elif Stype == 'Smag':
2088                            for j in range(6):
2089                                name = names['Smag'][j]
2090                                namstr += '%12s'%(names['Smag'][j])
2091                                valstr += '%12.4f'%(wave[0][j])
2092                                if name+stiw in wavesSig:
2093                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2094                                else:
2095                                    sigstr += 12*' '
2096                               
2097                    pFile.write(namstr+'\n')
2098                    pFile.write(valstr+'\n')
2099                    pFile.write(sigstr+'\n')
2100       
2101               
2102    def PrintRBObjPOAndSig(rbfx,rbsx):
2103        for i in WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
2104            pFile.write(i+'\n')
2105       
2106    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
2107        for i in WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
2108            pFile.write(i)
2109       
2110    def PrintRBObjTorAndSig(rbsx):
2111        nTors = len(RBObj['Torsions'])
2112        if nTors:
2113            for i in WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
2114                pFile.write(i)
2115               
2116    def PrintSHtextureAndSig(textureData,SHtextureSig):
2117        Tindx = 1.0
2118        Tvar = 0.0
2119        pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
2120        names = ['omega','chi','phi']
2121        namstr = '  names :'
2122        ptstr =  '  values:'
2123        sigstr = '  esds  :'
2124        for name in names:
2125            namstr += '%12s'%(name)
2126            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2127            if 'Sample '+name in SHtextureSig:
2128                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2129            else:
2130                sigstr += 12*' '
2131        pFile.write(namstr+'\n')
2132        pFile.write(ptstr+'\n')
2133        pFile.write(sigstr+'\n')
2134        pFile.write('\n Texture coefficients:\n')
2135        SHcoeff = textureData['SH Coeff'][1]
2136        SHkeys = list(SHcoeff.keys())
2137        nCoeff = len(SHcoeff)
2138        nBlock = nCoeff//10+1
2139        iBeg = 0
2140        iFin = min(iBeg+10,nCoeff)
2141        for block in range(nBlock):
2142            namstr = '  names :'
2143            ptstr =  '  values:'
2144            sigstr = '  esds  :'
2145            for name in SHkeys[iBeg:iFin]:
2146                namstr += '%12s'%(name)
2147                ptstr += '%12.3f'%(SHcoeff[name])
2148                l = 2.0*eval(name.strip('C'))[0]+1
2149                Tindx += SHcoeff[name]**2/l
2150                if name in SHtextureSig:
2151                    Tvar += (2.*SHcoeff[name]*SHtextureSig[name]/l)**2
2152                    sigstr += '%12.3f'%(SHtextureSig[name])
2153                else:
2154                    sigstr += 12*' '
2155            pFile.write(namstr+'\n')
2156            pFile.write(ptstr+'\n')
2157            pFile.write(sigstr+'\n')
2158            iBeg += 10
2159            iFin = min(iBeg+10,nCoeff)
2160        pFile.write(' Texture index J = %.3f(%d)'%(Tindx,int(1000*np.sqrt(Tvar))))
2161           
2162    ##########################################################################
2163    # SetPhaseData starts here
2164    if pFile: pFile.write('\n Phases:\n')
2165    for phase in Phases:
2166        if pFile: pFile.write(' Result for phase: %s\n'%phase)
2167        if pFile: pFile.write(135*'='+'\n')
2168        Phase = Phases[phase]
2169        General = Phase['General']
2170        SGData = General['SGData']
2171        Atoms = Phase['Atoms']
2172        AtLookup = []
2173        if Atoms and not General.get('doPawley'):
2174            cx,ct,cs,cia = General['AtomPtrs']
2175            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
2176        cell = General['Cell']
2177        pId = Phase['pId']
2178        pfx = str(pId)+'::'
2179        if cell[0]:
2180            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
2181            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
2182            if pFile: pFile.write(' Reciprocal metric tensor: \n')
2183            ptfmt = "%15.9f"
2184            names = ['A11','A22','A33','A12','A13','A23']
2185            namstr = '  names :'
2186            ptstr =  '  values:'
2187            sigstr = '  esds  :'
2188            for name,a,siga in zip(names,A,sigA):
2189                namstr += '%15s'%(name)
2190                ptstr += ptfmt%(a)
2191                if siga:
2192                    sigstr += ptfmt%(siga)
2193                else:
2194                    sigstr += 15*' '
2195            if pFile: pFile.write(namstr+'\n')
2196            if pFile: pFile.write(ptstr+'\n')
2197            if pFile: pFile.write(sigstr+'\n')
2198            cell[1:7] = G2lat.A2cell(A)
2199            cell[7] = G2lat.calc_V(A)
2200            if pFile: pFile.write(' New unit cell:\n')
2201            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
2202            names = ['a','b','c','alpha','beta','gamma','Volume']
2203            namstr = '  names :'
2204            ptstr =  '  values:'
2205            sigstr = '  esds  :'
2206            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
2207                namstr += '%12s'%(name)
2208                ptstr += fmt%(a)
2209                if siga:
2210                    sigstr += fmt%(siga)
2211                else:
2212                    sigstr += 12*' '
2213            if pFile: pFile.write(namstr+'\n')
2214            if pFile: pFile.write(ptstr+'\n')
2215            if pFile: pFile.write(sigstr+'\n')
2216        ik = 6  #for Pawley stuff below
2217        if General.get('Modulated',False):
2218            ik = 7
2219            Vec,vRef,maxH = General['SuperVec']
2220            if vRef:
2221                if pFile: pFile.write(' New modulation vector:\n')
2222                namstr = '  names :'
2223                ptstr =  '  values:'
2224                sigstr = '  esds  :'
2225                for iv,var in enumerate(['mV0','mV1','mV2']):
2226                    namstr += '%12s'%(pfx+var)
2227                    ptstr += '%12.6f'%(parmDict[pfx+var])
2228                    if pfx+var in sigDict:
2229                        Vec[iv] = parmDict[pfx+var]
2230                        sigstr += '%12.6f'%(sigDict[pfx+var])
2231                    else:
2232                        sigstr += 12*' '
2233                if pFile: pFile.write(namstr+'\n')
2234                if pFile: pFile.write(ptstr+'\n')
2235                if pFile: pFile.write(sigstr+'\n')
2236           
2237        General['Mass'] = 0.
2238        if Phase['General'].get('doPawley'):
2239            pawleyRef = Phase['Pawley ref']
2240            for i,refl in enumerate(pawleyRef):
2241                key = pfx+'PWLref:'+str(i)
2242                refl[ik] = parmDict[key]
2243                if key in sigDict:
2244                    refl[ik+1] = sigDict[key]
2245                else:
2246                    refl[ik+1] = 0
2247        else:
2248            VRBIds = RBIds['Vector']
2249            RRBIds = RBIds['Residue']
2250            RBModels = Phase['RBModels']
2251            if pFile:
2252                for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2253                    jrb = VRBIds.index(RBObj['RBId'])
2254                    rbsx = str(irb)+':'+str(jrb)
2255                    pFile.write(' Vector rigid body parameters:\n')
2256                    PrintRBObjPOAndSig('RBV',rbsx)
2257                    PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2258                for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2259                    jrb = RRBIds.index(RBObj['RBId'])
2260                    rbsx = str(irb)+':'+str(jrb)
2261                    pFile.write(' Residue rigid body parameters:\n')
2262                    PrintRBObjPOAndSig('RBR',rbsx)
2263                    PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2264                    PrintRBObjTorAndSig(rbsx)
2265            atomsSig = {}
2266            wavesSig = {}
2267            cx,ct,cs,cia = General['AtomPtrs']
2268            for i,at in enumerate(Atoms):
2269                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2270                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2271                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2272                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2273                for ind in range(cx,cx+4):
2274                    at[ind] = parmDict[names[ind]]
2275                    if ind in range(cx,cx+3):
2276                        name = names[ind].replace('A','dA')
2277                    else:
2278                        name = names[ind]
2279                    if name in sigDict:
2280                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2281                if at[cia] == 'I':
2282                    at[cia+1] = parmDict[names[cia+1]]
2283                    if names[cia+1] in sigDict:
2284                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2285                else:
2286                    for ind in range(cia+2,cia+8):
2287                        at[ind] = parmDict[names[ind]]
2288                        if names[ind] in sigDict:
2289                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2290                if General['Type'] == 'magnetic':
2291                    for ind in range(cx+4,cx+7):
2292                        at[ind] = parmDict[names[ind]]
2293                        if names[ind] in sigDict:
2294                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2295                ind = General['AtomTypes'].index(at[ct])
2296                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2297                if General.get('Modulated',False):
2298                    AtomSS = at[-1]['SS1']
2299                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2300                        Waves = AtomSS[Stype]
2301                        if len(Waves):
2302                            waveType = Waves[0]
2303                        else:
2304                            continue
2305                        for iw,wave in enumerate(Waves[1:]):
2306                            stiw = str(i)+':'+str(iw)
2307                            if Stype == 'Spos':
2308                                if waveType in ['ZigZag','Block',] and not iw:
2309                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2310                                else:
2311                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2312                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2313                            elif Stype == 'Sadp':
2314                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2315                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2316                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2317                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2318                            elif Stype == 'Sfrac':
2319                                if 'Crenel' in waveType and not iw:
2320                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2321                                else:
2322                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2323                            elif Stype == 'Smag':
2324                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2325                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2326                            for iname,name in enumerate(names):
2327                                AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2328                                if pfx+name in sigDict:
2329                                    wavesSig[name] = sigDict[pfx+name]
2330
2331            if pFile: PrintAtomsAndSig(General,Atoms,atomsSig)
2332            if pFile and General['Type'] == 'magnetic':
2333                PrintMomentsAndSig(General,Atoms,atomsSig)
2334            if pFile and General.get('Modulated',False):
2335                PrintWavesAndSig(General,Atoms,wavesSig)
2336               
2337            density = G2mth.getDensity(General)[0]
2338            if pFile: pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density))
2339           
2340       
2341        textureData = General['SH Texture']   
2342        if textureData['Order']:
2343            SHtextureSig = {}
2344            for name in ['omega','chi','phi']:
2345                aname = pfx+'SH '+name
2346                textureData['Sample '+name][1] = parmDict[aname]
2347                if aname in sigDict:
2348                    SHtextureSig['Sample '+name] = sigDict[aname]
2349            for name in textureData['SH Coeff'][1]:
2350                aname = pfx+name
2351                textureData['SH Coeff'][1][name] = parmDict[aname]
2352                if aname in sigDict:
2353                    SHtextureSig[name] = sigDict[aname]
2354            PrintSHtextureAndSig(textureData,SHtextureSig)
2355        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2356            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2357                textureData,RestraintDict[phase],pFile)
2358                   
2359################################################################################
2360##### Histogram & Phase data
2361################################################################################       
2362                   
2363def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2364    '''Loads the HAP histogram/phase information into dicts
2365
2366    :param dict Phases: phase information
2367    :param dict Histograms: Histogram information
2368    :param bool Print: prints information as it is read
2369    :param file pFile: file object to print to (the default, None causes printing to the console)
2370    :param bool resetRefList: Should the contents of the Reflection List be initialized
2371      on loading. The default, True, initializes the Reflection List as it is loaded.
2372
2373    :returns: (hapVary,hapDict,controlDict)
2374      * hapVary: list of refined variables
2375      * hapDict: dict with refined variables and their values
2376      * controlDict: dict with fixed parameters
2377    '''
2378   
2379    def PrintSize(hapData):
2380        if hapData[0] in ['isotropic','uniaxial']:
2381            line = '\n Size model    : %9s'%(hapData[0])
2382            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2383            if hapData[0] == 'uniaxial':
2384                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2385            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2386            pFile.write(line+'\n')
2387        else:
2388            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2389                (hapData[0],hapData[1][2],hapData[2][2]))
2390            Snames = ['S11','S22','S33','S12','S13','S23']
2391            ptlbls = ' names :'
2392            ptstr =  ' values:'
2393            varstr = ' refine:'
2394            for i,name in enumerate(Snames):
2395                ptlbls += '%12s' % (name)
2396                ptstr += '%12.3f' % (hapData[4][i])
2397                varstr += '%12s' % (str(hapData[5][i]))
2398            pFile.write(ptlbls+'\n')
2399            pFile.write(ptstr+'\n')
2400            pFile.write(varstr+'\n')
2401       
2402    def PrintMuStrain(hapData,SGData):
2403        if hapData[0] in ['isotropic','uniaxial']:
2404            line = '\n Mustrain model: %9s'%(hapData[0])
2405            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2406            if hapData[0] == 'uniaxial':
2407                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2408            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2409            pFile.write(line+'\n')
2410        else:
2411            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2412                (hapData[0],hapData[1][2],hapData[2][2]))
2413            Snames = G2spc.MustrainNames(SGData)
2414            ptlbls = ' names :'
2415            ptstr =  ' values:'
2416            varstr = ' refine:'
2417            for i,name in enumerate(Snames):
2418                ptlbls += '%12s' % (name)
2419                ptstr += '%12.1f' % (hapData[4][i])
2420                varstr += '%12s' % (str(hapData[5][i]))
2421            pFile.write(ptlbls+'\n')
2422            pFile.write(ptstr+'\n')
2423            pFile.write(varstr+'\n')
2424
2425    def PrintHStrain(hapData,SGData):
2426        pFile.write('\n Hydrostatic/elastic strain:\n')
2427        Hsnames = G2spc.HStrainNames(SGData)
2428        ptlbls = ' names :'
2429        ptstr =  ' values:'
2430        varstr = ' refine:'
2431        for i,name in enumerate(Hsnames):
2432            ptlbls += '%12s' % (name)
2433            ptstr += '%12.4g' % (hapData[0][i])
2434            varstr += '%12s' % (str(hapData[1][i]))
2435        pFile.write(ptlbls+'\n')
2436        pFile.write(ptstr+'\n')
2437        pFile.write(varstr+'\n')
2438
2439    def PrintSHPO(hapData):
2440        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2441        ptlbls = ' names :'
2442        ptstr =  ' values:'
2443        for item in hapData[5]:
2444            ptlbls += '%12s'%(item)
2445            ptstr += '%12.3f'%(hapData[5][item]) 
2446        pFile.write(ptlbls+'\n')
2447        pFile.write(ptstr+'\n')
2448   
2449    def PrintBabinet(hapData):
2450        pFile.write('\n Babinet form factor modification:\n')
2451        ptlbls = ' names :'
2452        ptstr =  ' values:'
2453        varstr = ' refine:'
2454        for name in ['BabA','BabU']:
2455            ptlbls += '%12s' % (name)
2456            ptstr += '%12.6f' % (hapData[name][0])
2457            varstr += '%12s' % (str(hapData[name][1]))
2458        pFile.write(ptlbls+'\n')
2459        pFile.write(ptstr+'\n')
2460        pFile.write(varstr+'\n')
2461       
2462    hapDict = {}
2463    hapVary = []
2464    controlDict = {}
2465   
2466    for phase in Phases:
2467        HistoPhase = Phases[phase]['Histograms']
2468        SGData = Phases[phase]['General']['SGData']
2469        cell = Phases[phase]['General']['Cell'][1:7]
2470        A = G2lat.cell2A(cell)
2471        if Phases[phase]['General'].get('Modulated',False):
2472            SSGData = Phases[phase]['General']['SSGData']
2473            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2474        pId = Phases[phase]['pId']
2475        for histogram in Histograms:
2476            if histogram not in HistoPhase and phase in Histograms[histogram]['Reflection Lists']:
2477                #remove previously created reflection list if histogram is removed from phase
2478                #print("removing ",phase,"from",histogram)
2479                del Histograms[histogram]['Reflection Lists'][phase]
2480        histoList = list(HistoPhase.keys())
2481        histoList.sort()
2482        for histogram in histoList:
2483            try:
2484                Histogram = Histograms[histogram]
2485            except KeyError:                       
2486                #skip if histogram not included e.g. in a sequential refinement
2487                continue
2488            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2489                if phase in Histograms[histogram]['Reflection Lists']:
2490                    del Histograms[histogram]['Reflection Lists'][phase]
2491                continue
2492            hapData = HistoPhase[histogram]
2493            hId = Histogram['hId']
2494            if 'PWDR' in histogram:
2495                limits = Histogram['Limits'][1]
2496                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2497                if 'C' in inst['Type'][1]:
2498                    try:
2499                        wave = inst['Lam'][1]
2500                    except KeyError:
2501                        wave = inst['Lam1'][1]
2502                    dmin = wave/(2.0*sind(limits[1]/2.0))
2503                elif 'T' in inst['Type'][0]:
2504                    dmin = limits[0]/inst['difC'][1]
2505                else:
2506                    wave = inst['Lam'][1]
2507                    dmin = wave/(2.0*sind(limits[1]/2.0))
2508                pfx = str(pId)+':'+str(hId)+':'
2509                if Phases[phase]['General']['doPawley']:
2510                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2511                    hapDict[pfx+'newLeBail'] = True
2512                    Tmin = G2lat.Dsp2pos(inst,dmin)
2513                    if 'T' in inst['Type'][1]:
2514                        limits[0] = max(limits[0],Tmin)
2515                    else:
2516                        limits[1] = min(limits[1],Tmin)
2517                else:
2518                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2519                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2520                if Phases[phase]['General']['Type'] == 'magnetic':
2521                    dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.))
2522                for item in ['Scale','Extinction']:
2523                    hapDict[pfx+item] = hapData[item][0]
2524                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2525                        hapVary.append(pfx+item)
2526                names = G2spc.HStrainNames(SGData)
2527                HSvals = []
2528                for i,name in enumerate(names):
2529                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2530                    HSvals.append(hapDict[pfx+name])
2531                    if hapData['HStrain'][1][i]:
2532                        hapVary.append(pfx+name)
2533                if 'Layer Disp' in hapData:
2534                    hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0]
2535                    if hapData['Layer Disp'][1]:
2536                        hapVary.append(pfx+'LayerDisp')
2537                else:
2538                    hapDict[pfx+'LayerDisp'] = 0.0
2539                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2540                if hapData['Pref.Ori.'][0] == 'MD':
2541                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2542                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2543                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2544                        hapVary.append(pfx+'MD')
2545                else:                           #'SH' spherical harmonics
2546                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2547                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2548                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2549                    controlDict[pfx+'SHhkl'] = []
2550                    try: #patch for old Pref.Ori. items
2551                        controlDict[pfx+'SHtoler'] = 0.1
2552                        if hapData['Pref.Ori.'][6][0] != '':
2553                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2554                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2555                    except IndexError:
2556                        pass
2557                    for item in hapData['Pref.Ori.'][5]:
2558                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2559                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2560                            hapVary.append(pfx+item)
2561                for item in ['Mustrain','Size']:
2562                    controlDict[pfx+item+'Type'] = hapData[item][0]
2563                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2564                    if hapData[item][2][2]:
2565                        hapVary.append(pfx+item+';mx')
2566                    if hapData[item][0] in ['isotropic','uniaxial']:
2567                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2568                        if hapData[item][2][0]:
2569                            hapVary.append(pfx+item+';i')
2570                        if hapData[item][0] == 'uniaxial':
2571                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2572                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2573                            if hapData[item][2][1]:
2574                                hapVary.append(pfx+item+';a')
2575                    else:       #generalized for mustrain or ellipsoidal for size
2576                        Nterms = len(hapData[item][4])
2577                        if item == 'Mustrain':
2578                            names = G2spc.MustrainNames(SGData)
2579                            pwrs = []
2580                            for name in names:
2581                                h,k,l = name[1:]
2582                                pwrs.append([int(h),int(k),int(l)])
2583                            controlDict[pfx+'MuPwrs'] = pwrs
2584                        for i in range(Nterms):
2585                            sfx = ';'+str(i)
2586                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2587                            if hapData[item][5][i]:
2588                                hapVary.append(pfx+item+sfx)
2589                if Phases[phase]['General']['Type'] != 'magnetic':
2590                    for bab in ['BabA','BabU']:
2591                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2592                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2593                            hapVary.append(pfx+bab)
2594                               
2595                if Print: 
2596                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2597                    pFile.write(135*'='+'\n')
2598                    if hapDict[pfx+'LeBail']:
2599                        pFile.write(' Perform LeBail extraction\n')                     
2600                    else:
2601                        pFile.write(' Phase fraction  : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2602                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2603                        if hapData['Pref.Ori.'][0] == 'MD':
2604                            Ax = hapData['Pref.Ori.'][3]
2605                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2606                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2607                        else: #'SH' for spherical harmonics
2608                            PrintSHPO(hapData['Pref.Ori.'])
2609                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2610                    PrintSize(hapData['Size'])
2611                    PrintMuStrain(hapData['Mustrain'],SGData)
2612                    PrintHStrain(hapData['HStrain'],SGData)
2613                    if 'Layer Disp' in hapData:
2614                        pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1]))
2615                    if Phases[phase]['General']['Type'] != 'magnetic':
2616                        if hapData['Babinet']['BabA'][0]:
2617                            PrintBabinet(hapData['Babinet'])
2618                if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2619                    hapData['newLeBail'] = True
2620                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2621                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2622                        hapData['newLeBail'] = False
2623                    refList = []
2624#                    Uniq = []
2625#                    Phi = []
2626                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2627                    if Phases[phase]['General'].get('Modulated',False):
2628                        ifSuper = True
2629                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2630                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2631                        for h,k,l,m,d in HKLd:
2632                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2633                            mul *= 2      # for powder overlap of Friedel pairs
2634                            if m or not ext or useExt:
2635                                if 'C' in inst['Type'][0]:
2636                                    pos = G2lat.Dsp2pos(inst,d)
2637                                    if limits[0] < pos < limits[1]:
2638                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,1., 0.0,0.0,1.0,1.0,1.0])
2639                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2640#                                        Uniq.append(uniq)
2641#                                        Phi.append(phi)
2642                                elif 'T' in inst['Type'][0]:
2643                                    pos = G2lat.Dsp2pos(inst,d)
2644                                    if limits[0] < pos < limits[1]:
2645                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2646                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,1., 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
2647                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2648                                        #TODO - if tabulated put alp & bet in here
2649#                                        Uniq.append(uniq)
2650#                                        Phi.append(phi)
2651                                elif 'B' in inst['Type'][0]:
2652                                    pos = G2lat.Dsp2pos(inst,d)
2653                                    if limits[0] < pos < limits[1]:
2654                                        refList.append([h,k,l,m,mul,d, pos,0.0,0.0,0.0,1., 0.0,0.0,0.0,0.0, 1.0,1.0,1.0])
2655                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2656                    else:
2657                        ifSuper = False
2658                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2659                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2660                        for h,k,l,d in HKLd:
2661                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2662                            if ext and 'N' in inst['Type'][0] and  'MagSpGrp' in SGData:
2663                                ext = G2spc.checkMagextc([h,k,l],SGData)
2664                            mul *= 2      # for powder overlap of Friedel pairs
2665                            if ext and not useExt:
2666                                continue
2667                            if 'C' in inst['Type'][0]:
2668                                pos = G2lat.Dsp2pos(inst,d)
2669                                if limits[0] < pos < limits[1]:
2670                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,1., 0.0,0.0,1.0,1.0,1.0])
2671                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2672#                                    Uniq.append(uniq)
2673#                                    Phi.append(phi)
2674                            elif 'T' in inst['Type'][0]:
2675                                pos = G2lat.Dsp2pos(inst,d)
2676                                if limits[0] < pos < limits[1]:
2677                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2678                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,1., 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
2679                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2680#                                    Uniq.append(uniq)
2681#                                    Phi.append(phi)
2682                            elif 'B' in inst['Type'][0]:
2683                                pos = G2lat.Dsp2pos(inst,d)
2684                                if limits[0] < pos < limits[1]:
2685                                    refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,1., 0.0,0.0,0.0,0.0, 1.0,1.0,1.0])
2686                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2687                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2688            elif 'HKLF' in histogram:
2689                inst = Histogram['Instrument Parameters'][0]
2690                hId = Histogram['hId']
2691                hfx = ':%d:'%(hId)
2692                for item in inst:
2693                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2694                    hapDict[hfx+item] = inst[item][1]
2695                pfx = str(pId)+':'+str(hId)+':'
2696                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2697                if hapData['Scale'][1]:
2698                    hapVary.append(pfx+'Scale')
2699                               
2700                extApprox,extType,extParms = hapData['Extinction']
2701                controlDict[pfx+'EType'] = extType
2702                controlDict[pfx+'EApprox'] = extApprox
2703                if 'C' in inst['Type'][0]:
2704                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2705                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2706                if 'Primary' in extType:
2707                    Ekey = ['Ep',]
2708                elif 'I & II' in extType:
2709                    Ekey = ['Eg','Es']
2710                elif 'Secondary Type II' == extType:
2711                    Ekey = ['Es',]
2712                elif 'Secondary Type I' == extType:
2713                    Ekey = ['Eg',]
2714                else:   #'None'
2715                    Ekey = []
2716                for eKey in Ekey:
2717                    hapDict[pfx+eKey] = extParms[eKey][0]
2718                    if extParms[eKey][1]:
2719                        hapVary.append(pfx+eKey)
2720                for bab in ['BabA','BabU']:
2721                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2722                    if hapData['Babinet'][bab][1]:
2723                        hapVary.append(pfx+bab)
2724                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2725                if len(Twins) == 1:
2726                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2727                    if hapData.get('Flack',[0,False])[1]:
2728                        hapVary.append(pfx+'Flack')
2729                sumTwFr = 0.
2730                controlDict[pfx+'TwinLaw'] = []
2731                controlDict[pfx+'TwinInv'] = []
2732                NTL = 0           
2733                for it,twin in enumerate(Twins):
2734                    if 'bool' in str(type(twin[0])):
2735                        controlDict[pfx+'TwinInv'].append(twin[0])
2736                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2737                    else:
2738                        NTL += 1
2739                        controlDict[pfx+'TwinInv'].append(False)
2740                        controlDict[pfx+'TwinLaw'].append(twin[0])
2741                    if it:
2742                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2743                        sumTwFr += twin[1]
2744                    else:
2745                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2746                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2747                    if Twins[0][1][1]:
2748                        hapVary.append(pfx+'TwinFr:'+str(it))
2749                controlDict[pfx+'NTL'] = NTL
2750                #need to make constraint on TwinFr
2751                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2752                if len(Twins) > 1:    #force sum to unity
2753                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2754                if Print: 
2755                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2756                    pFile.write(135*'='+'\n')
2757                    pFile.write(' Scale factor     : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2758                    if extType != 'None':
2759                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2760                        text = ' Parameters       :'
2761                        for eKey in Ekey:
2762                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2763                        pFile.write(text+'\n')
2764                    if hapData['Babinet']['BabA'][0]:
2765                        PrintBabinet(hapData['Babinet'])
2766                    if not SGData['SGInv'] and len(Twins) == 1:
2767                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2768                    if len(Twins) > 1:
2769                        for it,twin in enumerate(Twins):
2770                            if 'bool' in str(type(twin[0])):
2771                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2772                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2773                            else:
2774                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2775                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2776                       
2777                Histogram['Reflection Lists'] = phase       
2778               
2779    return hapVary,hapDict,controlDict
2780   
2781def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,Print=True,pFile=None):
2782    'needs a doc string'
2783   
2784    def PrintSizeAndSig(hapData,sizeSig):
2785        line = '\n Size model:     %9s'%(hapData[0])
2786        refine = False
2787        if hapData[0] in ['isotropic','uniaxial']:
2788            line += ' equatorial:%12.5f'%(hapData[1][0])
2789            if sizeSig[0][0]:
2790                line += ', sig:%8.4f'%(sizeSig[0][0])
2791                refine = True
2792            if hapData[0] == 'uniaxial':
2793                line += ' axial:%12.4f'%(hapData[1][1])
2794                if sizeSig[0][1]:
2795                    refine = True
2796                    line += ', sig:%8.4f'%(sizeSig[0][1])
2797            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2798            if sizeSig[0][2]:
2799                refine = True
2800                line += ', sig:%8.4f'%(sizeSig[0][2])
2801            if refine:
2802                pFile.write(line+'\n')
2803        else:
2804            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2805            if sizeSig[0][2]:
2806                refine = True
2807                line += ', sig:%8.4f'%(sizeSig[0][2])
2808            Snames = ['S11','S22','S33','S12','S13','S23']
2809            ptlbls = ' name  :'
2810            ptstr =  ' value :'
2811            sigstr = ' sig   :'
2812            for i,name in enumerate(Snames):
2813                ptlbls += '%12s' % (name)
2814                ptstr += '%12.6f' % (hapData[4][i])
2815                if sizeSig[1][i]:
2816                    refine = True
2817                    sigstr += '%12.6f' % (sizeSig[1][i])
2818                else:
2819                    sigstr += 12*' '
2820            if refine:
2821                pFile.write(line+'\n')
2822                pFile.write(ptlbls+'\n')
2823                pFile.write(ptstr+'\n')
2824                pFile.write(sigstr+'\n')
2825       
2826    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2827        line = '\n Mustrain model: %9s\n'%(hapData[0])
2828        refine = False
2829        if hapData[0] in ['isotropic','uniaxial']:
2830            line += ' equatorial:%12.1f'%(hapData[1][0])
2831            if mustrainSig[0][0]:
2832                line += ', sig:%8.1f'%(mustrainSig[0][0])
2833                refine = True
2834            if hapData[0] == 'uniaxial':
2835                line += ' axial:%12.1f'%(hapData[1][1])
2836                if mustrainSig[0][1]:
2837                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2838            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2839            if mustrainSig[0][2]:
2840                refine = True
2841                line += ', sig:%8.3f'%(mustrainSig[0][2])
2842            if refine:
2843                pFile.write(line+'\n')
2844        else:
2845            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2846            if mustrainSig[0][2]:
2847                refine = True
2848                line += ', sig:%8.3f'%(mustrainSig[0][2])
2849            Snames = G2spc.MustrainNames(SGData)
2850            ptlbls = ' name  :'
2851            ptstr =  ' value :'
2852            sigstr = ' sig   :'
2853            for i,name in enumerate(Snames):
2854                ptlbls += '%12s' % (name)
2855                ptstr += '%12.1f' % (hapData[4][i])
2856                if mustrainSig[1][i]:
2857                    refine = True
2858                    sigstr += '%12.1f' % (mustrainSig[1][i])
2859                else:
2860                    sigstr += 12*' '
2861            if refine:
2862                pFile.write(line+'\n')
2863                pFile.write(ptlbls+'\n')
2864                pFile.write(ptstr+'\n')
2865                pFile.write(sigstr+'\n')
2866           
2867    def PrintHStrainAndSig(hapData,strainSig,SGData):
2868        Hsnames = G2spc.HStrainNames(SGData)
2869        ptlbls = ' name  :'
2870        ptstr =  ' value :'
2871        sigstr = ' sig   :'
2872        refine = False
2873        for i,name in enumerate(Hsnames):
2874            ptlbls += '%12s' % (name)
2875            ptstr += '%12.4g' % (hapData[0][i])
2876            if name in strainSig:
2877                refine = True
2878                sigstr += '%12.4g' % (strainSig[name])
2879            else:
2880                sigstr += 12*' '
2881        if refine:
2882            pFile.write('\n Hydrostatic/elastic strain:\n')
2883            pFile.write(ptlbls+'\n')
2884            pFile.write(ptstr+'\n')
2885            pFile.write(sigstr+'\n')
2886       
2887    def PrintSHPOAndSig(pfx,hapData,POsig):
2888        Tindx = 1.0
2889        Tvar = 0.0
2890        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2891        ptlbls = ' names :'
2892        ptstr =  ' values:'
2893        sigstr = ' sig   :'
2894        for item in hapData[5]:
2895            ptlbls += '%12s'%(item)
2896            ptstr += '%12.3f'%(hapData[5][item])
2897            l = 2.0*eval(item.strip('C'))[0]+1
2898            Tindx += hapData[5][item]**2/l
2899            if pfx+item in POsig:
2900                Tvar += (2.*hapData[5][item]*POsig[pfx+item]/l)**2
2901                sigstr += '%12.3f'%(POsig[pfx+item])
2902            else:
2903                sigstr += 12*' ' 
2904        pFile.write(ptlbls+'\n')
2905        pFile.write(ptstr+'\n')
2906        pFile.write(sigstr+'\n')
2907        pFile.write('\n Texture index J = %.3f(%d)\n'%(Tindx,int(1000*np.sqrt(Tvar))))
2908       
2909    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2910        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2911        text = ''
2912        for item in hapData[2]:
2913            if pfx+item in ScalExtSig:
2914                text += '       %s: '%(item)
2915                text += '%12.2e'%(hapData[2][item][0])
2916                if pfx+item in ScalExtSig:
2917                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2918        pFile.write(text+'\n')   
2919       
2920    def PrintBabinetAndSig(pfx,hapData,BabSig):
2921        pFile.write('\n Babinet form factor modification:\n')
2922        ptlbls = ' names :'
2923        ptstr =  ' values:'
2924        sigstr = ' sig   :'
2925        for item in hapData:
2926            ptlbls += '%12s'%(item)
2927            ptstr += '%12.3f'%(hapData[item][0])
2928            if pfx+item in BabSig:
2929                sigstr += '%12.3f'%(BabSig[pfx+item])
2930            else:
2931                sigstr += 12*' ' 
2932        pFile.write(ptlbls+'\n')
2933        pFile.write(ptstr+'\n')
2934        pFile.write(sigstr+'\n')
2935       
2936    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2937        pFile.write('\n Twin Law fractions :\n')
2938        ptlbls = ' names :'
2939        ptstr =  ' values:'
2940        sigstr = ' sig   :'
2941        for it,item in enumerate(twinData):
2942            ptlbls += '     twin #%d'%(it)
2943            if it:
2944                ptstr += '%12.3f'%(item[1])
2945            else:
2946                ptstr += '%12.3f'%(item[1][0])
2947            if pfx+'TwinFr:'+str(it) in TwinSig:
2948                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2949            else:
2950                sigstr += 12*' ' 
2951        pFile.write(ptlbls+'\n')
2952        pFile.write(ptstr+'\n')
2953        pFile.write(sigstr+'\n')
2954       
2955   
2956    PhFrExtPOSig = {}
2957    SizeMuStrSig = {}
2958    ScalExtSig = {}
2959    BabSig = {}
2960    TwinFrSig = {}
2961    wtFrSum = {}
2962    for phase in Phases:
2963        HistoPhase = Phases[phase]['Histograms']
2964        General = Phases[phase]['General']
2965        SGData = General['SGData']
2966        pId = Phases[phase]['pId']
2967        histoList = list(HistoPhase.keys())
2968        histoList.sort()
2969        for histogram in histoList:
2970            try:
2971                Histogram = Histograms[histogram]
2972            except KeyError:                       
2973                #skip if histogram not included e.g. in a sequential refinement
2974                continue
2975            if not Phases[phase]['Histograms'][histogram]['Use']:
2976                #skip if phase absent from this histogram
2977                continue
2978            hapData = HistoPhase[histogram]
2979            hId = Histogram['hId']
2980            pfx = str(pId)+':'+str(hId)+':'
2981            if hId not in wtFrSum:
2982                wtFrSum[hId] = 0.
2983            if 'PWDR' in histogram:
2984                for item in ['Scale','Extinction']:
2985                    hapData[item][0] = parmDict[pfx+item]
2986                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2987                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2988                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2989                if hapData['Pref.Ori.'][0] == 'MD':
2990                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2991                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2992                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2993                else:                           #'SH' spherical harmonics
2994                    for item in hapData['Pref.Ori.'][5]:
2995                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2996                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2997                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2998                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2999                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
3000                    pfx+'HStrain':{}})                 
3001                for item in ['Mustrain','Size']:
3002                    hapData[item][1][2] = parmDict[pfx+item+';mx']
3003#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
3004                    if pfx+item+';mx' in sigDict:
3005                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
3006                    if hapData[item][0] in ['isotropic','uniaxial']:                   
3007                        hapData[item][1][0] = parmDict[pfx+item+';i']
3008                        if item == 'Size':
3009                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
3010                        if pfx+item+';i' in sigDict: 
3011                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
3012                        if hapData[item][0] == 'uniaxial':
3013                            hapData[item][1][1] = parmDict[pfx+item+';a']
3014                            if item == 'Size':
3015                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
3016                            if pfx+item+';a' in sigDict:
3017                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
3018                    else:       #generalized for mustrain or ellipsoidal for size
3019                        Nterms = len(hapData[item][4])
3020                        for i in range(Nterms):
3021                            sfx = ';'+str(i)
3022                            hapData[item][4][i] = parmDict[pfx+item+sfx]
3023                            if pfx+item+sfx in sigDict:
3024                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
3025                names = G2spc.HStrainNames(SGData)
3026                for i,name in enumerate(names):
3027                    hapData['HStrain'][0][i] = parmDict[pfx+name]
3028                    if pfx+name in sigDict:
3029                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
3030                if 'Layer Disp' in hapData:
3031                    hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp']
3032                    if pfx+'LayerDisp' in sigDict:
3033                        SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp']
3034                if Phases[phase]['General']['Type'] != 'magnetic':
3035                    for name in ['BabA','BabU']:
3036                        hapData['Babinet'][name][0] = parmDict[pfx+name]
3037                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
3038                            BabSig[pfx+name] = sigDict[pfx+name]               
3039               
3040            elif 'HKLF' in histogram:
3041                for item in ['Scale','Flack']:
3042                    if parmDict.get(pfx+item):
3043                        hapData[item][0] = parmDict[pfx+item]
3044                        if pfx+item in sigDict:
3045                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3046                for item in ['Ep','Eg','Es']:
3047                    if parmDict.get(pfx+item):
3048                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
3049                        if pfx+item in sigDict:
3050                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3051                for item in ['BabA','BabU']:
3052                    hapData['Babinet'][item][0] = parmDict[pfx+item]
3053                    if pfx+item in sigDict:
3054                        BabSig[pfx+item] = sigDict[pfx+item]
3055                if 'Twins' in hapData:
3056                    it = 1
3057                    sumTwFr = 0.
3058                    while True:
3059                        try:
3060                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
3061                            if pfx+'TwinFr:'+str(it) in sigDict:
3062                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
3063                            if it:
3064                                sumTwFr += hapData['Twins'][it][1]
3065                            it += 1
3066                        except KeyError:
3067                            break
3068                    hapData['Twins'][0][1][0] = 1.-sumTwFr
3069
3070    if Print:
3071        for phase in Phases:
3072            HistoPhase = Phases[phase]['Histograms']
3073            General = Phases[phase]['General']
3074            SGData = General['SGData']
3075            pId = Phases[phase]['pId']
3076            histoList = list(HistoPhase.keys())
3077            histoList.sort()
3078            for histogram in histoList:
3079                try:
3080                    Histogram = Histograms[histogram]
3081                except KeyError:                       
3082                    #skip if histogram not included e.g. in a sequential refinement
3083                    continue
3084                hapData = HistoPhase[histogram]
3085                hId = Histogram['hId']
3086                Histogram['Residuals'][str(pId)+'::Name'] = phase
3087                pfx = str(pId)+':'+str(hId)+':'
3088                hfx = ':%s:'%(hId)
3089                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
3090                    continue
3091                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
3092                pFile.write(135*'='+'\n')
3093                Inst = Histogram['Instrument Parameters'][0]
3094                if 'PWDR' in histogram:
3095                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
3096                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
3097                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
3098                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
3099                   
3100                    if parmDict[pfx+'LeBail']:
3101                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
3102                    else:
3103                        # if calcControls != None:    #skipped in seqRefine
3104                        #     if 'X'in Inst['Type'][0]:
3105                        #         PrintFprime(calcControls['FFtables'],hfx,pFile)
3106                        #     elif 'NC' in Inst['Type'][0]:
3107                        #         PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3108                        if pfx+'Scale' in PhFrExtPOSig:
3109                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
3110                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
3111                            pFile.write(' Phase fraction  : %10.5g, sig %10.5g Weight fraction  : %8.5f, sig %10.5f\n'%
3112                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
3113                        if pfx+'Extinction' in PhFrExtPOSig:
3114                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
3115                        if hapData['Pref.Ori.'][0] == 'MD':
3116                            if pfx+'MD' in PhFrExtPOSig:
3117                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
3118                        else:
3119                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
3120                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
3121                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
3122                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
3123                    if pfx+'LayerDisp' in SizeMuStrSig:
3124                        pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp']))           
3125                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
3126                        if len(BabSig):
3127                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3128                   
3129                elif 'HKLF' in histogram:
3130                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
3131                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
3132                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
3133                    if calcControls != None:    #skipped in seqRefine
3134                        if 'X'in Inst['Type'][0]:
3135                            PrintFprime(calcControls['FFtables'],hfx,pFile)
3136                        elif 'NC' in Inst['Type'][0]:
3137                            PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3138                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
3139                    if pfx+'Scale' in ScalExtSig:
3140                        pFile.write(' Scale factor : %10.4g, sig %10.4g\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
3141                    if hapData['Extinction'][0] != 'None':
3142                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
3143                    if len(BabSig):
3144                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3145                    if pfx+'Flack' in ScalExtSig:
3146                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
3147                    if len(TwinFrSig):
3148                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
3149
3150################################################################################
3151##### Histogram data
3152################################################################################       
3153                   
3154def GetHistogramData(Histograms,Print=True,pFile=None):
3155    'needs a doc string'
3156   
3157    def GetBackgroundParms(hId,Background):
3158        Back = Background[0]
3159        DebyePeaks = Background[1]
3160        bakType,bakFlag = Back[:2]
3161        backVals = Back[3:]
3162        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
3163        backDict = dict(zip(backNames,backVals))
3164        backVary = []
3165        if bakFlag:
3166            backVary = backNames
3167        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
3168        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
3169        debyeDict = {}
3170        debyeList = []
3171        for i in range(DebyePeaks['nDebye']):
3172            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
3173            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
3174            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
3175        debyeVary = []
3176        for item in debyeList:
3177            if item[1]:
3178                debyeVary.append(item[0])
3179        backDict.update(debyeDict)
3180        backVary += debyeVary
3181        peakDict = {}
3182        peakList = []
3183        for i in range(DebyePeaks['nPeaks']):
3184            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
3185                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
3186            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
3187            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
3188        peakVary = []
3189        for item in peakList:
3190            if item[1]:
3191                peakVary.append(item[0])
3192        backDict.update(peakDict)
3193        backVary += peakVary
3194        if 'background PWDR' in Background[1]:
3195            backDict[':'+str(hId)+':Back File'] = Background[1]['background PWDR'][0]
3196            backDict[':'+str(hId)+':BF mult'] = Background[1]['background PWDR'][1]
3197            try:
3198                if Background[1]['background PWDR'][0] and Background[1]['background PWDR'][2]:
3199                    backVary.append(':'+str(hId)+':BF mult')
3200            except IndexError:  # old version without refine flag
3201                pass
3202        return bakType,backDict,backVary           
3203       
3204    def GetInstParms(hId,Inst):
3205        #patch
3206        if 'Z' not in Inst:
3207            Inst['Z'] = [0.0,0.0,False]
3208        dataType = Inst['Type'][0]
3209        instDict = {}
3210        insVary = []
3211        pfx = ':'+str(hId)+':'
3212        insKeys = list(Inst.keys())
3213        insKeys.sort()
3214        for item in insKeys:
3215            if item in 'Azimuth':
3216                continue
3217            insName = pfx+item
3218            instDict[insName] = Inst[item][1]
3219            if len(Inst[item]) > 2 and Inst[item][2]:
3220                insVary.append(insName)
3221        if 'C' in dataType:
3222            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
3223        elif 'T' in dataType:   #trap zero alp, bet coeff.
3224            if not instDict[pfx+'alpha']:
3225                instDict[pfx+'alpha'] = 1.0
3226            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3227                instDict[pfx+'beta-1'] = 1.0
3228        elif 'B' in dataType:   #trap zero alp, bet coeff.
3229            if not instDict[pfx+'alpha-0'] and not instDict[pfx+'alpha-1']:
3230                instDict[pfx+'alpha-1'] = 1.0
3231            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3232                instDict[pfx+'beta-1'] = 1.0
3233        return dataType,instDict,insVary
3234       
3235    def GetSampleParms(hId,Sample):
3236        sampVary = []
3237        hfx = ':'+str(hId)+':'       
3238        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
3239            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi'],hfx+'Azimuth':Sample['Azimuth']}
3240        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
3241            if key in Sample:
3242                sampDict[hfx+key] = Sample[key]
3243        Type = Sample['Type']
3244        if 'Bragg' in Type:             #Bragg-Brentano
3245            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3246                sampDict[hfx+item] = Sample[item][0]
3247                if Sample[item][1]:
3248                    sampVary.append(hfx+item)
3249        elif 'Debye' in Type:        #Debye-Scherrer
3250            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3251                sampDict[hfx+item] = Sample[item][0]
3252                if Sample[item][1]:
3253                    sampVary.append(hfx+item)
3254        return Type,sampDict,sampVary
3255       
3256    def PrintBackground(Background):
3257        Back = Background[0]
3258        DebyePeaks = Background[1]
3259        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
3260        line = ' Coefficients: '
3261        for i,back in enumerate(Back[3:]):
3262            line += '%10.3f'%(back)
3263            if i and not i%10:
3264                line += '\n'+15*' '
3265        pFile.write(line+'\n')
3266        if DebyePeaks['nDebye']:
3267            pFile.write('\n Debye diffuse scattering coefficients\n')
3268            parms = ['DebyeA','DebyeR','DebyeU']
3269            line = ' names :  '
3270            for parm in parms:
3271                line += '%8s refine?'%(parm)
3272            pFile.write(line+'\n')
3273            for j,term in enumerate(DebyePeaks['debyeTerms']):
3274                line = ' term'+'%2d'%(j)+':'
3275                for i in range(3):
3276                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3277                pFile.write(line+'\n')
3278        if DebyePeaks['nPeaks']:
3279            pFile.write('\n Single peak coefficients\n')
3280            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3281            line = ' names :  '
3282            for parm in parms:
3283                line += '%8s refine?'%(parm)
3284            pFile.write(line+'\n')
3285            for j,term in enumerate(DebyePeaks['peaksList']):
3286                line = ' peak'+'%2d'%(j)+':'
3287                for i in range(4):
3288                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3289                pFile.write(line+'\n')
3290        if 'background PWDR' in DebyePeaks:
3291            try:
3292                pFile.write(' Fixed background file: %s; mult: %.3f  Refine? %s\n'%(DebyePeaks['background PWDR'][0],
3293                    DebyePeaks['background PWDR'][1],DebyePeaks['background PWDR'][2]))
3294            except IndexError:  #old version without refine flag
3295                pass
3296       
3297    def PrintInstParms(Inst):
3298        pFile.write('\n Instrument Parameters:\n')
3299        insKeys = list(Inst.keys())
3300        insKeys.sort()
3301        iBeg = 0
3302        Ok = True
3303        while Ok:
3304            ptlbls = ' name  :'
3305            ptstr =  ' value :'
3306            varstr = ' refine:'
3307            iFin = min(iBeg+9,len(insKeys))
3308            for item in insKeys[iBeg:iFin]:
3309                if item not in ['Type','Source']:
3310                    ptlbls += '%12s' % (item)
3311                    ptstr += '%12.6f' % (Inst[item][1])
3312                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3313                        varstr += 12*' '
3314                    else:
3315                        varstr += '%12s' % (str(bool(Inst[item][2])))
3316            pFile.write(ptlbls+'\n')
3317            pFile.write(ptstr+'\n')
3318            pFile.write(varstr+'\n')
3319            iBeg = iFin
3320            if iBeg == len(insKeys):
3321                Ok = False
3322            else:
3323                pFile.write('\n')
3324       
3325    def PrintSampleParms(Sample):
3326        pFile.write('\n Sample Parameters:\n')
3327        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3328            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3329        ptlbls = ' name  :'
3330        ptstr =  ' value :'
3331        varstr = ' refine:'
3332        if 'Bragg' in Sample['Type']:
3333            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3334                ptlbls += '%14s'%(item)
3335                ptstr += '%14.4f'%(Sample[item][0])
3336                varstr += '%14s'%(str(bool(Sample[item][1])))
3337           
3338        elif 'Debye' in Type:        #Debye-Scherrer
3339            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3340                ptlbls += '%14s'%(item)
3341                ptstr += '%14.4f'%(Sample[item][0])
3342                varstr += '%14s'%(str(bool(Sample[item][1])))
3343
3344        pFile.write(ptlbls+'\n')
3345        pFile.write(ptstr+'\n')
3346        pFile.write(varstr+'\n')
3347       
3348    histDict = {}
3349    histVary = []
3350    controlDict = {}
3351    histoList = list(Histograms.keys())
3352    histoList.sort()
3353    for histogram in histoList:
3354        if 'PWDR' in histogram:
3355            Histogram = Histograms[histogram]
3356            hId = Histogram['hId']
3357            pfx = ':'+str(hId)+':'
3358            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3359            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3360            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3361            for excl in controlDict[pfx+'Exclude']:
3362                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3363            if controlDict[pfx+'Exclude']:
3364                ma.mask_rows(Histogram['Data'])
3365            Background = Histogram['Background']
3366            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3367            controlDict[pfx+'bakType'] = Type
3368            histDict.update(bakDict)
3369            histVary += bakVary
3370           
3371            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3372            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3373                G2fil.G2Print ('Warning: tabulated profile coefficients are ignored')
3374            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3375            controlDict[pfx+'histType'] = Type
3376            if 'XC' in Type or 'XB' in Type:
3377                if pfx+'Lam1' in instDict:
3378                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3379                else:
3380                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3381            histDict.update(instDict)
3382            histVary += insVary
3383           
3384            Sample = Histogram['Sample Parameters']
3385            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3386            controlDict[pfx+'instType'] = Type
3387            histDict.update(sampDict)
3388            histVary += sampVary
3389           
3390   
3391            if Print: 
3392                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3393                pFile.write(135*'='+'\n')
3394                Units = {'C':' deg','T':' msec','B':' deg'}
3395                units = Units[controlDict[pfx+'histType'][2]]
3396                Limits = controlDict[pfx+'Limits']
3397                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3398                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3399                if len(controlDict[pfx+'Exclude']):
3400                    excls = controlDict[pfx+'Exclude']
3401                    for excl in excls:
3402                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3403                PrintSampleParms(Sample)
3404                PrintInstParms(Inst[0])
3405                PrintBackground(Background)
3406        elif 'HKLF' in histogram:
3407            Histogram = Histograms[histogram]
3408            hId = Histogram['hId']
3409            pfx = ':'+str(hId)+':'
3410            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3411            Inst = Histogram['Instrument Parameters'][0]
3412            controlDict[pfx+'histType'] = Inst['Type'][0]
3413            if 'X' in Inst['Type'][0]:
3414                histDict[pfx+'Lam'] = Inst['Lam'][1]
3415                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3416            elif 'NC' in Inst['Type'][0] or 'NB' in Inst['Type'][0]:                   
3417                histDict[pfx+'Lam'] = Inst['Lam'][1]
3418    return histVary,histDict,controlDict
3419   
3420def SetHistogramData(parmDict,sigDict,Histograms,calcControls,Print=True,pFile=None,seq=False):
3421    'Shows histogram data after a refinement'
3422   
3423    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3424        Back = Background[0]
3425        DebyePeaks = Background[1]
3426        lenBack = len(Back[3:])
3427        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3428        for i in range(lenBack):
3429            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3430            if pfx+'Back;'+str(i) in sigDict:
3431                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3432        if DebyePeaks['nDebye']:
3433            for i in range(DebyePeaks['nDebye']):
3434                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3435                for j,name in enumerate(names):
3436                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3437                    if name in sigDict:
3438                        backSig[lenBack+3*i+j] = sigDict[name]           
3439        if DebyePeaks['nPeaks']:
3440            for i in range(DebyePeaks['nPeaks']):
3441                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3442                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3443                for j,name in enumerate(names):
3444                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3445                    if name in sigDict:
3446                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3447        if pfx+'BF mult' in sigDict:
3448            DebyePeaks['background PWDR'][1] = parmDict[pfx+'BF mult']
3449            backSig.append(sigDict[pfx+'BF mult'])
3450               
3451        return backSig
3452       
3453    def SetInstParms(pfx,Inst,parmDict,sigDict):
3454        instSig = {}
3455        insKeys = list(Inst.keys())
3456        insKeys.sort()
3457        for item in insKeys:
3458            insName = pfx+item
3459            Inst[item][1] = parmDict[insName]
3460            if insName in sigDict:
3461                instSig[item] = sigDict[insName]
3462            else:
3463                instSig[item] = 0
3464        return instSig
3465       
3466    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3467        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3468            sampSig = [0 for i in range(5)]
3469            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3470                Sample[item][0] = parmDict[pfx+item]
3471                if pfx+item in sigDict:
3472                    sampSig[i] = sigDict[pfx+item]
3473        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3474            sampSig = [0 for i in range(4)]
3475            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3476                Sample[item][0] = parmDict[pfx+item]
3477                if pfx+item in sigDict:
3478                    sampSig[i] = sigDict[pfx+item]
3479        return sampSig
3480       
3481    def PrintBackgroundSig(Background,backSig):
3482        Back = Background[0]
3483        DebyePeaks = Background[1]
3484        valstr = ' value : '
3485        sigstr = ' sig   : '
3486        refine = False
3487        for i,back in enumerate(Back[3:]):
3488            valstr += '%10.4g'%(back)
3489            if Back[1]:
3490                refine = True
3491                sigstr += '%10.4g'%(backSig[i])
3492            else:
3493                sigstr += 10*' '
3494        if refine:
3495            pFile.write('\n Background function: %s\n'%Back[0])
3496            pFile.write(valstr+'\n')
3497            pFile.write(sigstr+'\n')
3498        if DebyePeaks['nDebye']:
3499            ifAny = False
3500            ptfmt = "%12.3f"
3501            names =  ' names :'
3502            ptstr =  ' values:'
3503            sigstr = ' esds  :'
3504            for item in sigDict:
3505                if 'Debye' in item:
3506                    ifAny = True
3507                    names += '%12s'%(item)
3508                    ptstr += ptfmt%(parmDict[item])
3509                    sigstr += ptfmt%(sigDict[item])
3510            if ifAny:
3511                pFile.write('\n Debye diffuse scattering coefficients\n')
3512                pFile.write(names+'\n')
3513                pFile.write(ptstr+'\n')
3514                pFile.write(sigstr+'\n')
3515        if DebyePeaks['nPeaks']:
3516            pFile.write('\n Single peak coefficients:\n')
3517            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3518            line = ' peak no. '
3519            for parm in parms:
3520                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3521            pFile.write(line+'\n')
3522            for ip in range(DebyePeaks['nPeaks']):
3523                ptstr = ' %4d '%(ip)
3524                for parm in parms:
3525                    fmt = '%14.3f'
3526                    efmt = '%12.3f'
3527                    if parm == 'BkPkpos':
3528                        fmt = '%14.4f'
3529                        efmt = '%12.4f'
3530                    name = pfx+parm+';%d'%(ip)
3531                    ptstr += fmt%(parmDict[name])
3532                    if name in sigDict:
3533                        ptstr += efmt%(sigDict[name])
3534                    else:
3535                        ptstr += 12*' '
3536                pFile.write(ptstr+'\n')
3537        if 'background PWDR' in DebyePeaks and DebyePeaks['background PWDR'][2]:
3538            pFile.write(' Fixed background scale: %.3f(%d)\n'%(DebyePeaks['background PWDR'][1],int(1000*backSig[-1])))
3539        sumBk = np.array(Histogram['sumBk'])
3540        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3541            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3542       
3543    def PrintInstParmsSig(Inst,instSig):
3544        refine = False
3545        insKeys = list(instSig.keys())
3546        insKeys.sort()
3547        iBeg = 0
3548        Ok = True
3549        while Ok:
3550            ptlbls = ' names :'
3551            ptstr =  ' value :'
3552            sigstr = ' sig   :'
3553            iFin = min(iBeg+9,len(insKeys))
3554            for name in insKeys[iBeg:iFin]:
3555                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3556                    ptlbls += '%12s' % (name)
3557                    ptstr += '%12.6f' % (Inst[name][1])
3558                    if instSig[name]:
3559                        refine = True
3560                        sigstr += '%12.6f' % (instSig[name])
3561                    else:
3562                        sigstr += 12*' '
3563            if refine:
3564                pFile.write('\n Instrument Parameters:\n')
3565                pFile.write(ptlbls+'\n')
3566                pFile.write(ptstr+'\n')
3567                pFile.write(sigstr+'\n')
3568            iBeg = iFin
3569            if iBeg == len(insKeys):
3570                Ok = False
3571       
3572    def PrintSampleParmsSig(Sample,sampleSig):
3573        ptlbls = ' names :'
3574        ptstr =  ' values:'
3575        sigstr = ' sig   :'
3576        refine = False
3577        if 'Bragg' in Sample['Type']:
3578            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3579                ptlbls += '%14s'%(item)
3580                ptstr += '%14.4f'%(Sample[item][0])
3581                if sampleSig[i]:
3582                    refine = True
3583                    sigstr += '%14.4f'%(sampleSig[i])
3584                else:
3585                    sigstr += 14*' '
3586           
3587        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3588            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3589                ptlbls += '%14s'%(item)
3590                ptstr += '%14.4f'%(Sample[item][0])
3591                if sampleSig[i]:
3592                    refine = True
3593                    sigstr += '%14.4f'%(sampleSig[i])
3594                else:
3595                    sigstr += 14*' '
3596
3597        if refine:
3598            pFile.write('\n Sample Parameters:\n')
3599            pFile.write(ptlbls+'\n')
3600            pFile.write(ptstr+'\n')
3601            pFile.write(sigstr+'\n')
3602       
3603    histoList = list(Histograms.keys())
3604    histoList.sort()
3605    for histogram in histoList:
3606        if 'PWDR' in histogram:
3607            Histogram = Histograms[histogram]
3608            hId = Histogram['hId']
3609            pfx = ':'+str(hId)+':'
3610            Background = Histogram['Background']
3611            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3612           
3613            Inst = Histogram['Instrument Parameters'][0]
3614            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3615       
3616            Sample = Histogram['Sample Parameters']
3617            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3618
3619            if Print and not seq:
3620                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3621                pFile.write(135*'='+'\n')
3622            if Print:
3623                pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3624                pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3625                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3626                pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3627                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3628                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3629                if calcControls != None:    #skipped in seqRefine
3630                    if 'X' in Inst['Type'][0]:
3631                        PrintFprime(calcControls['FFtables'],pfx,pFile)
3632                    elif 'NC' in Inst['Type'][0]:
3633                        PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3634                PrintSampleParmsSig(Sample,sampSig)
3635                PrintInstParmsSig(Inst,instSig)
3636                PrintBackgroundSig(Background,backSig)
3637               
3638def WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
3639    '''Cribbed version of PrintRBObjPOAndSig but returns lists of strings.
3640    Moved so it can be used in ExportCIF
3641    '''
3642    namstr = '  names :'
3643    valstr = '  values:'
3644    sigstr = '  esds  :'
3645    for i,px in enumerate(['Px:','Py:','Pz:']):
3646        name = pfx+rbfx+px+rbsx
3647        namstr += '%12s'%('Pos '+px[1])
3648        valstr += '%12.5f'%(parmDict[name])
3649        if name in sigDict:
3650            sigstr += '%12.5f'%(sigDict[name])
3651        else:
3652            sigstr += 12*' '
3653    for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
3654        name = pfx+rbfx+po+rbsx
3655        namstr += '%12s'%('Ori '+po[1])
3656        valstr += '%12.5f'%(parmDict[name])
3657        if name in sigDict:
3658            sigstr += '%12.5f'%(sigDict[name])
3659        else:
3660            sigstr += 12*' '
3661    name = pfx+rbfx+'f:'+rbsx
3662    namstr += '%12s'%('Frac')
3663    valstr += '%12.5f'%(parmDict[name])
3664    if name in sigDict:
3665        sigstr += '%12.5f'%(sigDict[name])
3666    else:
3667        sigstr += 12*' '
3668    return (namstr,valstr,sigstr)
3669
3670def WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
3671    '''Cribbed version of PrintRBObjTLSAndSig but returns lists of strings.
3672    Moved so it can be used in ExportCIF
3673    '''
3674    out = []
3675    namstr = '  names :'
3676    valstr = '  values:'
3677    sigstr = '  esds  :'
3678    if 'N' not in TLS:
3679        out.append(' Thermal motion:\n')
3680    if 'T' in TLS:
3681        for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
3682            name = pfx+rbfx+pt+rbsx
3683            namstr += '%12s'%(pt[:3])
3684            valstr += '%12.5f'%(parmDict[name])
3685            if name in sigDict:
3686                sigstr += '%12.5f'%(sigDict[name])
3687            else:
3688                sigstr += 12*' '
3689        out.append(namstr+'\n')
3690        out.append(valstr+'\n')
3691        out.append(sigstr+'\n')
3692    if 'L' in TLS:
3693        namstr = '  names :'
3694        valstr = '  values:'
3695        sigstr = '  esds  :'
3696        for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
3697            name = pfx+rbfx+pt+rbsx
3698            namstr += '%12s'%(pt[:3])
3699            valstr += '%12.3f'%(parmDict[name])
3700            if name in sigDict:
3701                sigstr += '%12.3f'%(sigDict[name])
3702            else:
3703                sigstr += 12*' '
3704        out.append(namstr+'\n')
3705        out.append(valstr+'\n')
3706        out.append(sigstr+'\n')
3707    if 'S' in TLS:
3708        namstr = '  names :'
3709        valstr = '  values:'
3710        sigstr = '  esds  :'
3711        for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
3712            name = pfx+rbfx+pt+rbsx
3713            namstr += '%12s'%(pt[:3])
3714            valstr += '%12.4f'%(parmDict[name])
3715            if name in sigDict:
3716                sigstr += '%12.4f'%(sigDict[name])
3717            else:
3718                sigstr += 12*' '
3719        out.append(namstr+'\n')
3720        out.append(valstr+'\n')
3721        out.append(sigstr+'\n')
3722    if 'U' in TLS:
3723        name = pfx+rbfx+'U:'+rbsx
3724        namstr = '  names :'+'%12s'%('Uiso')
3725        valstr = '  values:'+'%12.5f'%(parmDict[name])
3726        if name in sigDict:
3727            sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
3728        else:
3729            sigstr = '  esds  :'+12*' '
3730        out.append(namstr+'\n')
3731        out.append(valstr+'\n')
3732        out.append(sigstr+'\n')
3733    return out
3734
3735def WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
3736    '''Cribbed version of PrintRBObjTorAndSig but returns lists of strings.
3737    Moved so it can be used in ExportCIF
3738    '''
3739    out = []
3740    namstr = '  names :'
3741    valstr = '  values:'
3742    sigstr = '  esds  :'
3743    out.append(' Torsions:\n')
3744    for it in range(nTors):
3745        name = pfx+'RBRTr;'+str(it)+':'+rbsx
3746        namstr += '%12s'%('Tor'+str(it))
3747        valstr += '%12.4f'%(parmDict[name])
3748        if name in sigDict:
3749            sigstr += '%12.4f'%(sigDict[name])
3750    out.append(namstr+'\n')
3751    out.append(valstr+'\n')
3752    out.append(sigstr+'\n')
3753    return out
3754
3755def WriteResRBModel(RBModel):
3756    '''Write description of a residue rigid body. Code shifted from
3757    PrintResRBModel to make usable from G2export_CIF
3758    '''
3759    out = []
3760    atNames = RBModel['atNames']
3761    rbRef = RBModel['rbRef']
3762    rbSeq = RBModel['rbSeq']
3763    out.append('    At name       x          y          z\n')
3764    for name,xyz in zip(atNames,RBModel['rbXYZ']):
3765        out.append(%8s %10.4f %10.4f %10.4f\n'%(name,xyz[0],xyz[1],xyz[2]))
3766    out.append('  Orientation defined by: %s -> %s & %s -> %s\n'%
3767        (atNames[rbRef[0]],atNames[rbRef[1]],atNames[rbRef[0]],atNames[rbRef[2]]))
3768    if rbSeq:
3769        for i,rbseq in enumerate(rbSeq):
3770#                nameLst = [atNames[i] for i in rbseq[3]]
3771            out.append('  Torsion sequence %d Bond: %s - %s riding: %s\n'%
3772                    (i,atNames[rbseq[0]],atNames[rbseq[1]],str([atNames[i] for i in rbseq[3]])))
3773    return out
3774
3775def WriteVecRBModel(RBModel,sigDict={},irb=None):
3776    '''Write description of a vector rigid body. Code shifted from
3777    PrintVecRBModel to make usable from G2export_CIF
3778    '''
3779    out = []
3780    rbRef = RBModel['rbRef']
3781    atTypes = RBModel['rbTypes']
3782    for i in range(len(RBModel['VectMag'])):
3783        if sigDict and irb is not None:
3784            name = '::RBV;'+str(i)+':'+str(irb)
3785            s = G2mth.ValEsd(RBModel['VectMag'][i],sigDict.get(name,-.0001))
3786            out.append('Vector no.: %d Magnitude: %s\n'%(i,s))
3787        else:
3788            out.append('Vector no.: %d Magnitude: %8.4f Refine? %s\n'%
3789                           (i,RBModel['VectMag'][i],RBModel['VectRef'][i]))
3790        out.append('  No. Type     vx         vy         vz\n')
3791        for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
3792            out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(j,name,xyz[0],xyz[1],xyz[2]))
3793    out.append('  No. Type      x          y          z\n')
3794    for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
3795        out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(i,name,xyz[0],xyz[1],xyz[2]))
3796    return out
Note: See TracBrowser for help on using the repository browser.