source: trunk/GSASIIstrIO.py @ 4813

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

Add 'V' option (vector only) for rigid body fitting

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 171.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-02-12 21:11:08 +0000 (Fri, 12 Feb 2021) $
4# $Author: toby $
5# $Revision: 4813 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 4813 2021-02-12 21:11:08Z toby $
8########### SVN repository information ###################
9'''
10*GSASIIstrIO: structure I/O routines*
11-------------------------------------
12
13Contains routines for reading from GPX files and printing to the .LST file.
14Used for refinements and in G2scriptable.
15
16Should not contain any wxpython references as this should be able to be used
17in non-GUI settings.
18
19'''
20from __future__ import division, print_function
21import platform
22import os
23import os.path as ospath
24import time
25import math
26import copy
27if '2' in platform.python_version_tuple()[0]:
28    import cPickle
29else:
30    import pickle as cPickle
31import numpy as np
32import numpy.ma as ma
33import GSASIIpath
34GSASIIpath.SetVersionNumber("$Revision: 4813 $")
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        A,V = G2mth.Q2AV(RB['Orig'][0])
1372        fixAxis = [0, np.abs(V).argmax()+1]
1373        for i in range(4):
1374            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1375            phaseDict[name] = RB['Orient'][0][i]
1376            if RB['Orient'][1] == 'AV' and i:
1377                phaseVary += [name,]
1378            elif RB['Orient'][1] == 'A' and not i:
1379                phaseVary += [name,]
1380            elif RB['Orient'][1] == 'V' and i not in fixAxis:
1381                phaseVary += [name,]
1382                print('refine',name)
1383    def MakeRBThermals(rbKey,phaseVary,phaseDict):
1384        rbid = str(rbids.index(RB['RBId']))
1385        tlstr = ['11','22','33','12','13','23']
1386        sstr = ['12','13','21','23','31','32','AA','BB']
1387        if 'T' in RB['ThermalMotion'][0]:
1388            pfxRB = pfx+'RB'+rbKey+'T'
1389            for i in range(6):
1390                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1391                phaseDict[name] = RB['ThermalMotion'][1][i]
1392                if RB['ThermalMotion'][2][i]:
1393                    phaseVary += [name,]
1394        if 'L' in RB['ThermalMotion'][0]:
1395            pfxRB = pfx+'RB'+rbKey+'L'
1396            for i in range(6):
1397                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1398                phaseDict[name] = RB['ThermalMotion'][1][i+6]
1399                if RB['ThermalMotion'][2][i+6]:
1400                    phaseVary += [name,]
1401        if 'S' in RB['ThermalMotion'][0]:
1402            pfxRB = pfx+'RB'+rbKey+'S'
1403            for i in range(8):
1404                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1405                phaseDict[name] = RB['ThermalMotion'][1][i+12]
1406                if RB['ThermalMotion'][2][i+12]:
1407                    phaseVary += [name,]
1408        if 'U' in RB['ThermalMotion'][0]:
1409            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1410            phaseDict[name] = RB['ThermalMotion'][1][0]
1411            if RB['ThermalMotion'][2][0]:
1412                phaseVary += [name,]
1413               
1414    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1415        rbid = str(rbids.index(RB['RBId']))
1416        pfxRB = pfx+'RB'+rbKey+'Tr;'
1417        for i,torsion in enumerate(RB['Torsions']):
1418            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1419            phaseDict[name] = torsion[0]
1420            if torsion[1]:
1421                phaseVary += [name,]
1422                   
1423    if Print:
1424        pFile.write('\n Phases:\n')
1425    phaseVary = []
1426    phaseDict = {}
1427    pawleyLookup = {}
1428    FFtables = {}                   #scattering factors - xrays
1429    MFtables = {}                   #Mag. form factors
1430    BLtables = {}                   # neutrons
1431    Natoms = {}
1432    maxSSwave = {}
1433    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1434    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1435    atomIndx = {}
1436    for name in PhaseData:
1437        General = PhaseData[name]['General']
1438        pId = PhaseData[name]['pId']
1439        pfx = str(pId)+'::'
1440        FFtable = G2el.GetFFtable(General['AtomTypes'])
1441        BLtable = G2el.GetBLtable(General)
1442        FFtables.update(FFtable)
1443        BLtables.update(BLtable)
1444        phaseDict[pfx+'isMag'] = False
1445        SGData = General['SGData']
1446        SGtext,SGtable = G2spc.SGPrint(SGData)
1447        if General['Type'] == 'magnetic':
1448            MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g'])
1449            MFtables.update(MFtable)
1450            phaseDict[pfx+'isMag'] = True
1451            SpnFlp = SGData['SpnFlp']
1452        Atoms = PhaseData[name]['Atoms']
1453        if Atoms and not General.get('doPawley'):
1454            cx,ct,cs,cia = General['AtomPtrs']
1455            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1456        PawleyRef = PhaseData[name].get('Pawley ref',[])
1457        cell = General['Cell']
1458        A = G2lat.cell2A(cell[1:7])
1459        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1460            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1461        if cell[0]:
1462            phaseVary += cellVary(pfx,SGData)       #also fills in symmetry required constraints
1463        SSGtext = []    #no superstructure
1464        im = 0
1465        if General.get('Modulated',False):
1466            im = 1      #refl offset
1467            Vec,vRef,maxH = General['SuperVec']
1468            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1469            SSGData = General['SSGData']
1470            SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1471            if vRef:
1472                phaseVary += modVary(pfx,SSGData)       
1473        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1474        if resRBData:
1475            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
1476            for iRB,RB in enumerate(resRBData):
1477                MakeRBParms('R',phaseVary,phaseDict)
1478                MakeRBThermals('R',phaseVary,phaseDict)
1479                MakeRBTorsions('R',phaseVary,phaseDict)
1480       
1481        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1482        if vecRBData:
1483            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
1484            for iRB,RB in enumerate(vecRBData):
1485                MakeRBParms('V',phaseVary,phaseDict)
1486                MakeRBThermals('V',phaseVary,phaseDict)
1487                   
1488        Natoms[pfx] = 0
1489        maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1490        if Atoms and not General.get('doPawley'):
1491            cx,ct,cs,cia = General['AtomPtrs']
1492            Natoms[pfx] = len(Atoms)
1493            for i,at in enumerate(Atoms):
1494                atomIndx[at[cia+8]] = [pfx,i]      #lookup table for restraints
1495                phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1496                    pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1497                    pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1498                    pfx+'AI/A:'+str(i):at[cia],})
1499                if at[cia] == 'I':
1500                    phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1501                else:
1502                    phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1503                        pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1504                if General['Type'] == 'magnetic':
1505                    phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]})
1506                if 'F' in at[ct+1]:
1507                    phaseVary.append(pfx+'Afrac:'+str(i))
1508                if 'X' in at[ct+1]:
1509                    try:    #patch for sytsym name changes
1510                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1511                    except KeyError:
1512                        Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1513                        at[cs] = Sytsym
1514                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1515                    names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1516                    equivs = [[],[],[]]
1517                    for j in range(3):
1518                        if xId[j] > 0:                               
1519                            phaseVary.append(names[j])
1520                            equivs[xId[j]-1].append([names[j],xCoef[j]])
1521                    for equiv in equivs:
1522                        if len(equiv) > 1:
1523                            name = equiv[0][0]
1524                            coef = equiv[0][1]
1525                            for eqv in equiv[1:]:
1526                                eqv[1] /= coef
1527                                G2mv.StoreEquivalence(name,(eqv,))
1528                if 'U' in at[ct+1]:
1529                    if at[cia] == 'I':
1530                        phaseVary.append(pfx+'AUiso:'+str(i))
1531                    else:
1532                        try:    #patch for sytsym name changes
1533                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1534                        except KeyError:
1535                            Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1536                            at[cs] = Sytsym
1537                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1538                        names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1539                            pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1540                        equivs = [[],[],[],[],[],[]]
1541                        for j in range(6):
1542                            if uId[j] > 0:                               
1543                                phaseVary.append(names[j])
1544                                equivs[uId[j]-1].append([names[j],uCoef[j]])
1545                        for equiv in equivs:
1546                            if len(equiv) > 1:
1547                                name = equiv[0][0]
1548                                coef = equiv[0][1]
1549                                for eqv in equiv[1:]:
1550                                    eqv[1] /= coef
1551                                    G2mv.StoreEquivalence(name,(eqv,))
1552                if 'M' in at[ct+1]:
1553                    SytSym,Mul,Nop,dupDir = G2spc.SytSym(at[cx:cx+3],SGData)
1554                    mId,mCoef = G2spc.GetCSpqinel(SpnFlp,dupDir)
1555                    names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1556                    equivs = [[],[],[]]
1557                    for j in range(3):
1558                        if mId[j] > 0:
1559                            phaseVary.append(names[j])
1560                            equivs[mId[j]-1].append([names[j],mCoef[j]])
1561                    for equiv in equivs:
1562                        if len(equiv) > 1:
1563                            name = equiv[0][0]
1564                            coef = equiv[0][1]
1565                            for eqv in equiv[1:]:
1566                                eqv[1] /= coef
1567                                G2mv.StoreEquivalence(name,(eqv,))
1568                if General.get('Modulated',False):
1569                    AtomSS = at[-1]['SS1']
1570                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1571                        Waves = AtomSS[Stype]
1572                        if len(Waves):
1573                            waveType = Waves[0]
1574                        else:
1575                            continue
1576                        phaseDict[pfx+Stype[1].upper()+'waveType:'+str(i)] = waveType
1577                        nx = 0
1578                        for iw,wave in enumerate(Waves[1:]):
1579                            if not iw:
1580                                if waveType in ['ZigZag','Block']:
1581                                    nx = 1
1582                                CSI = G2spc.GetSSfxuinel(waveType,Stype,1,at[cx:cx+3],SGData,SSGData)
1583                            else:
1584                                CSI = G2spc.GetSSfxuinel('Fourier',Stype,iw+1-nx,at[cx:cx+3],SGData,SSGData)
1585                            uId,uCoef = CSI[0]
1586                            stiw = str(i)+':'+str(iw)
1587                            if Stype == 'Spos':
1588                                if waveType in ['ZigZag','Block',] and not iw:
1589                                    names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1590                                    equivs = [[],[], [],[],[]]
1591                                else:
1592                                    names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1593                                        pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1594                                    equivs = [[],[],[], [],[],[]]
1595                            elif Stype == 'Sadp':
1596                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1597                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1598                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1599                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1600                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1601                            elif Stype == 'Sfrac':
1602                                equivs = [[],[]]
1603                                if 'Crenel' in waveType and not iw:
1604                                    names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1605                                else:
1606                                    names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1607                            elif Stype == 'Smag':
1608                                equivs = [[],[],[], [],[],[]]
1609                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1610                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1611                            phaseDict.update(dict(zip(names,wave[0])))
1612                            if wave[1]: #what do we do here for multiple terms in modulation constraints?
1613                                for j in range(len(equivs)):
1614                                    if uId[j][0] > 0:                               
1615                                        phaseVary.append(names[j])
1616                                        equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1617                                for equiv in equivs:
1618                                    if len(equiv) > 1:
1619                                        name = equiv[0][0]
1620                                        coef = equiv[0][1]
1621                                        for eqv in equiv[1:]:
1622                                            eqv[1] /= coef
1623                                            G2mv.StoreEquivalence(name,(eqv,))
1624                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1625            textureData = General['SH Texture']
1626            if textureData['Order'] and not seqRef:
1627                phaseDict[pfx+'SHorder'] = textureData['Order']
1628                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1629                for item in ['omega','chi','phi']:
1630                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1631                    if textureData['Sample '+item][0]:
1632                        phaseVary.append(pfx+'SH '+item)
1633                for item in textureData['SH Coeff'][1]:
1634                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1635                    if textureData['SH Coeff'][0]:
1636                        phaseVary.append(pfx+item)
1637               
1638            if Print:
1639                pFile.write('\n Phase name: %s\n'%General['Name'])
1640                pFile.write(135*'='+'\n')
1641                PrintFFtable(FFtable)
1642                PrintBLtable(BLtable)
1643                if General['Type'] == 'magnetic':
1644                    PrintMFtable(MFtable)
1645                pFile.write('\n')
1646                #how do we print magnetic symmetry table? TBD
1647                if len(SSGtext):    #if superstructure
1648                    for line in SSGtext: pFile.write(line+'\n')
1649                    if len(SSGtable):
1650                        for item in SSGtable:
1651                            line = ' %s '%(item)
1652                            pFile.write(line+'\n') 
1653                    else:
1654                        pFile.write(' ( 1)    %s\n'%(SSGtable[0]))
1655                else:
1656                    for line in SGtext: pFile.write(line+'\n')
1657                    if len(SGtable):
1658                        for item in SGtable:
1659                            line = ' %s '%(item)
1660                            pFile.write(line+'\n') 
1661                    else:
1662                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1663                PrintRBObjects(resRBData,vecRBData)
1664                PrintAtoms(General,Atoms)
1665                if General['Type'] == 'magnetic':
1666                    PrintMoments(General,Atoms)
1667                if General.get('Modulated',False):
1668                    PrintWaves(General,Atoms)
1669                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1670                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1671                if len(SSGtext):    #if superstructure
1672                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1673                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1674                if not seqRef:
1675                    PrintTexture(textureData)
1676                if name in RestraintDict:
1677                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1678                        textureData,RestraintDict[name],pFile)
1679                   
1680        elif PawleyRef:
1681            if Print:
1682                pFile.write('\n Phase name: %s\n'%General['Name'])
1683                pFile.write(135*'='+'\n')
1684                pFile.write('\n')
1685                if len(SSGtext):    #if superstructure
1686                    for line in SSGtext: pFile.write(line+'\n')
1687                    if len(SSGtable):
1688                        for item in SSGtable:
1689                            line = ' %s '%(item)
1690                            pFile.write(line+'\n') 
1691                    else:
1692                        pFile.write(' ( 1)    %s\n'%SSGtable[0])
1693                else:
1694                    for line in SGtext: pFile.write(line+'\n')
1695                    if len(SGtable):
1696                        for item in SGtable:
1697                            line = ' %s '%(item)
1698                            pFile.write(line+'\n')
1699                    else:
1700                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1701                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1702                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1703                if len(SSGtext):    #if superstructure
1704                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1705                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1706            pawleyVary = []
1707            for i,refl in enumerate(PawleyRef):
1708                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1709                if im:
1710                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1711                else:
1712                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1713                if refl[5+im]:
1714                    pawleyVary.append(pfx+'PWLref:'+str(i))
1715            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1716            phaseVary += pawleyVary
1717               
1718    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1719   
1720def cellFill(pfx,SGData,parmDict,sigDict): 
1721    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1722    from the parameter and sig dictionaries.
1723
1724    :param str pfx: parameter prefix ("n::", where n is a phase number)
1725    :param dict SGdata: a symmetry object
1726    :param dict parmDict: a dictionary of parameters
1727    :param dict sigDict:  a dictionary of uncertainties on parameters
1728
1729    :returns: A,sigA where each is a list of six terms with the A terms
1730    '''
1731    if SGData['SGLaue'] in ['-1',]:
1732        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1733            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1734    elif SGData['SGLaue'] in ['2/m',]:
1735        if SGData['SGUniq'] == 'a':
1736            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1737                0,0,parmDict[pfx+'A5']]
1738        elif SGData['SGUniq'] == 'b':
1739            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1740                0,parmDict[pfx+'A4'],0]
1741        else:
1742            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1743                parmDict[pfx+'A3'],0,0]
1744    elif SGData['SGLaue'] in ['mmm',]:
1745        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1746    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1747        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1748    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1749        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1750            parmDict[pfx+'A0'],0,0]
1751    elif SGData['SGLaue'] in ['3R', '3mR']:
1752        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1753            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1754    elif SGData['SGLaue'] in ['m3m','m3']:
1755        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1756
1757    try:
1758        if SGData['SGLaue'] in ['-1',]:
1759            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1760                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1761        elif SGData['SGLaue'] in ['2/m',]:
1762            if SGData['SGUniq'] == 'a':
1763                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1764                    0,0,sigDict[pfx+'A5']]
1765            elif SGData['SGUniq'] == 'b':
1766                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1767                    0,sigDict[pfx+'A4'],0]
1768            else:
1769                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1770                    sigDict[pfx+'A3'],0,0]
1771        elif SGData['SGLaue'] in ['mmm',]:
1772            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1773        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1774            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1775        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1776            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1777        elif SGData['SGLaue'] in ['3R', '3mR']:
1778            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1779        elif SGData['SGLaue'] in ['m3m','m3']:
1780            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1781    except KeyError:
1782        sigA = [0,0,0,0,0,0]
1783    return A,sigA
1784       
1785def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1786    'needs a doc string'
1787    if phaseRest:
1788        Amat = G2lat.cell2AB(cell)[0]
1789        cx,ct,cs = AtPtrs[:3]
1790        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1791            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1792            ['ChemComp','Sites'],['Texture','HKLs']]
1793        for name,rest in names:
1794            itemRest = phaseRest[name]
1795            if rest in itemRest and itemRest[rest] and itemRest['Use']:
1796                pFile.write('\n  %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use'])))
1797                if name in ['Bond','Angle','Plane','Chiral']:
1798                    pFile.write('     calc       obs      sig   delt/sig  atoms(symOp): \n')
1799                    for indx,ops,obs,esd in itemRest[rest]:
1800                        try:
1801                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1802                            AtName = ''
1803                            for i,Aname in enumerate(AtNames):
1804                                AtName += Aname
1805                                if ops[i] == '1':
1806                                    AtName += '-'
1807                                else:
1808                                    AtName += '+('+ops[i]+')-'
1809                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1810                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1811                            if name == 'Bond':
1812                                calc = G2mth.getRestDist(XYZ,Amat)
1813                            elif name == 'Angle':
1814                                calc = G2mth.getRestAngle(XYZ,Amat)
1815                            elif name == 'Plane':
1816                                calc = G2mth.getRestPlane(XYZ,Amat)
1817                            elif name == 'Chiral':
1818                                calc = G2mth.getRestChiral(XYZ,Amat)
1819                            pFile.write(' %9.3f %9.3f %8.3f %8.3f   %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1]))
1820                        except KeyError:
1821                            del itemRest[rest]
1822                elif name in ['Torsion','Rama']:
1823                    pFile.write('  atoms(symOp)  calc  obs  sig  delt/sig  torsions: \n')
1824                    coeffDict = itemRest['Coeff']
1825                    for indx,ops,cofName,esd in itemRest[rest]:
1826                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1827                        AtName = ''
1828                        for i,Aname in enumerate(AtNames):
1829                            AtName += Aname+'+('+ops[i]+')-'
1830                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1831                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1832                        if name == 'Torsion':
1833                            tor = G2mth.getRestTorsion(XYZ,Amat)
1834                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1835                            pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1]))
1836                        else:
1837                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1838                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1839                            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]))
1840                elif name == 'ChemComp':
1841                    pFile.write('     atoms   mul*frac  factor     prod\n')
1842                    for indx,factors,obs,esd in itemRest[rest]:
1843                        try:
1844                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1845                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1846                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1847                            mulfrac = mul*frac
1848                            calcs = mul*frac*factors
1849                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1850                                pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc))
1851                            pFile.write(' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd))
1852                        except KeyError:
1853                            del itemRest[rest]
1854                elif name == 'Texture' and textureData['Order']:
1855                    Start = False
1856                    SHCoef = textureData['SH Coeff'][1]
1857                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1858                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1859                    pFile.write ('    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd \n')
1860                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1861                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1862                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1863                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1864                        Z = ma.masked_greater(Z,0.0)
1865                        num = ma.count(Z)
1866                        sum = 0
1867                        if num:
1868                            sum = np.sum(Z)
1869                        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))
1870       
1871def getCellEsd(pfx,SGData,A,covData):
1872    'needs a doc string'
1873    rVsq = G2lat.calc_rVsq(A)
1874    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1875    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1876    varyList = covData['varyList']
1877    covMatrix = covData['covMatrix']
1878    if len(covMatrix):
1879        vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1880        if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']:
1881            vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0]
1882            vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0]
1883            vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2]
1884        elif SGData['SGLaue'] in ['m3','m3m']:
1885            vcov[0:3,0:3] = vcov[0,0]
1886        elif SGData['SGLaue'] in ['4/m', '4/mmm']:
1887            vcov[0:2,0:2] = vcov[0,0]
1888            vcov[1,2] = vcov[2,1] = vcov[0,2]
1889        elif SGData['SGLaue'] in ['3R','3mR']:
1890            vcov[0:3,0:3] = vcov[0,0]
1891    #        vcov[4,4] = vcov[5,5] = vcov[3,3]
1892            vcov[3:6,3:6] = vcov[3,3]
1893            vcov[0:3,3:6] = vcov[0,3]
1894            vcov[3:6,0:3] = vcov[3,0]
1895    else:
1896        vcov = np.eye(6)
1897    delt = 1.e-9
1898    drVdA = np.zeros(6)
1899    for i in range(6):
1900        A[i] += delt
1901        drVdA[i] = G2lat.calc_rVsq(A)
1902        A[i] -= 2*delt
1903        drVdA[i] -= G2lat.calc_rVsq(A)
1904        A[i] += delt
1905    drVdA /= 2.*delt   
1906    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1907    Vol = 1/np.sqrt(rVsq)
1908    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1909   
1910    dcdA = np.zeros((6,6))
1911    for i in range(6):
1912        pdcdA =np.zeros(6)
1913        A[i] += delt
1914        pdcdA += G2lat.A2cell(A)
1915        A[i] -= 2*delt
1916        pdcdA -= G2lat.A2cell(A)
1917        A[i] += delt
1918        dcdA[i] = pdcdA/(2.*delt)
1919   
1920    sigMat = np.inner(dcdA,np.inner(dcdA,vcov))
1921    var = np.diag(sigMat)
1922    CS = np.where(var>0.,np.sqrt(var),0.)
1923    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']:
1924        CS[3:6] = 0.0
1925    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1926   
1927def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1928    '''Called after a refinement to transfer parameters from the parameter dict to
1929    the phase(s) information read from a GPX file. Also prints values to the .lst file
1930    '''
1931   
1932    def PrintAtomsAndSig(General,Atoms,atomsSig):
1933        pFile.write('\n Atoms:\n')
1934        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1935        if General['Type'] == 'macromolecular':
1936            line = ' res no residue chain '+line
1937        cx,ct,cs,cia = General['AtomPtrs']
1938        pFile.write(line+'\n')
1939        pFile.write(135*'-'+'\n')
1940        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1941            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1942        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1943        for atyp in General['AtomTypes']:       #zero composition
1944            General['NoAtoms'][atyp] = 0.
1945        for i,at in enumerate(Atoms):
1946            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1947            if General['Type'] == 'macromolecular':
1948                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1949                valstr = ' values:          '
1950                sigstr = ' sig   :          '
1951            else:
1952                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1953                valstr = ' values:'
1954                sigstr = ' sig   :'
1955            for ind in range(cx,cx+4):
1956                sigind = str(i)+':'+str(ind)
1957                valstr += fmt[ind]%(at[ind])                   
1958                if sigind in atomsSig:
1959                    sigstr += fmt[ind]%(atomsSig[sigind])
1960                else:
1961                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1962            if at[cia] == 'I':
1963                valstr += fmt[cia+1]%(at[cia+1])
1964                if '%d:%d'%(i,cia+1) in atomsSig:
1965                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1966                else:
1967                    sigstr += 8*' '
1968            else:
1969                valstr += 8*' '
1970                sigstr += 8*' '
1971                for ind in range(cia+2,cia+8):
1972                    sigind = str(i)+':'+str(ind)
1973                    valstr += fmt[ind]%(at[ind])
1974                    if sigind in atomsSig:                       
1975                        sigstr += fmt[ind]%(atomsSig[sigind])
1976                    else:
1977                        sigstr += 8*' '
1978            pFile.write(name+'\n')
1979            pFile.write(valstr+'\n')
1980            pFile.write(sigstr+'\n')
1981           
1982    def PrintMomentsAndSig(General,Atoms,atomsSig):
1983        cell = General['Cell'][1:7]
1984        G = G2lat.fillgmat(cell)
1985        ast = np.sqrt(np.diag(G))
1986        GS = G/np.outer(ast,ast)
1987        pFile.write('\n Magnetic Moments:\n')    #add magnitude & angle, etc.? TBD
1988        line = '   name       Mx        My        Mz       |Mag|'
1989        cx,ct,cs,cia = General['AtomPtrs']
1990        cmx = cx+4
1991        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1992        pFile.write(line+'\n')
1993        pFile.write(135*'-'+'\n')
1994        fmt = {0:'%7s',ct:'%7s',cmx:'%10.3f',cmx+1:'%10.3f',cmx+2:'%10.3f'}
1995        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
1996        for i,at in enumerate(Atoms):
1997            if AtInfo[at[ct]]:
1998                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1999                valstr = ' values:'
2000                sigstr = ' sig   :'
2001                for ind in range(cmx,cmx+3):
2002                    sigind = str(i)+':'+str(ind)
2003                    valstr += fmt[ind]%(at[ind])                   
2004                    if sigind in atomsSig:
2005                        sigstr += fmt[ind]%(atomsSig[sigind])
2006                    else:
2007                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
2008                mag = np.array(at[cmx:cmx+3])
2009                Mag = np.sqrt(np.inner(mag,np.inner(mag,GS)))
2010                valstr += '%10.3f'%Mag
2011                sigstr += 10*' '
2012                pFile.write(name+'\n')
2013                pFile.write(valstr+'\n')
2014                pFile.write(sigstr+'\n')
2015           
2016    def PrintWavesAndSig(General,Atoms,wavesSig):
2017        cx,ct,cs,cia = General['AtomPtrs']
2018        pFile.write('\n Modulation waves\n')
2019        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
2020            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
2021            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
2022        pFile.write(135*'-'+'\n')
2023        for i,at in enumerate(Atoms):
2024            AtomSS = at[-1]['SS1']
2025            for Stype in ['Sfrac','Spos','Sadp','Smag']:
2026                Waves = AtomSS[Stype]
2027                if len(Waves) > 1:
2028                    waveType = Waves[0]
2029                else:
2030                    continue
2031                if len(Waves):
2032                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
2033                        (at[ct-1],at[cs],Stype,waveType))
2034                    for iw,wave in enumerate(Waves[1:]):
2035                        stiw = ':'+str(i)+':'+str(iw)
2036                        namstr = '  names :'
2037                        valstr = '  values:'
2038                        sigstr = '  esds  :'
2039                        if Stype == 'Spos':
2040                            nt = 6
2041                            ot = 0
2042                            if waveType in ['ZigZag','Block',] and not iw:
2043                                nt = 5
2044                                ot = 6
2045                            for j in range(nt):
2046                                name = names['Spos'][j+ot]
2047                                namstr += '%12s'%(name)
2048                                valstr += '%12.4f'%(wave[0][j])
2049                                if name+stiw in wavesSig:
2050                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2051                                else:
2052                                    sigstr += 12*' '
2053                        elif Stype == 'Sfrac':
2054                            ot = 0
2055                            if 'Crenel' in waveType and not iw:
2056                                ot = 2
2057                            for j in range(2):
2058                                name = names['Sfrac'][j+ot]
2059                                namstr += '%12s'%(names['Sfrac'][j+ot])
2060                                valstr += '%12.4f'%(wave[0][j])
2061                                if name+stiw in wavesSig:
2062                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2063                                else:
2064                                    sigstr += 12*' '
2065                        elif Stype == 'Sadp':
2066                            for j in range(12):
2067                                name = names['Sadp'][j]
2068                                namstr += '%10s'%(names['Sadp'][j])
2069                                valstr += '%10.6f'%(wave[0][j])
2070                                if name+stiw in wavesSig:
2071                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
2072                                else:
2073                                    sigstr += 10*' '
2074                        elif Stype == 'Smag':
2075                            for j in range(6):
2076                                name = names['Smag'][j]
2077                                namstr += '%12s'%(names['Smag'][j])
2078                                valstr += '%12.4f'%(wave[0][j])
2079                                if name+stiw in wavesSig:
2080                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2081                                else:
2082                                    sigstr += 12*' '
2083                               
2084                    pFile.write(namstr+'\n')
2085                    pFile.write(valstr+'\n')
2086                    pFile.write(sigstr+'\n')
2087       
2088               
2089    def PrintRBObjPOAndSig(rbfx,rbsx):
2090        for i in WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
2091            pFile.write(i+'\n')
2092       
2093    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
2094        for i in WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
2095            pFile.write(i)
2096       
2097    def PrintRBObjTorAndSig(rbsx):
2098        nTors = len(RBObj['Torsions'])
2099        if nTors:
2100            for i in WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
2101                pFile.write(i)
2102               
2103    def PrintSHtextureAndSig(textureData,SHtextureSig):
2104        Tindx = 1.0
2105        Tvar = 0.0
2106        pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
2107        names = ['omega','chi','phi']
2108        namstr = '  names :'
2109        ptstr =  '  values:'
2110        sigstr = '  esds  :'
2111        for name in names:
2112            namstr += '%12s'%(name)
2113            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2114            if 'Sample '+name in SHtextureSig:
2115                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2116            else:
2117                sigstr += 12*' '
2118        pFile.write(namstr+'\n')
2119        pFile.write(ptstr+'\n')
2120        pFile.write(sigstr+'\n')
2121        pFile.write('\n Texture coefficients:\n')
2122        SHcoeff = textureData['SH Coeff'][1]
2123        SHkeys = list(SHcoeff.keys())
2124        nCoeff = len(SHcoeff)
2125        nBlock = nCoeff//10+1
2126        iBeg = 0
2127        iFin = min(iBeg+10,nCoeff)
2128        for block in range(nBlock):
2129            namstr = '  names :'
2130            ptstr =  '  values:'
2131            sigstr = '  esds  :'
2132            for name in SHkeys[iBeg:iFin]:
2133                namstr += '%12s'%(name)
2134                ptstr += '%12.3f'%(SHcoeff[name])
2135                l = 2.0*eval(name.strip('C'))[0]+1
2136                Tindx += SHcoeff[name]**2/l
2137                if name in SHtextureSig:
2138                    Tvar += (2.*SHcoeff[name]*SHtextureSig[name]/l)**2
2139                    sigstr += '%12.3f'%(SHtextureSig[name])
2140                else:
2141                    sigstr += 12*' '
2142            pFile.write(namstr+'\n')
2143            pFile.write(ptstr+'\n')
2144            pFile.write(sigstr+'\n')
2145            iBeg += 10
2146            iFin = min(iBeg+10,nCoeff)
2147        pFile.write(' Texture index J = %.3f(%d)'%(Tindx,int(1000*np.sqrt(Tvar))))
2148           
2149    ##########################################################################
2150    # SetPhaseData starts here
2151    pFile.write('\n Phases:\n')
2152    for phase in Phases:
2153        pFile.write(' Result for phase: %s\n'%phase)
2154        pFile.write(135*'='+'\n')
2155        Phase = Phases[phase]
2156        General = Phase['General']
2157        SGData = General['SGData']
2158        Atoms = Phase['Atoms']
2159        AtLookup = []
2160        if Atoms and not General.get('doPawley'):
2161            cx,ct,cs,cia = General['AtomPtrs']
2162            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
2163        cell = General['Cell']
2164        pId = Phase['pId']
2165        pfx = str(pId)+'::'
2166        if cell[0]:
2167            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
2168            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
2169            pFile.write(' Reciprocal metric tensor: \n')
2170            ptfmt = "%15.9f"
2171            names = ['A11','A22','A33','A12','A13','A23']
2172            namstr = '  names :'
2173            ptstr =  '  values:'
2174            sigstr = '  esds  :'
2175            for name,a,siga in zip(names,A,sigA):
2176                namstr += '%15s'%(name)
2177                ptstr += ptfmt%(a)
2178                if siga:
2179                    sigstr += ptfmt%(siga)
2180                else:
2181                    sigstr += 15*' '
2182            pFile.write(namstr+'\n')
2183            pFile.write(ptstr+'\n')
2184            pFile.write(sigstr+'\n')
2185            cell[1:7] = G2lat.A2cell(A)
2186            cell[7] = G2lat.calc_V(A)
2187            pFile.write(' New unit cell:\n')
2188            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
2189            names = ['a','b','c','alpha','beta','gamma','Volume']
2190            namstr = '  names :'
2191            ptstr =  '  values:'
2192            sigstr = '  esds  :'
2193            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
2194                namstr += '%12s'%(name)
2195                ptstr += fmt%(a)
2196                if siga:
2197                    sigstr += fmt%(siga)
2198                else:
2199                    sigstr += 12*' '
2200            pFile.write(namstr+'\n')
2201            pFile.write(ptstr+'\n')
2202            pFile.write(sigstr+'\n')
2203        ik = 6  #for Pawley stuff below
2204        if General.get('Modulated',False):
2205            ik = 7
2206            Vec,vRef,maxH = General['SuperVec']
2207            if vRef:
2208                pFile.write(' New modulation vector:\n')
2209                namstr = '  names :'
2210                ptstr =  '  values:'
2211                sigstr = '  esds  :'
2212                for iv,var in enumerate(['mV0','mV1','mV2']):
2213                    namstr += '%12s'%(pfx+var)
2214                    ptstr += '%12.6f'%(parmDict[pfx+var])
2215                    if pfx+var in sigDict:
2216                        Vec[iv] = parmDict[pfx+var]
2217                        sigstr += '%12.6f'%(sigDict[pfx+var])
2218                    else:
2219                        sigstr += 12*' '
2220                pFile.write(namstr+'\n')
2221                pFile.write(ptstr+'\n')
2222                pFile.write(sigstr+'\n')
2223           
2224        General['Mass'] = 0.
2225        if Phase['General'].get('doPawley'):
2226            pawleyRef = Phase['Pawley ref']
2227            for i,refl in enumerate(pawleyRef):
2228                key = pfx+'PWLref:'+str(i)
2229                refl[ik] = parmDict[key]
2230                if key in sigDict:
2231                    refl[ik+1] = sigDict[key]
2232                else:
2233                    refl[ik+1] = 0
2234        else:
2235            VRBIds = RBIds['Vector']
2236            RRBIds = RBIds['Residue']
2237            RBModels = Phase['RBModels']
2238            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2239                jrb = VRBIds.index(RBObj['RBId'])
2240                rbsx = str(irb)+':'+str(jrb)
2241                pFile.write(' Vector rigid body parameters:\n')
2242                PrintRBObjPOAndSig('RBV',rbsx)
2243                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2244            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2245                jrb = RRBIds.index(RBObj['RBId'])
2246                rbsx = str(irb)+':'+str(jrb)
2247                pFile.write(' Residue rigid body parameters:\n')
2248                PrintRBObjPOAndSig('RBR',rbsx)
2249                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2250                PrintRBObjTorAndSig(rbsx)
2251            atomsSig = {}
2252            wavesSig = {}
2253            cx,ct,cs,cia = General['AtomPtrs']
2254            for i,at in enumerate(Atoms):
2255                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2256                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2257                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2258                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2259                for ind in range(cx,cx+4):
2260                    at[ind] = parmDict[names[ind]]
2261                    if ind in range(cx,cx+3):
2262                        name = names[ind].replace('A','dA')
2263                    else:
2264                        name = names[ind]
2265                    if name in sigDict:
2266                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2267                if at[cia] == 'I':
2268                    at[cia+1] = parmDict[names[cia+1]]
2269                    if names[cia+1] in sigDict:
2270                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2271                else:
2272                    for ind in range(cia+2,cia+8):
2273                        at[ind] = parmDict[names[ind]]
2274                        if names[ind] in sigDict:
2275                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2276                if General['Type'] == 'magnetic':
2277                    for ind in range(cx+4,cx+7):
2278                        at[ind] = parmDict[names[ind]]
2279                        if names[ind] in sigDict:
2280                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2281                ind = General['AtomTypes'].index(at[ct])
2282                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2283                if General.get('Modulated',False):
2284                    AtomSS = at[-1]['SS1']
2285                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2286                        Waves = AtomSS[Stype]
2287                        if len(Waves):
2288                            waveType = Waves[0]
2289                        else:
2290                            continue
2291                        for iw,wave in enumerate(Waves[1:]):
2292                            stiw = str(i)+':'+str(iw)
2293                            if Stype == 'Spos':
2294                                if waveType in ['ZigZag','Block',] and not iw:
2295                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2296                                else:
2297                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2298                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2299                            elif Stype == 'Sadp':
2300                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2301                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2302                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2303                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2304                            elif Stype == 'Sfrac':
2305                                if 'Crenel' in waveType and not iw:
2306                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2307                                else:
2308                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2309                            elif Stype == 'Smag':
2310                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2311                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2312                            for iname,name in enumerate(names):
2313                                AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2314                                if pfx+name in sigDict:
2315                                    wavesSig[name] = sigDict[pfx+name]
2316                   
2317            PrintAtomsAndSig(General,Atoms,atomsSig)
2318            if General['Type'] == 'magnetic':
2319                PrintMomentsAndSig(General,Atoms,atomsSig)
2320            if General.get('Modulated',False):
2321                PrintWavesAndSig(General,Atoms,wavesSig)
2322               
2323            density = G2mth.getDensity(General)[0]
2324            pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density))
2325           
2326       
2327        textureData = General['SH Texture']   
2328        if textureData['Order']:
2329            SHtextureSig = {}
2330            for name in ['omega','chi','phi']:
2331                aname = pfx+'SH '+name
2332                textureData['Sample '+name][1] = parmDict[aname]
2333                if aname in sigDict:
2334                    SHtextureSig['Sample '+name] = sigDict[aname]
2335            for name in textureData['SH Coeff'][1]:
2336                aname = pfx+name
2337                textureData['SH Coeff'][1][name] = parmDict[aname]
2338                if aname in sigDict:
2339                    SHtextureSig[name] = sigDict[aname]
2340            PrintSHtextureAndSig(textureData,SHtextureSig)
2341        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2342            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2343                textureData,RestraintDict[phase],pFile)
2344                   
2345################################################################################
2346##### Histogram & Phase data
2347################################################################################       
2348                   
2349def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2350    '''Loads the HAP histogram/phase information into dicts
2351
2352    :param dict Phases: phase information
2353    :param dict Histograms: Histogram information
2354    :param bool Print: prints information as it is read
2355    :param file pFile: file object to print to (the default, None causes printing to the console)
2356    :param bool resetRefList: Should the contents of the Reflection List be initialized
2357      on loading. The default, True, initializes the Reflection List as it is loaded.
2358
2359    :returns: (hapVary,hapDict,controlDict)
2360      * hapVary: list of refined variables
2361      * hapDict: dict with refined variables and their values
2362      * controlDict: dict with fixed parameters
2363    '''
2364   
2365    def PrintSize(hapData):
2366        if hapData[0] in ['isotropic','uniaxial']:
2367            line = '\n Size model    : %9s'%(hapData[0])
2368            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2369            if hapData[0] == 'uniaxial':
2370                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2371            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2372            pFile.write(line+'\n')
2373        else:
2374            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2375                (hapData[0],hapData[1][2],hapData[2][2]))
2376            Snames = ['S11','S22','S33','S12','S13','S23']
2377            ptlbls = ' names :'
2378            ptstr =  ' values:'
2379            varstr = ' refine:'
2380            for i,name in enumerate(Snames):
2381                ptlbls += '%12s' % (name)
2382                ptstr += '%12.3f' % (hapData[4][i])
2383                varstr += '%12s' % (str(hapData[5][i]))
2384            pFile.write(ptlbls+'\n')
2385            pFile.write(ptstr+'\n')
2386            pFile.write(varstr+'\n')
2387       
2388    def PrintMuStrain(hapData,SGData):
2389        if hapData[0] in ['isotropic','uniaxial']:
2390            line = '\n Mustrain model: %9s'%(hapData[0])
2391            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2392            if hapData[0] == 'uniaxial':
2393                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2394            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2395            pFile.write(line+'\n')
2396        else:
2397            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2398                (hapData[0],hapData[1][2],hapData[2][2]))
2399            Snames = G2spc.MustrainNames(SGData)
2400            ptlbls = ' names :'
2401            ptstr =  ' values:'
2402            varstr = ' refine:'
2403            for i,name in enumerate(Snames):
2404                ptlbls += '%12s' % (name)
2405                ptstr += '%12.1f' % (hapData[4][i])
2406                varstr += '%12s' % (str(hapData[5][i]))
2407            pFile.write(ptlbls+'\n')
2408            pFile.write(ptstr+'\n')
2409            pFile.write(varstr+'\n')
2410
2411    def PrintHStrain(hapData,SGData):
2412        pFile.write('\n Hydrostatic/elastic strain:\n')
2413        Hsnames = G2spc.HStrainNames(SGData)
2414        ptlbls = ' names :'
2415        ptstr =  ' values:'
2416        varstr = ' refine:'
2417        for i,name in enumerate(Hsnames):
2418            ptlbls += '%12s' % (name)
2419            ptstr += '%12.4g' % (hapData[0][i])
2420            varstr += '%12s' % (str(hapData[1][i]))
2421        pFile.write(ptlbls+'\n')
2422        pFile.write(ptstr+'\n')
2423        pFile.write(varstr+'\n')
2424
2425    def PrintSHPO(hapData):
2426        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2427        ptlbls = ' names :'
2428        ptstr =  ' values:'
2429        for item in hapData[5]:
2430            ptlbls += '%12s'%(item)
2431            ptstr += '%12.3f'%(hapData[5][item]) 
2432        pFile.write(ptlbls+'\n')
2433        pFile.write(ptstr+'\n')
2434   
2435    def PrintBabinet(hapData):
2436        pFile.write('\n Babinet form factor modification:\n')
2437        ptlbls = ' names :'
2438        ptstr =  ' values:'
2439        varstr = ' refine:'
2440        for name in ['BabA','BabU']:
2441            ptlbls += '%12s' % (name)
2442            ptstr += '%12.6f' % (hapData[name][0])
2443            varstr += '%12s' % (str(hapData[name][1]))
2444        pFile.write(ptlbls+'\n')
2445        pFile.write(ptstr+'\n')
2446        pFile.write(varstr+'\n')
2447       
2448    hapDict = {}
2449    hapVary = []
2450    controlDict = {}
2451   
2452    for phase in Phases:
2453        HistoPhase = Phases[phase]['Histograms']
2454        SGData = Phases[phase]['General']['SGData']
2455        cell = Phases[phase]['General']['Cell'][1:7]
2456        A = G2lat.cell2A(cell)
2457        if Phases[phase]['General'].get('Modulated',False):
2458            SSGData = Phases[phase]['General']['SSGData']
2459            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2460        pId = Phases[phase]['pId']
2461        for histogram in Histograms:
2462            if histogram not in HistoPhase and phase in Histograms[histogram]['Reflection Lists']:
2463                #remove previously created reflection list if histogram is removed from phase
2464                #print("removing ",phase,"from",histogram)
2465                del Histograms[histogram]['Reflection Lists'][phase]
2466        histoList = list(HistoPhase.keys())
2467        histoList.sort()
2468        for histogram in histoList:
2469            try:
2470                Histogram = Histograms[histogram]
2471            except KeyError:                       
2472                #skip if histogram not included e.g. in a sequential refinement
2473                continue
2474            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2475                if phase in Histograms[histogram]['Reflection Lists']:
2476                    del Histograms[histogram]['Reflection Lists'][phase]
2477                continue
2478            hapData = HistoPhase[histogram]
2479            hId = Histogram['hId']
2480            if 'PWDR' in histogram:
2481                limits = Histogram['Limits'][1]
2482                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2483                if 'C' in inst['Type'][1]:
2484                    try:
2485                        wave = inst['Lam'][1]
2486                    except KeyError:
2487                        wave = inst['Lam1'][1]
2488                    dmin = wave/(2.0*sind(limits[1]/2.0))
2489                elif 'T' in inst['Type'][0]:
2490                    dmin = limits[0]/inst['difC'][1]
2491                else:
2492                    wave = inst['Lam'][1]
2493                    dmin = wave/(2.0*sind(limits[1]/2.0))
2494                pfx = str(pId)+':'+str(hId)+':'
2495                if Phases[phase]['General']['doPawley']:
2496                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2497                    hapDict[pfx+'newLeBail'] = True
2498                    Tmin = G2lat.Dsp2pos(inst,dmin)
2499                    if 'T' in inst['Type'][1]:
2500                        limits[0] = max(limits[0],Tmin)
2501                    else:
2502                        limits[1] = min(limits[1],Tmin)
2503                else:
2504                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2505                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2506                if Phases[phase]['General']['Type'] == 'magnetic':
2507                    dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.))
2508                for item in ['Scale','Extinction']:
2509                    hapDict[pfx+item] = hapData[item][0]
2510                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2511                        hapVary.append(pfx+item)
2512                names = G2spc.HStrainNames(SGData)
2513                HSvals = []
2514                for i,name in enumerate(names):
2515                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2516                    HSvals.append(hapDict[pfx+name])
2517                    if hapData['HStrain'][1][i]:
2518                        hapVary.append(pfx+name)
2519                if 'Layer Disp' in hapData:
2520                    hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0]
2521                    if hapData['Layer Disp'][1]:
2522                        hapVary.append(pfx+'LayerDisp')
2523                else:
2524                    hapDict[pfx+'LayerDisp'] = 0.0
2525                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2526                if hapData['Pref.Ori.'][0] == 'MD':
2527                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2528                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2529                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2530                        hapVary.append(pfx+'MD')
2531                else:                           #'SH' spherical harmonics
2532                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2533                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2534                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2535                    controlDict[pfx+'SHhkl'] = []
2536                    try: #patch for old Pref.Ori. items
2537                        controlDict[pfx+'SHtoler'] = 0.1
2538                        if hapData['Pref.Ori.'][6][0] != '':
2539                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2540                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2541                    except IndexError:
2542                        pass
2543                    for item in hapData['Pref.Ori.'][5]:
2544                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2545                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2546                            hapVary.append(pfx+item)
2547                for item in ['Mustrain','Size']:
2548                    controlDict[pfx+item+'Type'] = hapData[item][0]
2549                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2550                    if hapData[item][2][2]:
2551                        hapVary.append(pfx+item+';mx')
2552                    if hapData[item][0] in ['isotropic','uniaxial']:
2553                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2554                        if hapData[item][2][0]:
2555                            hapVary.append(pfx+item+';i')
2556                        if hapData[item][0] == 'uniaxial':
2557                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2558                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2559                            if hapData[item][2][1]:
2560                                hapVary.append(pfx+item+';a')
2561                    else:       #generalized for mustrain or ellipsoidal for size
2562                        Nterms = len(hapData[item][4])
2563                        if item == 'Mustrain':
2564                            names = G2spc.MustrainNames(SGData)
2565                            pwrs = []
2566                            for name in names:
2567                                h,k,l = name[1:]
2568                                pwrs.append([int(h),int(k),int(l)])
2569                            controlDict[pfx+'MuPwrs'] = pwrs
2570                        for i in range(Nterms):
2571                            sfx = ';'+str(i)
2572                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2573                            if hapData[item][5][i]:
2574                                hapVary.append(pfx+item+sfx)
2575                if Phases[phase]['General']['Type'] != 'magnetic':
2576                    for bab in ['BabA','BabU']:
2577                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2578                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2579                            hapVary.append(pfx+bab)
2580                               
2581                if Print: 
2582                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2583                    pFile.write(135*'='+'\n')
2584                    if hapDict[pfx+'LeBail']:
2585                        pFile.write(' Perform LeBail extraction\n')                     
2586                    else:
2587                        pFile.write(' Phase fraction  : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2588                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2589                        if hapData['Pref.Ori.'][0] == 'MD':
2590                            Ax = hapData['Pref.Ori.'][3]
2591                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2592                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2593                        else: #'SH' for spherical harmonics
2594                            PrintSHPO(hapData['Pref.Ori.'])
2595                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2596                    PrintSize(hapData['Size'])
2597                    PrintMuStrain(hapData['Mustrain'],SGData)
2598                    PrintHStrain(hapData['HStrain'],SGData)
2599                    if 'Layer Disp' in hapData:
2600                        pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1]))
2601                    if Phases[phase]['General']['Type'] != 'magnetic':
2602                        if hapData['Babinet']['BabA'][0]:
2603                            PrintBabinet(hapData['Babinet'])
2604                if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2605                    hapData['newLeBail'] = True
2606                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2607                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2608                        hapData['newLeBail'] = False
2609                    refList = []
2610#                    Uniq = []
2611#                    Phi = []
2612                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2613                    if Phases[phase]['General'].get('Modulated',False):
2614                        ifSuper = True
2615                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2616                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2617                        for h,k,l,m,d in HKLd:
2618                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2619                            mul *= 2      # for powder overlap of Friedel pairs
2620                            if m or not ext or useExt:
2621                                if 'C' in inst['Type'][0]:
2622                                    pos = G2lat.Dsp2pos(inst,d)
2623                                    if limits[0] < pos < limits[1]:
2624                                        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])
2625                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2626#                                        Uniq.append(uniq)
2627#                                        Phi.append(phi)
2628                                elif 'T' in inst['Type'][0]:
2629                                    pos = G2lat.Dsp2pos(inst,d)
2630                                    if limits[0] < pos < limits[1]:
2631                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2632                                        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])
2633                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2634                                        #TODO - if tabulated put alp & bet in here
2635#                                        Uniq.append(uniq)
2636#                                        Phi.append(phi)
2637                                elif 'B' in inst['Type'][0]:
2638                                    pos = G2lat.Dsp2pos(inst,d)
2639                                    if limits[0] < pos < limits[1]:
2640                                        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])
2641                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2642                    else:
2643                        ifSuper = False
2644                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2645                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2646                        for h,k,l,d in HKLd:
2647                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2648                            if ext and 'N' in inst['Type'][0] and  'MagSpGrp' in SGData:
2649                                ext = G2spc.checkMagextc([h,k,l],SGData)
2650                            mul *= 2      # for powder overlap of Friedel pairs
2651                            if ext and not useExt:
2652                                continue
2653                            if 'C' in inst['Type'][0]:
2654                                pos = G2lat.Dsp2pos(inst,d)
2655                                if limits[0] < pos < limits[1]:
2656                                    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])
2657                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2658#                                    Uniq.append(uniq)
2659#                                    Phi.append(phi)
2660                            elif 'T' in inst['Type'][0]:
2661                                pos = G2lat.Dsp2pos(inst,d)
2662                                if limits[0] < pos < limits[1]:
2663                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2664                                    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])
2665                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2666#                                    Uniq.append(uniq)
2667#                                    Phi.append(phi)
2668                            elif 'B' in inst['Type'][0]:
2669                                pos = G2lat.Dsp2pos(inst,d)
2670                                if limits[0] < pos < limits[1]:
2671                                    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])
2672                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2673                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2674            elif 'HKLF' in histogram:
2675                inst = Histogram['Instrument Parameters'][0]
2676                hId = Histogram['hId']
2677                hfx = ':%d:'%(hId)
2678                for item in inst:
2679                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2680                    hapDict[hfx+item] = inst[item][1]
2681                pfx = str(pId)+':'+str(hId)+':'
2682                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2683                if hapData['Scale'][1]:
2684                    hapVary.append(pfx+'Scale')
2685                               
2686                extApprox,extType,extParms = hapData['Extinction']
2687                controlDict[pfx+'EType'] = extType
2688                controlDict[pfx+'EApprox'] = extApprox
2689                if 'C' in inst['Type'][0]:
2690                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2691                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2692                if 'Primary' in extType:
2693                    Ekey = ['Ep',]
2694                elif 'I & II' in extType:
2695                    Ekey = ['Eg','Es']
2696                elif 'Secondary Type II' == extType:
2697                    Ekey = ['Es',]
2698                elif 'Secondary Type I' == extType:
2699                    Ekey = ['Eg',]
2700                else:   #'None'
2701                    Ekey = []
2702                for eKey in Ekey:
2703                    hapDict[pfx+eKey] = extParms[eKey][0]
2704                    if extParms[eKey][1]:
2705                        hapVary.append(pfx+eKey)
2706                for bab in ['BabA','BabU']:
2707                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2708                    if hapData['Babinet'][bab][1]:
2709                        hapVary.append(pfx+bab)
2710                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2711                if len(Twins) == 1:
2712                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2713                    if hapData.get('Flack',[0,False])[1]:
2714                        hapVary.append(pfx+'Flack')
2715                sumTwFr = 0.
2716                controlDict[pfx+'TwinLaw'] = []
2717                controlDict[pfx+'TwinInv'] = []
2718                NTL = 0           
2719                for it,twin in enumerate(Twins):
2720                    if 'bool' in str(type(twin[0])):
2721                        controlDict[pfx+'TwinInv'].append(twin[0])
2722                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2723                    else:
2724                        NTL += 1
2725                        controlDict[pfx+'TwinInv'].append(False)
2726                        controlDict[pfx+'TwinLaw'].append(twin[0])
2727                    if it:
2728                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2729                        sumTwFr += twin[1]
2730                    else:
2731                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2732                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2733                    if Twins[0][1][1]:
2734                        hapVary.append(pfx+'TwinFr:'+str(it))
2735                controlDict[pfx+'NTL'] = NTL
2736                #need to make constraint on TwinFr
2737                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2738                if len(Twins) > 1:    #force sum to unity
2739                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2740                if Print: 
2741                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2742                    pFile.write(135*'='+'\n')
2743                    pFile.write(' Scale factor     : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2744                    if extType != 'None':
2745                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2746                        text = ' Parameters       :'
2747                        for eKey in Ekey:
2748                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2749                        pFile.write(text+'\n')
2750                    if hapData['Babinet']['BabA'][0]:
2751                        PrintBabinet(hapData['Babinet'])
2752                    if not SGData['SGInv'] and len(Twins) == 1:
2753                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2754                    if len(Twins) > 1:
2755                        for it,twin in enumerate(Twins):
2756                            if 'bool' in str(type(twin[0])):
2757                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2758                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2759                            else:
2760                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2761                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2762                       
2763                Histogram['Reflection Lists'] = phase       
2764               
2765    return hapVary,hapDict,controlDict
2766   
2767def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,Print=True,pFile=None):
2768    'needs a doc string'
2769   
2770    def PrintSizeAndSig(hapData,sizeSig):
2771        line = '\n Size model:     %9s'%(hapData[0])
2772        refine = False
2773        if hapData[0] in ['isotropic','uniaxial']:
2774            line += ' equatorial:%12.5f'%(hapData[1][0])
2775            if sizeSig[0][0]:
2776                line += ', sig:%8.4f'%(sizeSig[0][0])
2777                refine = True
2778            if hapData[0] == 'uniaxial':
2779                line += ' axial:%12.4f'%(hapData[1][1])
2780                if sizeSig[0][1]:
2781                    refine = True
2782                    line += ', sig:%8.4f'%(sizeSig[0][1])
2783            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2784            if sizeSig[0][2]:
2785                refine = True
2786                line += ', sig:%8.4f'%(sizeSig[0][2])
2787            if refine:
2788                pFile.write(line+'\n')
2789        else:
2790            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2791            if sizeSig[0][2]:
2792                refine = True
2793                line += ', sig:%8.4f'%(sizeSig[0][2])
2794            Snames = ['S11','S22','S33','S12','S13','S23']
2795            ptlbls = ' name  :'
2796            ptstr =  ' value :'
2797            sigstr = ' sig   :'
2798            for i,name in enumerate(Snames):
2799                ptlbls += '%12s' % (name)
2800                ptstr += '%12.6f' % (hapData[4][i])
2801                if sizeSig[1][i]:
2802                    refine = True
2803                    sigstr += '%12.6f' % (sizeSig[1][i])
2804                else:
2805                    sigstr += 12*' '
2806            if refine:
2807                pFile.write(line+'\n')
2808                pFile.write(ptlbls+'\n')
2809                pFile.write(ptstr+'\n')
2810                pFile.write(sigstr+'\n')
2811       
2812    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2813        line = '\n Mustrain model: %9s\n'%(hapData[0])
2814        refine = False
2815        if hapData[0] in ['isotropic','uniaxial']:
2816            line += ' equatorial:%12.1f'%(hapData[1][0])
2817            if mustrainSig[0][0]:
2818                line += ', sig:%8.1f'%(mustrainSig[0][0])
2819                refine = True
2820            if hapData[0] == 'uniaxial':
2821                line += ' axial:%12.1f'%(hapData[1][1])
2822                if mustrainSig[0][1]:
2823                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2824            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2825            if mustrainSig[0][2]:
2826                refine = True
2827                line += ', sig:%8.3f'%(mustrainSig[0][2])
2828            if refine:
2829                pFile.write(line+'\n')
2830        else:
2831            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2832            if mustrainSig[0][2]:
2833                refine = True
2834                line += ', sig:%8.3f'%(mustrainSig[0][2])
2835            Snames = G2spc.MustrainNames(SGData)
2836            ptlbls = ' name  :'
2837            ptstr =  ' value :'
2838            sigstr = ' sig   :'
2839            for i,name in enumerate(Snames):
2840                ptlbls += '%12s' % (name)
2841                ptstr += '%12.1f' % (hapData[4][i])
2842                if mustrainSig[1][i]:
2843                    refine = True
2844                    sigstr += '%12.1f' % (mustrainSig[1][i])
2845                else:
2846                    sigstr += 12*' '
2847            if refine:
2848                pFile.write(line+'\n')
2849                pFile.write(ptlbls+'\n')
2850                pFile.write(ptstr+'\n')
2851                pFile.write(sigstr+'\n')
2852           
2853    def PrintHStrainAndSig(hapData,strainSig,SGData):
2854        Hsnames = G2spc.HStrainNames(SGData)
2855        ptlbls = ' name  :'
2856        ptstr =  ' value :'
2857        sigstr = ' sig   :'
2858        refine = False
2859        for i,name in enumerate(Hsnames):
2860            ptlbls += '%12s' % (name)
2861            ptstr += '%12.4g' % (hapData[0][i])
2862            if name in strainSig:
2863                refine = True
2864                sigstr += '%12.4g' % (strainSig[name])
2865            else:
2866                sigstr += 12*' '
2867        if refine:
2868            pFile.write('\n Hydrostatic/elastic strain:\n')
2869            pFile.write(ptlbls+'\n')
2870            pFile.write(ptstr+'\n')
2871            pFile.write(sigstr+'\n')
2872       
2873    def PrintSHPOAndSig(pfx,hapData,POsig):
2874        Tindx = 1.0
2875        Tvar = 0.0
2876        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2877        ptlbls = ' names :'
2878        ptstr =  ' values:'
2879        sigstr = ' sig   :'
2880        for item in hapData[5]:
2881            ptlbls += '%12s'%(item)
2882            ptstr += '%12.3f'%(hapData[5][item])
2883            l = 2.0*eval(item.strip('C'))[0]+1
2884            Tindx += hapData[5][item]**2/l
2885            if pfx+item in POsig:
2886                Tvar += (2.*hapData[5][item]*POsig[pfx+item]/l)**2
2887                sigstr += '%12.3f'%(POsig[pfx+item])
2888            else:
2889                sigstr += 12*' ' 
2890        pFile.write(ptlbls+'\n')
2891        pFile.write(ptstr+'\n')
2892        pFile.write(sigstr+'\n')
2893        pFile.write('\n Texture index J = %.3f(%d)\n'%(Tindx,int(1000*np.sqrt(Tvar))))
2894       
2895    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2896        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2897        text = ''
2898        for item in hapData[2]:
2899            if pfx+item in ScalExtSig:
2900                text += '       %s: '%(item)
2901                text += '%12.2e'%(hapData[2][item][0])
2902                if pfx+item in ScalExtSig:
2903                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2904        pFile.write(text+'\n')   
2905       
2906    def PrintBabinetAndSig(pfx,hapData,BabSig):
2907        pFile.write('\n Babinet form factor modification:\n')
2908        ptlbls = ' names :'
2909        ptstr =  ' values:'
2910        sigstr = ' sig   :'
2911        for item in hapData:
2912            ptlbls += '%12s'%(item)
2913            ptstr += '%12.3f'%(hapData[item][0])
2914            if pfx+item in BabSig:
2915                sigstr += '%12.3f'%(BabSig[pfx+item])
2916            else:
2917                sigstr += 12*' ' 
2918        pFile.write(ptlbls+'\n')
2919        pFile.write(ptstr+'\n')
2920        pFile.write(sigstr+'\n')
2921       
2922    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2923        pFile.write('\n Twin Law fractions :\n')
2924        ptlbls = ' names :'
2925        ptstr =  ' values:'
2926        sigstr = ' sig   :'
2927        for it,item in enumerate(twinData):
2928            ptlbls += '     twin #%d'%(it)
2929            if it:
2930                ptstr += '%12.3f'%(item[1])
2931            else:
2932                ptstr += '%12.3f'%(item[1][0])
2933            if pfx+'TwinFr:'+str(it) in TwinSig:
2934                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2935            else:
2936                sigstr += 12*' ' 
2937        pFile.write(ptlbls+'\n')
2938        pFile.write(ptstr+'\n')
2939        pFile.write(sigstr+'\n')
2940       
2941   
2942    PhFrExtPOSig = {}
2943    SizeMuStrSig = {}
2944    ScalExtSig = {}
2945    BabSig = {}
2946    TwinFrSig = {}
2947    wtFrSum = {}
2948    for phase in Phases:
2949        HistoPhase = Phases[phase]['Histograms']
2950        General = Phases[phase]['General']
2951        SGData = General['SGData']
2952        pId = Phases[phase]['pId']
2953        histoList = list(HistoPhase.keys())
2954        histoList.sort()
2955        for histogram in histoList:
2956            try:
2957                Histogram = Histograms[histogram]
2958            except KeyError:                       
2959                #skip if histogram not included e.g. in a sequential refinement
2960                continue
2961            if not Phases[phase]['Histograms'][histogram]['Use']:
2962                #skip if phase absent from this histogram
2963                continue
2964            hapData = HistoPhase[histogram]
2965            hId = Histogram['hId']
2966            pfx = str(pId)+':'+str(hId)+':'
2967            if hId not in wtFrSum:
2968                wtFrSum[hId] = 0.
2969            if 'PWDR' in histogram:
2970                for item in ['Scale','Extinction']:
2971                    hapData[item][0] = parmDict[pfx+item]
2972                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2973                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2974                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2975                if hapData['Pref.Ori.'][0] == 'MD':
2976                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2977                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2978                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2979                else:                           #'SH' spherical harmonics
2980                    for item in hapData['Pref.Ori.'][5]:
2981                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2982                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2983                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2984                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2985                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2986                    pfx+'HStrain':{}})                 
2987                for item in ['Mustrain','Size']:
2988                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2989#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2990                    if pfx+item+';mx' in sigDict:
2991                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2992                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2993                        hapData[item][1][0] = parmDict[pfx+item+';i']
2994                        if item == 'Size':
2995                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2996                        if pfx+item+';i' in sigDict: 
2997                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2998                        if hapData[item][0] == 'uniaxial':
2999                            hapData[item][1][1] = parmDict[pfx+item+';a']
3000                            if item == 'Size':
3001                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
3002                            if pfx+item+';a' in sigDict:
3003                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
3004                    else:       #generalized for mustrain or ellipsoidal for size
3005                        Nterms = len(hapData[item][4])
3006                        for i in range(Nterms):
3007                            sfx = ';'+str(i)
3008                            hapData[item][4][i] = parmDict[pfx+item+sfx]
3009                            if pfx+item+sfx in sigDict:
3010                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
3011                names = G2spc.HStrainNames(SGData)
3012                for i,name in enumerate(names):
3013                    hapData['HStrain'][0][i] = parmDict[pfx+name]
3014                    if pfx+name in sigDict:
3015                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
3016                if 'Layer Disp' in hapData:
3017                    hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp']
3018                    if pfx+'LayerDisp' in sigDict:
3019                        SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp']
3020                if Phases[phase]['General']['Type'] != 'magnetic':
3021                    for name in ['BabA','BabU']:
3022                        hapData['Babinet'][name][0] = parmDict[pfx+name]
3023                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
3024                            BabSig[pfx+name] = sigDict[pfx+name]               
3025               
3026            elif 'HKLF' in histogram:
3027                for item in ['Scale','Flack']:
3028                    if parmDict.get(pfx+item):
3029                        hapData[item][0] = parmDict[pfx+item]
3030                        if pfx+item in sigDict:
3031                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3032                for item in ['Ep','Eg','Es']:
3033                    if parmDict.get(pfx+item):
3034                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
3035                        if pfx+item in sigDict:
3036                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3037                for item in ['BabA','BabU']:
3038                    hapData['Babinet'][item][0] = parmDict[pfx+item]
3039                    if pfx+item in sigDict:
3040                        BabSig[pfx+item] = sigDict[pfx+item]
3041                if 'Twins' in hapData:
3042                    it = 1
3043                    sumTwFr = 0.
3044                    while True:
3045                        try:
3046                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
3047                            if pfx+'TwinFr:'+str(it) in sigDict:
3048                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
3049                            if it:
3050                                sumTwFr += hapData['Twins'][it][1]
3051                            it += 1
3052                        except KeyError:
3053                            break
3054                    hapData['Twins'][0][1][0] = 1.-sumTwFr
3055
3056    if Print:
3057        for phase in Phases:
3058            HistoPhase = Phases[phase]['Histograms']
3059            General = Phases[phase]['General']
3060            SGData = General['SGData']
3061            pId = Phases[phase]['pId']
3062            histoList = list(HistoPhase.keys())
3063            histoList.sort()
3064            for histogram in histoList:
3065                try:
3066                    Histogram = Histograms[histogram]
3067                except KeyError:                       
3068                    #skip if histogram not included e.g. in a sequential refinement
3069                    continue
3070                hapData = HistoPhase[histogram]
3071                hId = Histogram['hId']
3072                Histogram['Residuals'][str(pId)+'::Name'] = phase
3073                pfx = str(pId)+':'+str(hId)+':'
3074                hfx = ':%s:'%(hId)
3075                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
3076                    continue
3077                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
3078                pFile.write(135*'='+'\n')
3079                Inst = Histogram['Instrument Parameters'][0]
3080                if 'PWDR' in histogram:
3081                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
3082                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
3083                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
3084                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
3085                   
3086                    if parmDict[pfx+'LeBail']:
3087                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
3088                    else:
3089                        # if calcControls != None:    #skipped in seqRefine
3090                        #     if 'X'in Inst['Type'][0]:
3091                        #         PrintFprime(calcControls['FFtables'],hfx,pFile)
3092                        #     elif 'NC' in Inst['Type'][0]:
3093                        #         PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3094                        if pfx+'Scale' in PhFrExtPOSig:
3095                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
3096                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
3097                            pFile.write(' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f\n'%
3098                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
3099                        if pfx+'Extinction' in PhFrExtPOSig:
3100                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
3101                        if hapData['Pref.Ori.'][0] == 'MD':
3102                            if pfx+'MD' in PhFrExtPOSig:
3103                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
3104                        else:
3105                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
3106                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
3107                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
3108                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
3109                    if pfx+'LayerDisp' in SizeMuStrSig:
3110                        pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp']))           
3111                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
3112                        if len(BabSig):
3113                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3114                   
3115                elif 'HKLF' in histogram:
3116                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
3117                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
3118                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
3119                    if calcControls != None:    #skipped in seqRefine
3120                        if 'X'in Inst['Type'][0]:
3121                            PrintFprime(calcControls['FFtables'],hfx,pFile)
3122                        elif 'NC' in Inst['Type'][0]:
3123                            PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3124                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
3125                    if pfx+'Scale' in ScalExtSig:
3126                        pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
3127                    if hapData['Extinction'][0] != 'None':
3128                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
3129                    if len(BabSig):
3130                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3131                    if pfx+'Flack' in ScalExtSig:
3132                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
3133                    if len(TwinFrSig):
3134                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
3135
3136################################################################################
3137##### Histogram data
3138################################################################################       
3139                   
3140def GetHistogramData(Histograms,Print=True,pFile=None):
3141    'needs a doc string'
3142   
3143    def GetBackgroundParms(hId,Background):
3144        Back = Background[0]
3145        DebyePeaks = Background[1]
3146        bakType,bakFlag = Back[:2]
3147        backVals = Back[3:]
3148        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
3149        backDict = dict(zip(backNames,backVals))
3150        backVary = []
3151        if bakFlag:
3152            backVary = backNames
3153        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
3154        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
3155        debyeDict = {}
3156        debyeList = []
3157        for i in range(DebyePeaks['nDebye']):
3158            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
3159            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
3160            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
3161        debyeVary = []
3162        for item in debyeList:
3163            if item[1]:
3164                debyeVary.append(item[0])
3165        backDict.update(debyeDict)
3166        backVary += debyeVary
3167        peakDict = {}
3168        peakList = []
3169        for i in range(DebyePeaks['nPeaks']):
3170            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
3171                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
3172            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
3173            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
3174        peakVary = []
3175        for item in peakList:
3176            if item[1]:
3177                peakVary.append(item[0])
3178        backDict.update(peakDict)
3179        backVary += peakVary
3180        if 'background PWDR' in Background[1]:
3181            backDict[':'+str(hId)+':Back File'] = Background[1]['background PWDR'][0]
3182            backDict[':'+str(hId)+':BF mult'] = Background[1]['background PWDR'][1]
3183            try:
3184                if Background[1]['background PWDR'][0] and Background[1]['background PWDR'][2]:
3185                    backVary.append(':'+str(hId)+':BF mult')
3186            except IndexError:  # old version without refine flag
3187                pass
3188        return bakType,backDict,backVary           
3189       
3190    def GetInstParms(hId,Inst):
3191        #patch
3192        if 'Z' not in Inst:
3193            Inst['Z'] = [0.0,0.0,False]
3194        dataType = Inst['Type'][0]
3195        instDict = {}
3196        insVary = []
3197        pfx = ':'+str(hId)+':'
3198        insKeys = list(Inst.keys())
3199        insKeys.sort()
3200        for item in insKeys:
3201            if item in 'Azimuth':
3202                continue
3203            insName = pfx+item
3204            instDict[insName] = Inst[item][1]
3205            if len(Inst[item]) > 2 and Inst[item][2]:
3206                insVary.append(insName)
3207        if 'C' in dataType:
3208            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
3209        elif 'T' in dataType:   #trap zero alp, bet coeff.
3210            if not instDict[pfx+'alpha']:
3211                instDict[pfx+'alpha'] = 1.0
3212            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3213                instDict[pfx+'beta-1'] = 1.0
3214        elif 'B' in dataType:   #trap zero alp, bet coeff.
3215            if not instDict[pfx+'alpha-0'] and not instDict[pfx+'alpha-1']:
3216                instDict[pfx+'alpha-1'] = 1.0
3217            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3218                instDict[pfx+'beta-1'] = 1.0
3219        return dataType,instDict,insVary
3220       
3221    def GetSampleParms(hId,Sample):
3222        sampVary = []
3223        hfx = ':'+str(hId)+':'       
3224        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
3225            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi'],hfx+'Azimuth':Sample['Azimuth']}
3226        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
3227            if key in Sample:
3228                sampDict[hfx+key] = Sample[key]
3229        Type = Sample['Type']
3230        if 'Bragg' in Type:             #Bragg-Brentano
3231            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3232                sampDict[hfx+item] = Sample[item][0]
3233                if Sample[item][1]:
3234                    sampVary.append(hfx+item)
3235        elif 'Debye' in Type:        #Debye-Scherrer
3236            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3237                sampDict[hfx+item] = Sample[item][0]
3238                if Sample[item][1]:
3239                    sampVary.append(hfx+item)
3240        return Type,sampDict,sampVary
3241       
3242    def PrintBackground(Background):
3243        Back = Background[0]
3244        DebyePeaks = Background[1]
3245        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
3246        line = ' Coefficients: '
3247        for i,back in enumerate(Back[3:]):
3248            line += '%10.3f'%(back)
3249            if i and not i%10:
3250                line += '\n'+15*' '
3251        pFile.write(line+'\n')
3252        if DebyePeaks['nDebye']:
3253            pFile.write('\n Debye diffuse scattering coefficients\n')
3254            parms = ['DebyeA','DebyeR','DebyeU']
3255            line = ' names :  '
3256            for parm in parms:
3257                line += '%8s refine?'%(parm)
3258            pFile.write(line+'\n')
3259            for j,term in enumerate(DebyePeaks['debyeTerms']):
3260                line = ' term'+'%2d'%(j)+':'
3261                for i in range(3):
3262                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3263                pFile.write(line+'\n')
3264        if DebyePeaks['nPeaks']:
3265            pFile.write('\n Single peak coefficients\n')
3266            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3267            line = ' names :  '
3268            for parm in parms:
3269                line += '%8s refine?'%(parm)
3270            pFile.write(line+'\n')
3271            for j,term in enumerate(DebyePeaks['peaksList']):
3272                line = ' peak'+'%2d'%(j)+':'
3273                for i in range(4):
3274                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3275                pFile.write(line+'\n')
3276        if 'background PWDR' in DebyePeaks:
3277            try:
3278                pFile.write(' Fixed background file: %s; mult: %.3f  Refine? %s\n'%(DebyePeaks['background PWDR'][0],
3279                    DebyePeaks['background PWDR'][1],DebyePeaks['background PWDR'][2]))
3280            except IndexError:  #old version without refine flag
3281                pass
3282       
3283    def PrintInstParms(Inst):
3284        pFile.write('\n Instrument Parameters:\n')
3285        insKeys = list(Inst.keys())
3286        insKeys.sort()
3287        iBeg = 0
3288        Ok = True
3289        while Ok:
3290            ptlbls = ' name  :'
3291            ptstr =  ' value :'
3292            varstr = ' refine:'
3293            iFin = min(iBeg+9,len(insKeys))
3294            for item in insKeys[iBeg:iFin]:
3295                if item not in ['Type','Source']:
3296                    ptlbls += '%12s' % (item)
3297                    ptstr += '%12.6f' % (Inst[item][1])
3298                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3299                        varstr += 12*' '
3300                    else:
3301                        varstr += '%12s' % (str(bool(Inst[item][2])))
3302            pFile.write(ptlbls+'\n')
3303            pFile.write(ptstr+'\n')
3304            pFile.write(varstr+'\n')
3305            iBeg = iFin
3306            if iBeg == len(insKeys):
3307                Ok = False
3308            else:
3309                pFile.write('\n')
3310       
3311    def PrintSampleParms(Sample):
3312        pFile.write('\n Sample Parameters:\n')
3313        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3314            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3315        ptlbls = ' name  :'
3316        ptstr =  ' value :'
3317        varstr = ' refine:'
3318        if 'Bragg' in Sample['Type']:
3319            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3320                ptlbls += '%14s'%(item)
3321                ptstr += '%14.4f'%(Sample[item][0])
3322                varstr += '%14s'%(str(bool(Sample[item][1])))
3323           
3324        elif 'Debye' in Type:        #Debye-Scherrer
3325            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3326                ptlbls += '%14s'%(item)
3327                ptstr += '%14.4f'%(Sample[item][0])
3328                varstr += '%14s'%(str(bool(Sample[item][1])))
3329
3330        pFile.write(ptlbls+'\n')
3331        pFile.write(ptstr+'\n')
3332        pFile.write(varstr+'\n')
3333       
3334    histDict = {}
3335    histVary = []
3336    controlDict = {}
3337    histoList = list(Histograms.keys())
3338    histoList.sort()
3339    for histogram in histoList:
3340        if 'PWDR' in histogram:
3341            Histogram = Histograms[histogram]
3342            hId = Histogram['hId']
3343            pfx = ':'+str(hId)+':'
3344            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3345            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3346            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3347            for excl in controlDict[pfx+'Exclude']:
3348                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3349            if controlDict[pfx+'Exclude']:
3350                ma.mask_rows(Histogram['Data'])
3351            Background = Histogram['Background']
3352            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3353            controlDict[pfx+'bakType'] = Type
3354            histDict.update(bakDict)
3355            histVary += bakVary
3356           
3357            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3358            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3359                G2fil.G2Print ('Warning: tabulated profile coefficients are ignored')
3360            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3361            controlDict[pfx+'histType'] = Type
3362            if 'XC' in Type or 'XB' in Type:
3363                if pfx+'Lam1' in instDict:
3364                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3365                else:
3366                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3367            histDict.update(instDict)
3368            histVary += insVary
3369           
3370            Sample = Histogram['Sample Parameters']
3371            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3372            controlDict[pfx+'instType'] = Type
3373            histDict.update(sampDict)
3374            histVary += sampVary
3375           
3376   
3377            if Print: 
3378                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3379                pFile.write(135*'='+'\n')
3380                Units = {'C':' deg','T':' msec','B':' deg'}
3381                units = Units[controlDict[pfx+'histType'][2]]
3382                Limits = controlDict[pfx+'Limits']
3383                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3384                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3385                if len(controlDict[pfx+'Exclude']):
3386                    excls = controlDict[pfx+'Exclude']
3387                    for excl in excls:
3388                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3389                PrintSampleParms(Sample)
3390                PrintInstParms(Inst[0])
3391                PrintBackground(Background)
3392        elif 'HKLF' in histogram:
3393            Histogram = Histograms[histogram]
3394            hId = Histogram['hId']
3395            pfx = ':'+str(hId)+':'
3396            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3397            Inst = Histogram['Instrument Parameters'][0]
3398            controlDict[pfx+'histType'] = Inst['Type'][0]
3399            if 'X' in Inst['Type'][0]:
3400                histDict[pfx+'Lam'] = Inst['Lam'][1]
3401                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3402            elif 'NC' in Inst['Type'][0] or 'NB' in Inst['Type'][0]:                   
3403                histDict[pfx+'Lam'] = Inst['Lam'][1]
3404    return histVary,histDict,controlDict
3405   
3406def SetHistogramData(parmDict,sigDict,Histograms,calcControls,Print=True,pFile=None,seq=False):
3407    'Shows histogram data after a refinement'
3408   
3409    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3410        Back = Background[0]
3411        DebyePeaks = Background[1]
3412        lenBack = len(Back[3:])
3413        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3414        for i in range(lenBack):
3415            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3416            if pfx+'Back;'+str(i) in sigDict:
3417                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3418        if DebyePeaks['nDebye']:
3419            for i in range(DebyePeaks['nDebye']):
3420                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3421                for j,name in enumerate(names):
3422                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3423                    if name in sigDict:
3424                        backSig[lenBack+3*i+j] = sigDict[name]           
3425        if DebyePeaks['nPeaks']:
3426            for i in range(DebyePeaks['nPeaks']):
3427                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3428                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3429                for j,name in enumerate(names):
3430                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3431                    if name in sigDict:
3432                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3433        if pfx+'BF mult' in sigDict:
3434            DebyePeaks['background PWDR'][1] = parmDict[pfx+'BF mult']
3435            backSig.append(sigDict[pfx+'BF mult'])
3436               
3437        return backSig
3438       
3439    def SetInstParms(pfx,Inst,parmDict,sigDict):
3440        instSig = {}
3441        insKeys = list(Inst.keys())
3442        insKeys.sort()
3443        for item in insKeys:
3444            insName = pfx+item
3445            Inst[item][1] = parmDict[insName]
3446            if insName in sigDict:
3447                instSig[item] = sigDict[insName]
3448            else:
3449                instSig[item] = 0
3450        return instSig
3451       
3452    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3453        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3454            sampSig = [0 for i in range(5)]
3455            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3456                Sample[item][0] = parmDict[pfx+item]
3457                if pfx+item in sigDict:
3458                    sampSig[i] = sigDict[pfx+item]
3459        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3460            sampSig = [0 for i in range(4)]
3461            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3462                Sample[item][0] = parmDict[pfx+item]
3463                if pfx+item in sigDict:
3464                    sampSig[i] = sigDict[pfx+item]
3465        return sampSig
3466       
3467    def PrintBackgroundSig(Background,backSig):
3468        Back = Background[0]
3469        DebyePeaks = Background[1]
3470        valstr = ' value : '
3471        sigstr = ' sig   : '
3472        refine = False
3473        for i,back in enumerate(Back[3:]):
3474            valstr += '%10.4g'%(back)
3475            if Back[1]:
3476                refine = True
3477                sigstr += '%10.4g'%(backSig[i])
3478            else:
3479                sigstr += 10*' '
3480        if refine:
3481            pFile.write('\n Background function: %s\n'%Back[0])
3482            pFile.write(valstr+'\n')
3483            pFile.write(sigstr+'\n')
3484        if DebyePeaks['nDebye']:
3485            ifAny = False
3486            ptfmt = "%12.3f"
3487            names =  ' names :'
3488            ptstr =  ' values:'
3489            sigstr = ' esds  :'
3490            for item in sigDict:
3491                if 'Debye' in item:
3492                    ifAny = True
3493                    names += '%12s'%(item)
3494                    ptstr += ptfmt%(parmDict[item])
3495                    sigstr += ptfmt%(sigDict[item])
3496            if ifAny:
3497                pFile.write('\n Debye diffuse scattering coefficients\n')
3498                pFile.write(names+'\n')
3499                pFile.write(ptstr+'\n')
3500                pFile.write(sigstr+'\n')
3501        if DebyePeaks['nPeaks']:
3502            pFile.write('\n Single peak coefficients:\n')
3503            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3504            line = ' peak no. '
3505            for parm in parms:
3506                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3507            pFile.write(line+'\n')
3508            for ip in range(DebyePeaks['nPeaks']):
3509                ptstr = ' %4d '%(ip)
3510                for parm in parms:
3511                    fmt = '%14.3f'
3512                    efmt = '%12.3f'
3513                    if parm == 'BkPkpos':
3514                        fmt = '%14.4f'
3515                        efmt = '%12.4f'
3516                    name = pfx+parm+';%d'%(ip)
3517                    ptstr += fmt%(parmDict[name])
3518                    if name in sigDict:
3519                        ptstr += efmt%(sigDict[name])
3520                    else:
3521                        ptstr += 12*' '
3522                pFile.write(ptstr+'\n')
3523        if 'background PWDR' in DebyePeaks and DebyePeaks['background PWDR'][2]:
3524            pFile.write(' Fixed background scale: %.3f(%d)\n'%(DebyePeaks['background PWDR'][1],int(1000*backSig[-1])))
3525        sumBk = np.array(Histogram['sumBk'])
3526        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3527            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3528       
3529    def PrintInstParmsSig(Inst,instSig):
3530        refine = False
3531        insKeys = list(instSig.keys())
3532        insKeys.sort()
3533        iBeg = 0
3534        Ok = True
3535        while Ok:
3536            ptlbls = ' names :'
3537            ptstr =  ' value :'
3538            sigstr = ' sig   :'
3539            iFin = min(iBeg+9,len(insKeys))
3540            for name in insKeys[iBeg:iFin]:
3541                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3542                    ptlbls += '%12s' % (name)
3543                    ptstr += '%12.6f' % (Inst[name][1])
3544                    if instSig[name]:
3545                        refine = True
3546                        sigstr += '%12.6f' % (instSig[name])
3547                    else:
3548                        sigstr += 12*' '
3549            if refine:
3550                pFile.write('\n Instrument Parameters:\n')
3551                pFile.write(ptlbls+'\n')
3552                pFile.write(ptstr+'\n')
3553                pFile.write(sigstr+'\n')
3554            iBeg = iFin
3555            if iBeg == len(insKeys):
3556                Ok = False
3557       
3558    def PrintSampleParmsSig(Sample,sampleSig):
3559        ptlbls = ' names :'
3560        ptstr =  ' values:'
3561        sigstr = ' sig   :'
3562        refine = False
3563        if 'Bragg' in Sample['Type']:
3564            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3565                ptlbls += '%14s'%(item)
3566                ptstr += '%14.4f'%(Sample[item][0])
3567                if sampleSig[i]:
3568                    refine = True
3569                    sigstr += '%14.4f'%(sampleSig[i])
3570                else:
3571                    sigstr += 14*' '
3572           
3573        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3574            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3575                ptlbls += '%14s'%(item)
3576                ptstr += '%14.4f'%(Sample[item][0])
3577                if sampleSig[i]:
3578                    refine = True
3579                    sigstr += '%14.4f'%(sampleSig[i])
3580                else:
3581                    sigstr += 14*' '
3582
3583        if refine:
3584            pFile.write('\n Sample Parameters:\n')
3585            pFile.write(ptlbls+'\n')
3586            pFile.write(ptstr+'\n')
3587            pFile.write(sigstr+'\n')
3588       
3589    histoList = list(Histograms.keys())
3590    histoList.sort()
3591    for histogram in histoList:
3592        if 'PWDR' in histogram:
3593            Histogram = Histograms[histogram]
3594            hId = Histogram['hId']
3595            pfx = ':'+str(hId)+':'
3596            Background = Histogram['Background']
3597            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3598           
3599            Inst = Histogram['Instrument Parameters'][0]
3600            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3601       
3602            Sample = Histogram['Sample Parameters']
3603            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3604
3605            if not seq:
3606                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3607                pFile.write(135*'='+'\n')
3608            pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3609            pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3610                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3611            pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3612                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3613            if Print:
3614                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3615                if calcControls != None:    #skipped in seqRefine
3616                    if 'X' in Inst['Type'][0]:
3617                        PrintFprime(calcControls['FFtables'],pfx,pFile)
3618                    elif 'NC' in Inst['Type'][0]:
3619                        PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3620                PrintSampleParmsSig(Sample,sampSig)
3621                PrintInstParmsSig(Inst,instSig)
3622                PrintBackgroundSig(Background,backSig)
3623               
3624def WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
3625    '''Cribbed version of PrintRBObjPOAndSig but returns lists of strings.
3626    Moved so it can be used in ExportCIF
3627    '''
3628    namstr = '  names :'
3629    valstr = '  values:'
3630    sigstr = '  esds  :'
3631    for i,px in enumerate(['Px:','Py:','Pz:']):
3632        name = pfx+rbfx+px+rbsx
3633        namstr += '%12s'%('Pos '+px[1])
3634        valstr += '%12.5f'%(parmDict[name])
3635        if name in sigDict:
3636            sigstr += '%12.5f'%(sigDict[name])
3637        else:
3638            sigstr += 12*' '
3639    for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
3640        name = pfx+rbfx+po+rbsx
3641        namstr += '%12s'%('Ori '+po[1])
3642        valstr += '%12.5f'%(parmDict[name])
3643        if name in sigDict:
3644            sigstr += '%12.5f'%(sigDict[name])
3645        else:
3646            sigstr += 12*' '
3647    return (namstr,valstr,sigstr)
3648
3649def WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
3650    '''Cribbed version of PrintRBObjTLSAndSig but returns lists of strings.
3651    Moved so it can be used in ExportCIF
3652    '''
3653    out = []
3654    namstr = '  names :'
3655    valstr = '  values:'
3656    sigstr = '  esds  :'
3657    if 'N' not in TLS:
3658        out.append(' Thermal motion:\n')
3659    if 'T' in TLS:
3660        for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
3661            name = pfx+rbfx+pt+rbsx
3662            namstr += '%12s'%(pt[:3])
3663            valstr += '%12.5f'%(parmDict[name])
3664            if name in sigDict:
3665                sigstr += '%12.5f'%(sigDict[name])
3666            else:
3667                sigstr += 12*' '
3668        out.append(namstr+'\n')
3669        out.append(valstr+'\n')
3670        out.append(sigstr+'\n')
3671    if 'L' in TLS:
3672        namstr = '  names :'
3673        valstr = '  values:'
3674        sigstr = '  esds  :'
3675        for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
3676            name = pfx+rbfx+pt+rbsx
3677            namstr += '%12s'%(pt[:3])
3678            valstr += '%12.3f'%(parmDict[name])
3679            if name in sigDict:
3680                sigstr += '%12.3f'%(sigDict[name])
3681            else:
3682                sigstr += 12*' '
3683        out.append(namstr+'\n')
3684        out.append(valstr+'\n')
3685        out.append(sigstr+'\n')
3686    if 'S' in TLS:
3687        namstr = '  names :'
3688        valstr = '  values:'
3689        sigstr = '  esds  :'
3690        for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
3691            name = pfx+rbfx+pt+rbsx
3692            namstr += '%12s'%(pt[:3])
3693            valstr += '%12.4f'%(parmDict[name])
3694            if name in sigDict:
3695                sigstr += '%12.4f'%(sigDict[name])
3696            else:
3697                sigstr += 12*' '
3698        out.append(namstr+'\n')
3699        out.append(valstr+'\n')
3700        out.append(sigstr+'\n')
3701    if 'U' in TLS:
3702        name = pfx+rbfx+'U:'+rbsx
3703        namstr = '  names :'+'%12s'%('Uiso')
3704        valstr = '  values:'+'%12.5f'%(parmDict[name])
3705        if name in sigDict:
3706            sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
3707        else:
3708            sigstr = '  esds  :'+12*' '
3709        out.append(namstr+'\n')
3710        out.append(valstr+'\n')
3711        out.append(sigstr+'\n')
3712    return out
3713
3714def WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
3715    '''Cribbed version of PrintRBObjTorAndSig but returns lists of strings.
3716    Moved so it can be used in ExportCIF
3717    '''
3718    out = []
3719    namstr = '  names :'
3720    valstr = '  values:'
3721    sigstr = '  esds  :'
3722    out.append(' Torsions:\n')
3723    for it in range(nTors):
3724        name = pfx+'RBRTr;'+str(it)+':'+rbsx
3725        namstr += '%12s'%('Tor'+str(it))
3726        valstr += '%12.4f'%(parmDict[name])
3727        if name in sigDict:
3728            sigstr += '%12.4f'%(sigDict[name])
3729    out.append(namstr+'\n')
3730    out.append(valstr+'\n')
3731    out.append(sigstr+'\n')
3732    return out
3733
3734def WriteResRBModel(RBModel):
3735    '''Write description of a residue rigid body. Code shifted from
3736    PrintResRBModel to make usable from G2export_CIF
3737    '''
3738    out = []
3739    atNames = RBModel['atNames']
3740    rbRef = RBModel['rbRef']
3741    rbSeq = RBModel['rbSeq']
3742    out.append('    At name       x          y          z\n')
3743    for name,xyz in zip(atNames,RBModel['rbXYZ']):
3744        out.append(%8s %10.4f %10.4f %10.4f\n'%(name,xyz[0],xyz[1],xyz[2]))
3745    out.append('  Orientation defined by: %s -> %s & %s -> %s\n'%
3746        (atNames[rbRef[0]],atNames[rbRef[1]],atNames[rbRef[0]],atNames[rbRef[2]]))
3747    if rbSeq:
3748        for i,rbseq in enumerate(rbSeq):
3749#                nameLst = [atNames[i] for i in rbseq[3]]
3750            out.append('  Torsion sequence %d Bond: %s - %s riding: %s\n'%
3751                    (i,atNames[rbseq[0]],atNames[rbseq[1]],str([atNames[i] for i in rbseq[3]])))
3752    return out
3753
3754def WriteVecRBModel(RBModel,sigDict={},irb=None):
3755    '''Write description of a vector rigid body. Code shifted from
3756    PrintVecRBModel to make usable from G2export_CIF
3757    '''
3758    out = []
3759    rbRef = RBModel['rbRef']
3760    atTypes = RBModel['rbTypes']
3761    for i in range(len(RBModel['VectMag'])):
3762        if sigDict and irb is not None:
3763            name = '::RBV;'+str(i)+':'+str(irb)
3764            s = G2mth.ValEsd(RBModel['VectMag'][i],sigDict.get(name,-.0001))
3765            out.append('Vector no.: %d Magnitude: %s\n'%(i,s))
3766        else:
3767            out.append('Vector no.: %d Magnitude: %8.4f Refine? %s\n'%
3768                           (i,RBModel['VectMag'][i],RBModel['VectRef'][i]))
3769        out.append('  No. Type     vx         vy         vz\n')
3770        for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
3771            out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(j,name,xyz[0],xyz[1],xyz[2]))
3772    out.append('  No. Type      x          y          z\n')
3773    for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
3774        out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(i,name,xyz[0],xyz[1],xyz[2]))
3775    return out
Note: See TracBrowser for help on using the repository browser.