source: trunk/GSASIIstrIO.py @ 3373

Last change on this file since 3373 was 3373, checked in by toby, 3 years ago

make sequential refinement setting work like a global flag; add indicator to g2frame legend; clean up constrain selection by showing only wildcards in seq; and not showing them in for non-eq.; speed up some calls of G2stIO.GetHistogramPhaseData? but providing only current histogram (in seq. ref.)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 149.8 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2018-05-06 19:13:25 +0000 (Sun, 06 May 2018) $
4# $Author: toby $
5# $Revision: 3373 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 3373 2018-05-06 19:13:25Z toby $
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: 3373 $")
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            if phase not in Phases:
442                pId = phaseNames.index(phase)
443                Phase['pId'] = pId
444                Phases[phase] = Phase
445            for hist in Phase['Histograms']:
446                if 'Use' not in Phase['Histograms'][hist]:      #patch
447                    Phase['Histograms'][hist]['Use'] = True         
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    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1565    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']:
1566        vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0]
1567        vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0]
1568        vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2]
1569    elif SGData['SGLaue'] in ['m3','m3m']:
1570        vcov[0:3,0:3] = vcov[0,0]
1571    elif SGData['SGLaue'] in ['4/m', '4/mmm']:
1572        vcov[0:2,0:2] = vcov[0,0]
1573        vcov[1,2] = vcov[2,1] = vcov[0,2]
1574    elif SGData['SGLaue'] in ['3R','3mR']:
1575        vcov[0:3,0:3] = vcov[0,0]
1576#        vcov[4,4] = vcov[5,5] = vcov[3,3]
1577        vcov[3:6,3:6] = vcov[3,3]
1578        vcov[0:3,3:6] = vcov[0,3]
1579        vcov[3:6,0:3] = vcov[3,0]
1580    delt = 1.e-9
1581    drVdA = np.zeros(6)
1582    for i in range(6):
1583        A[i] += delt
1584        drVdA[i] = G2lat.calc_rVsq(A)
1585        A[i] -= 2*delt
1586        drVdA[i] -= G2lat.calc_rVsq(A)
1587        A[i] += delt
1588    drVdA /= 2.*delt   
1589    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1590    Vol = 1/np.sqrt(rVsq)
1591    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1592   
1593    dcdA = np.zeros((6,6))
1594    for i in range(6):
1595        pdcdA =np.zeros(6)
1596        A[i] += delt
1597        pdcdA += G2lat.A2cell(A)
1598        A[i] -= 2*delt
1599        pdcdA -= G2lat.A2cell(A)
1600        A[i] += delt
1601        dcdA[i] = pdcdA/(2.*delt)
1602   
1603    sigMat = np.inner(dcdA,np.inner(dcdA,vcov))
1604    var = np.diag(sigMat)
1605    CS = np.where(var>0.,np.sqrt(var),0.)
1606    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']:
1607        CS[3:6] = 0.0
1608    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1609   
1610def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1611    '''Called after a refinement to transfer parameters from the parameter dict to
1612    the phase(s) information read from a GPX file. Also prints values to the .lst file
1613    '''
1614   
1615    def PrintAtomsAndSig(General,Atoms,atomsSig):
1616        pFile.write('\n Atoms:\n')
1617        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1618        if General['Type'] == 'macromolecular':
1619            line = ' res no residue chain '+line
1620        cx,ct,cs,cia = General['AtomPtrs']
1621        pFile.write(line+'\n')
1622        pFile.write(135*'-'+'\n')
1623        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1624            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1625        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1626        for atyp in General['AtomTypes']:       #zero composition
1627            General['NoAtoms'][atyp] = 0.
1628        for i,at in enumerate(Atoms):
1629            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1630            if General['Type'] == 'macromolecular':
1631                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1632                valstr = ' values:          '
1633                sigstr = ' sig   :          '
1634            else:
1635                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1636                valstr = ' values:'
1637                sigstr = ' sig   :'
1638            for ind in range(cx,cx+4):
1639                sigind = str(i)+':'+str(ind)
1640                valstr += fmt[ind]%(at[ind])                   
1641                if sigind in atomsSig:
1642                    sigstr += fmt[ind]%(atomsSig[sigind])
1643                else:
1644                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1645            if at[cia] == 'I':
1646                valstr += fmt[cia+1]%(at[cia+1])
1647                if '%d:%d'%(i,cia+1) in atomsSig:
1648                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1649                else:
1650                    sigstr += 8*' '
1651            else:
1652                valstr += 8*' '
1653                sigstr += 8*' '
1654                for ind in range(cia+2,cia+8):
1655                    sigind = str(i)+':'+str(ind)
1656                    valstr += fmt[ind]%(at[ind])
1657                    if sigind in atomsSig:                       
1658                        sigstr += fmt[ind]%(atomsSig[sigind])
1659                    else:
1660                        sigstr += 8*' '
1661            pFile.write(name+'\n')
1662            pFile.write(valstr+'\n')
1663            pFile.write(sigstr+'\n')
1664           
1665    def PrintMomentsAndSig(General,Atoms,atomsSig):
1666        cell = General['Cell'][1:7]
1667        G = G2lat.fillgmat(cell)
1668        ast = np.sqrt(np.diag(G))
1669        GS = G/np.outer(ast,ast)
1670        pFile.write('\n Magnetic Moments:\n')    #add magnitude & angle, etc.? TBD
1671        line = '   name       Mx        My        Mz       |Mag|'
1672        cx,ct,cs,cia = General['AtomPtrs']
1673        cmx = cx+4
1674        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1675        pFile.write(line+'\n')
1676        pFile.write(135*'-'+'\n')
1677        fmt = {0:'%7s',ct:'%7s',cmx:'%10.3f',cmx+1:'%10.3f',cmx+2:'%10.3f'}
1678        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
1679        for i,at in enumerate(Atoms):
1680            if AtInfo[at[ct]]:
1681                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1682                valstr = ' values:'
1683                sigstr = ' sig   :'
1684                for ind in range(cmx,cmx+3):
1685                    sigind = str(i)+':'+str(ind)
1686                    valstr += fmt[ind]%(at[ind])                   
1687                    if sigind in atomsSig:
1688                        sigstr += fmt[ind]%(atomsSig[sigind])
1689                    else:
1690                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1691                mag = np.array(at[cmx:cmx+3])
1692                Mag = np.sqrt(np.inner(mag,np.inner(mag,GS)))
1693                valstr += '%10.3f'%Mag
1694                sigstr += 10*' '
1695                pFile.write(name+'\n')
1696                pFile.write(valstr+'\n')
1697                pFile.write(sigstr+'\n')
1698           
1699    def PrintWavesAndSig(General,Atoms,wavesSig):
1700        cx,ct,cs,cia = General['AtomPtrs']
1701        pFile.write('\n Modulation waves\n')
1702        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1703            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1704            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1705        pFile.write(135*'-'+'\n')
1706        for i,at in enumerate(Atoms):
1707            AtomSS = at[-1]['SS1']
1708            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1709                Waves = AtomSS[Stype]
1710                if len(Waves):
1711                    waveType = Waves[0]
1712                else:
1713                    continue
1714                if len(Waves):
1715                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
1716                        (at[ct-1],at[cs],Stype,waveType))
1717                    for iw,wave in enumerate(Waves[1:]):
1718                        stiw = ':'+str(i)+':'+str(iw)
1719                        namstr = '  names :'
1720                        valstr = '  values:'
1721                        sigstr = '  esds  :'
1722                        if Stype == 'Spos':
1723                            nt = 6
1724                            ot = 0
1725                            if waveType in ['ZigZag','Block',] and not iw:
1726                                nt = 5
1727                                ot = 6
1728                            for j in range(nt):
1729                                name = names['Spos'][j+ot]
1730                                namstr += '%12s'%(name)
1731                                valstr += '%12.4f'%(wave[0][j])
1732                                if name+stiw in wavesSig:
1733                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1734                                else:
1735                                    sigstr += 12*' '
1736                        elif Stype == 'Sfrac':
1737                            ot = 0
1738                            if 'Crenel' in waveType and not iw:
1739                                ot = 2
1740                            for j in range(2):
1741                                name = names['Sfrac'][j+ot]
1742                                namstr += '%12s'%(names['Sfrac'][j+ot])
1743                                valstr += '%12.4f'%(wave[0][j])
1744                                if name+stiw in wavesSig:
1745                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1746                                else:
1747                                    sigstr += 12*' '
1748                        elif Stype == 'Sadp':
1749                            for j in range(12):
1750                                name = names['Sadp'][j]
1751                                namstr += '%10s'%(names['Sadp'][j])
1752                                valstr += '%10.6f'%(wave[0][j])
1753                                if name+stiw in wavesSig:
1754                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
1755                                else:
1756                                    sigstr += 10*' '
1757                        elif Stype == 'Smag':
1758                            for j in range(6):
1759                                name = names['Smag'][j]
1760                                namstr += '%12s'%(names['Smag'][j])
1761                                valstr += '%12.4f'%(wave[0][j])
1762                                if name+stiw in wavesSig:
1763                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1764                                else:
1765                                    sigstr += 12*' '
1766                               
1767                    pFile.write(namstr+'\n')
1768                    pFile.write(valstr+'\n')
1769                    pFile.write(sigstr+'\n')
1770       
1771               
1772    def PrintRBObjPOAndSig(rbfx,rbsx):
1773        namstr = '  names :'
1774        valstr = '  values:'
1775        sigstr = '  esds  :'
1776        for i,px in enumerate(['Px:','Py:','Pz:']):
1777            name = pfx+rbfx+px+rbsx
1778            namstr += '%12s'%('Pos '+px[1])
1779            valstr += '%12.5f'%(parmDict[name])
1780            if name in sigDict:
1781                sigstr += '%12.5f'%(sigDict[name])
1782            else:
1783                sigstr += 12*' '
1784        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1785            name = pfx+rbfx+po+rbsx
1786            namstr += '%12s'%('Ori '+po[1])
1787            valstr += '%12.5f'%(parmDict[name])
1788            if name in sigDict:
1789                sigstr += '%12.5f'%(sigDict[name])
1790            else:
1791                sigstr += 12*' '
1792        pFile.write(namstr+'\n')
1793        pFile.write(valstr+'\n')
1794        pFile.write(sigstr+'\n')
1795       
1796    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1797        namstr = '  names :'
1798        valstr = '  values:'
1799        sigstr = '  esds  :'
1800        if 'N' not in TLS:
1801            pFile.write(' Thermal motion:\n')
1802        if 'T' in TLS:
1803            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1804                name = pfx+rbfx+pt+rbsx
1805                namstr += '%12s'%(pt[:3])
1806                valstr += '%12.5f'%(parmDict[name])
1807                if name in sigDict:
1808                    sigstr += '%12.5f'%(sigDict[name])
1809                else:
1810                    sigstr += 12*' '
1811            pFile.write(namstr+'\n')
1812            pFile.write(valstr+'\n')
1813            pFile.write(sigstr+'\n')
1814        if 'L' in TLS:
1815            namstr = '  names :'
1816            valstr = '  values:'
1817            sigstr = '  esds  :'
1818            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1819                name = pfx+rbfx+pt+rbsx
1820                namstr += '%12s'%(pt[:3])
1821                valstr += '%12.3f'%(parmDict[name])
1822                if name in sigDict:
1823                    sigstr += '%12.3f'%(sigDict[name])
1824                else:
1825                    sigstr += 12*' '
1826            pFile.write(namstr+'\n')
1827            pFile.write(valstr+'\n')
1828            pFile.write(sigstr+'\n')
1829        if 'S' in TLS:
1830            namstr = '  names :'
1831            valstr = '  values:'
1832            sigstr = '  esds  :'
1833            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1834                name = pfx+rbfx+pt+rbsx
1835                namstr += '%12s'%(pt[:3])
1836                valstr += '%12.4f'%(parmDict[name])
1837                if name in sigDict:
1838                    sigstr += '%12.4f'%(sigDict[name])
1839                else:
1840                    sigstr += 12*' '
1841            pFile.write(namstr+'\n')
1842            pFile.write(valstr+'\n')
1843            pFile.write(sigstr+'\n')
1844        if 'U' in TLS:
1845            name = pfx+rbfx+'U:'+rbsx
1846            namstr = '  names :'+'%12s'%('Uiso')
1847            valstr = '  values:'+'%12.5f'%(parmDict[name])
1848            if name in sigDict:
1849                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1850            else:
1851                sigstr = '  esds  :'+12*' '
1852            pFile.write(namstr+'\n')
1853            pFile.write(valstr+'\n')
1854            pFile.write(sigstr+'\n')
1855       
1856    def PrintRBObjTorAndSig(rbsx):
1857        namstr = '  names :'
1858        valstr = '  values:'
1859        sigstr = '  esds  :'
1860        nTors = len(RBObj['Torsions'])
1861        if nTors:
1862            pFile.write(' Torsions:\n')
1863            for it in range(nTors):
1864                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1865                namstr += '%12s'%('Tor'+str(it))
1866                valstr += '%12.4f'%(parmDict[name])
1867                if name in sigDict:
1868                    sigstr += '%12.4f'%(sigDict[name])
1869            pFile.write(namstr+'\n')
1870            pFile.write(valstr+'\n')
1871            pFile.write(sigstr+'\n')
1872               
1873    def PrintSHtextureAndSig(textureData,SHtextureSig):
1874        pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
1875        names = ['omega','chi','phi']
1876        namstr = '  names :'
1877        ptstr =  '  values:'
1878        sigstr = '  esds  :'
1879        for name in names:
1880            namstr += '%12s'%(name)
1881            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1882            if 'Sample '+name in SHtextureSig:
1883                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1884            else:
1885                sigstr += 12*' '
1886        pFile.write(namstr+'\n')
1887        pFile.write(ptstr+'\n')
1888        pFile.write(sigstr+'\n')
1889        pFile.write('\n Texture coefficients:\n')
1890        SHcoeff = textureData['SH Coeff'][1]
1891        SHkeys = list(SHcoeff.keys())
1892        nCoeff = len(SHcoeff)
1893        nBlock = nCoeff//10+1
1894        iBeg = 0
1895        iFin = min(iBeg+10,nCoeff)
1896        for block in range(nBlock):
1897            namstr = '  names :'
1898            ptstr =  '  values:'
1899            sigstr = '  esds  :'
1900            for name in SHkeys[iBeg:iFin]:
1901                namstr += '%12s'%(name)
1902                ptstr += '%12.3f'%(SHcoeff[name])
1903                if name in SHtextureSig:
1904                    sigstr += '%12.3f'%(SHtextureSig[name])
1905                else:
1906                    sigstr += 12*' '
1907            pFile.write(namstr+'\n')
1908            pFile.write(ptstr+'\n')
1909            pFile.write(sigstr+'\n')
1910            iBeg += 10
1911            iFin = min(iBeg+10,nCoeff)
1912           
1913    ##########################################################################
1914    # SetPhaseData starts here
1915    pFile.write('\n Phases:\n')
1916    for phase in Phases:
1917        pFile.write(' Result for phase: %s\n'%phase)
1918        pFile.write(135*'='+'\n')
1919        Phase = Phases[phase]
1920        General = Phase['General']
1921        SGData = General['SGData']
1922        Atoms = Phase['Atoms']
1923        AtLookup = []
1924        if Atoms and not General.get('doPawley'):
1925            cx,ct,cs,cia = General['AtomPtrs']
1926            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1927        cell = General['Cell']
1928        pId = Phase['pId']
1929        pfx = str(pId)+'::'
1930        if cell[0]:
1931            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1932            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1933            pFile.write(' Reciprocal metric tensor: \n')
1934            ptfmt = "%15.9f"
1935            names = ['A11','A22','A33','A12','A13','A23']
1936            namstr = '  names :'
1937            ptstr =  '  values:'
1938            sigstr = '  esds  :'
1939            for name,a,siga in zip(names,A,sigA):
1940                namstr += '%15s'%(name)
1941                ptstr += ptfmt%(a)
1942                if siga:
1943                    sigstr += ptfmt%(siga)
1944                else:
1945                    sigstr += 15*' '
1946            pFile.write(namstr+'\n')
1947            pFile.write(ptstr+'\n')
1948            pFile.write(sigstr+'\n')
1949            cell[1:7] = G2lat.A2cell(A)
1950            cell[7] = G2lat.calc_V(A)
1951            pFile.write(' New unit cell:\n')
1952            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1953            names = ['a','b','c','alpha','beta','gamma','Volume']
1954            namstr = '  names :'
1955            ptstr =  '  values:'
1956            sigstr = '  esds  :'
1957            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1958                namstr += '%12s'%(name)
1959                ptstr += fmt%(a)
1960                if siga:
1961                    sigstr += fmt%(siga)
1962                else:
1963                    sigstr += 12*' '
1964            pFile.write(namstr+'\n')
1965            pFile.write(ptstr+'\n')
1966            pFile.write(sigstr+'\n')
1967        ik = 6  #for Pawley stuff below
1968        if General.get('Modulated',False):
1969            ik = 7
1970            Vec,vRef,maxH = General['SuperVec']
1971            if vRef:
1972                pFile.write(' New modulation vector:\n')
1973                namstr = '  names :'
1974                ptstr =  '  values:'
1975                sigstr = '  esds  :'
1976                for var in ['mV0','mV1','mV2']:
1977                    namstr += '%12s'%(pfx+var)
1978                    ptstr += '%12.6f'%(parmDict[pfx+var])
1979                    if pfx+var in sigDict:
1980                        sigstr += '%12.6f'%(sigDict[pfx+var])
1981                    else:
1982                        sigstr += 12*' '
1983                pFile.write(namstr+'\n')
1984                pFile.write(ptstr+'\n')
1985                pFile.write(sigstr+'\n')
1986           
1987        General['Mass'] = 0.
1988        if Phase['General'].get('doPawley'):
1989            pawleyRef = Phase['Pawley ref']
1990            for i,refl in enumerate(pawleyRef):
1991                key = pfx+'PWLref:'+str(i)
1992                refl[ik] = parmDict[key]
1993                if key in sigDict:
1994                    refl[ik+1] = sigDict[key]
1995                else:
1996                    refl[ik+1] = 0
1997        else:
1998            VRBIds = RBIds['Vector']
1999            RRBIds = RBIds['Residue']
2000            RBModels = Phase['RBModels']
2001            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
2002                jrb = VRBIds.index(RBObj['RBId'])
2003                rbsx = str(irb)+':'+str(jrb)
2004                pFile.write(' Vector rigid body parameters:\n')
2005                PrintRBObjPOAndSig('RBV',rbsx)
2006                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
2007            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
2008                jrb = RRBIds.index(RBObj['RBId'])
2009                rbsx = str(irb)+':'+str(jrb)
2010                pFile.write(' Residue rigid body parameters:\n')
2011                PrintRBObjPOAndSig('RBR',rbsx)
2012                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
2013                PrintRBObjTorAndSig(rbsx)
2014            atomsSig = {}
2015            wavesSig = {}
2016            cx,ct,cs,cia = General['AtomPtrs']
2017            for i,at in enumerate(Atoms):
2018                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
2019                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
2020                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
2021                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
2022                for ind in range(cx,cx+4):
2023                    at[ind] = parmDict[names[ind]]
2024                    if ind in range(cx,cx+3):
2025                        name = names[ind].replace('A','dA')
2026                    else:
2027                        name = names[ind]
2028                    if name in sigDict:
2029                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2030                if at[cia] == 'I':
2031                    at[cia+1] = parmDict[names[cia+1]]
2032                    if names[cia+1] in sigDict:
2033                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2034                else:
2035                    for ind in range(cia+2,cia+8):
2036                        at[ind] = parmDict[names[ind]]
2037                        if names[ind] in sigDict:
2038                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2039                if General['Type'] == 'magnetic':
2040                    for ind in range(cx+4,cx+7):
2041                        at[ind] = parmDict[names[ind]]
2042                        if names[ind] in sigDict:
2043                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2044                ind = General['AtomTypes'].index(at[ct])
2045                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2046                if General.get('Modulated',False):
2047                    AtomSS = at[-1]['SS1']
2048                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2049                        Waves = AtomSS[Stype]
2050                        if len(Waves):
2051                            waveType = Waves[0]
2052                        else:
2053                            continue
2054                        for iw,wave in enumerate(Waves[1:]):
2055                            stiw = str(i)+':'+str(iw)
2056                            if Stype == 'Spos':
2057                                if waveType in ['ZigZag','Block',] and not iw:
2058                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2059                                else:
2060                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2061                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2062                            elif Stype == 'Sadp':
2063                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2064                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2065                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2066                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2067                            elif Stype == 'Sfrac':
2068                                if 'Crenel' in waveType and not iw:
2069                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2070                                else:
2071                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2072                            elif Stype == 'Smag':
2073                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2074                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2075                            for iname,name in enumerate(names):
2076                                AtomSS[Stype][iw+1][0][iname] = parmDict[pfx+name]
2077                                if pfx+name in sigDict:
2078                                    wavesSig[name] = sigDict[pfx+name]
2079                   
2080            PrintAtomsAndSig(General,Atoms,atomsSig)
2081            if General['Type'] == 'magnetic':
2082                PrintMomentsAndSig(General,Atoms,atomsSig)
2083            if General.get('Modulated',False):
2084                PrintWavesAndSig(General,Atoms,wavesSig)
2085           
2086       
2087        textureData = General['SH Texture']   
2088        if textureData['Order']:
2089            SHtextureSig = {}
2090            for name in ['omega','chi','phi']:
2091                aname = pfx+'SH '+name
2092                textureData['Sample '+name][1] = parmDict[aname]
2093                if aname in sigDict:
2094                    SHtextureSig['Sample '+name] = sigDict[aname]
2095            for name in textureData['SH Coeff'][1]:
2096                aname = pfx+name
2097                textureData['SH Coeff'][1][name] = parmDict[aname]
2098                if aname in sigDict:
2099                    SHtextureSig[name] = sigDict[aname]
2100            PrintSHtextureAndSig(textureData,SHtextureSig)
2101        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2102            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2103                textureData,RestraintDict[phase],pFile)
2104                   
2105################################################################################
2106##### Histogram & Phase data
2107################################################################################       
2108                   
2109def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2110    '''Loads the HAP histogram/phase information into dicts
2111
2112    :param dict Phases: phase information
2113    :param dict Histograms: Histogram information
2114    :param bool Print: prints information as it is read
2115    :param file pFile: file object to print to (the default, None causes printing to the console)
2116    :param bool resetRefList: Should the contents of the Reflection List be initialized
2117      on loading. The default, True, initializes the Reflection List as it is loaded.
2118
2119    :returns: (hapVary,hapDict,controlDict)
2120      * hapVary: list of refined variables
2121      * hapDict: dict with refined variables and their values
2122      * controlDict: dict with fixed parameters
2123    '''
2124   
2125    def PrintSize(hapData):
2126        if hapData[0] in ['isotropic','uniaxial']:
2127            line = '\n Size model    : %9s'%(hapData[0])
2128            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2129            if hapData[0] == 'uniaxial':
2130                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2131            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2132            pFile.write(line+'\n')
2133        else:
2134            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2135                (hapData[0],hapData[1][2],hapData[2][2]))
2136            Snames = ['S11','S22','S33','S12','S13','S23']
2137            ptlbls = ' names :'
2138            ptstr =  ' values:'
2139            varstr = ' refine:'
2140            for i,name in enumerate(Snames):
2141                ptlbls += '%12s' % (name)
2142                ptstr += '%12.3f' % (hapData[4][i])
2143                varstr += '%12s' % (str(hapData[5][i]))
2144            pFile.write(ptlbls+'\n')
2145            pFile.write(ptstr+'\n')
2146            pFile.write(varstr+'\n')
2147       
2148    def PrintMuStrain(hapData,SGData):
2149        if hapData[0] in ['isotropic','uniaxial']:
2150            line = '\n Mustrain model: %9s'%(hapData[0])
2151            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2152            if hapData[0] == 'uniaxial':
2153                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2154            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2155            pFile.write(line+'\n')
2156        else:
2157            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2158                (hapData[0],hapData[1][2],hapData[2][2]))
2159            Snames = G2spc.MustrainNames(SGData)
2160            ptlbls = ' names :'
2161            ptstr =  ' values:'
2162            varstr = ' refine:'
2163            for i,name in enumerate(Snames):
2164                ptlbls += '%12s' % (name)
2165                ptstr += '%12.1f' % (hapData[4][i])
2166                varstr += '%12s' % (str(hapData[5][i]))
2167            pFile.write(ptlbls+'\n')
2168            pFile.write(ptstr+'\n')
2169            pFile.write(varstr+'\n')
2170
2171    def PrintHStrain(hapData,SGData):
2172        pFile.write('\n Hydrostatic/elastic strain:\n')
2173        Hsnames = G2spc.HStrainNames(SGData)
2174        ptlbls = ' names :'
2175        ptstr =  ' values:'
2176        varstr = ' refine:'
2177        for i,name in enumerate(Hsnames):
2178            ptlbls += '%12s' % (name)
2179            ptstr += '%12.4g' % (hapData[0][i])
2180            varstr += '%12s' % (str(hapData[1][i]))
2181        pFile.write(ptlbls+'\n')
2182        pFile.write(ptstr+'\n')
2183        pFile.write(varstr+'\n')
2184
2185    def PrintSHPO(hapData):
2186        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2187        ptlbls = ' names :'
2188        ptstr =  ' values:'
2189        for item in hapData[5]:
2190            ptlbls += '%12s'%(item)
2191            ptstr += '%12.3f'%(hapData[5][item]) 
2192        pFile.write(ptlbls+'\n')
2193        pFile.write(ptstr+'\n')
2194   
2195    def PrintBabinet(hapData):
2196        pFile.write('\n Babinet form factor modification:\n')
2197        ptlbls = ' names :'
2198        ptstr =  ' values:'
2199        varstr = ' refine:'
2200        for name in ['BabA','BabU']:
2201            ptlbls += '%12s' % (name)
2202            ptstr += '%12.6f' % (hapData[name][0])
2203            varstr += '%12s' % (str(hapData[name][1]))
2204        pFile.write(ptlbls+'\n')
2205        pFile.write(ptstr+'\n')
2206        pFile.write(varstr+'\n')
2207       
2208    hapDict = {}
2209    hapVary = []
2210    controlDict = {}
2211   
2212    for phase in Phases:
2213        HistoPhase = Phases[phase]['Histograms']
2214        SGData = Phases[phase]['General']['SGData']
2215        cell = Phases[phase]['General']['Cell'][1:7]
2216        A = G2lat.cell2A(cell)
2217        if Phases[phase]['General'].get('Modulated',False):
2218            SSGData = Phases[phase]['General']['SSGData']
2219            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2220        pId = Phases[phase]['pId']
2221        histoList = list(HistoPhase.keys())
2222        histoList.sort()
2223        for histogram in histoList:
2224            try:
2225                Histogram = Histograms[histogram]
2226            except KeyError:                       
2227                #skip if histogram not included e.g. in a sequential refinement
2228                continue
2229            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2230                if phase in Histograms[histogram]['Reflection Lists']:
2231                    del Histograms[histogram]['Reflection Lists'][phase]
2232                continue
2233            hapData = HistoPhase[histogram]
2234            hId = Histogram['hId']
2235            if 'PWDR' in histogram:
2236                limits = Histogram['Limits'][1]
2237                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2238                if 'C' in inst['Type'][1]:
2239                    try:
2240                        wave = inst['Lam'][1]
2241                    except KeyError:
2242                        wave = inst['Lam1'][1]
2243                    dmin = wave/(2.0*sind(limits[1]/2.0))
2244                elif 'T' in inst['Type'][0]:
2245                    dmin = limits[0]/inst['difC'][1]
2246                pfx = str(pId)+':'+str(hId)+':'
2247                if Phases[phase]['General']['doPawley']:
2248                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2249                    hapDict[pfx+'newLeBail'] = True
2250                    Tmin = G2lat.Dsp2pos(inst,dmin)
2251                    if 'C' in inst['Type'][1]:
2252                        limits[1] = min(limits[1],Tmin)
2253                    else:
2254                        limits[0] = max(limits[0],Tmin)
2255                else:
2256                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2257                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2258                if Phases[phase]['General']['Type'] == 'magnetic':
2259                    dmin = max(dmin,Phases[phase]['General']['MagDmin'])
2260                for item in ['Scale','Extinction']:
2261                    hapDict[pfx+item] = hapData[item][0]
2262                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2263                        hapVary.append(pfx+item)
2264                names = G2spc.HStrainNames(SGData)
2265                HSvals = []
2266                for i,name in enumerate(names):
2267                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2268                    HSvals.append(hapDict[pfx+name])
2269                    if hapData['HStrain'][1][i]:
2270#                    if hapData['HStrain'][1][i] and not hapDict[pfx+'LeBail']:
2271                        hapVary.append(pfx+name)
2272                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2273                if hapData['Pref.Ori.'][0] == 'MD':
2274                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2275                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2276                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2277                        hapVary.append(pfx+'MD')
2278                else:                           #'SH' spherical harmonics
2279                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2280                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2281                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2282                    controlDict[pfx+'SHhkl'] = []
2283                    try: #patch for old Pref.Ori. items
2284                        controlDict[pfx+'SHtoler'] = 0.1
2285                        if hapData['Pref.Ori.'][6][0] != '':
2286                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2287                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2288                    except IndexError:
2289                        pass
2290                    for item in hapData['Pref.Ori.'][5]:
2291                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2292                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2293                            hapVary.append(pfx+item)
2294                for item in ['Mustrain','Size']:
2295                    controlDict[pfx+item+'Type'] = hapData[item][0]
2296                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2297                    if hapData[item][2][2]:
2298                        hapVary.append(pfx+item+';mx')
2299                    if hapData[item][0] in ['isotropic','uniaxial']:
2300                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2301                        if hapData[item][2][0]:
2302                            hapVary.append(pfx+item+';i')
2303                        if hapData[item][0] == 'uniaxial':
2304                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2305                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2306                            if hapData[item][2][1]:
2307                                hapVary.append(pfx+item+';a')
2308                    else:       #generalized for mustrain or ellipsoidal for size
2309                        Nterms = len(hapData[item][4])
2310                        if item == 'Mustrain':
2311                            names = G2spc.MustrainNames(SGData)
2312                            pwrs = []
2313                            for name in names:
2314                                h,k,l = name[1:]
2315                                pwrs.append([int(h),int(k),int(l)])
2316                            controlDict[pfx+'MuPwrs'] = pwrs
2317                        for i in range(Nterms):
2318                            sfx = ';'+str(i)
2319                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2320                            if hapData[item][5][i]:
2321                                hapVary.append(pfx+item+sfx)
2322                if Phases[phase]['General']['Type'] != 'magnetic':
2323                    for bab in ['BabA','BabU']:
2324                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2325                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2326                            hapVary.append(pfx+bab)
2327                               
2328                if Print: 
2329                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2330                    pFile.write(135*'='+'\n')
2331                    if hapDict[pfx+'LeBail']:
2332                        pFile.write(' Perform LeBail extraction\n')                     
2333                    else:
2334                        pFile.write(' Phase fraction  : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2335                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2336                        if hapData['Pref.Ori.'][0] == 'MD':
2337                            Ax = hapData['Pref.Ori.'][3]
2338                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2339                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2340                        else: #'SH' for spherical harmonics
2341                            PrintSHPO(hapData['Pref.Ori.'])
2342                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2343                    PrintSize(hapData['Size'])
2344                    PrintMuStrain(hapData['Mustrain'],SGData)
2345                    PrintHStrain(hapData['HStrain'],SGData)
2346                    if Phases[phase]['General']['Type'] != 'magnetic':
2347                        if hapData['Babinet']['BabA'][0]:
2348                            PrintBabinet(hapData['Babinet'])                       
2349                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData['LeBail'] and hapData['newLeBail'])):
2350                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2351                        hapData['newLeBail'] = False
2352                    refList = []
2353                    Uniq = []
2354                    Phi = []
2355                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2356                    if Phases[phase]['General'].get('Modulated',False):
2357                        ifSuper = True
2358                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2359                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2360                        for h,k,l,m,d in HKLd:
2361                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2362                            mul *= 2      # for powder overlap of Friedel pairs
2363                            if m or not ext or useExt:
2364                                if 'C' in inst['Type'][0]:
2365                                    pos = G2lat.Dsp2pos(inst,d)
2366                                    if limits[0] < pos < limits[1]:
2367                                        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])
2368                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2369                                        Uniq.append(uniq)
2370                                        Phi.append(phi)
2371                                elif 'T' in inst['Type'][0]:
2372                                    pos = G2lat.Dsp2pos(inst,d)
2373                                    if limits[0] < pos < limits[1]:
2374                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2375                                        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])
2376                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2377                                        #TODO - if tabulated put alp & bet in here
2378                                        Uniq.append(uniq)
2379                                        Phi.append(phi)
2380                    else:
2381                        ifSuper = False
2382                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2383                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2384                        for h,k,l,d in HKLd:
2385                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2386                            mul *= 2      # for powder overlap of Friedel pairs
2387                            if ext and not useExt:
2388                                continue
2389                            if 'C' in inst['Type'][0]:
2390                                pos = G2lat.Dsp2pos(inst,d)
2391                                if limits[0] < pos < limits[1]:
2392                                    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])
2393                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2394                                    Uniq.append(uniq)
2395                                    Phi.append(phi)
2396                            elif 'T' in inst['Type'][0]:
2397                                pos = G2lat.Dsp2pos(inst,d)
2398                                if limits[0] < pos < limits[1]:
2399                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2400                                    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])
2401                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2402                                    Uniq.append(uniq)
2403                                    Phi.append(phi)
2404                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2405            elif 'HKLF' in histogram:
2406                inst = Histogram['Instrument Parameters'][0]
2407                hId = Histogram['hId']
2408                hfx = ':%d:'%(hId)
2409                for item in inst:
2410                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2411                    hapDict[hfx+item] = inst[item][1]
2412                pfx = str(pId)+':'+str(hId)+':'
2413                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2414                if hapData['Scale'][1]:
2415                    hapVary.append(pfx+'Scale')
2416                               
2417                extApprox,extType,extParms = hapData['Extinction']
2418                controlDict[pfx+'EType'] = extType
2419                controlDict[pfx+'EApprox'] = extApprox
2420                if 'C' in inst['Type'][0]:
2421                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2422                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2423                if 'Primary' in extType:
2424                    Ekey = ['Ep',]
2425                elif 'I & II' in extType:
2426                    Ekey = ['Eg','Es']
2427                elif 'Secondary Type II' == extType:
2428                    Ekey = ['Es',]
2429                elif 'Secondary Type I' == extType:
2430                    Ekey = ['Eg',]
2431                else:   #'None'
2432                    Ekey = []
2433                for eKey in Ekey:
2434                    hapDict[pfx+eKey] = extParms[eKey][0]
2435                    if extParms[eKey][1]:
2436                        hapVary.append(pfx+eKey)
2437                for bab in ['BabA','BabU']:
2438                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2439                    if hapData['Babinet'][bab][1]:
2440                        hapVary.append(pfx+bab)
2441                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2442                if len(Twins) == 1:
2443                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2444                    if hapData.get('Flack',[0,False])[1]:
2445                        hapVary.append(pfx+'Flack')
2446                sumTwFr = 0.
2447                controlDict[pfx+'TwinLaw'] = []
2448                controlDict[pfx+'TwinInv'] = []
2449                NTL = 0           
2450                for it,twin in enumerate(Twins):
2451                    if 'bool' in str(type(twin[0])):
2452                        controlDict[pfx+'TwinInv'].append(twin[0])
2453                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2454                    else:
2455                        NTL += 1
2456                        controlDict[pfx+'TwinInv'].append(False)
2457                        controlDict[pfx+'TwinLaw'].append(twin[0])
2458                    if it:
2459                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2460                        sumTwFr += twin[1]
2461                    else:
2462                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2463                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2464                    if Twins[0][1][1]:
2465                        hapVary.append(pfx+'TwinFr:'+str(it))
2466                controlDict[pfx+'NTL'] = NTL
2467                #need to make constraint on TwinFr
2468                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2469                if len(Twins) > 1:    #force sum to unity
2470                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2471                if Print: 
2472                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2473                    pFile.write(135*'='+'\n')
2474                    pFile.write(' Scale factor     : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2475                    if extType != 'None':
2476                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2477                        text = ' Parameters       :'
2478                        for eKey in Ekey:
2479                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2480                        pFile.write(text+'\n')
2481                    if hapData['Babinet']['BabA'][0]:
2482                        PrintBabinet(hapData['Babinet'])
2483                    if not SGData['SGInv'] and len(Twins) == 1:
2484                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2485                    if len(Twins) > 1:
2486                        for it,twin in enumerate(Twins):
2487                            if 'bool' in str(type(twin[0])):
2488                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine?\n'%
2489                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1]) 
2490                            else:
2491                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2492                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],Twins[0][1][1]))
2493                       
2494                Histogram['Reflection Lists'] = phase       
2495               
2496    return hapVary,hapDict,controlDict
2497   
2498def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2499    'needs a doc string'
2500   
2501    def PrintSizeAndSig(hapData,sizeSig):
2502        line = '\n Size model:     %9s'%(hapData[0])
2503        refine = False
2504        if hapData[0] in ['isotropic','uniaxial']:
2505            line += ' equatorial:%12.5f'%(hapData[1][0])
2506            if sizeSig[0][0]:
2507                line += ', sig:%8.4f'%(sizeSig[0][0])
2508                refine = True
2509            if hapData[0] == 'uniaxial':
2510                line += ' axial:%12.4f'%(hapData[1][1])
2511                if sizeSig[0][1]:
2512                    refine = True
2513                    line += ', sig:%8.4f'%(sizeSig[0][1])
2514            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2515            if sizeSig[0][2]:
2516                refine = True
2517                line += ', sig:%8.4f'%(sizeSig[0][2])
2518            if refine:
2519                pFile.write(line+'\n')
2520        else:
2521            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2522            if sizeSig[0][2]:
2523                refine = True
2524                line += ', sig:%8.4f'%(sizeSig[0][2])
2525            Snames = ['S11','S22','S33','S12','S13','S23']
2526            ptlbls = ' name  :'
2527            ptstr =  ' value :'
2528            sigstr = ' sig   :'
2529            for i,name in enumerate(Snames):
2530                ptlbls += '%12s' % (name)
2531                ptstr += '%12.6f' % (hapData[4][i])
2532                if sizeSig[1][i]:
2533                    refine = True
2534                    sigstr += '%12.6f' % (sizeSig[1][i])
2535                else:
2536                    sigstr += 12*' '
2537            if refine:
2538                pFile.write(line+'\n')
2539                pFile.write(ptlbls+'\n')
2540                pFile.write(ptstr+'\n')
2541                pFile.write(sigstr+'\n')
2542       
2543    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2544        line = '\n Mustrain model: %9s\n'%(hapData[0])
2545        refine = False
2546        if hapData[0] in ['isotropic','uniaxial']:
2547            line += ' equatorial:%12.1f'%(hapData[1][0])
2548            if mustrainSig[0][0]:
2549                line += ', sig:%8.1f'%(mustrainSig[0][0])
2550                refine = True
2551            if hapData[0] == 'uniaxial':
2552                line += ' axial:%12.1f'%(hapData[1][1])
2553                if mustrainSig[0][1]:
2554                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2555            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2556            if mustrainSig[0][2]:
2557                refine = True
2558                line += ', sig:%8.3f'%(mustrainSig[0][2])
2559            if refine:
2560                pFile.write(line+'\n')
2561        else:
2562            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2563            if mustrainSig[0][2]:
2564                refine = True
2565                line += ', sig:%8.3f'%(mustrainSig[0][2])
2566            Snames = G2spc.MustrainNames(SGData)
2567            ptlbls = ' name  :'
2568            ptstr =  ' value :'
2569            sigstr = ' sig   :'
2570            for i,name in enumerate(Snames):
2571                ptlbls += '%12s' % (name)
2572                ptstr += '%12.1f' % (hapData[4][i])
2573                if mustrainSig[1][i]:
2574                    refine = True
2575                    sigstr += '%12.1f' % (mustrainSig[1][i])
2576                else:
2577                    sigstr += 12*' '
2578            if refine:
2579                pFile.write(line+'\n')
2580                pFile.write(ptlbls+'\n')
2581                pFile.write(ptstr+'\n')
2582                pFile.write(sigstr+'\n')
2583           
2584    def PrintHStrainAndSig(hapData,strainSig,SGData):
2585        Hsnames = G2spc.HStrainNames(SGData)
2586        ptlbls = ' name  :'
2587        ptstr =  ' value :'
2588        sigstr = ' sig   :'
2589        refine = False
2590        for i,name in enumerate(Hsnames):
2591            ptlbls += '%12s' % (name)
2592            ptstr += '%12.4g' % (hapData[0][i])
2593            if name in strainSig:
2594                refine = True
2595                sigstr += '%12.4g' % (strainSig[name])
2596            else:
2597                sigstr += 12*' '
2598        if refine:
2599            pFile.write('\n Hydrostatic/elastic strain:\n')
2600            pFile.write(ptlbls+'\n')
2601            pFile.write(ptstr+'\n')
2602            pFile.write(sigstr+'\n')
2603       
2604    def PrintSHPOAndSig(pfx,hapData,POsig):
2605        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2606        ptlbls = ' names :'
2607        ptstr =  ' values:'
2608        sigstr = ' sig   :'
2609        for item in hapData[5]:
2610            ptlbls += '%12s'%(item)
2611            ptstr += '%12.3f'%(hapData[5][item])
2612            if pfx+item in POsig:
2613                sigstr += '%12.3f'%(POsig[pfx+item])
2614            else:
2615                sigstr += 12*' ' 
2616        pFile.write(ptlbls+'\n')
2617        pFile.write(ptstr+'\n')
2618        pFile.write(sigstr+'\n')
2619       
2620    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2621        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2622        text = ''
2623        for item in hapData[2]:
2624            if pfx+item in ScalExtSig:
2625                text += '       %s: '%(item)
2626                text += '%12.2e'%(hapData[2][item][0])
2627                if pfx+item in ScalExtSig:
2628                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2629        pFile.write(text+'\n')   
2630       
2631    def PrintBabinetAndSig(pfx,hapData,BabSig):
2632        pFile.write('\n Babinet form factor modification:\n')
2633        ptlbls = ' names :'
2634        ptstr =  ' values:'
2635        sigstr = ' sig   :'
2636        for item in hapData:
2637            ptlbls += '%12s'%(item)
2638            ptstr += '%12.3f'%(hapData[item][0])
2639            if pfx+item in BabSig:
2640                sigstr += '%12.3f'%(BabSig[pfx+item])
2641            else:
2642                sigstr += 12*' ' 
2643        pFile.write(ptlbls+'\n')
2644        pFile.write(ptstr+'\n')
2645        pFile.write(sigstr+'\n')
2646       
2647    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2648        pFile.write('\n Twin Law fractions :\n')
2649        ptlbls = ' names :'
2650        ptstr =  ' values:'
2651        sigstr = ' sig   :'
2652        for it,item in enumerate(twinData):
2653            ptlbls += '     twin #%d'%(it)
2654            if it:
2655                ptstr += '%12.3f'%(item[1])
2656            else:
2657                ptstr += '%12.3f'%(item[1][0])
2658            if pfx+'TwinFr:'+str(it) in TwinSig:
2659                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2660            else:
2661                sigstr += 12*' ' 
2662        pFile.write(ptlbls+'\n')
2663        pFile.write(ptstr+'\n')
2664        pFile.write(sigstr+'\n')
2665       
2666   
2667    PhFrExtPOSig = {}
2668    SizeMuStrSig = {}
2669    ScalExtSig = {}
2670    BabSig = {}
2671    TwinFrSig = {}
2672    wtFrSum = {}
2673    for phase in Phases:
2674        HistoPhase = Phases[phase]['Histograms']
2675        General = Phases[phase]['General']
2676        SGData = General['SGData']
2677        pId = Phases[phase]['pId']
2678        histoList = list(HistoPhase.keys())
2679        histoList.sort()
2680        for histogram in histoList:
2681            try:
2682                Histogram = Histograms[histogram]
2683            except KeyError:                       
2684                #skip if histogram not included e.g. in a sequential refinement
2685                continue
2686            if not Phases[phase]['Histograms'][histogram]['Use']:
2687                #skip if phase absent from this histogram
2688                continue
2689            hapData = HistoPhase[histogram]
2690            hId = Histogram['hId']
2691            pfx = str(pId)+':'+str(hId)+':'
2692            if hId not in wtFrSum:
2693                wtFrSum[hId] = 0.
2694            if 'PWDR' in histogram:
2695                for item in ['Scale','Extinction']:
2696                    hapData[item][0] = parmDict[pfx+item]
2697                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2698                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2699                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2700                if hapData['Pref.Ori.'][0] == 'MD':
2701                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2702                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2703                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2704                else:                           #'SH' spherical harmonics
2705                    for item in hapData['Pref.Ori.'][5]:
2706                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2707                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2708                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2709                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2710                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2711                    pfx+'HStrain':{}})                 
2712                for item in ['Mustrain','Size']:
2713                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2714#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2715                    if pfx+item+';mx' in sigDict:
2716                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2717                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2718                        hapData[item][1][0] = parmDict[pfx+item+';i']
2719                        if item == 'Size':
2720                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2721                        if pfx+item+';i' in sigDict: 
2722                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2723                        if hapData[item][0] == 'uniaxial':
2724                            hapData[item][1][1] = parmDict[pfx+item+';a']
2725                            if item == 'Size':
2726                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2727                            if pfx+item+';a' in sigDict:
2728                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2729                    else:       #generalized for mustrain or ellipsoidal for size
2730                        Nterms = len(hapData[item][4])
2731                        for i in range(Nterms):
2732                            sfx = ';'+str(i)
2733                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2734                            if pfx+item+sfx in sigDict:
2735                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2736                names = G2spc.HStrainNames(SGData)
2737                for i,name in enumerate(names):
2738                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2739                    if pfx+name in sigDict:
2740                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2741                if Phases[phase]['General']['Type'] != 'magnetic':
2742                    for name in ['BabA','BabU']:
2743                        hapData['Babinet'][name][0] = parmDict[pfx+name]
2744                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
2745                            BabSig[pfx+name] = sigDict[pfx+name]               
2746               
2747            elif 'HKLF' in histogram:
2748                for item in ['Scale','Flack']:
2749                    if parmDict.get(pfx+item):
2750                        hapData[item][0] = parmDict[pfx+item]
2751                        if pfx+item in sigDict:
2752                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2753                for item in ['Ep','Eg','Es']:
2754                    if parmDict.get(pfx+item):
2755                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2756                        if pfx+item in sigDict:
2757                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2758                for item in ['BabA','BabU']:
2759                    hapData['Babinet'][item][0] = parmDict[pfx+item]
2760                    if pfx+item in sigDict:
2761                        BabSig[pfx+item] = sigDict[pfx+item]
2762                if 'Twins' in hapData:
2763                    it = 1
2764                    sumTwFr = 0.
2765                    while True:
2766                        try:
2767                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2768                            if pfx+'TwinFr:'+str(it) in sigDict:
2769                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2770                            if it:
2771                                sumTwFr += hapData['Twins'][it][1]
2772                            it += 1
2773                        except KeyError:
2774                            break
2775                    hapData['Twins'][0][1][0] = 1.-sumTwFr
2776
2777    if Print:
2778        for phase in Phases:
2779            HistoPhase = Phases[phase]['Histograms']
2780            General = Phases[phase]['General']
2781            SGData = General['SGData']
2782            pId = Phases[phase]['pId']
2783            histoList = list(HistoPhase.keys())
2784            histoList.sort()
2785            for histogram in histoList:
2786                try:
2787                    Histogram = Histograms[histogram]
2788                except KeyError:                       
2789                    #skip if histogram not included e.g. in a sequential refinement
2790                    continue
2791                hapData = HistoPhase[histogram]
2792                hId = Histogram['hId']
2793                Histogram['Residuals'][str(pId)+'::Name'] = phase
2794                pfx = str(pId)+':'+str(hId)+':'
2795                hfx = ':%s:'%(hId)
2796                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
2797                    continue
2798                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2799                pFile.write(135*'='+'\n')
2800                if 'PWDR' in histogram:
2801                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
2802                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
2803                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
2804                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
2805                   
2806                    if parmDict[pfx+'LeBail']:
2807                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
2808                    else:
2809                        if pfx+'Scale' in PhFrExtPOSig:
2810                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2811                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2812                            pFile.write(' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f\n'%
2813                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
2814                        if pfx+'Extinction' in PhFrExtPOSig:
2815                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
2816                        if hapData['Pref.Ori.'][0] == 'MD':
2817                            if pfx+'MD' in PhFrExtPOSig:
2818                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
2819                        else:
2820                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2821                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2822                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2823                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2824                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
2825                        if len(BabSig):
2826                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2827                   
2828                elif 'HKLF' in histogram:
2829                    Inst = Histogram['Instrument Parameters'][0]
2830                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
2831                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2832                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
2833                    if FFtables != None and 'N' not in Inst['Type'][0]:
2834                        PrintFprime(FFtables,hfx,pFile)
2835                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
2836                    if pfx+'Scale' in ScalExtSig:
2837                        pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
2838                    if hapData['Extinction'][0] != 'None':
2839                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2840                    if len(BabSig):
2841                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2842                    if pfx+'Flack' in ScalExtSig:
2843                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
2844                    if len(TwinFrSig):
2845                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2846
2847################################################################################
2848##### Histogram data
2849################################################################################       
2850                   
2851def GetHistogramData(Histograms,Print=True,pFile=None):
2852    'needs a doc string'
2853   
2854    def GetBackgroundParms(hId,Background):
2855        Back = Background[0]
2856        DebyePeaks = Background[1]
2857        bakType,bakFlag = Back[:2]
2858        backVals = Back[3:]
2859        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2860        backDict = dict(zip(backNames,backVals))
2861        backVary = []
2862        if bakFlag:
2863            backVary = backNames
2864        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2865        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2866        debyeDict = {}
2867        debyeList = []
2868        for i in range(DebyePeaks['nDebye']):
2869            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2870            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2871            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2872        debyeVary = []
2873        for item in debyeList:
2874            if item[1]:
2875                debyeVary.append(item[0])
2876        backDict.update(debyeDict)
2877        backVary += debyeVary
2878        peakDict = {}
2879        peakList = []
2880        for i in range(DebyePeaks['nPeaks']):
2881            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2882                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2883            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2884            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2885        peakVary = []
2886        for item in peakList:
2887            if item[1]:
2888                peakVary.append(item[0])
2889        backDict.update(peakDict)
2890        backVary += peakVary
2891        return bakType,backDict,backVary           
2892       
2893    def GetInstParms(hId,Inst):
2894        #patch
2895        if 'Z' not in Inst:
2896            Inst['Z'] = [0.0,0.0,False]
2897        dataType = Inst['Type'][0]
2898        instDict = {}
2899        insVary = []
2900        pfx = ':'+str(hId)+':'
2901        insKeys = list(Inst.keys())
2902        insKeys.sort()
2903        for item in insKeys:
2904            insName = pfx+item
2905            instDict[insName] = Inst[item][1]
2906            if len(Inst[item]) > 2 and Inst[item][2]:
2907                insVary.append(insName)
2908        if 'C' in dataType:
2909            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2910        elif 'T' in dataType:   #trap zero alp, bet coeff.
2911            if not instDict[pfx+'alpha']:
2912                instDict[pfx+'alpha'] = 1.0
2913            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
2914                instDict[pfx+'beta-1'] = 1.0
2915        return dataType,instDict,insVary
2916       
2917    def GetSampleParms(hId,Sample):
2918        sampVary = []
2919        hfx = ':'+str(hId)+':'       
2920        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2921            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2922        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2923            if key in Sample:
2924                sampDict[hfx+key] = Sample[key]
2925        Type = Sample['Type']
2926        if 'Bragg' in Type:             #Bragg-Brentano
2927            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2928                sampDict[hfx+item] = Sample[item][0]
2929                if Sample[item][1]:
2930                    sampVary.append(hfx+item)
2931        elif 'Debye' in Type:        #Debye-Scherrer
2932            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2933                sampDict[hfx+item] = Sample[item][0]
2934                if Sample[item][1]:
2935                    sampVary.append(hfx+item)
2936        return Type,sampDict,sampVary
2937       
2938    def PrintBackground(Background):
2939        Back = Background[0]
2940        DebyePeaks = Background[1]
2941        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
2942        line = ' Coefficients: '
2943        for i,back in enumerate(Back[3:]):
2944            line += '%10.3f'%(back)
2945            if i and not i%10:
2946                line += '\n'+15*' '
2947        pFile.write(line+'\n')
2948        if DebyePeaks['nDebye']:
2949            pFile.write('\n Debye diffuse scattering coefficients\n')
2950            parms = ['DebyeA','DebyeR','DebyeU']
2951            line = ' names :  '
2952            for parm in parms:
2953                line += '%8s refine?'%(parm)
2954            pFile.write(line+'\n')
2955            for j,term in enumerate(DebyePeaks['debyeTerms']):
2956                line = ' term'+'%2d'%(j)+':'
2957                for i in range(3):
2958                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2959                pFile.write(line+'\n')
2960        if DebyePeaks['nPeaks']:
2961            pFile.write('\n Single peak coefficients\n')
2962            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2963            line = ' names :  '
2964            for parm in parms:
2965                line += '%8s refine?'%(parm)
2966            pFile.write(line+'\n')
2967            for j,term in enumerate(DebyePeaks['peaksList']):
2968                line = ' peak'+'%2d'%(j)+':'
2969                for i in range(4):
2970                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2971                pFile.write(line+'\n')
2972       
2973    def PrintInstParms(Inst):
2974        pFile.write('\n Instrument Parameters:\n')
2975        insKeys = list(Inst.keys())
2976        insKeys.sort()
2977        iBeg = 0
2978        Ok = True
2979        while Ok:
2980            ptlbls = ' name  :'
2981            ptstr =  ' value :'
2982            varstr = ' refine:'
2983            iFin = min(iBeg+9,len(insKeys))
2984            for item in insKeys[iBeg:iFin]:
2985                if item not in ['Type','Source']:
2986                    ptlbls += '%12s' % (item)
2987                    ptstr += '%12.6f' % (Inst[item][1])
2988                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2989                        varstr += 12*' '
2990                    else:
2991                        varstr += '%12s' % (str(bool(Inst[item][2])))
2992            pFile.write(ptlbls+'\n')
2993            pFile.write(ptstr+'\n')
2994            pFile.write(varstr+'\n')
2995            iBeg = iFin
2996            if iBeg == len(insKeys):
2997                Ok = False
2998            else:
2999                pFile.write('\n')
3000       
3001    def PrintSampleParms(Sample):
3002        pFile.write('\n Sample Parameters:\n')
3003        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
3004            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
3005        ptlbls = ' name  :'
3006        ptstr =  ' value :'
3007        varstr = ' refine:'
3008        if 'Bragg' in Sample['Type']:
3009            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
3010                ptlbls += '%14s'%(item)
3011                ptstr += '%14.4f'%(Sample[item][0])
3012                varstr += '%14s'%(str(bool(Sample[item][1])))
3013           
3014        elif 'Debye' in Type:        #Debye-Scherrer
3015            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
3016                ptlbls += '%14s'%(item)
3017                ptstr += '%14.4f'%(Sample[item][0])
3018                varstr += '%14s'%(str(bool(Sample[item][1])))
3019
3020        pFile.write(ptlbls+'\n')
3021        pFile.write(ptstr+'\n')
3022        pFile.write(varstr+'\n')
3023       
3024    histDict = {}
3025    histVary = []
3026    controlDict = {}
3027    histoList = list(Histograms.keys())
3028    histoList.sort()
3029    for histogram in histoList:
3030        if 'PWDR' in histogram:
3031            Histogram = Histograms[histogram]
3032            hId = Histogram['hId']
3033            pfx = ':'+str(hId)+':'
3034            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3035            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3036            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3037            for excl in controlDict[pfx+'Exclude']:
3038                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3039            if controlDict[pfx+'Exclude']:
3040                ma.mask_rows(Histogram['Data'])
3041            Background = Histogram['Background']
3042            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3043            controlDict[pfx+'bakType'] = Type
3044            histDict.update(bakDict)
3045            histVary += bakVary
3046           
3047            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3048            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3049                print ('Warning: tabulated profile coefficients are ignored')
3050            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3051            controlDict[pfx+'histType'] = Type
3052            if 'XC' in Type:
3053                if pfx+'Lam1' in instDict:
3054                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3055                else:
3056                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3057            histDict.update(instDict)
3058            histVary += insVary
3059           
3060            Sample = Histogram['Sample Parameters']
3061            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3062            controlDict[pfx+'instType'] = Type
3063            histDict.update(sampDict)
3064            histVary += sampVary
3065           
3066   
3067            if Print: 
3068                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3069                pFile.write(135*'='+'\n')
3070                Units = {'C':' deg','T':' msec'}
3071                units = Units[controlDict[pfx+'histType'][2]]
3072                Limits = controlDict[pfx+'Limits']
3073                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3074                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3075                if len(controlDict[pfx+'Exclude']):
3076                    excls = controlDict[pfx+'Exclude']
3077                    for excl in excls:
3078                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3079                PrintSampleParms(Sample)
3080                PrintInstParms(Inst[0])
3081                PrintBackground(Background)
3082        elif 'HKLF' in histogram:
3083            Histogram = Histograms[histogram]
3084            hId = Histogram['hId']
3085            pfx = ':'+str(hId)+':'
3086            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3087            Inst = Histogram['Instrument Parameters'][0]
3088            controlDict[pfx+'histType'] = Inst['Type'][0]
3089            if 'X' in Inst['Type'][0]:
3090                histDict[pfx+'Lam'] = Inst['Lam'][1]
3091                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
3092    return histVary,histDict,controlDict
3093   
3094def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3095    'needs a doc string'
3096   
3097    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3098        Back = Background[0]
3099        DebyePeaks = Background[1]
3100        lenBack = len(Back[3:])
3101        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3102        for i in range(lenBack):
3103            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3104            if pfx+'Back;'+str(i) in sigDict:
3105                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3106        if DebyePeaks['nDebye']:
3107            for i in range(DebyePeaks['nDebye']):
3108                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3109                for j,name in enumerate(names):
3110                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3111                    if name in sigDict:
3112                        backSig[lenBack+3*i+j] = sigDict[name]           
3113        if DebyePeaks['nPeaks']:
3114            for i in range(DebyePeaks['nPeaks']):
3115                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3116                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3117                for j,name in enumerate(names):
3118                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3119                    if name in sigDict:
3120                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3121        return backSig
3122       
3123    def SetInstParms(pfx,Inst,parmDict,sigDict):
3124        instSig = {}
3125        insKeys = list(Inst.keys())
3126        insKeys.sort()
3127        for item in insKeys:
3128            insName = pfx+item
3129            Inst[item][1] = parmDict[insName]
3130            if insName in sigDict:
3131                instSig[item] = sigDict[insName]
3132            else:
3133                instSig[item] = 0
3134        return instSig
3135       
3136    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3137        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3138            sampSig = [0 for i in range(5)]
3139            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3140                Sample[item][0] = parmDict[pfx+item]
3141                if pfx+item in sigDict:
3142                    sampSig[i] = sigDict[pfx+item]
3143        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3144            sampSig = [0 for i in range(4)]
3145            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3146                Sample[item][0] = parmDict[pfx+item]
3147                if pfx+item in sigDict:
3148                    sampSig[i] = sigDict[pfx+item]
3149        return sampSig
3150       
3151    def PrintBackgroundSig(Background,backSig):
3152        Back = Background[0]
3153        DebyePeaks = Background[1]
3154        valstr = ' value : '
3155        sigstr = ' sig   : '
3156        refine = False
3157        for i,back in enumerate(Back[3:]):
3158            valstr += '%10.4g'%(back)
3159            if Back[1]:
3160                refine = True
3161                sigstr += '%10.4g'%(backSig[i])
3162            else:
3163                sigstr += 10*' '
3164        if refine:
3165            pFile.write('\n Background function: %s\n'%Back[0])
3166            pFile.write(valstr+'\n')
3167            pFile.write(sigstr+'\n')
3168        if DebyePeaks['nDebye']:
3169            ifAny = False
3170            ptfmt = "%12.3f"
3171            names =  ' names :'
3172            ptstr =  ' values:'
3173            sigstr = ' esds  :'
3174            for item in sigDict:
3175                if 'Debye' in item:
3176                    ifAny = True
3177                    names += '%12s'%(item)
3178                    ptstr += ptfmt%(parmDict[item])
3179                    sigstr += ptfmt%(sigDict[item])
3180            if ifAny:
3181                pFile.write('\n Debye diffuse scattering coefficients\n')
3182                pFile.write(names+'\n')
3183                pFile.write(ptstr+'\n')
3184                pFile.write(sigstr+'\n')
3185        if DebyePeaks['nPeaks']:
3186            pFile.write('\n Single peak coefficients:\n')
3187            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3188            line = ' peak no. '
3189            for parm in parms:
3190                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3191            pFile.write(line+'\n')
3192            for ip in range(DebyePeaks['nPeaks']):
3193                ptstr = ' %4d '%(ip)
3194                for parm in parms:
3195                    fmt = '%14.3f'
3196                    efmt = '%12.3f'
3197                    if parm == 'BkPkpos':
3198                        fmt = '%14.4f'
3199                        efmt = '%12.4f'
3200                    name = pfx+parm+';%d'%(ip)
3201                    ptstr += fmt%(parmDict[name])
3202                    if name in sigDict:
3203                        ptstr += efmt%(sigDict[name])
3204                    else:
3205                        ptstr += 12*' '
3206                pFile.write(ptstr+'\n')
3207        sumBk = np.array(Histogram['sumBk'])
3208        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3209            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3210       
3211    def PrintInstParmsSig(Inst,instSig):
3212        refine = False
3213        insKeys = list(instSig.keys())
3214        insKeys.sort()
3215        iBeg = 0
3216        Ok = True
3217        while Ok:
3218            ptlbls = ' names :'
3219            ptstr =  ' value :'
3220            sigstr = ' sig   :'
3221            iFin = min(iBeg+9,len(insKeys))
3222            for name in insKeys[iBeg:iFin]:
3223                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3224                    ptlbls += '%12s' % (name)
3225                    ptstr += '%12.6f' % (Inst[name][1])
3226                    if instSig[name]:
3227                        refine = True
3228                        sigstr += '%12.6f' % (instSig[name])
3229                    else:
3230                        sigstr += 12*' '
3231            if refine:
3232                pFile.write('\n Instrument Parameters:\n')
3233                pFile.write(ptlbls+'\n')
3234                pFile.write(ptstr+'\n')
3235                pFile.write(sigstr+'\n')
3236            iBeg = iFin
3237            if iBeg == len(insKeys):
3238                Ok = False
3239       
3240    def PrintSampleParmsSig(Sample,sampleSig):
3241        ptlbls = ' names :'
3242        ptstr =  ' values:'
3243        sigstr = ' sig   :'
3244        refine = False
3245        if 'Bragg' in Sample['Type']:
3246            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3247                ptlbls += '%14s'%(item)
3248                ptstr += '%14.4f'%(Sample[item][0])
3249                if sampleSig[i]:
3250                    refine = True
3251                    sigstr += '%14.4f'%(sampleSig[i])
3252                else:
3253                    sigstr += 14*' '
3254           
3255        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3256            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3257                ptlbls += '%14s'%(item)
3258                ptstr += '%14.4f'%(Sample[item][0])
3259                if sampleSig[i]:
3260                    refine = True
3261                    sigstr += '%14.4f'%(sampleSig[i])
3262                else:
3263                    sigstr += 14*' '
3264
3265        if refine:
3266            pFile.write('\n Sample Parameters:\n')
3267            pFile.write(ptlbls+'\n')
3268            pFile.write(ptstr+'\n')
3269            pFile.write(sigstr+'\n')
3270       
3271    histoList = list(Histograms.keys())
3272    histoList.sort()
3273    for histogram in histoList:
3274        if 'PWDR' in histogram:
3275            Histogram = Histograms[histogram]
3276            hId = Histogram['hId']
3277            pfx = ':'+str(hId)+':'
3278            Background = Histogram['Background']
3279            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3280           
3281            Inst = Histogram['Instrument Parameters'][0]
3282            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3283       
3284            Sample = Histogram['Sample Parameters']
3285            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3286
3287            pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3288            pFile.write(135*'='+'\n')
3289            pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3290            pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3291                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3292            pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3293                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3294            if Print:
3295                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3296                if FFtables != None and 'N' not in Inst['Type'][0]:
3297                    PrintFprime(FFtables,pfx,pFile)
3298                PrintSampleParmsSig(Sample,sampSig)
3299                PrintInstParmsSig(Inst,instSig)
3300                PrintBackgroundSig(Background,backSig)
3301               
Note: See TracBrowser for help on using the repository browser.