source: trunk/GSASIIstrIO.py @ 4370

Last change on this file since 4370 was 4370, checked in by vondreele, 20 months ago

Add new parameter to Data tab - Layer DispSizerThis? is intended to cover multilayered materials of different phases with layers perpendicular to beam - it might have other uses.

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 161.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2020-03-15 07:57:46 +0000 (Sun, 15 Mar 2020) $
4# $Author: vondreele $
5# $Revision: 4370 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 4370 2020-03-15 07:57:46Z 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: 4370 $")
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                        hapVary.append(pfx+name)
2506                if 'Layer Disp' in hapData:
2507                    hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0]
2508                    if hapData['Layer Disp'][1]:
2509                        hapVary.append(pfx+'LayerDisp')
2510                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2511                if hapData['Pref.Ori.'][0] == 'MD':
2512                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2513                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2514                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2515                        hapVary.append(pfx+'MD')
2516                else:                           #'SH' spherical harmonics
2517                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2518                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2519                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2520                    controlDict[pfx+'SHhkl'] = []
2521                    try: #patch for old Pref.Ori. items
2522                        controlDict[pfx+'SHtoler'] = 0.1
2523                        if hapData['Pref.Ori.'][6][0] != '':
2524                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2525                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2526                    except IndexError:
2527                        pass
2528                    for item in hapData['Pref.Ori.'][5]:
2529                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2530                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2531                            hapVary.append(pfx+item)
2532                for item in ['Mustrain','Size']:
2533                    controlDict[pfx+item+'Type'] = hapData[item][0]
2534                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2535                    if hapData[item][2][2]:
2536                        hapVary.append(pfx+item+';mx')
2537                    if hapData[item][0] in ['isotropic','uniaxial']:
2538                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2539                        if hapData[item][2][0]:
2540                            hapVary.append(pfx+item+';i')
2541                        if hapData[item][0] == 'uniaxial':
2542                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2543                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2544                            if hapData[item][2][1]:
2545                                hapVary.append(pfx+item+';a')
2546                    else:       #generalized for mustrain or ellipsoidal for size
2547                        Nterms = len(hapData[item][4])
2548                        if item == 'Mustrain':
2549                            names = G2spc.MustrainNames(SGData)
2550                            pwrs = []
2551                            for name in names:
2552                                h,k,l = name[1:]
2553                                pwrs.append([int(h),int(k),int(l)])
2554                            controlDict[pfx+'MuPwrs'] = pwrs
2555                        for i in range(Nterms):
2556                            sfx = ';'+str(i)
2557                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2558                            if hapData[item][5][i]:
2559                                hapVary.append(pfx+item+sfx)
2560                if Phases[phase]['General']['Type'] != 'magnetic':
2561                    for bab in ['BabA','BabU']:
2562                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2563                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2564                            hapVary.append(pfx+bab)
2565                               
2566                if Print: 
2567                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2568                    pFile.write(135*'='+'\n')
2569                    if hapDict[pfx+'LeBail']:
2570                        pFile.write(' Perform LeBail extraction\n')                     
2571                    else:
2572                        pFile.write(' Phase fraction  : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2573                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2574                        if hapData['Pref.Ori.'][0] == 'MD':
2575                            Ax = hapData['Pref.Ori.'][3]
2576                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2577                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2578                        else: #'SH' for spherical harmonics
2579                            PrintSHPO(hapData['Pref.Ori.'])
2580                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2581                    PrintSize(hapData['Size'])
2582                    PrintMuStrain(hapData['Mustrain'],SGData)
2583                    PrintHStrain(hapData['HStrain'],SGData)
2584                    if 'Layer Disp' in hapData:
2585                        pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1]))
2586                    if Phases[phase]['General']['Type'] != 'magnetic':
2587                        if hapData['Babinet']['BabA'][0]:
2588                            PrintBabinet(hapData['Babinet'])
2589                if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2590                    hapData['newLeBail'] = True
2591                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2592                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2593                        hapData['newLeBail'] = False
2594                    refList = []
2595#                    Uniq = []
2596#                    Phi = []
2597                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2598                    if Phases[phase]['General'].get('Modulated',False):
2599                        ifSuper = True
2600                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2601                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2602                        for h,k,l,m,d in HKLd:
2603                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2604                            mul *= 2      # for powder overlap of Friedel pairs
2605                            if m or not ext or useExt:
2606                                if 'C' in inst['Type'][0]:
2607                                    pos = G2lat.Dsp2pos(inst,d)
2608                                    if limits[0] < pos < limits[1]:
2609                                        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])
2610                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2611#                                        Uniq.append(uniq)
2612#                                        Phi.append(phi)
2613                                elif 'T' in inst['Type'][0]:
2614                                    pos = G2lat.Dsp2pos(inst,d)
2615                                    if limits[0] < pos < limits[1]:
2616                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2617                                        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])
2618                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2619                                        #TODO - if tabulated put alp & bet in here
2620#                                        Uniq.append(uniq)
2621#                                        Phi.append(phi)
2622                    else:
2623                        ifSuper = False
2624                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2625                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2626                        for h,k,l,d in HKLd:
2627                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2628                            if ext and 'N' in inst['Type'][0] and  'MagSpGrp' in SGData:
2629                                ext = G2spc.checkMagextc([h,k,l],SGData)
2630                            mul *= 2      # for powder overlap of Friedel pairs
2631                            if ext and not useExt:
2632                                continue
2633                            if 'C' in inst['Type'][0]:
2634                                pos = G2lat.Dsp2pos(inst,d)
2635                                if limits[0] < pos < limits[1]:
2636                                    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])
2637                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2638#                                    Uniq.append(uniq)
2639#                                    Phi.append(phi)
2640                            elif 'T' in inst['Type'][0]:
2641                                pos = G2lat.Dsp2pos(inst,d)
2642                                if limits[0] < pos < limits[1]:
2643                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2644                                    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])
2645                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2646#                                    Uniq.append(uniq)
2647#                                    Phi.append(phi)
2648                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2649            elif 'HKLF' in histogram:
2650                inst = Histogram['Instrument Parameters'][0]
2651                hId = Histogram['hId']
2652                hfx = ':%d:'%(hId)
2653                for item in inst:
2654                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2655                    hapDict[hfx+item] = inst[item][1]
2656                pfx = str(pId)+':'+str(hId)+':'
2657                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2658                if hapData['Scale'][1]:
2659                    hapVary.append(pfx+'Scale')
2660                               
2661                extApprox,extType,extParms = hapData['Extinction']
2662                controlDict[pfx+'EType'] = extType
2663                controlDict[pfx+'EApprox'] = extApprox
2664                if 'C' in inst['Type'][0]:
2665                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2666                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2667                if 'Primary' in extType:
2668                    Ekey = ['Ep',]
2669                elif 'I & II' in extType:
2670                    Ekey = ['Eg','Es']
2671                elif 'Secondary Type II' == extType:
2672                    Ekey = ['Es',]
2673                elif 'Secondary Type I' == extType:
2674                    Ekey = ['Eg',]
2675                else:   #'None'
2676                    Ekey = []
2677                for eKey in Ekey:
2678                    hapDict[pfx+eKey] = extParms[eKey][0]
2679                    if extParms[eKey][1]:
2680                        hapVary.append(pfx+eKey)
2681                for bab in ['BabA','BabU']:
2682                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2683                    if hapData['Babinet'][bab][1]:
2684                        hapVary.append(pfx+bab)
2685                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2686                if len(Twins) == 1:
2687                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2688                    if hapData.get('Flack',[0,False])[1]:
2689                        hapVary.append(pfx+'Flack')
2690                sumTwFr = 0.
2691                controlDict[pfx+'TwinLaw'] = []
2692                controlDict[pfx+'TwinInv'] = []
2693                NTL = 0           
2694                for it,twin in enumerate(Twins):
2695                    if 'bool' in str(type(twin[0])):
2696                        controlDict[pfx+'TwinInv'].append(twin[0])
2697                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2698                    else:
2699                        NTL += 1
2700                        controlDict[pfx+'TwinInv'].append(False)
2701                        controlDict[pfx+'TwinLaw'].append(twin[0])
2702                    if it:
2703                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2704                        sumTwFr += twin[1]
2705                    else:
2706                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2707                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2708                    if Twins[0][1][1]:
2709                        hapVary.append(pfx+'TwinFr:'+str(it))
2710                controlDict[pfx+'NTL'] = NTL
2711                #need to make constraint on TwinFr
2712                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2713                if len(Twins) > 1:    #force sum to unity
2714                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2715                if Print: 
2716                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2717                    pFile.write(135*'='+'\n')
2718                    pFile.write(' Scale factor     : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2719                    if extType != 'None':
2720                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2721                        text = ' Parameters       :'
2722                        for eKey in Ekey:
2723                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2724                        pFile.write(text+'\n')
2725                    if hapData['Babinet']['BabA'][0]:
2726                        PrintBabinet(hapData['Babinet'])
2727                    if not SGData['SGInv'] and len(Twins) == 1:
2728                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2729                    if len(Twins) > 1:
2730                        for it,twin in enumerate(Twins):
2731                            if 'bool' in str(type(twin[0])):
2732                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2733                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2734                            else:
2735                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2736                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2737                       
2738                Histogram['Reflection Lists'] = phase       
2739               
2740    return hapVary,hapDict,controlDict
2741   
2742def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2743    'needs a doc string'
2744   
2745    def PrintSizeAndSig(hapData,sizeSig):
2746        line = '\n Size model:     %9s'%(hapData[0])
2747        refine = False
2748        if hapData[0] in ['isotropic','uniaxial']:
2749            line += ' equatorial:%12.5f'%(hapData[1][0])
2750            if sizeSig[0][0]:
2751                line += ', sig:%8.4f'%(sizeSig[0][0])
2752                refine = True
2753            if hapData[0] == 'uniaxial':
2754                line += ' axial:%12.4f'%(hapData[1][1])
2755                if sizeSig[0][1]:
2756                    refine = True
2757                    line += ', sig:%8.4f'%(sizeSig[0][1])
2758            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2759            if sizeSig[0][2]:
2760                refine = True
2761                line += ', sig:%8.4f'%(sizeSig[0][2])
2762            if refine:
2763                pFile.write(line+'\n')
2764        else:
2765            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2766            if sizeSig[0][2]:
2767                refine = True
2768                line += ', sig:%8.4f'%(sizeSig[0][2])
2769            Snames = ['S11','S22','S33','S12','S13','S23']
2770            ptlbls = ' name  :'
2771            ptstr =  ' value :'
2772            sigstr = ' sig   :'
2773            for i,name in enumerate(Snames):
2774                ptlbls += '%12s' % (name)
2775                ptstr += '%12.6f' % (hapData[4][i])
2776                if sizeSig[1][i]:
2777                    refine = True
2778                    sigstr += '%12.6f' % (sizeSig[1][i])
2779                else:
2780                    sigstr += 12*' '
2781            if refine:
2782                pFile.write(line+'\n')
2783                pFile.write(ptlbls+'\n')
2784                pFile.write(ptstr+'\n')
2785                pFile.write(sigstr+'\n')
2786       
2787    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2788        line = '\n Mustrain model: %9s\n'%(hapData[0])
2789        refine = False
2790        if hapData[0] in ['isotropic','uniaxial']:
2791            line += ' equatorial:%12.1f'%(hapData[1][0])
2792            if mustrainSig[0][0]:
2793                line += ', sig:%8.1f'%(mustrainSig[0][0])
2794                refine = True
2795            if hapData[0] == 'uniaxial':
2796                line += ' axial:%12.1f'%(hapData[1][1])
2797                if mustrainSig[0][1]:
2798                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2799            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2800            if mustrainSig[0][2]:
2801                refine = True
2802                line += ', sig:%8.3f'%(mustrainSig[0][2])
2803            if refine:
2804                pFile.write(line+'\n')
2805        else:
2806            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2807            if mustrainSig[0][2]:
2808                refine = True
2809                line += ', sig:%8.3f'%(mustrainSig[0][2])
2810            Snames = G2spc.MustrainNames(SGData)
2811            ptlbls = ' name  :'
2812            ptstr =  ' value :'
2813            sigstr = ' sig   :'
2814            for i,name in enumerate(Snames):
2815                ptlbls += '%12s' % (name)
2816                ptstr += '%12.1f' % (hapData[4][i])
2817                if mustrainSig[1][i]:
2818                    refine = True
2819                    sigstr += '%12.1f' % (mustrainSig[1][i])
2820                else:
2821                    sigstr += 12*' '
2822            if refine:
2823                pFile.write(line+'\n')
2824                pFile.write(ptlbls+'\n')
2825                pFile.write(ptstr+'\n')
2826                pFile.write(sigstr+'\n')
2827           
2828    def PrintHStrainAndSig(hapData,strainSig,SGData):
2829        Hsnames = G2spc.HStrainNames(SGData)
2830        ptlbls = ' name  :'
2831        ptstr =  ' value :'
2832        sigstr = ' sig   :'
2833        refine = False
2834        for i,name in enumerate(Hsnames):
2835            ptlbls += '%12s' % (name)
2836            ptstr += '%12.4g' % (hapData[0][i])
2837            if name in strainSig:
2838                refine = True
2839                sigstr += '%12.4g' % (strainSig[name])
2840            else:
2841                sigstr += 12*' '
2842        if refine:
2843            pFile.write('\n Hydrostatic/elastic strain:\n')
2844            pFile.write(ptlbls+'\n')
2845            pFile.write(ptstr+'\n')
2846            pFile.write(sigstr+'\n')
2847       
2848    def PrintSHPOAndSig(pfx,hapData,POsig):
2849        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2850        ptlbls = ' names :'
2851        ptstr =  ' values:'
2852        sigstr = ' sig   :'
2853        for item in hapData[5]:
2854            ptlbls += '%12s'%(item)
2855            ptstr += '%12.3f'%(hapData[5][item])
2856            if pfx+item in POsig:
2857                sigstr += '%12.3f'%(POsig[pfx+item])
2858            else:
2859                sigstr += 12*' ' 
2860        pFile.write(ptlbls+'\n')
2861        pFile.write(ptstr+'\n')
2862        pFile.write(sigstr+'\n')
2863       
2864    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2865        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2866        text = ''
2867        for item in hapData[2]:
2868            if pfx+item in ScalExtSig:
2869                text += '       %s: '%(item)
2870                text += '%12.2e'%(hapData[2][item][0])
2871                if pfx+item in ScalExtSig:
2872                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2873        pFile.write(text+'\n')   
2874       
2875    def PrintBabinetAndSig(pfx,hapData,BabSig):
2876        pFile.write('\n Babinet form factor modification:\n')
2877        ptlbls = ' names :'
2878        ptstr =  ' values:'
2879        sigstr = ' sig   :'
2880        for item in hapData:
2881            ptlbls += '%12s'%(item)
2882            ptstr += '%12.3f'%(hapData[item][0])
2883            if pfx+item in BabSig:
2884                sigstr += '%12.3f'%(BabSig[pfx+item])
2885            else:
2886                sigstr += 12*' ' 
2887        pFile.write(ptlbls+'\n')
2888        pFile.write(ptstr+'\n')
2889        pFile.write(sigstr+'\n')
2890       
2891    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2892        pFile.write('\n Twin Law fractions :\n')
2893        ptlbls = ' names :'
2894        ptstr =  ' values:'
2895        sigstr = ' sig   :'
2896        for it,item in enumerate(twinData):
2897            ptlbls += '     twin #%d'%(it)
2898            if it:
2899                ptstr += '%12.3f'%(item[1])
2900            else:
2901                ptstr += '%12.3f'%(item[1][0])
2902            if pfx+'TwinFr:'+str(it) in TwinSig:
2903                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2904            else:
2905                sigstr += 12*' ' 
2906        pFile.write(ptlbls+'\n')
2907        pFile.write(ptstr+'\n')
2908        pFile.write(sigstr+'\n')
2909       
2910   
2911    PhFrExtPOSig = {}
2912    SizeMuStrSig = {}
2913    ScalExtSig = {}
2914    BabSig = {}
2915    TwinFrSig = {}
2916    wtFrSum = {}
2917    for phase in Phases:
2918        HistoPhase = Phases[phase]['Histograms']
2919        General = Phases[phase]['General']
2920        SGData = General['SGData']
2921        pId = Phases[phase]['pId']
2922        histoList = list(HistoPhase.keys())
2923        histoList.sort()
2924        for histogram in histoList:
2925            try:
2926                Histogram = Histograms[histogram]
2927            except KeyError:                       
2928                #skip if histogram not included e.g. in a sequential refinement
2929                continue
2930            if not Phases[phase]['Histograms'][histogram]['Use']:
2931                #skip if phase absent from this histogram
2932                continue
2933            hapData = HistoPhase[histogram]
2934            hId = Histogram['hId']
2935            pfx = str(pId)+':'+str(hId)+':'
2936            if hId not in wtFrSum:
2937                wtFrSum[hId] = 0.
2938            if 'PWDR' in histogram:
2939                for item in ['Scale','Extinction']:
2940                    hapData[item][0] = parmDict[pfx+item]
2941                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2942                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2943                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2944                if hapData['Pref.Ori.'][0] == 'MD':
2945                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2946                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2947                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2948                else:                           #'SH' spherical harmonics
2949                    for item in hapData['Pref.Ori.'][5]:
2950                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2951                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2952                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2953                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2954                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2955                    pfx+'HStrain':{}})                 
2956                for item in ['Mustrain','Size']:
2957                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2958#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2959                    if pfx+item+';mx' in sigDict:
2960                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2961                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2962                        hapData[item][1][0] = parmDict[pfx+item+';i']
2963                        if item == 'Size':
2964                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2965                        if pfx+item+';i' in sigDict: 
2966                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2967                        if hapData[item][0] == 'uniaxial':
2968                            hapData[item][1][1] = parmDict[pfx+item+';a']
2969                            if item == 'Size':
2970                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2971                            if pfx+item+';a' in sigDict:
2972                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2973                    else:       #generalized for mustrain or ellipsoidal for size
2974                        Nterms = len(hapData[item][4])
2975                        for i in range(Nterms):
2976                            sfx = ';'+str(i)
2977                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2978                            if pfx+item+sfx in sigDict:
2979                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2980                names = G2spc.HStrainNames(SGData)
2981                for i,name in enumerate(names):
2982                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2983                    if pfx+name in sigDict:
2984                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2985                if 'Layer Disp' in hapData:
2986                    hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp']
2987                    if pfx+'LayerDisp' in sigDict:
2988                        SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp']
2989                if Phases[phase]['General']['Type'] != 'magnetic':
2990                    for name in ['BabA','BabU']:
2991                        hapData['Babinet'][name][0] = parmDict[pfx+name]
2992                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
2993                            BabSig[pfx+name] = sigDict[pfx+name]               
2994               
2995            elif 'HKLF' in histogram:
2996                for item in ['Scale','Flack']:
2997                    if parmDict.get(pfx+item):
2998                        hapData[item][0] = parmDict[pfx+item]
2999                        if pfx+item in sigDict:
3000                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3001                for item in ['Ep','Eg','Es']:
3002                    if parmDict.get(pfx+item):
3003                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
3004                        if pfx+item in sigDict:
3005                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3006                for item in ['BabA','BabU']:
3007                    hapData['Babinet'][item][0] = parmDict[pfx+item]
3008                    if pfx+item in sigDict:
3009                        BabSig[pfx+item] = sigDict[pfx+item]
3010                if 'Twins' in hapData:
3011                    it = 1
3012                    sumTwFr = 0.
3013                    while True:
3014                        try:
3015                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
3016                            if pfx+'TwinFr:'+str(it) in sigDict:
3017                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
3018                            if it:
3019                                sumTwFr += hapData['Twins'][it][1]
3020                            it += 1
3021                        except KeyError:
3022                            break
3023                    hapData['Twins'][0][1][0] = 1.-sumTwFr
3024
3025    if Print:
3026        for phase in Phases:
3027            HistoPhase = Phases[phase]['Histograms']
3028            General = Phases[phase]['General']
3029            SGData = General['SGData']
3030            pId = Phases[phase]['pId']
3031            histoList = list(HistoPhase.keys())
3032            histoList.sort()
3033            for histogram in histoList:
3034                try:
3035                    Histogram = Histograms[histogram]
3036                except KeyError:                       
3037                    #skip if histogram not included e.g. in a sequential refinement
3038                    continue
3039                hapData = HistoPhase[histogram]
3040                hId = Histogram['hId']
3041                Histogram['Residuals'][str(pId)+'::Name'] = phase
3042                pfx = str(pId)+':'+str(hId)+':'
3043                hfx = ':%s:'%(hId)
3044                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
3045                    continue
3046                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
3047                pFile.write(135*'='+'\n')
3048                if 'PWDR' in histogram:
3049                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
3050                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
3051                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
3052                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
3053                   
3054                    if parmDict[pfx+'LeBail']:
3055                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
3056                    else:
3057                        if pfx+'Scale' in PhFrExtPOSig:
3058                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
3059                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
3060                            pFile.write(' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f\n'%
3061                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
3062                        if pfx+'Extinction' in PhFrExtPOSig:
3063                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
3064                        if hapData['Pref.Ori.'][0] == 'MD':
3065                            if pfx+'MD' in PhFrExtPOSig:
3066                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
3067                        else:
3068                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
3069                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
3070                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
3071                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
3072                    if pfx+'LayerDisp' in SizeMuStrSig:
3073                        pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp']))           
3074                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
3075                        if len(BabSig):
3076                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3077                   
3078                elif 'HKLF' in histogram:
3079                    Inst = Histogram['Instrument Parameters'][0]
3080                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
3081                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
3082                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
3083                    if FFtables != None and 'N' not in Inst['Type'][0]:
3084                        PrintFprime(FFtables,hfx,pFile)
3085                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
3086                    if pfx+'Scale' in ScalExtSig:
3087                        pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
3088                    if hapData['Extinction'][0] != 'None':
3089                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
3090                    if len(BabSig):
3091                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3092                    if pfx+'Flack' in ScalExtSig:
3093                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
3094                    if len(TwinFrSig):
3095                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
3096
3097################################################################################
3098##### Histogram data
3099################################################################################       
3100                   
3101def GetHistogramData(Histograms,Print=True,pFile=None):
3102    'needs a doc string'
3103   
3104    def GetBackgroundParms(hId,Background):
3105        Back = Background[0]
3106        DebyePeaks = Background[1]
3107        bakType,bakFlag = Back[:2]
3108        backVals = Back[3:]
3109        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
3110        backDict = dict(zip(backNames,backVals))
3111        backVary = []
3112        if bakFlag:
3113            backVary = backNames
3114        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
3115        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
3116        debyeDict = {}
3117        debyeList = []
3118        for i in range(DebyePeaks['nDebye']):
3119            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
3120            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
3121            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
3122        debyeVary = []
3123        for item in debyeList:
3124            if item[1]:
3125                debyeVary.append(item[0])
3126        backDict.update(debyeDict)
3127        backVary += debyeVary
3128        peakDict = {}
3129        peakList = []
3130        for i in range(DebyePeaks['nPeaks']):
3131            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
3132                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
3133            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
3134            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
3135        peakVary = []
3136        for item in peakList:
3137            if item[1]:
3138                peakVary.append(item[0])
3139        backDict.update(peakDict)
3140        backVary += peakVary
3141        return bakType,backDict,backVary           
3142       
3143    def GetInstParms(hId,Inst):
3144        #patch
3145        if 'Z' not in Inst:
3146            Inst['Z'] = [0.0,0.0,False]
3147        dataType = Inst['Type'][0]
3148        instDict = {}
3149        insVary = []
3150        pfx = ':'+str(hId)+':'
3151        insKeys = list(Inst.keys())
3152        insKeys.sort()
3153        for item in insKeys:
3154            insName = pfx+item
3155            instDict[insName] = Inst[item][1]
3156            if len(Inst[item]) > 2 and Inst[item][2]:
3157                insVary.append(insName)
3158        if 'C' in dataType:
3159            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
3160        elif 'T' in dataType:   #trap zero alp, bet coeff.
3161            if not instDict[pfx+'alpha']:
3162                instDict[pfx+'alpha'] = 1.0
3163            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3164                instDict[pfx+'beta-1'] = 1.0
3165        return dataType,instDict,insVary
3166       
3167    def GetSampleParms(hId,Sample):
3168        sampVary = []
3169        hfx = ':'+str(hId)+':'       
3170        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
3171            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
3172        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
3173            if key in Sample:
3174                sampDict[hfx+key] = Sample[key]
3175        Type = Sample['Type']
3176        if 'Bragg' in Type:             #Bragg-Brentano
3177            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3178                sampDict[hfx+item] = Sample[item][0]
3179                if Sample[item][1]:
3180                    sampVary.append(hfx+item)
3181        elif 'Debye' in Type:        #Debye-Scherrer
3182            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3183                sampDict[hfx+item] = Sample[item][0]
3184                if Sample[item][1]:
3185                    sampVary.append(hfx+item)
3186        return Type,sampDict,sampVary
3187       
3188    def PrintBackground(Background):
3189        Back = Background[0]
3190        DebyePeaks = Background[1]
3191        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
3192        line = ' Coefficients: '
3193        for i,back in enumerate(Back[3:]):
3194            line += '%10.3f'%(back)
3195            if i and not i%10:
3196                line += '\n'+15*' '
3197        pFile.write(line+'\n')
3198        if DebyePeaks['nDebye']:
3199            pFile.write('\n Debye diffuse scattering coefficients\n')
3200            parms = ['DebyeA','DebyeR','DebyeU']
3201            line = ' names :  '
3202            for parm in parms:
3203                line += '%8s refine?'%(parm)
3204            pFile.write(line+'\n')
3205            for j,term in enumerate(DebyePeaks['debyeTerms']):
3206                line = ' term'+'%2d'%(j)+':'
3207                for i in range(3):
3208                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3209                pFile.write(line+'\n')
3210        if DebyePeaks['nPeaks']:
3211            pFile.write('\n Single peak coefficients\n')
3212            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3213            line = ' names :  '
3214            for parm in parms:
3215                line += '%8s refine?'%(parm)
3216            pFile.write(line+'\n')
3217            for j,term in enumerate(DebyePeaks['peaksList']):
3218                line = ' peak'+'%2d'%(j)+':'
3219                for i in range(4):
3220                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3221                pFile.write(line+'\n')
3222       
3223    def PrintInstParms(Inst):
3224        pFile.write('\n Instrument Parameters:\n')
3225        insKeys = list(Inst.keys())
3226        insKeys.sort()
3227        iBeg = 0
3228        Ok = True
3229        while Ok:
3230            ptlbls = ' name  :'
3231            ptstr =  ' value :'
3232            varstr = ' refine:'
3233            iFin = min(iBeg+9,len(insKeys))
3234            for item in insKeys[iBeg:iFin]:
3235                if item not in ['Type','Source']:
3236                    ptlbls += '%12s' % (item)
3237                    ptstr += '%12.6f' % (Inst[item][1])
3238                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3239                        varstr += 12*' '
3240                    else:
3241                        varstr += '%12s' % (str(bool(Inst[item][2])))
3242            pFile.write(ptlbls+'\n')
3243            pFile.write(ptstr+'\n')
3244            pFile.write(varstr+'\n')
3245            iBeg = iFin
3246            if iBeg == len(insKeys):
3247                Ok = False
3248            else:
3249                pFile.write('\n')
3250       
3251    def PrintSampleParms(Sample):
3252        pFile.write('\n Sample Parameters:\n')
3253        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3254            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3255        ptlbls = ' name  :'
3256        ptstr =  ' value :'
3257        varstr = ' refine:'
3258        if 'Bragg' in Sample['Type']:
3259            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3260                ptlbls += '%14s'%(item)
3261                ptstr += '%14.4f'%(Sample[item][0])
3262                varstr += '%14s'%(str(bool(Sample[item][1])))
3263           
3264        elif 'Debye' in Type:        #Debye-Scherrer
3265            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3266                ptlbls += '%14s'%(item)
3267                ptstr += '%14.4f'%(Sample[item][0])
3268                varstr += '%14s'%(str(bool(Sample[item][1])))
3269
3270        pFile.write(ptlbls+'\n')
3271        pFile.write(ptstr+'\n')
3272        pFile.write(varstr+'\n')
3273       
3274    histDict = {}
3275    histVary = []
3276    controlDict = {}
3277    histoList = list(Histograms.keys())
3278    histoList.sort()
3279    for histogram in histoList:
3280        if 'PWDR' in histogram:
3281            Histogram = Histograms[histogram]
3282            hId = Histogram['hId']
3283            pfx = ':'+str(hId)+':'
3284            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3285            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3286            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3287            for excl in controlDict[pfx+'Exclude']:
3288                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3289            if controlDict[pfx+'Exclude']:
3290                ma.mask_rows(Histogram['Data'])
3291            Background = Histogram['Background']
3292            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3293            controlDict[pfx+'bakType'] = Type
3294            histDict.update(bakDict)
3295            histVary += bakVary
3296           
3297            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3298            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3299                G2fil.G2Print ('Warning: tabulated profile coefficients are ignored')
3300            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3301            controlDict[pfx+'histType'] = Type
3302            if 'XC' in Type:
3303                if pfx+'Lam1' in instDict:
3304                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3305                else:
3306                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3307            histDict.update(instDict)
3308            histVary += insVary
3309           
3310            Sample = Histogram['Sample Parameters']
3311            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3312            controlDict[pfx+'instType'] = Type
3313            histDict.update(sampDict)
3314            histVary += sampVary
3315           
3316   
3317            if Print: 
3318                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3319                pFile.write(135*'='+'\n')
3320                Units = {'C':' deg','T':' msec'}
3321                units = Units[controlDict[pfx+'histType'][2]]
3322                Limits = controlDict[pfx+'Limits']
3323                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3324                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3325                if len(controlDict[pfx+'Exclude']):
3326                    excls = controlDict[pfx+'Exclude']
3327                    for excl in excls:
3328                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3329                PrintSampleParms(Sample)
3330                PrintInstParms(Inst[0])
3331                PrintBackground(Background)
3332        elif 'HKLF' in histogram:
3333            Histogram = Histograms[histogram]
3334            hId = Histogram['hId']
3335            pfx = ':'+str(hId)+':'
3336            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3337            Inst = Histogram['Instrument Parameters'][0]
3338            controlDict[pfx+'histType'] = Inst['Type'][0]
3339            if 'X' in Inst['Type'][0]:
3340                histDict[pfx+'Lam'] = Inst['Lam'][1]
3341                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3342            elif 'NC' in Inst['Type'][0]:                   
3343                histDict[pfx+'Lam'] = Inst['Lam'][1]
3344    return histVary,histDict,controlDict
3345   
3346def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3347    'needs a doc string'
3348   
3349    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3350        Back = Background[0]
3351        DebyePeaks = Background[1]
3352        lenBack = len(Back[3:])
3353        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3354        for i in range(lenBack):
3355            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3356            if pfx+'Back;'+str(i) in sigDict:
3357                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3358        if DebyePeaks['nDebye']:
3359            for i in range(DebyePeaks['nDebye']):
3360                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3361                for j,name in enumerate(names):
3362                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3363                    if name in sigDict:
3364                        backSig[lenBack+3*i+j] = sigDict[name]           
3365        if DebyePeaks['nPeaks']:
3366            for i in range(DebyePeaks['nPeaks']):
3367                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3368                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3369                for j,name in enumerate(names):
3370                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3371                    if name in sigDict:
3372                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3373        return backSig
3374       
3375    def SetInstParms(pfx,Inst,parmDict,sigDict):
3376        instSig = {}
3377        insKeys = list(Inst.keys())
3378        insKeys.sort()
3379        for item in insKeys:
3380            insName = pfx+item
3381            Inst[item][1] = parmDict[insName]
3382            if insName in sigDict:
3383                instSig[item] = sigDict[insName]
3384            else:
3385                instSig[item] = 0
3386        return instSig
3387       
3388    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3389        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3390            sampSig = [0 for i in range(5)]
3391            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3392                Sample[item][0] = parmDict[pfx+item]
3393                if pfx+item in sigDict:
3394                    sampSig[i] = sigDict[pfx+item]
3395        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3396            sampSig = [0 for i in range(4)]
3397            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3398                Sample[item][0] = parmDict[pfx+item]
3399                if pfx+item in sigDict:
3400                    sampSig[i] = sigDict[pfx+item]
3401        return sampSig
3402       
3403    def PrintBackgroundSig(Background,backSig):
3404        Back = Background[0]
3405        DebyePeaks = Background[1]
3406        valstr = ' value : '
3407        sigstr = ' sig   : '
3408        refine = False
3409        for i,back in enumerate(Back[3:]):
3410            valstr += '%10.4g'%(back)
3411            if Back[1]:
3412                refine = True
3413                sigstr += '%10.4g'%(backSig[i])
3414            else:
3415                sigstr += 10*' '
3416        if refine:
3417            pFile.write('\n Background function: %s\n'%Back[0])
3418            pFile.write(valstr+'\n')
3419            pFile.write(sigstr+'\n')
3420        if DebyePeaks['nDebye']:
3421            ifAny = False
3422            ptfmt = "%12.3f"
3423            names =  ' names :'
3424            ptstr =  ' values:'
3425            sigstr = ' esds  :'
3426            for item in sigDict:
3427                if 'Debye' in item:
3428                    ifAny = True
3429                    names += '%12s'%(item)
3430                    ptstr += ptfmt%(parmDict[item])
3431                    sigstr += ptfmt%(sigDict[item])
3432            if ifAny:
3433                pFile.write('\n Debye diffuse scattering coefficients\n')
3434                pFile.write(names+'\n')
3435                pFile.write(ptstr+'\n')
3436                pFile.write(sigstr+'\n')
3437        if DebyePeaks['nPeaks']:
3438            pFile.write('\n Single peak coefficients:\n')
3439            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3440            line = ' peak no. '
3441            for parm in parms:
3442                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3443            pFile.write(line+'\n')
3444            for ip in range(DebyePeaks['nPeaks']):
3445                ptstr = ' %4d '%(ip)
3446                for parm in parms:
3447                    fmt = '%14.3f'
3448                    efmt = '%12.3f'
3449                    if parm == 'BkPkpos':
3450                        fmt = '%14.4f'
3451                        efmt = '%12.4f'
3452                    name = pfx+parm+';%d'%(ip)
3453                    ptstr += fmt%(parmDict[name])
3454                    if name in sigDict:
3455                        ptstr += efmt%(sigDict[name])
3456                    else:
3457                        ptstr += 12*' '
3458                pFile.write(ptstr+'\n')
3459        sumBk = np.array(Histogram['sumBk'])
3460        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3461            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3462       
3463    def PrintInstParmsSig(Inst,instSig):
3464        refine = False
3465        insKeys = list(instSig.keys())
3466        insKeys.sort()
3467        iBeg = 0
3468        Ok = True
3469        while Ok:
3470            ptlbls = ' names :'
3471            ptstr =  ' value :'
3472            sigstr = ' sig   :'
3473            iFin = min(iBeg+9,len(insKeys))
3474            for name in insKeys[iBeg:iFin]:
3475                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3476                    ptlbls += '%12s' % (name)
3477                    ptstr += '%12.6f' % (Inst[name][1])
3478                    if instSig[name]:
3479                        refine = True
3480                        sigstr += '%12.6f' % (instSig[name])
3481                    else:
3482                        sigstr += 12*' '
3483            if refine:
3484                pFile.write('\n Instrument Parameters:\n')
3485                pFile.write(ptlbls+'\n')
3486                pFile.write(ptstr+'\n')
3487                pFile.write(sigstr+'\n')
3488            iBeg = iFin
3489            if iBeg == len(insKeys):
3490                Ok = False
3491       
3492    def PrintSampleParmsSig(Sample,sampleSig):
3493        ptlbls = ' names :'
3494        ptstr =  ' values:'
3495        sigstr = ' sig   :'
3496        refine = False
3497        if 'Bragg' in Sample['Type']:
3498            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3499                ptlbls += '%14s'%(item)
3500                ptstr += '%14.4f'%(Sample[item][0])
3501                if sampleSig[i]:
3502                    refine = True
3503                    sigstr += '%14.4f'%(sampleSig[i])
3504                else:
3505                    sigstr += 14*' '
3506           
3507        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3508            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3509                ptlbls += '%14s'%(item)
3510                ptstr += '%14.4f'%(Sample[item][0])
3511                if sampleSig[i]:
3512                    refine = True
3513                    sigstr += '%14.4f'%(sampleSig[i])
3514                else:
3515                    sigstr += 14*' '
3516
3517        if refine:
3518            pFile.write('\n Sample Parameters:\n')
3519            pFile.write(ptlbls+'\n')
3520            pFile.write(ptstr+'\n')
3521            pFile.write(sigstr+'\n')
3522       
3523    histoList = list(Histograms.keys())
3524    histoList.sort()
3525    for histogram in histoList:
3526        if 'PWDR' in histogram:
3527            Histogram = Histograms[histogram]
3528            hId = Histogram['hId']
3529            pfx = ':'+str(hId)+':'
3530            Background = Histogram['Background']
3531            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3532           
3533            Inst = Histogram['Instrument Parameters'][0]
3534            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3535       
3536            Sample = Histogram['Sample Parameters']
3537            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3538
3539            pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3540            pFile.write(135*'='+'\n')
3541            pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3542            pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3543                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3544            pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3545                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3546            if Print:
3547                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3548                if FFtables != None and 'N' not in Inst['Type'][0]:
3549                    PrintFprime(FFtables,pfx,pFile)
3550                PrintSampleParmsSig(Sample,sampSig)
3551                PrintInstParmsSig(Inst,instSig)
3552                PrintBackgroundSig(Background,backSig)
3553               
Note: See TracBrowser for help on using the repository browser.