source: trunk/GSASIIstrIO.py @ 3855

Last change on this file since 3855 was 3855, checked in by toby, 3 years ago

major speed up with very big sequential fits

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