source: trunk/GSASIIstrIO.py

Last change on this file was 5494, checked in by vondreele, 3 days ago

fixes to spin RB stuff - now will refine (no derivatives yet)

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