source: trunk/GSASIIstrIO.py @ 4839

Last change on this file since 4839 was 4839, checked in by vondreele, 8 months ago

change formatting for scale/phase fraction in G2ddataGUI & G2strIO now 'g' format
fix array problem in G2math/findOffset
missed a DrawAtomsReplaceByID reference to G2mth
add progress bar to Map eaks move peaks & remove the BusyCursor? in OnPeaksUnique?

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 172.6 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-03-04 18:42:09 +0000 (Thu, 04 Mar 2021) $
4# $Author: vondreele $
5# $Revision: 4839 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 4839 2021-03-04 18:42:09Z vondreele $
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: 4839 $")
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    rbVary = []
1079    rbDict = {}
1080    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
1081    if len(rbIds['Vector']):
1082        for irb,item in enumerate(rbIds['Vector']):
1083            if rigidbodyDict['Vector'][item]['useCount']:
1084                RBmags = rigidbodyDict['Vector'][item]['VectMag']
1085                RBrefs = rigidbodyDict['Vector'][item]['VectRef']
1086                for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
1087                    pid = '::RBV;'+str(i)+':'+str(irb)
1088                    rbDict[pid] = mag
1089                    if ref:
1090                        rbVary.append(pid)
1091                if Print:
1092                    pFile.write('\nVector rigid body model:\n')
1093                    PrintVecRBModel(rigidbodyDict['Vector'][item])
1094    if len(rbIds['Residue']):
1095        for item in rbIds['Residue']:
1096            if rigidbodyDict['Residue'][item]['useCount']:
1097                if Print:
1098                    pFile.write('\nResidue rigid body model:\n')
1099                    PrintResRBModel(rigidbodyDict['Residue'][item])
1100    return rbVary,rbDict
1101   
1102def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
1103    'needs a doc string'
1104   
1105    def PrintRBVectandSig(VectRB,VectSig):
1106        pFile.write('\n Rigid body vector magnitudes for %s:\n'%VectRB['RBname'])
1107        namstr = '  names :'
1108        valstr = '  values:'
1109        sigstr = '  esds  :'
1110        for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
1111            namstr += '%12s'%('Vect '+str(i))
1112            valstr += '%12.4f'%(val)
1113            if sig:
1114                sigstr += '%12.4f'%(sig)
1115            else:
1116                sigstr += 12*' '
1117        pFile.write(namstr+'\n')
1118        pFile.write(valstr+'\n')
1119        pFile.write(sigstr+'\n')       
1120       
1121    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
1122    if not RBIds['Vector']:
1123        return
1124    for irb,item in enumerate(RBIds['Vector']):
1125        if rigidbodyDict['Vector'][item]['useCount']:
1126            VectSig = []
1127            RBmags = rigidbodyDict['Vector'][item]['VectMag']
1128            for i,mag in enumerate(RBmags):
1129                name = '::RBV;'+str(i)+':'+str(irb)
1130                if name in sigDict:
1131                    VectSig.append(sigDict[name])
1132            PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)   
1133       
1134################################################################################
1135##### Phase data
1136################################################################################                   
1137def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None,seqRef=False):
1138    '''Setup the phase information for a structural refinement, used for
1139    regular and sequential refinements, optionally printing information
1140    to the .lst file (if Print is True)
1141    '''
1142           
1143    def PrintFFtable(FFtable):
1144        pFile.write('\n X-ray scattering factors:\n')
1145        pFile.write('   Symbol     fa                                      fb                                      fc\n')
1146        pFile.write(99*'-'+'\n')
1147        for Ename in FFtable:
1148            ffdata = FFtable[Ename]
1149            fa = ffdata['fa']
1150            fb = ffdata['fb']
1151            pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1152                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc']))
1153               
1154    def PrintMFtable(MFtable):
1155        pFile.write('\n <j0> Magnetic scattering factors:\n')
1156        pFile.write('   Symbol     mfa                                    mfb                                     mfc\n')
1157        pFile.write(99*'-'+'\n')
1158        for Ename in MFtable:
1159            mfdata = MFtable[Ename]
1160            fa = mfdata['mfa']
1161            fb = mfdata['mfb']
1162            pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1163                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['mfc']))
1164        pFile.write('\n <j2> Magnetic scattering factors:\n')
1165        pFile.write('   Symbol     nfa                                    nfb                                     nfc\n')
1166        pFile.write(99*'-'+'\n')
1167        for Ename in MFtable:
1168            mfdata = MFtable[Ename]
1169            fa = mfdata['nfa']
1170            fb = mfdata['nfb']
1171            pFile.write(' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f\n'%
1172                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],mfdata['nfc']))
1173               
1174    def PrintBLtable(BLtable):
1175        pFile.write('\n Neutron scattering factors:\n')
1176        pFile.write('   Symbol   isotope       mass       b       resonant terms\n')
1177        pFile.write(99*'-'+'\n')
1178        for Ename in BLtable:
1179            bldata = BLtable[Ename]
1180            isotope = bldata[0]
1181            mass = bldata[1]['Mass']
1182            if 'BW-LS' in bldata[1]:
1183                bres = bldata[1]['BW-LS']
1184                blen = 0
1185            else:
1186                blen = bldata[1]['SL'][0]
1187                bres = []
1188            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
1189            for item in bres:
1190                line += '%10.5g'%(item)
1191            pFile.write(line+'\n')
1192           
1193    def PrintRBObjects(resRBData,vecRBData):
1194               
1195        def PrintRBThermals():
1196            tlstr = ['11','22','33','12','13','23']
1197            sstr = ['12','13','21','23','31','32','AA','BB']
1198            TLS = RB['ThermalMotion'][1]
1199            TLSvar = RB['ThermalMotion'][2]
1200            if 'T' in RB['ThermalMotion'][0]:
1201                pFile.write('TLS data\n')
1202                text = ''
1203                for i in range(6):
1204                    text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
1205                pFile.write(text+'\n')
1206                if 'L' in RB['ThermalMotion'][0]: 
1207                    text = ''
1208                    for i in range(6,12):
1209                        text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
1210                    pFile.write(text+'\n')
1211                if 'S' in RB['ThermalMotion'][0]:
1212                    text = ''
1213                    for i in range(12,20):
1214                        text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
1215                    pFile.write(text+'\n')
1216            if 'U' in RB['ThermalMotion'][0]:
1217                pFile.write('Uiso data\n')
1218                text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])           
1219                pFile.write(text+'\n')
1220           
1221        if len(resRBData):
1222            for RB in resRBData:
1223                Oxyz = RB['Orig'][0]
1224                Qrijk = RB['Orient'][0]
1225                Angle = 2.0*acosd(Qrijk[0])
1226                pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'%
1227                    (RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1]))
1228                pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'%
1229                    (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1]))
1230                pFile.write('Atom site frac: %10.3f Refine? %s\n'%(RB['AtomFrac'][0],RB['AtomFrac'][1]))
1231                Torsions = RB['Torsions']
1232                if len(Torsions):
1233                    text = 'Torsions: '
1234                    for torsion in Torsions:
1235                        text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
1236                    pFile.write(text+'\n')
1237                PrintRBThermals()
1238        if len(vecRBData):
1239            for RB in vecRBData:
1240                Oxyz = RB['Orig'][0]
1241                Qrijk = RB['Orient'][0]
1242                Angle = 2.0*acosd(Qrijk[0])
1243                pFile.write('\nRBObject %s at %10.4f %10.4f %10.4f Refine? %s\n'%
1244                    (RB['RBname'],Oxyz[0],Oxyz[1],Oxyz[2],RB['Orig'][1]))
1245                pFile.write('Orientation angle,vector: %10.3f %10.4f %10.4f %10.4f Refine? %s\n'%
1246                    (Angle,Qrijk[1],Qrijk[2],Qrijk[3],RB['Orient'][1]))
1247                pFile.write('Atom site frac: %10.3f Refine? %s\n'%(RB['AtomFrac'][0],RB['AtomFrac'][1]))
1248                PrintRBThermals()
1249               
1250    def PrintAtoms(General,Atoms):
1251        cx,ct,cs,cia = General['AtomPtrs']
1252        pFile.write('\n Atoms:\n')
1253        line = '   name    type  refine?   x         y         z    '+ \
1254            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
1255        if General['Type'] == 'macromolecular':
1256            line = ' res no residue chain'+line
1257        pFile.write(line+'\n')
1258        if General['Type'] in ['nuclear','magnetic','faulted',]:
1259            pFile.write(135*'-'+'\n')
1260            for i,at in enumerate(Atoms):
1261                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
1262                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
1263                if at[cia] == 'I':
1264                    line += '%8.5f'%(at[cia+1])+48*' '
1265                else:
1266                    line += 8*' '
1267                    for j in range(6):
1268                        line += '%8.5f'%(at[cia+2+j])
1269                pFile.write(line+'\n')
1270        elif General['Type'] == 'macromolecular':
1271            pFile.write(135*'-'+'\n')           
1272            for i,at in enumerate(Atoms):
1273                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])+ \
1274                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
1275                if at[cia] == 'I':
1276                    line += '%8.4f'%(at[cia+1])+48*' '
1277                else:
1278                    line += 8*' '
1279                    for j in range(6):
1280                        line += '%8.4f'%(at[cia+2+j])
1281                pFile.write(line+'\n')
1282               
1283    def PrintMoments(General,Atoms):
1284        cx,ct,cs,cia = General['AtomPtrs']
1285        cmx = cx+4
1286        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1287        pFile.write('\n Magnetic moments:\n')
1288        line = '   name    type  refine?  Mx        My        Mz    '
1289        pFile.write(line+'\n')
1290        pFile.write(135*'-'+'\n')
1291        for i,at in enumerate(Atoms):
1292            if AtInfo[at[ct]]:
1293                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cmx])+'%10.5f'%(at[cmx+1])+ \
1294                    '%10.5f'%(at[cmx+2])
1295                pFile.write(line+'\n')
1296       
1297               
1298    def PrintWaves(General,Atoms):
1299        cx,ct,cs,cia = General['AtomPtrs']
1300        pFile.write('\n Modulation waves\n')
1301        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1302            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1303            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1304        pFile.write(135*'-'+'\n')
1305        for i,at in enumerate(Atoms):
1306            AtomSS = at[-1]['SS1']
1307            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1308                Waves = AtomSS[Stype]
1309                if len(Waves):
1310                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
1311                        (at[ct-1],at[cs],Stype,Waves[0]))
1312                else:
1313                    continue
1314                for iw,wave in enumerate(Waves[1:]):                   
1315                    line = ''
1316                    if Waves[0] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
1317                        for item in names[Stype][6:]:
1318                            line += '%8s '%(item)                       
1319                    else:
1320                        if Stype == 'Spos':
1321                            for item in names[Stype][:6]:
1322                                line += '%8s '%(item)
1323                        else:
1324                            for item in names[Stype]:
1325                                line += '%8s '%(item)
1326                    pFile.write(line+'\n')
1327                    line = ''
1328                    for item in wave[0]:
1329                        line += '%8.4f '%(item)
1330                    line += ' Refine? '+str(wave[1])
1331                    pFile.write(line+'\n')
1332       
1333    def PrintTexture(textureData):
1334        topstr = '\n Spherical harmonics texture: Order:' + \
1335            str(textureData['Order'])
1336        if textureData['Order']:
1337            pFile.write('%s Refine? %s\n'%(topstr,textureData['SH Coeff'][0]))
1338        else:
1339            pFile.write(topstr+'\n')
1340            return
1341        names = ['omega','chi','phi']
1342        line = '\n'
1343        for name in names:
1344            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1345        pFile.write(line+'\n')
1346        pFile.write('\n Texture coefficients:\n')
1347        SHcoeff = textureData['SH Coeff'][1]
1348        SHkeys = list(SHcoeff.keys())
1349        nCoeff = len(SHcoeff)
1350        nBlock = nCoeff//10+1
1351        iBeg = 0
1352        iFin = min(iBeg+10,nCoeff)
1353        for block in range(nBlock):
1354            ptlbls = ' names :'
1355            ptstr =  ' values:'
1356            for item in SHkeys[iBeg:iFin]:
1357                ptlbls += '%12s'%(item)
1358                ptstr += '%12.4f'%(SHcoeff[item]) 
1359            pFile.write(ptlbls+'\n')
1360            pFile.write(ptstr+'\n')
1361            iBeg += 10
1362            iFin = min(iBeg+10,nCoeff)
1363       
1364    def MakeRBParms(rbKey,phaseVary,phaseDict):
1365        #### patch 2/24/21 BHT: new param, AtomFrac in RB
1366        if 'AtomFrac' not in RB: raise Exception('out of date RB: edit in RB Models')
1367        # end patch
1368        rbid = str(rbids.index(RB['RBId']))
1369        pfxRB = pfx+'RB'+rbKey+'P'
1370        pstr = ['x','y','z']
1371        ostr = ['a','i','j','k']
1372        for i in range(3):
1373            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1374            phaseDict[name] = RB['Orig'][0][i]
1375            if RB['Orig'][1]:
1376                phaseVary += [name,]
1377        pfxRB = pfx+'RB'+rbKey+'O'
1378        A,V = G2mth.Q2AV(RB['Orig'][0])
1379        fixAxis = [0, np.abs(V).argmax()+1]
1380        for i in range(4):
1381            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1382            phaseDict[name] = RB['Orient'][0][i]
1383            if RB['Orient'][1] == 'AV' and i:
1384                phaseVary += [name,]
1385            elif RB['Orient'][1] == 'A' and not i:
1386                phaseVary += [name,]
1387            elif RB['Orient'][1] == 'V' and i not in fixAxis:
1388                phaseVary += [name,]
1389        name = pfx+'RB'+rbKey+'f:'+str(iRB)+':'+rbid
1390        phaseDict[name] = RB['AtomFrac'][0]
1391        if RB['AtomFrac'][1]:
1392            phaseVary += [name,]
1393               
1394    def MakeRBThermals(rbKey,phaseVary,phaseDict):
1395        rbid = str(rbids.index(RB['RBId']))
1396        tlstr = ['11','22','33','12','13','23']
1397        sstr = ['12','13','21','23','31','32','AA','BB']
1398        if 'T' in RB['ThermalMotion'][0]:
1399            pfxRB = pfx+'RB'+rbKey+'T'
1400            for i in range(6):
1401                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1402                phaseDict[name] = RB['ThermalMotion'][1][i]
1403                if RB['ThermalMotion'][2][i]:
1404                    phaseVary += [name,]
1405        if 'L' in RB['ThermalMotion'][0]:
1406            pfxRB = pfx+'RB'+rbKey+'L'
1407            for i in range(6):
1408                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1409                phaseDict[name] = RB['ThermalMotion'][1][i+6]
1410                if RB['ThermalMotion'][2][i+6]:
1411                    phaseVary += [name,]
1412        if 'S' in RB['ThermalMotion'][0]:
1413            pfxRB = pfx+'RB'+rbKey+'S'
1414            for i in range(8):
1415                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1416                phaseDict[name] = RB['ThermalMotion'][1][i+12]
1417                if RB['ThermalMotion'][2][i+12]:
1418                    phaseVary += [name,]
1419        if 'U' in RB['ThermalMotion'][0]:
1420            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1421            phaseDict[name] = RB['ThermalMotion'][1][0]
1422            if RB['ThermalMotion'][2][0]:
1423                phaseVary += [name,]
1424               
1425    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1426        rbid = str(rbids.index(RB['RBId']))
1427        pfxRB = pfx+'RB'+rbKey+'Tr;'
1428        for i,torsion in enumerate(RB['Torsions']):
1429            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1430            phaseDict[name] = torsion[0]
1431            if torsion[1]:
1432                phaseVary += [name,]
1433                   
1434    if Print:
1435        pFile.write('\n Phases:\n')
1436    phaseVary = []
1437    phaseDict = {}
1438    pawleyLookup = {}
1439    FFtables = {}                   #scattering factors - xrays
1440    MFtables = {}                   #Mag. form factors
1441    BLtables = {}                   # neutrons
1442    Natoms = {}
1443    maxSSwave = {}
1444    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1445    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1446    atomIndx = {}
1447    for name in PhaseData:
1448        General = PhaseData[name]['General']
1449        pId = PhaseData[name]['pId']
1450        pfx = str(pId)+'::'
1451        FFtable = G2el.GetFFtable(General['AtomTypes'])
1452        BLtable = G2el.GetBLtable(General)
1453        FFtables.update(FFtable)
1454        BLtables.update(BLtable)
1455        phaseDict[pfx+'isMag'] = False
1456        SGData = General['SGData']
1457        SGtext,SGtable = G2spc.SGPrint(SGData)
1458        if General['Type'] == 'magnetic':
1459            MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g'])
1460            MFtables.update(MFtable)
1461            phaseDict[pfx+'isMag'] = True
1462            SpnFlp = SGData['SpnFlp']
1463        Atoms = PhaseData[name]['Atoms']
1464        if Atoms and not General.get('doPawley'):
1465            cx,ct,cs,cia = General['AtomPtrs']
1466            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1467        PawleyRef = PhaseData[name].get('Pawley ref',[])
1468        cell = General['Cell']
1469        A = G2lat.cell2A(cell[1:7])
1470        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1471            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1472        if cell[0]:
1473            phaseVary += cellVary(pfx,SGData)       #also fills in symmetry required constraints
1474        SSGtext = []    #no superstructure
1475        im = 0
1476        if General.get('Modulated',False):
1477            im = 1      #refl offset
1478            Vec,vRef,maxH = General['SuperVec']
1479            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1480            SSGData = General['SSGData']
1481            SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1482            if vRef:
1483                phaseVary += modVary(pfx,SSGData)       
1484        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1485        if resRBData:
1486            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
1487            for iRB,RB in enumerate(resRBData):
1488                MakeRBParms('R',phaseVary,phaseDict)
1489                MakeRBThermals('R',phaseVary,phaseDict)
1490                MakeRBTorsions('R',phaseVary,phaseDict)
1491       
1492        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1493        if vecRBData:
1494            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
1495            for iRB,RB in enumerate(vecRBData):
1496                MakeRBParms('V',phaseVary,phaseDict)
1497                MakeRBThermals('V',phaseVary,phaseDict)
1498                   
1499        Natoms[pfx] = 0
1500        maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1501        if Atoms and not General.get('doPawley'):
1502            cx,ct,cs,cia = General['AtomPtrs']
1503            Natoms[pfx] = len(Atoms)
1504            for i,at in enumerate(Atoms):
1505                atomIndx[at[cia+8]] = [pfx,i]      #lookup table for restraints
1506                phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1507                    pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1508                    pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1509                    pfx+'AI/A:'+str(i):at[cia],})
1510                if at[cia] == 'I':
1511                    phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1512                else:
1513                    phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1514                        pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1515                if General['Type'] == 'magnetic':
1516                    phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]})
1517                if 'F' in at[ct+1]:
1518                    phaseVary.append(pfx+'Afrac:'+str(i))
1519                if 'X' in at[ct+1]:
1520                    try:    #patch for sytsym name changes
1521                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1522                    except KeyError:
1523                        Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1524                        at[cs] = Sytsym
1525                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1526                    names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1527                    equivs = [[],[],[]]
1528                    for j in range(3):
1529                        if xId[j] > 0:                               
1530                            phaseVary.append(names[j])
1531                            equivs[xId[j]-1].append([names[j],xCoef[j]])
1532                    for equiv in equivs:
1533                        if len(equiv) > 1:
1534                            name = equiv[0][0]
1535                            coef = equiv[0][1]
1536                            for eqv in equiv[1:]:
1537                                eqv[1] /= coef
1538                                G2mv.StoreEquivalence(name,(eqv,))
1539                if 'U' in at[ct+1]:
1540                    if at[cia] == 'I':
1541                        phaseVary.append(pfx+'AUiso:'+str(i))
1542                    else:
1543                        try:    #patch for sytsym name changes
1544                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1545                        except KeyError:
1546                            Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1547                            at[cs] = Sytsym
1548                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1549                        names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1550                            pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1551                        equivs = [[],[],[],[],[],[]]
1552                        for j in range(6):
1553                            if uId[j] > 0:                               
1554                                phaseVary.append(names[j])
1555                                equivs[uId[j]-1].append([names[j],uCoef[j]])
1556                        for equiv in equivs:
1557                            if len(equiv) > 1:
1558                                name = equiv[0][0]
1559                                coef = equiv[0][1]
1560                                for eqv in equiv[1:]:
1561                                    eqv[1] /= coef
1562                                    G2mv.StoreEquivalence(name,(eqv,))
1563                if 'M' in at[ct+1]:
1564                    SytSym,Mul,Nop,dupDir = G2spc.SytSym(at[cx:cx+3],SGData)
1565                    mId,mCoef = G2spc.GetCSpqinel(SpnFlp,dupDir)
1566                    names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1567                    equivs = [[],[],[]]
1568                    for j in range(3):
1569                        if mId[j] > 0:
1570                            phaseVary.append(names[j])
1571                            equivs[mId[j]-1].append([names[j],mCoef[j]])
1572                    for equiv in equivs:
1573                        if len(equiv) > 1:
1574                            name = equiv[0][0]
1575                            coef = equiv[0][1]
1576                            for eqv in equiv[1:]:
1577                                eqv[1] /= coef
1578                                G2mv.StoreEquivalence(name,(eqv,))
1579                if General.get('Modulated',False):
1580                    AtomSS = at[-1]['SS1']
1581                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1582                        Waves = AtomSS[Stype]
1583                        if len(Waves):
1584                            waveType = Waves[0]
1585                        else:
1586                            continue
1587                        phaseDict[pfx+Stype[1].upper()+'waveType:'+str(i)] = waveType
1588                        nx = 0
1589                        for iw,wave in enumerate(Waves[1:]):
1590                            if not iw:
1591                                if waveType in ['ZigZag','Block']:
1592                                    nx = 1
1593                                CSI = G2spc.GetSSfxuinel(waveType,Stype,1,at[cx:cx+3],SGData,SSGData)
1594                            else:
1595                                CSI = G2spc.GetSSfxuinel('Fourier',Stype,iw+1-nx,at[cx:cx+3],SGData,SSGData)
1596                            uId,uCoef = CSI[0]
1597                            stiw = str(i)+':'+str(iw)
1598                            if Stype == 'Spos':
1599                                if waveType in ['ZigZag','Block',] and not iw:
1600                                    names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1601                                    equivs = [[],[], [],[],[]]
1602                                else:
1603                                    names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1604                                        pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1605                                    equivs = [[],[],[], [],[],[]]
1606                            elif Stype == 'Sadp':
1607                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1608                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1609                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1610                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1611                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1612                            elif Stype == 'Sfrac':
1613                                equivs = [[],[]]
1614                                if 'Crenel' in waveType and not iw:
1615                                    names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1616                                else:
1617                                    names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1618                            elif Stype == 'Smag':
1619                                equivs = [[],[],[], [],[],[]]
1620                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1621                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1622                            phaseDict.update(dict(zip(names,wave[0])))
1623                            if wave[1]: #what do we do here for multiple terms in modulation constraints?
1624                                for j in range(len(equivs)):
1625                                    if uId[j][0] > 0:                               
1626                                        phaseVary.append(names[j])
1627                                        equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1628                                for equiv in equivs:
1629                                    if len(equiv) > 1:
1630                                        name = equiv[0][0]
1631                                        coef = equiv[0][1]
1632                                        for eqv in equiv[1:]:
1633                                            eqv[1] /= coef
1634                                            G2mv.StoreEquivalence(name,(eqv,))
1635                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1636            textureData = General['SH Texture']
1637            if textureData['Order'] and not seqRef:
1638                phaseDict[pfx+'SHorder'] = textureData['Order']
1639                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1640                for item in ['omega','chi','phi']:
1641                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1642                    if textureData['Sample '+item][0]:
1643                        phaseVary.append(pfx+'SH '+item)
1644                for item in textureData['SH Coeff'][1]:
1645                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1646                    if textureData['SH Coeff'][0]:
1647                        phaseVary.append(pfx+item)
1648               
1649            if Print:
1650                pFile.write('\n Phase name: %s\n'%General['Name'])
1651                pFile.write(135*'='+'\n')
1652                PrintFFtable(FFtable)
1653                PrintBLtable(BLtable)
1654                if General['Type'] == 'magnetic':
1655                    PrintMFtable(MFtable)
1656                pFile.write('\n')
1657                #how do we print magnetic symmetry table? TBD
1658                if len(SSGtext):    #if superstructure
1659                    for line in SSGtext: pFile.write(line+'\n')
1660                    if len(SSGtable):
1661                        for item in SSGtable:
1662                            line = ' %s '%(item)
1663                            pFile.write(line+'\n') 
1664                    else:
1665                        pFile.write(' ( 1)    %s\n'%(SSGtable[0]))
1666                else:
1667                    for line in SGtext: pFile.write(line+'\n')
1668                    if len(SGtable):
1669                        for item in SGtable:
1670                            line = ' %s '%(item)
1671                            pFile.write(line+'\n') 
1672                    else:
1673                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1674                PrintRBObjects(resRBData,vecRBData)
1675                PrintAtoms(General,Atoms)
1676                if General['Type'] == 'magnetic':
1677                    PrintMoments(General,Atoms)
1678                if General.get('Modulated',False):
1679                    PrintWaves(General,Atoms)
1680                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1681                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1682                if len(SSGtext):    #if superstructure
1683                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1684                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1685                if not seqRef:
1686                    PrintTexture(textureData)
1687                if name in RestraintDict:
1688                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1689                        textureData,RestraintDict[name],pFile)
1690                   
1691        elif PawleyRef:
1692            if Print:
1693                pFile.write('\n Phase name: %s\n'%General['Name'])
1694                pFile.write(135*'='+'\n')
1695                pFile.write('\n')
1696                if len(SSGtext):    #if superstructure
1697                    for line in SSGtext: pFile.write(line+'\n')
1698                    if len(SSGtable):
1699                        for item in SSGtable:
1700                            line = ' %s '%(item)
1701                            pFile.write(line+'\n') 
1702                    else:
1703                        pFile.write(' ( 1)    %s\n'%SSGtable[0])
1704                else:
1705                    for line in SGtext: pFile.write(line+'\n')
1706                    if len(SGtable):
1707                        for item in SGtable:
1708                            line = ' %s '%(item)
1709                            pFile.write(line+'\n')
1710                    else:
1711                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1712                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1713                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1714                if len(SSGtext):    #if superstructure
1715                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1716                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1717            pawleyVary = []
1718            for i,refl in enumerate(PawleyRef):
1719                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1720                if im:
1721                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1722                else:
1723                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1724                if refl[5+im]:
1725                    pawleyVary.append(pfx+'PWLref:'+str(i))
1726            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1727            phaseVary += pawleyVary
1728               
1729    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1730   
1731def cellFill(pfx,SGData,parmDict,sigDict): 
1732    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1733    from the parameter and sig dictionaries.
1734
1735    :param str pfx: parameter prefix ("n::", where n is a phase number)
1736    :param dict SGdata: a symmetry object
1737    :param dict parmDict: a dictionary of parameters
1738    :param dict sigDict:  a dictionary of uncertainties on parameters
1739
1740    :returns: A,sigA where each is a list of six terms with the A terms
1741    '''
1742    if SGData['SGLaue'] in ['-1',]:
1743        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1744            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1745    elif SGData['SGLaue'] in ['2/m',]:
1746        if SGData['SGUniq'] == 'a':
1747            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1748                0,0,parmDict[pfx+'A5']]
1749        elif SGData['SGUniq'] == 'b':
1750            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1751                0,parmDict[pfx+'A4'],0]
1752        else:
1753            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1754                parmDict[pfx+'A3'],0,0]
1755    elif SGData['SGLaue'] in ['mmm',]:
1756        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1757    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1758        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1759    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1760        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1761            parmDict[pfx+'A0'],0,0]
1762    elif SGData['SGLaue'] in ['3R', '3mR']:
1763        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1764            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1765    elif SGData['SGLaue'] in ['m3m','m3']:
1766        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1767
1768    try:
1769        if SGData['SGLaue'] in ['-1',]:
1770            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1771                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1772        elif SGData['SGLaue'] in ['2/m',]:
1773            if SGData['SGUniq'] == 'a':
1774                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1775                    0,0,sigDict[pfx+'A5']]
1776            elif SGData['SGUniq'] == 'b':
1777                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1778                    0,sigDict[pfx+'A4'],0]
1779            else:
1780                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1781                    sigDict[pfx+'A3'],0,0]
1782        elif SGData['SGLaue'] in ['mmm',]:
1783            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1784        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1785            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1786        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1787            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1788        elif SGData['SGLaue'] in ['3R', '3mR']:
1789            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1790        elif SGData['SGLaue'] in ['m3m','m3']:
1791            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1792    except KeyError:
1793        sigA = [0,0,0,0,0,0]
1794    return A,sigA
1795       
1796def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1797    'needs a doc string'
1798    if phaseRest:
1799        Amat = G2lat.cell2AB(cell)[0]
1800        cx,ct,cs = AtPtrs[:3]
1801        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1802            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1803            ['ChemComp','Sites'],['Texture','HKLs']]
1804        for name,rest in names:
1805            itemRest = phaseRest[name]
1806            if rest in itemRest and itemRest[rest] and itemRest['Use']:
1807                pFile.write('\n  %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use'])))
1808                if name in ['Bond','Angle','Plane','Chiral']:
1809                    pFile.write('     calc       obs      sig   delt/sig  atoms(symOp): \n')
1810                    for indx,ops,obs,esd in itemRest[rest]:
1811                        try:
1812                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1813                            AtName = ''
1814                            for i,Aname in enumerate(AtNames):
1815                                AtName += Aname
1816                                if ops[i] == '1':
1817                                    AtName += '-'
1818                                else:
1819                                    AtName += '+('+ops[i]+')-'
1820                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1821                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1822                            if name == 'Bond':
1823                                calc = G2mth.getRestDist(XYZ,Amat)
1824                            elif name == 'Angle':
1825                                calc = G2mth.getRestAngle(XYZ,Amat)
1826                            elif name == 'Plane':
1827                                calc = G2mth.getRestPlane(XYZ,Amat)
1828                            elif name == 'Chiral':
1829                                calc = G2mth.getRestChiral(XYZ,Amat)
1830                            pFile.write(' %9.3f %9.3f %8.3f %8.3f   %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1]))
1831                        except KeyError:
1832                            del itemRest[rest]
1833                elif name in ['Torsion','Rama']:
1834                    pFile.write('  atoms(symOp)  calc  obs  sig  delt/sig  torsions: \n')
1835                    coeffDict = itemRest['Coeff']
1836                    for indx,ops,cofName,esd in itemRest[rest]:
1837                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1838                        AtName = ''
1839                        for i,Aname in enumerate(AtNames):
1840                            AtName += Aname+'+('+ops[i]+')-'
1841                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1842                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1843                        if name == 'Torsion':
1844                            tor = G2mth.getRestTorsion(XYZ,Amat)
1845                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1846                            pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1]))
1847                        else:
1848                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1849                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1850                            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]))
1851                elif name == 'ChemComp':
1852                    pFile.write('     atoms   mul*frac  factor     prod\n')
1853                    for indx,factors,obs,esd in itemRest[rest]:
1854                        try:
1855                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1856                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1857                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1858                            mulfrac = mul*frac
1859                            calcs = mul*frac*factors
1860                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1861                                pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc))
1862                            pFile.write(' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd))
1863                        except KeyError:
1864                            del itemRest[rest]
1865                elif name == 'Texture' and textureData['Order']:
1866                    Start = False
1867                    SHCoef = textureData['SH Coeff'][1]
1868                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1869                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1870                    pFile.write ('    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd \n')
1871                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1872                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1873                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1874                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1875                        Z = ma.masked_greater(Z,0.0)
1876                        num = ma.count(Z)
1877                        sum = 0
1878                        if num:
1879                            sum = np.sum(Z)
1880                        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))
1881       
1882def getCellEsd(pfx,SGData,A,covData):
1883    'needs a doc string'
1884    rVsq = G2lat.calc_rVsq(A)
1885    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1886    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1887    varyList = covData['varyList']
1888    covMatrix = covData['covMatrix']
1889    if len(covMatrix):
1890        vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1891        if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']:
1892            vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0]
1893            vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0]
1894            vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2]
1895        elif SGData['SGLaue'] in ['m3','m3m']:
1896            vcov[0:3,0:3] = vcov[0,0]
1897        elif SGData['SGLaue'] in ['4/m', '4/mmm']:
1898            vcov[0:2,0:2] = vcov[0,0]
1899            vcov[1,2] = vcov[2,1] = vcov[0,2]
1900        elif SGData['SGLaue'] in ['3R','3mR']:
1901            vcov[0:3,0:3] = vcov[0,0]
1902    #        vcov[4,4] = vcov[5,5] = vcov[3,3]
1903            vcov[3:6,3:6] = vcov[3,3]
1904            vcov[0:3,3:6] = vcov[0,3]
1905            vcov[3:6,0:3] = vcov[3,0]
1906    else:
1907        vcov = np.eye(6)
1908    delt = 1.e-9
1909    drVdA = np.zeros(6)
1910    for i in range(6):
1911        A[i] += delt
1912        drVdA[i] = G2lat.calc_rVsq(A)
1913        A[i] -= 2*delt
1914        drVdA[i] -= G2lat.calc_rVsq(A)
1915        A[i] += delt
1916    drVdA /= 2.*delt   
1917    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1918    Vol = 1/np.sqrt(rVsq)
1919    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1920   
1921    dcdA = np.zeros((6,6))
1922    for i in range(6):
1923        pdcdA =np.zeros(6)
1924        A[i] += delt
1925        pdcdA += G2lat.A2cell(A)
1926        A[i] -= 2*delt
1927        pdcdA -= G2lat.A2cell(A)
1928        A[i] += delt
1929        dcdA[i] = pdcdA/(2.*delt)
1930   
1931    sigMat = np.inner(dcdA,np.inner(dcdA,vcov))
1932    var = np.diag(sigMat)
1933    CS = np.where(var>0.,np.sqrt(var),0.)
1934    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']:
1935        CS[3:6] = 0.0
1936    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1937   
1938def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1939    '''Called after a refinement to transfer parameters from the parameter dict to
1940    the phase(s) information read from a GPX file. Also prints values to the .lst file
1941    '''
1942   
1943    def PrintAtomsAndSig(General,Atoms,atomsSig):
1944        pFile.write('\n Atoms:\n')
1945        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1946        if General['Type'] == 'macromolecular':
1947            line = ' res no residue chain '+line
1948        cx,ct,cs,cia = General['AtomPtrs']
1949        pFile.write(line+'\n')
1950        pFile.write(135*'-'+'\n')
1951        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1952            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1953        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1954        for atyp in General['AtomTypes']:       #zero composition
1955            General['NoAtoms'][atyp] = 0.
1956        for i,at in enumerate(Atoms):
1957            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1958            if General['Type'] == 'macromolecular':
1959                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1960                valstr = ' values:          '
1961                sigstr = ' sig   :          '
1962            else:
1963                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1964                valstr = ' values:'
1965                sigstr = ' sig   :'
1966            for ind in range(cx,cx+4):
1967                sigind = str(i)+':'+str(ind)
1968                valstr += fmt[ind]%(at[ind])                   
1969                if sigind in atomsSig:
1970                    sigstr += fmt[ind]%(atomsSig[sigind])
1971                else:
1972                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1973            if at[cia] == 'I':
1974                valstr += fmt[cia+1]%(at[cia+1])
1975                if '%d:%d'%(i,cia+1) in atomsSig:
1976                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1977                else:
1978                    sigstr += 8*' '
1979            else:
1980                valstr += 8*' '
1981                sigstr += 8*' '
1982                for ind in range(cia+2,cia+8):
1983                    sigind = str(i)+':'+str(ind)
1984                    valstr += fmt[ind]%(at[ind])
1985                    if sigind in atomsSig:                       
1986                        sigstr += fmt[ind]%(atomsSig[sigind])
1987                    else:
1988                        sigstr += 8*' '
1989            pFile.write(name+'\n')
1990            pFile.write(valstr+'\n')
1991            pFile.write(sigstr+'\n')
1992           
1993    def PrintMomentsAndSig(General,Atoms,atomsSig):
1994        cell = General['Cell'][1:7]
1995        G = G2lat.fillgmat(cell)
1996        ast = np.sqrt(np.diag(G))
1997        GS = G/np.outer(ast,ast)
1998        pFile.write('\n Magnetic Moments:\n')    #add magnitude & angle, etc.? TBD
1999        line = '   name       Mx        My        Mz       |Mag|'
2000        cx,ct,cs,cia = General['AtomPtrs']
2001        cmx = cx+4
2002        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
2003        pFile.write(line+'\n')
2004        pFile.write(135*'-'+'\n')
2005        fmt = {0:'%7s',ct:'%7s',cmx:'%10.3f',cmx+1:'%10.3f',cmx+2:'%10.3f'}
2006        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
2007        for i,at in enumerate(Atoms):
2008            if AtInfo[at[ct]]:
2009                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
2010                valstr = ' values:'
2011                sigstr = ' sig   :'
2012                for ind in range(cmx,cmx+3):
2013                    sigind = str(i)+':'+str(ind)
2014                    valstr += fmt[ind]%(at[ind])                   
2015                    if sigind in atomsSig:
2016                        sigstr += fmt[ind]%(atomsSig[sigind])
2017                    else:
2018                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
2019                mag = np.array(at[cmx:cmx+3])
2020                Mag = np.sqrt(np.inner(mag,np.inner(mag,GS)))
2021                valstr += '%10.3f'%Mag
2022                sigstr += 10*' '
2023                pFile.write(name+'\n')
2024                pFile.write(valstr+'\n')
2025                pFile.write(sigstr+'\n')
2026           
2027    def PrintWavesAndSig(General,Atoms,wavesSig):
2028        cx,ct,cs,cia = General['AtomPtrs']
2029        pFile.write('\n Modulation waves\n')
2030        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
2031            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
2032            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
2033        pFile.write(135*'-'+'\n')
2034        for i,at in enumerate(Atoms):
2035            AtomSS = at[-1]['SS1']
2036            for Stype in ['Sfrac','Spos','Sadp','Smag']:
2037                Waves = AtomSS[Stype]
2038                if len(Waves) > 1:
2039                    waveType = Waves[0]
2040                else:
2041                    continue
2042                if len(Waves):
2043                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
2044                        (at[ct-1],at[cs],Stype,waveType))
2045                    for iw,wave in enumerate(Waves[1:]):
2046                        stiw = ':'+str(i)+':'+str(iw)
2047                        namstr = '  names :'
2048                        valstr = '  values:'
2049                        sigstr = '  esds  :'
2050                        if Stype == 'Spos':
2051                            nt = 6
2052                            ot = 0
2053                            if waveType in ['ZigZag','Block',] and not iw:
2054                                nt = 5
2055                                ot = 6
2056                            for j in range(nt):
2057                                name = names['Spos'][j+ot]
2058                                namstr += '%12s'%(name)
2059                                valstr += '%12.4f'%(wave[0][j])
2060                                if name+stiw in wavesSig:
2061                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2062                                else:
2063                                    sigstr += 12*' '
2064                        elif Stype == 'Sfrac':
2065                            ot = 0
2066                            if 'Crenel' in waveType and not iw:
2067                                ot = 2
2068                            for j in range(2):
2069                                name = names['Sfrac'][j+ot]
2070                                namstr += '%12s'%(names['Sfrac'][j+ot])
2071                                valstr += '%12.4f'%(wave[0][j])
2072                                if name+stiw in wavesSig:
2073                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2074                                else:
2075                                    sigstr += 12*' '
2076                        elif Stype == 'Sadp':
2077                            for j in range(12):
2078                                name = names['Sadp'][j]
2079                                namstr += '%10s'%(names['Sadp'][j])
2080                                valstr += '%10.6f'%(wave[0][j])
2081                                if name+stiw in wavesSig:
2082                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
2083                                else:
2084                                    sigstr += 10*' '
2085                        elif Stype == 'Smag':
2086                            for j in range(6):
2087                                name = names['Smag'][j]
2088                                namstr += '%12s'%(names['Smag'][j])
2089                                valstr += '%12.4f'%(wave[0][j])
2090                                if name+stiw in wavesSig:
2091                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2092                                else:
2093                                    sigstr += 12*' '
2094                               
2095                    pFile.write(namstr+'\n')
2096                    pFile.write(valstr+'\n')
2097                    pFile.write(sigstr+'\n')
2098       
2099               
2100    def PrintRBObjPOAndSig(rbfx,rbsx):
2101        for i in WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
2102            pFile.write(i+'\n')
2103       
2104    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
2105        for i in WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
2106            pFile.write(i)
2107       
2108    def PrintRBObjTorAndSig(rbsx):
2109        nTors = len(RBObj['Torsions'])
2110        if nTors:
2111            for i in WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
2112                pFile.write(i)
2113               
2114    def PrintSHtextureAndSig(textureData,SHtextureSig):
2115        Tindx = 1.0
2116        Tvar = 0.0
2117        pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
2118        names = ['omega','chi','phi']
2119        namstr = '  names :'
2120        ptstr =  '  values:'
2121        sigstr = '  esds  :'
2122        for name in names:
2123            namstr += '%12s'%(name)
2124            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2125            if 'Sample '+name in SHtextureSig:
2126                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2127            else:
2128                sigstr += 12*' '
2129        pFile.write(namstr+'\n')
2130        pFile.write(ptstr+'\n')
2131        pFile.write(sigstr+'\n')
2132        pFile.write('\n Texture coefficients:\n')
2133        SHcoeff = textureData['SH Coeff'][1]
2134        SHkeys = list(SHcoeff.keys())
2135        nCoeff = len(SHcoeff)
2136        nBlock = nCoeff//10+1
2137        iBeg = 0
2138        iFin = min(iBeg+10,nCoeff)
2139        for block in range(nBlock):
2140            namstr = '  names :'
2141            ptstr =  '  values:'
2142            sigstr = '  esds  :'
2143            for name in SHkeys[iBeg:iFin]:
2144                namstr += '%12s'%(name)
2145                ptstr += '%12.3f'%(SHcoeff[name])
2146                l = 2.0*eval(name.strip('C'))[0]+1
2147                Tindx += SHcoeff[name]**2/l
2148                if name in SHtextureSig:
2149                    Tvar += (2.*SHcoeff[name]*SHtextureSig[name]/l)**2
2150                    sigstr += '%12.3f'%(SHtextureSig[name])
2151                else:
2152                    sigstr += 12*' '
2153            pFile.write(namstr+'\n')
2154            pFile.write(ptstr+'\n')
2155            pFile.write(sigstr+'\n')
2156            iBeg += 10
2157            iFin = min(iBeg+10,nCoeff)
2158        pFile.write(' Texture index J = %.3f(%d)'%(Tindx,int(1000*np.sqrt(Tvar))))
2159           
2160    ##########################################################################
2161    # SetPhaseData starts here
2162    pFile.write('\n Phases:\n')
2163    for phase in Phases:
2164        pFile.write(' Result for phase: %s\n'%phase)
2165        pFile.write(135*'='+'\n')
2166        Phase = Phases[phase]
2167        General = Phase['General']
2168        SGData = General['SGData']
2169        Atoms = Phase['Atoms']
2170        AtLookup = []
2171        if Atoms and not General.get('doPawley'):
2172            cx,ct,cs,cia = General['AtomPtrs']
2173            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
2174        cell = General['Cell']
2175        pId = Phase['pId']
2176        pfx = str(pId)+'::'
2177        if cell[0]:
2178            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
2179            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
2180            pFile.write(' Reciprocal metric tensor: \n')
2181            ptfmt = "%15.9f"
2182            names = ['A11','A22','A33','A12','A13','A23']
2183            namstr = '  names :'
2184            ptstr =  '  values:'
2185            sigstr = '  esds  :'
2186            for name,a,siga in zip(names,A,sigA):
2187                namstr += '%15s'%(name)
2188                ptstr += ptfmt%(a)
2189                if siga:
2190                    sigstr += ptfmt%(siga)
2191                else:
2192                    sigstr += 15*' '
2193            pFile.write(namstr+'\n')
2194            pFile.write(ptstr+'\n')
2195            pFile.write(sigstr+'\n')
2196            cell[1:7] = G2lat.A2cell(A)
2197            cell[7] = G2lat.calc_V(A)
2198            pFile.write(' New unit cell:\n')
2199            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
2200            names = ['a','b','c','alpha','beta','gamma','Volume']
2201            namstr = '  names :'
2202            ptstr =  '  values:'
2203            sigstr = '  esds  :'
2204            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
2205                namstr += '%12s'%(name)
2206                ptstr += fmt%(a)
2207                if siga:
2208                    sigstr += fmt%(siga)
2209                else:
2210                    sigstr += 12*' '
2211            pFile.write(namstr+'\n')
2212            pFile.write(ptstr+'\n')
2213            pFile.write(sigstr+'\n')
2214        ik = 6  #for Pawley stuff below
2215        if General.get('Modulated',False):
2216            ik = 7
2217            Vec,vRef,maxH = General['SuperVec']
2218            if vRef:
2219                pFile.write(' New modulation vector:\n')
2220                namstr = '  names :'
2221                ptstr =  '  values:'
2222                sigstr = '  esds  :'
2223                for iv,var in enumerate(['mV0','mV1','mV2']):
2224                    namstr += '%12s'%(pfx+var)
2225                    ptstr += '%12.6f'%(parmDict[pfx+var])
2226                    if pfx+var in sigDict:
2227                        Vec[iv] = parmDict[pfx+var]
2228                        sigstr += '%12.6f'%(sigDict[pfx+var])
2229                    else:
2230                        sigstr += 12*' '
2231                pFile.write(namstr+'\n')
2232                pFile.write(ptstr+'\n')
2233                pFile.write(sigstr+'\n')
2234           
2235        General['Mass'] = 0.
2236        if Phase['General'].get('doPawley'):
2237            pawleyRef = Phase['Pawley ref']
2238            for i,refl in enumerate(pawleyRef):
2239                key = pfx+'PWLref:'+str(i)
2240                refl[ik] = parmDict[key]
2241                if key in sigDict:
2242                    refl[ik+1] = sigDict[key]
2243                else:
2244                    refl[ik+1] = 0
2245        else:
2246            VRBIds = RBIds['Vector']
2247            RRBIds = RBIds['Residue']
2248            RBModels = Phase['RBModels']
2249            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2250                jrb = VRBIds.index(RBObj['RBId'])
2251                rbsx = str(irb)+':'+str(jrb)
2252                pFile.write(' Vector rigid body parameters:\n')
2253                PrintRBObjPOAndSig('RBV',rbsx)
2254                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2255            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2256                jrb = RRBIds.index(RBObj['RBId'])
2257                rbsx = str(irb)+':'+str(jrb)
2258                pFile.write(' Residue rigid body parameters:\n')
2259                PrintRBObjPOAndSig('RBR',rbsx)
2260                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2261                PrintRBObjTorAndSig(rbsx)
2262            atomsSig = {}
2263            wavesSig = {}
2264            cx,ct,cs,cia = General['AtomPtrs']
2265            for i,at in enumerate(Atoms):
2266                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2267                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2268                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2269                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2270                for ind in range(cx,cx+4):
2271                    at[ind] = parmDict[names[ind]]
2272                    if ind in range(cx,cx+3):
2273                        name = names[ind].replace('A','dA')
2274                    else:
2275                        name = names[ind]
2276                    if name in sigDict:
2277                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2278                if at[cia] == 'I':
2279                    at[cia+1] = parmDict[names[cia+1]]
2280                    if names[cia+1] in sigDict:
2281                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2282                else:
2283                    for ind in range(cia+2,cia+8):
2284                        at[ind] = parmDict[names[ind]]
2285                        if names[ind] in sigDict:
2286                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2287                if General['Type'] == 'magnetic':
2288                    for ind in range(cx+4,cx+7):
2289                        at[ind] = parmDict[names[ind]]
2290                        if names[ind] in sigDict:
2291                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2292                ind = General['AtomTypes'].index(at[ct])
2293                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2294                if General.get('Modulated',False):
2295                    AtomSS = at[-1]['SS1']
2296                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2297                        Waves = AtomSS[Stype]
2298                        if len(Waves):
2299                            waveType = Waves[0]
2300                        else:
2301                            continue
2302                        for iw,wave in enumerate(Waves[1:]):
2303                            stiw = str(i)+':'+str(iw)
2304                            if Stype == 'Spos':
2305                                if waveType in ['ZigZag','Block',] and not iw:
2306                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2307                                else:
2308                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2309                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2310                            elif Stype == 'Sadp':
2311                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2312                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2313                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2314                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2315                            elif Stype == 'Sfrac':
2316                                if 'Crenel' in waveType and not iw:
2317                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2318                                else:
2319                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2320                            elif Stype == 'Smag':
2321                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2322                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2323                            for iname,name in enumerate(names):
2324                                AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2325                                if pfx+name in sigDict:
2326                                    wavesSig[name] = sigDict[pfx+name]
2327                   
2328            PrintAtomsAndSig(General,Atoms,atomsSig)
2329            if General['Type'] == 'magnetic':
2330                PrintMomentsAndSig(General,Atoms,atomsSig)
2331            if General.get('Modulated',False):
2332                PrintWavesAndSig(General,Atoms,wavesSig)
2333               
2334            density = G2mth.getDensity(General)[0]
2335            pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density))
2336           
2337       
2338        textureData = General['SH Texture']   
2339        if textureData['Order']:
2340            SHtextureSig = {}
2341            for name in ['omega','chi','phi']:
2342                aname = pfx+'SH '+name
2343                textureData['Sample '+name][1] = parmDict[aname]
2344                if aname in sigDict:
2345                    SHtextureSig['Sample '+name] = sigDict[aname]
2346            for name in textureData['SH Coeff'][1]:
2347                aname = pfx+name
2348                textureData['SH Coeff'][1][name] = parmDict[aname]
2349                if aname in sigDict:
2350                    SHtextureSig[name] = sigDict[aname]
2351            PrintSHtextureAndSig(textureData,SHtextureSig)
2352        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2353            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2354                textureData,RestraintDict[phase],pFile)
2355                   
2356################################################################################
2357##### Histogram & Phase data
2358################################################################################       
2359                   
2360def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2361    '''Loads the HAP histogram/phase information into dicts
2362
2363    :param dict Phases: phase information
2364    :param dict Histograms: Histogram information
2365    :param bool Print: prints information as it is read
2366    :param file pFile: file object to print to (the default, None causes printing to the console)
2367    :param bool resetRefList: Should the contents of the Reflection List be initialized
2368      on loading. The default, True, initializes the Reflection List as it is loaded.
2369
2370    :returns: (hapVary,hapDict,controlDict)
2371      * hapVary: list of refined variables
2372      * hapDict: dict with refined variables and their values
2373      * controlDict: dict with fixed parameters
2374    '''
2375   
2376    def PrintSize(hapData):
2377        if hapData[0] in ['isotropic','uniaxial']:
2378            line = '\n Size model    : %9s'%(hapData[0])
2379            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2380            if hapData[0] == 'uniaxial':
2381                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2382            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2383            pFile.write(line+'\n')
2384        else:
2385            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2386                (hapData[0],hapData[1][2],hapData[2][2]))
2387            Snames = ['S11','S22','S33','S12','S13','S23']
2388            ptlbls = ' names :'
2389            ptstr =  ' values:'
2390            varstr = ' refine:'
2391            for i,name in enumerate(Snames):
2392                ptlbls += '%12s' % (name)
2393                ptstr += '%12.3f' % (hapData[4][i])
2394                varstr += '%12s' % (str(hapData[5][i]))
2395            pFile.write(ptlbls+'\n')
2396            pFile.write(ptstr+'\n')
2397            pFile.write(varstr+'\n')
2398       
2399    def PrintMuStrain(hapData,SGData):
2400        if hapData[0] in ['isotropic','uniaxial']:
2401            line = '\n Mustrain model: %9s'%(hapData[0])
2402            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2403            if hapData[0] == 'uniaxial':
2404                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2405            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2406            pFile.write(line+'\n')
2407        else:
2408            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2409                (hapData[0],hapData[1][2],hapData[2][2]))
2410            Snames = G2spc.MustrainNames(SGData)
2411            ptlbls = ' names :'
2412            ptstr =  ' values:'
2413            varstr = ' refine:'
2414            for i,name in enumerate(Snames):
2415                ptlbls += '%12s' % (name)
2416                ptstr += '%12.1f' % (hapData[4][i])
2417                varstr += '%12s' % (str(hapData[5][i]))
2418            pFile.write(ptlbls+'\n')
2419            pFile.write(ptstr+'\n')
2420            pFile.write(varstr+'\n')
2421
2422    def PrintHStrain(hapData,SGData):
2423        pFile.write('\n Hydrostatic/elastic strain:\n')
2424        Hsnames = G2spc.HStrainNames(SGData)
2425        ptlbls = ' names :'
2426        ptstr =  ' values:'
2427        varstr = ' refine:'
2428        for i,name in enumerate(Hsnames):
2429            ptlbls += '%12s' % (name)
2430            ptstr += '%12.4g' % (hapData[0][i])
2431            varstr += '%12s' % (str(hapData[1][i]))
2432        pFile.write(ptlbls+'\n')
2433        pFile.write(ptstr+'\n')
2434        pFile.write(varstr+'\n')
2435
2436    def PrintSHPO(hapData):
2437        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2438        ptlbls = ' names :'
2439        ptstr =  ' values:'
2440        for item in hapData[5]:
2441            ptlbls += '%12s'%(item)
2442            ptstr += '%12.3f'%(hapData[5][item]) 
2443        pFile.write(ptlbls+'\n')
2444        pFile.write(ptstr+'\n')
2445   
2446    def PrintBabinet(hapData):
2447        pFile.write('\n Babinet form factor modification:\n')
2448        ptlbls = ' names :'
2449        ptstr =  ' values:'
2450        varstr = ' refine:'
2451        for name in ['BabA','BabU']:
2452            ptlbls += '%12s' % (name)
2453            ptstr += '%12.6f' % (hapData[name][0])
2454            varstr += '%12s' % (str(hapData[name][1]))
2455        pFile.write(ptlbls+'\n')
2456        pFile.write(ptstr+'\n')
2457        pFile.write(varstr+'\n')
2458       
2459    hapDict = {}
2460    hapVary = []
2461    controlDict = {}
2462   
2463    for phase in Phases:
2464        HistoPhase = Phases[phase]['Histograms']
2465        SGData = Phases[phase]['General']['SGData']
2466        cell = Phases[phase]['General']['Cell'][1:7]
2467        A = G2lat.cell2A(cell)
2468        if Phases[phase]['General'].get('Modulated',False):
2469            SSGData = Phases[phase]['General']['SSGData']
2470            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2471        pId = Phases[phase]['pId']
2472        for histogram in Histograms:
2473            if histogram not in HistoPhase and phase in Histograms[histogram]['Reflection Lists']:
2474                #remove previously created reflection list if histogram is removed from phase
2475                #print("removing ",phase,"from",histogram)
2476                del Histograms[histogram]['Reflection Lists'][phase]
2477        histoList = list(HistoPhase.keys())
2478        histoList.sort()
2479        for histogram in histoList:
2480            try:
2481                Histogram = Histograms[histogram]
2482            except KeyError:                       
2483                #skip if histogram not included e.g. in a sequential refinement
2484                continue
2485            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2486                if phase in Histograms[histogram]['Reflection Lists']:
2487                    del Histograms[histogram]['Reflection Lists'][phase]
2488                continue
2489            hapData = HistoPhase[histogram]
2490            hId = Histogram['hId']
2491            if 'PWDR' in histogram:
2492                limits = Histogram['Limits'][1]
2493                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2494                if 'C' in inst['Type'][1]:
2495                    try:
2496                        wave = inst['Lam'][1]
2497                    except KeyError:
2498                        wave = inst['Lam1'][1]
2499                    dmin = wave/(2.0*sind(limits[1]/2.0))
2500                elif 'T' in inst['Type'][0]:
2501                    dmin = limits[0]/inst['difC'][1]
2502                else:
2503                    wave = inst['Lam'][1]
2504                    dmin = wave/(2.0*sind(limits[1]/2.0))
2505                pfx = str(pId)+':'+str(hId)+':'
2506                if Phases[phase]['General']['doPawley']:
2507                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2508                    hapDict[pfx+'newLeBail'] = True
2509                    Tmin = G2lat.Dsp2pos(inst,dmin)
2510                    if 'T' in inst['Type'][1]:
2511                        limits[0] = max(limits[0],Tmin)
2512                    else:
2513                        limits[1] = min(limits[1],Tmin)
2514                else:
2515                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2516                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2517                if Phases[phase]['General']['Type'] == 'magnetic':
2518                    dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.))
2519                for item in ['Scale','Extinction']:
2520                    hapDict[pfx+item] = hapData[item][0]
2521                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2522                        hapVary.append(pfx+item)
2523                names = G2spc.HStrainNames(SGData)
2524                HSvals = []
2525                for i,name in enumerate(names):
2526                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2527                    HSvals.append(hapDict[pfx+name])
2528                    if hapData['HStrain'][1][i]:
2529                        hapVary.append(pfx+name)
2530                if 'Layer Disp' in hapData:
2531                    hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0]
2532                    if hapData['Layer Disp'][1]:
2533                        hapVary.append(pfx+'LayerDisp')
2534                else:
2535                    hapDict[pfx+'LayerDisp'] = 0.0
2536                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2537                if hapData['Pref.Ori.'][0] == 'MD':
2538                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2539                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2540                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2541                        hapVary.append(pfx+'MD')
2542                else:                           #'SH' spherical harmonics
2543                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2544                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2545                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2546                    controlDict[pfx+'SHhkl'] = []
2547                    try: #patch for old Pref.Ori. items
2548                        controlDict[pfx+'SHtoler'] = 0.1
2549                        if hapData['Pref.Ori.'][6][0] != '':
2550                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2551                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2552                    except IndexError:
2553                        pass
2554                    for item in hapData['Pref.Ori.'][5]:
2555                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2556                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2557                            hapVary.append(pfx+item)
2558                for item in ['Mustrain','Size']:
2559                    controlDict[pfx+item+'Type'] = hapData[item][0]
2560                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2561                    if hapData[item][2][2]:
2562                        hapVary.append(pfx+item+';mx')
2563                    if hapData[item][0] in ['isotropic','uniaxial']:
2564                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2565                        if hapData[item][2][0]:
2566                            hapVary.append(pfx+item+';i')
2567                        if hapData[item][0] == 'uniaxial':
2568                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2569                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2570                            if hapData[item][2][1]:
2571                                hapVary.append(pfx+item+';a')
2572                    else:       #generalized for mustrain or ellipsoidal for size
2573                        Nterms = len(hapData[item][4])
2574                        if item == 'Mustrain':
2575                            names = G2spc.MustrainNames(SGData)
2576                            pwrs = []
2577                            for name in names:
2578                                h,k,l = name[1:]
2579                                pwrs.append([int(h),int(k),int(l)])
2580                            controlDict[pfx+'MuPwrs'] = pwrs
2581                        for i in range(Nterms):
2582                            sfx = ';'+str(i)
2583                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2584                            if hapData[item][5][i]:
2585                                hapVary.append(pfx+item+sfx)
2586                if Phases[phase]['General']['Type'] != 'magnetic':
2587                    for bab in ['BabA','BabU']:
2588                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2589                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2590                            hapVary.append(pfx+bab)
2591                               
2592                if Print: 
2593                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2594                    pFile.write(135*'='+'\n')
2595                    if hapDict[pfx+'LeBail']:
2596                        pFile.write(' Perform LeBail extraction\n')                     
2597                    else:
2598                        pFile.write(' Phase fraction  : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2599                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2600                        if hapData['Pref.Ori.'][0] == 'MD':
2601                            Ax = hapData['Pref.Ori.'][3]
2602                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2603                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2604                        else: #'SH' for spherical harmonics
2605                            PrintSHPO(hapData['Pref.Ori.'])
2606                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2607                    PrintSize(hapData['Size'])
2608                    PrintMuStrain(hapData['Mustrain'],SGData)
2609                    PrintHStrain(hapData['HStrain'],SGData)
2610                    if 'Layer Disp' in hapData:
2611                        pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1]))
2612                    if Phases[phase]['General']['Type'] != 'magnetic':
2613                        if hapData['Babinet']['BabA'][0]:
2614                            PrintBabinet(hapData['Babinet'])
2615                if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2616                    hapData['newLeBail'] = True
2617                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2618                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2619                        hapData['newLeBail'] = False
2620                    refList = []
2621#                    Uniq = []
2622#                    Phi = []
2623                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2624                    if Phases[phase]['General'].get('Modulated',False):
2625                        ifSuper = True
2626                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2627                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2628                        for h,k,l,m,d in HKLd:
2629                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2630                            mul *= 2      # for powder overlap of Friedel pairs
2631                            if m or not ext or useExt:
2632                                if 'C' in inst['Type'][0]:
2633                                    pos = G2lat.Dsp2pos(inst,d)
2634                                    if limits[0] < pos < limits[1]:
2635                                        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])
2636                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2637#                                        Uniq.append(uniq)
2638#                                        Phi.append(phi)
2639                                elif 'T' in inst['Type'][0]:
2640                                    pos = G2lat.Dsp2pos(inst,d)
2641                                    if limits[0] < pos < limits[1]:
2642                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2643                                        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])
2644                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2645                                        #TODO - if tabulated put alp & bet in here
2646#                                        Uniq.append(uniq)
2647#                                        Phi.append(phi)
2648                                elif 'B' in inst['Type'][0]:
2649                                    pos = G2lat.Dsp2pos(inst,d)
2650                                    if limits[0] < pos < limits[1]:
2651                                        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])
2652                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2653                    else:
2654                        ifSuper = False
2655                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2656                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2657                        for h,k,l,d in HKLd:
2658                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2659                            if ext and 'N' in inst['Type'][0] and  'MagSpGrp' in SGData:
2660                                ext = G2spc.checkMagextc([h,k,l],SGData)
2661                            mul *= 2      # for powder overlap of Friedel pairs
2662                            if ext and not useExt:
2663                                continue
2664                            if 'C' in inst['Type'][0]:
2665                                pos = G2lat.Dsp2pos(inst,d)
2666                                if limits[0] < pos < limits[1]:
2667                                    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])
2668                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2669#                                    Uniq.append(uniq)
2670#                                    Phi.append(phi)
2671                            elif 'T' in inst['Type'][0]:
2672                                pos = G2lat.Dsp2pos(inst,d)
2673                                if limits[0] < pos < limits[1]:
2674                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2675                                    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])
2676                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2677#                                    Uniq.append(uniq)
2678#                                    Phi.append(phi)
2679                            elif 'B' in inst['Type'][0]:
2680                                pos = G2lat.Dsp2pos(inst,d)
2681                                if limits[0] < pos < limits[1]:
2682                                    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])
2683                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2684                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2685            elif 'HKLF' in histogram:
2686                inst = Histogram['Instrument Parameters'][0]
2687                hId = Histogram['hId']
2688                hfx = ':%d:'%(hId)
2689                for item in inst:
2690                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2691                    hapDict[hfx+item] = inst[item][1]
2692                pfx = str(pId)+':'+str(hId)+':'
2693                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2694                if hapData['Scale'][1]:
2695                    hapVary.append(pfx+'Scale')
2696                               
2697                extApprox,extType,extParms = hapData['Extinction']
2698                controlDict[pfx+'EType'] = extType
2699                controlDict[pfx+'EApprox'] = extApprox
2700                if 'C' in inst['Type'][0]:
2701                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2702                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2703                if 'Primary' in extType:
2704                    Ekey = ['Ep',]
2705                elif 'I & II' in extType:
2706                    Ekey = ['Eg','Es']
2707                elif 'Secondary Type II' == extType:
2708                    Ekey = ['Es',]
2709                elif 'Secondary Type I' == extType:
2710                    Ekey = ['Eg',]
2711                else:   #'None'
2712                    Ekey = []
2713                for eKey in Ekey:
2714                    hapDict[pfx+eKey] = extParms[eKey][0]
2715                    if extParms[eKey][1]:
2716                        hapVary.append(pfx+eKey)
2717                for bab in ['BabA','BabU']:
2718                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2719                    if hapData['Babinet'][bab][1]:
2720                        hapVary.append(pfx+bab)
2721                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2722                if len(Twins) == 1:
2723                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2724                    if hapData.get('Flack',[0,False])[1]:
2725                        hapVary.append(pfx+'Flack')
2726                sumTwFr = 0.
2727                controlDict[pfx+'TwinLaw'] = []
2728                controlDict[pfx+'TwinInv'] = []
2729                NTL = 0           
2730                for it,twin in enumerate(Twins):
2731                    if 'bool' in str(type(twin[0])):
2732                        controlDict[pfx+'TwinInv'].append(twin[0])
2733                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2734                    else:
2735                        NTL += 1
2736                        controlDict[pfx+'TwinInv'].append(False)
2737                        controlDict[pfx+'TwinLaw'].append(twin[0])
2738                    if it:
2739                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2740                        sumTwFr += twin[1]
2741                    else:
2742                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2743                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2744                    if Twins[0][1][1]:
2745                        hapVary.append(pfx+'TwinFr:'+str(it))
2746                controlDict[pfx+'NTL'] = NTL
2747                #need to make constraint on TwinFr
2748                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2749                if len(Twins) > 1:    #force sum to unity
2750                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2751                if Print: 
2752                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2753                    pFile.write(135*'='+'\n')
2754                    pFile.write(' Scale factor     : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2755                    if extType != 'None':
2756                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2757                        text = ' Parameters       :'
2758                        for eKey in Ekey:
2759                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2760                        pFile.write(text+'\n')
2761                    if hapData['Babinet']['BabA'][0]:
2762                        PrintBabinet(hapData['Babinet'])
2763                    if not SGData['SGInv'] and len(Twins) == 1:
2764                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2765                    if len(Twins) > 1:
2766                        for it,twin in enumerate(Twins):
2767                            if 'bool' in str(type(twin[0])):
2768                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2769                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2770                            else:
2771                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2772                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2773                       
2774                Histogram['Reflection Lists'] = phase       
2775               
2776    return hapVary,hapDict,controlDict
2777   
2778def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,Print=True,pFile=None):
2779    'needs a doc string'
2780   
2781    def PrintSizeAndSig(hapData,sizeSig):
2782        line = '\n Size model:     %9s'%(hapData[0])
2783        refine = False
2784        if hapData[0] in ['isotropic','uniaxial']:
2785            line += ' equatorial:%12.5f'%(hapData[1][0])
2786            if sizeSig[0][0]:
2787                line += ', sig:%8.4f'%(sizeSig[0][0])
2788                refine = True
2789            if hapData[0] == 'uniaxial':
2790                line += ' axial:%12.4f'%(hapData[1][1])
2791                if sizeSig[0][1]:
2792                    refine = True
2793                    line += ', sig:%8.4f'%(sizeSig[0][1])
2794            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2795            if sizeSig[0][2]:
2796                refine = True
2797                line += ', sig:%8.4f'%(sizeSig[0][2])
2798            if refine:
2799                pFile.write(line+'\n')
2800        else:
2801            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2802            if sizeSig[0][2]:
2803                refine = True
2804                line += ', sig:%8.4f'%(sizeSig[0][2])
2805            Snames = ['S11','S22','S33','S12','S13','S23']
2806            ptlbls = ' name  :'
2807            ptstr =  ' value :'
2808            sigstr = ' sig   :'
2809            for i,name in enumerate(Snames):
2810                ptlbls += '%12s' % (name)
2811                ptstr += '%12.6f' % (hapData[4][i])
2812                if sizeSig[1][i]:
2813                    refine = True
2814                    sigstr += '%12.6f' % (sizeSig[1][i])
2815                else:
2816                    sigstr += 12*' '
2817            if refine:
2818                pFile.write(line+'\n')
2819                pFile.write(ptlbls+'\n')
2820                pFile.write(ptstr+'\n')
2821                pFile.write(sigstr+'\n')
2822       
2823    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2824        line = '\n Mustrain model: %9s\n'%(hapData[0])
2825        refine = False
2826        if hapData[0] in ['isotropic','uniaxial']:
2827            line += ' equatorial:%12.1f'%(hapData[1][0])
2828            if mustrainSig[0][0]:
2829                line += ', sig:%8.1f'%(mustrainSig[0][0])
2830                refine = True
2831            if hapData[0] == 'uniaxial':
2832                line += ' axial:%12.1f'%(hapData[1][1])
2833                if mustrainSig[0][1]:
2834                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2835            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2836            if mustrainSig[0][2]:
2837                refine = True
2838                line += ', sig:%8.3f'%(mustrainSig[0][2])
2839            if refine:
2840                pFile.write(line+'\n')
2841        else:
2842            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2843            if mustrainSig[0][2]:
2844                refine = True
2845                line += ', sig:%8.3f'%(mustrainSig[0][2])
2846            Snames = G2spc.MustrainNames(SGData)
2847            ptlbls = ' name  :'
2848            ptstr =  ' value :'
2849            sigstr = ' sig   :'
2850            for i,name in enumerate(Snames):
2851                ptlbls += '%12s' % (name)
2852                ptstr += '%12.1f' % (hapData[4][i])
2853                if mustrainSig[1][i]:
2854                    refine = True
2855                    sigstr += '%12.1f' % (mustrainSig[1][i])
2856                else:
2857                    sigstr += 12*' '
2858            if refine:
2859                pFile.write(line+'\n')
2860                pFile.write(ptlbls+'\n')
2861                pFile.write(ptstr+'\n')
2862                pFile.write(sigstr+'\n')
2863           
2864    def PrintHStrainAndSig(hapData,strainSig,SGData):
2865        Hsnames = G2spc.HStrainNames(SGData)
2866        ptlbls = ' name  :'
2867        ptstr =  ' value :'
2868        sigstr = ' sig   :'
2869        refine = False
2870        for i,name in enumerate(Hsnames):
2871            ptlbls += '%12s' % (name)
2872            ptstr += '%12.4g' % (hapData[0][i])
2873            if name in strainSig:
2874                refine = True
2875                sigstr += '%12.4g' % (strainSig[name])
2876            else:
2877                sigstr += 12*' '
2878        if refine:
2879            pFile.write('\n Hydrostatic/elastic strain:\n')
2880            pFile.write(ptlbls+'\n')
2881            pFile.write(ptstr+'\n')
2882            pFile.write(sigstr+'\n')
2883       
2884    def PrintSHPOAndSig(pfx,hapData,POsig):
2885        Tindx = 1.0
2886        Tvar = 0.0
2887        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2888        ptlbls = ' names :'
2889        ptstr =  ' values:'
2890        sigstr = ' sig   :'
2891        for item in hapData[5]:
2892            ptlbls += '%12s'%(item)
2893            ptstr += '%12.3f'%(hapData[5][item])
2894            l = 2.0*eval(item.strip('C'))[0]+1
2895            Tindx += hapData[5][item]**2/l
2896            if pfx+item in POsig:
2897                Tvar += (2.*hapData[5][item]*POsig[pfx+item]/l)**2
2898                sigstr += '%12.3f'%(POsig[pfx+item])
2899            else:
2900                sigstr += 12*' ' 
2901        pFile.write(ptlbls+'\n')
2902        pFile.write(ptstr+'\n')
2903        pFile.write(sigstr+'\n')
2904        pFile.write('\n Texture index J = %.3f(%d)\n'%(Tindx,int(1000*np.sqrt(Tvar))))
2905       
2906    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2907        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2908        text = ''
2909        for item in hapData[2]:
2910            if pfx+item in ScalExtSig:
2911                text += '       %s: '%(item)
2912                text += '%12.2e'%(hapData[2][item][0])
2913                if pfx+item in ScalExtSig:
2914                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2915        pFile.write(text+'\n')   
2916       
2917    def PrintBabinetAndSig(pfx,hapData,BabSig):
2918        pFile.write('\n Babinet form factor modification:\n')
2919        ptlbls = ' names :'
2920        ptstr =  ' values:'
2921        sigstr = ' sig   :'
2922        for item in hapData:
2923            ptlbls += '%12s'%(item)
2924            ptstr += '%12.3f'%(hapData[item][0])
2925            if pfx+item in BabSig:
2926                sigstr += '%12.3f'%(BabSig[pfx+item])
2927            else:
2928                sigstr += 12*' ' 
2929        pFile.write(ptlbls+'\n')
2930        pFile.write(ptstr+'\n')
2931        pFile.write(sigstr+'\n')
2932       
2933    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2934        pFile.write('\n Twin Law fractions :\n')
2935        ptlbls = ' names :'
2936        ptstr =  ' values:'
2937        sigstr = ' sig   :'
2938        for it,item in enumerate(twinData):
2939            ptlbls += '     twin #%d'%(it)
2940            if it:
2941                ptstr += '%12.3f'%(item[1])
2942            else:
2943                ptstr += '%12.3f'%(item[1][0])
2944            if pfx+'TwinFr:'+str(it) in TwinSig:
2945                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2946            else:
2947                sigstr += 12*' ' 
2948        pFile.write(ptlbls+'\n')
2949        pFile.write(ptstr+'\n')
2950        pFile.write(sigstr+'\n')
2951       
2952   
2953    PhFrExtPOSig = {}
2954    SizeMuStrSig = {}
2955    ScalExtSig = {}
2956    BabSig = {}
2957    TwinFrSig = {}
2958    wtFrSum = {}
2959    for phase in Phases:
2960        HistoPhase = Phases[phase]['Histograms']
2961        General = Phases[phase]['General']
2962        SGData = General['SGData']
2963        pId = Phases[phase]['pId']
2964        histoList = list(HistoPhase.keys())
2965        histoList.sort()
2966        for histogram in histoList:
2967            try:
2968                Histogram = Histograms[histogram]
2969            except KeyError:                       
2970                #skip if histogram not included e.g. in a sequential refinement
2971                continue
2972            if not Phases[phase]['Histograms'][histogram]['Use']:
2973                #skip if phase absent from this histogram
2974                continue
2975            hapData = HistoPhase[histogram]
2976            hId = Histogram['hId']
2977            pfx = str(pId)+':'+str(hId)+':'
2978            if hId not in wtFrSum:
2979                wtFrSum[hId] = 0.
2980            if 'PWDR' in histogram:
2981                for item in ['Scale','Extinction']:
2982                    hapData[item][0] = parmDict[pfx+item]
2983                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2984                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2985                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2986                if hapData['Pref.Ori.'][0] == 'MD':
2987                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2988                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2989                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2990                else:                           #'SH' spherical harmonics
2991                    for item in hapData['Pref.Ori.'][5]:
2992                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2993                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2994                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2995                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2996                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2997                    pfx+'HStrain':{}})                 
2998                for item in ['Mustrain','Size']:
2999                    hapData[item][1][2] = parmDict[pfx+item+';mx']
3000#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
3001                    if pfx+item+';mx' in sigDict:
3002                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
3003                    if hapData[item][0] in ['isotropic','uniaxial']:                   
3004                        hapData[item][1][0] = parmDict[pfx+item+';i']
3005                        if item == 'Size':
3006                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
3007                        if pfx+item+';i' in sigDict: 
3008                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
3009                        if hapData[item][0] == 'uniaxial':
3010                            hapData[item][1][1] = parmDict[pfx+item+';a']
3011                            if item == 'Size':
3012                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
3013                            if pfx+item+';a' in sigDict:
3014                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
3015                    else:       #generalized for mustrain or ellipsoidal for size
3016                        Nterms = len(hapData[item][4])
3017                        for i in range(Nterms):
3018                            sfx = ';'+str(i)
3019                            hapData[item][4][i] = parmDict[pfx+item+sfx]
3020                            if pfx+item+sfx in sigDict:
3021                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
3022                names = G2spc.HStrainNames(SGData)
3023                for i,name in enumerate(names):
3024                    hapData['HStrain'][0][i] = parmDict[pfx+name]
3025                    if pfx+name in sigDict:
3026                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
3027                if 'Layer Disp' in hapData:
3028                    hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp']
3029                    if pfx+'LayerDisp' in sigDict:
3030                        SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp']
3031                if Phases[phase]['General']['Type'] != 'magnetic':
3032                    for name in ['BabA','BabU']:
3033                        hapData['Babinet'][name][0] = parmDict[pfx+name]
3034                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
3035                            BabSig[pfx+name] = sigDict[pfx+name]               
3036               
3037            elif 'HKLF' in histogram:
3038                for item in ['Scale','Flack']:
3039                    if parmDict.get(pfx+item):
3040                        hapData[item][0] = parmDict[pfx+item]
3041                        if pfx+item in sigDict:
3042                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3043                for item in ['Ep','Eg','Es']:
3044                    if parmDict.get(pfx+item):
3045                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
3046                        if pfx+item in sigDict:
3047                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3048                for item in ['BabA','BabU']:
3049                    hapData['Babinet'][item][0] = parmDict[pfx+item]
3050                    if pfx+item in sigDict:
3051                        BabSig[pfx+item] = sigDict[pfx+item]
3052                if 'Twins' in hapData:
3053                    it = 1
3054                    sumTwFr = 0.
3055                    while True:
3056                        try:
3057                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
3058                            if pfx+'TwinFr:'+str(it) in sigDict:
3059                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
3060                            if it:
3061                                sumTwFr += hapData['Twins'][it][1]
3062                            it += 1
3063                        except KeyError:
3064                            break
3065                    hapData['Twins'][0][1][0] = 1.-sumTwFr
3066
3067    if Print:
3068        for phase in Phases:
3069            HistoPhase = Phases[phase]['Histograms']
3070            General = Phases[phase]['General']
3071            SGData = General['SGData']
3072            pId = Phases[phase]['pId']
3073            histoList = list(HistoPhase.keys())
3074            histoList.sort()
3075            for histogram in histoList:
3076                try:
3077                    Histogram = Histograms[histogram]
3078                except KeyError:                       
3079                    #skip if histogram not included e.g. in a sequential refinement
3080                    continue
3081                hapData = HistoPhase[histogram]
3082                hId = Histogram['hId']
3083                Histogram['Residuals'][str(pId)+'::Name'] = phase
3084                pfx = str(pId)+':'+str(hId)+':'
3085                hfx = ':%s:'%(hId)
3086                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
3087                    continue
3088                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
3089                pFile.write(135*'='+'\n')
3090                Inst = Histogram['Instrument Parameters'][0]
3091                if 'PWDR' in histogram:
3092                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
3093                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
3094                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
3095                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
3096                   
3097                    if parmDict[pfx+'LeBail']:
3098                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
3099                    else:
3100                        # if calcControls != None:    #skipped in seqRefine
3101                        #     if 'X'in Inst['Type'][0]:
3102                        #         PrintFprime(calcControls['FFtables'],hfx,pFile)
3103                        #     elif 'NC' in Inst['Type'][0]:
3104                        #         PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3105                        if pfx+'Scale' in PhFrExtPOSig:
3106                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
3107                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
3108                            pFile.write(' Phase fraction  : %10.5g, sig %10.5g Weight fraction  : %8.5f, sig %10.5f\n'%
3109                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
3110                        if pfx+'Extinction' in PhFrExtPOSig:
3111                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
3112                        if hapData['Pref.Ori.'][0] == 'MD':
3113                            if pfx+'MD' in PhFrExtPOSig:
3114                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
3115                        else:
3116                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
3117                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
3118                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
3119                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
3120                    if pfx+'LayerDisp' in SizeMuStrSig:
3121                        pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp']))           
3122                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
3123                        if len(BabSig):
3124                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3125                   
3126                elif 'HKLF' in histogram:
3127                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
3128                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
3129                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
3130                    if calcControls != None:    #skipped in seqRefine
3131                        if 'X'in Inst['Type'][0]:
3132                            PrintFprime(calcControls['FFtables'],hfx,pFile)
3133                        elif 'NC' in Inst['Type'][0]:
3134                            PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3135                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
3136                    if pfx+'Scale' in ScalExtSig:
3137                        pFile.write(' Scale factor : %10.4g, sig %10.4g\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
3138                    if hapData['Extinction'][0] != 'None':
3139                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
3140                    if len(BabSig):
3141                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3142                    if pfx+'Flack' in ScalExtSig:
3143                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
3144                    if len(TwinFrSig):
3145                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
3146
3147################################################################################
3148##### Histogram data
3149################################################################################       
3150                   
3151def GetHistogramData(Histograms,Print=True,pFile=None):
3152    'needs a doc string'
3153   
3154    def GetBackgroundParms(hId,Background):
3155        Back = Background[0]
3156        DebyePeaks = Background[1]
3157        bakType,bakFlag = Back[:2]
3158        backVals = Back[3:]
3159        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
3160        backDict = dict(zip(backNames,backVals))
3161        backVary = []
3162        if bakFlag:
3163            backVary = backNames
3164        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
3165        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
3166        debyeDict = {}
3167        debyeList = []
3168        for i in range(DebyePeaks['nDebye']):
3169            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
3170            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
3171            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
3172        debyeVary = []
3173        for item in debyeList:
3174            if item[1]:
3175                debyeVary.append(item[0])
3176        backDict.update(debyeDict)
3177        backVary += debyeVary
3178        peakDict = {}
3179        peakList = []
3180        for i in range(DebyePeaks['nPeaks']):
3181            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
3182                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
3183            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
3184            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
3185        peakVary = []
3186        for item in peakList:
3187            if item[1]:
3188                peakVary.append(item[0])
3189        backDict.update(peakDict)
3190        backVary += peakVary
3191        if 'background PWDR' in Background[1]:
3192            backDict[':'+str(hId)+':Back File'] = Background[1]['background PWDR'][0]
3193            backDict[':'+str(hId)+':BF mult'] = Background[1]['background PWDR'][1]
3194            try:
3195                if Background[1]['background PWDR'][0] and Background[1]['background PWDR'][2]:
3196                    backVary.append(':'+str(hId)+':BF mult')
3197            except IndexError:  # old version without refine flag
3198                pass
3199        return bakType,backDict,backVary           
3200       
3201    def GetInstParms(hId,Inst):
3202        #patch
3203        if 'Z' not in Inst:
3204            Inst['Z'] = [0.0,0.0,False]
3205        dataType = Inst['Type'][0]
3206        instDict = {}
3207        insVary = []
3208        pfx = ':'+str(hId)+':'
3209        insKeys = list(Inst.keys())
3210        insKeys.sort()
3211        for item in insKeys:
3212            if item in 'Azimuth':
3213                continue
3214            insName = pfx+item
3215            instDict[insName] = Inst[item][1]
3216            if len(Inst[item]) > 2 and Inst[item][2]:
3217                insVary.append(insName)
3218        if 'C' in dataType:
3219            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
3220        elif 'T' in dataType:   #trap zero alp, bet coeff.
3221            if not instDict[pfx+'alpha']:
3222                instDict[pfx+'alpha'] = 1.0
3223            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3224                instDict[pfx+'beta-1'] = 1.0
3225        elif 'B' in dataType:   #trap zero alp, bet coeff.
3226            if not instDict[pfx+'alpha-0'] and not instDict[pfx+'alpha-1']:
3227                instDict[pfx+'alpha-1'] = 1.0
3228            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3229                instDict[pfx+'beta-1'] = 1.0
3230        return dataType,instDict,insVary
3231       
3232    def GetSampleParms(hId,Sample):
3233        sampVary = []
3234        hfx = ':'+str(hId)+':'       
3235        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
3236            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi'],hfx+'Azimuth':Sample['Azimuth']}
3237        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
3238            if key in Sample:
3239                sampDict[hfx+key] = Sample[key]
3240        Type = Sample['Type']
3241        if 'Bragg' in Type:             #Bragg-Brentano
3242            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3243                sampDict[hfx+item] = Sample[item][0]
3244                if Sample[item][1]:
3245                    sampVary.append(hfx+item)
3246        elif 'Debye' in Type:        #Debye-Scherrer
3247            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3248                sampDict[hfx+item] = Sample[item][0]
3249                if Sample[item][1]:
3250                    sampVary.append(hfx+item)
3251        return Type,sampDict,sampVary
3252       
3253    def PrintBackground(Background):
3254        Back = Background[0]
3255        DebyePeaks = Background[1]
3256        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
3257        line = ' Coefficients: '
3258        for i,back in enumerate(Back[3:]):
3259            line += '%10.3f'%(back)
3260            if i and not i%10:
3261                line += '\n'+15*' '
3262        pFile.write(line+'\n')
3263        if DebyePeaks['nDebye']:
3264            pFile.write('\n Debye diffuse scattering coefficients\n')
3265            parms = ['DebyeA','DebyeR','DebyeU']
3266            line = ' names :  '
3267            for parm in parms:
3268                line += '%8s refine?'%(parm)
3269            pFile.write(line+'\n')
3270            for j,term in enumerate(DebyePeaks['debyeTerms']):
3271                line = ' term'+'%2d'%(j)+':'
3272                for i in range(3):
3273                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3274                pFile.write(line+'\n')
3275        if DebyePeaks['nPeaks']:
3276            pFile.write('\n Single peak coefficients\n')
3277            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3278            line = ' names :  '
3279            for parm in parms:
3280                line += '%8s refine?'%(parm)
3281            pFile.write(line+'\n')
3282            for j,term in enumerate(DebyePeaks['peaksList']):
3283                line = ' peak'+'%2d'%(j)+':'
3284                for i in range(4):
3285                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3286                pFile.write(line+'\n')
3287        if 'background PWDR' in DebyePeaks:
3288            try:
3289                pFile.write(' Fixed background file: %s; mult: %.3f  Refine? %s\n'%(DebyePeaks['background PWDR'][0],
3290                    DebyePeaks['background PWDR'][1],DebyePeaks['background PWDR'][2]))
3291            except IndexError:  #old version without refine flag
3292                pass
3293       
3294    def PrintInstParms(Inst):
3295        pFile.write('\n Instrument Parameters:\n')
3296        insKeys = list(Inst.keys())
3297        insKeys.sort()
3298        iBeg = 0
3299        Ok = True
3300        while Ok:
3301            ptlbls = ' name  :'
3302            ptstr =  ' value :'
3303            varstr = ' refine:'
3304            iFin = min(iBeg+9,len(insKeys))
3305            for item in insKeys[iBeg:iFin]:
3306                if item not in ['Type','Source']:
3307                    ptlbls += '%12s' % (item)
3308                    ptstr += '%12.6f' % (Inst[item][1])
3309                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3310                        varstr += 12*' '
3311                    else:
3312                        varstr += '%12s' % (str(bool(Inst[item][2])))
3313            pFile.write(ptlbls+'\n')
3314            pFile.write(ptstr+'\n')
3315            pFile.write(varstr+'\n')
3316            iBeg = iFin
3317            if iBeg == len(insKeys):
3318                Ok = False
3319            else:
3320                pFile.write('\n')
3321       
3322    def PrintSampleParms(Sample):
3323        pFile.write('\n Sample Parameters:\n')
3324        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3325            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3326        ptlbls = ' name  :'
3327        ptstr =  ' value :'
3328        varstr = ' refine:'
3329        if 'Bragg' in Sample['Type']:
3330            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3331                ptlbls += '%14s'%(item)
3332                ptstr += '%14.4f'%(Sample[item][0])
3333                varstr += '%14s'%(str(bool(Sample[item][1])))
3334           
3335        elif 'Debye' in Type:        #Debye-Scherrer
3336            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3337                ptlbls += '%14s'%(item)
3338                ptstr += '%14.4f'%(Sample[item][0])
3339                varstr += '%14s'%(str(bool(Sample[item][1])))
3340
3341        pFile.write(ptlbls+'\n')
3342        pFile.write(ptstr+'\n')
3343        pFile.write(varstr+'\n')
3344       
3345    histDict = {}
3346    histVary = []
3347    controlDict = {}
3348    histoList = list(Histograms.keys())
3349    histoList.sort()
3350    for histogram in histoList:
3351        if 'PWDR' in histogram:
3352            Histogram = Histograms[histogram]
3353            hId = Histogram['hId']
3354            pfx = ':'+str(hId)+':'
3355            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3356            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3357            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3358            for excl in controlDict[pfx+'Exclude']:
3359                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3360            if controlDict[pfx+'Exclude']:
3361                ma.mask_rows(Histogram['Data'])
3362            Background = Histogram['Background']
3363            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3364            controlDict[pfx+'bakType'] = Type
3365            histDict.update(bakDict)
3366            histVary += bakVary
3367           
3368            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3369            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3370                G2fil.G2Print ('Warning: tabulated profile coefficients are ignored')
3371            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3372            controlDict[pfx+'histType'] = Type
3373            if 'XC' in Type or 'XB' in Type:
3374                if pfx+'Lam1' in instDict:
3375                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3376                else:
3377                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3378            histDict.update(instDict)
3379            histVary += insVary
3380           
3381            Sample = Histogram['Sample Parameters']
3382            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3383            controlDict[pfx+'instType'] = Type
3384            histDict.update(sampDict)
3385            histVary += sampVary
3386           
3387   
3388            if Print: 
3389                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3390                pFile.write(135*'='+'\n')
3391                Units = {'C':' deg','T':' msec','B':' deg'}
3392                units = Units[controlDict[pfx+'histType'][2]]
3393                Limits = controlDict[pfx+'Limits']
3394                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3395                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3396                if len(controlDict[pfx+'Exclude']):
3397                    excls = controlDict[pfx+'Exclude']
3398                    for excl in excls:
3399                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3400                PrintSampleParms(Sample)
3401                PrintInstParms(Inst[0])
3402                PrintBackground(Background)
3403        elif 'HKLF' in histogram:
3404            Histogram = Histograms[histogram]
3405            hId = Histogram['hId']
3406            pfx = ':'+str(hId)+':'
3407            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3408            Inst = Histogram['Instrument Parameters'][0]
3409            controlDict[pfx+'histType'] = Inst['Type'][0]
3410            if 'X' in Inst['Type'][0]:
3411                histDict[pfx+'Lam'] = Inst['Lam'][1]
3412                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3413            elif 'NC' in Inst['Type'][0] or 'NB' in Inst['Type'][0]:                   
3414                histDict[pfx+'Lam'] = Inst['Lam'][1]
3415    return histVary,histDict,controlDict
3416   
3417def SetHistogramData(parmDict,sigDict,Histograms,calcControls,Print=True,pFile=None,seq=False):
3418    'Shows histogram data after a refinement'
3419   
3420    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3421        Back = Background[0]
3422        DebyePeaks = Background[1]
3423        lenBack = len(Back[3:])
3424        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3425        for i in range(lenBack):
3426            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3427            if pfx+'Back;'+str(i) in sigDict:
3428                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3429        if DebyePeaks['nDebye']:
3430            for i in range(DebyePeaks['nDebye']):
3431                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3432                for j,name in enumerate(names):
3433                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3434                    if name in sigDict:
3435                        backSig[lenBack+3*i+j] = sigDict[name]           
3436        if DebyePeaks['nPeaks']:
3437            for i in range(DebyePeaks['nPeaks']):
3438                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3439                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3440                for j,name in enumerate(names):
3441                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3442                    if name in sigDict:
3443                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3444        if pfx+'BF mult' in sigDict:
3445            DebyePeaks['background PWDR'][1] = parmDict[pfx+'BF mult']
3446            backSig.append(sigDict[pfx+'BF mult'])
3447               
3448        return backSig
3449       
3450    def SetInstParms(pfx,Inst,parmDict,sigDict):
3451        instSig = {}
3452        insKeys = list(Inst.keys())
3453        insKeys.sort()
3454        for item in insKeys:
3455            insName = pfx+item
3456            Inst[item][1] = parmDict[insName]
3457            if insName in sigDict:
3458                instSig[item] = sigDict[insName]
3459            else:
3460                instSig[item] = 0
3461        return instSig
3462       
3463    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3464        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3465            sampSig = [0 for i in range(5)]
3466            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3467                Sample[item][0] = parmDict[pfx+item]
3468                if pfx+item in sigDict:
3469                    sampSig[i] = sigDict[pfx+item]
3470        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3471            sampSig = [0 for i in range(4)]
3472            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3473                Sample[item][0] = parmDict[pfx+item]
3474                if pfx+item in sigDict:
3475                    sampSig[i] = sigDict[pfx+item]
3476        return sampSig
3477       
3478    def PrintBackgroundSig(Background,backSig):
3479        Back = Background[0]
3480        DebyePeaks = Background[1]
3481        valstr = ' value : '
3482        sigstr = ' sig   : '
3483        refine = False
3484        for i,back in enumerate(Back[3:]):
3485            valstr += '%10.4g'%(back)
3486            if Back[1]:
3487                refine = True
3488                sigstr += '%10.4g'%(backSig[i])
3489            else:
3490                sigstr += 10*' '
3491        if refine:
3492            pFile.write('\n Background function: %s\n'%Back[0])
3493            pFile.write(valstr+'\n')
3494            pFile.write(sigstr+'\n')
3495        if DebyePeaks['nDebye']:
3496            ifAny = False
3497            ptfmt = "%12.3f"
3498            names =  ' names :'
3499            ptstr =  ' values:'
3500            sigstr = ' esds  :'
3501            for item in sigDict:
3502                if 'Debye' in item:
3503                    ifAny = True
3504                    names += '%12s'%(item)
3505                    ptstr += ptfmt%(parmDict[item])
3506                    sigstr += ptfmt%(sigDict[item])
3507            if ifAny:
3508                pFile.write('\n Debye diffuse scattering coefficients\n')
3509                pFile.write(names+'\n')
3510                pFile.write(ptstr+'\n')
3511                pFile.write(sigstr+'\n')
3512        if DebyePeaks['nPeaks']:
3513            pFile.write('\n Single peak coefficients:\n')
3514            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3515            line = ' peak no. '
3516            for parm in parms:
3517                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3518            pFile.write(line+'\n')
3519            for ip in range(DebyePeaks['nPeaks']):
3520                ptstr = ' %4d '%(ip)
3521                for parm in parms:
3522                    fmt = '%14.3f'
3523                    efmt = '%12.3f'
3524                    if parm == 'BkPkpos':
3525                        fmt = '%14.4f'
3526                        efmt = '%12.4f'
3527                    name = pfx+parm+';%d'%(ip)
3528                    ptstr += fmt%(parmDict[name])
3529                    if name in sigDict:
3530                        ptstr += efmt%(sigDict[name])
3531                    else:
3532                        ptstr += 12*' '
3533                pFile.write(ptstr+'\n')
3534        if 'background PWDR' in DebyePeaks and DebyePeaks['background PWDR'][2]:
3535            pFile.write(' Fixed background scale: %.3f(%d)\n'%(DebyePeaks['background PWDR'][1],int(1000*backSig[-1])))
3536        sumBk = np.array(Histogram['sumBk'])
3537        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3538            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3539       
3540    def PrintInstParmsSig(Inst,instSig):
3541        refine = False
3542        insKeys = list(instSig.keys())
3543        insKeys.sort()
3544        iBeg = 0
3545        Ok = True
3546        while Ok:
3547            ptlbls = ' names :'
3548            ptstr =  ' value :'
3549            sigstr = ' sig   :'
3550            iFin = min(iBeg+9,len(insKeys))
3551            for name in insKeys[iBeg:iFin]:
3552                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3553                    ptlbls += '%12s' % (name)
3554                    ptstr += '%12.6f' % (Inst[name][1])
3555                    if instSig[name]:
3556                        refine = True
3557                        sigstr += '%12.6f' % (instSig[name])
3558                    else:
3559                        sigstr += 12*' '
3560            if refine:
3561                pFile.write('\n Instrument Parameters:\n')
3562                pFile.write(ptlbls+'\n')
3563                pFile.write(ptstr+'\n')
3564                pFile.write(sigstr+'\n')
3565            iBeg = iFin
3566            if iBeg == len(insKeys):
3567                Ok = False
3568       
3569    def PrintSampleParmsSig(Sample,sampleSig):
3570        ptlbls = ' names :'
3571        ptstr =  ' values:'
3572        sigstr = ' sig   :'
3573        refine = False
3574        if 'Bragg' in Sample['Type']:
3575            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3576                ptlbls += '%14s'%(item)
3577                ptstr += '%14.4f'%(Sample[item][0])
3578                if sampleSig[i]:
3579                    refine = True
3580                    sigstr += '%14.4f'%(sampleSig[i])
3581                else:
3582                    sigstr += 14*' '
3583           
3584        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3585            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3586                ptlbls += '%14s'%(item)
3587                ptstr += '%14.4f'%(Sample[item][0])
3588                if sampleSig[i]:
3589                    refine = True
3590                    sigstr += '%14.4f'%(sampleSig[i])
3591                else:
3592                    sigstr += 14*' '
3593
3594        if refine:
3595            pFile.write('\n Sample Parameters:\n')
3596            pFile.write(ptlbls+'\n')
3597            pFile.write(ptstr+'\n')
3598            pFile.write(sigstr+'\n')
3599       
3600    histoList = list(Histograms.keys())
3601    histoList.sort()
3602    for histogram in histoList:
3603        if 'PWDR' in histogram:
3604            Histogram = Histograms[histogram]
3605            hId = Histogram['hId']
3606            pfx = ':'+str(hId)+':'
3607            Background = Histogram['Background']
3608            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3609           
3610            Inst = Histogram['Instrument Parameters'][0]
3611            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3612       
3613            Sample = Histogram['Sample Parameters']
3614            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3615
3616            if not seq:
3617                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3618                pFile.write(135*'='+'\n')
3619            pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3620            pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3621                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3622            pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3623                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3624            if Print:
3625                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3626                if calcControls != None:    #skipped in seqRefine
3627                    if 'X' in Inst['Type'][0]:
3628                        PrintFprime(calcControls['FFtables'],pfx,pFile)
3629                    elif 'NC' in Inst['Type'][0]:
3630                        PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3631                PrintSampleParmsSig(Sample,sampSig)
3632                PrintInstParmsSig(Inst,instSig)
3633                PrintBackgroundSig(Background,backSig)
3634               
3635def WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
3636    '''Cribbed version of PrintRBObjPOAndSig but returns lists of strings.
3637    Moved so it can be used in ExportCIF
3638    '''
3639    namstr = '  names :'
3640    valstr = '  values:'
3641    sigstr = '  esds  :'
3642    for i,px in enumerate(['Px:','Py:','Pz:']):
3643        name = pfx+rbfx+px+rbsx
3644        namstr += '%12s'%('Pos '+px[1])
3645        valstr += '%12.5f'%(parmDict[name])
3646        if name in sigDict:
3647            sigstr += '%12.5f'%(sigDict[name])
3648        else:
3649            sigstr += 12*' '
3650    for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
3651        name = pfx+rbfx+po+rbsx
3652        namstr += '%12s'%('Ori '+po[1])
3653        valstr += '%12.5f'%(parmDict[name])
3654        if name in sigDict:
3655            sigstr += '%12.5f'%(sigDict[name])
3656        else:
3657            sigstr += 12*' '
3658    name = pfx+rbfx+'f:'+rbsx
3659    namstr += '%12s'%('Frac')
3660    valstr += '%12.5f'%(parmDict[name])
3661    if name in sigDict:
3662        sigstr += '%12.5f'%(sigDict[name])
3663    else:
3664        sigstr += 12*' '
3665    return (namstr,valstr,sigstr)
3666
3667def WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
3668    '''Cribbed version of PrintRBObjTLSAndSig but returns lists of strings.
3669    Moved so it can be used in ExportCIF
3670    '''
3671    out = []
3672    namstr = '  names :'
3673    valstr = '  values:'
3674    sigstr = '  esds  :'
3675    if 'N' not in TLS:
3676        out.append(' Thermal motion:\n')
3677    if 'T' in TLS:
3678        for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
3679            name = pfx+rbfx+pt+rbsx
3680            namstr += '%12s'%(pt[:3])
3681            valstr += '%12.5f'%(parmDict[name])
3682            if name in sigDict:
3683                sigstr += '%12.5f'%(sigDict[name])
3684            else:
3685                sigstr += 12*' '
3686        out.append(namstr+'\n')
3687        out.append(valstr+'\n')
3688        out.append(sigstr+'\n')
3689    if 'L' in TLS:
3690        namstr = '  names :'
3691        valstr = '  values:'
3692        sigstr = '  esds  :'
3693        for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
3694            name = pfx+rbfx+pt+rbsx
3695            namstr += '%12s'%(pt[:3])
3696            valstr += '%12.3f'%(parmDict[name])
3697            if name in sigDict:
3698                sigstr += '%12.3f'%(sigDict[name])
3699            else:
3700                sigstr += 12*' '
3701        out.append(namstr+'\n')
3702        out.append(valstr+'\n')
3703        out.append(sigstr+'\n')
3704    if 'S' in TLS:
3705        namstr = '  names :'
3706        valstr = '  values:'
3707        sigstr = '  esds  :'
3708        for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
3709            name = pfx+rbfx+pt+rbsx
3710            namstr += '%12s'%(pt[:3])
3711            valstr += '%12.4f'%(parmDict[name])
3712            if name in sigDict:
3713                sigstr += '%12.4f'%(sigDict[name])
3714            else:
3715                sigstr += 12*' '
3716        out.append(namstr+'\n')
3717        out.append(valstr+'\n')
3718        out.append(sigstr+'\n')
3719    if 'U' in TLS:
3720        name = pfx+rbfx+'U:'+rbsx
3721        namstr = '  names :'+'%12s'%('Uiso')
3722        valstr = '  values:'+'%12.5f'%(parmDict[name])
3723        if name in sigDict:
3724            sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
3725        else:
3726            sigstr = '  esds  :'+12*' '
3727        out.append(namstr+'\n')
3728        out.append(valstr+'\n')
3729        out.append(sigstr+'\n')
3730    return out
3731
3732def WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
3733    '''Cribbed version of PrintRBObjTorAndSig but returns lists of strings.
3734    Moved so it can be used in ExportCIF
3735    '''
3736    out = []
3737    namstr = '  names :'
3738    valstr = '  values:'
3739    sigstr = '  esds  :'
3740    out.append(' Torsions:\n')
3741    for it in range(nTors):
3742        name = pfx+'RBRTr;'+str(it)+':'+rbsx
3743        namstr += '%12s'%('Tor'+str(it))
3744        valstr += '%12.4f'%(parmDict[name])
3745        if name in sigDict:
3746            sigstr += '%12.4f'%(sigDict[name])
3747    out.append(namstr+'\n')
3748    out.append(valstr+'\n')
3749    out.append(sigstr+'\n')
3750    return out
3751
3752def WriteResRBModel(RBModel):
3753    '''Write description of a residue rigid body. Code shifted from
3754    PrintResRBModel to make usable from G2export_CIF
3755    '''
3756    out = []
3757    atNames = RBModel['atNames']
3758    rbRef = RBModel['rbRef']
3759    rbSeq = RBModel['rbSeq']
3760    out.append('    At name       x          y          z\n')
3761    for name,xyz in zip(atNames,RBModel['rbXYZ']):
3762        out.append(%8s %10.4f %10.4f %10.4f\n'%(name,xyz[0],xyz[1],xyz[2]))
3763    out.append('  Orientation defined by: %s -> %s & %s -> %s\n'%
3764        (atNames[rbRef[0]],atNames[rbRef[1]],atNames[rbRef[0]],atNames[rbRef[2]]))
3765    if rbSeq:
3766        for i,rbseq in enumerate(rbSeq):
3767#                nameLst = [atNames[i] for i in rbseq[3]]
3768            out.append('  Torsion sequence %d Bond: %s - %s riding: %s\n'%
3769                    (i,atNames[rbseq[0]],atNames[rbseq[1]],str([atNames[i] for i in rbseq[3]])))
3770    return out
3771
3772def WriteVecRBModel(RBModel,sigDict={},irb=None):
3773    '''Write description of a vector rigid body. Code shifted from
3774    PrintVecRBModel to make usable from G2export_CIF
3775    '''
3776    out = []
3777    rbRef = RBModel['rbRef']
3778    atTypes = RBModel['rbTypes']
3779    for i in range(len(RBModel['VectMag'])):
3780        if sigDict and irb is not None:
3781            name = '::RBV;'+str(i)+':'+str(irb)
3782            s = G2mth.ValEsd(RBModel['VectMag'][i],sigDict.get(name,-.0001))
3783            out.append('Vector no.: %d Magnitude: %s\n'%(i,s))
3784        else:
3785            out.append('Vector no.: %d Magnitude: %8.4f Refine? %s\n'%
3786                           (i,RBModel['VectMag'][i],RBModel['VectRef'][i]))
3787        out.append('  No. Type     vx         vy         vz\n')
3788        for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
3789            out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(j,name,xyz[0],xyz[1],xyz[2]))
3790    out.append('  No. Type      x          y          z\n')
3791    for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
3792        out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(i,name,xyz[0],xyz[1],xyz[2]))
3793    return out
Note: See TracBrowser for help on using the repository browser.