source: trunk/GSASIIstrIO.py @ 4629

Last change on this file since 4629 was 4629, checked in by vondreele, 12 months ago

GetXsectionCoeff? - remove a 'U' from an open - deprecated
close infile & outfile in SetUsedHistogramsAndPhases?

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