source: trunk/GSASIIstrIO.py @ 3925

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

bugfix for fixed background: conflict for single xtal

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