source: trunk/GSASIIstrIO.py @ 3887

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

Fix restore of seq results after crash

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