source: trunk/GSASIIstrIO.py @ 5147

Last change on this file since 5147 was 5147, checked in by vondreele, 11 months ago

Add new routine to G2lat - AplusDij? & use it in a couple of places
add a new column to the ISODISTORT mode display - shows values obtained from ISODISTORT
remove calls to GetHstrainShift? for all powder data cases - wrong thing to do.

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