source: trunk/GSASIIstrIO.py

Last change on this file was 5025, checked in by vondreele, 3 days ago

Add shell plots to 3D reflection plots.
Fix escape character issues in CopyRietveld2Origin
Revise LeBail? operation - now works in sequential refinements. Remove "Fit LeBail?" from Calculate menu.
OnLeBail? does 10 cycles instead of 3 in the "zero cycle" step.
Initial FOSQ set to Scale*Phase fr. instead of 1.0 - improves initial LeBaill? fit.
Put newLeBail flag in Controls - as more global.
Allow LeBail? refinements in sequential refinements.
Remove a wx.Yield() - obsolete
Disable deletion of selected parameters in a LeBail? refinement & let the LS handle any singularities that may ensue
Derivatives no longer controlled by LeBail? flag.

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