source: trunk/GSASIIstrIO.py @ 3747

Last change on this file since 3747 was 3747, checked in by vondreele, 3 years ago

provide fix for missing reflection sets in LeBail? sequential refinements

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