source: trunk/GSASIIstrIO.py @ 4591

Last change on this file since 4591 was 4591, checked in by toby, 13 months ago

fix limit bug

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