source: trunk/GSASIIstrIO.py @ 1409

Last change on this file since 1409 was 1409, checked in by toby, 8 years ago

more problems with older GPX files

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