source: trunk/GSASIIstrIO.py @ 3873

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

trap error if unable to delete seq files

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