source: trunk/GSASIIstrIO.py @ 4849

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

apply symmetry constraints on rigid body origin; show position vars that are fixed on sym-gen tab for constraints; minor bug fixes for missing keys

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