source: trunk/GSASIIstrIO.py @ 4534

Last change on this file since 4534 was 4534, checked in by toby, 15 months ago

implement variable limits; show cell under Dij vals

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