source: trunk/GSASIIstrIO.py @ 3885

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

correct seq ref hist error

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