source: trunk/GSASIIstrIO.py @ 4399

Last change on this file since 4399 was 4399, checked in by toby, 19 months ago

major refactoring & corrections to ISODISTORT mode import and computations

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 164.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2020-04-13 01:46:57 +0000 (Mon, 13 Apr 2020) $
4# $Author: toby $
5# $Revision: 4399 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 4399 2020-04-13 01:46:57Z toby $
8########### SVN repository information ###################
9'''
10*GSASIIstrIO: structure I/O routines*
11-------------------------------------
12
13Contains routines for reading from GPX files and printing to the .LST file.
14Used for refinements and in G2scriptable.
15
16Should not contain any wxpython references as this should be able to be used
17in non-GUI settings.
18
19'''
20from __future__ import division, print_function
21import platform
22import os
23import os.path as ospath
24import time
25import math
26import copy
27if '2' in platform.python_version_tuple()[0]:
28    import cPickle
29else:
30    import pickle as cPickle
31import numpy as np
32import numpy.ma as ma
33import GSASIIpath
34GSASIIpath.SetVersionNumber("$Revision: 4399 $")
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        pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
2168        names = ['omega','chi','phi']
2169        namstr = '  names :'
2170        ptstr =  '  values:'
2171        sigstr = '  esds  :'
2172        for name in names:
2173            namstr += '%12s'%(name)
2174            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2175            if 'Sample '+name in SHtextureSig:
2176                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2177            else:
2178                sigstr += 12*' '
2179        pFile.write(namstr+'\n')
2180        pFile.write(ptstr+'\n')
2181        pFile.write(sigstr+'\n')
2182        pFile.write('\n Texture coefficients:\n')
2183        SHcoeff = textureData['SH Coeff'][1]
2184        SHkeys = list(SHcoeff.keys())
2185        nCoeff = len(SHcoeff)
2186        nBlock = nCoeff//10+1
2187        iBeg = 0
2188        iFin = min(iBeg+10,nCoeff)
2189        for block in range(nBlock):
2190            namstr = '  names :'
2191            ptstr =  '  values:'
2192            sigstr = '  esds  :'
2193            for name in SHkeys[iBeg:iFin]:
2194                namstr += '%12s'%(name)
2195                ptstr += '%12.3f'%(SHcoeff[name])
2196                if name in SHtextureSig:
2197                    sigstr += '%12.3f'%(SHtextureSig[name])
2198                else:
2199                    sigstr += 12*' '
2200            pFile.write(namstr+'\n')
2201            pFile.write(ptstr+'\n')
2202            pFile.write(sigstr+'\n')
2203            iBeg += 10
2204            iFin = min(iBeg+10,nCoeff)
2205           
2206    ##########################################################################
2207    # SetPhaseData starts here
2208    pFile.write('\n Phases:\n')
2209    for phase in Phases:
2210        pFile.write(' Result for phase: %s\n'%phase)
2211        pFile.write(135*'='+'\n')
2212        Phase = Phases[phase]
2213        General = Phase['General']
2214        SGData = General['SGData']
2215        Atoms = Phase['Atoms']
2216        AtLookup = []
2217        if Atoms and not General.get('doPawley'):
2218            cx,ct,cs,cia = General['AtomPtrs']
2219            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
2220        cell = General['Cell']
2221        pId = Phase['pId']
2222        pfx = str(pId)+'::'
2223        if cell[0]:
2224            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
2225            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
2226            pFile.write(' Reciprocal metric tensor: \n')
2227            ptfmt = "%15.9f"
2228            names = ['A11','A22','A33','A12','A13','A23']
2229            namstr = '  names :'
2230            ptstr =  '  values:'
2231            sigstr = '  esds  :'
2232            for name,a,siga in zip(names,A,sigA):
2233                namstr += '%15s'%(name)
2234                ptstr += ptfmt%(a)
2235                if siga:
2236                    sigstr += ptfmt%(siga)
2237                else:
2238                    sigstr += 15*' '
2239            pFile.write(namstr+'\n')
2240            pFile.write(ptstr+'\n')
2241            pFile.write(sigstr+'\n')
2242            cell[1:7] = G2lat.A2cell(A)
2243            cell[7] = G2lat.calc_V(A)
2244            pFile.write(' New unit cell:\n')
2245            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
2246            names = ['a','b','c','alpha','beta','gamma','Volume']
2247            namstr = '  names :'
2248            ptstr =  '  values:'
2249            sigstr = '  esds  :'
2250            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
2251                namstr += '%12s'%(name)
2252                ptstr += fmt%(a)
2253                if siga:
2254                    sigstr += fmt%(siga)
2255                else:
2256                    sigstr += 12*' '
2257            pFile.write(namstr+'\n')
2258            pFile.write(ptstr+'\n')
2259            pFile.write(sigstr+'\n')
2260        ik = 6  #for Pawley stuff below
2261        if General.get('Modulated',False):
2262            ik = 7
2263            Vec,vRef,maxH = General['SuperVec']
2264            if vRef:
2265                pFile.write(' New modulation vector:\n')
2266                namstr = '  names :'
2267                ptstr =  '  values:'
2268                sigstr = '  esds  :'
2269                for iv,var in enumerate(['mV0','mV1','mV2']):
2270                    namstr += '%12s'%(pfx+var)
2271                    ptstr += '%12.6f'%(parmDict[pfx+var])
2272                    if pfx+var in sigDict:
2273                        Vec[iv] = parmDict[pfx+var]
2274                        sigstr += '%12.6f'%(sigDict[pfx+var])
2275                    else:
2276                        sigstr += 12*' '
2277                pFile.write(namstr+'\n')
2278                pFile.write(ptstr+'\n')
2279                pFile.write(sigstr+'\n')
2280           
2281        General['Mass'] = 0.
2282        if Phase['General'].get('doPawley'):
2283            pawleyRef = Phase['Pawley ref']
2284            for i,refl in enumerate(pawleyRef):
2285                key = pfx+'PWLref:'+str(i)
2286                refl[ik] = parmDict[key]
2287                if key in sigDict:
2288                    refl[ik+1] = sigDict[key]
2289                else:
2290                    refl[ik+1] = 0
2291        else:
2292            VRBIds = RBIds['Vector']
2293            RRBIds = RBIds['Residue']
2294            RBModels = Phase['RBModels']
2295            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2296                jrb = VRBIds.index(RBObj['RBId'])
2297                rbsx = str(irb)+':'+str(jrb)
2298                pFile.write(' Vector rigid body parameters:\n')
2299                PrintRBObjPOAndSig('RBV',rbsx)
2300                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2301            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2302                jrb = RRBIds.index(RBObj['RBId'])
2303                rbsx = str(irb)+':'+str(jrb)
2304                pFile.write(' Residue rigid body parameters:\n')
2305                PrintRBObjPOAndSig('RBR',rbsx)
2306                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2307                PrintRBObjTorAndSig(rbsx)
2308            atomsSig = {}
2309            wavesSig = {}
2310            cx,ct,cs,cia = General['AtomPtrs']
2311            for i,at in enumerate(Atoms):
2312                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2313                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2314                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2315                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2316                for ind in range(cx,cx+4):
2317                    at[ind] = parmDict[names[ind]]
2318                    if ind in range(cx,cx+3):
2319                        name = names[ind].replace('A','dA')
2320                    else:
2321                        name = names[ind]
2322                    if name in sigDict:
2323                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2324                if at[cia] == 'I':
2325                    at[cia+1] = parmDict[names[cia+1]]
2326                    if names[cia+1] in sigDict:
2327                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2328                else:
2329                    for ind in range(cia+2,cia+8):
2330                        at[ind] = parmDict[names[ind]]
2331                        if names[ind] in sigDict:
2332                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2333                if General['Type'] == 'magnetic':
2334                    for ind in range(cx+4,cx+7):
2335                        at[ind] = parmDict[names[ind]]
2336                        if names[ind] in sigDict:
2337                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2338                ind = General['AtomTypes'].index(at[ct])
2339                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2340                if General.get('Modulated',False):
2341                    AtomSS = at[-1]['SS1']
2342                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2343                        Waves = AtomSS[Stype]
2344                        if len(Waves):
2345                            waveType = Waves[0]
2346                        else:
2347                            continue
2348                        for iw,wave in enumerate(Waves[1:]):
2349                            stiw = str(i)+':'+str(iw)
2350                            if Stype == 'Spos':
2351                                if waveType in ['ZigZag','Block',] and not iw:
2352                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2353                                else:
2354                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2355                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2356                            elif Stype == 'Sadp':
2357                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2358                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2359                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2360                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2361                            elif Stype == 'Sfrac':
2362                                if 'Crenel' in waveType and not iw:
2363                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2364                                else:
2365                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2366                            elif Stype == 'Smag':
2367                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2368                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2369                            for iname,name in enumerate(names):
2370                                AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2371                                if pfx+name in sigDict:
2372                                    wavesSig[name] = sigDict[pfx+name]
2373                   
2374            PrintAtomsAndSig(General,Atoms,atomsSig)
2375            if General['Type'] == 'magnetic':
2376                PrintMomentsAndSig(General,Atoms,atomsSig)
2377            if General.get('Modulated',False):
2378                PrintWavesAndSig(General,Atoms,wavesSig)
2379               
2380            density = G2mth.getDensity(General)[0]
2381            pFile.write('\n Density: %f.4 g/cm**3\n'%density)
2382           
2383       
2384        textureData = General['SH Texture']   
2385        if textureData['Order']:
2386            SHtextureSig = {}
2387            for name in ['omega','chi','phi']:
2388                aname = pfx+'SH '+name
2389                textureData['Sample '+name][1] = parmDict[aname]
2390                if aname in sigDict:
2391                    SHtextureSig['Sample '+name] = sigDict[aname]
2392            for name in textureData['SH Coeff'][1]:
2393                aname = pfx+name
2394                textureData['SH Coeff'][1][name] = parmDict[aname]
2395                if aname in sigDict:
2396                    SHtextureSig[name] = sigDict[aname]
2397            PrintSHtextureAndSig(textureData,SHtextureSig)
2398        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2399            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2400                textureData,RestraintDict[phase],pFile)
2401                   
2402################################################################################
2403##### Histogram & Phase data
2404################################################################################       
2405                   
2406def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2407    '''Loads the HAP histogram/phase information into dicts
2408
2409    :param dict Phases: phase information
2410    :param dict Histograms: Histogram information
2411    :param bool Print: prints information as it is read
2412    :param file pFile: file object to print to (the default, None causes printing to the console)
2413    :param bool resetRefList: Should the contents of the Reflection List be initialized
2414      on loading. The default, True, initializes the Reflection List as it is loaded.
2415
2416    :returns: (hapVary,hapDict,controlDict)
2417      * hapVary: list of refined variables
2418      * hapDict: dict with refined variables and their values
2419      * controlDict: dict with fixed parameters
2420    '''
2421   
2422    def PrintSize(hapData):
2423        if hapData[0] in ['isotropic','uniaxial']:
2424            line = '\n Size model    : %9s'%(hapData[0])
2425            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2426            if hapData[0] == 'uniaxial':
2427                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2428            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2429            pFile.write(line+'\n')
2430        else:
2431            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2432                (hapData[0],hapData[1][2],hapData[2][2]))
2433            Snames = ['S11','S22','S33','S12','S13','S23']
2434            ptlbls = ' names :'
2435            ptstr =  ' values:'
2436            varstr = ' refine:'
2437            for i,name in enumerate(Snames):
2438                ptlbls += '%12s' % (name)
2439                ptstr += '%12.3f' % (hapData[4][i])
2440                varstr += '%12s' % (str(hapData[5][i]))
2441            pFile.write(ptlbls+'\n')
2442            pFile.write(ptstr+'\n')
2443            pFile.write(varstr+'\n')
2444       
2445    def PrintMuStrain(hapData,SGData):
2446        if hapData[0] in ['isotropic','uniaxial']:
2447            line = '\n Mustrain model: %9s'%(hapData[0])
2448            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2449            if hapData[0] == 'uniaxial':
2450                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2451            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2452            pFile.write(line+'\n')
2453        else:
2454            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2455                (hapData[0],hapData[1][2],hapData[2][2]))
2456            Snames = G2spc.MustrainNames(SGData)
2457            ptlbls = ' names :'
2458            ptstr =  ' values:'
2459            varstr = ' refine:'
2460            for i,name in enumerate(Snames):
2461                ptlbls += '%12s' % (name)
2462                ptstr += '%12.1f' % (hapData[4][i])
2463                varstr += '%12s' % (str(hapData[5][i]))
2464            pFile.write(ptlbls+'\n')
2465            pFile.write(ptstr+'\n')
2466            pFile.write(varstr+'\n')
2467
2468    def PrintHStrain(hapData,SGData):
2469        pFile.write('\n Hydrostatic/elastic strain:\n')
2470        Hsnames = G2spc.HStrainNames(SGData)
2471        ptlbls = ' names :'
2472        ptstr =  ' values:'
2473        varstr = ' refine:'
2474        for i,name in enumerate(Hsnames):
2475            ptlbls += '%12s' % (name)
2476            ptstr += '%12.4g' % (hapData[0][i])
2477            varstr += '%12s' % (str(hapData[1][i]))
2478        pFile.write(ptlbls+'\n')
2479        pFile.write(ptstr+'\n')
2480        pFile.write(varstr+'\n')
2481
2482    def PrintSHPO(hapData):
2483        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2484        ptlbls = ' names :'
2485        ptstr =  ' values:'
2486        for item in hapData[5]:
2487            ptlbls += '%12s'%(item)
2488            ptstr += '%12.3f'%(hapData[5][item]) 
2489        pFile.write(ptlbls+'\n')
2490        pFile.write(ptstr+'\n')
2491   
2492    def PrintBabinet(hapData):
2493        pFile.write('\n Babinet form factor modification:\n')
2494        ptlbls = ' names :'
2495        ptstr =  ' values:'
2496        varstr = ' refine:'
2497        for name in ['BabA','BabU']:
2498            ptlbls += '%12s' % (name)
2499            ptstr += '%12.6f' % (hapData[name][0])
2500            varstr += '%12s' % (str(hapData[name][1]))
2501        pFile.write(ptlbls+'\n')
2502        pFile.write(ptstr+'\n')
2503        pFile.write(varstr+'\n')
2504       
2505    hapDict = {}
2506    hapVary = []
2507    controlDict = {}
2508   
2509    for phase in Phases:
2510        HistoPhase = Phases[phase]['Histograms']
2511        SGData = Phases[phase]['General']['SGData']
2512        cell = Phases[phase]['General']['Cell'][1:7]
2513        A = G2lat.cell2A(cell)
2514        if Phases[phase]['General'].get('Modulated',False):
2515            SSGData = Phases[phase]['General']['SSGData']
2516            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2517        pId = Phases[phase]['pId']
2518        for histogram in Histograms:
2519            if histogram not in HistoPhase and phase in Histograms[histogram]['Reflection Lists']:
2520                #remove previously created reflection list if histogram is removed from phase
2521                #print("removing ",phase,"from",histogram)
2522                del Histograms[histogram]['Reflection Lists'][phase]
2523        histoList = list(HistoPhase.keys())
2524        histoList.sort()
2525        for histogram in histoList:
2526            try:
2527                Histogram = Histograms[histogram]
2528            except KeyError:                       
2529                #skip if histogram not included e.g. in a sequential refinement
2530                continue
2531            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2532                if phase in Histograms[histogram]['Reflection Lists']:
2533                    del Histograms[histogram]['Reflection Lists'][phase]
2534                continue
2535            hapData = HistoPhase[histogram]
2536            hId = Histogram['hId']
2537            if 'PWDR' in histogram:
2538                limits = Histogram['Limits'][1]
2539                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2540                if 'C' in inst['Type'][1]:
2541                    try:
2542                        wave = inst['Lam'][1]
2543                    except KeyError:
2544                        wave = inst['Lam1'][1]
2545                    dmin = wave/(2.0*sind(limits[1]/2.0))
2546                elif 'T' in inst['Type'][0]:
2547                    dmin = limits[0]/inst['difC'][1]
2548                pfx = str(pId)+':'+str(hId)+':'
2549                if Phases[phase]['General']['doPawley']:
2550                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2551                    hapDict[pfx+'newLeBail'] = True
2552                    Tmin = G2lat.Dsp2pos(inst,dmin)
2553                    if 'C' in inst['Type'][1]:
2554                        limits[1] = min(limits[1],Tmin)
2555                    else:
2556                        limits[0] = max(limits[0],Tmin)
2557                else:
2558                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2559                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2560                if Phases[phase]['General']['Type'] == 'magnetic':
2561                    dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.))
2562                for item in ['Scale','Extinction']:
2563                    hapDict[pfx+item] = hapData[item][0]
2564                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2565                        hapVary.append(pfx+item)
2566                names = G2spc.HStrainNames(SGData)
2567                HSvals = []
2568                for i,name in enumerate(names):
2569                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2570                    HSvals.append(hapDict[pfx+name])
2571                    if hapData['HStrain'][1][i]:
2572                        hapVary.append(pfx+name)
2573                if 'Layer Disp' in hapData:
2574                    hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0]
2575                    if hapData['Layer Disp'][1]:
2576                        hapVary.append(pfx+'LayerDisp')
2577                else:
2578                    hapDict[pfx+'LayerDisp'] = 0.0
2579                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2580                if hapData['Pref.Ori.'][0] == 'MD':
2581                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2582                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2583                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2584                        hapVary.append(pfx+'MD')
2585                else:                           #'SH' spherical harmonics
2586                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2587                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2588                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2589                    controlDict[pfx+'SHhkl'] = []
2590                    try: #patch for old Pref.Ori. items
2591                        controlDict[pfx+'SHtoler'] = 0.1
2592                        if hapData['Pref.Ori.'][6][0] != '':
2593                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2594                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2595                    except IndexError:
2596                        pass
2597                    for item in hapData['Pref.Ori.'][5]:
2598                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2599                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2600                            hapVary.append(pfx+item)
2601                for item in ['Mustrain','Size']:
2602                    controlDict[pfx+item+'Type'] = hapData[item][0]
2603                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2604                    if hapData[item][2][2]:
2605                        hapVary.append(pfx+item+';mx')
2606                    if hapData[item][0] in ['isotropic','uniaxial']:
2607                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2608                        if hapData[item][2][0]:
2609                            hapVary.append(pfx+item+';i')
2610                        if hapData[item][0] == 'uniaxial':
2611                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2612                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2613                            if hapData[item][2][1]:
2614                                hapVary.append(pfx+item+';a')
2615                    else:       #generalized for mustrain or ellipsoidal for size
2616                        Nterms = len(hapData[item][4])
2617                        if item == 'Mustrain':
2618                            names = G2spc.MustrainNames(SGData)
2619                            pwrs = []
2620                            for name in names:
2621                                h,k,l = name[1:]
2622                                pwrs.append([int(h),int(k),int(l)])
2623                            controlDict[pfx+'MuPwrs'] = pwrs
2624                        for i in range(Nterms):
2625                            sfx = ';'+str(i)
2626                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2627                            if hapData[item][5][i]:
2628                                hapVary.append(pfx+item+sfx)
2629                if Phases[phase]['General']['Type'] != 'magnetic':
2630                    for bab in ['BabA','BabU']:
2631                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2632                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2633                            hapVary.append(pfx+bab)
2634                               
2635                if Print: 
2636                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2637                    pFile.write(135*'='+'\n')
2638                    if hapDict[pfx+'LeBail']:
2639                        pFile.write(' Perform LeBail extraction\n')                     
2640                    else:
2641                        pFile.write(' Phase fraction  : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2642                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2643                        if hapData['Pref.Ori.'][0] == 'MD':
2644                            Ax = hapData['Pref.Ori.'][3]
2645                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2646                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2647                        else: #'SH' for spherical harmonics
2648                            PrintSHPO(hapData['Pref.Ori.'])
2649                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2650                    PrintSize(hapData['Size'])
2651                    PrintMuStrain(hapData['Mustrain'],SGData)
2652                    PrintHStrain(hapData['HStrain'],SGData)
2653                    if 'Layer Disp' in hapData:
2654                        pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1]))
2655                    if Phases[phase]['General']['Type'] != 'magnetic':
2656                        if hapData['Babinet']['BabA'][0]:
2657                            PrintBabinet(hapData['Babinet'])
2658                if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2659                    hapData['newLeBail'] = True
2660                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2661                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2662                        hapData['newLeBail'] = False
2663                    refList = []
2664#                    Uniq = []
2665#                    Phi = []
2666                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2667                    if Phases[phase]['General'].get('Modulated',False):
2668                        ifSuper = True
2669                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2670                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2671                        for h,k,l,m,d in HKLd:
2672                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2673                            mul *= 2      # for powder overlap of Friedel pairs
2674                            if m or not ext or useExt:
2675                                if 'C' in inst['Type'][0]:
2676                                    pos = G2lat.Dsp2pos(inst,d)
2677                                    if limits[0] < pos < limits[1]:
2678                                        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])
2679                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2680#                                        Uniq.append(uniq)
2681#                                        Phi.append(phi)
2682                                elif 'T' in inst['Type'][0]:
2683                                    pos = G2lat.Dsp2pos(inst,d)
2684                                    if limits[0] < pos < limits[1]:
2685                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2686                                        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])
2687                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2688                                        #TODO - if tabulated put alp & bet in here
2689#                                        Uniq.append(uniq)
2690#                                        Phi.append(phi)
2691                    else:
2692                        ifSuper = False
2693                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2694                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2695                        for h,k,l,d in HKLd:
2696                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2697                            if ext and 'N' in inst['Type'][0] and  'MagSpGrp' in SGData:
2698                                ext = G2spc.checkMagextc([h,k,l],SGData)
2699                            mul *= 2      # for powder overlap of Friedel pairs
2700                            if ext and not useExt:
2701                                continue
2702                            if 'C' in inst['Type'][0]:
2703                                pos = G2lat.Dsp2pos(inst,d)
2704                                if limits[0] < pos < limits[1]:
2705                                    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])
2706                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2707#                                    Uniq.append(uniq)
2708#                                    Phi.append(phi)
2709                            elif 'T' in inst['Type'][0]:
2710                                pos = G2lat.Dsp2pos(inst,d)
2711                                if limits[0] < pos < limits[1]:
2712                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2713                                    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])
2714                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2715#                                    Uniq.append(uniq)
2716#                                    Phi.append(phi)
2717                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2718            elif 'HKLF' in histogram:
2719                inst = Histogram['Instrument Parameters'][0]
2720                hId = Histogram['hId']
2721                hfx = ':%d:'%(hId)
2722                for item in inst:
2723                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2724                    hapDict[hfx+item] = inst[item][1]
2725                pfx = str(pId)+':'+str(hId)+':'
2726                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2727                if hapData['Scale'][1]:
2728                    hapVary.append(pfx+'Scale')
2729                               
2730                extApprox,extType,extParms = hapData['Extinction']
2731                controlDict[pfx+'EType'] = extType
2732                controlDict[pfx+'EApprox'] = extApprox
2733                if 'C' in inst['Type'][0]:
2734                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2735                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2736                if 'Primary' in extType:
2737                    Ekey = ['Ep',]
2738                elif 'I & II' in extType:
2739                    Ekey = ['Eg','Es']
2740                elif 'Secondary Type II' == extType:
2741                    Ekey = ['Es',]
2742                elif 'Secondary Type I' == extType:
2743                    Ekey = ['Eg',]
2744                else:   #'None'
2745                    Ekey = []
2746                for eKey in Ekey:
2747                    hapDict[pfx+eKey] = extParms[eKey][0]
2748                    if extParms[eKey][1]:
2749                        hapVary.append(pfx+eKey)
2750                for bab in ['BabA','BabU']:
2751                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2752                    if hapData['Babinet'][bab][1]:
2753                        hapVary.append(pfx+bab)
2754                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2755                if len(Twins) == 1:
2756                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2757                    if hapData.get('Flack',[0,False])[1]:
2758                        hapVary.append(pfx+'Flack')
2759                sumTwFr = 0.
2760                controlDict[pfx+'TwinLaw'] = []
2761                controlDict[pfx+'TwinInv'] = []
2762                NTL = 0           
2763                for it,twin in enumerate(Twins):
2764                    if 'bool' in str(type(twin[0])):
2765                        controlDict[pfx+'TwinInv'].append(twin[0])
2766                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2767                    else:
2768                        NTL += 1
2769                        controlDict[pfx+'TwinInv'].append(False)
2770                        controlDict[pfx+'TwinLaw'].append(twin[0])
2771                    if it:
2772                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2773                        sumTwFr += twin[1]
2774                    else:
2775                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2776                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2777                    if Twins[0][1][1]:
2778                        hapVary.append(pfx+'TwinFr:'+str(it))
2779                controlDict[pfx+'NTL'] = NTL
2780                #need to make constraint on TwinFr
2781                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2782                if len(Twins) > 1:    #force sum to unity
2783                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2784                if Print: 
2785                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2786                    pFile.write(135*'='+'\n')
2787                    pFile.write(' Scale factor     : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2788                    if extType != 'None':
2789                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2790                        text = ' Parameters       :'
2791                        for eKey in Ekey:
2792                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2793                        pFile.write(text+'\n')
2794                    if hapData['Babinet']['BabA'][0]:
2795                        PrintBabinet(hapData['Babinet'])
2796                    if not SGData['SGInv'] and len(Twins) == 1:
2797                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2798                    if len(Twins) > 1:
2799                        for it,twin in enumerate(Twins):
2800                            if 'bool' in str(type(twin[0])):
2801                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2802                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2803                            else:
2804                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2805                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2806                       
2807                Histogram['Reflection Lists'] = phase       
2808               
2809    return hapVary,hapDict,controlDict
2810   
2811def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2812    'needs a doc string'
2813   
2814    def PrintSizeAndSig(hapData,sizeSig):
2815        line = '\n Size model:     %9s'%(hapData[0])
2816        refine = False
2817        if hapData[0] in ['isotropic','uniaxial']:
2818            line += ' equatorial:%12.5f'%(hapData[1][0])
2819            if sizeSig[0][0]:
2820                line += ', sig:%8.4f'%(sizeSig[0][0])
2821                refine = True
2822            if hapData[0] == 'uniaxial':
2823                line += ' axial:%12.4f'%(hapData[1][1])
2824                if sizeSig[0][1]:
2825                    refine = True
2826                    line += ', sig:%8.4f'%(sizeSig[0][1])
2827            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2828            if sizeSig[0][2]:
2829                refine = True
2830                line += ', sig:%8.4f'%(sizeSig[0][2])
2831            if refine:
2832                pFile.write(line+'\n')
2833        else:
2834            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2835            if sizeSig[0][2]:
2836                refine = True
2837                line += ', sig:%8.4f'%(sizeSig[0][2])
2838            Snames = ['S11','S22','S33','S12','S13','S23']
2839            ptlbls = ' name  :'
2840            ptstr =  ' value :'
2841            sigstr = ' sig   :'
2842            for i,name in enumerate(Snames):
2843                ptlbls += '%12s' % (name)
2844                ptstr += '%12.6f' % (hapData[4][i])
2845                if sizeSig[1][i]:
2846                    refine = True
2847                    sigstr += '%12.6f' % (sizeSig[1][i])
2848                else:
2849                    sigstr += 12*' '
2850            if refine:
2851                pFile.write(line+'\n')
2852                pFile.write(ptlbls+'\n')
2853                pFile.write(ptstr+'\n')
2854                pFile.write(sigstr+'\n')
2855       
2856    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2857        line = '\n Mustrain model: %9s\n'%(hapData[0])
2858        refine = False
2859        if hapData[0] in ['isotropic','uniaxial']:
2860            line += ' equatorial:%12.1f'%(hapData[1][0])
2861            if mustrainSig[0][0]:
2862                line += ', sig:%8.1f'%(mustrainSig[0][0])
2863                refine = True
2864            if hapData[0] == 'uniaxial':
2865                line += ' axial:%12.1f'%(hapData[1][1])
2866                if mustrainSig[0][1]:
2867                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2868            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2869            if mustrainSig[0][2]:
2870                refine = True
2871                line += ', sig:%8.3f'%(mustrainSig[0][2])
2872            if refine:
2873                pFile.write(line+'\n')
2874        else:
2875            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2876            if mustrainSig[0][2]:
2877                refine = True
2878                line += ', sig:%8.3f'%(mustrainSig[0][2])
2879            Snames = G2spc.MustrainNames(SGData)
2880            ptlbls = ' name  :'
2881            ptstr =  ' value :'
2882            sigstr = ' sig   :'
2883            for i,name in enumerate(Snames):
2884                ptlbls += '%12s' % (name)
2885                ptstr += '%12.1f' % (hapData[4][i])
2886                if mustrainSig[1][i]:
2887                    refine = True
2888                    sigstr += '%12.1f' % (mustrainSig[1][i])
2889                else:
2890                    sigstr += 12*' '
2891            if refine:
2892                pFile.write(line+'\n')
2893                pFile.write(ptlbls+'\n')
2894                pFile.write(ptstr+'\n')
2895                pFile.write(sigstr+'\n')
2896           
2897    def PrintHStrainAndSig(hapData,strainSig,SGData):
2898        Hsnames = G2spc.HStrainNames(SGData)
2899        ptlbls = ' name  :'
2900        ptstr =  ' value :'
2901        sigstr = ' sig   :'
2902        refine = False
2903        for i,name in enumerate(Hsnames):
2904            ptlbls += '%12s' % (name)
2905            ptstr += '%12.4g' % (hapData[0][i])
2906            if name in strainSig:
2907                refine = True
2908                sigstr += '%12.4g' % (strainSig[name])
2909            else:
2910                sigstr += 12*' '
2911        if refine:
2912            pFile.write('\n Hydrostatic/elastic strain:\n')
2913            pFile.write(ptlbls+'\n')
2914            pFile.write(ptstr+'\n')
2915            pFile.write(sigstr+'\n')
2916       
2917    def PrintSHPOAndSig(pfx,hapData,POsig):
2918        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2919        ptlbls = ' names :'
2920        ptstr =  ' values:'
2921        sigstr = ' sig   :'
2922        for item in hapData[5]:
2923            ptlbls += '%12s'%(item)
2924            ptstr += '%12.3f'%(hapData[5][item])
2925            if pfx+item in POsig:
2926                sigstr += '%12.3f'%(POsig[pfx+item])
2927            else:
2928                sigstr += 12*' ' 
2929        pFile.write(ptlbls+'\n')
2930        pFile.write(ptstr+'\n')
2931        pFile.write(sigstr+'\n')
2932       
2933    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2934        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2935        text = ''
2936        for item in hapData[2]:
2937            if pfx+item in ScalExtSig:
2938                text += '       %s: '%(item)
2939                text += '%12.2e'%(hapData[2][item][0])
2940                if pfx+item in ScalExtSig:
2941                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2942        pFile.write(text+'\n')   
2943       
2944    def PrintBabinetAndSig(pfx,hapData,BabSig):
2945        pFile.write('\n Babinet form factor modification:\n')
2946        ptlbls = ' names :'
2947        ptstr =  ' values:'
2948        sigstr = ' sig   :'
2949        for item in hapData:
2950            ptlbls += '%12s'%(item)
2951            ptstr += '%12.3f'%(hapData[item][0])
2952            if pfx+item in BabSig:
2953                sigstr += '%12.3f'%(BabSig[pfx+item])
2954            else:
2955                sigstr += 12*' ' 
2956        pFile.write(ptlbls+'\n')
2957        pFile.write(ptstr+'\n')
2958        pFile.write(sigstr+'\n')
2959       
2960    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2961        pFile.write('\n Twin Law fractions :\n')
2962        ptlbls = ' names :'
2963        ptstr =  ' values:'
2964        sigstr = ' sig   :'
2965        for it,item in enumerate(twinData):
2966            ptlbls += '     twin #%d'%(it)
2967            if it:
2968                ptstr += '%12.3f'%(item[1])
2969            else:
2970                ptstr += '%12.3f'%(item[1][0])
2971            if pfx+'TwinFr:'+str(it) in TwinSig:
2972                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2973            else:
2974                sigstr += 12*' ' 
2975        pFile.write(ptlbls+'\n')
2976        pFile.write(ptstr+'\n')
2977        pFile.write(sigstr+'\n')
2978       
2979   
2980    PhFrExtPOSig = {}
2981    SizeMuStrSig = {}
2982    ScalExtSig = {}
2983    BabSig = {}
2984    TwinFrSig = {}
2985    wtFrSum = {}
2986    for phase in Phases:
2987        HistoPhase = Phases[phase]['Histograms']
2988        General = Phases[phase]['General']
2989        SGData = General['SGData']
2990        pId = Phases[phase]['pId']
2991        histoList = list(HistoPhase.keys())
2992        histoList.sort()
2993        for histogram in histoList:
2994            try:
2995                Histogram = Histograms[histogram]
2996            except KeyError:                       
2997                #skip if histogram not included e.g. in a sequential refinement
2998                continue
2999            if not Phases[phase]['Histograms'][histogram]['Use']:
3000                #skip if phase absent from this histogram
3001                continue
3002            hapData = HistoPhase[histogram]
3003            hId = Histogram['hId']
3004            pfx = str(pId)+':'+str(hId)+':'
3005            if hId not in wtFrSum:
3006                wtFrSum[hId] = 0.
3007            if 'PWDR' in histogram:
3008                for item in ['Scale','Extinction']:
3009                    hapData[item][0] = parmDict[pfx+item]
3010                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
3011                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
3012                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
3013                if hapData['Pref.Ori.'][0] == 'MD':
3014                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
3015                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
3016                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
3017                else:                           #'SH' spherical harmonics
3018                    for item in hapData['Pref.Ori.'][5]:
3019                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
3020                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
3021                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
3022                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
3023                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
3024                    pfx+'HStrain':{}})                 
3025                for item in ['Mustrain','Size']:
3026                    hapData[item][1][2] = parmDict[pfx+item+';mx']
3027#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
3028                    if pfx+item+';mx' in sigDict:
3029                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
3030                    if hapData[item][0] in ['isotropic','uniaxial']:                   
3031                        hapData[item][1][0] = parmDict[pfx+item+';i']
3032                        if item == 'Size':
3033                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
3034                        if pfx+item+';i' in sigDict: 
3035                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
3036                        if hapData[item][0] == 'uniaxial':
3037                            hapData[item][1][1] = parmDict[pfx+item+';a']
3038                            if item == 'Size':
3039                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
3040                            if pfx+item+';a' in sigDict:
3041                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
3042                    else:       #generalized for mustrain or ellipsoidal for size
3043                        Nterms = len(hapData[item][4])
3044                        for i in range(Nterms):
3045                            sfx = ';'+str(i)
3046                            hapData[item][4][i] = parmDict[pfx+item+sfx]
3047                            if pfx+item+sfx in sigDict:
3048                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
3049                names = G2spc.HStrainNames(SGData)
3050                for i,name in enumerate(names):
3051                    hapData['HStrain'][0][i] = parmDict[pfx+name]
3052                    if pfx+name in sigDict:
3053                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
3054                if 'Layer Disp' in hapData:
3055                    hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp']
3056                    if pfx+'LayerDisp' in sigDict:
3057                        SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp']
3058                if Phases[phase]['General']['Type'] != 'magnetic':
3059                    for name in ['BabA','BabU']:
3060                        hapData['Babinet'][name][0] = parmDict[pfx+name]
3061                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
3062                            BabSig[pfx+name] = sigDict[pfx+name]               
3063               
3064            elif 'HKLF' in histogram:
3065                for item in ['Scale','Flack']:
3066                    if parmDict.get(pfx+item):
3067                        hapData[item][0] = parmDict[pfx+item]
3068                        if pfx+item in sigDict:
3069                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3070                for item in ['Ep','Eg','Es']:
3071                    if parmDict.get(pfx+item):
3072                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
3073                        if pfx+item in sigDict:
3074                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3075                for item in ['BabA','BabU']:
3076                    hapData['Babinet'][item][0] = parmDict[pfx+item]
3077                    if pfx+item in sigDict:
3078                        BabSig[pfx+item] = sigDict[pfx+item]
3079                if 'Twins' in hapData:
3080                    it = 1
3081                    sumTwFr = 0.
3082                    while True:
3083                        try:
3084                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
3085                            if pfx+'TwinFr:'+str(it) in sigDict:
3086                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
3087                            if it:
3088                                sumTwFr += hapData['Twins'][it][1]
3089                            it += 1
3090                        except KeyError:
3091                            break
3092                    hapData['Twins'][0][1][0] = 1.-sumTwFr
3093
3094    if Print:
3095        for phase in Phases:
3096            HistoPhase = Phases[phase]['Histograms']
3097            General = Phases[phase]['General']
3098            SGData = General['SGData']
3099            pId = Phases[phase]['pId']
3100            histoList = list(HistoPhase.keys())
3101            histoList.sort()
3102            for histogram in histoList:
3103                try:
3104                    Histogram = Histograms[histogram]
3105                except KeyError:                       
3106                    #skip if histogram not included e.g. in a sequential refinement
3107                    continue
3108                hapData = HistoPhase[histogram]
3109                hId = Histogram['hId']
3110                Histogram['Residuals'][str(pId)+'::Name'] = phase
3111                pfx = str(pId)+':'+str(hId)+':'
3112                hfx = ':%s:'%(hId)
3113                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
3114                    continue
3115                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
3116                pFile.write(135*'='+'\n')
3117                if 'PWDR' in histogram:
3118                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
3119                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
3120                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
3121                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
3122                   
3123                    if parmDict[pfx+'LeBail']:
3124                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
3125                    else:
3126                        if pfx+'Scale' in PhFrExtPOSig:
3127                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
3128                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
3129                            pFile.write(' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f\n'%
3130                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
3131                        if pfx+'Extinction' in PhFrExtPOSig:
3132                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
3133                        if hapData['Pref.Ori.'][0] == 'MD':
3134                            if pfx+'MD' in PhFrExtPOSig:
3135                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
3136                        else:
3137                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
3138                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
3139                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
3140                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
3141                    if pfx+'LayerDisp' in SizeMuStrSig:
3142                        pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp']))           
3143                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
3144                        if len(BabSig):
3145                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3146                   
3147                elif 'HKLF' in histogram:
3148                    Inst = Histogram['Instrument Parameters'][0]
3149                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
3150                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
3151                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
3152                    if FFtables != None and 'N' not in Inst['Type'][0]:
3153                        PrintFprime(FFtables,hfx,pFile)
3154                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
3155                    if pfx+'Scale' in ScalExtSig:
3156                        pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
3157                    if hapData['Extinction'][0] != 'None':
3158                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
3159                    if len(BabSig):
3160                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3161                    if pfx+'Flack' in ScalExtSig:
3162                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
3163                    if len(TwinFrSig):
3164                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
3165
3166################################################################################
3167##### Histogram data
3168################################################################################       
3169                   
3170def GetHistogramData(Histograms,Print=True,pFile=None):
3171    'needs a doc string'
3172   
3173    def GetBackgroundParms(hId,Background):
3174        Back = Background[0]
3175        DebyePeaks = Background[1]
3176        bakType,bakFlag = Back[:2]
3177        backVals = Back[3:]
3178        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
3179        backDict = dict(zip(backNames,backVals))
3180        backVary = []
3181        if bakFlag:
3182            backVary = backNames
3183        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
3184        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
3185        debyeDict = {}
3186        debyeList = []
3187        for i in range(DebyePeaks['nDebye']):
3188            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
3189            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
3190            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
3191        debyeVary = []
3192        for item in debyeList:
3193            if item[1]:
3194                debyeVary.append(item[0])
3195        backDict.update(debyeDict)
3196        backVary += debyeVary
3197        peakDict = {}
3198        peakList = []
3199        for i in range(DebyePeaks['nPeaks']):
3200            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
3201                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
3202            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
3203            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
3204        peakVary = []
3205        for item in peakList:
3206            if item[1]:
3207                peakVary.append(item[0])
3208        backDict.update(peakDict)
3209        backVary += peakVary
3210        return bakType,backDict,backVary           
3211       
3212    def GetInstParms(hId,Inst):
3213        #patch
3214        if 'Z' not in Inst:
3215            Inst['Z'] = [0.0,0.0,False]
3216        dataType = Inst['Type'][0]
3217        instDict = {}
3218        insVary = []
3219        pfx = ':'+str(hId)+':'
3220        insKeys = list(Inst.keys())
3221        insKeys.sort()
3222        for item in insKeys:
3223            insName = pfx+item
3224            instDict[insName] = Inst[item][1]
3225            if len(Inst[item]) > 2 and Inst[item][2]:
3226                insVary.append(insName)
3227        if 'C' in dataType:
3228            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
3229        elif 'T' in dataType:   #trap zero alp, bet coeff.
3230            if not instDict[pfx+'alpha']:
3231                instDict[pfx+'alpha'] = 1.0
3232            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3233                instDict[pfx+'beta-1'] = 1.0
3234        return dataType,instDict,insVary
3235       
3236    def GetSampleParms(hId,Sample):
3237        sampVary = []
3238        hfx = ':'+str(hId)+':'       
3239        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
3240            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
3241        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
3242            if key in Sample:
3243                sampDict[hfx+key] = Sample[key]
3244        Type = Sample['Type']
3245        if 'Bragg' in Type:             #Bragg-Brentano
3246            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3247                sampDict[hfx+item] = Sample[item][0]
3248                if Sample[item][1]:
3249                    sampVary.append(hfx+item)
3250        elif 'Debye' in Type:        #Debye-Scherrer
3251            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3252                sampDict[hfx+item] = Sample[item][0]
3253                if Sample[item][1]:
3254                    sampVary.append(hfx+item)
3255        return Type,sampDict,sampVary
3256       
3257    def PrintBackground(Background):
3258        Back = Background[0]
3259        DebyePeaks = Background[1]
3260        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
3261        line = ' Coefficients: '
3262        for i,back in enumerate(Back[3:]):
3263            line += '%10.3f'%(back)
3264            if i and not i%10:
3265                line += '\n'+15*' '
3266        pFile.write(line+'\n')
3267        if DebyePeaks['nDebye']:
3268            pFile.write('\n Debye diffuse scattering coefficients\n')
3269            parms = ['DebyeA','DebyeR','DebyeU']
3270            line = ' names :  '
3271            for parm in parms:
3272                line += '%8s refine?'%(parm)
3273            pFile.write(line+'\n')
3274            for j,term in enumerate(DebyePeaks['debyeTerms']):
3275                line = ' term'+'%2d'%(j)+':'
3276                for i in range(3):
3277                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3278                pFile.write(line+'\n')
3279        if DebyePeaks['nPeaks']:
3280            pFile.write('\n Single peak coefficients\n')
3281            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
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['peaksList']):
3287                line = ' peak'+'%2d'%(j)+':'
3288                for i in range(4):
3289                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3290                pFile.write(line+'\n')
3291       
3292    def PrintInstParms(Inst):
3293        pFile.write('\n Instrument Parameters:\n')
3294        insKeys = list(Inst.keys())
3295        insKeys.sort()
3296        iBeg = 0
3297        Ok = True
3298        while Ok:
3299            ptlbls = ' name  :'
3300            ptstr =  ' value :'
3301            varstr = ' refine:'
3302            iFin = min(iBeg+9,len(insKeys))
3303            for item in insKeys[iBeg:iFin]:
3304                if item not in ['Type','Source']:
3305                    ptlbls += '%12s' % (item)
3306                    ptstr += '%12.6f' % (Inst[item][1])
3307                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3308                        varstr += 12*' '
3309                    else:
3310                        varstr += '%12s' % (str(bool(Inst[item][2])))
3311            pFile.write(ptlbls+'\n')
3312            pFile.write(ptstr+'\n')
3313            pFile.write(varstr+'\n')
3314            iBeg = iFin
3315            if iBeg == len(insKeys):
3316                Ok = False
3317            else:
3318                pFile.write('\n')
3319       
3320    def PrintSampleParms(Sample):
3321        pFile.write('\n Sample Parameters:\n')
3322        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3323            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3324        ptlbls = ' name  :'
3325        ptstr =  ' value :'
3326        varstr = ' refine:'
3327        if 'Bragg' in Sample['Type']:
3328            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3329                ptlbls += '%14s'%(item)
3330                ptstr += '%14.4f'%(Sample[item][0])
3331                varstr += '%14s'%(str(bool(Sample[item][1])))
3332           
3333        elif 'Debye' in Type:        #Debye-Scherrer
3334            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3335                ptlbls += '%14s'%(item)
3336                ptstr += '%14.4f'%(Sample[item][0])
3337                varstr += '%14s'%(str(bool(Sample[item][1])))
3338
3339        pFile.write(ptlbls+'\n')
3340        pFile.write(ptstr+'\n')
3341        pFile.write(varstr+'\n')
3342       
3343    histDict = {}
3344    histVary = []
3345    controlDict = {}
3346    histoList = list(Histograms.keys())
3347    histoList.sort()
3348    for histogram in histoList:
3349        if 'PWDR' in histogram:
3350            Histogram = Histograms[histogram]
3351            hId = Histogram['hId']
3352            pfx = ':'+str(hId)+':'
3353            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3354            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3355            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3356            for excl in controlDict[pfx+'Exclude']:
3357                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3358            if controlDict[pfx+'Exclude']:
3359                ma.mask_rows(Histogram['Data'])
3360            Background = Histogram['Background']
3361            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3362            controlDict[pfx+'bakType'] = Type
3363            histDict.update(bakDict)
3364            histVary += bakVary
3365           
3366            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3367            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3368                G2fil.G2Print ('Warning: tabulated profile coefficients are ignored')
3369            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3370            controlDict[pfx+'histType'] = Type
3371            if 'XC' in Type:
3372                if pfx+'Lam1' in instDict:
3373                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3374                else:
3375                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3376            histDict.update(instDict)
3377            histVary += insVary
3378           
3379            Sample = Histogram['Sample Parameters']
3380            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3381            controlDict[pfx+'instType'] = Type
3382            histDict.update(sampDict)
3383            histVary += sampVary
3384           
3385   
3386            if Print: 
3387                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3388                pFile.write(135*'='+'\n')
3389                Units = {'C':' deg','T':' msec'}
3390                units = Units[controlDict[pfx+'histType'][2]]
3391                Limits = controlDict[pfx+'Limits']
3392                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3393                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3394                if len(controlDict[pfx+'Exclude']):
3395                    excls = controlDict[pfx+'Exclude']
3396                    for excl in excls:
3397                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3398                PrintSampleParms(Sample)
3399                PrintInstParms(Inst[0])
3400                PrintBackground(Background)
3401        elif 'HKLF' in histogram:
3402            Histogram = Histograms[histogram]
3403            hId = Histogram['hId']
3404            pfx = ':'+str(hId)+':'
3405            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3406            Inst = Histogram['Instrument Parameters'][0]
3407            controlDict[pfx+'histType'] = Inst['Type'][0]
3408            if 'X' in Inst['Type'][0]:
3409                histDict[pfx+'Lam'] = Inst['Lam'][1]
3410                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3411            elif 'NC' in Inst['Type'][0]:                   
3412                histDict[pfx+'Lam'] = Inst['Lam'][1]
3413    return histVary,histDict,controlDict
3414   
3415def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3416    'needs a doc string'
3417   
3418    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3419        Back = Background[0]
3420        DebyePeaks = Background[1]
3421        lenBack = len(Back[3:])
3422        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3423        for i in range(lenBack):
3424            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3425            if pfx+'Back;'+str(i) in sigDict:
3426                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3427        if DebyePeaks['nDebye']:
3428            for i in range(DebyePeaks['nDebye']):
3429                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3430                for j,name in enumerate(names):
3431                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3432                    if name in sigDict:
3433                        backSig[lenBack+3*i+j] = sigDict[name]           
3434        if DebyePeaks['nPeaks']:
3435            for i in range(DebyePeaks['nPeaks']):
3436                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3437                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3438                for j,name in enumerate(names):
3439                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3440                    if name in sigDict:
3441                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3442        return backSig
3443       
3444    def SetInstParms(pfx,Inst,parmDict,sigDict):
3445        instSig = {}
3446        insKeys = list(Inst.keys())
3447        insKeys.sort()
3448        for item in insKeys:
3449            insName = pfx+item
3450            Inst[item][1] = parmDict[insName]
3451            if insName in sigDict:
3452                instSig[item] = sigDict[insName]
3453            else:
3454                instSig[item] = 0
3455        return instSig
3456       
3457    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3458        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3459            sampSig = [0 for i in range(5)]
3460            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3461                Sample[item][0] = parmDict[pfx+item]
3462                if pfx+item in sigDict:
3463                    sampSig[i] = sigDict[pfx+item]
3464        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3465            sampSig = [0 for i in range(4)]
3466            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3467                Sample[item][0] = parmDict[pfx+item]
3468                if pfx+item in sigDict:
3469                    sampSig[i] = sigDict[pfx+item]
3470        return sampSig
3471       
3472    def PrintBackgroundSig(Background,backSig):
3473        Back = Background[0]
3474        DebyePeaks = Background[1]
3475        valstr = ' value : '
3476        sigstr = ' sig   : '
3477        refine = False
3478        for i,back in enumerate(Back[3:]):
3479            valstr += '%10.4g'%(back)
3480            if Back[1]:
3481                refine = True
3482                sigstr += '%10.4g'%(backSig[i])
3483            else:
3484                sigstr += 10*' '
3485        if refine:
3486            pFile.write('\n Background function: %s\n'%Back[0])
3487            pFile.write(valstr+'\n')
3488            pFile.write(sigstr+'\n')
3489        if DebyePeaks['nDebye']:
3490            ifAny = False
3491            ptfmt = "%12.3f"
3492            names =  ' names :'
3493            ptstr =  ' values:'
3494            sigstr = ' esds  :'
3495            for item in sigDict:
3496                if 'Debye' in item:
3497                    ifAny = True
3498                    names += '%12s'%(item)
3499                    ptstr += ptfmt%(parmDict[item])
3500                    sigstr += ptfmt%(sigDict[item])
3501            if ifAny:
3502                pFile.write('\n Debye diffuse scattering coefficients\n')
3503                pFile.write(names+'\n')
3504                pFile.write(ptstr+'\n')
3505                pFile.write(sigstr+'\n')
3506        if DebyePeaks['nPeaks']:
3507            pFile.write('\n Single peak coefficients:\n')
3508            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3509            line = ' peak no. '
3510            for parm in parms:
3511                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3512            pFile.write(line+'\n')
3513            for ip in range(DebyePeaks['nPeaks']):
3514                ptstr = ' %4d '%(ip)
3515                for parm in parms:
3516                    fmt = '%14.3f'
3517                    efmt = '%12.3f'
3518                    if parm == 'BkPkpos':
3519                        fmt = '%14.4f'
3520                        efmt = '%12.4f'
3521                    name = pfx+parm+';%d'%(ip)
3522                    ptstr += fmt%(parmDict[name])
3523                    if name in sigDict:
3524                        ptstr += efmt%(sigDict[name])
3525                    else:
3526                        ptstr += 12*' '
3527                pFile.write(ptstr+'\n')
3528        sumBk = np.array(Histogram['sumBk'])
3529        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3530            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3531       
3532    def PrintInstParmsSig(Inst,instSig):
3533        refine = False
3534        insKeys = list(instSig.keys())
3535        insKeys.sort()
3536        iBeg = 0
3537        Ok = True
3538        while Ok:
3539            ptlbls = ' names :'
3540            ptstr =  ' value :'
3541            sigstr = ' sig   :'
3542            iFin = min(iBeg+9,len(insKeys))
3543            for name in insKeys[iBeg:iFin]:
3544                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3545                    ptlbls += '%12s' % (name)
3546                    ptstr += '%12.6f' % (Inst[name][1])
3547                    if instSig[name]:
3548                        refine = True
3549                        sigstr += '%12.6f' % (instSig[name])
3550                    else:
3551                        sigstr += 12*' '
3552            if refine:
3553                pFile.write('\n Instrument Parameters:\n')
3554                pFile.write(ptlbls+'\n')
3555                pFile.write(ptstr+'\n')
3556                pFile.write(sigstr+'\n')
3557            iBeg = iFin
3558            if iBeg == len(insKeys):
3559                Ok = False
3560       
3561    def PrintSampleParmsSig(Sample,sampleSig):
3562        ptlbls = ' names :'
3563        ptstr =  ' values:'
3564        sigstr = ' sig   :'
3565        refine = False
3566        if 'Bragg' in Sample['Type']:
3567            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3568                ptlbls += '%14s'%(item)
3569                ptstr += '%14.4f'%(Sample[item][0])
3570                if sampleSig[i]:
3571                    refine = True
3572                    sigstr += '%14.4f'%(sampleSig[i])
3573                else:
3574                    sigstr += 14*' '
3575           
3576        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3577            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3578                ptlbls += '%14s'%(item)
3579                ptstr += '%14.4f'%(Sample[item][0])
3580                if sampleSig[i]:
3581                    refine = True
3582                    sigstr += '%14.4f'%(sampleSig[i])
3583                else:
3584                    sigstr += 14*' '
3585
3586        if refine:
3587            pFile.write('\n Sample Parameters:\n')
3588            pFile.write(ptlbls+'\n')
3589            pFile.write(ptstr+'\n')
3590            pFile.write(sigstr+'\n')
3591       
3592    histoList = list(Histograms.keys())
3593    histoList.sort()
3594    for histogram in histoList:
3595        if 'PWDR' in histogram:
3596            Histogram = Histograms[histogram]
3597            hId = Histogram['hId']
3598            pfx = ':'+str(hId)+':'
3599            Background = Histogram['Background']
3600            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3601           
3602            Inst = Histogram['Instrument Parameters'][0]
3603            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3604       
3605            Sample = Histogram['Sample Parameters']
3606            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3607
3608            pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3609            pFile.write(135*'='+'\n')
3610            pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3611            pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3612                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3613            pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3614                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3615            if Print:
3616                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3617                if FFtables != None and 'N' not in Inst['Type'][0]:
3618                    PrintFprime(FFtables,pfx,pFile)
3619                PrintSampleParmsSig(Sample,sampSig)
3620                PrintInstParmsSig(Inst,instSig)
3621                PrintBackgroundSig(Background,backSig)
3622               
Note: See TracBrowser for help on using the repository browser.