source: trunk/GSASIIstrIO.py @ 3883

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

correct seq ref hist error

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