source: trunk/GSASIIstrIO.py @ 1460

Last change on this file since 1460 was 1460, checked in by vondreele, 7 years ago

fix missing radius Factor for MC/SA plots
add Prfo, Abs & xt to reflection list items
fix most TOF derivatives (remove a multiply by 100!)

  • Property svn:eol-style set to native
  • Property svn:keywords set to Date Author Revision URL Id
File size: 112.7 KB
Line 
1# -*- coding: utf-8 -*-
2########### SVN repository information ###################
3# $Date: 2014-08-12 21:08:52 +0000 (Tue, 12 Aug 2014) $
4# $Author: vondreele $
5# $Revision: 1460 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 1460 2014-08-12 21:08:52Z vondreele $
8########### SVN repository information ###################
9'''
10*GSASIIstrIO: structure I/O routines*
11-------------------------------------
12
13'''
14import sys
15import os
16import os.path as ospath
17import time
18import math
19import copy
20import random
21import cPickle
22import numpy as np
23import numpy.ma as ma
24import numpy.linalg as nl
25import scipy.optimize as so
26import GSASIIpath
27GSASIIpath.SetVersionNumber("$Revision: 1460 $")
28import GSASIIElem as G2el
29import GSASIIgrid as G2gd
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 GetControls(GPXfile):
46    ''' Returns dictionary of control items found in GSASII gpx file
47
48    :param str GPXfile: full .gpx file name
49    :return: dictionary of control items
50    '''
51    Controls = copy.copy(G2obj.DefaultControls)
52    fl = open(GPXfile,'rb')
53    while True:
54        try:
55            data = cPickle.load(fl)
56        except EOFError:
57            break
58        datum = data[0]
59        if datum[0] == 'Controls':
60            Controls.update(datum[1])
61    fl.close()
62    return Controls
63   
64def GetConstraints(GPXfile):
65    '''Read the constraints from the GPX file and interpret them
66
67    called in :func:`CheckConstraints`, :func:`GSASIIstrMain.Refine`
68    and :func:`GSASIIstrMain.SeqRefine`.
69    '''
70    fl = open(GPXfile,'rb')
71    while True:
72        try:
73            data = cPickle.load(fl)
74        except EOFError:
75            break
76        datum = data[0]
77        if datum[0] == 'Constraints':
78            constList = []
79            for item in datum[1]:
80                if item.startswith('_'): continue
81                constList += datum[1][item]
82            fl.close()
83            constDict,fixedList,ignored = ProcessConstraints(constList)
84            if ignored:
85                print ignored,'Constraints were rejected. Was a constrained phase, histogram or atom deleted?'
86            return constDict,fixedList
87    fl.close()
88    raise Exception,"No constraints in GPX file"
89   
90def ProcessConstraints(constList):
91    """Interpret the constraints in the constList input into a dictionary, etc.
92    All :class:`GSASIIobj.G2VarObj` objects are mapped to the appropriate
93    phase/hist/atoms based on the object internals (random Ids). If this can't be
94    done (if a phase has been deleted, etc.), the variable is ignored.
95    If the constraint cannot be used due to too many dropped variables,
96    it is counted as ignored.
97   
98    :param list constList: a list of lists where each item in the outer list
99      specifies a constraint of some form, as described in the :mod:`GSASIIobj`
100      :ref:`Constraint definition <Constraint_definitions_table>`.
101
102    :returns:  a tuple of (constDict,fixedList,ignored) where:
103     
104      * constDict (list of dicts) contains the constraint relationships
105      * fixedList (list) contains the fixed values for each type
106        of constraint.
107      * ignored (int) counts the number of invalid constraint items
108        (should always be zero!)
109    """
110    constDict = []
111    fixedList = []
112    ignored = 0
113    for item in constList:
114        if item[-1] == 'h':
115            # process a hold
116            fixedList.append('0')
117            var = str(item[0][1])
118            if '?' not in var:
119                constDict.append({var:0.0})
120            else:
121                ignored += 1
122        elif item[-1] == 'f':
123            # process a new variable
124            fixedList.append(None)
125            D = {}
126            varyFlag = item[-2]
127            varname = item[-3]
128            for term in item[:-3]:
129                var = str(term[1])
130                if '?' not in var:
131                    D[var] = term[0]
132            if len(D) > 1:
133                # add extra dict terms for input variable name and vary flag
134                if varname is not None:                   
135                    if varname.startswith('::'):
136                        varname = varname[2:].replace(':',';')
137                    else:
138                        varname = varname.replace(':',';')
139                    D['_name'] = '::' + varname
140                D['_vary'] = varyFlag == True # force to bool
141                constDict.append(D)
142            else:
143                ignored += 1
144            #constFlag[-1] = ['Vary']
145        elif item[-1] == 'c': 
146            # process a contraint relationship
147            D = {}
148            for term in item[:-3]:
149                var = str(term[1])
150                if '?' not in var:
151                    D[var] = term[0]
152            if len(D) >= 1:
153                fixedList.append(str(item[-3]))
154                constDict.append(D)
155            else:
156                ignored += 1
157        elif item[-1] == 'e':
158            # process an equivalence
159            firstmult = None
160            eqlist = []
161            for term in item[:-3]:
162                if term[0] == 0: term[0] = 1.0
163                var = str(term[1])
164                if '?' in var: continue
165                if firstmult is None:
166                    firstmult = term[0]
167                    firstvar = var
168                else:
169                    eqlist.append([var,firstmult/term[0]])
170            if len(eqlist) > 0:
171                G2mv.StoreEquivalence(firstvar,eqlist)
172            else:
173                ignored += 1
174        else:
175            ignored += 1
176    return constDict,fixedList,ignored
177
178def CheckConstraints(GPXfile):
179    '''Load constraints and related info and return any error or warning messages'''
180    # init constraints
181    G2mv.InitVars()   
182    # get variables
183    Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile)
184    if not Phases:
185        return 'Error: No Phases!',''
186    if not Histograms:
187        return 'Error: no diffraction data',''
188    rigidbodyDict = GetRigidBodies(GPXfile)
189    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
190    rbVary,rbDict = GetRigidBodyModels(rigidbodyDict,Print=False)
191    Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,RestraintDict=None,rbIds=rbIds,Print=False)
192    hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms,Print=False)
193    histVary,histDict,controlDict = GetHistogramData(Histograms,Print=False)
194    varyList = rbVary+phaseVary+hapVary+histVary
195    constrDict,fixedList = GetConstraints(GPXfile)
196    return G2mv.CheckConstraints(varyList,constrDict,fixedList)
197   
198def GetRestraints(GPXfile):
199    '''Read the restraints from the GPX file.
200    Throws an exception if not found in the .GPX file
201    '''
202    fl = open(GPXfile,'rb')
203    while True:
204        try:
205            data = cPickle.load(fl)
206        except EOFError:
207            break
208        datum = data[0]
209        if datum[0] == 'Restraints':
210            restraintDict = datum[1]
211    fl.close()
212    return restraintDict
213   
214def GetRigidBodies(GPXfile):
215    '''Read the rigid body models from the GPX file
216    '''
217    fl = open(GPXfile,'rb')
218    while True:
219        try:
220            data = cPickle.load(fl)
221        except EOFError:
222            break
223        datum = data[0]
224        if datum[0] == 'Rigid bodies':
225            rigidbodyDict = datum[1]
226    fl.close()
227    return rigidbodyDict
228       
229def GetFprime(controlDict,Histograms):
230    'Needs a doc string'
231    FFtables = controlDict['FFtables']
232    if not FFtables:
233        return
234    histoList = Histograms.keys()
235    histoList.sort()
236    for histogram in histoList:
237        if histogram[:4] in ['PWDR','HKLF']:
238            Histogram = Histograms[histogram]
239            hId = Histogram['hId']
240            hfx = ':%d:'%(hId)
241            if 'X' in controlDict[hfx+'histType']:
242                keV = controlDict[hfx+'keV']
243                for El in FFtables:
244                    Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0])
245                    FP,FPP,Mu = G2el.FPcalc(Orbs, keV)
246                    FFtables[El][hfx+'FP'] = FP
247                    FFtables[El][hfx+'FPP'] = FPP               
248           
249def GetPhaseNames(GPXfile):
250    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
251
252    :param str GPXfile: full .gpx file name
253    :return: list of phase names
254    '''
255    fl = open(GPXfile,'rb')
256    PhaseNames = []
257    while True:
258        try:
259            data = cPickle.load(fl)
260        except EOFError:
261            break
262        datum = data[0]
263        if 'Phases' == datum[0]:
264            for datus in data[1:]:
265                PhaseNames.append(datus[0])
266    fl.close()
267    return PhaseNames
268
269def GetAllPhaseData(GPXfile,PhaseName):
270    ''' Returns the entire dictionary for PhaseName from GSASII gpx file
271
272    :param str GPXfile: full .gpx file name
273    :param str PhaseName: phase name
274    :return: phase dictionary
275    '''       
276    fl = open(GPXfile,'rb')
277    General = {}
278    Atoms = []
279    while True:
280        try:
281            data = cPickle.load(fl)
282        except EOFError:
283            break
284        datum = data[0]
285        if 'Phases' == datum[0]:
286            for datus in data[1:]:
287                if datus[0] == PhaseName:
288                    break
289    fl.close()
290    return datus[1]
291   
292def GetHistograms(GPXfile,hNames):
293    """ Returns a dictionary of histograms found in GSASII gpx file
294
295    :param str GPXfile: full .gpx file name
296    :param str hNames: list of histogram names
297    :return: dictionary of histograms (types = PWDR & HKLF)
298
299    """
300    fl = open(GPXfile,'rb')
301    Histograms = {}
302    while True:
303        try:
304            data = cPickle.load(fl)
305        except EOFError:
306            break
307        datum = data[0]
308        hist = datum[0]
309        if hist in hNames:
310            if 'PWDR' in hist[:4]:
311                PWDRdata = {}
312                PWDRdata.update(datum[1][0])        #weight factor
313                PWDRdata['Data'] = ma.array(ma.getdata(datum[1][1]))          #masked powder data arrays/clear previous masks
314                PWDRdata[data[2][0]] = data[2][1]       #Limits & excluded regions (if any)
315                PWDRdata[data[3][0]] = data[3][1]       #Background
316                PWDRdata[data[4][0]] = data[4][1]       #Instrument parameters
317                PWDRdata[data[5][0]] = data[5][1]       #Sample parameters
318                try:
319                    PWDRdata[data[9][0]] = data[9][1]       #Reflection lists might be missing
320                except IndexError:
321                    PWDRdata['Reflection Lists'] = {}
322                PWDRdata['Residuals'] = {}
323   
324                Histograms[hist] = PWDRdata
325            elif 'HKLF' in hist[:4]:
326                HKLFdata = {}
327                HKLFdata.update(datum[1][0])        #weight factor
328#patch
329                if isinstance(datum[1][1],list):
330                    RefData = {'RefList':[],'FF':{}}
331                    for ref in datum[1][1]:
332                        RefData['RefList'].append(ref[:11]+[ref[13],])
333                    RefData['RefList'] = np.array(RefData['RefList'])
334                    datum[1][1] = RefData
335#end patch
336                datum[1][1]['FF'] = {}
337                HKLFdata['Data'] = datum[1][1]
338                HKLFdata[data[1][0]] = data[1][1]       #Instrument parameters
339                HKLFdata['Reflection Lists'] = None
340                HKLFdata['Residuals'] = {}
341                Histograms[hist] = HKLFdata           
342    fl.close()
343    return Histograms
344   
345def GetHistogramNames(GPXfile,hType):
346    """ Returns a list of histogram names found in GSASII gpx file
347
348    :param str GPXfile: full .gpx file name
349    :param str hType: list of histogram types
350    :return: list of histogram names (types = PWDR & HKLF)
351
352    """
353    fl = open(GPXfile,'rb')
354    HistogramNames = []
355    while True:
356        try:
357            data = cPickle.load(fl)
358        except EOFError:
359            break
360        datum = data[0]
361        if datum[0][:4] in hType:
362            HistogramNames.append(datum[0])
363    fl.close()
364    return HistogramNames
365   
366def GetUsedHistogramsAndPhases(GPXfile):
367    ''' Returns all histograms that are found in any phase
368    and any phase that uses a histogram. This also
369    assigns numbers to used phases and histograms by the
370    order they appear in the file.
371
372    :param str GPXfile: full .gpx file name
373    :returns: (Histograms,Phases)
374
375     * Histograms = dictionary of histograms as {name:data,...}
376     * Phases = dictionary of phases that use histograms
377
378    '''
379    phaseNames = GetPhaseNames(GPXfile)
380    histoList = GetHistogramNames(GPXfile,['PWDR','HKLF'])
381    allHistograms = GetHistograms(GPXfile,histoList)
382    phaseData = {}
383    for name in phaseNames: 
384        phaseData[name] =  GetAllPhaseData(GPXfile,name)
385    Histograms = {}
386    Phases = {}
387    for phase in phaseData:
388        Phase = phaseData[phase]
389        if Phase['Histograms']:
390            if phase not in Phases:
391                pId = phaseNames.index(phase)
392                Phase['pId'] = pId
393                Phases[phase] = Phase
394            for hist in Phase['Histograms']:
395                if 'Use' not in Phase['Histograms'][hist]:      #patch
396                    Phase['Histograms'][hist]['Use'] = True         
397                if hist not in Histograms and Phase['Histograms'][hist]['Use']:
398                    try:
399                        Histograms[hist] = allHistograms[hist]
400                        hId = histoList.index(hist)
401                        Histograms[hist]['hId'] = hId
402                    except KeyError: # would happen if a referenced histogram were
403                        # renamed or deleted
404                        print('For phase "'+str(phase)+
405                              '" unresolved reference to histogram "'+str(hist)+'"')
406    G2obj.IndexAllIds(Histograms=Histograms,Phases=Phases)
407    return Histograms,Phases
408   
409def getBackupName(GPXfile,makeBack):
410    '''
411    Get the name for the backup .gpx file name
412   
413    :param str GPXfile: full .gpx file name
414    :param bool makeBack: if True the name of a new file is returned, if
415      False the name of the last file that exists is returned
416    :returns: the name of a backup file
417   
418    '''
419    GPXpath,GPXname = ospath.split(GPXfile)
420    if GPXpath == '': GPXpath = '.'
421    Name = ospath.splitext(GPXname)[0]
422    files = os.listdir(GPXpath)
423    last = 0
424    for name in files:
425        name = name.split('.')
426        if len(name) == 3 and name[0] == Name and 'bak' in name[1]:
427            if makeBack:
428                last = max(last,int(name[1].strip('bak'))+1)
429            else:
430                last = max(last,int(name[1].strip('bak')))
431    GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx')
432    return GPXback   
433       
434def GPXBackup(GPXfile,makeBack=True):
435    '''
436    makes a backup of the current .gpx file (?)
437   
438    :param str GPXfile: full .gpx file name
439    :param bool makeBack: if True (default), the backup is written to
440      a new file; if False, the last backup is overwritten
441    :returns: the name of the backup file that was written
442    '''
443    import distutils.file_util as dfu
444    GPXback = getBackupName(GPXfile,makeBack)
445    dfu.copy_file(GPXfile,GPXback)
446    return GPXback
447
448def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,RigidBodies,CovData,makeBack=True):
449    ''' Updates gpxfile from all histograms that are found in any phase
450    and any phase that used a histogram. Also updates rigid body definitions.
451
452
453    :param str GPXfile: full .gpx file name
454    :param dict Histograms: dictionary of histograms as {name:data,...}
455    :param dict Phases: dictionary of phases that use histograms
456    :param dict RigidBodies: dictionary of rigid bodies
457    :param dict CovData: dictionary of refined variables, varyList, & covariance matrix
458    :param bool makeBack: True if new backup of .gpx file is to be made; else use the last one made
459
460    '''
461                       
462    GPXback = GPXBackup(GPXfile,makeBack)
463    print 'Read from file:',GPXback
464    print 'Save to file  :',GPXfile
465    infile = open(GPXback,'rb')
466    outfile = open(GPXfile,'wb')
467    while True:
468        try:
469            data = cPickle.load(infile)
470        except EOFError:
471            break
472        datum = data[0]
473#        print 'read: ',datum[0]
474        if datum[0] == 'Phases':
475            for iphase in range(len(data)):
476                if data[iphase][0] in Phases:
477                    phaseName = data[iphase][0]
478                    data[iphase][1].update(Phases[phaseName])
479        elif datum[0] == 'Covariance':
480            data[0][1] = CovData
481        elif datum[0] == 'Rigid bodies':
482            data[0][1] = RigidBodies
483        try:
484            histogram = Histograms[datum[0]]
485#            print 'found ',datum[0]
486            data[0][1][1] = histogram['Data']
487            data[0][1][0].update(histogram['Residuals'])
488            for datus in data[1:]:
489#                print '    read: ',datus[0]
490                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
491                    datus[1] = histogram[datus[0]]
492        except KeyError:
493            pass
494                               
495        cPickle.dump(data,outfile,1)
496    infile.close()
497    outfile.close()
498    print 'GPX file save successful'
499   
500def SetSeqResult(GPXfile,Histograms,SeqResult):
501    '''
502    Needs doc string
503   
504    :param str GPXfile: full .gpx file name
505    '''
506    GPXback = GPXBackup(GPXfile)
507    print 'Read from file:',GPXback
508    print 'Save to file  :',GPXfile
509    infile = open(GPXback,'rb')
510    outfile = open(GPXfile,'wb')
511    while True:
512        try:
513            data = cPickle.load(infile)
514        except EOFError:
515            break
516        datum = data[0]
517        if datum[0] == 'Sequential results':
518            data[0][1] = SeqResult
519        # reset the Copy Next flag, since it should not be needed twice in a row
520        if datum[0] == 'Controls':
521            data[0][1]['Copy2Next'] = False
522        try:
523            histogram = Histograms[datum[0]]
524            data[0][1][1] = list(histogram['Data'])
525            for datus in data[1:]:
526                if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']:
527                    datus[1] = histogram[datus[0]]
528        except KeyError:
529            pass
530                               
531        cPickle.dump(data,outfile,1)
532    infile.close()
533    outfile.close()
534    print 'GPX file save successful'
535                       
536def ShowBanner(pFile=None):
537    'Print authorship, copyright and citation notice'
538    print >>pFile,80*'*'
539    print >>pFile,'   General Structure Analysis System-II Crystal Structure Refinement'
540    print >>pFile,'              by Robert B. Von Dreele & Brian H. Toby'
541    print >>pFile,'                Argonne National Laboratory(C), 2010'
542    print >>pFile,' This product includes software developed by the UChicago Argonne, LLC,' 
543    print >>pFile,'            as Operator of Argonne National Laboratory.'
544    print >>pFile,'                          Please cite:'
545    print >>pFile,'   B.H. Toby & R.B. Von Dreele, J. Appl. Cryst. 46, 544-549 (2013)'
546
547    print >>pFile,80*'*','\n'
548
549def ShowControls(Controls,pFile=None,SeqRef=False):
550    'Print controls information'
551    print >>pFile,' Least squares controls:'
552    print >>pFile,' Refinement type: ',Controls['deriv type']
553    if 'Hessian' in Controls['deriv type']:
554        print >>pFile,' Maximum number of cycles:',Controls['max cyc']
555    else:
556        print >>pFile,' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M'])
557    print >>pFile,' Initial shift factor: ','%.3f'%(Controls['shift factor'])
558    if SeqRef:
559        print >>pFile,' Sequential refinement controls:'
560        print >>pFile,' Copy of histogram results to next: ',Controls['Copy2Next']
561        print >>pFile,' Process histograms in reverse order: ',Controls['Reverse Seq']
562   
563def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary):
564    'needs a doc string'
565#    if SGLaue in ['-1','2/m','mmm']:
566#        return                      #no Pawley symmetry required constraints
567    eqvDict = {}
568    for i,varyI in enumerate(pawleyVary):
569        eqvDict[varyI] = []
570        refI = int(varyI.split(':')[-1])
571        ih,ik,il = PawleyRef[refI][:3]
572        dspI = PawleyRef[refI][4]
573        for varyJ in pawleyVary[i+1:]:
574            refJ = int(varyJ.split(':')[-1])
575            jh,jk,jl = PawleyRef[refJ][:3]
576            dspJ = PawleyRef[refJ][4]
577            if SGLaue in ['4/m','4/mmm']:
578                isum = ih**2+ik**2
579                jsum = jh**2+jk**2
580                if abs(il) == abs(jl) and isum == jsum:
581                    eqvDict[varyI].append(varyJ) 
582            elif SGLaue in ['3R','3mR']:
583                isum = ih**2+ik**2+il**2
584                jsum = jh**2+jk**2+jl**2
585                isum2 = ih*ik+ih*il+ik*il
586                jsum2 = jh*jk+jh*jl+jk*jl
587                if isum == jsum and isum2 == jsum2:
588                    eqvDict[varyI].append(varyJ) 
589            elif SGLaue in ['3','3m1','31m','6/m','6/mmm']:
590                isum = ih**2+ik**2+ih*ik
591                jsum = jh**2+jk**2+jh*jk
592                if abs(il) == abs(jl) and isum == jsum:
593                    eqvDict[varyI].append(varyJ) 
594            elif SGLaue in ['m3','m3m']:
595                isum = ih**2+ik**2+il**2
596                jsum = jh**2+jk**2+jl**2
597                if isum == jsum:
598                    eqvDict[varyI].append(varyJ)
599            elif abs(dspI-dspJ)/dspI < 1.e-4:
600                eqvDict[varyI].append(varyJ) 
601    for item in pawleyVary:
602        if eqvDict[item]:
603            for item2 in pawleyVary:
604                if item2 in eqvDict[item]:
605                    eqvDict[item2] = []
606            G2mv.StoreEquivalence(item,eqvDict[item])
607                   
608def cellVary(pfx,SGData): 
609    'needs a doc string'
610    if SGData['SGLaue'] in ['-1',]:
611        return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
612    elif SGData['SGLaue'] in ['2/m',]:
613        if SGData['SGUniq'] == 'a':
614            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3']
615        elif SGData['SGUniq'] == 'b':
616            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4']
617        else:
618            return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5']
619    elif SGData['SGLaue'] in ['mmm',]:
620        return [pfx+'A0',pfx+'A1',pfx+'A2']
621    elif SGData['SGLaue'] in ['4/m','4/mmm']:
622        return [pfx+'A0',pfx+'A2']
623    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
624        return [pfx+'A0',pfx+'A2']
625    elif SGData['SGLaue'] in ['3R', '3mR']:
626        return [pfx+'A0',pfx+'A3']                       
627    elif SGData['SGLaue'] in ['m3m','m3']:
628        return [pfx+'A0',]
629       
630################################################################################
631##### Rigid Body Models and not General.get('doPawley')
632################################################################################
633       
634def GetRigidBodyModels(rigidbodyDict,Print=True,pFile=None):
635    'needs a doc string'
636   
637    def PrintResRBModel(RBModel):
638        atNames = RBModel['atNames']
639        rbRef = RBModel['rbRef']
640        rbSeq = RBModel['rbSeq']
641        print >>pFile,'Residue RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
642            'No. times used: ',RBModel['useCount']
643        print >>pFile,'    At name       x          y          z'
644        for name,xyz in zip(atNames,RBModel['rbXYZ']):
645            print >>pFile,%8s %10.4f %10.4f %10.4f'%(name,xyz[0],xyz[1],xyz[2])
646        print >>pFile,'Orientation defined by:',atNames[rbRef[0]],' -> ',atNames[rbRef[1]], \
647            ' & ',atNames[rbRef[0]],' -> ',atNames[rbRef[2]]
648        if rbSeq:
649            for i,rbseq in enumerate(rbSeq):
650                print >>pFile,'Torsion sequence ',i,' Bond: '+atNames[rbseq[0]],' - ', \
651                    atNames[rbseq[1]],' riding: ',[atNames[i] for i in rbseq[3]]
652       
653    def PrintVecRBModel(RBModel):
654        rbRef = RBModel['rbRef']
655        atTypes = RBModel['rbTypes']
656        print >>pFile,'Vector RB name: ',RBModel['RBname'],' no.atoms: ',len(RBModel['rbTypes']), \
657            'No. times used: ',RBModel['useCount']
658        for i in range(len(RBModel['VectMag'])):
659            print >>pFile,'Vector no.: ',i,' Magnitude: ', \
660                '%8.4f'%(RBModel['VectMag'][i]),' Refine? ',RBModel['VectRef'][i]
661            print >>pFile,'  No. Type     vx         vy         vz'
662            for j,[name,xyz] in enumerate(zip(atTypes,RBModel['rbVect'][i])):
663                print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(j,name,xyz[0],xyz[1],xyz[2])
664        print >>pFile,'  No. Type      x          y          z'
665        for i,[name,xyz] in enumerate(zip(atTypes,RBModel['rbXYZ'])):
666            print >>pFile,%d   %2s %10.4f %10.4f %10.4f'%(i,name,xyz[0],xyz[1],xyz[2])
667        print >>pFile,'Orientation defined by: atom ',rbRef[0],' -> atom ',rbRef[1], \
668            ' & atom ',rbRef[0],' -> atom ',rbRef[2]
669    rbVary = []
670    rbDict = {}
671    rbIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})
672    if len(rbIds['Vector']):
673        for irb,item in enumerate(rbIds['Vector']):
674            if rigidbodyDict['Vector'][item]['useCount']:
675                RBmags = rigidbodyDict['Vector'][item]['VectMag']
676                RBrefs = rigidbodyDict['Vector'][item]['VectRef']
677                for i,[mag,ref] in enumerate(zip(RBmags,RBrefs)):
678                    pid = '::RBV;'+str(i)+':'+str(irb)
679                    rbDict[pid] = mag
680                    if ref:
681                        rbVary.append(pid)
682                if Print:
683                    print >>pFile,'\nVector rigid body model:'
684                    PrintVecRBModel(rigidbodyDict['Vector'][item])
685    if len(rbIds['Residue']):
686        for item in rbIds['Residue']:
687            if rigidbodyDict['Residue'][item]['useCount']:
688                if Print:
689                    print >>pFile,'\nResidue rigid body model:'
690                    PrintResRBModel(rigidbodyDict['Residue'][item])
691    return rbVary,rbDict
692   
693def SetRigidBodyModels(parmDict,sigDict,rigidbodyDict,pFile=None):
694    'needs a doc string'
695   
696    def PrintRBVectandSig(VectRB,VectSig):
697        print >>pFile,'\n Rigid body vector magnitudes for '+VectRB['RBname']+':'
698        namstr = '  names :'
699        valstr = '  values:'
700        sigstr = '  esds  :'
701        for i,[val,sig] in enumerate(zip(VectRB['VectMag'],VectSig)):
702            namstr += '%12s'%('Vect '+str(i))
703            valstr += '%12.4f'%(val)
704            if sig:
705                sigstr += '%12.4f'%(sig)
706            else:
707                sigstr += 12*' '
708        print >>pFile,namstr
709        print >>pFile,valstr
710        print >>pFile,sigstr       
711       
712    RBIds = rigidbodyDict.get('RBIds',{'Vector':[],'Residue':[]})  #these are lists of rbIds
713    if not RBIds['Vector']:
714        return
715    for irb,item in enumerate(RBIds['Vector']):
716        if rigidbodyDict['Vector'][item]['useCount']:
717            VectSig = []
718            RBmags = rigidbodyDict['Vector'][item]['VectMag']
719            for i,mag in enumerate(RBmags):
720                name = '::RBV;'+str(i)+':'+str(irb)
721                mag = parmDict[name]
722                if name in sigDict:
723                    VectSig.append(sigDict[name])
724            PrintRBVectandSig(rigidbodyDict['Vector'][item],VectSig)   
725       
726################################################################################
727##### Phase data
728################################################################################       
729                   
730def GetPhaseData(PhaseData,RestraintDict={},rbIds={},Print=True,pFile=None):
731    'needs a doc string'
732           
733    def PrintFFtable(FFtable):
734        print >>pFile,'\n X-ray scattering factors:'
735        print >>pFile,'   Symbol     fa                                      fb                                      fc'
736        print >>pFile,99*'-'
737        for Ename in FFtable:
738            ffdata = FFtable[Ename]
739            fa = ffdata['fa']
740            fb = ffdata['fb']
741            print >>pFile,' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' %  \
742                (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc'])
743               
744    def PrintBLtable(BLtable):
745        print >>pFile,'\n Neutron scattering factors:'
746        print >>pFile,'   Symbol   isotope       mass       b       resonant terms'
747        print >>pFile,99*'-'
748        for Ename in BLtable:
749            bldata = BLtable[Ename]
750            isotope = bldata[0]
751            mass = bldata[1]['Mass']
752            if 'SL' in bldata[1]:
753                blen = bldata[1]['SL'][0]
754                bres = []
755            else:
756                blen = 0
757                bres = bldata[1]['BW-LS']
758            line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen)
759            for item in bres:
760                line += '%10.5g'%(item)
761            print >>pFile,line
762           
763    def PrintRBObjects(resRBData,vecRBData):
764       
765        def PrintRBThermals():
766            tlstr = ['11','22','33','12','13','23']
767            sstr = ['12','13','21','23','31','32','AA','BB']
768            TLS = RB['ThermalMotion'][1]
769            TLSvar = RB['ThermalMotion'][2]
770            if 'T' in RB['ThermalMotion'][0]:
771                print >>pFile,'TLS data'
772                text = ''
773                for i in range(6):
774                    text += 'T'+tlstr[i]+' %8.4f %s '%(TLS[i],str(TLSvar[i])[0])
775                print >>pFile,text
776                if 'L' in RB['ThermalMotion'][0]: 
777                    text = ''
778                    for i in range(6,12):
779                        text += 'L'+tlstr[i-6]+' %8.2f %s '%(TLS[i],str(TLSvar[i])[0])
780                    print >>pFile,text
781                if 'S' in RB['ThermalMotion'][0]:
782                    text = ''
783                    for i in range(12,20):
784                        text += 'S'+sstr[i-12]+' %8.3f %s '%(TLS[i],str(TLSvar[i])[0])
785                    print >>pFile,text
786            if 'U' in RB['ThermalMotion'][0]:
787                print >>pFile,'Uiso data'
788                text = 'Uiso'+' %10.3f %s'%(TLS[0],str(TLSvar[0])[0])           
789           
790        if len(resRBData):
791            for RB in resRBData:
792                Oxyz = RB['Orig'][0]
793                Qrijk = RB['Orient'][0]
794                Angle = 2.0*acosd(Qrijk[0])
795                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
796                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]
797                print >>pFile,'Orientation angle,vector:',      \
798                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
799                Torsions = RB['Torsions']
800                if len(Torsions):
801                    text = 'Torsions: '
802                    for torsion in Torsions:
803                        text += '%10.4f Refine? %s'%(torsion[0],torsion[1])
804                    print >>pFile,text
805                PrintRBThermals()
806        if len(vecRBData):
807            for RB in vecRBData:
808                Oxyz = RB['Orig'][0]
809                Qrijk = RB['Orient'][0]
810                Angle = 2.0*acosd(Qrijk[0])
811                print >>pFile,'\nRBObject ',RB['RBname'],' at ',      \
812                    '%10.4f %10.4f %10.4f'%(Oxyz[0],Oxyz[1],Oxyz[2]),' Refine?',RB['Orig'][1]           
813                print >>pFile,'Orientation angle,vector:',      \
814                    '%10.3f %10.4f %10.4f %10.4f'%(Angle,Qrijk[1],Qrijk[2],Qrijk[3]),' Refine? ',RB['Orient'][1]
815                PrintRBThermals()
816               
817    def PrintAtoms(General,Atoms):
818        cx,ct,cs,cia = General['AtomPtrs']
819        print >>pFile,'\n Atoms:'
820        line = '   name    type  refine?   x         y         z    '+ \
821            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
822        if General['Type'] == 'magnetic':
823            line += '   Mx     My     Mz'
824        elif General['Type'] == 'macromolecular':
825            line = ' res no residue chain'+line
826        print >>pFile,line
827        if General['Type'] == 'nuclear':
828            print >>pFile,135*'-'
829            for i,at in enumerate(Atoms):
830                line = '%7s'%(at[ct-1])+'%7s'%(at[ct])+'%7s'%(at[ct+1])+'%10.5f'%(at[cx])+'%10.5f'%(at[cx+1])+ \
831                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
832                if at[cia] == 'I':
833                    line += '%8.4f'%(at[cia+1])+48*' '
834                else:
835                    line += 8*' '
836                    for j in range(6):
837                        line += '%8.4f'%(at[cia+1+j])
838                print >>pFile,line
839        elif General['Type'] == 'macromolecular':
840            print >>pFile,135*'-'           
841            for i,at in enumerate(Atoms):
842                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])+ \
843                    '%10.5f'%(at[cx+2])+'%8.3f'%(at[cx+3])+'%7s'%(at[cs])+'%5d'%(at[cs+1])+'%5s'%(at[cia])
844                if at[cia] == 'I':
845                    line += '%8.4f'%(at[cia+1])+48*' '
846                else:
847                    line += 8*' '
848                    for j in range(6):
849                        line += '%8.4f'%(at[cia+1+j])
850                print >>pFile,line
851       
852    def PrintTexture(textureData):
853        topstr = '\n Spherical harmonics texture: Order:' + \
854            str(textureData['Order'])
855        if textureData['Order']:
856            print >>pFile,topstr+' Refine? '+str(textureData['SH Coeff'][0])
857        else:
858            print >>pFile,topstr
859            return
860        names = ['omega','chi','phi']
861        line = '\n'
862        for name in names:
863            line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0])
864        print >>pFile,line
865        print >>pFile,'\n Texture coefficients:'
866        ptlbls = ' names :'
867        ptstr =  ' values:'
868        SHcoeff = textureData['SH Coeff'][1]
869        for item in SHcoeff:
870            ptlbls += '%12s'%(item)
871            ptstr += '%12.4f'%(SHcoeff[item]) 
872        print >>pFile,ptlbls
873        print >>pFile,ptstr
874       
875    def MakeRBParms(rbKey,phaseVary,phaseDict):
876        rbid = str(rbids.index(RB['RBId']))
877        pfxRB = pfx+'RB'+rbKey+'P'
878        pstr = ['x','y','z']
879        ostr = ['a','i','j','k']
880        for i in range(3):
881            name = pfxRB+pstr[i]+':'+str(iRB)+':'+rbid
882            phaseDict[name] = RB['Orig'][0][i]
883            if RB['Orig'][1]:
884                phaseVary += [name,]
885        pfxRB = pfx+'RB'+rbKey+'O'
886        for i in range(4):
887            name = pfxRB+ostr[i]+':'+str(iRB)+':'+rbid
888            phaseDict[name] = RB['Orient'][0][i]
889            if RB['Orient'][1] == 'AV' and i:
890                phaseVary += [name,]
891            elif RB['Orient'][1] == 'A' and not i:
892                phaseVary += [name,]
893           
894    def MakeRBThermals(rbKey,phaseVary,phaseDict):
895        rbid = str(rbids.index(RB['RBId']))
896        tlstr = ['11','22','33','12','13','23']
897        sstr = ['12','13','21','23','31','32','AA','BB']
898        if 'T' in RB['ThermalMotion'][0]:
899            pfxRB = pfx+'RB'+rbKey+'T'
900            for i in range(6):
901                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
902                phaseDict[name] = RB['ThermalMotion'][1][i]
903                if RB['ThermalMotion'][2][i]:
904                    phaseVary += [name,]
905        if 'L' in RB['ThermalMotion'][0]:
906            pfxRB = pfx+'RB'+rbKey+'L'
907            for i in range(6):
908                name = pfxRB+tlstr[i]+':'+str(iRB)+':'+rbid
909                phaseDict[name] = RB['ThermalMotion'][1][i+6]
910                if RB['ThermalMotion'][2][i+6]:
911                    phaseVary += [name,]
912        if 'S' in RB['ThermalMotion'][0]:
913            pfxRB = pfx+'RB'+rbKey+'S'
914            for i in range(8):
915                name = pfxRB+sstr[i]+':'+str(iRB)+':'+rbid
916                phaseDict[name] = RB['ThermalMotion'][1][i+12]
917                if RB['ThermalMotion'][2][i+12]:
918                    phaseVary += [name,]
919        if 'U' in RB['ThermalMotion'][0]:
920            name = pfx+'RB'+rbKey+'U:'+str(iRB)+':'+rbid
921            phaseDict[name] = RB['ThermalMotion'][1][0]
922            if RB['ThermalMotion'][2][0]:
923                phaseVary += [name,]
924               
925    def MakeRBTorsions(rbKey,phaseVary,phaseDict):
926        rbid = str(rbids.index(RB['RBId']))
927        pfxRB = pfx+'RB'+rbKey+'Tr;'
928        for i,torsion in enumerate(RB['Torsions']):
929            name = pfxRB+str(i)+':'+str(iRB)+':'+rbid
930            phaseDict[name] = torsion[0]
931            if torsion[1]:
932                phaseVary += [name,]
933                   
934    if Print:
935        print  >>pFile,'\n Phases:'
936    phaseVary = []
937    phaseDict = {}
938    phaseConstr = {}
939    pawleyLookup = {}
940    FFtables = {}                   #scattering factors - xrays
941    BLtables = {}                   # neutrons
942    Natoms = {}
943    AtMults = {}
944    AtIA = {}
945    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
946    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
947    atomIndx = {}
948    for name in PhaseData:
949        General = PhaseData[name]['General']
950        pId = PhaseData[name]['pId']
951        pfx = str(pId)+'::'
952        FFtable = G2el.GetFFtable(General['AtomTypes'])
953        BLtable = G2el.GetBLtable(General)
954        FFtables.update(FFtable)
955        BLtables.update(BLtable)
956        Atoms = PhaseData[name]['Atoms']
957        AtLookup = G2mth.FillAtomLookUp(Atoms)
958        PawleyRef = PhaseData[name].get('Pawley ref',[])
959        SGData = General['SGData']
960        SGtext = G2spc.SGPrint(SGData)
961        cell = General['Cell']
962        A = G2lat.cell2A(cell[1:7])
963        phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],
964            pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5],pfx+'Vol':G2lat.calc_V(A)})
965        if cell[0]:
966            phaseVary += cellVary(pfx,SGData)
967        resRBData = PhaseData[name]['RBModels'].get('Residue',[])
968        if resRBData:
969            rbids = rbIds['Residue']    #NB: used in the MakeRB routines
970            for iRB,RB in enumerate(resRBData):
971                MakeRBParms('R',phaseVary,phaseDict)
972                MakeRBThermals('R',phaseVary,phaseDict)
973                MakeRBTorsions('R',phaseVary,phaseDict)
974       
975        vecRBData = PhaseData[name]['RBModels'].get('Vector',[])
976        if vecRBData:
977            rbids = rbIds['Vector']    #NB: used in the MakeRB routines
978            for iRB,RB in enumerate(vecRBData):
979                MakeRBParms('V',phaseVary,phaseDict)
980                MakeRBThermals('V',phaseVary,phaseDict)
981                   
982        Natoms[pfx] = 0
983        if Atoms and not General.get('doPawley'):
984            cx,ct,cs,cia = General['AtomPtrs']
985            if General['Type'] in ['nuclear','macromolecular']:
986                Natoms[pfx] = len(Atoms)
987                for i,at in enumerate(Atoms):
988                    atomIndx[at[-1]] = [pfx,i]      #lookup table for restraints
989                    phaseDict.update({pfx+'Atype:'+str(i):at[ct],pfx+'Afrac:'+str(i):at[cx+3],pfx+'Amul:'+str(i):at[cs+1],
990                        pfx+'Ax:'+str(i):at[cx],pfx+'Ay:'+str(i):at[cx+1],pfx+'Az:'+str(i):at[cx+2],
991                        pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0.,         #refined shifts for x,y,z
992                        pfx+'AI/A:'+str(i):at[cia],})
993                    if at[cia] == 'I':
994                        phaseDict[pfx+'AUiso:'+str(i)] = at[cia+1]
995                    else:
996                        phaseDict.update({pfx+'AU11:'+str(i):at[cia+2],pfx+'AU22:'+str(i):at[cia+3],pfx+'AU33:'+str(i):at[cia+4],
997                            pfx+'AU12:'+str(i):at[cia+5],pfx+'AU13:'+str(i):at[cia+6],pfx+'AU23:'+str(i):at[cia+7]})
998                    if 'F' in at[ct+1]:
999                        phaseVary.append(pfx+'Afrac:'+str(i))
1000                    if 'X' in at[ct+1]:
1001                        xId,xCoef = G2spc.GetCSxinel(at[cs])
1002                        names = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)]
1003                        equivs = [[],[],[]]
1004                        for j in range(3):
1005                            if xId[j] > 0:                               
1006                                phaseVary.append(names[j])
1007                                equivs[xId[j]-1].append([names[j],xCoef[j]])
1008                        for equiv in equivs:
1009                            if len(equiv) > 1:
1010                                name = equiv[0][0]
1011                                for eqv in equiv[1:]:
1012                                    G2mv.StoreEquivalence(name,(eqv,))
1013                    if 'U' in at[ct+1]:
1014                        if at[cia] == 'I':
1015                            phaseVary.append(pfx+'AUiso:'+str(i))
1016                        else:
1017                            uId,uCoef = G2spc.GetCSuinel(at[cs])[:2]
1018                            names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i),
1019                                pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)]
1020                            equivs = [[],[],[],[],[],[]]
1021                            for j in range(6):
1022                                if uId[j] > 0:                               
1023                                    phaseVary.append(names[j])
1024                                    equivs[uId[j]-1].append([names[j],uCoef[j]])
1025                            for equiv in equivs:
1026                                if len(equiv) > 1:
1027                                    name = equiv[0][0]
1028                                    for eqv in equiv[1:]:
1029                                        G2mv.StoreEquivalence(name,(eqv,))
1030#            elif General['Type'] == 'magnetic':
1031#            elif General['Type'] == 'macromolecular':
1032            textureData = General['SH Texture']
1033            if textureData['Order']:
1034                phaseDict[pfx+'SHorder'] = textureData['Order']
1035                phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']]
1036                for item in ['omega','chi','phi']:
1037                    phaseDict[pfx+'SH '+item] = textureData['Sample '+item][1]
1038                    if textureData['Sample '+item][0]:
1039                        phaseVary.append(pfx+'SH '+item)
1040                for item in textureData['SH Coeff'][1]:
1041                    phaseDict[pfx+item] = textureData['SH Coeff'][1][item]
1042                    if textureData['SH Coeff'][0]:
1043                        phaseVary.append(pfx+item)
1044               
1045            if Print:
1046                print >>pFile,'\n Phase name: ',General['Name']
1047                print >>pFile,135*'-'
1048                PrintFFtable(FFtable)
1049                PrintBLtable(BLtable)
1050                print >>pFile,''
1051                for line in SGtext: print >>pFile,line
1052                PrintRBObjects(resRBData,vecRBData)
1053                PrintAtoms(General,Atoms)
1054                print >>pFile,'\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
1055                    ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \
1056                    '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0]
1057                PrintTexture(textureData)
1058                if name in RestraintDict:
1059                    PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1060                        textureData,RestraintDict[name],pFile)
1061                   
1062        elif PawleyRef:
1063            pawleyVary = []
1064            for i,refl in enumerate(PawleyRef):
1065                phaseDict[pfx+'PWLref:'+str(i)] = refl[6]
1066                pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i
1067                if refl[5]:
1068                    pawleyVary.append(pfx+'PWLref:'+str(i))
1069            GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary)      #does G2mv.StoreEquivalence
1070            phaseVary += pawleyVary
1071               
1072    return Natoms,atomIndx,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables
1073   
1074def cellFill(pfx,SGData,parmDict,sigDict): 
1075    '''Returns the filled-out reciprocal cell (A) terms and their uncertainties
1076    from the parameter and sig dictionaries.
1077
1078    :param str pfx: parameter prefix ("n::", where n is a phase number)
1079    :param dict SGdata: a symmetry object
1080    :param dict parmDict: a dictionary of parameters
1081    :param dict sigDict:  a dictionary of uncertainties on parameters
1082
1083    :returns: A,sigA where each is a list of six terms with the A terms
1084    '''
1085    if SGData['SGLaue'] in ['-1',]:
1086        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1087            parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']]
1088    elif SGData['SGLaue'] in ['2/m',]:
1089        if SGData['SGUniq'] == 'a':
1090            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1091                parmDict[pfx+'A3'],0,0]
1092        elif SGData['SGUniq'] == 'b':
1093            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1094                0,parmDict[pfx+'A4'],0]
1095        else:
1096            A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],
1097                0,0,parmDict[pfx+'A5']]
1098    elif SGData['SGLaue'] in ['mmm',]:
1099        A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0]
1100    elif SGData['SGLaue'] in ['4/m','4/mmm']:
1101        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0]
1102    elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1103        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],
1104            parmDict[pfx+'A0'],0,0]
1105    elif SGData['SGLaue'] in ['3R', '3mR']:
1106        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],
1107            parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']]
1108    elif SGData['SGLaue'] in ['m3m','m3']:
1109        A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0]
1110
1111    try:
1112        if SGData['SGLaue'] in ['-1',]:
1113            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1114                sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']]
1115        elif SGData['SGLaue'] in ['2/m',]:
1116            if SGData['SGUniq'] == 'a':
1117                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1118                    sigDict[pfx+'A3'],0,0]
1119            elif SGData['SGUniq'] == 'b':
1120                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1121                    0,sigDict[pfx+'A4'],0]
1122            else:
1123                sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],
1124                    0,0,sigDict[pfx+'A5']]
1125        elif SGData['SGLaue'] in ['mmm',]:
1126            sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0]
1127        elif SGData['SGLaue'] in ['4/m','4/mmm']:
1128            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1129        elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']:
1130            sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0]
1131        elif SGData['SGLaue'] in ['3R', '3mR']:
1132            sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0]
1133        elif SGData['SGLaue'] in ['m3m','m3']:
1134            sigA = [sigDict[pfx+'A0'],0,0,0,0,0]
1135    except KeyError:
1136        sigA = [0,0,0,0,0,0]
1137
1138    return A,sigA
1139       
1140def PrintRestraints(cell,SGData,AtPtrs,Atoms,AtLookup,textureData,phaseRest,pFile):
1141    'needs a doc string'
1142    if phaseRest:
1143        Amat = G2lat.cell2AB(cell)[0]
1144        cx,ct,cs = AtPtrs[:3]
1145        names = [['Bond','Bonds'],['Angle','Angles'],['Plane','Planes'],
1146            ['Chiral','Volumes'],['Torsion','Torsions'],['Rama','Ramas'],
1147            ['ChemComp','Sites'],['Texture','HKLs']]
1148        for name,rest in names:
1149            itemRest = phaseRest[name]
1150            if itemRest[rest] and itemRest['Use']:
1151                print >>pFile,'\n  %s %10.3f Use: %s'%(name+' restraint weight factor',itemRest['wtFactor'],str(itemRest['Use']))
1152                if name in ['Bond','Angle','Plane','Chiral']:
1153                    print >>pFile,'     calc       obs      sig   delt/sig  atoms(symOp): '
1154                    for indx,ops,obs,esd in itemRest[rest]:
1155                        try:
1156                            AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1157                            AtName = ''
1158                            for i,Aname in enumerate(AtNames):
1159                                AtName += Aname
1160                                if ops[i] == '1':
1161                                    AtName += '-'
1162                                else:
1163                                    AtName += '+('+ops[i]+')-'
1164                            XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1165                            XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1166                            if name == 'Bond':
1167                                calc = G2mth.getRestDist(XYZ,Amat)
1168                            elif name == 'Angle':
1169                                calc = G2mth.getRestAngle(XYZ,Amat)
1170                            elif name == 'Plane':
1171                                calc = G2mth.getRestPlane(XYZ,Amat)
1172                            elif name == 'Chiral':
1173                                calc = G2mth.getRestChiral(XYZ,Amat)
1174                            print >>pFile,' %9.3f %9.3f %8.3f %8.3f   %s'%(calc,obs,esd,(obs-calc)/esd,AtName[:-1])
1175                        except KeyError:
1176                            del itemRest[rest]
1177                elif name in ['Torsion','Rama']:
1178                    print >>pFile,'  atoms(symOp)  calc  obs  sig  delt/sig  torsions: '
1179                    coeffDict = itemRest['Coeff']
1180                    for indx,ops,cofName,esd in itemRest[rest]:
1181                        AtNames = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1182                        AtName = ''
1183                        for i,Aname in enumerate(AtNames):
1184                            AtName += Aname+'+('+ops[i]+')-'
1185                        XYZ = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cx,3))
1186                        XYZ = G2mth.getSyXYZ(XYZ,ops,SGData)
1187                        if name == 'Torsion':
1188                            tor = G2mth.getRestTorsion(XYZ,Amat)
1189                            restr,calc = G2mth.calcTorsionEnergy(tor,coeffDict[cofName])
1190                            print >>pFile,' %8.3f %8.3f %.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,tor,AtName[:-1])
1191                        else:
1192                            phi,psi = G2mth.getRestRama(XYZ,Amat)
1193                            restr,calc = G2mth.calcRamaEnergy(phi,psi,coeffDict[cofName])                               
1194                            print >>pFile,' %8.3f %8.3f %8.3f %8.3f %8.3f %8.3f %s'%(calc,obs,esd,(obs-calc)/esd,phi,psi,AtName[:-1])
1195                elif name == 'ChemComp':
1196                    print >>pFile,'     atoms   mul*frac  factor     prod'
1197                    for indx,factors,obs,esd in itemRest[rest]:
1198                        try:
1199                            atoms = G2mth.GetAtomItemsById(Atoms,AtLookup,indx,ct-1)
1200                            mul = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs+1))
1201                            frac = np.array(G2mth.GetAtomItemsById(Atoms,AtLookup,indx,cs-1))
1202                            mulfrac = mul*frac
1203                            calcs = mul*frac*factors
1204                            for iatm,[atom,mf,fr,clc] in enumerate(zip(atoms,mulfrac,factors,calcs)):
1205                                print >>pFile,' %10s %8.3f %8.3f %8.3f'%(atom,mf,fr,clc)
1206                            print >>pFile,' Sum:                   calc: %8.3f obs: %8.3f esd: %8.3f'%(np.sum(calcs),obs,esd)
1207                        except KeyError:
1208                            del itemRest[rest]
1209                elif name == 'Texture' and textureData['Order']:
1210                    Start = False
1211                    SHCoef = textureData['SH Coeff'][1]
1212                    shModels = ['cylindrical','none','shear - 2/m','rolling - mmm']
1213                    SamSym = dict(zip(shModels,['0','-1','2/m','mmm']))
1214                    print '    HKL  grid  neg esd  sum neg  num neg use unit?  unit esd  '
1215                    for hkl,grid,esd1,ifesd2,esd2 in itemRest[rest]:
1216                        PH = np.array(hkl)
1217                        phi,beta = G2lat.CrsAng(np.array(hkl),cell,SGData)
1218                        ODFln = G2lat.Flnh(Start,SHCoef,phi,beta,SGData)
1219                        R,P,Z = G2mth.getRestPolefig(ODFln,SamSym[textureData['Model']],grid)
1220                        Z = ma.masked_greater(Z,0.0)
1221                        num = ma.count(Z)
1222                        sum = 0
1223                        if num:
1224                            sum = np.sum(Z)
1225                        print '   %d %d %d  %d %8.3f %8.3f %8d   %s    %8.3f'%(hkl[0],hkl[1],hkl[2],grid,esd1,sum,num,str(ifesd2),esd2)
1226       
1227def getCellEsd(pfx,SGData,A,covData):
1228    'needs a doc string'
1229    dpr = 180./np.pi
1230    rVsq = G2lat.calc_rVsq(A)
1231    G,g = G2lat.A2Gmat(A)       #get recip. & real metric tensors
1232    cell = np.array(G2lat.Gmat2cell(g))   #real cell
1233    cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell
1234    scos = cosd(cellst[3:6])
1235    ssin = sind(cellst[3:6])
1236    scot = scos/ssin
1237    rcos = cosd(cell[3:6])
1238    rsin = sind(cell[3:6])
1239    rcot = rcos/rsin
1240    RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5']
1241    varyList = covData['varyList']
1242    covMatrix = covData['covMatrix']
1243    vcov = G2mth.getVCov(RMnames,varyList,covMatrix)
1244    Ax = np.array(A)
1245    Ax[3:] /= 2.
1246    drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2,
1247        Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]])
1248    srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T))
1249    Vol = 1/np.sqrt(rVsq)
1250    sigVol = Vol**3*np.sqrt(srcvlsq)/2.
1251    R123 = Ax[0]*Ax[1]*Ax[2]
1252    dsasdg = np.zeros((3,6))
1253    dadg = np.zeros((6,6))
1254    for i0 in range(3):         #0  1   2
1255        i1 = (i0+1)%3           #1  2   0
1256        i2 = (i1+1)%3           #2  0   1
1257        i3 = 5-i2               #3  5   4
1258        i4 = 5-i1               #4  3   5
1259        i5 = 5-i0               #5  4   3
1260        dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1]
1261        dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2]
1262        dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2])
1263        denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2)
1264        denom = np.sqrt(denmsq)
1265        dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2)
1266        dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2)
1267        dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2)
1268        dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2)
1269        dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4])
1270        dadg[i5][i5] = -Ax[i0]/denom
1271    for i0 in range(3):
1272        i1 = (i0+1)%3
1273        i2 = (i1+1)%3
1274        i3 = 5-i2
1275        for ij in range(6):
1276            dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1])
1277            if ij == i0:
1278                dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0]
1279            dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr
1280    sigMat = np.inner(dadg,np.inner(vcov,dadg.T))
1281    var = np.diag(sigMat)
1282    CS = np.where(var>0.,np.sqrt(var),0.)
1283    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
1284   
1285def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1286    'needs a doc string'
1287   
1288    def PrintAtomsAndSig(General,Atoms,atomsSig):
1289        print >>pFile,'\n Atoms:'
1290        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1291        if General['Type'] == 'magnetic':
1292            line += '   Mx     My     Mz'
1293        elif General['Type'] == 'macromolecular':
1294            line = ' res no  residue  chain '+line
1295        print >>pFile,line
1296        if General['Type'] == 'nuclear':
1297            print >>pFile,135*'-'
1298            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1299                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1300            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
1301            for atyp in General['AtomTypes']:       #zero composition
1302                General['NoAtoms'][atyp] = 0.
1303            for i,at in enumerate(Atoms):
1304                General['NoAtoms'][at[1]] += at[6]*float(at[8])     #new composition
1305                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1306                valstr = ' values:'
1307                sigstr = ' sig   :'
1308                for ind in [3,4,5,6]:
1309                    sigind = str(i)+':'+str(ind)
1310                    valstr += fmt[ind]%(at[ind])                   
1311                    if sigind in atomsSig:
1312                        sigstr += fmt[ind]%(atomsSig[sigind])
1313                    else:
1314                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1315                if at[9] == 'I':
1316                    valstr += fmt[10]%(at[10])
1317                    if str(i)+':10' in atomsSig:
1318                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1319                    else:
1320                        sigstr += 8*' '
1321                else:
1322                    valstr += 8*' '
1323                    sigstr += 8*' '
1324                    for ind in [11,12,13,14,15,16]:
1325                        sigind = str(i)+':'+str(ind)
1326                        valstr += fmt[ind]%(at[ind])
1327                        if sigind in atomsSig:                       
1328                            sigstr += fmt[ind]%(atomsSig[sigind])
1329                        else:
1330                            sigstr += 8*' '
1331                print >>pFile,name
1332                print >>pFile,valstr
1333                print >>pFile,sigstr
1334               
1335    def PrintRBObjPOAndSig(rbfx,rbsx):
1336        namstr = '  names :'
1337        valstr = '  values:'
1338        sigstr = '  esds  :'
1339        for i,px in enumerate(['Px:','Py:','Pz:']):
1340            name = pfx+rbfx+px+rbsx
1341            namstr += '%12s'%('Pos '+px[1])
1342            valstr += '%12.5f'%(parmDict[name])
1343            if name in sigDict:
1344                sigstr += '%12.5f'%(sigDict[name])
1345            else:
1346                sigstr += 12*' '
1347        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1348            name = pfx+rbfx+po+rbsx
1349            namstr += '%12s'%('Ori '+po[1])
1350            valstr += '%12.5f'%(parmDict[name])
1351            if name in sigDict:
1352                sigstr += '%12.5f'%(sigDict[name])
1353            else:
1354                sigstr += 12*' '
1355        print >>pFile,namstr
1356        print >>pFile,valstr
1357        print >>pFile,sigstr
1358       
1359    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1360        namstr = '  names :'
1361        valstr = '  values:'
1362        sigstr = '  esds  :'
1363        if 'N' not in TLS:
1364            print >>pFile,' Thermal motion:'
1365        if 'T' in TLS:
1366            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1367                name = pfx+rbfx+pt+rbsx
1368                namstr += '%12s'%(pt[:3])
1369                valstr += '%12.5f'%(parmDict[name])
1370                if name in sigDict:
1371                    sigstr += '%12.5f'%(sigDict[name])
1372                else:
1373                    sigstr += 12*' '
1374            print >>pFile,namstr
1375            print >>pFile,valstr
1376            print >>pFile,sigstr
1377        if 'L' in TLS:
1378            namstr = '  names :'
1379            valstr = '  values:'
1380            sigstr = '  esds  :'
1381            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1382                name = pfx+rbfx+pt+rbsx
1383                namstr += '%12s'%(pt[:3])
1384                valstr += '%12.3f'%(parmDict[name])
1385                if name in sigDict:
1386                    sigstr += '%12.3f'%(sigDict[name])
1387                else:
1388                    sigstr += 12*' '
1389            print >>pFile,namstr
1390            print >>pFile,valstr
1391            print >>pFile,sigstr
1392        if 'S' in TLS:
1393            namstr = '  names :'
1394            valstr = '  values:'
1395            sigstr = '  esds  :'
1396            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1397                name = pfx+rbfx+pt+rbsx
1398                namstr += '%12s'%(pt[:3])
1399                valstr += '%12.4f'%(parmDict[name])
1400                if name in sigDict:
1401                    sigstr += '%12.4f'%(sigDict[name])
1402                else:
1403                    sigstr += 12*' '
1404            print >>pFile,namstr
1405            print >>pFile,valstr
1406            print >>pFile,sigstr
1407        if 'U' in TLS:
1408            name = pfx+rbfx+'U:'+rbsx
1409            namstr = '  names :'+'%12s'%('Uiso')
1410            valstr = '  values:'+'%12.5f'%(parmDict[name])
1411            if name in sigDict:
1412                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1413            else:
1414                sigstr = '  esds  :'+12*' '
1415            print >>pFile,namstr
1416            print >>pFile,valstr
1417            print >>pFile,sigstr
1418       
1419    def PrintRBObjTorAndSig(rbsx):
1420        namstr = '  names :'
1421        valstr = '  values:'
1422        sigstr = '  esds  :'
1423        nTors = len(RBObj['Torsions'])
1424        if nTors:
1425            print >>pFile,' Torsions:'
1426            for it in range(nTors):
1427                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1428                namstr += '%12s'%('Tor'+str(it))
1429                valstr += '%12.4f'%(parmDict[name])
1430                if name in sigDict:
1431                    sigstr += '%12.4f'%(sigDict[name])
1432            print >>pFile,namstr
1433            print >>pFile,valstr
1434            print >>pFile,sigstr
1435               
1436    def PrintSHtextureAndSig(textureData,SHtextureSig):
1437        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1438        names = ['omega','chi','phi']
1439        namstr = '  names :'
1440        ptstr =  '  values:'
1441        sigstr = '  esds  :'
1442        for name in names:
1443            namstr += '%12s'%(name)
1444            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1445            if 'Sample '+name in SHtextureSig:
1446                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1447            else:
1448                sigstr += 12*' '
1449        print >>pFile,namstr
1450        print >>pFile,ptstr
1451        print >>pFile,sigstr
1452        print >>pFile,'\n Texture coefficients:'
1453        namstr = '  names :'
1454        ptstr =  '  values:'
1455        sigstr = '  esds  :'
1456        SHcoeff = textureData['SH Coeff'][1]
1457        for name in SHcoeff:
1458            namstr += '%12s'%(name)
1459            ptstr += '%12.3f'%(SHcoeff[name])
1460            if name in SHtextureSig:
1461                sigstr += '%12.3f'%(SHtextureSig[name])
1462            else:
1463                sigstr += 12*' '
1464        print >>pFile,namstr
1465        print >>pFile,ptstr
1466        print >>pFile,sigstr
1467           
1468    print >>pFile,'\n Phases:'
1469    for phase in Phases:
1470        print >>pFile,' Result for phase: ',phase
1471        Phase = Phases[phase]
1472        General = Phase['General']
1473        SGData = General['SGData']
1474        Atoms = Phase['Atoms']
1475        AtLookup = G2mth.FillAtomLookUp(Atoms)
1476        cell = General['Cell']
1477        pId = Phase['pId']
1478        pfx = str(pId)+'::'
1479        if cell[0]:
1480            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1481            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1482            print >>pFile,' Reciprocal metric tensor: '
1483            ptfmt = "%15.9f"
1484            names = ['A11','A22','A33','A12','A13','A23']
1485            namstr = '  names :'
1486            ptstr =  '  values:'
1487            sigstr = '  esds  :'
1488            for name,a,siga in zip(names,A,sigA):
1489                namstr += '%15s'%(name)
1490                ptstr += ptfmt%(a)
1491                if siga:
1492                    sigstr += ptfmt%(siga)
1493                else:
1494                    sigstr += 15*' '
1495            print >>pFile,namstr
1496            print >>pFile,ptstr
1497            print >>pFile,sigstr
1498            cell[1:7] = G2lat.A2cell(A)
1499            cell[7] = G2lat.calc_V(A)
1500            print >>pFile,' New unit cell:'
1501            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1502            names = ['a','b','c','alpha','beta','gamma','Volume']
1503            namstr = '  names :'
1504            ptstr =  '  values:'
1505            sigstr = '  esds  :'
1506            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1507                namstr += '%12s'%(name)
1508                ptstr += fmt%(a)
1509                if siga:
1510                    sigstr += fmt%(siga)
1511                else:
1512                    sigstr += 12*' '
1513            print >>pFile,namstr
1514            print >>pFile,ptstr
1515            print >>pFile,sigstr
1516           
1517        General['Mass'] = 0.
1518        if Phase['General'].get('doPawley'):
1519            pawleyRef = Phase['Pawley ref']
1520            for i,refl in enumerate(pawleyRef):
1521                key = pfx+'PWLref:'+str(i)
1522                refl[6] = parmDict[key]
1523                if key in sigDict:
1524                    refl[7] = sigDict[key]
1525                else:
1526                    refl[7] = 0
1527        else:
1528            VRBIds = RBIds['Vector']
1529            RRBIds = RBIds['Residue']
1530            RBModels = Phase['RBModels']
1531            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1532                jrb = VRBIds.index(RBObj['RBId'])
1533                rbsx = str(irb)+':'+str(jrb)
1534                print >>pFile,' Vector rigid body parameters:'
1535                PrintRBObjPOAndSig('RBV',rbsx)
1536                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1537            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1538                jrb = RRBIds.index(RBObj['RBId'])
1539                rbsx = str(irb)+':'+str(jrb)
1540                print >>pFile,' Residue rigid body parameters:'
1541                PrintRBObjPOAndSig('RBR',rbsx)
1542                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1543                PrintRBObjTorAndSig(rbsx)
1544            atomsSig = {}
1545            if General['Type'] == 'nuclear':        #this needs macromolecular variant!
1546                for i,at in enumerate(Atoms):
1547                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1548                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1549                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1550                    for ind in [3,4,5,6]:
1551                        at[ind] = parmDict[names[ind]]
1552                        if ind in [3,4,5]:
1553                            name = names[ind].replace('A','dA')
1554                        else:
1555                            name = names[ind]
1556                        if name in sigDict:
1557                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1558                    if at[9] == 'I':
1559                        at[10] = parmDict[names[10]]
1560                        if names[10] in sigDict:
1561                            atomsSig[str(i)+':10'] = sigDict[names[10]]
1562                    else:
1563                        for ind in [11,12,13,14,15,16]:
1564                            at[ind] = parmDict[names[ind]]
1565                            if names[ind] in sigDict:
1566                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1567                    ind = General['AtomTypes'].index(at[1])
1568                    General['Mass'] += General['AtomMass'][ind]*at[6]*at[8]
1569            PrintAtomsAndSig(General,Atoms,atomsSig)
1570       
1571        textureData = General['SH Texture']   
1572        if textureData['Order']:
1573            SHtextureSig = {}
1574            for name in ['omega','chi','phi']:
1575                aname = pfx+'SH '+name
1576                textureData['Sample '+name][1] = parmDict[aname]
1577                if aname in sigDict:
1578                    SHtextureSig['Sample '+name] = sigDict[aname]
1579            for name in textureData['SH Coeff'][1]:
1580                aname = pfx+name
1581                textureData['SH Coeff'][1][name] = parmDict[aname]
1582                if aname in sigDict:
1583                    SHtextureSig[name] = sigDict[aname]
1584            PrintSHtextureAndSig(textureData,SHtextureSig)
1585        if phase in RestraintDict:
1586            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1587                textureData,RestraintDict[phase],pFile)
1588                   
1589################################################################################
1590##### Histogram & Phase data
1591################################################################################       
1592                   
1593def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1594    '''Loads the HAP histogram/phase information into dicts
1595
1596    :param dict Phases: phase information
1597    :param dict Histograms: Histogram information
1598    :param bool Print: prints information as it is read
1599    :param file pFile: file object to print to (the default, None causes printing to the console)
1600    :param bool resetRefList: Should the contents of the Reflection List be initialized
1601      on loading. The default, True, initializes the Reflection List as it is loaded.
1602
1603    :returns: (hapVary,hapDict,controlDict)
1604      * hapVary: list of refined variables
1605      * hapDict: dict with refined variables and their values
1606      * controlDict: dict with computation controls (?)
1607    '''
1608   
1609    def PrintSize(hapData):
1610        if hapData[0] in ['isotropic','uniaxial']:
1611            line = '\n Size model    : %9s'%(hapData[0])
1612            line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1613            if hapData[0] == 'uniaxial':
1614                line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1615            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1616            print >>pFile,line
1617        else:
1618            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
1619                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1620            Snames = ['S11','S22','S33','S12','S13','S23']
1621            ptlbls = ' names :'
1622            ptstr =  ' values:'
1623            varstr = ' refine:'
1624            for i,name in enumerate(Snames):
1625                ptlbls += '%12s' % (name)
1626                ptstr += '%12.6f' % (hapData[4][i])
1627                varstr += '%12s' % (str(hapData[5][i]))
1628            print >>pFile,ptlbls
1629            print >>pFile,ptstr
1630            print >>pFile,varstr
1631       
1632    def PrintMuStrain(hapData,SGData):
1633        if hapData[0] in ['isotropic','uniaxial']:
1634            line = '\n Mustrain model: %9s'%(hapData[0])
1635            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1636            if hapData[0] == 'uniaxial':
1637                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1638            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1639            print >>pFile,line
1640        else:
1641            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1642                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1643            Snames = G2spc.MustrainNames(SGData)
1644            ptlbls = ' names :'
1645            ptstr =  ' values:'
1646            varstr = ' refine:'
1647            for i,name in enumerate(Snames):
1648                ptlbls += '%12s' % (name)
1649                ptstr += '%12.6f' % (hapData[4][i])
1650                varstr += '%12s' % (str(hapData[5][i]))
1651            print >>pFile,ptlbls
1652            print >>pFile,ptstr
1653            print >>pFile,varstr
1654
1655    def PrintHStrain(hapData,SGData):
1656        print >>pFile,'\n Hydrostatic/elastic strain: '
1657        Hsnames = G2spc.HStrainNames(SGData)
1658        ptlbls = ' names :'
1659        ptstr =  ' values:'
1660        varstr = ' refine:'
1661        for i,name in enumerate(Hsnames):
1662            ptlbls += '%12s' % (name)
1663            ptstr += '%12.6f' % (hapData[0][i])
1664            varstr += '%12s' % (str(hapData[1][i]))
1665        print >>pFile,ptlbls
1666        print >>pFile,ptstr
1667        print >>pFile,varstr
1668
1669    def PrintSHPO(hapData):
1670        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1671            str(hapData[4])+' Refine? '+str(hapData[2])
1672        ptlbls = ' names :'
1673        ptstr =  ' values:'
1674        for item in hapData[5]:
1675            ptlbls += '%12s'%(item)
1676            ptstr += '%12.3f'%(hapData[5][item]) 
1677        print >>pFile,ptlbls
1678        print >>pFile,ptstr
1679   
1680    def PrintBabinet(hapData):
1681        print >>pFile,'\n Babinet form factor modification: '
1682        ptlbls = ' names :'
1683        ptstr =  ' values:'
1684        varstr = ' refine:'
1685        for name in ['BabA','BabU']:
1686            ptlbls += '%12s' % (name)
1687            ptstr += '%12.6f' % (hapData[name][0])
1688            varstr += '%12s' % (str(hapData[name][1]))
1689        print >>pFile,ptlbls
1690        print >>pFile,ptstr
1691        print >>pFile,varstr
1692       
1693   
1694    hapDict = {}
1695    hapVary = []
1696    controlDict = {}
1697    poType = {}
1698    poAxes = {}
1699    spAxes = {}
1700    spType = {}
1701   
1702    for phase in Phases:
1703        HistoPhase = Phases[phase]['Histograms']
1704        SGData = Phases[phase]['General']['SGData']
1705        cell = Phases[phase]['General']['Cell'][1:7]
1706        A = G2lat.cell2A(cell)
1707        pId = Phases[phase]['pId']
1708        histoList = HistoPhase.keys()
1709        histoList.sort()
1710        for histogram in histoList:
1711            try:
1712                Histogram = Histograms[histogram]
1713            except KeyError:                       
1714                #skip if histogram not included e.g. in a sequential refinement
1715                continue
1716            hapData = HistoPhase[histogram]
1717            hId = Histogram['hId']
1718            if 'PWDR' in histogram:
1719                limits = Histogram['Limits'][1]
1720                inst = Histogram['Instrument Parameters'][0]
1721                Zero = inst['Zero'][0]
1722                if 'C' in inst['Type'][1]:
1723                    try:
1724                        wave = inst['Lam'][1]
1725                    except KeyError:
1726                        wave = inst['Lam1'][1]
1727                    dmin = wave/(2.0*sind(limits[1]/2.0))
1728                elif 'T' in inst['Type'][0]:
1729                    dmin = limits[0]/inst['difC'][1]
1730                pfx = str(pId)+':'+str(hId)+':'
1731                for item in ['Scale','Extinction']:
1732                    hapDict[pfx+item] = hapData[item][0]
1733                    if hapData[item][1]:
1734                        hapVary.append(pfx+item)
1735                names = G2spc.HStrainNames(SGData)
1736                for i,name in enumerate(names):
1737                    hapDict[pfx+name] = hapData['HStrain'][0][i]
1738                    if hapData['HStrain'][1][i]:
1739                        hapVary.append(pfx+name)
1740                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1741                if hapData['Pref.Ori.'][0] == 'MD':
1742                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1743                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1744                    if hapData['Pref.Ori.'][2]:
1745                        hapVary.append(pfx+'MD')
1746                else:                           #'SH' spherical harmonics
1747                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1748                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1749                    for item in hapData['Pref.Ori.'][5]:
1750                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1751                        if hapData['Pref.Ori.'][2]:
1752                            hapVary.append(pfx+item)
1753                for item in ['Mustrain','Size']:
1754                    controlDict[pfx+item+'Type'] = hapData[item][0]
1755                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
1756                    if hapData[item][2][2]:
1757                        hapVary.append(pfx+item+';mx')
1758                    if hapData[item][0] in ['isotropic','uniaxial']:
1759                        hapDict[pfx+item+';i'] = hapData[item][1][0]
1760                        if hapData[item][2][0]:
1761                            hapVary.append(pfx+item+';i')
1762                        if hapData[item][0] == 'uniaxial':
1763                            controlDict[pfx+item+'Axis'] = hapData[item][3]
1764                            hapDict[pfx+item+';a'] = hapData[item][1][1]
1765                            if hapData[item][2][1]:
1766                                hapVary.append(pfx+item+';a')
1767                    else:       #generalized for mustrain or ellipsoidal for size
1768                        Nterms = len(hapData[item][4])
1769                        if item == 'Mustrain':
1770                            names = G2spc.MustrainNames(SGData)
1771                            pwrs = []
1772                            for name in names:
1773                                h,k,l = name[1:]
1774                                pwrs.append([int(h),int(k),int(l)])
1775                            controlDict[pfx+'MuPwrs'] = pwrs
1776                        for i in range(Nterms):
1777                            sfx = ':'+str(i)
1778                            hapDict[pfx+item+sfx] = hapData[item][4][i]
1779                            if hapData[item][5][i]:
1780                                hapVary.append(pfx+item+sfx)
1781                for bab in ['BabA','BabU']:
1782                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1783                    if hapData['Babinet'][bab][1]:
1784                        hapVary.append(pfx+bab)
1785                               
1786                if Print: 
1787                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1788                    print >>pFile,135*'-'
1789                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1790                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1791                    if hapData['Pref.Ori.'][0] == 'MD':
1792                        Ax = hapData['Pref.Ori.'][3]
1793                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1794                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1795                    else: #'SH' for spherical harmonics
1796                        PrintSHPO(hapData['Pref.Ori.'])
1797                    PrintSize(hapData['Size'])
1798                    PrintMuStrain(hapData['Mustrain'],SGData)
1799                    PrintHStrain(hapData['HStrain'],SGData)
1800                    if hapData['Babinet']['BabA'][0]:
1801                        PrintBabinet(hapData['Babinet'])
1802                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1803                if resetRefList:
1804                    refList = []
1805                    Uniq = []
1806                    Phi = []
1807                    for h,k,l,d in HKLd:
1808                        ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1809                        mul *= 2      # for powder overlap of Friedel pairs
1810                        if ext:
1811                            continue
1812                        if 'C' in inst['Type'][0]:
1813                            pos = 2.0*asind(wave/(2.0*d))+Zero
1814                            if limits[0] < pos < limits[1]:
1815                                refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,1.0,1.0,1.0])
1816                                #... sig,gam,fotsq,fctsq, phase,icorr,prfo,abs,ext
1817                                Uniq.append(uniq)
1818                                Phi.append(phi)
1819                        elif 'T' in inst['Type'][0]:
1820                            pos = inst['difC'][1]*d+inst['difA'][1]*d**2+inst['difB'][1]*d**3+Zero
1821                            if limits[0] < pos < limits[1]:
1822                                wave = inst['difC'][1]*d/(252.816*inst['fltPath'][0])
1823                                refList.append([h,k,l,mul,d, pos,0.0,0.0,0.0,0.0, 0.0,0.0,0.0,0.0,wave, 1.0,1.0,1.0])
1824                                # ... sig,gam,fotsq,fctsq, phase,icorr,alp,bet,wave, prfo,abs,ext
1825                                Uniq.append(uniq)
1826                                Phi.append(phi)
1827                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{}}
1828            elif 'HKLF' in histogram:
1829                inst = Histogram['Instrument Parameters'][0]
1830                hId = Histogram['hId']
1831                hfx = ':%d:'%(hId)
1832                for item in inst:
1833                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
1834                    hapDict[hfx+item] = inst[item][1]
1835                pfx = str(pId)+':'+str(hId)+':'
1836                hapDict[pfx+'Scale'] = hapData['Scale'][0]
1837                if hapData['Scale'][1]:
1838                    hapVary.append(pfx+'Scale')
1839                               
1840                extApprox,extType,extParms = hapData['Extinction']
1841                controlDict[pfx+'EType'] = extType
1842                controlDict[pfx+'EApprox'] = extApprox
1843                controlDict[pfx+'Tbar'] = extParms['Tbar']
1844                controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1845                if 'Primary' in extType:
1846                    Ekey = ['Ep',]
1847                elif 'I & II' in extType:
1848                    Ekey = ['Eg','Es']
1849                elif 'Secondary Type II' == extType:
1850                    Ekey = ['Es',]
1851                elif 'Secondary Type I' == extType:
1852                    Ekey = ['Eg',]
1853                else:   #'None'
1854                    Ekey = []
1855                for eKey in Ekey:
1856                    hapDict[pfx+eKey] = extParms[eKey][0]
1857                    if extParms[eKey][1]:
1858                        hapVary.append(pfx+eKey)
1859                for bab in ['BabA','BabU']:
1860                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1861                    if hapData['Babinet'][bab][1]:
1862                        hapVary.append(pfx+bab)
1863                if Print: 
1864                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1865                    print >>pFile,135*'-'
1866                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1867                    if extType != 'None':
1868                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1869                        text = ' Parameters       :'
1870                        for eKey in Ekey:
1871                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1872                        print >>pFile,text
1873                    if hapData['Babinet']['BabA'][0]:
1874                        PrintBabinet(hapData['Babinet'])
1875                Histogram['Reflection Lists'] = phase       
1876               
1877    return hapVary,hapDict,controlDict
1878   
1879def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1880    'needs a doc string'
1881   
1882    def PrintSizeAndSig(hapData,sizeSig):
1883        line = '\n Size model:     %9s'%(hapData[0])
1884        refine = False
1885        if hapData[0] in ['isotropic','uniaxial']:
1886            line += ' equatorial:%12.5f'%(hapData[1][0])
1887            if sizeSig[0][0]:
1888                line += ', sig:%8.4f'%(sizeSig[0][0])
1889                refine = True
1890            if hapData[0] == 'uniaxial':
1891                line += ' axial:%12.4f'%(hapData[1][1])
1892                if sizeSig[0][1]:
1893                    refine = True
1894                    line += ', sig:%8.4f'%(sizeSig[0][1])
1895            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1896            if sizeSig[0][2]:
1897                refine = True
1898                line += ', sig:%8.4f'%(sizeSig[0][2])
1899            if refine:
1900                print >>pFile,line
1901        else:
1902            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1903            if sizeSig[0][2]:
1904                refine = True
1905                line += ', sig:%8.4f'%(sizeSig[0][2])
1906            Snames = ['S11','S22','S33','S12','S13','S23']
1907            ptlbls = ' name  :'
1908            ptstr =  ' value :'
1909            sigstr = ' sig   :'
1910            for i,name in enumerate(Snames):
1911                ptlbls += '%12s' % (name)
1912                ptstr += '%12.6f' % (hapData[4][i])
1913                if sizeSig[1][i]:
1914                    refine = True
1915                    sigstr += '%12.6f' % (sizeSig[1][i])
1916                else:
1917                    sigstr += 12*' '
1918            if refine:
1919                print >>pFile,line
1920                print >>pFile,ptlbls
1921                print >>pFile,ptstr
1922                print >>pFile,sigstr
1923       
1924    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1925        line = '\n Mustrain model: %9s'%(hapData[0])
1926        refine = False
1927        if hapData[0] in ['isotropic','uniaxial']:
1928            line += ' equatorial:%12.1f'%(hapData[1][0])
1929            if mustrainSig[0][0]:
1930                line += ', sig:%8.1f'%(mustrainSig[0][0])
1931                refine = True
1932            if hapData[0] == 'uniaxial':
1933                line += ' axial:%12.1f'%(hapData[1][1])
1934                if mustrainSig[0][1]:
1935                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1936            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1937            if mustrainSig[0][2]:
1938                refine = True
1939                line += ', sig:%8.3f'%(mustrainSig[0][2])
1940            if refine:
1941                print >>pFile,line
1942        else:
1943            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1944            if mustrainSig[0][2]:
1945                refine = True
1946                line += ', sig:%8.3f'%(mustrainSig[0][2])
1947            Snames = G2spc.MustrainNames(SGData)
1948            ptlbls = ' name  :'
1949            ptstr =  ' value :'
1950            sigstr = ' sig   :'
1951            for i,name in enumerate(Snames):
1952                ptlbls += '%12s' % (name)
1953                ptstr += '%12.6f' % (hapData[4][i])
1954                if mustrainSig[1][i]:
1955                    refine = True
1956                    sigstr += '%12.6f' % (mustrainSig[1][i])
1957                else:
1958                    sigstr += 12*' '
1959            if refine:
1960                print >>pFile,line
1961                print >>pFile,ptlbls
1962                print >>pFile,ptstr
1963                print >>pFile,sigstr
1964           
1965    def PrintHStrainAndSig(hapData,strainSig,SGData):
1966        Hsnames = G2spc.HStrainNames(SGData)
1967        ptlbls = ' name  :'
1968        ptstr =  ' value :'
1969        sigstr = ' sig   :'
1970        refine = False
1971        for i,name in enumerate(Hsnames):
1972            ptlbls += '%12s' % (name)
1973            ptstr += '%12.6g' % (hapData[0][i])
1974            if name in strainSig:
1975                refine = True
1976                sigstr += '%12.6g' % (strainSig[name])
1977            else:
1978                sigstr += 12*' '
1979        if refine:
1980            print >>pFile,'\n Hydrostatic/elastic strain: '
1981            print >>pFile,ptlbls
1982            print >>pFile,ptstr
1983            print >>pFile,sigstr
1984       
1985    def PrintSHPOAndSig(pfx,hapData,POsig):
1986        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1987        ptlbls = ' names :'
1988        ptstr =  ' values:'
1989        sigstr = ' sig   :'
1990        for item in hapData[5]:
1991            ptlbls += '%12s'%(item)
1992            ptstr += '%12.3f'%(hapData[5][item])
1993            if pfx+item in POsig:
1994                sigstr += '%12.3f'%(POsig[pfx+item])
1995            else:
1996                sigstr += 12*' ' 
1997        print >>pFile,ptlbls
1998        print >>pFile,ptstr
1999        print >>pFile,sigstr
2000       
2001    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2002        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2003        text = ''
2004        for item in hapData[2]:
2005            if pfx+item in ScalExtSig:
2006                text += '       %s: '%(item)
2007                text += '%12.2e'%(hapData[2][item][0])
2008                if pfx+item in ScalExtSig:
2009                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2010        print >>pFile,text       
2011       
2012    def PrintBabinetAndSig(pfx,hapData,BabSig):
2013        print >>pFile,'\n Babinet form factor modification: '
2014        ptlbls = ' names :'
2015        ptstr =  ' values:'
2016        sigstr = ' sig   :'
2017        for item in hapData:
2018            ptlbls += '%12s'%(item)
2019            ptstr += '%12.3f'%(hapData[item][0])
2020            if pfx+item in BabSig:
2021                sigstr += '%12.3f'%(BabSig[pfx+item])
2022            else:
2023                sigstr += 12*' ' 
2024        print >>pFile,ptlbls
2025        print >>pFile,ptstr
2026        print >>pFile,sigstr
2027   
2028    PhFrExtPOSig = {}
2029    SizeMuStrSig = {}
2030    ScalExtSig = {}
2031    BabSig = {}
2032    wtFrSum = {}
2033    for phase in Phases:
2034        HistoPhase = Phases[phase]['Histograms']
2035        General = Phases[phase]['General']
2036        SGData = General['SGData']
2037        pId = Phases[phase]['pId']
2038        histoList = HistoPhase.keys()
2039        histoList.sort()
2040        for histogram in histoList:
2041            try:
2042                Histogram = Histograms[histogram]
2043            except KeyError:                       
2044                #skip if histogram not included e.g. in a sequential refinement
2045                continue
2046            hapData = HistoPhase[histogram]
2047            hId = Histogram['hId']
2048            pfx = str(pId)+':'+str(hId)+':'
2049            if hId not in wtFrSum:
2050                wtFrSum[hId] = 0.
2051            if 'PWDR' in histogram:
2052                for item in ['Scale','Extinction']:
2053                    hapData[item][0] = parmDict[pfx+item]
2054                    if pfx+item in sigDict:
2055                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2056                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2057                if hapData['Pref.Ori.'][0] == 'MD':
2058                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2059                    if pfx+'MD' in sigDict:
2060                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2061                else:                           #'SH' spherical harmonics
2062                    for item in hapData['Pref.Ori.'][5]:
2063                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2064                        if pfx+item in sigDict:
2065                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2066                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2067                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2068                    pfx+'HStrain':{}})                 
2069                for item in ['Mustrain','Size']:
2070                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2071                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
2072                    if pfx+item+';mx' in sigDict:
2073                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2074                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2075                        hapData[item][1][0] = parmDict[pfx+item+';i']
2076                        if item == 'Size':
2077                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2078                        if pfx+item+';i' in sigDict: 
2079                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2080                        if hapData[item][0] == 'uniaxial':
2081                            hapData[item][1][1] = parmDict[pfx+item+';a']
2082                            if item == 'Size':
2083                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2084                            if pfx+item+';a' in sigDict:
2085                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2086                    else:       #generalized for mustrain or ellipsoidal for size
2087                        Nterms = len(hapData[item][4])
2088                        for i in range(Nterms):
2089                            sfx = ':'+str(i)
2090                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2091                            if pfx+item+sfx in sigDict:
2092                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2093                names = G2spc.HStrainNames(SGData)
2094                for i,name in enumerate(names):
2095                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2096                    if pfx+name in sigDict:
2097                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2098                for name in ['BabA','BabU']:
2099                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2100                    if pfx+name in sigDict:
2101                        BabSig[pfx+name] = sigDict[pfx+name]               
2102               
2103            elif 'HKLF' in histogram:
2104                for item in ['Scale',]:
2105                    if parmDict.get(pfx+item):
2106                        hapData[item][0] = parmDict[pfx+item]
2107                        if pfx+item in sigDict:
2108                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2109                for item in ['Ep','Eg','Es']:
2110                    if parmDict.get(pfx+item):
2111                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2112                        if pfx+item in sigDict:
2113                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2114                for name in ['BabA','BabU']:
2115                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2116                    if pfx+name in sigDict:
2117                        BabSig[pfx+name] = sigDict[pfx+name]               
2118
2119    if Print:
2120        for phase in Phases:
2121            HistoPhase = Phases[phase]['Histograms']
2122            General = Phases[phase]['General']
2123            SGData = General['SGData']
2124            pId = Phases[phase]['pId']
2125            histoList = HistoPhase.keys()
2126            histoList.sort()
2127            for histogram in histoList:
2128                try:
2129                    Histogram = Histograms[histogram]
2130                except KeyError:                       
2131                    #skip if histogram not included e.g. in a sequential refinement
2132                    continue
2133                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2134                print >>pFile,130*'-'
2135                hapData = HistoPhase[histogram]
2136                hId = Histogram['hId']
2137                Histogram['Residuals'][str(pId)+'::Name'] = phase
2138                pfx = str(pId)+':'+str(hId)+':'
2139                if 'PWDR' in histogram:
2140                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2141                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2142               
2143                    if pfx+'Scale' in PhFrExtPOSig:
2144                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2145                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2146                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2147                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2148                    if pfx+'Extinction' in PhFrExtPOSig:
2149                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2150                    if hapData['Pref.Ori.'][0] == 'MD':
2151                        if pfx+'MD' in PhFrExtPOSig:
2152                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2153                    else:
2154                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2155                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2156                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2157                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2158                    if len(BabSig):
2159                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2160                   
2161                elif 'HKLF' in histogram:
2162                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2163                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2164                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2165                    if pfx+'Scale' in ScalExtSig:
2166                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2167                    if hapData['Extinction'][0] != 'None':
2168                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2169                    if len(BabSig):
2170                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2171
2172################################################################################
2173##### Histogram data
2174################################################################################       
2175                   
2176def GetHistogramData(Histograms,Print=True,pFile=None):
2177    'needs a doc string'
2178   
2179    def GetBackgroundParms(hId,Background):
2180        Back = Background[0]
2181        DebyePeaks = Background[1]
2182        bakType,bakFlag = Back[:2]
2183        backVals = Back[3:]
2184        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2185        backDict = dict(zip(backNames,backVals))
2186        backVary = []
2187        if bakFlag:
2188            backVary = backNames
2189        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2190        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2191        debyeDict = {}
2192        debyeList = []
2193        for i in range(DebyePeaks['nDebye']):
2194            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2195            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2196            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2197        debyeVary = []
2198        for item in debyeList:
2199            if item[1]:
2200                debyeVary.append(item[0])
2201        backDict.update(debyeDict)
2202        backVary += debyeVary
2203        peakDict = {}
2204        peakList = []
2205        for i in range(DebyePeaks['nPeaks']):
2206            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2207                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2208            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2209            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2210        peakVary = []
2211        for item in peakList:
2212            if item[1]:
2213                peakVary.append(item[0])
2214        backDict.update(peakDict)
2215        backVary += peakVary
2216        return bakType,backDict,backVary           
2217       
2218    def GetInstParms(hId,Inst):     
2219        dataType = Inst['Type'][0]
2220        instDict = {}
2221        insVary = []
2222        pfx = ':'+str(hId)+':'
2223        insKeys = Inst.keys()
2224        insKeys.sort()
2225        for item in insKeys:
2226            insName = pfx+item
2227            instDict[insName] = Inst[item][1]
2228            if len(Inst[item]) > 2 and Inst[item][2]:
2229                insVary.append(insName)
2230        if 'C' in dataType:
2231            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2232        return dataType,instDict,insVary
2233       
2234    def GetSampleParms(hId,Sample):
2235        sampVary = []
2236        hfx = ':'+str(hId)+':'       
2237        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2238            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2239        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2240            if key in Sample:
2241                sampDict[hfx+key] = Sample[key]
2242        Type = Sample['Type']
2243        if 'Bragg' in Type:             #Bragg-Brentano
2244            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2245                sampDict[hfx+item] = Sample[item][0]
2246                if Sample[item][1]:
2247                    sampVary.append(hfx+item)
2248        elif 'Debye' in Type:        #Debye-Scherrer
2249            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2250                sampDict[hfx+item] = Sample[item][0]
2251                if Sample[item][1]:
2252                    sampVary.append(hfx+item)
2253        return Type,sampDict,sampVary
2254       
2255    def PrintBackground(Background):
2256        Back = Background[0]
2257        DebyePeaks = Background[1]
2258        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2259        line = ' Coefficients: '
2260        for i,back in enumerate(Back[3:]):
2261            line += '%10.3f'%(back)
2262            if i and not i%10:
2263                line += '\n'+15*' '
2264        print >>pFile,line
2265        if DebyePeaks['nDebye']:
2266            print >>pFile,'\n Debye diffuse scattering coefficients'
2267            parms = ['DebyeA','DebyeR','DebyeU']
2268            line = ' names :  '
2269            for parm in parms:
2270                line += '%8s refine?'%(parm)
2271            print >>pFile,line
2272            for j,term in enumerate(DebyePeaks['debyeTerms']):
2273                line = ' term'+'%2d'%(j)+':'
2274                for i in range(3):
2275                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2276                print >>pFile,line
2277        if DebyePeaks['nPeaks']:
2278            print >>pFile,'\n Single peak coefficients'
2279            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2280            line = ' names :  '
2281            for parm in parms:
2282                line += '%8s refine?'%(parm)
2283            print >>pFile,line
2284            for j,term in enumerate(DebyePeaks['peaksList']):
2285                line = ' peak'+'%2d'%(j)+':'
2286                for i in range(4):
2287                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2288                print >>pFile,line
2289       
2290    def PrintInstParms(Inst):
2291        print >>pFile,'\n Instrument Parameters:'
2292        insKeys = Inst.keys()
2293        insKeys.sort()
2294        iBeg = 0
2295        Ok = True
2296        while Ok:
2297            ptlbls = ' name  :'
2298            ptstr =  ' value :'
2299            varstr = ' refine:'
2300            iFin = min(iBeg+9,len(insKeys))
2301            for item in insKeys[iBeg:iFin]:
2302                if item not in ['Type','Source']:
2303                    ptlbls += '%12s' % (item)
2304                    ptstr += '%12.6f' % (Inst[item][1])
2305                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2306                        varstr += 12*' '
2307                    else:
2308                        varstr += '%12s' % (str(bool(Inst[item][2])))
2309            print >>pFile,ptlbls
2310            print >>pFile,ptstr
2311            print >>pFile,varstr
2312            iBeg = iFin
2313            if iBeg == len(insKeys):
2314                Ok = False
2315            else:
2316                print >>pFile,'\n'
2317       
2318    def PrintSampleParms(Sample):
2319        print >>pFile,'\n Sample Parameters:'
2320        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2321            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2322        ptlbls = ' name  :'
2323        ptstr =  ' value :'
2324        varstr = ' refine:'
2325        if 'Bragg' in Sample['Type']:
2326            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2327                ptlbls += '%14s'%(item)
2328                ptstr += '%14.4f'%(Sample[item][0])
2329                varstr += '%14s'%(str(bool(Sample[item][1])))
2330           
2331        elif 'Debye' in Type:        #Debye-Scherrer
2332            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2333                ptlbls += '%14s'%(item)
2334                ptstr += '%14.4f'%(Sample[item][0])
2335                varstr += '%14s'%(str(bool(Sample[item][1])))
2336
2337        print >>pFile,ptlbls
2338        print >>pFile,ptstr
2339        print >>pFile,varstr
2340       
2341    histDict = {}
2342    histVary = []
2343    controlDict = {}
2344    histoList = Histograms.keys()
2345    histoList.sort()
2346    for histogram in histoList:
2347        if 'PWDR' in histogram:
2348            Histogram = Histograms[histogram]
2349            hId = Histogram['hId']
2350            pfx = ':'+str(hId)+':'
2351            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2352            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2353            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2354            for excl in controlDict[pfx+'Exclude']:
2355                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2356            if controlDict[pfx+'Exclude']:
2357                ma.mask_rows(Histogram['Data'])
2358            Background = Histogram['Background']
2359            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2360            controlDict[pfx+'bakType'] = Type
2361            histDict.update(bakDict)
2362            histVary += bakVary
2363           
2364            Inst = Histogram['Instrument Parameters'][0]
2365            Type,instDict,insVary = GetInstParms(hId,Inst)
2366            controlDict[pfx+'histType'] = Type
2367            if 'XC' in Type:
2368                if pfx+'Lam1' in instDict:
2369                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2370                else:
2371                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2372            histDict.update(instDict)
2373            histVary += insVary
2374           
2375            Sample = Histogram['Sample Parameters']
2376            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2377            controlDict[pfx+'instType'] = Type
2378            histDict.update(sampDict)
2379            histVary += sampVary
2380           
2381   
2382            if Print: 
2383                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2384                print >>pFile,135*'-'
2385                Units = {'C':' deg','T':' msec'}
2386                units = Units[controlDict[pfx+'histType'][2]]
2387                Limits = controlDict[pfx+'Limits']
2388                print >>pFile,' Instrument type: ',Sample['Type']
2389                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2390                if len(controlDict[pfx+'Exclude']):
2391                    excls = controlDict[pfx+'Exclude']
2392                    for excl in excls:
2393                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
2394                PrintSampleParms(Sample)
2395                PrintInstParms(Inst)
2396                PrintBackground(Background)
2397        elif 'HKLF' in histogram:
2398            Histogram = Histograms[histogram]
2399            hId = Histogram['hId']
2400            pfx = ':'+str(hId)+':'
2401            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2402            Inst = Histogram['Instrument Parameters'][0]
2403            controlDict[pfx+'histType'] = Inst['Type'][0]
2404            if 'X' in Inst['Type'][0]:
2405                histDict[pfx+'Lam'] = Inst['Lam'][1]
2406                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2407    return histVary,histDict,controlDict
2408   
2409def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2410    'needs a doc string'
2411   
2412    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2413        Back = Background[0]
2414        DebyePeaks = Background[1]
2415        lenBack = len(Back[3:])
2416        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2417        for i in range(lenBack):
2418            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2419            if pfx+'Back:'+str(i) in sigDict:
2420                backSig[i] = sigDict[pfx+'Back:'+str(i)]
2421        if DebyePeaks['nDebye']:
2422            for i in range(DebyePeaks['nDebye']):
2423                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2424                for j,name in enumerate(names):
2425                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2426                    if name in sigDict:
2427                        backSig[lenBack+3*i+j] = sigDict[name]           
2428        if DebyePeaks['nPeaks']:
2429            for i in range(DebyePeaks['nPeaks']):
2430                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2431                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2432                for j,name in enumerate(names):
2433                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2434                    if name in sigDict:
2435                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2436        return backSig
2437       
2438    def SetInstParms(pfx,Inst,parmDict,sigDict):
2439        instSig = {}
2440        insKeys = Inst.keys()
2441        insKeys.sort()
2442        for item in insKeys:
2443            insName = pfx+item
2444            Inst[item][1] = parmDict[insName]
2445            if insName in sigDict:
2446                instSig[item] = sigDict[insName]
2447            else:
2448                instSig[item] = 0
2449        return instSig
2450       
2451    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2452        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2453            sampSig = [0 for i in range(5)]
2454            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2455                Sample[item][0] = parmDict[pfx+item]
2456                if pfx+item in sigDict:
2457                    sampSig[i] = sigDict[pfx+item]
2458        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2459            sampSig = [0 for i in range(4)]
2460            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2461                Sample[item][0] = parmDict[pfx+item]
2462                if pfx+item in sigDict:
2463                    sampSig[i] = sigDict[pfx+item]
2464        return sampSig
2465       
2466    def PrintBackgroundSig(Background,backSig):
2467        Back = Background[0]
2468        DebyePeaks = Background[1]
2469        lenBack = len(Back[3:])
2470        valstr = ' value : '
2471        sigstr = ' sig   : '
2472        refine = False
2473        for i,back in enumerate(Back[3:]):
2474            valstr += '%10.4g'%(back)
2475            if Back[1]:
2476                refine = True
2477                sigstr += '%10.4g'%(backSig[i])
2478            else:
2479                sigstr += 10*' '
2480        if refine:
2481            print >>pFile,'\n Background function: ',Back[0]
2482            print >>pFile,valstr
2483            print >>pFile,sigstr
2484        if DebyePeaks['nDebye']:
2485            ifAny = False
2486            ptfmt = "%12.3f"
2487            names =  ' names :'
2488            ptstr =  ' values:'
2489            sigstr = ' esds  :'
2490            for item in sigDict:
2491                if 'Debye' in item:
2492                    ifAny = True
2493                    names += '%12s'%(item)
2494                    ptstr += ptfmt%(parmDict[item])
2495                    sigstr += ptfmt%(sigDict[item])
2496            if ifAny:
2497                print >>pFile,'\n Debye diffuse scattering coefficients'
2498                print >>pFile,names
2499                print >>pFile,ptstr
2500                print >>pFile,sigstr
2501        if DebyePeaks['nPeaks']:
2502            ifAny = False
2503            ptfmt = "%14.3f"
2504            names =  ' names :'
2505            ptstr =  ' values:'
2506            sigstr = ' esds  :'
2507            for item in sigDict:
2508                if 'BkPk' in item:
2509                    ifAny = True
2510                    names += '%14s'%(item)
2511                    ptstr += ptfmt%(parmDict[item])
2512                    sigstr += ptfmt%(sigDict[item])
2513            if ifAny:
2514                print >>pFile,'\n Single peak coefficients'
2515                print >>pFile,names
2516                print >>pFile,ptstr
2517                print >>pFile,sigstr
2518       
2519    def PrintInstParmsSig(Inst,instSig):
2520        refine = False
2521        insKeys = instSig.keys()
2522        insKeys.sort()
2523        iBeg = 0
2524        Ok = True
2525        while Ok:
2526            ptlbls = ' names :'
2527            ptstr =  ' value :'
2528            sigstr = ' sig   :'
2529            iFin = min(iBeg+9,len(insKeys))
2530            for name in insKeys[iBeg:iFin]:
2531                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
2532                    ptlbls += '%12s' % (name)
2533                    ptstr += '%12.6f' % (Inst[name][1])
2534                    if instSig[name]:
2535                        refine = True
2536                        sigstr += '%12.6f' % (instSig[name])
2537                    else:
2538                        sigstr += 12*' '
2539            if refine:
2540                print >>pFile,'\n Instrument Parameters:'
2541                print >>pFile,ptlbls
2542                print >>pFile,ptstr
2543                print >>pFile,sigstr
2544            iBeg = iFin
2545            if iBeg == len(insKeys):
2546                Ok = False
2547       
2548    def PrintSampleParmsSig(Sample,sampleSig):
2549        ptlbls = ' names :'
2550        ptstr =  ' values:'
2551        sigstr = ' sig   :'
2552        refine = False
2553        if 'Bragg' in Sample['Type']:
2554            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2555                ptlbls += '%14s'%(item)
2556                ptstr += '%14.4f'%(Sample[item][0])
2557                if sampleSig[i]:
2558                    refine = True
2559                    sigstr += '%14.4f'%(sampleSig[i])
2560                else:
2561                    sigstr += 14*' '
2562           
2563        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2564            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2565                ptlbls += '%14s'%(item)
2566                ptstr += '%14.4f'%(Sample[item][0])
2567                if sampleSig[i]:
2568                    refine = True
2569                    sigstr += '%14.4f'%(sampleSig[i])
2570                else:
2571                    sigstr += 14*' '
2572
2573        if refine:
2574            print >>pFile,'\n Sample Parameters:'
2575            print >>pFile,ptlbls
2576            print >>pFile,ptstr
2577            print >>pFile,sigstr
2578       
2579    histoList = Histograms.keys()
2580    histoList.sort()
2581    for histogram in histoList:
2582        if 'PWDR' in histogram:
2583            Histogram = Histograms[histogram]
2584            hId = Histogram['hId']
2585            pfx = ':'+str(hId)+':'
2586            Background = Histogram['Background']
2587            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2588           
2589            Inst = Histogram['Instrument Parameters'][0]
2590            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2591       
2592            Sample = Histogram['Sample Parameters']
2593            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2594
2595            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2596            print >>pFile,135*'-'
2597            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2598            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2599            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2600                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2601            if Print:
2602                print >>pFile,' Instrument type: ',Sample['Type']
2603                PrintSampleParmsSig(Sample,sampSig)
2604                PrintInstParmsSig(Inst,instSig)
2605                PrintBackgroundSig(Background,backSig)
2606               
Note: See TracBrowser for help on using the repository browser.