source: trunk/GSASIIstrIO.py @ 1242

Last change on this file since 1242 was 1242, checked in by vondreele, 8 years ago

minor mods to sequential refinement print controls output

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