source: trunk/GSASIIstrIO.py @ 4671

Last change on this file since 4671 was 4671, checked in by vondreele, 10 months ago

add refinable multiplier for a fixed background; multiplier is now normally > 0
usable for peak fitting and Rietved refinement
add the fixed background entry to all background defaults
fix bug in GetDetectorXY fo when cursor outside image - returns [0,0] not None; changes elsewhere to use this
GetTthAzmDsp? now returns explicit list not assumed tuple
put the abs in the nl.qr test for singularities in HessianLSQ
Add 'BF mult' to name list in G2obj
put a try - except TypeError? around setting plot style stuff in PlotPatterns?
Remove the setting of Pattern[0]BackFile? - this was redundant for PWDR
remove picker/pickradius from linescan plot - failed
remove the alternate fixed background definition ('_fixedVary', etc.)
clear the PhaseReaderClass?Drawing? dictionary upon import of phase from a gpx file

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