source: trunk/GSASIIstrIO.py @ 5142

Last change on this file since 5142 was 5142, checked in by toby, 11 months ago

revise CIF output to include only unique su values

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