source: trunk/GSASIIstrIO.py @ 3392

Last change on this file since 3392 was 3392, checked in by vondreele, 4 years ago

fix issues with multiphase refinements & multiphase plots, etc.
move plot OffSet? parameters to Page.plotStyle - required some reorganization of PlotPatterns?
correct colors for peak fit lines (make always blue)
make plot limits behave

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