source: trunk/GSASIIstrIO.py @ 3711

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

major constraint update: move conflicting equivs to constraints; allow a formula for multiplier; update docs extensively. New routine EvaluateMultipliers? needed to evaluate formulae

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