source: trunk/GSASIIstrIO.py @ 3157

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

add Lorentzian term 'Z' to CW profile instrument parameters - constant term
change sig-q function to d*sig-q from sig-q/d2 for TOF profile instrument parameters
add filter to SumDialog? for PWDR & IMG data

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 148.4 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2017-11-21 19:52:20 +0000 (Tue, 21 Nov 2017) $
4# $Author: vondreele $
5# $Revision: 3157 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 3157 2017-11-21 19:52:20Z 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: 3157 $")
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 'SL' in bldata[1]:
863                blen = bldata[1]['SL'][0]
864                bres = []
865            else:
866                blen = 0
867                bres = bldata[1]['BW-LS']
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,AtomSS['waveType']))
990                for iw,wave in enumerate(Waves):                   
991                    line = ''
992                    if AtomSS['waveType'] in ['Block','ZigZag'] and Stype == 'Spos' and not iw:
993                        for item in names[Stype][6:]:
994                            line += '%8s '%(item)                       
995                    else:
996                        if Stype == 'Spos':
997                            for item in names[Stype][:6]:
998                                line += '%8s '%(item)
999                        else:
1000                            for item in names[Stype]:
1001                                line += '%8s '%(item)
1002                    pFile.write(line+'\n')
1003                    line = ''
1004                    for item in wave[0]:
1005                        line += '%8.4f '%(item)
1006                    line += ' Refine? '+str(wave[1])
1007                    pFile.write(line+'\n')
1008       
1009    def PrintTexture(textureData):
1010        topstr = '\n Spherical harmonics texture: Order:' + \
1011            str(textureData['Order'])
1012        if textureData['Order']:
1013            pFile.write('%s Refine? %s\n'%(topstr,textureData['SH Coeff'][0]))
1014        else:
1015            pFile.write(topstr+'\n')
1016            return
1017        names = ['omega','chi','phi']
1018        line = '\n'
1019        for name in names:
1020            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
1021        pFile.write(line+'\n')
1022        pFile.write('\n Texture coefficients:\n')
1023        SHcoeff = textureData['SH Coeff'][1]
1024        SHkeys = list(SHcoeff.keys())
1025        nCoeff = len(SHcoeff)
1026        nBlock = nCoeff//10+1
1027        iBeg = 0
1028        iFin = min(iBeg+10,nCoeff)
1029        for block in range(nBlock):
1030            ptlbls = ' names :'
1031            ptstr =  ' values:'
1032            for item in SHkeys[iBeg:iFin]:
1033                ptlbls += '%12s'%(item)
1034                ptstr += '%12.4f'%(SHcoeff[item]) 
1035            pFile.write(ptlbls+'\n')
1036            pFile.write(ptstr+'\n')
1037            iBeg += 10
1038            iFin = min(iBeg+10,nCoeff)
1039       
1040    def MakeRBParms(rbKey,phaseVary,phaseDict):
1041        rbid = str(rbids.index(RB['RBId']))
1042        pfxRB = pfx+'RB'+rbKey+'P'
1043        pstr = ['x','y','z']
1044        ostr = ['a','i','j','k']
1045        for i in range(3):
1046            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
1047            phaseDict[name] = RB['Orig'][0][i]
1048            if RB['Orig'][1]:
1049                phaseVary += [name,]
1050        pfxRB = pfx+'RB'+rbKey+'O'
1051        for i in range(4):
1052            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
1053            phaseDict[name] = RB['Orient'][0][i]
1054            if RB['Orient'][1] == 'AV' and i:
1055                phaseVary += [name,]
1056            elif RB['Orient'][1] == 'A' and not i:
1057                phaseVary += [name,]
1058           
1059    def MakeRBThermals(rbKey,phaseVary,phaseDict):
1060        rbid = str(rbids.index(RB['RBId']))
1061        tlstr = ['11','22','33','12','13','23']
1062        sstr = ['12','13','21','23','31','32','AA','BB']
1063        if 'T' in RB['ThermalMotion'][0]:
1064            pfxRB = pfx+'RB'+rbKey+'T'
1065            for i in range(6):
1066                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1067                phaseDict[name] = RB['ThermalMotion'][1][i]
1068                if RB['ThermalMotion'][2][i]:
1069                    phaseVary += [name,]
1070        if 'L' in RB['ThermalMotion'][0]:
1071            pfxRB = pfx+'RB'+rbKey+'L'
1072            for i in range(6):
1073                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
1074                phaseDict[name] = RB['ThermalMotion'][1][i+6]
1075                if RB['ThermalMotion'][2][i+6]:
1076                    phaseVary += [name,]
1077        if 'S' in RB['ThermalMotion'][0]:
1078            pfxRB = pfx+'RB'+rbKey+'S'
1079            for i in range(8):
1080                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
1081                phaseDict[name] = RB['ThermalMotion'][1][i+12]
1082                if RB['ThermalMotion'][2][i+12]:
1083                    phaseVary += [name,]
1084        if 'U' in RB['ThermalMotion'][0]:
1085            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
1086            phaseDict[name] = RB['ThermalMotion'][1][0]
1087            if RB['ThermalMotion'][2][0]:
1088                phaseVary += [name,]
1089               
1090    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
1091        rbid = str(rbids.index(RB['RBId']))
1092        pfxRB = pfx+'RB'+rbKey+'Tr;'
1093        for i,torsion in enumerate(RB['Torsions']):
1094            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
1095            phaseDict[name] = torsion[0]
1096            if torsion[1]:
1097                phaseVary += [name,]
1098                   
1099    if Print:
1100        pFile.write('\n Phases:\n')
1101    phaseVary = []
1102    phaseDict = {}
1103    pawleyLookup = {}
1104    FFtables = {}                   #scattering factors - xrays
1105    MFtables = {}                   #Mag. form factors
1106    BLtables = {}                   # neutrons
1107    Natoms = {}
1108    maxSSwave = {}
1109    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1110    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1111    atomIndx = {}
1112    for name in PhaseData:
1113        General = PhaseData[name]['General']
1114        pId = PhaseData[name]['pId']
1115        pfx = str(pId)+'::'
1116        FFtable = G2el.GetFFtable(General['AtomTypes'])
1117        BLtable = G2el.GetBLtable(General)
1118        FFtables.update(FFtable)
1119        BLtables.update(BLtable)
1120        phaseDict[pfx+'isMag'] = False
1121        SGData = General['SGData']
1122        SGtext,SGtable = G2spc.SGPrint(SGData)
1123        if General['Type'] == 'magnetic':
1124            MFtable = G2el.GetMFtable(General['AtomTypes'],General['Lande g'])
1125            MFtables.update(MFtable)
1126            phaseDict[pfx+'isMag'] = True
1127            SpnFlp = SGData['SpnFlp']
1128        Atoms = PhaseData[name]['Atoms']
1129        if Atoms and not General.get('doPawley'):
1130            cx,ct,cs,cia = General['AtomPtrs']
1131            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1132        PawleyRef = PhaseData[name].get('Pawley ref',[])
1133        cell = General['Cell']
1134        A = G2lat.cell2A(cell[1:7])
1135        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
1136            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
1137        if cell[0]:
1138            phaseVary += cellVary(pfx,SGData)       #also fills in symmetry required constraints
1139        SSGtext = []    #no superstructure
1140        im = 0
1141        if General.get('Modulated',False):
1142            im = 1      #refl offset
1143            Vec,vRef,maxH = General['SuperVec']
1144            phaseDict.update({pfx+'mV0':Vec[0],pfx+'mV1':Vec[1],pfx+'mV2':Vec[2]})
1145            SSGData = General['SSGData']
1146            SSGtext,SSGtable = G2spc.SSGPrint(SGData,SSGData)
1147            if vRef:
1148                phaseVary += modVary(pfx,SSGData)       
1149        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
1150        if resRBData:
1151            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
1152            for iRB,RB in enumerate(resRBData):
1153                MakeRBParms('R',phaseVary,phaseDict)
1154                MakeRBThermals('R',phaseVary,phaseDict)
1155                MakeRBTorsions('R',phaseVary,phaseDict)
1156       
1157        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
1158        if vecRBData:
1159            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
1160            for iRB,RB in enumerate(vecRBData):
1161                MakeRBParms('V',phaseVary,phaseDict)
1162                MakeRBThermals('V',phaseVary,phaseDict)
1163                   
1164        Natoms[pfx] = 0
1165        maxSSwave[pfx] = {'Sfrac':0,'Spos':0,'Sadp':0,'Smag':0}
1166        if Atoms and not General.get('doPawley'):
1167            cx,ct,cs,cia = General['AtomPtrs']
1168            Natoms[pfx] = len(Atoms)
1169            for i,at in enumerate(Atoms):
1170                atomIndx[at[cia+8]] = [pfx,i]      #lookup table for restraints
1171                phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
1172                    pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
1173                    pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
1174                    pfx+'AI/A:'+str(i):at[cia],})
1175                if at[cia] == 'I':
1176                    phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
1177                else:
1178                    phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
1179                        pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
1180                if General['Type'] == 'magnetic':
1181                    phaseDict.update({pfx+'AMx:'+str(i):at[cx+4],pfx+'AMy:'+str(i):at[cx+5],pfx+'AMz:'+str(i):at[cx+6]})
1182                if 'F' in at[ct+1]:
1183                    phaseVary.append(pfx+'Afrac:'+str(i))
1184                if 'X' in at[ct+1]:
1185                    try:    #patch for sytsym name changes
1186                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1187                    except KeyError:
1188                        Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1189                        at[cs] = Sytsym
1190                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1191                    names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1192                    equivs = [[],[],[]]
1193                    for j in range(3):
1194                        if xId[j] > 0:                               
1195                            phaseVary.append(names[j])
1196                            equivs[xId[j]-1].append([names[j],xCoef[j]])
1197                    for equiv in equivs:
1198                        if len(equiv) > 1:
1199                            name = equiv[0][0]
1200                            coef = equiv[0][1]
1201                            for eqv in equiv[1:]:
1202                                eqv[1] /= coef
1203                                G2mv.StoreEquivalence(name,(eqv,))
1204                if 'U' in at[ct+1]:
1205                    if at[cia] == 'I':
1206                        phaseVary.append(pfx+'AUiso:'+str(i))
1207                    else:
1208                        try:    #patch for sytsym name changes
1209                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1210                        except KeyError:
1211                            Sytsym = G2spc.SytSym(at[cx:cx+3],SGData)[0]
1212                            at[cs] = Sytsym
1213                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1214                        names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1215                            pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1216                        equivs = [[],[],[],[],[],[]]
1217                        for j in range(6):
1218                            if uId[j] > 0:                               
1219                                phaseVary.append(names[j])
1220                                equivs[uId[j]-1].append([names[j],uCoef[j]])
1221                        for equiv in equivs:
1222                            if len(equiv) > 1:
1223                                name = equiv[0][0]
1224                                coef = equiv[0][1]
1225                                for eqv in equiv[1:]:
1226                                    eqv[1] /= coef
1227                                    G2mv.StoreEquivalence(name,(eqv,))
1228                if 'M' in at[ct+1]:
1229                    SytSym,Mul,Nop,dupDir = G2spc.SytSym(at[cx:cx+3],SGData)
1230                    mId,mCoef = G2spc.GetCSpqinel(SytSym,SpnFlp,dupDir)
1231                    names = [pfx+'AMx:'+str(i),pfx+'AMy:'+str(i),pfx+'AMz:'+str(i)]
1232                    equivs = [[],[],[]]
1233                    for j in range(3):
1234                        if mId[j] > 0:
1235                            phaseVary.append(names[j])
1236                            equivs[mId[j]-1].append([names[j],mCoef[j]])
1237                    for equiv in equivs:
1238                        if len(equiv) > 1:
1239                            name = equiv[0][0]
1240                            coef = equiv[0][1]
1241                            for eqv in equiv[1:]:
1242                                eqv[1] /= coef
1243                                G2mv.StoreEquivalence(name,(eqv,))
1244                if General.get('Modulated',False):
1245                    AtomSS = at[-1]['SS1']
1246                    waveType = AtomSS['waveType']
1247                    phaseDict[pfx+'waveType:'+str(i)] = waveType
1248                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
1249                        Waves = AtomSS[Stype]
1250                        nx = 0
1251                        for iw,wave in enumerate(Waves):
1252                            if not iw:
1253                                if waveType in ['ZigZag','Block']:
1254                                    nx = 1
1255                                CSI = G2spc.GetSSfxuinel(waveType,1,at[cx:cx+3],SGData,SSGData)
1256                            else:
1257                                CSI = G2spc.GetSSfxuinel('Fourier',iw+1-nx,at[cx:cx+3],SGData,SSGData)
1258                            uId,uCoef = CSI[Stype]
1259                            stiw = str(i)+':'+str(iw)
1260                            if Stype == 'Spos':
1261                                if waveType in ['ZigZag','Block',] and not iw:
1262                                    names = [pfx+'Tmin:'+stiw,pfx+'Tmax:'+stiw,pfx+'Xmax:'+stiw,pfx+'Ymax:'+stiw,pfx+'Zmax:'+stiw]
1263                                    equivs = [[],[], [],[],[]]
1264                                else:
1265                                    names = [pfx+'Xsin:'+stiw,pfx+'Ysin:'+stiw,pfx+'Zsin:'+stiw,
1266                                        pfx+'Xcos:'+stiw,pfx+'Ycos:'+stiw,pfx+'Zcos:'+stiw]
1267                                    equivs = [[],[],[], [],[],[]]
1268                            elif Stype == 'Sadp':
1269                                names = [pfx+'U11sin:'+stiw,pfx+'U22sin:'+stiw,pfx+'U33sin:'+stiw,
1270                                    pfx+'U12sin:'+stiw,pfx+'U13sin:'+stiw,pfx+'U23sin:'+stiw,
1271                                    pfx+'U11cos:'+stiw,pfx+'U22cos:'+stiw,pfx+'U33cos:'+stiw,
1272                                    pfx+'U12cos:'+stiw,pfx+'U13cos:'+stiw,pfx+'U23cos:'+stiw]
1273                                equivs = [[],[],[],[],[],[], [],[],[],[],[],[]]
1274                            elif Stype == 'Sfrac':
1275                                equivs = [[],[]]
1276                                if 'Crenel' in waveType and not iw:
1277                                    names = [pfx+'Fzero:'+stiw,pfx+'Fwid:'+stiw]
1278                                else:
1279                                    names = [pfx+'Fsin:'+stiw,pfx+'Fcos:'+stiw]
1280                            elif Stype == 'Smag':
1281                                equivs = [[],[],[], [],[],[]]
1282                                names = [pfx+'MXsin:'+stiw,pfx+'MYsin:'+stiw,pfx+'MZsin:'+stiw,
1283                                    pfx+'MXcos:'+stiw,pfx+'MYcos:'+stiw,pfx+'MZcos:'+stiw]
1284                            phaseDict.update(dict(zip(names,wave[0])))
1285                            if wave[1]: #what do we do here for multiple terms in modulation constraints?
1286                                for j in range(len(equivs)):
1287                                    if uId[j][0] > 0:                               
1288                                        phaseVary.append(names[j])
1289                                        equivs[uId[j][0]-1].append([names[j],uCoef[j][0]])
1290                                for equiv in equivs:
1291                                    if len(equiv) > 1:
1292                                        name = equiv[0][0]
1293                                        coef = equiv[0][1]
1294                                        for eqv in equiv[1:]:
1295                                            eqv[1] /= coef
1296                                            G2mv.StoreEquivalence(name,(eqv,))
1297                            maxSSwave[pfx][Stype] = max(maxSSwave[pfx][Stype],iw+1)
1298            textureData = General['SH Texture']
1299            if textureData['Order'] and not seqRef:
1300                phaseDict[pfx+'SHorder'] = textureData['Order']
1301                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1302                for item in ['omega','chi','phi']:
1303                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1304                    if textureData['Sample '+item][0]:
1305                        phaseVary.append(pfx+'SH '+item)
1306                for item in textureData['SH Coeff'][1]:
1307                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1308                    if textureData['SH Coeff'][0]:
1309                        phaseVary.append(pfx+item)
1310               
1311            if Print:
1312                pFile.write('\n Phase name: %s\n'%General['Name'])
1313                pFile.write(135*'='+'\n')
1314                PrintFFtable(FFtable)
1315                PrintBLtable(BLtable)
1316                if General['Type'] == 'magnetic':
1317                    PrintMFtable(MFtable)
1318                pFile.write('\n')
1319                #how do we print magnetic symmetry table? TBD
1320                if len(SSGtext):    #if superstructure
1321                    for line in SSGtext: pFile.write(line+'\n')
1322                    if len(SSGtable):
1323                        for item in SSGtable:
1324                            line = ' %s '%(item)
1325                            pFile.write(line+'\n') 
1326                    else:
1327                        pFile.write(' ( 1)    %s\n'%(SSGtable[0]))
1328                else:
1329                    for line in SGtext: pFile.write(line+'\n')
1330                    if len(SGtable):
1331                        for item in SGtable:
1332                            line = ' %s '%(item)
1333                            pFile.write(line+'\n') 
1334                    else:
1335                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1336                PrintRBObjects(resRBData,vecRBData)
1337                PrintAtoms(General,Atoms)
1338                if General['Type'] == 'magnetic':
1339                    PrintMoments(General,Atoms)
1340                if General.get('Modulated',False):
1341                    PrintWaves(General,Atoms)
1342                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1343                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1344                if len(SSGtext):    #if superstructure
1345                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1346                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1347                if not seqRef:
1348                    PrintTexture(textureData)
1349                if name in RestraintDict:
1350                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1351                        textureData,RestraintDict[name],pFile)
1352                   
1353        elif PawleyRef:
1354            if Print:
1355                pFile.write('\n Phase name: %s\n'%General['Name'])
1356                pFile.write(135*'='+'\n')
1357                pFile.write('\n')
1358                if len(SSGtext):    #if superstructure
1359                    for line in SSGtext: pFile.write(line+'\n')
1360                    if len(SSGtable):
1361                        for item in SSGtable:
1362                            line = ' %s '%(item)
1363                            pFile.write(line+'\n') 
1364                    else:
1365                        pFile.write(' ( 1)    %s\n'%SSGtable[0])
1366                else:
1367                    for line in SGtext: pFile.write(line+'\n')
1368                    if len(SGtable):
1369                        for item in SGtable:
1370                            line = ' %s '%(item)
1371                            pFile.write(line+'\n')
1372                    else:
1373                        pFile.write(' ( 1)    %s\n'%(SGtable[0]))
1374                pFile.write('\n Unit cell: a = %.5f b = %.5f c = %.5f alpha = %.3f beta = %.3f gamma = %.3f volume = %.3f Refine? %s\n'%
1375                    (cell[1],cell[2],cell[3],cell[4],cell[5],cell[6],cell[7],cell[0]))
1376                if len(SSGtext):    #if superstructure
1377                    pFile.write('\n Modulation vector: mV0 = %.4f mV1 = %.4f mV2 = %.4f max mod. index = %d Refine? %s\n'%
1378                        (Vec[0],Vec[1],Vec[2],maxH,vRef))
1379            pawleyVary = []
1380            for i,refl in enumerate(PawleyRef):
1381                phaseDict[pfx+'PWLref:'+str(i)] = refl[6+im]
1382                if im:
1383                    pawleyLookup[pfx+'%d,%d,%d,%d'%(refl[0],refl[1],refl[2],refl[3])] = i
1384                else:
1385                    pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1386                if refl[5+im]:
1387                    pawleyVary.append(pfx+'PWLref:'+str(i))
1388            GetPawleyConstr(SGData['SGLaue'],PawleyRef,im,pawleyVary)      #does G2mv.StoreEquivalence
1389            phaseVary += pawleyVary
1390               
1391    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables,MFtables,maxSSwave
1392   
1393def cellFill(pfx,SGData,parmDict,sigDict): 
1394    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1395    from the parameter and sig dictionaries.
1396
1397    :param str pfx: parameter prefix ("n::", where n is a phase number)
1398    :param dict SGdata: a symmetry object
1399    :param dict parmDict: a dictionary of parameters
1400    :param dict sigDict:  a dictionary of uncertainties on parameters
1401
1402    :returns: A,sigA where each is a list of six terms with the A terms
1403    '''
1404    if SGData['SGLaue'] in ['-1',]:
1405        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1406            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1407    elif SGData['SGLaue'] in ['2/m',]:
1408        if SGData['SGUniq'] == 'a':
1409            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1410                0,0,parmDict[pfx+'A5']]
1411        elif SGData['SGUniq'] == 'b':
1412            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1413                0,parmDict[pfx+'A4'],0]
1414        else:
1415            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1416                parmDict[pfx+'A3'],0,0]
1417    elif SGData['SGLaue'] in ['mmm',]:
1418        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1419    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1420        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1421    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1422        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1423            parmDict[pfx+'A0'],0,0]
1424    elif SGData['SGLaue'] in ['3R', '3mR']:
1425        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1426            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1427    elif SGData['SGLaue'] in ['m3m','m3']:
1428        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1429
1430    try:
1431        if SGData['SGLaue'] in ['-1',]:
1432            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1433                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1434        elif SGData['SGLaue'] in ['2/m',]:
1435            if SGData['SGUniq'] == 'a':
1436                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1437                    0,0,sigDict[pfx+'A5']]
1438            elif SGData['SGUniq'] == 'b':
1439                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1440                    0,sigDict[pfx+'A4'],0]
1441            else:
1442                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1443                    sigDict[pfx+'A3'],0,0]
1444        elif SGData['SGLaue'] in ['mmm',]:
1445            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1446        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1447            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1448        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1449            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1450        elif SGData['SGLaue'] in ['3R', '3mR']:
1451            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1452        elif SGData['SGLaue'] in ['m3m','m3']:
1453            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1454    except KeyError:
1455        sigA = [0,0,0,0,0,0]
1456    return A,sigA
1457       
1458def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1459    'needs a doc string'
1460    if phaseRest:
1461        Amat = G2lat.cell2AB(cell)[0]
1462        cx,ct,cs = AtPtrs[:3]
1463        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1464            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1465            ['ChemComp','Sites'],['Texture','HKLs']]
1466        for name,rest in names:
1467            itemRest = phaseRest[name]
1468            if rest in itemRest and itemRest[rest] and itemRest['Use']:
1469                pFile.write('\n  %s restraint weight factor %10.3f Use: %s\n'%(name,itemRest['wtFactor'],str(itemRest['Use'])))
1470                if name in ['Bond','Angle','Plane','Chiral']:
1471                    pFile.write('     calc       obs      sig   delt/sig  atoms(symOp): \n')
1472                    for indx,ops,obs,esd in itemRest[rest]:
1473                        try:
1474                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1475                            AtName = ''
1476                            for i,Aname in enumerate(AtNames):
1477                                AtName += Aname
1478                                if ops[i] == '1':
1479                                    AtName += '-'
1480                                else:
1481                                    AtName += '+('+ops[i]+')-'
1482                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1483                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1484                            if name == 'Bond':
1485                                calc = G2mth.getRestDist(XYZ,Amat)
1486                            elif name == 'Angle':
1487                                calc = G2mth.getRestAngle(XYZ,Amat)
1488                            elif name == 'Plane':
1489                                calc = G2mth.getRestPlane(XYZ,Amat)
1490                            elif name == 'Chiral':
1491                                calc = G2mth.getRestChiral(XYZ,Amat)
1492                            pFile.write(' %9.3f %9.3f %8.3f %8.3f   %s\n'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1]))
1493                        except KeyError:
1494                            del itemRest[rest]
1495                elif name in ['Torsion','Rama']:
1496                    pFile.write('  atoms(symOp)  calc  obs  sig  delt/sig  torsions: \n')
1497                    coeffDict = itemRest['Coeff']
1498                    for indx,ops,cofName,esd in itemRest[rest]:
1499                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1500                        AtName = ''
1501                        for i,Aname in enumerate(AtNames):
1502                            AtName += Aname+'+('+ops[i]+')-'
1503                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1504                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1505                        if name == 'Torsion':
1506                            tor = G2mth.getRestTorsion(XYZ,Amat)
1507                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1508                            pFile.write(' %8.3f %8.3f %.3f %8.3f %8.3f %s\n'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1]))
1509                        else:
1510                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1511                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1512                            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]))
1513                elif name == 'ChemComp':
1514                    pFile.write('     atoms   mul*frac  factor     prod\n')
1515                    for indx,factors,obs,esd in itemRest[rest]:
1516                        try:
1517                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1518                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1519                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1520                            mulfrac = mul*frac
1521                            calcs = mul*frac*factors
1522                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1523                                pFile.write(' %10s %8.3f %8.3f %8.3f\n'%(atom,mf,fr,clc))
1524                            pFile.write(' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f\n'%(np.sum(calcs),obs,esd))
1525                        except KeyError:
1526                            del itemRest[rest]
1527                elif name == 'Texture' and textureData['Order']:
1528                    Start = False
1529                    SHCoef = textureData['SH Coeff'][1]
1530                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1531                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1532                    pFile.write ('    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd \n')
1533                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1534                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1535                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1536                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1537                        Z = ma.masked_greater(Z,0.0)
1538                        num = ma.count(Z)
1539                        sum = 0
1540                        if num:
1541                            sum = np.sum(Z)
1542                        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))
1543       
1544def getCellEsd(pfx,SGData,A,covData):
1545    'needs a doc string'
1546    rVsq = G2lat.calc_rVsq(A)
1547    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1548    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1549    varyList = covData['varyList']
1550    covMatrix = covData['covMatrix']
1551    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1552    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm']:
1553        vcov[1,1] = vcov[3,3] = vcov[0,1] = vcov[1,0] = vcov[0,0]
1554        vcov[1,3] = vcov[3,1] = vcov[0,3] = vcov[3,0] = vcov[0,0]
1555        vcov[1,2] = vcov[2,1] = vcov[2,3] = vcov[3,2] = vcov[0,2]
1556    elif SGData['SGLaue'] in ['m3','m3m']:
1557        vcov[0:3,0:3] = vcov[0,0]
1558    elif SGData['SGLaue'] in ['4/m', '4/mmm']:
1559        vcov[0:2,0:2] = vcov[0,0]
1560        vcov[1,2] = vcov[2,1] = vcov[0,2]
1561    elif SGData['SGLaue'] in ['3R','3mR']:
1562        vcov[0:3,0:3] = vcov[0,0]
1563#        vcov[4,4] = vcov[5,5] = vcov[3,3]
1564        vcov[3:6,3:6] = vcov[3,3]
1565        vcov[0:3,3:6] = vcov[0,3]
1566        vcov[3:6,0:3] = vcov[3,0]
1567    delt = 1.e-9
1568    drVdA = np.zeros(6)
1569    for i in range(6):
1570        A[i] += delt
1571        drVdA[i] = G2lat.calc_rVsq(A)
1572        A[i] -= 2*delt
1573        drVdA[i] -= G2lat.calc_rVsq(A)
1574        A[i] += delt
1575    drVdA /= 2.*delt   
1576    srcvlsq = np.inner(drVdA,np.inner(drVdA,vcov))
1577    Vol = 1/np.sqrt(rVsq)
1578    sigVol = Vol**3*np.sqrt(srcvlsq)/2.         #ok - checks with GSAS
1579   
1580    dcdA = np.zeros((6,6))
1581    for i in range(6):
1582        pdcdA =np.zeros(6)
1583        A[i] += delt
1584        pdcdA += G2lat.A2cell(A)
1585        A[i] -= 2*delt
1586        pdcdA -= G2lat.A2cell(A)
1587        A[i] += delt
1588        dcdA[i] = pdcdA/(2.*delt)
1589   
1590    sigMat = np.inner(dcdA,np.inner(dcdA,vcov))
1591    var = np.diag(sigMat)
1592    CS = np.where(var>0.,np.sqrt(var),0.)
1593    if SGData['SGLaue'] in ['3', '3m1', '31m', '6/m', '6/mmm','m3','m3m','4/m','4/mmm']:
1594        CS[3:6] = 0.0
1595    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]
1596   
1597def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1598    'needs a doc string'
1599   
1600    def PrintAtomsAndSig(General,Atoms,atomsSig):
1601        pFile.write('\n Atoms:\n')
1602        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1603        if General['Type'] == 'macromolecular':
1604            line = ' res no residue chain '+line
1605        cx,ct,cs,cia = General['AtomPtrs']
1606        pFile.write(line+'\n')
1607        pFile.write(135*'-'+'\n')
1608        fmt = {0:'%7s',ct:'%7s',cx:'%10.5f',cx+1:'%10.5f',cx+2:'%10.5f',cx+3:'%8.3f',cia+1:'%8.5f',
1609            cia+2:'%8.5f',cia+3:'%8.5f',cia+4:'%8.5f',cia+5:'%8.5f',cia+6:'%8.5f',cia+7:'%8.5f'}
1610        noFXsig = {cx:[10*' ','%10s'],cx+1:[10*' ','%10s'],cx+2:[10*' ','%10s'],cx+3:[8*' ','%8s']}
1611        for atyp in General['AtomTypes']:       #zero composition
1612            General['NoAtoms'][atyp] = 0.
1613        for i,at in enumerate(Atoms):
1614            General['NoAtoms'][at[ct]] += at[cx+3]*float(at[cx+5])     #new composition
1615            if General['Type'] == 'macromolecular':
1616                name = ' %s %s %s %s:'%(at[0],at[1],at[2],at[3])
1617                valstr = ' values:          '
1618                sigstr = ' sig   :          '
1619            else:
1620                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1621                valstr = ' values:'
1622                sigstr = ' sig   :'
1623            for ind in range(cx,cx+4):
1624                sigind = str(i)+':'+str(ind)
1625                valstr += fmt[ind]%(at[ind])                   
1626                if sigind in atomsSig:
1627                    sigstr += fmt[ind]%(atomsSig[sigind])
1628                else:
1629                    sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1630            if at[cia] == 'I':
1631                valstr += fmt[cia+1]%(at[cia+1])
1632                if '%d:%d'%(i,cia+1) in atomsSig:
1633                    sigstr += fmt[cia+1]%(atomsSig['%d:%d'%(i,cia+1)])
1634                else:
1635                    sigstr += 8*' '
1636            else:
1637                valstr += 8*' '
1638                sigstr += 8*' '
1639                for ind in range(cia+2,cia+8):
1640                    sigind = str(i)+':'+str(ind)
1641                    valstr += fmt[ind]%(at[ind])
1642                    if sigind in atomsSig:                       
1643                        sigstr += fmt[ind]%(atomsSig[sigind])
1644                    else:
1645                        sigstr += 8*' '
1646            pFile.write(name+'\n')
1647            pFile.write(valstr+'\n')
1648            pFile.write(sigstr+'\n')
1649           
1650    def PrintMomentsAndSig(General,Atoms,atomsSig):
1651        pFile.write('\n Magnetic Moments:\n')    #add magnitude & angle, etc.? TBD
1652        line = '   name      Mx        My        Mz'
1653        cx,ct,cs,cia = General['AtomPtrs']
1654        cmx = cx+4
1655        AtInfo = dict(zip(General['AtomTypes'],General['Lande g']))
1656        pFile.write(line+'\n')
1657        pFile.write(135*'-'+'\n')
1658        fmt = {0:'%7s',ct:'%7s',cmx:'%10.5f',cmx+1:'%10.5f',cmx+2:'%10.5f'}
1659        noFXsig = {cmx:[10*' ','%10s'],cmx+1:[10*' ','%10s'],cmx+2:[10*' ','%10s']}
1660        for i,at in enumerate(Atoms):
1661            if AtInfo[at[ct]]:
1662                name = fmt[0]%(at[ct-1])+fmt[1]%(at[ct])+':'
1663                valstr = ' values:'
1664                sigstr = ' sig   :'
1665                for ind in range(cmx,cmx+3):
1666                    sigind = str(i)+':'+str(ind)
1667                    valstr += fmt[ind]%(at[ind])                   
1668                    if sigind in atomsSig:
1669                        sigstr += fmt[ind]%(atomsSig[sigind])
1670                    else:
1671                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1672                pFile.write(name+'\n')
1673                pFile.write(valstr+'\n')
1674                pFile.write(sigstr+'\n')
1675           
1676    def PrintWavesAndSig(General,Atoms,wavesSig):
1677        cx,ct,cs,cia = General['AtomPtrs']
1678        pFile.write('\n Modulation waves\n')
1679        names = {'Sfrac':['Fsin','Fcos','Fzero','Fwid'],'Spos':['Xsin','Ysin','Zsin','Xcos','Ycos','Zcos','Tmin','Tmax','Xmax','Ymax','Zmax'],
1680            'Sadp':['U11sin','U22sin','U33sin','U12sin','U13sin','U23sin','U11cos','U22cos',
1681            'U33cos','U12cos','U13cos','U23cos'],'Smag':['MXsin','MYsin','MZsin','MXcos','MYcos','MZcos']}
1682        pFile.write(135*'-'+'\n')
1683        for i,at in enumerate(Atoms):
1684            AtomSS = at[-1]['SS1']
1685            waveType = AtomSS['waveType']
1686            for Stype in ['Sfrac','Spos','Sadp','Smag']:
1687                Waves = AtomSS[Stype]
1688                if len(Waves):
1689                    pFile.write(' atom: %s, site sym: %s, type: %s wave type: %s:\n'%
1690                        (at[ct-1],at[cs],Stype,waveType))
1691                    for iw,wave in enumerate(Waves):
1692                        stiw = ':'+str(i)+':'+str(iw)
1693                        namstr = '  names :'
1694                        valstr = '  values:'
1695                        sigstr = '  esds  :'
1696                        if Stype == 'Spos':
1697                            nt = 6
1698                            ot = 0
1699                            if waveType in ['ZigZag','Block',] and not iw:
1700                                nt = 5
1701                                ot = 6
1702                            for j in range(nt):
1703                                name = names['Spos'][j+ot]
1704                                namstr += '%12s'%(name)
1705                                valstr += '%12.4f'%(wave[0][j])
1706                                if name+stiw in wavesSig:
1707                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1708                                else:
1709                                    sigstr += 12*' '
1710                        elif Stype == 'Sfrac':
1711                            ot = 0
1712                            if 'Crenel' in waveType and not iw:
1713                                ot = 2
1714                            for j in range(2):
1715                                name = names['Sfrac'][j+ot]
1716                                namstr += '%12s'%(names['Sfrac'][j+ot])
1717                                valstr += '%12.4f'%(wave[0][j])
1718                                if name+stiw in wavesSig:
1719                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1720                                else:
1721                                    sigstr += 12*' '
1722                        elif Stype == 'Sadp':
1723                            for j in range(12):
1724                                name = names['Sadp'][j]
1725                                namstr += '%10s'%(names['Sadp'][j])
1726                                valstr += '%10.6f'%(wave[0][j])
1727                                if name+stiw in wavesSig:
1728                                    sigstr += '%10.6f'%(wavesSig[name+stiw])
1729                                else:
1730                                    sigstr += 10*' '
1731                        elif Stype == 'Smag':
1732                            for j in range(6):
1733                                name = names['Smag'][j]
1734                                namstr += '%12s'%(names['Smag'][j])
1735                                valstr += '%12.4f'%(wave[0][j])
1736                                if name+stiw in wavesSig:
1737                                    sigstr += '%12.4f'%(wavesSig[name+stiw])
1738                                else:
1739                                    sigstr += 12*' '
1740                               
1741                    pFile.write(namstr+'\n')
1742                    pFile.write(valstr+'\n')
1743                    pFile.write(sigstr+'\n')
1744       
1745               
1746    def PrintRBObjPOAndSig(rbfx,rbsx):
1747        namstr = '  names :'
1748        valstr = '  values:'
1749        sigstr = '  esds  :'
1750        for i,px in enumerate(['Px:','Py:','Pz:']):
1751            name = pfx+rbfx+px+rbsx
1752            namstr += '%12s'%('Pos '+px[1])
1753            valstr += '%12.5f'%(parmDict[name])
1754            if name in sigDict:
1755                sigstr += '%12.5f'%(sigDict[name])
1756            else:
1757                sigstr += 12*' '
1758        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1759            name = pfx+rbfx+po+rbsx
1760            namstr += '%12s'%('Ori '+po[1])
1761            valstr += '%12.5f'%(parmDict[name])
1762            if name in sigDict:
1763                sigstr += '%12.5f'%(sigDict[name])
1764            else:
1765                sigstr += 12*' '
1766        pFile.write(namstr+'\n')
1767        pFile.write(valstr+'\n')
1768        pFile.write(sigstr+'\n')
1769       
1770    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1771        namstr = '  names :'
1772        valstr = '  values:'
1773        sigstr = '  esds  :'
1774        if 'N' not in TLS:
1775            pFile.write(' Thermal motion:\n')
1776        if 'T' in TLS:
1777            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1778                name = pfx+rbfx+pt+rbsx
1779                namstr += '%12s'%(pt[:3])
1780                valstr += '%12.5f'%(parmDict[name])
1781                if name in sigDict:
1782                    sigstr += '%12.5f'%(sigDict[name])
1783                else:
1784                    sigstr += 12*' '
1785            pFile.write(namstr+'\n')
1786            pFile.write(valstr+'\n')
1787            pFile.write(sigstr+'\n')
1788        if 'L' in TLS:
1789            namstr = '  names :'
1790            valstr = '  values:'
1791            sigstr = '  esds  :'
1792            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1793                name = pfx+rbfx+pt+rbsx
1794                namstr += '%12s'%(pt[:3])
1795                valstr += '%12.3f'%(parmDict[name])
1796                if name in sigDict:
1797                    sigstr += '%12.3f'%(sigDict[name])
1798                else:
1799                    sigstr += 12*' '
1800            pFile.write(namstr+'\n')
1801            pFile.write(valstr+'\n')
1802            pFile.write(sigstr+'\n')
1803        if 'S' in TLS:
1804            namstr = '  names :'
1805            valstr = '  values:'
1806            sigstr = '  esds  :'
1807            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1808                name = pfx+rbfx+pt+rbsx
1809                namstr += '%12s'%(pt[:3])
1810                valstr += '%12.4f'%(parmDict[name])
1811                if name in sigDict:
1812                    sigstr += '%12.4f'%(sigDict[name])
1813                else:
1814                    sigstr += 12*' '
1815            pFile.write(namstr+'\n')
1816            pFile.write(valstr+'\n')
1817            pFile.write(sigstr+'\n')
1818        if 'U' in TLS:
1819            name = pfx+rbfx+'U:'+rbsx
1820            namstr = '  names :'+'%12s'%('Uiso')
1821            valstr = '  values:'+'%12.5f'%(parmDict[name])
1822            if name in sigDict:
1823                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1824            else:
1825                sigstr = '  esds  :'+12*' '
1826            pFile.write(namstr+'\n')
1827            pFile.write(valstr+'\n')
1828            pFile.write(sigstr+'\n')
1829       
1830    def PrintRBObjTorAndSig(rbsx):
1831        namstr = '  names :'
1832        valstr = '  values:'
1833        sigstr = '  esds  :'
1834        nTors = len(RBObj['Torsions'])
1835        if nTors:
1836            pFile.write(' Torsions:\n')
1837            for it in range(nTors):
1838                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1839                namstr += '%12s'%('Tor'+str(it))
1840                valstr += '%12.4f'%(parmDict[name])
1841                if name in sigDict:
1842                    sigstr += '%12.4f'%(sigDict[name])
1843            pFile.write(namstr+'\n')
1844            pFile.write(valstr+'\n')
1845            pFile.write(sigstr+'\n')
1846               
1847    def PrintSHtextureAndSig(textureData,SHtextureSig):
1848        pFile.write('\n Spherical harmonics texture: Order: %d\n'%textureData['Order'])
1849        names = ['omega','chi','phi']
1850        namstr = '  names :'
1851        ptstr =  '  values:'
1852        sigstr = '  esds  :'
1853        for name in names:
1854            namstr += '%12s'%(name)
1855            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1856            if 'Sample '+name in SHtextureSig:
1857                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1858            else:
1859                sigstr += 12*' '
1860        pFile.write(namstr+'\n')
1861        pFile.write(ptstr+'\n')
1862        pFile.write(sigstr+'\n')
1863        pFile.write('\n Texture coefficients:\n')
1864        SHcoeff = textureData['SH Coeff'][1]
1865        SHkeys = list(SHcoeff.keys())
1866        nCoeff = len(SHcoeff)
1867        nBlock = nCoeff//10+1
1868        iBeg = 0
1869        iFin = min(iBeg+10,nCoeff)
1870        for block in range(nBlock):
1871            namstr = '  names :'
1872            ptstr =  '  values:'
1873            sigstr = '  esds  :'
1874            for name in SHkeys[iBeg:iFin]:
1875                namstr += '%12s'%(name)
1876                ptstr += '%12.3f'%(SHcoeff[name])
1877                if name in SHtextureSig:
1878                    sigstr += '%12.3f'%(SHtextureSig[name])
1879                else:
1880                    sigstr += 12*' '
1881            pFile.write(namstr+'\n')
1882            pFile.write(ptstr+'\n')
1883            pFile.write(sigstr+'\n')
1884            iBeg += 10
1885            iFin = min(iBeg+10,nCoeff)
1886           
1887    pFile.write('\n Phases:\n')
1888    for phase in Phases:
1889        pFile.write(' Result for phase: %s\n'%phase)
1890        pFile.write(135*'='+'\n')
1891        Phase = Phases[phase]
1892        General = Phase['General']
1893        SGData = General['SGData']
1894        Atoms = Phase['Atoms']
1895        AtLookup = []
1896        if Atoms and not General.get('doPawley'):
1897            cx,ct,cs,cia = General['AtomPtrs']
1898            AtLookup = G2mth.FillAtomLookUp(Atoms,cia+8)
1899        cell = General['Cell']
1900        pId = Phase['pId']
1901        pfx = str(pId)+'::'
1902        if cell[0]:
1903            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1904            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1905            pFile.write(' Reciprocal metric tensor: \n')
1906            ptfmt = "%15.9f"
1907            names = ['A11','A22','A33','A12','A13','A23']
1908            namstr = '  names :'
1909            ptstr =  '  values:'
1910            sigstr = '  esds  :'
1911            for name,a,siga in zip(names,A,sigA):
1912                namstr += '%15s'%(name)
1913                ptstr += ptfmt%(a)
1914                if siga:
1915                    sigstr += ptfmt%(siga)
1916                else:
1917                    sigstr += 15*' '
1918            pFile.write(namstr+'\n')
1919            pFile.write(ptstr+'\n')
1920            pFile.write(sigstr+'\n')
1921            cell[1:7] = G2lat.A2cell(A)
1922            cell[7] = G2lat.calc_V(A)
1923            pFile.write(' New unit cell:\n')
1924            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1925            names = ['a','b','c','alpha','beta','gamma','Volume']
1926            namstr = '  names :'
1927            ptstr =  '  values:'
1928            sigstr = '  esds  :'
1929            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1930                namstr += '%12s'%(name)
1931                ptstr += fmt%(a)
1932                if siga:
1933                    sigstr += fmt%(siga)
1934                else:
1935                    sigstr += 12*' '
1936            pFile.write(namstr+'\n')
1937            pFile.write(ptstr+'\n')
1938            pFile.write(sigstr+'\n')
1939        ik = 6  #for Pawley stuff below
1940        if General.get('Modulated',False):
1941            ik = 7
1942            Vec,vRef,maxH = General['SuperVec']
1943            if vRef:
1944                pFile.write(' New modulation vector:\n')
1945                namstr = '  names :'
1946                ptstr =  '  values:'
1947                sigstr = '  esds  :'
1948                for var in ['mV0','mV1','mV2']:
1949                    namstr += '%12s'%(pfx+var)
1950                    ptstr += '%12.6f'%(parmDict[pfx+var])
1951                    if pfx+var in sigDict:
1952                        sigstr += '%12.6f'%(sigDict[pfx+var])
1953                    else:
1954                        sigstr += 12*' '
1955                pFile.write(namstr+'\n')
1956                pFile.write(ptstr+'\n')
1957                pFile.write(sigstr+'\n')
1958           
1959        General['Mass'] = 0.
1960        if Phase['General'].get('doPawley'):
1961            pawleyRef = Phase['Pawley ref']
1962            for i,refl in enumerate(pawleyRef):
1963                key = pfx+'PWLref:'+str(i)
1964                refl[ik] = parmDict[key]
1965                if key in sigDict:
1966                    refl[ik+1] = sigDict[key]
1967                else:
1968                    refl[ik+1] = 0
1969        else:
1970            VRBIds = RBIds['Vector']
1971            RRBIds = RBIds['Residue']
1972            RBModels = Phase['RBModels']
1973            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1974                jrb = VRBIds.index(RBObj['RBId'])
1975                rbsx = str(irb)+':'+str(jrb)
1976                pFile.write(' Vector rigid body parameters:\n')
1977                PrintRBObjPOAndSig('RBV',rbsx)
1978                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1979            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1980                jrb = RRBIds.index(RBObj['RBId'])
1981                rbsx = str(irb)+':'+str(jrb)
1982                pFile.write(' Residue rigid body parameters:\n')
1983                PrintRBObjPOAndSig('RBR',rbsx)
1984                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1985                PrintRBObjTorAndSig(rbsx)
1986            atomsSig = {}
1987            wavesSig = {}
1988            cx,ct,cs,cia = General['AtomPtrs']
1989            for i,at in enumerate(Atoms):
1990                names = {cx:pfx+'Ax:'+str(i),cx+1:pfx+'Ay:'+str(i),cx+2:pfx+'Az:'+str(i),cx+3:pfx+'Afrac:'+str(i),
1991                    cia+1:pfx+'AUiso:'+str(i),cia+2:pfx+'AU11:'+str(i),cia+3:pfx+'AU22:'+str(i),cia+4:pfx+'AU33:'+str(i),
1992                    cia+5:pfx+'AU12:'+str(i),cia+6:pfx+'AU13:'+str(i),cia+7:pfx+'AU23:'+str(i),
1993                    cx+4:pfx+'AMx:'+str(i),cx+5:pfx+'AMy:'+str(i),cx+6:pfx+'AMz:'+str(i)}
1994                for ind in range(cx,cx+4):
1995                    at[ind] = parmDict[names[ind]]
1996                    if ind in range(cx,cx+3):
1997                        name = names[ind].replace('A','dA')
1998                    else:
1999                        name = names[ind]
2000                    if name in sigDict:
2001                        atomsSig[str(i)+':'+str(ind)] = sigDict[name]
2002                if at[cia] == 'I':
2003                    at[cia+1] = parmDict[names[cia+1]]
2004                    if names[cia+1] in sigDict:
2005                        atomsSig['%d:%d'%(i,cia+1)] = sigDict[names[cia+1]]
2006                else:
2007                    for ind in range(cia+2,cia+8):
2008                        at[ind] = parmDict[names[ind]]
2009                        if names[ind] in sigDict:
2010                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2011                if General['Type'] == 'magnetic':
2012                    for ind in range(cx+4,cx+7):
2013                        at[ind] = parmDict[names[ind]]
2014                        if names[ind] in sigDict:
2015                            atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
2016                ind = General['AtomTypes'].index(at[ct])
2017                General['Mass'] += General['AtomMass'][ind]*at[cx+3]*at[cx+5]
2018                if General.get('Modulated',False):
2019                    AtomSS = at[-1]['SS1']
2020                    waveType = AtomSS['waveType']
2021                    for Stype in ['Sfrac','Spos','Sadp','Smag']:
2022                        Waves = AtomSS[Stype]
2023                        for iw,wave in enumerate(Waves):
2024                            stiw = str(i)+':'+str(iw)
2025                            if Stype == 'Spos':
2026                                if waveType in ['ZigZag','Block',] and not iw:
2027                                    names = ['Tmin:'+stiw,'Tmax:'+stiw,'Xmax:'+stiw,'Ymax:'+stiw,'Zmax:'+stiw]
2028                                else:
2029                                    names = ['Xsin:'+stiw,'Ysin:'+stiw,'Zsin:'+stiw,
2030                                        'Xcos:'+stiw,'Ycos:'+stiw,'Zcos:'+stiw]
2031                            elif Stype == 'Sadp':
2032                                names = ['U11sin:'+stiw,'U22sin:'+stiw,'U33sin:'+stiw,
2033                                    'U12sin:'+stiw,'U13sin:'+stiw,'U23sin:'+stiw,
2034                                    'U11cos:'+stiw,'U22cos:'+stiw,'U33cos:'+stiw,
2035                                    'U12cos:'+stiw,'U13cos:'+stiw,'U23cos:'+stiw]
2036                            elif Stype == 'Sfrac':
2037                                if 'Crenel' in waveType and not iw:
2038                                    names = ['Fzero:'+stiw,'Fwid:'+stiw]
2039                                else:
2040                                    names = ['Fsin:'+stiw,'Fcos:'+stiw]
2041                            elif Stype == 'Smag':
2042                                names = ['MXsin:'+stiw,'MYsin:'+stiw,'MZsin:'+stiw,
2043                                    'MXcos:'+stiw,'MYcos:'+stiw,'MZcos:'+stiw]
2044                            for iname,name in enumerate(names):
2045                                AtomSS[Stype][iw][0][iname] = parmDict[pfx+name]
2046                                if pfx+name in sigDict:
2047                                    wavesSig[name] = sigDict[pfx+name]
2048                   
2049            PrintAtomsAndSig(General,Atoms,atomsSig)
2050            if General['Type'] == 'magnetic':
2051                PrintMomentsAndSig(General,Atoms,atomsSig)
2052            if General.get('Modulated',False):
2053                PrintWavesAndSig(General,Atoms,wavesSig)
2054           
2055       
2056        textureData = General['SH Texture']   
2057        if textureData['Order']:
2058            SHtextureSig = {}
2059            for name in ['omega','chi','phi']:
2060                aname = pfx+'SH '+name
2061                textureData['Sample '+name][1] = parmDict[aname]
2062                if aname in sigDict:
2063                    SHtextureSig['Sample '+name] = sigDict[aname]
2064            for name in textureData['SH Coeff'][1]:
2065                aname = pfx+name
2066                textureData['SH Coeff'][1][name] = parmDict[aname]
2067                if aname in sigDict:
2068                    SHtextureSig[name] = sigDict[aname]
2069            PrintSHtextureAndSig(textureData,SHtextureSig)
2070        if phase in RestraintDict and not Phase['General'].get('doPawley'):
2071            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
2072                textureData,RestraintDict[phase],pFile)
2073                   
2074################################################################################
2075##### Histogram & Phase data
2076################################################################################       
2077                   
2078def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
2079    '''Loads the HAP histogram/phase information into dicts
2080
2081    :param dict Phases: phase information
2082    :param dict Histograms: Histogram information
2083    :param bool Print: prints information as it is read
2084    :param file pFile: file object to print to (the default, None causes printing to the console)
2085    :param bool resetRefList: Should the contents of the Reflection List be initialized
2086      on loading. The default, True, initializes the Reflection List as it is loaded.
2087
2088    :returns: (hapVary,hapDict,controlDict)
2089      * hapVary: list of refined variables
2090      * hapDict: dict with refined variables and their values
2091      * controlDict: dict with fixed parameters
2092    '''
2093   
2094    def PrintSize(hapData):
2095        if hapData[0] in ['isotropic','uniaxial']:
2096            line = '\n Size model    : %9s'%(hapData[0])
2097            line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2098            if hapData[0] == 'uniaxial':
2099                line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2100            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2101            pFile.write(line+'\n')
2102        else:
2103            pFile.write('\n Size model    : %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2104                (hapData[0],hapData[1][2],hapData[2][2]))
2105            Snames = ['S11','S22','S33','S12','S13','S23']
2106            ptlbls = ' names :'
2107            ptstr =  ' values:'
2108            varstr = ' refine:'
2109            for i,name in enumerate(Snames):
2110                ptlbls += '%12s' % (name)
2111                ptstr += '%12.3f' % (hapData[4][i])
2112                varstr += '%12s' % (str(hapData[5][i]))
2113            pFile.write(ptlbls+'\n')
2114            pFile.write(ptstr+'\n')
2115            pFile.write(varstr+'\n')
2116       
2117    def PrintMuStrain(hapData,SGData):
2118        if hapData[0] in ['isotropic','uniaxial']:
2119            line = '\n Mustrain model: %9s'%(hapData[0])
2120            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
2121            if hapData[0] == 'uniaxial':
2122                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
2123            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
2124            pFile.write(line+'\n')
2125        else:
2126            pFile.write('\n Mustrain model: %s\n\t LG mixing coeff.:%12.4f Refine? %s\n'%
2127                (hapData[0],hapData[1][2],hapData[2][2]))
2128            Snames = G2spc.MustrainNames(SGData)
2129            ptlbls = ' names :'
2130            ptstr =  ' values:'
2131            varstr = ' refine:'
2132            for i,name in enumerate(Snames):
2133                ptlbls += '%12s' % (name)
2134                ptstr += '%12.1f' % (hapData[4][i])
2135                varstr += '%12s' % (str(hapData[5][i]))
2136            pFile.write(ptlbls+'\n')
2137            pFile.write(ptstr+'\n')
2138            pFile.write(varstr+'\n')
2139
2140    def PrintHStrain(hapData,SGData):
2141        pFile.write('\n Hydrostatic/elastic strain:\n')
2142        Hsnames = G2spc.HStrainNames(SGData)
2143        ptlbls = ' names :'
2144        ptstr =  ' values:'
2145        varstr = ' refine:'
2146        for i,name in enumerate(Hsnames):
2147            ptlbls += '%12s' % (name)
2148            ptstr += '%12.4g' % (hapData[0][i])
2149            varstr += '%12s' % (str(hapData[1][i]))
2150        pFile.write(ptlbls+'\n')
2151        pFile.write(ptstr+'\n')
2152        pFile.write(varstr+'\n')
2153
2154    def PrintSHPO(hapData):
2155        pFile.write('\n Spherical harmonics preferred orientation: Order: %d Refine? %s\n'%(hapData[4],hapData[2]))
2156        ptlbls = ' names :'
2157        ptstr =  ' values:'
2158        for item in hapData[5]:
2159            ptlbls += '%12s'%(item)
2160            ptstr += '%12.3f'%(hapData[5][item]) 
2161        pFile.write(ptlbls+'\n')
2162        pFile.write(ptstr+'\n')
2163   
2164    def PrintBabinet(hapData):
2165        pFile.write('\n Babinet form factor modification:\n')
2166        ptlbls = ' names :'
2167        ptstr =  ' values:'
2168        varstr = ' refine:'
2169        for name in ['BabA','BabU']:
2170            ptlbls += '%12s' % (name)
2171            ptstr += '%12.6f' % (hapData[name][0])
2172            varstr += '%12s' % (str(hapData[name][1]))
2173        pFile.write(ptlbls+'\n')
2174        pFile.write(ptstr+'\n')
2175        pFile.write(varstr+'\n')
2176       
2177    hapDict = {}
2178    hapVary = []
2179    controlDict = {}
2180   
2181    for phase in Phases:
2182        HistoPhase = Phases[phase]['Histograms']
2183        SGData = Phases[phase]['General']['SGData']
2184        cell = Phases[phase]['General']['Cell'][1:7]
2185        A = G2lat.cell2A(cell)
2186        if Phases[phase]['General'].get('Modulated',False):
2187            SSGData = Phases[phase]['General']['SSGData']
2188            Vec,x,maxH = Phases[phase]['General']['SuperVec']
2189        pId = Phases[phase]['pId']
2190        histoList = list(HistoPhase.keys())
2191        histoList.sort()
2192        for histogram in histoList:
2193            try:
2194                Histogram = Histograms[histogram]
2195            except KeyError:                       
2196                #skip if histogram not included e.g. in a sequential refinement
2197                continue
2198            if not HistoPhase[histogram]['Use']:        #remove previously created  & now unused phase reflection list
2199                if phase in Histograms[histogram]['Reflection Lists']:
2200                    del Histograms[histogram]['Reflection Lists'][phase]
2201                continue
2202            hapData = HistoPhase[histogram]
2203            hId = Histogram['hId']
2204            if 'PWDR' in histogram:
2205                limits = Histogram['Limits'][1]
2206                inst = Histogram['Instrument Parameters'][0]    #TODO - grab table here if present
2207                if 'C' in inst['Type'][1]:
2208                    try:
2209                        wave = inst['Lam'][1]
2210                    except KeyError:
2211                        wave = inst['Lam1'][1]
2212                    dmin = wave/(2.0*sind(limits[1]/2.0))
2213                elif 'T' in inst['Type'][0]:
2214                    dmin = limits[0]/inst['difC'][1]
2215                pfx = str(pId)+':'+str(hId)+':'
2216                if Phases[phase]['General']['doPawley']:
2217                    hapDict[pfx+'LeBail'] = False           #Pawley supercedes LeBail
2218                    hapDict[pfx+'newLeBail'] = True
2219                    Tmin = G2lat.Dsp2pos(inst,dmin)
2220                    if 'C' in inst['Type'][1]:
2221                        limits[1] = min(limits[1],Tmin)
2222                    else:
2223                        limits[0] = max(limits[0],Tmin)
2224                else:
2225                    hapDict[pfx+'LeBail'] = hapData.get('LeBail',False)
2226                    hapDict[pfx+'newLeBail'] = hapData.get('newLeBail',True)
2227                if Phases[phase]['General']['Type'] == 'magnetic':
2228                    dmin = max(dmin,Phases[phase]['General']['MagDmin'])
2229                for item in ['Scale','Extinction']:
2230                    hapDict[pfx+item] = hapData[item][0]
2231                    if hapData[item][1] and not hapDict[pfx+'LeBail']:
2232                        hapVary.append(pfx+item)
2233                names = G2spc.HStrainNames(SGData)
2234                HSvals = []
2235                for i,name in enumerate(names):
2236                    hapDict[pfx+name] = hapData['HStrain'][0][i]
2237                    HSvals.append(hapDict[pfx+name])
2238                    if hapData['HStrain'][1][i]:
2239#                    if hapData['HStrain'][1][i] and not hapDict[pfx+'LeBail']:
2240                        hapVary.append(pfx+name)
2241                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
2242                if hapData['Pref.Ori.'][0] == 'MD':
2243                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
2244                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
2245                    if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2246                        hapVary.append(pfx+'MD')
2247                else:                           #'SH' spherical harmonics
2248                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
2249                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
2250                    controlDict[pfx+'SHnames'] = G2lat.GenSHCoeff(SGData['SGLaue'],'0',controlDict[pfx+'SHord'],False)
2251                    controlDict[pfx+'SHhkl'] = []
2252                    try: #patch for old Pref.Ori. items
2253                        controlDict[pfx+'SHtoler'] = 0.1
2254                        if hapData['Pref.Ori.'][6][0] != '':
2255                            controlDict[pfx+'SHhkl'] = [eval(a.replace(' ',',')) for a in hapData['Pref.Ori.'][6]]
2256                        controlDict[pfx+'SHtoler'] = hapData['Pref.Ori.'][7]
2257                    except IndexError:
2258                        pass
2259                    for item in hapData['Pref.Ori.'][5]:
2260                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
2261                        if hapData['Pref.Ori.'][2] and not hapDict[pfx+'LeBail']:
2262                            hapVary.append(pfx+item)
2263                for item in ['Mustrain','Size']:
2264                    controlDict[pfx+item+'Type'] = hapData[item][0]
2265                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
2266                    if hapData[item][2][2]:
2267                        hapVary.append(pfx+item+';mx')
2268                    if hapData[item][0] in ['isotropic','uniaxial']:
2269                        hapDict[pfx+item+';i'] = hapData[item][1][0]
2270                        if hapData[item][2][0]:
2271                            hapVary.append(pfx+item+';i')
2272                        if hapData[item][0] == 'uniaxial':
2273                            controlDict[pfx+item+'Axis'] = hapData[item][3]
2274                            hapDict[pfx+item+';a'] = hapData[item][1][1]
2275                            if hapData[item][2][1]:
2276                                hapVary.append(pfx+item+';a')
2277                    else:       #generalized for mustrain or ellipsoidal for size
2278                        Nterms = len(hapData[item][4])
2279                        if item == 'Mustrain':
2280                            names = G2spc.MustrainNames(SGData)
2281                            pwrs = []
2282                            for name in names:
2283                                h,k,l = name[1:]
2284                                pwrs.append([int(h),int(k),int(l)])
2285                            controlDict[pfx+'MuPwrs'] = pwrs
2286                        for i in range(Nterms):
2287                            sfx = ';'+str(i)
2288                            hapDict[pfx+item+sfx] = hapData[item][4][i]
2289                            if hapData[item][5][i]:
2290                                hapVary.append(pfx+item+sfx)
2291                if Phases[phase]['General']['Type'] != 'magnetic':
2292                    for bab in ['BabA','BabU']:
2293                        hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2294                        if hapData['Babinet'][bab][1] and not hapDict[pfx+'LeBail']:
2295                            hapVary.append(pfx+bab)
2296                               
2297                if Print: 
2298                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2299                    pFile.write(135*'='+'\n')
2300                    if hapDict[pfx+'LeBail']:
2301                        pFile.write(' Perform LeBail extraction\n')                     
2302                    else:
2303                        pFile.write(' Phase fraction  : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2304                        pFile.write(' Extinction coeff: %10.4f Refine? %s\n'%(hapData['Extinction'][0],hapData['Extinction'][1]))
2305                        if hapData['Pref.Ori.'][0] == 'MD':
2306                            Ax = hapData['Pref.Ori.'][3]
2307                            pFile.write(' March-Dollase PO: %10.4f Refine? %s Axis: %d %d %d\n'%
2308                                (hapData['Pref.Ori.'][1],hapData['Pref.Ori.'][2],Ax[0],Ax[1],Ax[2]))
2309                        else: #'SH' for spherical harmonics
2310                            PrintSHPO(hapData['Pref.Ori.'])
2311                            pFile.write(' Penalty hkl list: %s tolerance: %.2f\n'%(controlDict[pfx+'SHhkl'],controlDict[pfx+'SHtoler']))
2312                    PrintSize(hapData['Size'])
2313                    PrintMuStrain(hapData['Mustrain'],SGData)
2314                    PrintHStrain(hapData['HStrain'],SGData)
2315                    if Phases[phase]['General']['Type'] != 'magnetic':
2316                        if hapData['Babinet']['BabA'][0]:
2317                            PrintBabinet(hapData['Babinet'])                       
2318                if resetRefList and (not hapDict[pfx+'LeBail'] or (hapData['LeBail'] and hapData['newLeBail'])):
2319                    if hapData.get('LeBail',True):         #stop regeneating reflections for LeBail
2320                        hapData['newLeBail'] = False
2321                    refList = []
2322                    Uniq = []
2323                    Phi = []
2324                    useExt = 'magnetic' in Phases[phase]['General']['Type'] and 'N' in inst['Type'][0]
2325                    if Phases[phase]['General'].get('Modulated',False):
2326                        ifSuper = True
2327                        HKLd = np.array(G2lat.GenSSHLaue(dmin,SGData,SSGData,Vec,maxH,A))
2328                        HKLd = G2mth.sortArray(HKLd,4,reverse=True)
2329                        for h,k,l,m,d in HKLd:
2330                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2331                            mul *= 2      # for powder overlap of Friedel pairs
2332                            if m or not ext or useExt:
2333                                if 'C' in inst['Type'][0]:
2334                                    pos = G2lat.Dsp2pos(inst,d)
2335                                    if limits[0] < pos < limits[1]:
2336                                        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])
2337                                        #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2338                                        Uniq.append(uniq)
2339                                        Phi.append(phi)
2340                                elif 'T' in inst['Type'][0]:
2341                                    pos = G2lat.Dsp2pos(inst,d)
2342                                    if limits[0] < pos < limits[1]:
2343                                        wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2344                                        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])
2345                                        # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2346                                        #TODO - if tabulated put alp & bet in here
2347                                        Uniq.append(uniq)
2348                                        Phi.append(phi)
2349                    else:
2350                        ifSuper = False
2351                        HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
2352                        HKLd = G2mth.sortArray(HKLd,3,reverse=True)
2353                        for h,k,l,d in HKLd:
2354                            ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
2355                            mul *= 2      # for powder overlap of Friedel pairs
2356                            if ext and not useExt:
2357                                continue
2358                            if 'C' in inst['Type'][0]:
2359                                pos = G2lat.Dsp2pos(inst,d)
2360                                if limits[0] < pos < limits[1]:
2361                                    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])
2362                                    #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
2363                                    Uniq.append(uniq)
2364                                    Phi.append(phi)
2365                            elif 'T' in inst['Type'][0]:
2366                                pos = G2lat.Dsp2pos(inst,d)
2367                                if limits[0] < pos < limits[1]:
2368                                    wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
2369                                    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])
2370                                    # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
2371                                    Uniq.append(uniq)
2372                                    Phi.append(phi)
2373                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{},'Type':inst['Type'][0],'Super':ifSuper}
2374            elif 'HKLF' in histogram:
2375                inst = Histogram['Instrument Parameters'][0]
2376                hId = Histogram['hId']
2377                hfx = ':%d:'%(hId)
2378                for item in inst:
2379                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
2380                    hapDict[hfx+item] = inst[item][1]
2381                pfx = str(pId)+':'+str(hId)+':'
2382                hapDict[pfx+'Scale'] = hapData['Scale'][0]
2383                if hapData['Scale'][1]:
2384                    hapVary.append(pfx+'Scale')
2385                               
2386                extApprox,extType,extParms = hapData['Extinction']
2387                controlDict[pfx+'EType'] = extType
2388                controlDict[pfx+'EApprox'] = extApprox
2389                if 'C' in inst['Type'][0]:
2390                    controlDict[pfx+'Tbar'] = extParms['Tbar']
2391                    controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
2392                if 'Primary' in extType:
2393                    Ekey = ['Ep',]
2394                elif 'I & II' in extType:
2395                    Ekey = ['Eg','Es']
2396                elif 'Secondary Type II' == extType:
2397                    Ekey = ['Es',]
2398                elif 'Secondary Type I' == extType:
2399                    Ekey = ['Eg',]
2400                else:   #'None'
2401                    Ekey = []
2402                for eKey in Ekey:
2403                    hapDict[pfx+eKey] = extParms[eKey][0]
2404                    if extParms[eKey][1]:
2405                        hapVary.append(pfx+eKey)
2406                for bab in ['BabA','BabU']:
2407                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
2408                    if hapData['Babinet'][bab][1]:
2409                        hapVary.append(pfx+bab)
2410                Twins = hapData.get('Twins',[[np.array([[1,0,0],[0,1,0],[0,0,1]]),[1.0,False,0]],])
2411                if len(Twins) == 1:
2412                    hapDict[pfx+'Flack'] = hapData.get('Flack',[0.,False])[0]
2413                    if hapData.get('Flack',[0,False])[1]:
2414                        hapVary.append(pfx+'Flack')
2415                sumTwFr = 0.
2416                controlDict[pfx+'TwinLaw'] = []
2417                controlDict[pfx+'TwinInv'] = []
2418                NTL = 0           
2419                for it,twin in enumerate(Twins):
2420                    if 'bool' in str(type(twin[0])):
2421                        controlDict[pfx+'TwinInv'].append(twin[0])
2422                        controlDict[pfx+'TwinLaw'].append(np.zeros((3,3)))
2423                    else:
2424                        NTL += 1
2425                        controlDict[pfx+'TwinInv'].append(False)
2426                        controlDict[pfx+'TwinLaw'].append(twin[0])
2427                    if it:
2428                        hapDict[pfx+'TwinFr:'+str(it)] = twin[1]
2429                        sumTwFr += twin[1]
2430                    else:
2431                        hapDict[pfx+'TwinFr:0'] = twin[1][0]
2432                        controlDict[pfx+'TwinNMN'] = twin[1][1]
2433                    if Twins[0][1][1]:
2434                        hapVary.append(pfx+'TwinFr:'+str(it))
2435                controlDict[pfx+'NTL'] = NTL
2436                #need to make constraint on TwinFr
2437                controlDict[pfx+'TwinLaw'] = np.array(controlDict[pfx+'TwinLaw'])
2438                if len(Twins) > 1:    #force sum to unity
2439                    hapDict[pfx+'TwinFr:0'] = 1.-sumTwFr
2440                if Print: 
2441                    pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2442                    pFile.write(135*'='+'\n')
2443                    pFile.write(' Scale factor     : %10.4f Refine? %s\n'%(hapData['Scale'][0],hapData['Scale'][1]))
2444                    if extType != 'None':
2445                        pFile.write(' Extinction  Type: %15s approx: %10s\n'%(extType,extApprox))
2446                        text = ' Parameters       :'
2447                        for eKey in Ekey:
2448                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
2449                        pFile.write(text+'\n')
2450                    if hapData['Babinet']['BabA'][0]:
2451                        PrintBabinet(hapData['Babinet'])
2452                    if not SGData['SGInv'] and len(Twins) == 1:
2453                        pFile.write(' Flack parameter: %10.3f Refine? %s\n'%(hapData['Flack'][0],hapData['Flack'][1]))
2454                    if len(Twins) > 1:
2455                        for it,twin in enumerate(Twins):
2456                            if 'bool' in str(type(twin[0])):
2457                                pFile.write(' Nonmerohedral twin fr.: %5.3f Inv? %s Refine?\n'%
2458                                    (hapDict[pfx+'TwinFr:'+str(it)],str(controlDict[pfx+'TwinInv'][it])),Twins[0][1][1]) 
2459                            else:
2460                                pFile.write(' Twin law: %s Twin fr.: %5.3f Refine? %s\n'%
2461                                    (str(twin[0]).replace('\n',','),hapDict[pfx+'TwinFr:'+str(it)],Twins[0][1][1]))
2462                       
2463                Histogram['Reflection Lists'] = phase       
2464               
2465    return hapVary,hapDict,controlDict
2466   
2467def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,FFtables,Print=True,pFile=None):
2468    'needs a doc string'
2469   
2470    def PrintSizeAndSig(hapData,sizeSig):
2471        line = '\n Size model:     %9s'%(hapData[0])
2472        refine = False
2473        if hapData[0] in ['isotropic','uniaxial']:
2474            line += ' equatorial:%12.5f'%(hapData[1][0])
2475            if sizeSig[0][0]:
2476                line += ', sig:%8.4f'%(sizeSig[0][0])
2477                refine = True
2478            if hapData[0] == 'uniaxial':
2479                line += ' axial:%12.4f'%(hapData[1][1])
2480                if sizeSig[0][1]:
2481                    refine = True
2482                    line += ', sig:%8.4f'%(sizeSig[0][1])
2483            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2484            if sizeSig[0][2]:
2485                refine = True
2486                line += ', sig:%8.4f'%(sizeSig[0][2])
2487            if refine:
2488                pFile.write(line+'\n')
2489        else:
2490            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2491            if sizeSig[0][2]:
2492                refine = True
2493                line += ', sig:%8.4f'%(sizeSig[0][2])
2494            Snames = ['S11','S22','S33','S12','S13','S23']
2495            ptlbls = ' name  :'
2496            ptstr =  ' value :'
2497            sigstr = ' sig   :'
2498            for i,name in enumerate(Snames):
2499                ptlbls += '%12s' % (name)
2500                ptstr += '%12.6f' % (hapData[4][i])
2501                if sizeSig[1][i]:
2502                    refine = True
2503                    sigstr += '%12.6f' % (sizeSig[1][i])
2504                else:
2505                    sigstr += 12*' '
2506            if refine:
2507                pFile.write(line+'\n')
2508                pFile.write(ptlbls+'\n')
2509                pFile.write(ptstr+'\n')
2510                pFile.write(sigstr+'\n')
2511       
2512    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
2513        line = '\n Mustrain model: %9s\n'%(hapData[0])
2514        refine = False
2515        if hapData[0] in ['isotropic','uniaxial']:
2516            line += ' equatorial:%12.1f'%(hapData[1][0])
2517            if mustrainSig[0][0]:
2518                line += ', sig:%8.1f'%(mustrainSig[0][0])
2519                refine = True
2520            if hapData[0] == 'uniaxial':
2521                line += ' axial:%12.1f'%(hapData[1][1])
2522                if mustrainSig[0][1]:
2523                     line += ', sig:%8.1f'%(mustrainSig[0][1])
2524            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2525            if mustrainSig[0][2]:
2526                refine = True
2527                line += ', sig:%8.3f'%(mustrainSig[0][2])
2528            if refine:
2529                pFile.write(line+'\n')
2530        else:
2531            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
2532            if mustrainSig[0][2]:
2533                refine = True
2534                line += ', sig:%8.3f'%(mustrainSig[0][2])
2535            Snames = G2spc.MustrainNames(SGData)
2536            ptlbls = ' name  :'
2537            ptstr =  ' value :'
2538            sigstr = ' sig   :'
2539            for i,name in enumerate(Snames):
2540                ptlbls += '%12s' % (name)
2541                ptstr += '%12.1f' % (hapData[4][i])
2542                if mustrainSig[1][i]:
2543                    refine = True
2544                    sigstr += '%12.1f' % (mustrainSig[1][i])
2545                else:
2546                    sigstr += 12*' '
2547            if refine:
2548                pFile.write(line+'\n')
2549                pFile.write(ptlbls+'\n')
2550                pFile.write(ptstr+'\n')
2551                pFile.write(sigstr+'\n')
2552           
2553    def PrintHStrainAndSig(hapData,strainSig,SGData):
2554        Hsnames = G2spc.HStrainNames(SGData)
2555        ptlbls = ' name  :'
2556        ptstr =  ' value :'
2557        sigstr = ' sig   :'
2558        refine = False
2559        for i,name in enumerate(Hsnames):
2560            ptlbls += '%12s' % (name)
2561            ptstr += '%12.4g' % (hapData[0][i])
2562            if name in strainSig:
2563                refine = True
2564                sigstr += '%12.4g' % (strainSig[name])
2565            else:
2566                sigstr += 12*' '
2567        if refine:
2568            pFile.write('\n Hydrostatic/elastic strain:\n')
2569            pFile.write(ptlbls+'\n')
2570            pFile.write(ptstr+'\n')
2571            pFile.write(sigstr+'\n')
2572       
2573    def PrintSHPOAndSig(pfx,hapData,POsig):
2574        pFile.write('\n Spherical harmonics preferred orientation: Order: %d\n'%hapData[4])
2575        ptlbls = ' names :'
2576        ptstr =  ' values:'
2577        sigstr = ' sig   :'
2578        for item in hapData[5]:
2579            ptlbls += '%12s'%(item)
2580            ptstr += '%12.3f'%(hapData[5][item])
2581            if pfx+item in POsig:
2582                sigstr += '%12.3f'%(POsig[pfx+item])
2583            else:
2584                sigstr += 12*' ' 
2585        pFile.write(ptlbls+'\n')
2586        pFile.write(ptstr+'\n')
2587        pFile.write(sigstr+'\n')
2588       
2589    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2590        pFile.write('\n Single crystal extinction: Type: %s Approx: %s\n'%(hapData[0],hapData[1]))
2591        text = ''
2592        for item in hapData[2]:
2593            if pfx+item in ScalExtSig:
2594                text += '       %s: '%(item)
2595                text += '%12.2e'%(hapData[2][item][0])
2596                if pfx+item in ScalExtSig:
2597                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2598        pFile.write(text+'\n')   
2599       
2600    def PrintBabinetAndSig(pfx,hapData,BabSig):
2601        pFile.write('\n Babinet form factor modification:\n')
2602        ptlbls = ' names :'
2603        ptstr =  ' values:'
2604        sigstr = ' sig   :'
2605        for item in hapData:
2606            ptlbls += '%12s'%(item)
2607            ptstr += '%12.3f'%(hapData[item][0])
2608            if pfx+item in BabSig:
2609                sigstr += '%12.3f'%(BabSig[pfx+item])
2610            else:
2611                sigstr += 12*' ' 
2612        pFile.write(ptlbls+'\n')
2613        pFile.write(ptstr+'\n')
2614        pFile.write(sigstr+'\n')
2615       
2616    def PrintTwinsAndSig(pfx,twinData,TwinSig):
2617        pFile.write('\n Twin Law fractions :\n')
2618        ptlbls = ' names :'
2619        ptstr =  ' values:'
2620        sigstr = ' sig   :'
2621        for it,item in enumerate(twinData):
2622            ptlbls += '     twin #%d'%(it)
2623            if it:
2624                ptstr += '%12.3f'%(item[1])
2625            else:
2626                ptstr += '%12.3f'%(item[1][0])
2627            if pfx+'TwinFr:'+str(it) in TwinSig:
2628                sigstr += '%12.3f'%(TwinSig[pfx+'TwinFr:'+str(it)])
2629            else:
2630                sigstr += 12*' ' 
2631        pFile.write(ptlbls+'\n')
2632        pFile.write(ptstr+'\n')
2633        pFile.write(sigstr+'\n')
2634       
2635   
2636    PhFrExtPOSig = {}
2637    SizeMuStrSig = {}
2638    ScalExtSig = {}
2639    BabSig = {}
2640    TwinFrSig = {}
2641    wtFrSum = {}
2642    for phase in Phases:
2643        HistoPhase = Phases[phase]['Histograms']
2644        General = Phases[phase]['General']
2645        SGData = General['SGData']
2646        pId = Phases[phase]['pId']
2647        histoList = list(HistoPhase.keys())
2648        histoList.sort()
2649        for histogram in histoList:
2650            try:
2651                Histogram = Histograms[histogram]
2652            except KeyError:                       
2653                #skip if histogram not included e.g. in a sequential refinement
2654                continue
2655            if not Phases[phase]['Histograms'][histogram]['Use']:
2656                #skip if phase absent from this histogram
2657                continue
2658            hapData = HistoPhase[histogram]
2659            hId = Histogram['hId']
2660            pfx = str(pId)+':'+str(hId)+':'
2661            if hId not in wtFrSum:
2662                wtFrSum[hId] = 0.
2663            if 'PWDR' in histogram:
2664                for item in ['Scale','Extinction']:
2665                    hapData[item][0] = parmDict[pfx+item]
2666                    if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2667                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2668                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2669                if hapData['Pref.Ori.'][0] == 'MD':
2670                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2671                    if pfx+'MD' in sigDict and not parmDict[pfx+'LeBail']:
2672                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2673                else:                           #'SH' spherical harmonics
2674                    for item in hapData['Pref.Ori.'][5]:
2675                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2676                        if pfx+item in sigDict and not parmDict[pfx+'LeBail']:
2677                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2678                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2679                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2680                    pfx+'HStrain':{}})                 
2681                for item in ['Mustrain','Size']:
2682                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2683#                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2684                    if pfx+item+';mx' in sigDict:
2685                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2686                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2687                        hapData[item][1][0] = parmDict[pfx+item+';i']
2688                        if item == 'Size':
2689                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2690                        if pfx+item+';i' in sigDict: 
2691                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2692                        if hapData[item][0] == 'uniaxial':
2693                            hapData[item][1][1] = parmDict[pfx+item+';a']
2694                            if item == 'Size':
2695                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2696                            if pfx+item+';a' in sigDict:
2697                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2698                    else:       #generalized for mustrain or ellipsoidal for size
2699                        Nterms = len(hapData[item][4])
2700                        for i in range(Nterms):
2701                            sfx = ';'+str(i)
2702                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2703                            if pfx+item+sfx in sigDict:
2704                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2705                names = G2spc.HStrainNames(SGData)
2706                for i,name in enumerate(names):
2707                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2708                    if pfx+name in sigDict:
2709                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2710                if Phases[phase]['General']['Type'] != 'magnetic':
2711                    for name in ['BabA','BabU']:
2712                        hapData['Babinet'][name][0] = parmDict[pfx+name]
2713                        if pfx+name in sigDict and not parmDict[pfx+'LeBail']:
2714                            BabSig[pfx+name] = sigDict[pfx+name]               
2715               
2716            elif 'HKLF' in histogram:
2717                for item in ['Scale','Flack']:
2718                    if parmDict.get(pfx+item):
2719                        hapData[item][0] = parmDict[pfx+item]
2720                        if pfx+item in sigDict:
2721                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2722                for item in ['Ep','Eg','Es']:
2723                    if parmDict.get(pfx+item):
2724                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2725                        if pfx+item in sigDict:
2726                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2727                for item in ['BabA','BabU']:
2728                    hapData['Babinet'][item][0] = parmDict[pfx+item]
2729                    if pfx+item in sigDict:
2730                        BabSig[pfx+item] = sigDict[pfx+item]
2731                if 'Twins' in hapData:
2732                    it = 1
2733                    sumTwFr = 0.
2734                    while True:
2735                        try:
2736                            hapData['Twins'][it][1] = parmDict[pfx+'TwinFr:'+str(it)]
2737                            if pfx+'TwinFr:'+str(it) in sigDict:
2738                                TwinFrSig[pfx+'TwinFr:'+str(it)] = sigDict[pfx+'TwinFr:'+str(it)]
2739                            if it:
2740                                sumTwFr += hapData['Twins'][it][1]
2741                            it += 1
2742                        except KeyError:
2743                            break
2744                    hapData['Twins'][0][1][0] = 1.-sumTwFr
2745
2746    if Print:
2747        for phase in Phases:
2748            HistoPhase = Phases[phase]['Histograms']
2749            General = Phases[phase]['General']
2750            SGData = General['SGData']
2751            pId = Phases[phase]['pId']
2752            histoList = list(HistoPhase.keys())
2753            histoList.sort()
2754            for histogram in histoList:
2755                try:
2756                    Histogram = Histograms[histogram]
2757                except KeyError:                       
2758                    #skip if histogram not included e.g. in a sequential refinement
2759                    continue
2760                hapData = HistoPhase[histogram]
2761                hId = Histogram['hId']
2762                Histogram['Residuals'][str(pId)+'::Name'] = phase
2763                pfx = str(pId)+':'+str(hId)+':'
2764                hfx = ':%s:'%(hId)
2765                if pfx+'Nref' not in Histogram['Residuals']:    #skip not used phase in histogram
2766                    continue
2767                pFile.write('\n Phase: %s in histogram: %s\n'%(phase,histogram))
2768                pFile.write(135*'='+'\n')
2769                if 'PWDR' in histogram:
2770                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections\n'%
2771                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref']))
2772                    pFile.write(' Durbin-Watson statistic = %.3f\n'%(Histogram['Residuals']['Durbin-Watson']))
2773                    pFile.write(' Bragg intensity sum = %.3g\n'%(Histogram['Residuals'][pfx+'sumInt']))
2774                   
2775                    if parmDict[pfx+'LeBail']:
2776                        pFile.write(' Performed LeBail extraction for phase %s in histogram %s\n'%(phase,histogram))
2777                    else:
2778                        if pfx+'Scale' in PhFrExtPOSig:
2779                            wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2780                            sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2781                            pFile.write(' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f\n'%
2782                                (hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr))
2783                        if pfx+'Extinction' in PhFrExtPOSig:
2784                            pFile.write(' Extinction coeff: %10.4f, sig %10.4f\n'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction']))
2785                        if hapData['Pref.Ori.'][0] == 'MD':
2786                            if pfx+'MD' in PhFrExtPOSig:
2787                                pFile.write(' March-Dollase PO: %10.4f, sig %10.4f\n'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD']))
2788                        else:
2789                            PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2790                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2791                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2792                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2793                    if Phases[phase]['General']['Type'] != 'magnetic' and not parmDict[pfx+'LeBail']:
2794                        if len(BabSig):
2795                            PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2796                   
2797                elif 'HKLF' in histogram:
2798                    Inst = Histogram['Instrument Parameters'][0]
2799                    pFile.write(' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections (%d user rejected, %d sp.gp.extinct)\n'%
2800                        (Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'],
2801                        Histogram['Residuals'][pfx+'Nrej'],Histogram['Residuals'][pfx+'Next']))
2802                    if FFtables != None and 'N' not in Inst['Type'][0]:
2803                        PrintFprime(FFtables,hfx,pFile)
2804                    pFile.write(' HKLF histogram weight factor = %.3f\n'%(Histogram['wtFactor']))
2805                    if pfx+'Scale' in ScalExtSig:
2806                        pFile.write(' Scale factor : %10.4f, sig %10.4f\n'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale']))
2807                    if hapData['Extinction'][0] != 'None':
2808                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2809                    if len(BabSig):
2810                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2811                    if pfx+'Flack' in ScalExtSig:
2812                        pFile.write(' Flack parameter : %10.3f, sig %10.3f\n'%(hapData['Flack'][0],ScalExtSig[pfx+'Flack']))
2813                    if len(TwinFrSig):
2814                        PrintTwinsAndSig(pfx,hapData['Twins'],TwinFrSig)
2815
2816################################################################################
2817##### Histogram data
2818################################################################################       
2819                   
2820def GetHistogramData(Histograms,Print=True,pFile=None):
2821    'needs a doc string'
2822   
2823    def GetBackgroundParms(hId,Background):
2824        Back = Background[0]
2825        DebyePeaks = Background[1]
2826        bakType,bakFlag = Back[:2]
2827        backVals = Back[3:]
2828        backNames = [':'+str(hId)+':Back;'+str(i) for i in range(len(backVals))]
2829        backDict = dict(zip(backNames,backVals))
2830        backVary = []
2831        if bakFlag:
2832            backVary = backNames
2833        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2834        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2835        debyeDict = {}
2836        debyeList = []
2837        for i in range(DebyePeaks['nDebye']):
2838            debyeNames = [':'+str(hId)+':DebyeA;'+str(i),':'+str(hId)+':DebyeR;'+str(i),':'+str(hId)+':DebyeU;'+str(i)]
2839            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2840            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2841        debyeVary = []
2842        for item in debyeList:
2843            if item[1]:
2844                debyeVary.append(item[0])
2845        backDict.update(debyeDict)
2846        backVary += debyeVary
2847        peakDict = {}
2848        peakList = []
2849        for i in range(DebyePeaks['nPeaks']):
2850            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2851                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2852            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2853            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2854        peakVary = []
2855        for item in peakList:
2856            if item[1]:
2857                peakVary.append(item[0])
2858        backDict.update(peakDict)
2859        backVary += peakVary
2860        return bakType,backDict,backVary           
2861       
2862    def GetInstParms(hId,Inst):
2863        #patch
2864        if 'Z' not in Inst:
2865            Inst['Z'] = [0.0,0.0,False]
2866        dataType = Inst['Type'][0]
2867        instDict = {}
2868        insVary = []
2869        pfx = ':'+str(hId)+':'
2870        insKeys = list(Inst.keys())
2871        insKeys.sort()
2872        for item in insKeys:
2873            insName = pfx+item
2874            instDict[insName] = Inst[item][1]
2875            if len(Inst[item]) > 2 and Inst[item][2]:
2876                insVary.append(insName)
2877        if 'C' in dataType:
2878            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2879        elif 'T' in dataType:   #trap zero alp, bet coeff.
2880            if not instDict[pfx+'alpha']:
2881                instDict[pfx+'alpha'] = 1.0
2882            if not instDict[pfx+'beta-0'] and not instDict[pfx+'beta-1']:
2883                instDict[pfx+'beta-1'] = 1.0
2884        return dataType,instDict,insVary
2885       
2886    def GetSampleParms(hId,Sample):
2887        sampVary = []
2888        hfx = ':'+str(hId)+':'       
2889        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2890            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2891        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2892            if key in Sample:
2893                sampDict[hfx+key] = Sample[key]
2894        Type = Sample['Type']
2895        if 'Bragg' in Type:             #Bragg-Brentano
2896            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2897                sampDict[hfx+item] = Sample[item][0]
2898                if Sample[item][1]:
2899                    sampVary.append(hfx+item)
2900        elif 'Debye' in Type:        #Debye-Scherrer
2901            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2902                sampDict[hfx+item] = Sample[item][0]
2903                if Sample[item][1]:
2904                    sampVary.append(hfx+item)
2905        return Type,sampDict,sampVary
2906       
2907    def PrintBackground(Background):
2908        Back = Background[0]
2909        DebyePeaks = Background[1]
2910        pFile.write('\n Background function: %s Refine? %s\n'%(Back[0],Back[1]))
2911        line = ' Coefficients: '
2912        for i,back in enumerate(Back[3:]):
2913            line += '%10.3f'%(back)
2914            if i and not i%10:
2915                line += '\n'+15*' '
2916        pFile.write(line+'\n')
2917        if DebyePeaks['nDebye']:
2918            pFile.write('\n Debye diffuse scattering coefficients\n')
2919            parms = ['DebyeA','DebyeR','DebyeU']
2920            line = ' names :  '
2921            for parm in parms:
2922                line += '%8s refine?'%(parm)
2923            pFile.write(line+'\n')
2924            for j,term in enumerate(DebyePeaks['debyeTerms']):
2925                line = ' term'+'%2d'%(j)+':'
2926                for i in range(3):
2927                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2928                pFile.write(line+'\n')
2929        if DebyePeaks['nPeaks']:
2930            pFile.write('\n Single peak coefficients\n')
2931            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2932            line = ' names :  '
2933            for parm in parms:
2934                line += '%8s refine?'%(parm)
2935            pFile.write(line+'\n')
2936            for j,term in enumerate(DebyePeaks['peaksList']):
2937                line = ' peak'+'%2d'%(j)+':'
2938                for i in range(4):
2939                    line += '%12.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2940                pFile.write(line+'\n')
2941       
2942    def PrintInstParms(Inst):
2943        pFile.write('\n Instrument Parameters:\n')
2944        insKeys = list(Inst.keys())
2945        insKeys.sort()
2946        iBeg = 0
2947        Ok = True
2948        while Ok:
2949            ptlbls = ' name  :'
2950            ptstr =  ' value :'
2951            varstr = ' refine:'
2952            iFin = min(iBeg+9,len(insKeys))
2953            for item in insKeys[iBeg:iFin]:
2954                if item not in ['Type','Source']:
2955                    ptlbls += '%12s' % (item)
2956                    ptstr += '%12.6f' % (Inst[item][1])
2957                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2958                        varstr += 12*' '
2959                    else:
2960                        varstr += '%12s' % (str(bool(Inst[item][2])))
2961            pFile.write(ptlbls+'\n')
2962            pFile.write(ptstr+'\n')
2963            pFile.write(varstr+'\n')
2964            iBeg = iFin
2965            if iBeg == len(insKeys):
2966                Ok = False
2967            else:
2968                pFile.write('\n')
2969       
2970    def PrintSampleParms(Sample):
2971        pFile.write('\n Sample Parameters:\n')
2972        pFile.write(' Goniometer omega = %.2f, chi = %.2f, phi = %.2f\n'%
2973            (Sample['Omega'],Sample['Chi'],Sample['Phi']))
2974        ptlbls = ' name  :'
2975        ptstr =  ' value :'
2976        varstr = ' refine:'
2977        if 'Bragg' in Sample['Type']:
2978            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2979                ptlbls += '%14s'%(item)
2980                ptstr += '%14.4f'%(Sample[item][0])
2981                varstr += '%14s'%(str(bool(Sample[item][1])))
2982           
2983        elif 'Debye' in Type:        #Debye-Scherrer
2984            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2985                ptlbls += '%14s'%(item)
2986                ptstr += '%14.4f'%(Sample[item][0])
2987                varstr += '%14s'%(str(bool(Sample[item][1])))
2988
2989        pFile.write(ptlbls+'\n')
2990        pFile.write(ptstr+'\n')
2991        pFile.write(varstr+'\n')
2992       
2993    histDict = {}
2994    histVary = []
2995    controlDict = {}
2996    histoList = list(Histograms.keys())
2997    histoList.sort()
2998    for histogram in histoList:
2999        if 'PWDR' in histogram:
3000            Histogram = Histograms[histogram]
3001            hId = Histogram['hId']
3002            pfx = ':'+str(hId)+':'
3003            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3004            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
3005            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
3006            for excl in controlDict[pfx+'Exclude']:
3007                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
3008            if controlDict[pfx+'Exclude']:
3009                ma.mask_rows(Histogram['Data'])
3010            Background = Histogram['Background']
3011            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
3012            controlDict[pfx+'bakType'] = Type
3013            histDict.update(bakDict)
3014            histVary += bakVary
3015           
3016            Inst = Histogram['Instrument Parameters']        #TODO ? ignores tabulated alp,bet & delt for TOF
3017            if 'T' in Type and len(Inst[1]):    #patch -  back-to-back exponential contribution to TOF line shape is removed
3018                print ('Warning: tabulated profile coefficients are ignored')
3019            Type,instDict,insVary = GetInstParms(hId,Inst[0])
3020            controlDict[pfx+'histType'] = Type
3021            if 'XC' in Type:
3022                if pfx+'Lam1' in instDict:
3023                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
3024                else:
3025                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
3026            histDict.update(instDict)
3027            histVary += insVary
3028           
3029            Sample = Histogram['Sample Parameters']
3030            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
3031            controlDict[pfx+'instType'] = Type
3032            histDict.update(sampDict)
3033            histVary += sampVary
3034           
3035   
3036            if Print: 
3037                pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3038                pFile.write(135*'='+'\n')
3039                Units = {'C':' deg','T':' msec'}
3040                units = Units[controlDict[pfx+'histType'][2]]
3041                Limits = controlDict[pfx+'Limits']
3042                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3043                pFile.write(' Histogram limits: %8.2f%s to %8.2f%s\n'%(Limits[0],units,Limits[1],units))
3044                if len(controlDict[pfx+'Exclude']):
3045                    excls = controlDict[pfx+'Exclude']
3046                    for excl in excls:
3047                        pFile.write(' Excluded region:  %8.2f%s to %8.2f%s\n'%(excl[0],units,excl[1],units)) 
3048                PrintSampleParms(Sample)
3049                PrintInstParms(Inst[0])
3050                PrintBackground(Background)
3051        elif 'HKLF' in histogram:
3052            Histogram = Histograms[histogram]
3053            hId = Histogram['hId']
3054            pfx = ':'+str(hId)+':'
3055            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
3056            Inst = Histogram['Instrument Parameters'][0]
3057            controlDict[pfx+'histType'] = Inst['Type'][0]
3058            if 'X' in Inst['Type'][0]:
3059                histDict[pfx+'Lam'] = Inst['Lam'][1]
3060                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
3061    return histVary,histDict,controlDict
3062   
3063def SetHistogramData(parmDict,sigDict,Histograms,FFtables,Print=True,pFile=None):
3064    'needs a doc string'
3065   
3066    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
3067        Back = Background[0]
3068        DebyePeaks = Background[1]
3069        lenBack = len(Back[3:])
3070        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
3071        for i in range(lenBack):
3072            Back[3+i] = parmDict[pfx+'Back;'+str(i)]
3073            if pfx+'Back;'+str(i) in sigDict:
3074                backSig[i] = sigDict[pfx+'Back;'+str(i)]
3075        if DebyePeaks['nDebye']:
3076            for i in range(DebyePeaks['nDebye']):
3077                names = [pfx+'DebyeA;'+str(i),pfx+'DebyeR;'+str(i),pfx+'DebyeU;'+str(i)]
3078                for j,name in enumerate(names):
3079                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
3080                    if name in sigDict:
3081                        backSig[lenBack+3*i+j] = sigDict[name]           
3082        if DebyePeaks['nPeaks']:
3083            for i in range(DebyePeaks['nPeaks']):
3084                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
3085                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
3086                for j,name in enumerate(names):
3087                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
3088                    if name in sigDict:
3089                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
3090        return backSig
3091       
3092    def SetInstParms(pfx,Inst,parmDict,sigDict):
3093        instSig = {}
3094        insKeys = list(Inst.keys())
3095        insKeys.sort()
3096        for item in insKeys:
3097            insName = pfx+item
3098            Inst[item][1] = parmDict[insName]
3099            if insName in sigDict:
3100                instSig[item] = sigDict[insName]
3101            else:
3102                instSig[item] = 0
3103        return instSig
3104       
3105    def SetSampleParms(pfx,Sample,parmDict,sigDict):
3106        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
3107            sampSig = [0 for i in range(5)]
3108            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3109                Sample[item][0] = parmDict[pfx+item]
3110                if pfx+item in sigDict:
3111                    sampSig[i] = sigDict[pfx+item]
3112        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3113            sampSig = [0 for i in range(4)]
3114            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3115                Sample[item][0] = parmDict[pfx+item]
3116                if pfx+item in sigDict:
3117                    sampSig[i] = sigDict[pfx+item]
3118        return sampSig
3119       
3120    def PrintBackgroundSig(Background,backSig):
3121        Back = Background[0]
3122        DebyePeaks = Background[1]
3123        valstr = ' value : '
3124        sigstr = ' sig   : '
3125        refine = False
3126        for i,back in enumerate(Back[3:]):
3127            valstr += '%10.4g'%(back)
3128            if Back[1]:
3129                refine = True
3130                sigstr += '%10.4g'%(backSig[i])
3131            else:
3132                sigstr += 10*' '
3133        if refine:
3134            pFile.write('\n Background function: %s\n'%Back[0])
3135            pFile.write(valstr+'\n')
3136            pFile.write(sigstr+'\n')
3137        if DebyePeaks['nDebye']:
3138            ifAny = False
3139            ptfmt = "%12.3f"
3140            names =  ' names :'
3141            ptstr =  ' values:'
3142            sigstr = ' esds  :'
3143            for item in sigDict:
3144                if 'Debye' in item:
3145                    ifAny = True
3146                    names += '%12s'%(item)
3147                    ptstr += ptfmt%(parmDict[item])
3148                    sigstr += ptfmt%(sigDict[item])
3149            if ifAny:
3150                pFile.write('\n Debye diffuse scattering coefficients\n')
3151                pFile.write(names+'\n')
3152                pFile.write(ptstr+'\n')
3153                pFile.write(sigstr+'\n')
3154        if DebyePeaks['nPeaks']:
3155            pFile.write('\n Single peak coefficients:\n')
3156            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
3157            line = ' peak no. '
3158            for parm in parms:
3159                line += '%14s%12s'%(parm.center(14),'esd'.center(12))
3160            pFile.write(line+'\n')
3161            for ip in range(DebyePeaks['nPeaks']):
3162                ptstr = ' %4d '%(ip)
3163                for parm in parms:
3164                    fmt = '%14.3f'
3165                    efmt = '%12.3f'
3166                    if parm == 'BkPkpos':
3167                        fmt = '%14.4f'
3168                        efmt = '%12.4f'
3169                    name = pfx+parm+';%d'%(ip)
3170                    ptstr += fmt%(parmDict[name])
3171                    if name in sigDict:
3172                        ptstr += efmt%(sigDict[name])
3173                    else:
3174                        ptstr += 12*' '
3175                pFile.write(ptstr+'\n')
3176        sumBk = np.array(Histogram['sumBk'])
3177        pFile.write(' Background sums: empirical %.3g, Debye %.3g, peaks %.3g, Total %.3g\n'%
3178            (sumBk[0],sumBk[1],sumBk[2],np.sum(sumBk)))
3179       
3180    def PrintInstParmsSig(Inst,instSig):
3181        refine = False
3182        insKeys = list(instSig.keys())
3183        insKeys.sort()
3184        iBeg = 0
3185        Ok = True
3186        while Ok:
3187            ptlbls = ' names :'
3188            ptstr =  ' value :'
3189            sigstr = ' sig   :'
3190            iFin = min(iBeg+9,len(insKeys))
3191            for name in insKeys[iBeg:iFin]:
3192                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
3193                    ptlbls += '%12s' % (name)
3194                    ptstr += '%12.6f' % (Inst[name][1])
3195                    if instSig[name]:
3196                        refine = True
3197                        sigstr += '%12.6f' % (instSig[name])
3198                    else:
3199                        sigstr += 12*' '
3200            if refine:
3201                pFile.write('\n Instrument Parameters:\n')
3202                pFile.write(ptlbls+'\n')
3203                pFile.write(ptstr+'\n')
3204                pFile.write(sigstr+'\n')
3205            iBeg = iFin
3206            if iBeg == len(insKeys):
3207                Ok = False
3208       
3209    def PrintSampleParmsSig(Sample,sampleSig):
3210        ptlbls = ' names :'
3211        ptstr =  ' values:'
3212        sigstr = ' sig   :'
3213        refine = False
3214        if 'Bragg' in Sample['Type']:
3215            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
3216                ptlbls += '%14s'%(item)
3217                ptstr += '%14.4f'%(Sample[item][0])
3218                if sampleSig[i]:
3219                    refine = True
3220                    sigstr += '%14.4f'%(sampleSig[i])
3221                else:
3222                    sigstr += 14*' '
3223           
3224        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
3225            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
3226                ptlbls += '%14s'%(item)
3227                ptstr += '%14.4f'%(Sample[item][0])
3228                if sampleSig[i]:
3229                    refine = True
3230                    sigstr += '%14.4f'%(sampleSig[i])
3231                else:
3232                    sigstr += 14*' '
3233
3234        if refine:
3235            pFile.write('\n Sample Parameters:\n')
3236            pFile.write(ptlbls+'\n')
3237            pFile.write(ptstr+'\n')
3238            pFile.write(sigstr+'\n')
3239       
3240    histoList = list(Histograms.keys())
3241    histoList.sort()
3242    for histogram in histoList:
3243        if 'PWDR' in histogram:
3244            Histogram = Histograms[histogram]
3245            hId = Histogram['hId']
3246            pfx = ':'+str(hId)+':'
3247            Background = Histogram['Background']
3248            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
3249           
3250            Inst = Histogram['Instrument Parameters'][0]
3251            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
3252       
3253            Sample = Histogram['Sample Parameters']
3254            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
3255
3256            pFile.write('\n Histogram: %s histogram Id: %d\n'%(histogram,hId))
3257            pFile.write(135*'='+'\n')
3258            pFile.write(' PWDR histogram weight factor = '+'%.3f\n'%(Histogram['wtFactor']))
3259            pFile.write(' Final refinement wR = %.2f%% on %d observations in this histogram\n'%
3260                (Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs']))
3261            pFile.write(' Other residuals: R = %.2f%%, R-bkg = %.2f%%, wR-bkg = %.2f%% wRmin = %.2f%%\n'%
3262                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wR'],Histogram['Residuals']['wRmin']))
3263            if Print:
3264                pFile.write(' Instrument type: %s\n'%Sample['Type'])
3265                if FFtables != None and 'N' not in Inst['Type'][0]:
3266                    PrintFprime(FFtables,pfx,pFile)
3267                PrintSampleParmsSig(Sample,sampSig)
3268                PrintInstParmsSig(Inst,instSig)
3269                PrintBackgroundSig(Background,backSig)
3270               
Note: See TracBrowser for help on using the repository browser.