source: trunk/GSASIIstrIO.py @ 5114

Last change on this file since 5114 was 5114, checked in by vondreele, 10 months ago

various fixes for EDX data fitting by LeBail?
various fixes for SAD analysis - plotting updates & importer
removed a couple of unused lambdas from G2sasd

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