source: trunk/GSASIIstrIO.py @ 4970

Last change on this file since 4970 was 4970, checked in by vondreele, 4 months ago

fix PlotAAProb bugs - leftover pickradius & an xpos error
fix scale factor bug & set floor as 1.e-12 for PWDR histograms

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