source: trunk/GSASIIstrIO.py @ 4588

Last change on this file since 4588 was 4588, checked in by toby, 2 years ago

use G2VarObj for param limits; add more info to seq. ref. done dialog; show Frozen in show LS parameters

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