source: trunk/GSASIIstrIO.py @ 4868

Last change on this file since 4868 was 4868, checked in by toby, 7 months ago

fix back print issue + some CIF fixes

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 174.1 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2021-03-30 22:45:34 +0000 (Tue, 30 Mar 2021) $
4# $Author: toby $
5# $Revision: 4868 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 4868 2021-03-30 22:45:34Z toby $
8########### SVN repository information ###################
9'''
10*GSASIIstrIO: structure I/O routines*
11-------------------------------------
12
13Contains routines for reading from GPX files and printing to the .LST file.
14Used for refinements and in G2scriptable.
15
16Should not contain any wxpython references as this should be able to be used
17in non-GUI settings.
18
19'''
20from __future__ import division, print_function
21import platform
22import os
23import os.path as ospath
24import time
25import math
26import copy
27if '2' in platform.python_version_tuple()[0]:
28    import cPickle
29else:
30    import pickle as cPickle
31import numpy as np
32import numpy.ma as ma
33import GSASIIpath
34GSASIIpath.SetVersionNumber("$Revision: 4868 $")
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    'needs a doc string'
1902    rVsq = G2lat.calc_rVsq(A)
1903    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1904    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1905    varyList = covData['varyList']
1906    covMatrix = covData['covMatrix']
1907    if len(covMatrix):
1908        vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1909        if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']:
1910            vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0]
1911            vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0]
1912            vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2]
1913        elif SGData['SGLaue'] in ['m3','m3m']:
1914            vcov[0:3,0:3] = vcov[0,0]
1915        elif SGData['SGLaue'] in ['4/m', '4/mmm']:
1916            vcov[0:2,0:2] = vcov[0,0]
1917            vcov[1,2] = vcov[2,1] = vcov[0,2]
1918        elif SGData['SGLaue'] in ['3R','3mR']:
1919            vcov[0:3,0:3] = vcov[0,0]
1920    #        vcov[4,4] = vcov[5,5] = vcov[3,3]
1921            vcov[3:6,3:6] = vcov[3,3]
1922            vcov[0:3,3:6] = vcov[0,3]
1923            vcov[3:6,0:3] = vcov[3,0]
1924    else:
1925        vcov = np.eye(6)
1926    delt = 1.e-9
1927    drVdA = np.zeros(6)
1928    for i in range(6):
1929        A[i] += delt
1930        drVdA[i] = G2lat.calc_rVsq(A)
1931        A[i] -= 2*delt
1932        drVdA[i] -= G2lat.calc_rVsq(A)
1933        A[i] += delt
1934    drVdA /= 2.*delt   
1935    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1936    Vol = 1/np.sqrt(rVsq)
1937    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1938   
1939    dcdA = np.zeros((6,6))
1940    for i in range(6):
1941        pdcdA =np.zeros(6)
1942        A[i] += delt
1943        pdcdA += G2lat.A2cell(A)
1944        A[i] -= 2*delt
1945        pdcdA -= G2lat.A2cell(A)
1946        A[i] += delt
1947        dcdA[i] = pdcdA/(2.*delt)
1948   
1949    sigMat = np.inner(dcdA,np.inner(dcdA,vcov))
1950    var = np.diag(sigMat)
1951    CS = np.where(var>0.,np.sqrt(var),0.)
1952    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']:
1953        CS[3:6] = 0.0
1954    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1955   
1956def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1957    '''Called after a refinement to transfer parameters from the parameter dict to
1958    the phase(s) information read from a GPX file. Also prints values to the .lst file
1959    '''
1960   
1961    def PrintAtomsAndSig(General,Atoms,atomsSig):
1962        pFile.write('\n Atoms:\n')
1963        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1964        if General['Type'] == 'macromolecular':
1965            line = ' res no residue chain '+line
1966        cx,ct,cs,cia = General['AtomPtrs']
1967        pFile.write(line+'\n')
1968        pFile.write(135*'-'+'\n')
1969        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1970            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1971        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1972        for atyp in General['AtomTypes']:       #zero composition
1973            General['NoAtoms'][atyp] = 0.
1974        for i,at in enumerate(Atoms):
1975            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1976            if General['Type'] == 'macromolecular':
1977                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1978                valstr = ' values:          '
1979                sigstr = ' sig   :          '
1980            else:
1981                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1982                valstr = ' values:'
1983                sigstr = ' sig   :'
1984            for ind in range(cx,cx+4):
1985                sigind = str(i)+':'+str(ind)
1986                valstr += fmt[ind]%(at[ind])                   
1987                if sigind in atomsSig:
1988                    sigstr += fmt[ind]%(atomsSig[sigind])
1989                else:
1990                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1991            if at[cia] == 'I':
1992                valstr += fmt[cia+1]%(at[cia+1])
1993                if '%d:%d'%(i,cia+1) in atomsSig:
1994                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1995                else:
1996                    sigstr += 8*' '
1997            else:
1998                valstr += 8*' '
1999                sigstr += 8*' '
2000                for ind in range(cia+2,cia+8):
2001                    sigind = str(i)+':'+str(ind)
2002                    valstr += fmt[ind]%(at[ind])
2003                    if sigind in atomsSig:                       
2004                        sigstr += fmt[ind]%(atomsSig[sigind])
2005                    else:
2006                        sigstr += 8*' '
2007            pFile.write(name+'\n')
2008            pFile.write(valstr+'\n')
2009            pFile.write(sigstr+'\n')
2010           
2011    def PrintMomentsAndSig(General,Atoms,atomsSig):
2012        cell = General['Cell'][1:7]
2013        G = G2lat.fillgmat(cell)
2014        ast = np.sqrt(np.diag(G))
2015        GS = G/np.outer(ast,ast)
2016        pFile.write('\n Magnetic Moments:\n')    #add magnitude & angle, etc.? TBD
2017        line = '   name       Mx        My        Mz       |Mag|'
2018        cx,ct,cs,cia = General['AtomPtrs']
2019        cmx = cx+4
2020        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
2021        pFile.write(line+'\n')
2022        pFile.write(135*'-'+'\n')
2023        fmt = {0:'%7s',ct:'%7s',cmx:'%10.3f',cmx+1:'%10.3f',cmx+2:'%10.3f'}
2024        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
2025        for i,at in enumerate(Atoms):
2026            if AtInfo[at[ct]]:
2027                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
2028                valstr = ' values:'
2029                sigstr = ' sig   :'
2030                for ind in range(cmx,cmx+3):
2031                    sigind = str(i)+':'+str(ind)
2032                    valstr += fmt[ind]%(at[ind])                   
2033                    if sigind in atomsSig:
2034                        sigstr += fmt[ind]%(atomsSig[sigind])
2035                    else:
2036                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
2037                mag = np.array(at[cmx:cmx+3])
2038                Mag = np.sqrt(np.inner(mag,np.inner(mag,GS)))
2039                valstr += '%10.3f'%Mag
2040                sigstr += 10*' '
2041                pFile.write(name+'\n')
2042                pFile.write(valstr+'\n')
2043                pFile.write(sigstr+'\n')
2044           
2045    def PrintWavesAndSig(General,Atoms,wavesSig):
2046        cx,ct,cs,cia = General['AtomPtrs']
2047        pFile.write('\n Modulation waves\n')
2048        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
2049            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
2050            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
2051        pFile.write(135*'-'+'\n')
2052        for i,at in enumerate(Atoms):
2053            AtomSS = at[-1]['SS1']
2054            for Stype in ['Sfrac','Spos','Sadp','Smag']:
2055                Waves = AtomSS[Stype]
2056                if len(Waves) > 1:
2057                    waveType = Waves[0]
2058                else:
2059                    continue
2060                if len(Waves):
2061                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
2062                        (at[ct-1],at[cs],Stype,waveType))
2063                    for iw,wave in enumerate(Waves[1:]):
2064                        stiw = ':'+str(i)+':'+str(iw)
2065                        namstr = '  names :'
2066                        valstr = '  values:'
2067                        sigstr = '  esds  :'
2068                        if Stype == 'Spos':
2069                            nt = 6
2070                            ot = 0
2071                            if waveType in ['ZigZag','Block',] and not iw:
2072                                nt = 5
2073                                ot = 6
2074                            for j in range(nt):
2075                                name = names['Spos'][j+ot]
2076                                namstr += '%12s'%(name)
2077                                valstr += '%12.4f'%(wave[0][j])
2078                                if name+stiw in wavesSig:
2079                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2080                                else:
2081                                    sigstr += 12*' '
2082                        elif Stype == 'Sfrac':
2083                            ot = 0
2084                            if 'Crenel' in waveType and not iw:
2085                                ot = 2
2086                            for j in range(2):
2087                                name = names['Sfrac'][j+ot]
2088                                namstr += '%12s'%(names['Sfrac'][j+ot])
2089                                valstr += '%12.4f'%(wave[0][j])
2090                                if name+stiw in wavesSig:
2091                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2092                                else:
2093                                    sigstr += 12*' '
2094                        elif Stype == 'Sadp':
2095                            for j in range(12):
2096                                name = names['Sadp'][j]
2097                                namstr += '%10s'%(names['Sadp'][j])
2098                                valstr += '%10.6f'%(wave[0][j])
2099                                if name+stiw in wavesSig:
2100                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
2101                                else:
2102                                    sigstr += 10*' '
2103                        elif Stype == 'Smag':
2104                            for j in range(6):
2105                                name = names['Smag'][j]
2106                                namstr += '%12s'%(names['Smag'][j])
2107                                valstr += '%12.4f'%(wave[0][j])
2108                                if name+stiw in wavesSig:
2109                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
2110                                else:
2111                                    sigstr += 12*' '
2112                               
2113                    pFile.write(namstr+'\n')
2114                    pFile.write(valstr+'\n')
2115                    pFile.write(sigstr+'\n')
2116       
2117               
2118    def PrintRBObjPOAndSig(rbfx,rbsx):
2119        for i in WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
2120            pFile.write(i+'\n')
2121       
2122    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
2123        for i in WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
2124            pFile.write(i)
2125       
2126    def PrintRBObjTorAndSig(rbsx):
2127        nTors = len(RBObj['Torsions'])
2128        if nTors:
2129            for i in WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
2130                pFile.write(i)
2131               
2132    def PrintSHtextureAndSig(textureData,SHtextureSig):
2133        Tindx = 1.0
2134        Tvar = 0.0
2135        pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
2136        names = ['omega','chi','phi']
2137        namstr = '  names :'
2138        ptstr =  '  values:'
2139        sigstr = '  esds  :'
2140        for name in names:
2141            namstr += '%12s'%(name)
2142            ptstr += '%12.3f'%(textureData['Sample '+name][1])
2143            if 'Sample '+name in SHtextureSig:
2144                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
2145            else:
2146                sigstr += 12*' '
2147        pFile.write(namstr+'\n')
2148        pFile.write(ptstr+'\n')
2149        pFile.write(sigstr+'\n')
2150        pFile.write('\n Texture coefficients:\n')
2151        SHcoeff = textureData['SH Coeff'][1]
2152        SHkeys = list(SHcoeff.keys())
2153        nCoeff = len(SHcoeff)
2154        nBlock = nCoeff//10+1
2155        iBeg = 0
2156        iFin = min(iBeg+10,nCoeff)
2157        for block in range(nBlock):
2158            namstr = '  names :'
2159            ptstr =  '  values:'
2160            sigstr = '  esds  :'
2161            for name in SHkeys[iBeg:iFin]:
2162                namstr += '%12s'%(name)
2163                ptstr += '%12.3f'%(SHcoeff[name])
2164                l = 2.0*eval(name.strip('C'))[0]+1
2165                Tindx += SHcoeff[name]**2/l
2166                if name in SHtextureSig:
2167                    Tvar += (2.*SHcoeff[name]*SHtextureSig[name]/l)**2
2168                    sigstr += '%12.3f'%(SHtextureSig[name])
2169                else:
2170                    sigstr += 12*' '
2171            pFile.write(namstr+'\n')
2172            pFile.write(ptstr+'\n')
2173            pFile.write(sigstr+'\n')
2174            iBeg += 10
2175            iFin = min(iBeg+10,nCoeff)
2176        pFile.write(' Texture index J = %.3f(%d)'%(Tindx,int(1000*np.sqrt(Tvar))))
2177           
2178    ##########################################################################
2179    # SetPhaseData starts here
2180    if pFile: pFile.write('\n Phases:\n')
2181    for phase in Phases:
2182        if pFile: pFile.write(' Result for phase: %s\n'%phase)
2183        if pFile: pFile.write(135*'='+'\n')
2184        Phase = Phases[phase]
2185        General = Phase['General']
2186        SGData = General['SGData']
2187        Atoms = Phase['Atoms']
2188        AtLookup = []
2189        if Atoms and not General.get('doPawley'):
2190            cx,ct,cs,cia = General['AtomPtrs']
2191            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
2192        cell = General['Cell']
2193        pId = Phase['pId']
2194        pfx = str(pId)+'::'
2195        if cell[0]:
2196            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
2197            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
2198            if pFile: pFile.write(' Reciprocal metric tensor: \n')
2199            ptfmt = "%15.9f"
2200            names = ['A11','A22','A33','A12','A13','A23']
2201            namstr = '  names :'
2202            ptstr =  '  values:'
2203            sigstr = '  esds  :'
2204            for name,a,siga in zip(names,A,sigA):
2205                namstr += '%15s'%(name)
2206                ptstr += ptfmt%(a)
2207                if siga:
2208                    sigstr += ptfmt%(siga)
2209                else:
2210                    sigstr += 15*' '
2211            if pFile: pFile.write(namstr+'\n')
2212            if pFile: pFile.write(ptstr+'\n')
2213            if pFile: pFile.write(sigstr+'\n')
2214            cell[1:7] = G2lat.A2cell(A)
2215            cell[7] = G2lat.calc_V(A)
2216            if pFile: pFile.write(' New unit cell:\n')
2217            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
2218            names = ['a','b','c','alpha','beta','gamma','Volume']
2219            namstr = '  names :'
2220            ptstr =  '  values:'
2221            sigstr = '  esds  :'
2222            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
2223                namstr += '%12s'%(name)
2224                ptstr += fmt%(a)
2225                if siga:
2226                    sigstr += fmt%(siga)
2227                else:
2228                    sigstr += 12*' '
2229            if pFile: pFile.write(namstr+'\n')
2230            if pFile: pFile.write(ptstr+'\n')
2231            if pFile: pFile.write(sigstr+'\n')
2232        ik = 6  #for Pawley stuff below
2233        if General.get('Modulated',False):
2234            ik = 7
2235            Vec,vRef,maxH = General['SuperVec']
2236            if vRef:
2237                if pFile: pFile.write(' New modulation vector:\n')
2238                namstr = '  names :'
2239                ptstr =  '  values:'
2240                sigstr = '  esds  :'
2241                for iv,var in enumerate(['mV0','mV1','mV2']):
2242                    namstr += '%12s'%(pfx+var)
2243                    ptstr += '%12.6f'%(parmDict[pfx+var])
2244                    if pfx+var in sigDict:
2245                        Vec[iv] = parmDict[pfx+var]
2246                        sigstr += '%12.6f'%(sigDict[pfx+var])
2247                    else:
2248                        sigstr += 12*' '
2249                if pFile: pFile.write(namstr+'\n')
2250                if pFile: pFile.write(ptstr+'\n')
2251                if pFile: pFile.write(sigstr+'\n')
2252           
2253        General['Mass'] = 0.
2254        if Phase['General'].get('doPawley'):
2255            pawleyRef = Phase['Pawley ref']
2256            for i,refl in enumerate(pawleyRef):
2257                key = pfx+'PWLref:'+str(i)
2258                refl[ik] = parmDict[key]
2259                if key in sigDict:
2260                    refl[ik+1] = sigDict[key]
2261                else:
2262                    refl[ik+1] = 0
2263        else:
2264            VRBIds = RBIds['Vector']
2265            RRBIds = RBIds['Residue']
2266            RBModels = Phase['RBModels']
2267            if pFile:
2268                for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2269                    jrb = VRBIds.index(RBObj['RBId'])
2270                    rbsx = str(irb)+':'+str(jrb)
2271                    pFile.write(' Vector rigid body parameters:\n')
2272                    PrintRBObjPOAndSig('RBV',rbsx)
2273                    PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2274                for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2275                    jrb = RRBIds.index(RBObj['RBId'])
2276                    rbsx = str(irb)+':'+str(jrb)
2277                    pFile.write(' Residue rigid body parameters:\n')
2278                    PrintRBObjPOAndSig('RBR',rbsx)
2279                    PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2280                    PrintRBObjTorAndSig(rbsx)
2281            atomsSig = {}
2282            wavesSig = {}
2283            cx,ct,cs,cia = General['AtomPtrs']
2284            for i,at in enumerate(Atoms):
2285                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2286                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2287                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2288                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2289                for ind in range(cx,cx+4):
2290                    at[ind] = parmDict[names[ind]]
2291                    if ind in range(cx,cx+3):
2292                        name = names[ind].replace('A','dA')
2293                    else:
2294                        name = names[ind]
2295                    if name in sigDict:
2296                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2297                if at[cia] == 'I':
2298                    at[cia+1] = parmDict[names[cia+1]]
2299                    if names[cia+1] in sigDict:
2300                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2301                else:
2302                    for ind in range(cia+2,cia+8):
2303                        at[ind] = parmDict[names[ind]]
2304                        if names[ind] in sigDict:
2305                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2306                if General['Type'] == 'magnetic':
2307                    for ind in range(cx+4,cx+7):
2308                        at[ind] = parmDict[names[ind]]
2309                        if names[ind] in sigDict:
2310                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2311                ind = General['AtomTypes'].index(at[ct])
2312                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2313                if General.get('Modulated',False):
2314                    AtomSS = at[-1]['SS1']
2315                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2316                        Waves = AtomSS[Stype]
2317                        if len(Waves):
2318                            waveType = Waves[0]
2319                        else:
2320                            continue
2321                        for iw,wave in enumerate(Waves[1:]):
2322                            stiw = str(i)+':'+str(iw)
2323                            if Stype == 'Spos':
2324                                if waveType in ['ZigZag','Block',] and not iw:
2325                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2326                                else:
2327                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2328                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2329                            elif Stype == 'Sadp':
2330                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2331                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2332                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2333                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2334                            elif Stype == 'Sfrac':
2335                                if 'Crenel' in waveType and not iw:
2336                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2337                                else:
2338                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2339                            elif Stype == 'Smag':
2340                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2341                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2342                            for iname,name in enumerate(names):
2343                                AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2344                                if pfx+name in sigDict:
2345                                    wavesSig[name] = sigDict[pfx+name]
2346
2347            if pFile: PrintAtomsAndSig(General,Atoms,atomsSig)
2348            if pFile and General['Type'] == 'magnetic':
2349                PrintMomentsAndSig(General,Atoms,atomsSig)
2350            if pFile and General.get('Modulated',False):
2351                PrintWavesAndSig(General,Atoms,wavesSig)
2352               
2353            density = G2mth.getDensity(General)[0]
2354            if pFile: pFile.write('\n Density: {:.4f} g/cm**3\n'.format(density))
2355           
2356       
2357        textureData = General['SH Texture']   
2358        if textureData['Order']:
2359            SHtextureSig = {}
2360            for name in ['omega','chi','phi']:
2361                aname = pfx+'SH '+name
2362                textureData['Sample '+name][1] = parmDict[aname]
2363                if aname in sigDict:
2364                    SHtextureSig['Sample '+name] = sigDict[aname]
2365            for name in textureData['SH Coeff'][1]:
2366                aname = pfx+name
2367                textureData['SH Coeff'][1][name] = parmDict[aname]
2368                if aname in sigDict:
2369                    SHtextureSig[name] = sigDict[aname]
2370            PrintSHtextureAndSig(textureData,SHtextureSig)
2371        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2372            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2373                textureData,RestraintDict[phase],pFile)
2374                   
2375################################################################################
2376##### Histogram & Phase data
2377################################################################################       
2378                   
2379def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2380    '''Loads the HAP histogram/phase information into dicts
2381
2382    :param dict Phases: phase information
2383    :param dict Histograms: Histogram information
2384    :param bool Print: prints information as it is read
2385    :param file pFile: file object to print to (the default, None causes printing to the console)
2386    :param bool resetRefList: Should the contents of the Reflection List be initialized
2387      on loading. The default, True, initializes the Reflection List as it is loaded.
2388
2389    :returns: (hapVary,hapDict,controlDict)
2390      * hapVary: list of refined variables
2391      * hapDict: dict with refined variables and their values
2392      * controlDict: dict with fixed parameters
2393    '''
2394   
2395    def PrintSize(hapData):
2396        if hapData[0] in ['isotropic','uniaxial']:
2397            line = '\n Size model    : %9s'%(hapData[0])
2398            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2399            if hapData[0] == 'uniaxial':
2400                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2401            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2402            pFile.write(line+'\n')
2403        else:
2404            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2405                (hapData[0],hapData[1][2],hapData[2][2]))
2406            Snames = ['S11','S22','S33','S12','S13','S23']
2407            ptlbls = ' names :'
2408            ptstr =  ' values:'
2409            varstr = ' refine:'
2410            for i,name in enumerate(Snames):
2411                ptlbls += '%12s' % (name)
2412                ptstr += '%12.3f' % (hapData[4][i])
2413                varstr += '%12s' % (str(hapData[5][i]))
2414            pFile.write(ptlbls+'\n')
2415            pFile.write(ptstr+'\n')
2416            pFile.write(varstr+'\n')
2417       
2418    def PrintMuStrain(hapData,SGData):
2419        if hapData[0] in ['isotropic','uniaxial']:
2420            line = '\n Mustrain model: %9s'%(hapData[0])
2421            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2422            if hapData[0] == 'uniaxial':
2423                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2424            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2425            pFile.write(line+'\n')
2426        else:
2427            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2428                (hapData[0],hapData[1][2],hapData[2][2]))
2429            Snames = G2spc.MustrainNames(SGData)
2430            ptlbls = ' names :'
2431            ptstr =  ' values:'
2432            varstr = ' refine:'
2433            for i,name in enumerate(Snames):
2434                ptlbls += '%12s' % (name)
2435                ptstr += '%12.1f' % (hapData[4][i])
2436                varstr += '%12s' % (str(hapData[5][i]))
2437            pFile.write(ptlbls+'\n')
2438            pFile.write(ptstr+'\n')
2439            pFile.write(varstr+'\n')
2440
2441    def PrintHStrain(hapData,SGData):
2442        pFile.write('\n Hydrostatic/elastic strain:\n')
2443        Hsnames = G2spc.HStrainNames(SGData)
2444        ptlbls = ' names :'
2445        ptstr =  ' values:'
2446        varstr = ' refine:'
2447        for i,name in enumerate(Hsnames):
2448            ptlbls += '%12s' % (name)
2449            ptstr += '%12.4g' % (hapData[0][i])
2450            varstr += '%12s' % (str(hapData[1][i]))
2451        pFile.write(ptlbls+'\n')
2452        pFile.write(ptstr+'\n')
2453        pFile.write(varstr+'\n')
2454
2455    def PrintSHPO(hapData):
2456        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2457        ptlbls = ' names :'
2458        ptstr =  ' values:'
2459        for item in hapData[5]:
2460            ptlbls += '%12s'%(item)
2461            ptstr += '%12.3f'%(hapData[5][item]) 
2462        pFile.write(ptlbls+'\n')
2463        pFile.write(ptstr+'\n')
2464   
2465    def PrintBabinet(hapData):
2466        pFile.write('\n Babinet form factor modification:\n')
2467        ptlbls = ' names :'
2468        ptstr =  ' values:'
2469        varstr = ' refine:'
2470        for name in ['BabA','BabU']:
2471            ptlbls += '%12s' % (name)
2472            ptstr += '%12.6f' % (hapData[name][0])
2473            varstr += '%12s' % (str(hapData[name][1]))
2474        pFile.write(ptlbls+'\n')
2475        pFile.write(ptstr+'\n')
2476        pFile.write(varstr+'\n')
2477       
2478    hapDict = {}
2479    hapVary = []
2480    controlDict = {}
2481   
2482    for phase in Phases:
2483        HistoPhase = Phases[phase]['Histograms']
2484        SGData = Phases[phase]['General']['SGData']
2485        cell = Phases[phase]['General']['Cell'][1:7]
2486        A = G2lat.cell2A(cell)
2487        if Phases[phase]['General'].get('Modulated',False):
2488            SSGData = Phases[phase]['General']['SSGData']
2489            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2490        pId = Phases[phase]['pId']
2491        for histogram in Histograms:
2492            if histogram not in HistoPhase and phase in Histograms[histogram]['Reflection Lists']:
2493                #remove previously created reflection list if histogram is removed from phase
2494                #print("removing ",phase,"from",histogram)
2495                del Histograms[histogram]['Reflection Lists'][phase]
2496        histoList = list(HistoPhase.keys())
2497        histoList.sort()
2498        for histogram in histoList:
2499            try:
2500                Histogram = Histograms[histogram]
2501            except KeyError:                       
2502                #skip if histogram not included e.g. in a sequential refinement
2503                continue
2504            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2505                if phase in Histograms[histogram]['Reflection Lists']:
2506                    del Histograms[histogram]['Reflection Lists'][phase]
2507                continue
2508            hapData = HistoPhase[histogram]
2509            hId = Histogram['hId']
2510            if 'PWDR' in histogram:
2511                limits = Histogram['Limits'][1]
2512                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2513                if 'C' in inst['Type'][1]:
2514                    try:
2515                        wave = inst['Lam'][1]
2516                    except KeyError:
2517                        wave = inst['Lam1'][1]
2518                    dmin = wave/(2.0*sind(limits[1]/2.0))
2519                elif 'T' in inst['Type'][0]:
2520                    dmin = limits[0]/inst['difC'][1]
2521                else:
2522                    wave = inst['Lam'][1]
2523                    dmin = wave/(2.0*sind(limits[1]/2.0))
2524                pfx = str(pId)+':'+str(hId)+':'
2525                if Phases[phase]['General']['doPawley']:
2526                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2527                    hapDict[pfx+'newLeBail'] = True
2528                    Tmin = G2lat.Dsp2pos(inst,dmin)
2529                    if 'T' in inst['Type'][1]:
2530                        limits[0] = max(limits[0],Tmin)
2531                    else:
2532                        limits[1] = min(limits[1],Tmin)
2533                else:
2534                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2535                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2536                if Phases[phase]['General']['Type'] == 'magnetic':
2537                    dmin = max(dmin,Phases[phase]['General'].get('MagDmin',0.))
2538                for item in ['Scale','Extinction']:
2539                    hapDict[pfx+item] = hapData[item][0]
2540                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2541                        hapVary.append(pfx+item)
2542                names = G2spc.HStrainNames(SGData)
2543                HSvals = []
2544                for i,name in enumerate(names):
2545                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2546                    HSvals.append(hapDict[pfx+name])
2547                    if hapData['HStrain'][1][i]:
2548                        hapVary.append(pfx+name)
2549                if 'Layer Disp' in hapData:
2550                    hapDict[pfx+'LayerDisp'] = hapData['Layer Disp'][0]
2551                    if hapData['Layer Disp'][1]:
2552                        hapVary.append(pfx+'LayerDisp')
2553                else:
2554                    hapDict[pfx+'LayerDisp'] = 0.0
2555                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2556                if hapData['Pref.Ori.'][0] == 'MD':
2557                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2558                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2559                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2560                        hapVary.append(pfx+'MD')
2561                else:                           #'SH' spherical harmonics
2562                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2563                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2564                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2565                    controlDict[pfx+'SHhkl'] = []
2566                    try: #patch for old Pref.Ori. items
2567                        controlDict[pfx+'SHtoler'] = 0.1
2568                        if hapData['Pref.Ori.'][6][0] != '':
2569                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2570                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2571                    except IndexError:
2572                        pass
2573                    for item in hapData['Pref.Ori.'][5]:
2574                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2575                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2576                            hapVary.append(pfx+item)
2577                for item in ['Mustrain','Size']:
2578                    controlDict[pfx+item+'Type'] = hapData[item][0]
2579                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2580                    if hapData[item][2][2]:
2581                        hapVary.append(pfx+item+';mx')
2582                    if hapData[item][0] in ['isotropic','uniaxial']:
2583                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2584                        if hapData[item][2][0]:
2585                            hapVary.append(pfx+item+';i')
2586                        if hapData[item][0] == 'uniaxial':
2587                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2588                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2589                            if hapData[item][2][1]:
2590                                hapVary.append(pfx+item+';a')
2591                    else:       #generalized for mustrain or ellipsoidal for size
2592                        Nterms = len(hapData[item][4])
2593                        if item == 'Mustrain':
2594                            names = G2spc.MustrainNames(SGData)
2595                            pwrs = []
2596                            for name in names:
2597                                h,k,l = name[1:]
2598                                pwrs.append([int(h),int(k),int(l)])
2599                            controlDict[pfx+'MuPwrs'] = pwrs
2600                        for i in range(Nterms):
2601                            sfx = ';'+str(i)
2602                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2603                            if hapData[item][5][i]:
2604                                hapVary.append(pfx+item+sfx)
2605                if Phases[phase]['General']['Type'] != 'magnetic':
2606                    for bab in ['BabA','BabU']:
2607                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2608                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2609                            hapVary.append(pfx+bab)
2610                               
2611                if Print: 
2612                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2613                    pFile.write(135*'='+'\n')
2614                    if hapDict[pfx+'LeBail']:
2615                        pFile.write(' Perform LeBail extraction\n')                     
2616                    else:
2617                        pFile.write(' Phase fraction  : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2618                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2619                        if hapData['Pref.Ori.'][0] == 'MD':
2620                            Ax = hapData['Pref.Ori.'][3]
2621                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2622                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2623                        else: #'SH' for spherical harmonics
2624                            PrintSHPO(hapData['Pref.Ori.'])
2625                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2626                    PrintSize(hapData['Size'])
2627                    PrintMuStrain(hapData['Mustrain'],SGData)
2628                    PrintHStrain(hapData['HStrain'],SGData)
2629                    if 'Layer Disp' in hapData:
2630                        pFile.write(' Layer Displacement: %10.3f Refine? %s\n'%(hapData['Layer Disp'][0],hapData['Layer Disp'][1]))
2631                    if Phases[phase]['General']['Type'] != 'magnetic':
2632                        if hapData['Babinet']['BabA'][0]:
2633                            PrintBabinet(hapData['Babinet'])
2634                if phase in Histogram['Reflection Lists'] and 'RefList' not in Histogram['Reflection Lists'][phase] and hapData.get('LeBail',False):
2635                    hapData['newLeBail'] = True
2636                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData.get('LeBail',False) and hapData['newLeBail'])):
2637                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2638                        hapData['newLeBail'] = False
2639                    refList = []
2640#                    Uniq = []
2641#                    Phi = []
2642                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2643                    if Phases[phase]['General'].get('Modulated',False):
2644                        ifSuper = True
2645                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2646                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2647                        for h,k,l,m,d in HKLd:
2648                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2649                            mul *= 2      # for powder overlap of Friedel pairs
2650                            if m or not ext or useExt:
2651                                if 'C' in inst['Type'][0]:
2652                                    pos = G2lat.Dsp2pos(inst,d)
2653                                    if limits[0] < pos < limits[1]:
2654                                        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])
2655                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2656#                                        Uniq.append(uniq)
2657#                                        Phi.append(phi)
2658                                elif 'T' in inst['Type'][0]:
2659                                    pos = G2lat.Dsp2pos(inst,d)
2660                                    if limits[0] < pos < limits[1]:
2661                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2662                                        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])
2663                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2664                                        #TODO - if tabulated put alp & bet in here
2665#                                        Uniq.append(uniq)
2666#                                        Phi.append(phi)
2667                                elif 'B' in inst['Type'][0]:
2668                                    pos = G2lat.Dsp2pos(inst,d)
2669                                    if limits[0] < pos < limits[1]:
2670                                        refList.append([h,k,l,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])
2671                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2672                    else:
2673                        ifSuper = False
2674                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2675                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2676                        for h,k,l,d in HKLd:
2677                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2678                            if ext and 'N' in inst['Type'][0] and  'MagSpGrp' in SGData:
2679                                ext = G2spc.checkMagextc([h,k,l],SGData)
2680                            mul *= 2      # for powder overlap of Friedel pairs
2681                            if ext and not useExt:
2682                                continue
2683                            if 'C' in inst['Type'][0]:
2684                                pos = G2lat.Dsp2pos(inst,d)
2685                                if limits[0] < pos < limits[1]:
2686                                    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])
2687                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2688#                                    Uniq.append(uniq)
2689#                                    Phi.append(phi)
2690                            elif 'T' in inst['Type'][0]:
2691                                pos = G2lat.Dsp2pos(inst,d)
2692                                if limits[0] < pos < limits[1]:
2693                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2694                                    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])
2695                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2696#                                    Uniq.append(uniq)
2697#                                    Phi.append(phi)
2698                            elif 'B' in inst['Type'][0]:
2699                                pos = G2lat.Dsp2pos(inst,d)
2700                                if limits[0] < pos < limits[1]:
2701                                    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])
2702                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet, prfo,abs,ext
2703                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2704            elif 'HKLF' in histogram:
2705                inst = Histogram['Instrument Parameters'][0]
2706                hId = Histogram['hId']
2707                hfx = ':%d:'%(hId)
2708                for item in inst:
2709                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2710                    hapDict[hfx+item] = inst[item][1]
2711                pfx = str(pId)+':'+str(hId)+':'
2712                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2713                if hapData['Scale'][1]:
2714                    hapVary.append(pfx+'Scale')
2715                               
2716                extApprox,extType,extParms = hapData['Extinction']
2717                controlDict[pfx+'EType'] = extType
2718                controlDict[pfx+'EApprox'] = extApprox
2719                if 'C' in inst['Type'][0]:
2720                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2721                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2722                if 'Primary' in extType:
2723                    Ekey = ['Ep',]
2724                elif 'I & II' in extType:
2725                    Ekey = ['Eg','Es']
2726                elif 'Secondary Type II' == extType:
2727                    Ekey = ['Es',]
2728                elif 'Secondary Type I' == extType:
2729                    Ekey = ['Eg',]
2730                else:   #'None'
2731                    Ekey = []
2732                for eKey in Ekey:
2733                    hapDict[pfx+eKey] = extParms[eKey][0]
2734                    if extParms[eKey][1]:
2735                        hapVary.append(pfx+eKey)
2736                for bab in ['BabA','BabU']:
2737                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2738                    if hapData['Babinet'][bab][1]:
2739                        hapVary.append(pfx+bab)
2740                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2741                if len(Twins) == 1:
2742                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2743                    if hapData.get('Flack',[0,False])[1]:
2744                        hapVary.append(pfx+'Flack')
2745                sumTwFr = 0.
2746                controlDict[pfx+'TwinLaw'] = []
2747                controlDict[pfx+'TwinInv'] = []
2748                NTL = 0           
2749                for it,twin in enumerate(Twins):
2750                    if 'bool' in str(type(twin[0])):
2751                        controlDict[pfx+'TwinInv'].append(twin[0])
2752                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2753                    else:
2754                        NTL += 1
2755                        controlDict[pfx+'TwinInv'].append(False)
2756                        controlDict[pfx+'TwinLaw'].append(twin[0])
2757                    if it:
2758                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2759                        sumTwFr += twin[1]
2760                    else:
2761                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2762                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2763                    if Twins[0][1][1]:
2764                        hapVary.append(pfx+'TwinFr:'+str(it))
2765                controlDict[pfx+'NTL'] = NTL
2766                #need to make constraint on TwinFr
2767                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2768                if len(Twins) > 1:    #force sum to unity
2769                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2770                if Print: 
2771                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2772                    pFile.write(135*'='+'\n')
2773                    pFile.write(' Scale factor     : %10.4g Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2774                    if extType != 'None':
2775                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2776                        text = ' Parameters       :'
2777                        for eKey in Ekey:
2778                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2779                        pFile.write(text+'\n')
2780                    if hapData['Babinet']['BabA'][0]:
2781                        PrintBabinet(hapData['Babinet'])
2782                    if not SGData['SGInv'] and len(Twins) == 1:
2783                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2784                    if len(Twins) > 1:
2785                        for it,twin in enumerate(Twins):
2786                            if 'bool' in str(type(twin[0])):
2787                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine? %s\n'%
2788                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it]),str(Twins[0][1][1])))
2789                            else:
2790                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2791                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],str(Twins[0][1][1])))
2792                       
2793                Histogram['Reflection Lists'] = phase       
2794               
2795    return hapVary,hapDict,controlDict
2796   
2797def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,calcControls,Print=True,pFile=None):
2798    'needs a doc string'
2799   
2800    def PrintSizeAndSig(hapData,sizeSig):
2801        line = '\n Size model:     %9s'%(hapData[0])
2802        refine = False
2803        if hapData[0] in ['isotropic','uniaxial']:
2804            line += ' equatorial:%12.5f'%(hapData[1][0])
2805            if sizeSig[0][0]:
2806                line += ', sig:%8.4f'%(sizeSig[0][0])
2807                refine = True
2808            if hapData[0] == 'uniaxial':
2809                line += ' axial:%12.4f'%(hapData[1][1])
2810                if sizeSig[0][1]:
2811                    refine = True
2812                    line += ', sig:%8.4f'%(sizeSig[0][1])
2813            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2814            if sizeSig[0][2]:
2815                refine = True
2816                line += ', sig:%8.4f'%(sizeSig[0][2])
2817            if refine:
2818                pFile.write(line+'\n')
2819        else:
2820            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2821            if sizeSig[0][2]:
2822                refine = True
2823                line += ', sig:%8.4f'%(sizeSig[0][2])
2824            Snames = ['S11','S22','S33','S12','S13','S23']
2825            ptlbls = ' name  :'
2826            ptstr =  ' value :'
2827            sigstr = ' sig   :'
2828            for i,name in enumerate(Snames):
2829                ptlbls += '%12s' % (name)
2830                ptstr += '%12.6f' % (hapData[4][i])
2831                if sizeSig[1][i]:
2832                    refine = True
2833                    sigstr += '%12.6f' % (sizeSig[1][i])
2834                else:
2835                    sigstr += 12*' '
2836            if refine:
2837                pFile.write(line+'\n')
2838                pFile.write(ptlbls+'\n')
2839                pFile.write(ptstr+'\n')
2840                pFile.write(sigstr+'\n')
2841       
2842    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2843        line = '\n Mustrain model: %9s\n'%(hapData[0])
2844        refine = False
2845        if hapData[0] in ['isotropic','uniaxial']:
2846            line += ' equatorial:%12.1f'%(hapData[1][0])
2847            if mustrainSig[0][0]:
2848                line += ', sig:%8.1f'%(mustrainSig[0][0])
2849                refine = True
2850            if hapData[0] == 'uniaxial':
2851                line += ' axial:%12.1f'%(hapData[1][1])
2852                if mustrainSig[0][1]:
2853                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2854            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2855            if mustrainSig[0][2]:
2856                refine = True
2857                line += ', sig:%8.3f'%(mustrainSig[0][2])
2858            if refine:
2859                pFile.write(line+'\n')
2860        else:
2861            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2862            if mustrainSig[0][2]:
2863                refine = True
2864                line += ', sig:%8.3f'%(mustrainSig[0][2])
2865            Snames = G2spc.MustrainNames(SGData)
2866            ptlbls = ' name  :'
2867            ptstr =  ' value :'
2868            sigstr = ' sig   :'
2869            for i,name in enumerate(Snames):
2870                ptlbls += '%12s' % (name)
2871                ptstr += '%12.1f' % (hapData[4][i])
2872                if mustrainSig[1][i]:
2873                    refine = True
2874                    sigstr += '%12.1f' % (mustrainSig[1][i])
2875                else:
2876                    sigstr += 12*' '
2877            if refine:
2878                pFile.write(line+'\n')
2879                pFile.write(ptlbls+'\n')
2880                pFile.write(ptstr+'\n')
2881                pFile.write(sigstr+'\n')
2882           
2883    def PrintHStrainAndSig(hapData,strainSig,SGData):
2884        Hsnames = G2spc.HStrainNames(SGData)
2885        ptlbls = ' name  :'
2886        ptstr =  ' value :'
2887        sigstr = ' sig   :'
2888        refine = False
2889        for i,name in enumerate(Hsnames):
2890            ptlbls += '%12s' % (name)
2891            ptstr += '%12.4g' % (hapData[0][i])
2892            if name in strainSig:
2893                refine = True
2894                sigstr += '%12.4g' % (strainSig[name])
2895            else:
2896                sigstr += 12*' '
2897        if refine:
2898            pFile.write('\n Hydrostatic/elastic strain:\n')
2899            pFile.write(ptlbls+'\n')
2900            pFile.write(ptstr+'\n')
2901            pFile.write(sigstr+'\n')
2902       
2903    def PrintSHPOAndSig(pfx,hapData,POsig):
2904        Tindx = 1.0
2905        Tvar = 0.0
2906        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2907        ptlbls = ' names :'
2908        ptstr =  ' values:'
2909        sigstr = ' sig   :'
2910        for item in hapData[5]:
2911            ptlbls += '%12s'%(item)
2912            ptstr += '%12.3f'%(hapData[5][item])
2913            l = 2.0*eval(item.strip('C'))[0]+1
2914            Tindx += hapData[5][item]**2/l
2915            if pfx+item in POsig:
2916                Tvar += (2.*hapData[5][item]*POsig[pfx+item]/l)**2
2917                sigstr += '%12.3f'%(POsig[pfx+item])
2918            else:
2919                sigstr += 12*' ' 
2920        pFile.write(ptlbls+'\n')
2921        pFile.write(ptstr+'\n')
2922        pFile.write(sigstr+'\n')
2923        pFile.write('\n Texture index J = %.3f(%d)\n'%(Tindx,int(1000*np.sqrt(Tvar))))
2924       
2925    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2926        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2927        text = ''
2928        for item in hapData[2]:
2929            if pfx+item in ScalExtSig:
2930                text += '       %s: '%(item)
2931                text += '%12.2e'%(hapData[2][item][0])
2932                if pfx+item in ScalExtSig:
2933                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2934        pFile.write(text+'\n')   
2935       
2936    def PrintBabinetAndSig(pfx,hapData,BabSig):
2937        pFile.write('\n Babinet form factor modification:\n')
2938        ptlbls = ' names :'
2939        ptstr =  ' values:'
2940        sigstr = ' sig   :'
2941        for item in hapData:
2942            ptlbls += '%12s'%(item)
2943            ptstr += '%12.3f'%(hapData[item][0])
2944            if pfx+item in BabSig:
2945                sigstr += '%12.3f'%(BabSig[pfx+item])
2946            else:
2947                sigstr += 12*' ' 
2948        pFile.write(ptlbls+'\n')
2949        pFile.write(ptstr+'\n')
2950        pFile.write(sigstr+'\n')
2951       
2952    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2953        pFile.write('\n Twin Law fractions :\n')
2954        ptlbls = ' names :'
2955        ptstr =  ' values:'
2956        sigstr = ' sig   :'
2957        for it,item in enumerate(twinData):
2958            ptlbls += '     twin #%d'%(it)
2959            if it:
2960                ptstr += '%12.3f'%(item[1])
2961            else:
2962                ptstr += '%12.3f'%(item[1][0])
2963            if pfx+'TwinFr:'+str(it) in TwinSig:
2964                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2965            else:
2966                sigstr += 12*' ' 
2967        pFile.write(ptlbls+'\n')
2968        pFile.write(ptstr+'\n')
2969        pFile.write(sigstr+'\n')
2970       
2971   
2972    PhFrExtPOSig = {}
2973    SizeMuStrSig = {}
2974    ScalExtSig = {}
2975    BabSig = {}
2976    TwinFrSig = {}
2977    wtFrSum = {}
2978    for phase in Phases:
2979        HistoPhase = Phases[phase]['Histograms']
2980        General = Phases[phase]['General']
2981        SGData = General['SGData']
2982        pId = Phases[phase]['pId']
2983        histoList = list(HistoPhase.keys())
2984        histoList.sort()
2985        for histogram in histoList:
2986            try:
2987                Histogram = Histograms[histogram]
2988            except KeyError:                       
2989                #skip if histogram not included e.g. in a sequential refinement
2990                continue
2991            if not Phases[phase]['Histograms'][histogram]['Use']:
2992                #skip if phase absent from this histogram
2993                continue
2994            hapData = HistoPhase[histogram]
2995            hId = Histogram['hId']
2996            pfx = str(pId)+':'+str(hId)+':'
2997            if hId not in wtFrSum:
2998                wtFrSum[hId] = 0.
2999            if 'PWDR' in histogram:
3000                for item in ['Scale','Extinction']:
3001                    hapData[item][0] = parmDict[pfx+item]
3002                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
3003                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
3004                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
3005                if hapData['Pref.Ori.'][0] == 'MD':
3006                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
3007                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
3008                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
3009                else:                           #'SH' spherical harmonics
3010                    for item in hapData['Pref.Ori.'][5]:
3011                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
3012                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
3013                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
3014                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
3015                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
3016                    pfx+'HStrain':{}})                 
3017                for item in ['Mustrain','Size']:
3018                    hapData[item][1][2] = parmDict[pfx+item+';mx']
3019#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
3020                    if pfx+item+';mx' in sigDict:
3021                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
3022                    if hapData[item][0] in ['isotropic','uniaxial']:                   
3023                        hapData[item][1][0] = parmDict[pfx+item+';i']
3024                        if item == 'Size':
3025                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
3026                        if pfx+item+';i' in sigDict: 
3027                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
3028                        if hapData[item][0] == 'uniaxial':
3029                            hapData[item][1][1] = parmDict[pfx+item+';a']
3030                            if item == 'Size':
3031                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
3032                            if pfx+item+';a' in sigDict:
3033                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
3034                    else:       #generalized for mustrain or ellipsoidal for size
3035                        Nterms = len(hapData[item][4])
3036                        for i in range(Nterms):
3037                            sfx = ';'+str(i)
3038                            hapData[item][4][i] = parmDict[pfx+item+sfx]
3039                            if pfx+item+sfx in sigDict:
3040                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
3041                names = G2spc.HStrainNames(SGData)
3042                for i,name in enumerate(names):
3043                    hapData['HStrain'][0][i] = parmDict[pfx+name]
3044                    if pfx+name in sigDict:
3045                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
3046                if 'Layer Disp' in hapData:
3047                    hapData['Layer Disp'][0] = parmDict[pfx+'LayerDisp']
3048                    if pfx+'LayerDisp' in sigDict:
3049                        SizeMuStrSig[pfx+'LayerDisp'] = sigDict[pfx+'LayerDisp']
3050                if Phases[phase]['General']['Type'] != 'magnetic':
3051                    for name in ['BabA','BabU']:
3052                        hapData['Babinet'][name][0] = parmDict[pfx+name]
3053                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
3054                            BabSig[pfx+name] = sigDict[pfx+name]               
3055               
3056            elif 'HKLF' in histogram:
3057                for item in ['Scale','Flack']:
3058                    if parmDict.get(pfx+item):
3059                        hapData[item][0] = parmDict[pfx+item]
3060                        if pfx+item in sigDict:
3061                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3062                for item in ['Ep','Eg','Es']:
3063                    if parmDict.get(pfx+item):
3064                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
3065                        if pfx+item in sigDict:
3066                            ScalExtSig[pfx+item] = sigDict[pfx+item]
3067                for item in ['BabA','BabU']:
3068                    hapData['Babinet'][item][0] = parmDict[pfx+item]
3069                    if pfx+item in sigDict:
3070                        BabSig[pfx+item] = sigDict[pfx+item]
3071                if 'Twins' in hapData:
3072                    it = 1
3073                    sumTwFr = 0.
3074                    while True:
3075                        try:
3076                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
3077                            if pfx+'TwinFr:'+str(it) in sigDict:
3078                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
3079                            if it:
3080                                sumTwFr += hapData['Twins'][it][1]
3081                            it += 1
3082                        except KeyError:
3083                            break
3084                    hapData['Twins'][0][1][0] = 1.-sumTwFr
3085
3086    if Print:
3087        for phase in Phases:
3088            HistoPhase = Phases[phase]['Histograms']
3089            General = Phases[phase]['General']
3090            SGData = General['SGData']
3091            pId = Phases[phase]['pId']
3092            histoList = list(HistoPhase.keys())
3093            histoList.sort()
3094            for histogram in histoList:
3095                try:
3096                    Histogram = Histograms[histogram]
3097                except KeyError:                       
3098                    #skip if histogram not included e.g. in a sequential refinement
3099                    continue
3100                hapData = HistoPhase[histogram]
3101                hId = Histogram['hId']
3102                Histogram['Residuals'][str(pId)+'::Name'] = phase
3103                pfx = str(pId)+':'+str(hId)+':'
3104                hfx = ':%s:'%(hId)
3105                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
3106                    continue
3107                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
3108                pFile.write(135*'='+'\n')
3109                Inst = Histogram['Instrument Parameters'][0]
3110                if 'PWDR' in histogram:
3111                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
3112                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
3113                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
3114                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
3115                   
3116                    if parmDict[pfx+'LeBail']:
3117                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
3118                    else:
3119                        # if calcControls != None:    #skipped in seqRefine
3120                        #     if 'X'in Inst['Type'][0]:
3121                        #         PrintFprime(calcControls['FFtables'],hfx,pFile)
3122                        #     elif 'NC' in Inst['Type'][0]:
3123                        #         PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3124                        if pfx+'Scale' in PhFrExtPOSig:
3125                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
3126                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
3127                            pFile.write(' Phase fraction  : %10.5g, sig %10.5g Weight fraction  : %8.5f, sig %10.5f\n'%
3128                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
3129                        if pfx+'Extinction' in PhFrExtPOSig:
3130                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
3131                        if hapData['Pref.Ori.'][0] == 'MD':
3132                            if pfx+'MD' in PhFrExtPOSig:
3133                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
3134                        else:
3135                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
3136                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
3137                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
3138                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
3139                    if pfx+'LayerDisp' in SizeMuStrSig:
3140                        pFile.write(' Layer displacement : %10.3f, sig %10.3f\n'%(hapData['Layer Disp'][0],SizeMuStrSig[pfx+'LayerDisp']))           
3141                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
3142                        if len(BabSig):
3143                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3144                   
3145                elif 'HKLF' in histogram:
3146                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
3147                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
3148                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
3149                    if calcControls != None:    #skipped in seqRefine
3150                        if 'X'in Inst['Type'][0]:
3151                            PrintFprime(calcControls['FFtables'],hfx,pFile)
3152                        elif 'NC' in Inst['Type'][0]:
3153                            PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3154                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
3155                    if pfx+'Scale' in ScalExtSig:
3156                        pFile.write(' Scale factor : %10.4g, sig %10.4g\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
3157                    if hapData['Extinction'][0] != 'None':
3158                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
3159                    if len(BabSig):
3160                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
3161                    if pfx+'Flack' in ScalExtSig:
3162                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
3163                    if len(TwinFrSig):
3164                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
3165
3166################################################################################
3167##### Histogram data
3168################################################################################       
3169                   
3170def GetHistogramData(Histograms,Print=True,pFile=None):
3171    'needs a doc string'
3172   
3173    def GetBackgroundParms(hId,Background):
3174        Back = Background[0]
3175        DebyePeaks = Background[1]
3176        bakType,bakFlag = Back[:2]
3177        backVals = Back[3:]
3178        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
3179        backDict = dict(zip(backNames,backVals))
3180        backVary = []
3181        if bakFlag:
3182            backVary = backNames
3183        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
3184        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
3185        debyeDict = {}
3186        debyeList = []
3187        for i in range(DebyePeaks['nDebye']):
3188            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
3189            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
3190            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
3191        debyeVary = []
3192        for item in debyeList:
3193            if item[1]:
3194                debyeVary.append(item[0])
3195        backDict.update(debyeDict)
3196        backVary += debyeVary
3197        peakDict = {}
3198        peakList = []
3199        for i in range(DebyePeaks['nPeaks']):
3200            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
3201                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
3202            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
3203            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
3204        peakVary = []
3205        for item in peakList:
3206            if item[1]:
3207                peakVary.append(item[0])
3208        backDict.update(peakDict)
3209        backVary += peakVary
3210        if 'background PWDR' in Background[1]:
3211            backDict[':'+str(hId)+':Back File'] = Background[1]['background PWDR'][0]
3212            backDict[':'+str(hId)+':BF mult'] = Background[1]['background PWDR'][1]
3213            try:
3214                if Background[1]['background PWDR'][0] and Background[1]['background PWDR'][2]:
3215                    backVary.append(':'+str(hId)+':BF mult')
3216            except IndexError:  # old version without refine flag
3217                pass
3218        return bakType,backDict,backVary           
3219       
3220    def GetInstParms(hId,Inst):
3221        #patch
3222        if 'Z' not in Inst:
3223            Inst['Z'] = [0.0,0.0,False]
3224        dataType = Inst['Type'][0]
3225        instDict = {}
3226        insVary = []
3227        pfx = ':'+str(hId)+':'
3228        insKeys = list(Inst.keys())
3229        insKeys.sort()
3230        for item in insKeys:
3231            if item in 'Azimuth':
3232                continue
3233            insName = pfx+item
3234            instDict[insName] = Inst[item][1]
3235            if len(Inst[item]) > 2 and Inst[item][2]:
3236                insVary.append(insName)
3237        if 'C' in dataType:
3238            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
3239        elif 'T' in dataType:   #trap zero alp, bet coeff.
3240            if not instDict[pfx+'alpha']:
3241                instDict[pfx+'alpha'] = 1.0
3242            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3243                instDict[pfx+'beta-1'] = 1.0
3244        elif 'B' in dataType:   #trap zero alp, bet coeff.
3245            if not instDict[pfx+'alpha-0'] and not instDict[pfx+'alpha-1']:
3246                instDict[pfx+'alpha-1'] = 1.0
3247            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
3248                instDict[pfx+'beta-1'] = 1.0
3249        return dataType,instDict,insVary
3250       
3251    def GetSampleParms(hId,Sample):
3252        sampVary = []
3253        hfx = ':'+str(hId)+':'       
3254        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
3255            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi'],hfx+'Azimuth':Sample['Azimuth']}
3256        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
3257            if key in Sample:
3258                sampDict[hfx+key] = Sample[key]
3259        Type = Sample['Type']
3260        if 'Bragg' in Type:             #Bragg-Brentano
3261            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3262                sampDict[hfx+item] = Sample[item][0]
3263                if Sample[item][1]:
3264                    sampVary.append(hfx+item)
3265        elif 'Debye' in Type:        #Debye-Scherrer
3266            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3267                sampDict[hfx+item] = Sample[item][0]
3268                if Sample[item][1]:
3269                    sampVary.append(hfx+item)
3270        return Type,sampDict,sampVary
3271       
3272    def PrintBackground(Background):
3273        Back = Background[0]
3274        DebyePeaks = Background[1]
3275        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
3276        line = ' Coefficients: '
3277        for i,back in enumerate(Back[3:]):
3278            line += '%10.3f'%(back)
3279            if i and not i%10:
3280                line += '\n'+15*' '
3281        pFile.write(line+'\n')
3282        if DebyePeaks['nDebye']:
3283            pFile.write('\n Debye diffuse scattering coefficients\n')
3284            parms = ['DebyeA','DebyeR','DebyeU']
3285            line = ' names :  '
3286            for parm in parms:
3287                line += '%8s refine?'%(parm)
3288            pFile.write(line+'\n')
3289            for j,term in enumerate(DebyePeaks['debyeTerms']):
3290                line = ' term'+'%2d'%(j)+':'
3291                for i in range(3):
3292                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3293                pFile.write(line+'\n')
3294        if DebyePeaks['nPeaks']:
3295            pFile.write('\n Single peak coefficients\n')
3296            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3297            line = ' names :  '
3298            for parm in parms:
3299                line += '%8s refine?'%(parm)
3300            pFile.write(line+'\n')
3301            for j,term in enumerate(DebyePeaks['peaksList']):
3302                line = ' peak'+'%2d'%(j)+':'
3303                for i in range(4):
3304                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
3305                pFile.write(line+'\n')
3306        if 'background PWDR' in DebyePeaks:
3307            try:
3308                pFile.write(' Fixed background file: %s; mult: %.3f  Refine? %s\n'%(DebyePeaks['background PWDR'][0],
3309                    DebyePeaks['background PWDR'][1],DebyePeaks['background PWDR'][2]))
3310            except IndexError:  #old version without refine flag
3311                pass
3312       
3313    def PrintInstParms(Inst):
3314        pFile.write('\n Instrument Parameters:\n')
3315        insKeys = list(Inst.keys())
3316        insKeys.sort()
3317        iBeg = 0
3318        Ok = True
3319        while Ok:
3320            ptlbls = ' name  :'
3321            ptstr =  ' value :'
3322            varstr = ' refine:'
3323            iFin = min(iBeg+9,len(insKeys))
3324            for item in insKeys[iBeg:iFin]:
3325                if item not in ['Type','Source']:
3326                    ptlbls += '%12s' % (item)
3327                    ptstr += '%12.6f' % (Inst[item][1])
3328                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
3329                        varstr += 12*' '
3330                    else:
3331                        varstr += '%12s' % (str(bool(Inst[item][2])))
3332            pFile.write(ptlbls+'\n')
3333            pFile.write(ptstr+'\n')
3334            pFile.write(varstr+'\n')
3335            iBeg = iFin
3336            if iBeg == len(insKeys):
3337                Ok = False
3338            else:
3339                pFile.write('\n')
3340       
3341    def PrintSampleParms(Sample):
3342        pFile.write('\n Sample Parameters:\n')
3343        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3344            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3345        ptlbls = ' name  :'
3346        ptstr =  ' value :'
3347        varstr = ' refine:'
3348        if 'Bragg' in Sample['Type']:
3349            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3350                ptlbls += '%14s'%(item)
3351                ptstr += '%14.4f'%(Sample[item][0])
3352                varstr += '%14s'%(str(bool(Sample[item][1])))
3353           
3354        elif 'Debye' in Type:        #Debye-Scherrer
3355            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3356                ptlbls += '%14s'%(item)
3357                ptstr += '%14.4f'%(Sample[item][0])
3358                varstr += '%14s'%(str(bool(Sample[item][1])))
3359
3360        pFile.write(ptlbls+'\n')
3361        pFile.write(ptstr+'\n')
3362        pFile.write(varstr+'\n')
3363       
3364    histDict = {}
3365    histVary = []
3366    controlDict = {}
3367    histoList = list(Histograms.keys())
3368    histoList.sort()
3369    for histogram in histoList:
3370        if 'PWDR' in histogram:
3371            Histogram = Histograms[histogram]
3372            hId = Histogram['hId']
3373            pfx = ':'+str(hId)+':'
3374            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3375            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3376            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3377            for excl in controlDict[pfx+'Exclude']:
3378                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3379            if controlDict[pfx+'Exclude']:
3380                ma.mask_rows(Histogram['Data'])
3381            Background = Histogram['Background']
3382            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3383            controlDict[pfx+'bakType'] = Type
3384            histDict.update(bakDict)
3385            histVary += bakVary
3386           
3387            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3388            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3389                G2fil.G2Print ('Warning: tabulated profile coefficients are ignored')
3390            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3391            controlDict[pfx+'histType'] = Type
3392            if 'XC' in Type or 'XB' in Type:
3393                if pfx+'Lam1' in instDict:
3394                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3395                else:
3396                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3397            histDict.update(instDict)
3398            histVary += insVary
3399           
3400            Sample = Histogram['Sample Parameters']
3401            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3402            controlDict[pfx+'instType'] = Type
3403            histDict.update(sampDict)
3404            histVary += sampVary
3405           
3406   
3407            if Print: 
3408                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3409                pFile.write(135*'='+'\n')
3410                Units = {'C':' deg','T':' msec','B':' deg'}
3411                units = Units[controlDict[pfx+'histType'][2]]
3412                Limits = controlDict[pfx+'Limits']
3413                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3414                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3415                if len(controlDict[pfx+'Exclude']):
3416                    excls = controlDict[pfx+'Exclude']
3417                    for excl in excls:
3418                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3419                PrintSampleParms(Sample)
3420                PrintInstParms(Inst[0])
3421                PrintBackground(Background)
3422        elif 'HKLF' in histogram:
3423            Histogram = Histograms[histogram]
3424            hId = Histogram['hId']
3425            pfx = ':'+str(hId)+':'
3426            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3427            Inst = Histogram['Instrument Parameters'][0]
3428            controlDict[pfx+'histType'] = Inst['Type'][0]
3429            if 'X' in Inst['Type'][0]:
3430                histDict[pfx+'Lam'] = Inst['Lam'][1]
3431                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']
3432            elif 'NC' in Inst['Type'][0] or 'NB' in Inst['Type'][0]:                   
3433                histDict[pfx+'Lam'] = Inst['Lam'][1]
3434    return histVary,histDict,controlDict
3435   
3436def SetHistogramData(parmDict,sigDict,Histograms,calcControls,Print=True,pFile=None,seq=False):
3437    'Shows histogram data after a refinement'
3438   
3439    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3440        Back = Background[0]
3441        DebyePeaks = Background[1]
3442        lenBack = len(Back[3:])
3443        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3444        for i in range(lenBack):
3445            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3446            if pfx+'Back;'+str(i) in sigDict:
3447                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3448        if DebyePeaks['nDebye']:
3449            for i in range(DebyePeaks['nDebye']):
3450                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3451                for j,name in enumerate(names):
3452                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3453                    if name in sigDict:
3454                        backSig[lenBack+3*i+j] = sigDict[name]           
3455        if DebyePeaks['nPeaks']:
3456            for i in range(DebyePeaks['nPeaks']):
3457                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3458                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3459                for j,name in enumerate(names):
3460                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3461                    if name in sigDict:
3462                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3463        if pfx+'BF mult' in sigDict:
3464            DebyePeaks['background PWDR'][1] = parmDict[pfx+'BF mult']
3465            backSig.append(sigDict[pfx+'BF mult'])
3466               
3467        return backSig
3468       
3469    def SetInstParms(pfx,Inst,parmDict,sigDict):
3470        instSig = {}
3471        insKeys = list(Inst.keys())
3472        insKeys.sort()
3473        for item in insKeys:
3474            insName = pfx+item
3475            Inst[item][1] = parmDict[insName]
3476            if insName in sigDict:
3477                instSig[item] = sigDict[insName]
3478            else:
3479                instSig[item] = 0
3480        return instSig
3481       
3482    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3483        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3484            sampSig = [0 for i in range(5)]
3485            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3486                Sample[item][0] = parmDict[pfx+item]
3487                if pfx+item in sigDict:
3488                    sampSig[i] = sigDict[pfx+item]
3489        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3490            sampSig = [0 for i in range(4)]
3491            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3492                Sample[item][0] = parmDict[pfx+item]
3493                if pfx+item in sigDict:
3494                    sampSig[i] = sigDict[pfx+item]
3495        return sampSig
3496       
3497    def PrintBackgroundSig(Background,backSig):
3498        Back = Background[0]
3499        DebyePeaks = Background[1]
3500        valstr = ' value : '
3501        sigstr = ' sig   : '
3502        refine = False
3503        for i,back in enumerate(Back[3:]):
3504            valstr += '%10.4g'%(back)
3505            if Back[1]:
3506                refine = True
3507                sigstr += '%10.4g'%(backSig[i])
3508            else:
3509                sigstr += 10*' '
3510        if refine:
3511            pFile.write('\n Background function: %s\n'%Back[0])
3512            pFile.write(valstr+'\n')
3513            pFile.write(sigstr+'\n')
3514        if DebyePeaks['nDebye']:
3515            ifAny = False
3516            ptfmt = "%12.3f"
3517            names =  ' names :'
3518            ptstr =  ' values:'
3519            sigstr = ' esds  :'
3520            for item in sigDict:
3521                if 'Debye' in item:
3522                    ifAny = True
3523                    names += '%12s'%(item)
3524                    ptstr += ptfmt%(parmDict[item])
3525                    sigstr += ptfmt%(sigDict[item])
3526            if ifAny:
3527                pFile.write('\n Debye diffuse scattering coefficients\n')
3528                pFile.write(names+'\n')
3529                pFile.write(ptstr+'\n')
3530                pFile.write(sigstr+'\n')
3531        if DebyePeaks['nPeaks']:
3532            pFile.write('\n Single peak coefficients:\n')
3533            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3534            line = ' peak no. '
3535            for parm in parms:
3536                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3537            pFile.write(line+'\n')
3538            for ip in range(DebyePeaks['nPeaks']):
3539                ptstr = ' %4d '%(ip)
3540                for parm in parms:
3541                    fmt = '%14.3f'
3542                    efmt = '%12.3f'
3543                    if parm == 'BkPkpos':
3544                        fmt = '%14.4f'
3545                        efmt = '%12.4f'
3546                    name = pfx+parm+';%d'%(ip)
3547                    ptstr += fmt%(parmDict[name])
3548                    if name in sigDict:
3549                        ptstr += efmt%(sigDict[name])
3550                    else:
3551                        ptstr += 12*' '
3552                pFile.write(ptstr+'\n')
3553        if ('background PWDR' in DebyePeaks and
3554                len(DebyePeaks['background PWDR']) >= 3 and
3555                DebyePeaks['background PWDR'][2]):
3556            pFile.write(' Fixed background scale: %.3f(%d)\n'%(DebyePeaks['background PWDR'][1],int(1000*backSig[-1])))
3557        sumBk = np.array(Histogram['sumBk'])
3558        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3559            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3560       
3561    def PrintInstParmsSig(Inst,instSig):
3562        refine = False
3563        insKeys = list(instSig.keys())
3564        insKeys.sort()
3565        iBeg = 0
3566        Ok = True
3567        while Ok:
3568            ptlbls = ' names :'
3569            ptstr =  ' value :'
3570            sigstr = ' sig   :'
3571            iFin = min(iBeg+9,len(insKeys))
3572            for name in insKeys[iBeg:iFin]:
3573                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3574                    ptlbls += '%12s' % (name)
3575                    ptstr += '%12.6f' % (Inst[name][1])
3576                    if instSig[name]:
3577                        refine = True
3578                        sigstr += '%12.6f' % (instSig[name])
3579                    else:
3580                        sigstr += 12*' '
3581            if refine:
3582                pFile.write('\n Instrument Parameters:\n')
3583                pFile.write(ptlbls+'\n')
3584                pFile.write(ptstr+'\n')
3585                pFile.write(sigstr+'\n')
3586            iBeg = iFin
3587            if iBeg == len(insKeys):
3588                Ok = False
3589       
3590    def PrintSampleParmsSig(Sample,sampleSig):
3591        ptlbls = ' names :'
3592        ptstr =  ' values:'
3593        sigstr = ' sig   :'
3594        refine = False
3595        if 'Bragg' in Sample['Type']:
3596            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3597                ptlbls += '%14s'%(item)
3598                ptstr += '%14.4f'%(Sample[item][0])
3599                if sampleSig[i]:
3600                    refine = True
3601                    sigstr += '%14.4f'%(sampleSig[i])
3602                else:
3603                    sigstr += 14*' '
3604           
3605        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3606            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3607                ptlbls += '%14s'%(item)
3608                ptstr += '%14.4f'%(Sample[item][0])
3609                if sampleSig[i]:
3610                    refine = True
3611                    sigstr += '%14.4f'%(sampleSig[i])
3612                else:
3613                    sigstr += 14*' '
3614
3615        if refine:
3616            pFile.write('\n Sample Parameters:\n')
3617            pFile.write(ptlbls+'\n')
3618            pFile.write(ptstr+'\n')
3619            pFile.write(sigstr+'\n')
3620       
3621    histoList = list(Histograms.keys())
3622    histoList.sort()
3623    for histogram in histoList:
3624        if 'PWDR' in histogram:
3625            Histogram = Histograms[histogram]
3626            hId = Histogram['hId']
3627            pfx = ':'+str(hId)+':'
3628            Background = Histogram['Background']
3629            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3630           
3631            Inst = Histogram['Instrument Parameters'][0]
3632            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3633       
3634            Sample = Histogram['Sample Parameters']
3635            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3636
3637            if Print and not seq:
3638                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3639                pFile.write(135*'='+'\n')
3640            if Print:
3641                pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3642                pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3643                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3644                pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3645                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3646                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3647                if calcControls != None:    #skipped in seqRefine
3648                    if 'X' in Inst['Type'][0]:
3649                        PrintFprime(calcControls['FFtables'],pfx,pFile)
3650                    elif 'NC' in Inst['Type'][0]:
3651                        PrintBlength(calcControls['BLtables'],Inst['Lam'][1],pFile)
3652                PrintSampleParmsSig(Sample,sampSig)
3653                PrintInstParmsSig(Inst,instSig)
3654                PrintBackgroundSig(Background,backSig)
3655               
3656def WriteRBObjPOAndSig(pfx,rbfx,rbsx,parmDict,sigDict):
3657    '''Cribbed version of PrintRBObjPOAndSig but returns lists of strings.
3658    Moved so it can be used in ExportCIF
3659    '''
3660    namstr = '  names :'
3661    valstr = '  values:'
3662    sigstr = '  esds  :'
3663    for i,px in enumerate(['Px:','Py:','Pz:']):
3664        name = pfx+rbfx+px+rbsx
3665        namstr += '%12s'%('Pos '+px[1])
3666        valstr += '%12.5f'%(parmDict[name])
3667        if name in sigDict:
3668            sigstr += '%12.5f'%(sigDict[name])
3669        else:
3670            sigstr += 12*' '
3671    for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
3672        name = pfx+rbfx+po+rbsx
3673        namstr += '%12s'%('Ori '+po[1])
3674        valstr += '%12.5f'%(parmDict[name])
3675        if name in sigDict:
3676            sigstr += '%12.5f'%(sigDict[name])
3677        else:
3678            sigstr += 12*' '
3679    name = pfx+rbfx+'f:'+rbsx
3680    namstr += '%12s'%('Frac')
3681    valstr += '%12.5f'%(parmDict[name])
3682    if name in sigDict:
3683        sigstr += '%12.5f'%(sigDict[name])
3684    else:
3685        sigstr += 12*' '
3686    return (namstr,valstr,sigstr)
3687
3688def WriteRBObjTLSAndSig(pfx,rbfx,rbsx,TLS,parmDict,sigDict):
3689    '''Cribbed version of PrintRBObjTLSAndSig but returns lists of strings.
3690    Moved so it can be used in ExportCIF
3691    '''
3692    out = []
3693    namstr = '  names :'
3694    valstr = '  values:'
3695    sigstr = '  esds  :'
3696    if 'N' not in TLS:
3697        out.append(' Thermal motion:\n')
3698    if 'T' in TLS:
3699        for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
3700            name = pfx+rbfx+pt+rbsx
3701            namstr += '%12s'%(pt[:3])
3702            valstr += '%12.5f'%(parmDict[name])
3703            if name in sigDict:
3704                sigstr += '%12.5f'%(sigDict[name])
3705            else:
3706                sigstr += 12*' '
3707        out.append(namstr+'\n')
3708        out.append(valstr+'\n')
3709        out.append(sigstr+'\n')
3710    if 'L' in TLS:
3711        namstr = '  names :'
3712        valstr = '  values:'
3713        sigstr = '  esds  :'
3714        for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
3715            name = pfx+rbfx+pt+rbsx
3716            namstr += '%12s'%(pt[:3])
3717            valstr += '%12.3f'%(parmDict[name])
3718            if name in sigDict:
3719                sigstr += '%12.3f'%(sigDict[name])
3720            else:
3721                sigstr += 12*' '
3722        out.append(namstr+'\n')
3723        out.append(valstr+'\n')
3724        out.append(sigstr+'\n')
3725    if 'S' in TLS:
3726        namstr = '  names :'
3727        valstr = '  values:'
3728        sigstr = '  esds  :'
3729        for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
3730            name = pfx+rbfx+pt+rbsx
3731            namstr += '%12s'%(pt[:3])
3732            valstr += '%12.4f'%(parmDict[name])
3733            if name in sigDict:
3734                sigstr += '%12.4f'%(sigDict[name])
3735            else:
3736                sigstr += 12*' '
3737        out.append(namstr+'\n')
3738        out.append(valstr+'\n')
3739        out.append(sigstr+'\n')
3740    if 'U' in TLS:
3741        name = pfx+rbfx+'U:'+rbsx
3742        namstr = '  names :'+'%12s'%('Uiso')
3743        valstr = '  values:'+'%12.5f'%(parmDict[name])
3744        if name in sigDict:
3745            sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
3746        else:
3747            sigstr = '  esds  :'+12*' '
3748        out.append(namstr+'\n')
3749        out.append(valstr+'\n')
3750        out.append(sigstr+'\n')
3751    return out
3752
3753def WriteRBObjTorAndSig(pfx,rbsx,parmDict,sigDict,nTors):
3754    '''Cribbed version of PrintRBObjTorAndSig but returns lists of strings.
3755    Moved so it can be used in ExportCIF
3756    '''
3757    out = []
3758    namstr = '  names :'
3759    valstr = '  values:'
3760    sigstr = '  esds  :'
3761    out.append(' Torsions:\n')
3762    for it in range(nTors):
3763        name = pfx+'RBRTr;'+str(it)+':'+rbsx
3764        namstr += '%12s'%('Tor'+str(it))
3765        valstr += '%12.4f'%(parmDict[name])
3766        if name in sigDict:
3767            sigstr += '%12.4f'%(sigDict[name])
3768    out.append(namstr+'\n')
3769    out.append(valstr+'\n')
3770    out.append(sigstr+'\n')
3771    return out
3772
3773def WriteResRBModel(RBModel):
3774    '''Write description of a residue rigid body. Code shifted from
3775    PrintResRBModel to make usable from G2export_CIF
3776    '''
3777    out = []
3778    atNames = RBModel['atNames']
3779    rbRef = RBModel['rbRef']
3780    rbSeq = RBModel['rbSeq']
3781    out.append('    At name       x          y          z\n')
3782    for name,xyz in zip(atNames,RBModel['rbXYZ']):
3783        out.append(%8s %10.4f %10.4f %10.4f\n'%(name,xyz[0],xyz[1],xyz[2]))
3784    out.append('  Orientation defined by: %s -> %s & %s -> %s\n'%
3785        (atNames[rbRef[0]],atNames[rbRef[1]],atNames[rbRef[0]],atNames[rbRef[2]]))
3786    if rbSeq:
3787        for i,rbseq in enumerate(rbSeq):
3788#                nameLst = [atNames[i] for i in rbseq[3]]
3789            out.append('  Torsion sequence %d Bond: %s - %s riding: %s\n'%
3790                    (i,atNames[rbseq[0]],atNames[rbseq[1]],str([atNames[i] for i in rbseq[3]])))
3791    return out
3792
3793def WriteVecRBModel(RBModel,sigDict={},irb=None):
3794    '''Write description of a vector rigid body. Code shifted from
3795    PrintVecRBModel to make usable from G2export_CIF
3796    '''
3797    out = []
3798    rbRef = RBModel['rbRef']
3799    atTypes = RBModel['rbTypes']
3800    for i in range(len(RBModel['VectMag'])):
3801        if sigDict and irb is not None:
3802            name = '::RBV;'+str(i)+':'+str(irb)
3803            s = G2mth.ValEsd(RBModel['VectMag'][i],sigDict.get(name,-.0001))
3804            out.append('Vector no.: %d Magnitude: %s\n'%(i,s))
3805        else:
3806            out.append('Vector no.: %d Magnitude: %8.4f Refine? %s\n'%
3807                           (i,RBModel['VectMag'][i],RBModel['VectRef'][i]))
3808        out.append('  No. Type     vx         vy         vz\n')
3809        for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
3810            out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(j,name,xyz[0],xyz[1],xyz[2]))
3811    out.append('  No. Type      x          y          z\n')
3812    for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
3813        out.append(%d   %2s %10.4f %10.4f %10.4f\n'%(i,name,xyz[0],xyz[1],xyz[2]))
3814    return out
Note: See TracBrowser for help on using the repository browser.