source: trunk/GSASIIstrIO.py @ 4823

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

fix naming of new variable constraints; label them & improve constraint help text

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