source: trunk/GSASIIstrIO.py @ 4814

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

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

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