source: trunk/GSASIIstrIO.py @ 1299

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

Set conf flags consistently for .py files

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