source: trunk/GSASIIstrIO.py @ 3266

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

change modulated wave type assignments so one for each kind of wave (pos, frac, etc).
implement optional display of PNT data set positions on pole figure plots. Picker displays histo. name.

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