source: trunk/GSASIIstrIO.py @ 4401

Last change on this file since 4401 was 4401, checked in by vondreele, 19 months ago

Add texture index & esd to texture results - all 3 different cases. Shows in lst file or on console.

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