source: trunk/GSASIIstrIO.py @ 4197

Last change on this file since 4197 was 4197, checked in by vondreele, 23 months ago

put reference to Cromer & Liberman in Absorb & Fprime as well as refinement LST output & atmdata.py

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