source: trunk/GSASIIstrIO.py @ 4666

Last change on this file since 4666 was 4666, checked in by vondreele, 11 months ago

Add output of resonant neutron scattering lengths for CW powder & single crystal histograms to .lst file

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