source: trunk/GSASIIstrIO.py @ 3438

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

g2strIO
fix magnetic hkl extinctions in reflection generation
remove Uniq & Phi - never used
change all config 'debug' out put to have 'DBG' in message - many places
fix phoenix decode problem for Comments
Display transformation math on Transform dialog box & labels on matrix & vectors
use a general exception in getSelection (G2plot line 7141) to avoid crashes

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