source: trunk/GSASIIstrIO.py @ 4859

Last change on this file since 4859 was 4859, checked in by vondreele, 7 months ago

use new symmetry constraint scheme in MakeRBParms

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