source: trunk/GSASIIstrIO.py @ 1312

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

Parametric fit; work on PseudoVars? (more 2 come); minor code cleanups

  • 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-05-01 02:42:51 +0000 (Thu, 01 May 2014) $
4# $Author: toby $
5# $Revision: 1312 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 1312 2014-05-01 02:42:51Z 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: 1312 $")
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    return [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol]  #exchange sig(alp) & sig(gam) to get in right order
1281   
1282def SetPhaseData(parmDict,sigDict,Phases,RBIds,covData,RestraintDict=None,pFile=None):
1283    'needs a doc string'
1284   
1285    def PrintAtomsAndSig(General,Atoms,atomsSig):
1286        print >>pFile,'\n Atoms:'
1287        line = '   name      x         y         z      frac   Uiso     U11     U22     U33     U12     U13     U23'
1288        if General['Type'] == 'magnetic':
1289            line += '   Mx     My     Mz'
1290        elif General['Type'] == 'macromolecular':
1291            line = ' res no  residue  chain '+line
1292        print >>pFile,line
1293        if General['Type'] == 'nuclear':
1294            print >>pFile,135*'-'
1295            fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f',
1296                11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'}
1297            noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']}
1298            for atyp in General['AtomTypes']:       #zero composition
1299                General['NoAtoms'][atyp] = 0.
1300            for i,at in enumerate(Atoms):
1301                General['NoAtoms'][at[1]] += at[6]*float(at[8])     #new composition
1302                name = fmt[0]%(at[0])+fmt[1]%(at[1])+':'
1303                valstr = ' values:'
1304                sigstr = ' sig   :'
1305                for ind in [3,4,5,6]:
1306                    sigind = str(i)+':'+str(ind)
1307                    valstr += fmt[ind]%(at[ind])                   
1308                    if sigind in atomsSig:
1309                        sigstr += fmt[ind]%(atomsSig[sigind])
1310                    else:
1311                        sigstr += noFXsig[ind][1]%(noFXsig[ind][0])
1312                if at[9] == 'I':
1313                    valstr += fmt[10]%(at[10])
1314                    if str(i)+':10' in atomsSig:
1315                        sigstr += fmt[10]%(atomsSig[str(i)+':10'])
1316                    else:
1317                        sigstr += 8*' '
1318                else:
1319                    valstr += 8*' '
1320                    sigstr += 8*' '
1321                    for ind in [11,12,13,14,15,16]:
1322                        sigind = str(i)+':'+str(ind)
1323                        valstr += fmt[ind]%(at[ind])
1324                        if sigind in atomsSig:                       
1325                            sigstr += fmt[ind]%(atomsSig[sigind])
1326                        else:
1327                            sigstr += 8*' '
1328                print >>pFile,name
1329                print >>pFile,valstr
1330                print >>pFile,sigstr
1331               
1332    def PrintRBObjPOAndSig(rbfx,rbsx):
1333        namstr = '  names :'
1334        valstr = '  values:'
1335        sigstr = '  esds  :'
1336        for i,px in enumerate(['Px:','Py:','Pz:']):
1337            name = pfx+rbfx+px+rbsx
1338            namstr += '%12s'%('Pos '+px[1])
1339            valstr += '%12.5f'%(parmDict[name])
1340            if name in sigDict:
1341                sigstr += '%12.5f'%(sigDict[name])
1342            else:
1343                sigstr += 12*' '
1344        for i,po in enumerate(['Oa:','Oi:','Oj:','Ok:']):
1345            name = pfx+rbfx+po+rbsx
1346            namstr += '%12s'%('Ori '+po[1])
1347            valstr += '%12.5f'%(parmDict[name])
1348            if name in sigDict:
1349                sigstr += '%12.5f'%(sigDict[name])
1350            else:
1351                sigstr += 12*' '
1352        print >>pFile,namstr
1353        print >>pFile,valstr
1354        print >>pFile,sigstr
1355       
1356    def PrintRBObjTLSAndSig(rbfx,rbsx,TLS):
1357        namstr = '  names :'
1358        valstr = '  values:'
1359        sigstr = '  esds  :'
1360        if 'N' not in TLS:
1361            print >>pFile,' Thermal motion:'
1362        if 'T' in TLS:
1363            for i,pt in enumerate(['T11:','T22:','T33:','T12:','T13:','T23:']):
1364                name = pfx+rbfx+pt+rbsx
1365                namstr += '%12s'%(pt[:3])
1366                valstr += '%12.5f'%(parmDict[name])
1367                if name in sigDict:
1368                    sigstr += '%12.5f'%(sigDict[name])
1369                else:
1370                    sigstr += 12*' '
1371            print >>pFile,namstr
1372            print >>pFile,valstr
1373            print >>pFile,sigstr
1374        if 'L' in TLS:
1375            namstr = '  names :'
1376            valstr = '  values:'
1377            sigstr = '  esds  :'
1378            for i,pt in enumerate(['L11:','L22:','L33:','L12:','L13:','L23:']):
1379                name = pfx+rbfx+pt+rbsx
1380                namstr += '%12s'%(pt[:3])
1381                valstr += '%12.3f'%(parmDict[name])
1382                if name in sigDict:
1383                    sigstr += '%12.3f'%(sigDict[name])
1384                else:
1385                    sigstr += 12*' '
1386            print >>pFile,namstr
1387            print >>pFile,valstr
1388            print >>pFile,sigstr
1389        if 'S' in TLS:
1390            namstr = '  names :'
1391            valstr = '  values:'
1392            sigstr = '  esds  :'
1393            for i,pt in enumerate(['S12:','S13:','S21:','S23:','S31:','S32:','SAA:','SBB:']):
1394                name = pfx+rbfx+pt+rbsx
1395                namstr += '%12s'%(pt[:3])
1396                valstr += '%12.4f'%(parmDict[name])
1397                if name in sigDict:
1398                    sigstr += '%12.4f'%(sigDict[name])
1399                else:
1400                    sigstr += 12*' '
1401            print >>pFile,namstr
1402            print >>pFile,valstr
1403            print >>pFile,sigstr
1404        if 'U' in TLS:
1405            name = pfx+rbfx+'U:'+rbsx
1406            namstr = '  names :'+'%12s'%('Uiso')
1407            valstr = '  values:'+'%12.5f'%(parmDict[name])
1408            if name in sigDict:
1409                sigstr = '  esds  :'+'%12.5f'%(sigDict[name])
1410            else:
1411                sigstr = '  esds  :'+12*' '
1412            print >>pFile,namstr
1413            print >>pFile,valstr
1414            print >>pFile,sigstr
1415       
1416    def PrintRBObjTorAndSig(rbsx):
1417        namstr = '  names :'
1418        valstr = '  values:'
1419        sigstr = '  esds  :'
1420        nTors = len(RBObj['Torsions'])
1421        if nTors:
1422            print >>pFile,' Torsions:'
1423            for it in range(nTors):
1424                name = pfx+'RBRTr;'+str(it)+':'+rbsx
1425                namstr += '%12s'%('Tor'+str(it))
1426                valstr += '%12.4f'%(parmDict[name])
1427                if name in sigDict:
1428                    sigstr += '%12.4f'%(sigDict[name])
1429            print >>pFile,namstr
1430            print >>pFile,valstr
1431            print >>pFile,sigstr
1432               
1433    def PrintSHtextureAndSig(textureData,SHtextureSig):
1434        print >>pFile,'\n Spherical harmonics texture: Order:' + str(textureData['Order'])
1435        names = ['omega','chi','phi']
1436        namstr = '  names :'
1437        ptstr =  '  values:'
1438        sigstr = '  esds  :'
1439        for name in names:
1440            namstr += '%12s'%(name)
1441            ptstr += '%12.3f'%(textureData['Sample '+name][1])
1442            if 'Sample '+name in SHtextureSig:
1443                sigstr += '%12.3f'%(SHtextureSig['Sample '+name])
1444            else:
1445                sigstr += 12*' '
1446        print >>pFile,namstr
1447        print >>pFile,ptstr
1448        print >>pFile,sigstr
1449        print >>pFile,'\n Texture coefficients:'
1450        namstr = '  names :'
1451        ptstr =  '  values:'
1452        sigstr = '  esds  :'
1453        SHcoeff = textureData['SH Coeff'][1]
1454        for name in SHcoeff:
1455            namstr += '%12s'%(name)
1456            ptstr += '%12.3f'%(SHcoeff[name])
1457            if name in SHtextureSig:
1458                sigstr += '%12.3f'%(SHtextureSig[name])
1459            else:
1460                sigstr += 12*' '
1461        print >>pFile,namstr
1462        print >>pFile,ptstr
1463        print >>pFile,sigstr
1464           
1465    print >>pFile,'\n Phases:'
1466    for phase in Phases:
1467        print >>pFile,' Result for phase: ',phase
1468        Phase = Phases[phase]
1469        General = Phase['General']
1470        SGData = General['SGData']
1471        Atoms = Phase['Atoms']
1472        AtLookup = G2mth.FillAtomLookUp(Atoms)
1473        cell = General['Cell']
1474        pId = Phase['pId']
1475        pfx = str(pId)+'::'
1476        if cell[0]:
1477            A,sigA = cellFill(pfx,SGData,parmDict,sigDict)
1478            cellSig = getCellEsd(pfx,SGData,A,covData)  #includes sigVol
1479            print >>pFile,' Reciprocal metric tensor: '
1480            ptfmt = "%15.9f"
1481            names = ['A11','A22','A33','A12','A13','A23']
1482            namstr = '  names :'
1483            ptstr =  '  values:'
1484            sigstr = '  esds  :'
1485            for name,a,siga in zip(names,A,sigA):
1486                namstr += '%15s'%(name)
1487                ptstr += ptfmt%(a)
1488                if siga:
1489                    sigstr += ptfmt%(siga)
1490                else:
1491                    sigstr += 15*' '
1492            print >>pFile,namstr
1493            print >>pFile,ptstr
1494            print >>pFile,sigstr
1495            cell[1:7] = G2lat.A2cell(A)
1496            cell[7] = G2lat.calc_V(A)
1497            print >>pFile,' New unit cell:'
1498            ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"]
1499            names = ['a','b','c','alpha','beta','gamma','Volume']
1500            namstr = '  names :'
1501            ptstr =  '  values:'
1502            sigstr = '  esds  :'
1503            for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig):
1504                namstr += '%12s'%(name)
1505                ptstr += fmt%(a)
1506                if siga:
1507                    sigstr += fmt%(siga)
1508                else:
1509                    sigstr += 12*' '
1510            print >>pFile,namstr
1511            print >>pFile,ptstr
1512            print >>pFile,sigstr
1513           
1514        General['Mass'] = 0.
1515        if Phase['General'].get('doPawley'):
1516            pawleyRef = Phase['Pawley ref']
1517            for i,refl in enumerate(pawleyRef):
1518                key = pfx+'PWLref:'+str(i)
1519                refl[6] = parmDict[key]
1520                if key in sigDict:
1521                    refl[7] = sigDict[key]
1522                else:
1523                    refl[7] = 0
1524        else:
1525            VRBIds = RBIds['Vector']
1526            RRBIds = RBIds['Residue']
1527            RBModels = Phase['RBModels']
1528            for irb,RBObj in enumerate(RBModels.get('Vector',[])):
1529                jrb = VRBIds.index(RBObj['RBId'])
1530                rbsx = str(irb)+':'+str(jrb)
1531                print >>pFile,' Vector rigid body parameters:'
1532                PrintRBObjPOAndSig('RBV',rbsx)
1533                PrintRBObjTLSAndSig('RBV',rbsx,RBObj['ThermalMotion'][0])
1534            for irb,RBObj in enumerate(RBModels.get('Residue',[])):
1535                jrb = RRBIds.index(RBObj['RBId'])
1536                rbsx = str(irb)+':'+str(jrb)
1537                print >>pFile,' Residue rigid body parameters:'
1538                PrintRBObjPOAndSig('RBR',rbsx)
1539                PrintRBObjTLSAndSig('RBR',rbsx,RBObj['ThermalMotion'][0])
1540                PrintRBObjTorAndSig(rbsx)
1541            atomsSig = {}
1542            if General['Type'] == 'nuclear':        #this needs macromolecular variant!
1543                for i,at in enumerate(Atoms):
1544                    names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i),
1545                        10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i),
1546                        14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)}
1547                    for ind in [3,4,5,6]:
1548                        at[ind] = parmDict[names[ind]]
1549                        if ind in [3,4,5]:
1550                            name = names[ind].replace('A','dA')
1551                        else:
1552                            name = names[ind]
1553                        if name in sigDict:
1554                            atomsSig[str(i)+':'+str(ind)] = sigDict[name]
1555                    if at[9] == 'I':
1556                        at[10] = parmDict[names[10]]
1557                        if names[10] in sigDict:
1558                            atomsSig[str(i)+':10'] = sigDict[names[10]]
1559                    else:
1560                        for ind in [11,12,13,14,15,16]:
1561                            at[ind] = parmDict[names[ind]]
1562                            if names[ind] in sigDict:
1563                                atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]]
1564                    ind = General['AtomTypes'].index(at[1])
1565                    General['Mass'] += General['AtomMass'][ind]*at[6]*at[8]
1566            PrintAtomsAndSig(General,Atoms,atomsSig)
1567       
1568        textureData = General['SH Texture']   
1569        if textureData['Order']:
1570            SHtextureSig = {}
1571            for name in ['omega','chi','phi']:
1572                aname = pfx+'SH '+name
1573                textureData['Sample '+name][1] = parmDict[aname]
1574                if aname in sigDict:
1575                    SHtextureSig['Sample '+name] = sigDict[aname]
1576            for name in textureData['SH Coeff'][1]:
1577                aname = pfx+name
1578                textureData['SH Coeff'][1][name] = parmDict[aname]
1579                if aname in sigDict:
1580                    SHtextureSig[name] = sigDict[aname]
1581            PrintSHtextureAndSig(textureData,SHtextureSig)
1582        if phase in RestraintDict:
1583            PrintRestraints(cell[1:7],SGData,General['AtomPtrs'],Atoms,AtLookup,
1584                textureData,RestraintDict[phase],pFile)
1585                   
1586################################################################################
1587##### Histogram & Phase data
1588################################################################################       
1589                   
1590def GetHistogramPhaseData(Phases,Histograms,Print=True,pFile=None,resetRefList=True):
1591    '''Loads the HAP histogram/phase information into dicts
1592
1593    :param dict Phases: phase information
1594    :param dict Histograms: Histogram information
1595    :param bool Print: prints information as it is read
1596    :param file pFile: file object to print to (the default, None causes printing to the console)
1597    :param bool resetRefList: Should the contents of the Reflection List be initialized
1598      on loading. The default, True, initializes the Reflection List as it is loaded.
1599
1600    :returns: (hapVary,hapDict,controlDict)
1601      * hapVary: list of refined variables
1602      * hapDict: dict with refined variables and their values
1603      * controlDict: dict with computation controls (?)
1604    '''
1605   
1606    def PrintSize(hapData):
1607        if hapData[0] in ['isotropic','uniaxial']:
1608            line = '\n Size model    : %9s'%(hapData[0])
1609            line += ' equatorial:'+'%12.5f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1610            if hapData[0] == 'uniaxial':
1611                line += ' axial:'+'%12.5f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1612            line += '\n\t LG mixing coeff.: %12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1613            print >>pFile,line
1614        else:
1615            print >>pFile,'\n Size model    : %s'%(hapData[0])+ \
1616                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1617            Snames = ['S11','S22','S33','S12','S13','S23']
1618            ptlbls = ' names :'
1619            ptstr =  ' values:'
1620            varstr = ' refine:'
1621            for i,name in enumerate(Snames):
1622                ptlbls += '%12s' % (name)
1623                ptstr += '%12.6f' % (hapData[4][i])
1624                varstr += '%12s' % (str(hapData[5][i]))
1625            print >>pFile,ptlbls
1626            print >>pFile,ptstr
1627            print >>pFile,varstr
1628       
1629    def PrintMuStrain(hapData,SGData):
1630        if hapData[0] in ['isotropic','uniaxial']:
1631            line = '\n Mustrain model: %9s'%(hapData[0])
1632            line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0])
1633            if hapData[0] == 'uniaxial':
1634                line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1])
1635            line +='\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1636            print >>pFile,line
1637        else:
1638            print >>pFile,'\n Mustrain model: %s'%(hapData[0])+ \
1639                '\n\t LG mixing coeff.:%12.4f'%(hapData[1][2])+' Refine? '+str(hapData[2][2])
1640            Snames = G2spc.MustrainNames(SGData)
1641            ptlbls = ' names :'
1642            ptstr =  ' values:'
1643            varstr = ' refine:'
1644            for i,name in enumerate(Snames):
1645                ptlbls += '%12s' % (name)
1646                ptstr += '%12.6f' % (hapData[4][i])
1647                varstr += '%12s' % (str(hapData[5][i]))
1648            print >>pFile,ptlbls
1649            print >>pFile,ptstr
1650            print >>pFile,varstr
1651
1652    def PrintHStrain(hapData,SGData):
1653        print >>pFile,'\n Hydrostatic/elastic strain: '
1654        Hsnames = G2spc.HStrainNames(SGData)
1655        ptlbls = ' names :'
1656        ptstr =  ' values:'
1657        varstr = ' refine:'
1658        for i,name in enumerate(Hsnames):
1659            ptlbls += '%12s' % (name)
1660            ptstr += '%12.6f' % (hapData[0][i])
1661            varstr += '%12s' % (str(hapData[1][i]))
1662        print >>pFile,ptlbls
1663        print >>pFile,ptstr
1664        print >>pFile,varstr
1665
1666    def PrintSHPO(hapData):
1667        print >>pFile,'\n Spherical harmonics preferred orientation: Order:' + \
1668            str(hapData[4])+' Refine? '+str(hapData[2])
1669        ptlbls = ' names :'
1670        ptstr =  ' values:'
1671        for item in hapData[5]:
1672            ptlbls += '%12s'%(item)
1673            ptstr += '%12.3f'%(hapData[5][item]) 
1674        print >>pFile,ptlbls
1675        print >>pFile,ptstr
1676   
1677    def PrintBabinet(hapData):
1678        print >>pFile,'\n Babinet form factor modification: '
1679        ptlbls = ' names :'
1680        ptstr =  ' values:'
1681        varstr = ' refine:'
1682        for name in ['BabA','BabU']:
1683            ptlbls += '%12s' % (name)
1684            ptstr += '%12.6f' % (hapData[name][0])
1685            varstr += '%12s' % (str(hapData[name][1]))
1686        print >>pFile,ptlbls
1687        print >>pFile,ptstr
1688        print >>pFile,varstr
1689       
1690   
1691    hapDict = {}
1692    hapVary = []
1693    controlDict = {}
1694    poType = {}
1695    poAxes = {}
1696    spAxes = {}
1697    spType = {}
1698   
1699    for phase in Phases:
1700        HistoPhase = Phases[phase]['Histograms']
1701        SGData = Phases[phase]['General']['SGData']
1702        cell = Phases[phase]['General']['Cell'][1:7]
1703        A = G2lat.cell2A(cell)
1704        pId = Phases[phase]['pId']
1705        histoList = HistoPhase.keys()
1706        histoList.sort()
1707        for histogram in histoList:
1708            try:
1709                Histogram = Histograms[histogram]
1710            except KeyError:                       
1711                #skip if histogram not included e.g. in a sequential refinement
1712                continue
1713            hapData = HistoPhase[histogram]
1714            hId = Histogram['hId']
1715            if 'PWDR' in histogram:
1716                limits = Histogram['Limits'][1]
1717                inst = Histogram['Instrument Parameters'][0]
1718                Zero = inst['Zero'][1]
1719                if 'C' in inst['Type'][1]:
1720                    try:
1721                        wave = inst['Lam'][1]
1722                    except KeyError:
1723                        wave = inst['Lam1'][1]
1724                    dmin = wave/(2.0*sind(limits[1]/2.0))
1725                pfx = str(pId)+':'+str(hId)+':'
1726                for item in ['Scale','Extinction']:
1727                    hapDict[pfx+item] = hapData[item][0]
1728                    if hapData[item][1]:
1729                        hapVary.append(pfx+item)
1730                names = G2spc.HStrainNames(SGData)
1731                for i,name in enumerate(names):
1732                    hapDict[pfx+name] = hapData['HStrain'][0][i]
1733                    if hapData['HStrain'][1][i]:
1734                        hapVary.append(pfx+name)
1735                controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0]
1736                if hapData['Pref.Ori.'][0] == 'MD':
1737                    hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1]
1738                    controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3]
1739                    if hapData['Pref.Ori.'][2]:
1740                        hapVary.append(pfx+'MD')
1741                else:                           #'SH' spherical harmonics
1742                    controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4]
1743                    controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5])
1744                    for item in hapData['Pref.Ori.'][5]:
1745                        hapDict[pfx+item] = hapData['Pref.Ori.'][5][item]
1746                        if hapData['Pref.Ori.'][2]:
1747                            hapVary.append(pfx+item)
1748                for item in ['Mustrain','Size']:
1749                    controlDict[pfx+item+'Type'] = hapData[item][0]
1750                    hapDict[pfx+item+';mx'] = hapData[item][1][2]
1751                    if hapData[item][2][2]:
1752                        hapVary.append(pfx+item+';mx')
1753                    if hapData[item][0] in ['isotropic','uniaxial']:
1754                        hapDict[pfx+item+';i'] = hapData[item][1][0]
1755                        if hapData[item][2][0]:
1756                            hapVary.append(pfx+item+';i')
1757                        if hapData[item][0] == 'uniaxial':
1758                            controlDict[pfx+item+'Axis'] = hapData[item][3]
1759                            hapDict[pfx+item+';a'] = hapData[item][1][1]
1760                            if hapData[item][2][1]:
1761                                hapVary.append(pfx+item+';a')
1762                    else:       #generalized for mustrain or ellipsoidal for size
1763                        Nterms = len(hapData[item][4])
1764                        if item == 'Mustrain':
1765                            names = G2spc.MustrainNames(SGData)
1766                            pwrs = []
1767                            for name in names:
1768                                h,k,l = name[1:]
1769                                pwrs.append([int(h),int(k),int(l)])
1770                            controlDict[pfx+'MuPwrs'] = pwrs
1771                        for i in range(Nterms):
1772                            sfx = ':'+str(i)
1773                            hapDict[pfx+item+sfx] = hapData[item][4][i]
1774                            if hapData[item][5][i]:
1775                                hapVary.append(pfx+item+sfx)
1776                for bab in ['BabA','BabU']:
1777                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1778                    if hapData['Babinet'][bab][1]:
1779                        hapVary.append(pfx+bab)
1780                               
1781                if Print: 
1782                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1783                    print >>pFile,135*'-'
1784                    print >>pFile,' Phase fraction  : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1785                    print >>pFile,' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1]
1786                    if hapData['Pref.Ori.'][0] == 'MD':
1787                        Ax = hapData['Pref.Ori.'][3]
1788                        print >>pFile,' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \
1789                            ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2])
1790                    else: #'SH' for spherical harmonics
1791                        PrintSHPO(hapData['Pref.Ori.'])
1792                    PrintSize(hapData['Size'])
1793                    PrintMuStrain(hapData['Mustrain'],SGData)
1794                    PrintHStrain(hapData['HStrain'],SGData)
1795                    if hapData['Babinet']['BabA'][0]:
1796                        PrintBabinet(hapData['Babinet'])
1797                HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A))
1798                if resetRefList:
1799                    refList = []
1800                    Uniq = []
1801                    Phi = []
1802                    for h,k,l,d in HKLd:
1803                        ext,mul,uniq,phi = G2spc.GenHKLf([h,k,l],SGData)
1804                        mul *= 2      # for powder overlap of Friedel pairs
1805                        if ext:
1806                            continue
1807                        if 'C' in inst['Type'][0]:
1808                            pos = 2.0*asind(wave/(2.0*d))+Zero
1809                            if limits[0] < pos < limits[1]:
1810                                refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,0.0])
1811                                Uniq.append(uniq)
1812                                Phi.append(phi)
1813                        else:
1814                            raise ValueError 
1815                    Histogram['Reflection Lists'][phase] = {'RefList':np.array(refList),'FF':{}}
1816            elif 'HKLF' in histogram:
1817                inst = Histogram['Instrument Parameters'][0]
1818                hId = Histogram['hId']
1819                hfx = ':%d:'%(hId)
1820                for item in inst:
1821                    if type(inst) is not list and item != 'Type': continue # skip over non-refined items (such as InstName)
1822                    hapDict[hfx+item] = inst[item][1]
1823                pfx = str(pId)+':'+str(hId)+':'
1824                hapDict[pfx+'Scale'] = hapData['Scale'][0]
1825                if hapData['Scale'][1]:
1826                    hapVary.append(pfx+'Scale')
1827                               
1828                extApprox,extType,extParms = hapData['Extinction']
1829                controlDict[pfx+'EType'] = extType
1830                controlDict[pfx+'EApprox'] = extApprox
1831                controlDict[pfx+'Tbar'] = extParms['Tbar']
1832                controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1833                if 'Primary' in extType:
1834                    Ekey = ['Ep',]
1835                elif 'I & II' in extType:
1836                    Ekey = ['Eg','Es']
1837                elif 'Secondary Type II' == extType:
1838                    Ekey = ['Es',]
1839                elif 'Secondary Type I' == extType:
1840                    Ekey = ['Eg',]
1841                else:   #'None'
1842                    Ekey = []
1843                for eKey in Ekey:
1844                    hapDict[pfx+eKey] = extParms[eKey][0]
1845                    if extParms[eKey][1]:
1846                        hapVary.append(pfx+eKey)
1847                for bab in ['BabA','BabU']:
1848                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1849                    if hapData['Babinet'][bab][1]:
1850                        hapVary.append(pfx+bab)
1851                if Print: 
1852                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1853                    print >>pFile,135*'-'
1854                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1855                    if extType != 'None':
1856                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1857                        text = ' Parameters       :'
1858                        for eKey in Ekey:
1859                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1860                        print >>pFile,text
1861                    if hapData['Babinet']['BabA'][0]:
1862                        PrintBabinet(hapData['Babinet'])
1863                Histogram['Reflection Lists'] = phase       
1864               
1865    return hapVary,hapDict,controlDict
1866   
1867def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1868    'needs a doc string'
1869   
1870    def PrintSizeAndSig(hapData,sizeSig):
1871        line = '\n Size model:     %9s'%(hapData[0])
1872        refine = False
1873        if hapData[0] in ['isotropic','uniaxial']:
1874            line += ' equatorial:%12.5f'%(hapData[1][0])
1875            if sizeSig[0][0]:
1876                line += ', sig:%8.4f'%(sizeSig[0][0])
1877                refine = True
1878            if hapData[0] == 'uniaxial':
1879                line += ' axial:%12.4f'%(hapData[1][1])
1880                if sizeSig[0][1]:
1881                    refine = True
1882                    line += ', sig:%8.4f'%(sizeSig[0][1])
1883            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1884            if sizeSig[0][2]:
1885                refine = True
1886                line += ', sig:%8.4f'%(sizeSig[0][2])
1887            if refine:
1888                print >>pFile,line
1889        else:
1890            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1891            if sizeSig[0][2]:
1892                refine = True
1893                line += ', sig:%8.4f'%(sizeSig[0][2])
1894            Snames = ['S11','S22','S33','S12','S13','S23']
1895            ptlbls = ' name  :'
1896            ptstr =  ' value :'
1897            sigstr = ' sig   :'
1898            for i,name in enumerate(Snames):
1899                ptlbls += '%12s' % (name)
1900                ptstr += '%12.6f' % (hapData[4][i])
1901                if sizeSig[1][i]:
1902                    refine = True
1903                    sigstr += '%12.6f' % (sizeSig[1][i])
1904                else:
1905                    sigstr += 12*' '
1906            if refine:
1907                print >>pFile,line
1908                print >>pFile,ptlbls
1909                print >>pFile,ptstr
1910                print >>pFile,sigstr
1911       
1912    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1913        line = '\n Mustrain model: %9s'%(hapData[0])
1914        refine = False
1915        if hapData[0] in ['isotropic','uniaxial']:
1916            line += ' equatorial:%12.1f'%(hapData[1][0])
1917            if mustrainSig[0][0]:
1918                line += ', sig:%8.1f'%(mustrainSig[0][0])
1919                refine = True
1920            if hapData[0] == 'uniaxial':
1921                line += ' axial:%12.1f'%(hapData[1][1])
1922                if mustrainSig[0][1]:
1923                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1924            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1925            if mustrainSig[0][2]:
1926                refine = True
1927                line += ', sig:%8.3f'%(mustrainSig[0][2])
1928            if refine:
1929                print >>pFile,line
1930        else:
1931            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1932            if mustrainSig[0][2]:
1933                refine = True
1934                line += ', sig:%8.3f'%(mustrainSig[0][2])
1935            Snames = G2spc.MustrainNames(SGData)
1936            ptlbls = ' name  :'
1937            ptstr =  ' value :'
1938            sigstr = ' sig   :'
1939            for i,name in enumerate(Snames):
1940                ptlbls += '%12s' % (name)
1941                ptstr += '%12.6f' % (hapData[4][i])
1942                if mustrainSig[1][i]:
1943                    refine = True
1944                    sigstr += '%12.6f' % (mustrainSig[1][i])
1945                else:
1946                    sigstr += 12*' '
1947            if refine:
1948                print >>pFile,line
1949                print >>pFile,ptlbls
1950                print >>pFile,ptstr
1951                print >>pFile,sigstr
1952           
1953    def PrintHStrainAndSig(hapData,strainSig,SGData):
1954        Hsnames = G2spc.HStrainNames(SGData)
1955        ptlbls = ' name  :'
1956        ptstr =  ' value :'
1957        sigstr = ' sig   :'
1958        refine = False
1959        for i,name in enumerate(Hsnames):
1960            ptlbls += '%12s' % (name)
1961            ptstr += '%12.6g' % (hapData[0][i])
1962            if name in strainSig:
1963                refine = True
1964                sigstr += '%12.6g' % (strainSig[name])
1965            else:
1966                sigstr += 12*' '
1967        if refine:
1968            print >>pFile,'\n Hydrostatic/elastic strain: '
1969            print >>pFile,ptlbls
1970            print >>pFile,ptstr
1971            print >>pFile,sigstr
1972       
1973    def PrintSHPOAndSig(pfx,hapData,POsig):
1974        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
1975        ptlbls = ' names :'
1976        ptstr =  ' values:'
1977        sigstr = ' sig   :'
1978        for item in hapData[5]:
1979            ptlbls += '%12s'%(item)
1980            ptstr += '%12.3f'%(hapData[5][item])
1981            if pfx+item in POsig:
1982                sigstr += '%12.3f'%(POsig[pfx+item])
1983            else:
1984                sigstr += 12*' ' 
1985        print >>pFile,ptlbls
1986        print >>pFile,ptstr
1987        print >>pFile,sigstr
1988       
1989    def PrintExtAndSig(pfx,hapData,ScalExtSig):
1990        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
1991        text = ''
1992        for item in hapData[2]:
1993            if pfx+item in ScalExtSig:
1994                text += '       %s: '%(item)
1995                text += '%12.2e'%(hapData[2][item][0])
1996                if pfx+item in ScalExtSig:
1997                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
1998        print >>pFile,text       
1999       
2000    def PrintBabinetAndSig(pfx,hapData,BabSig):
2001        print >>pFile,'\n Babinet form factor modification: '
2002        ptlbls = ' names :'
2003        ptstr =  ' values:'
2004        sigstr = ' sig   :'
2005        for item in hapData:
2006            ptlbls += '%12s'%(item)
2007            ptstr += '%12.3f'%(hapData[item][0])
2008            if pfx+item in BabSig:
2009                sigstr += '%12.3f'%(BabSig[pfx+item])
2010            else:
2011                sigstr += 12*' ' 
2012        print >>pFile,ptlbls
2013        print >>pFile,ptstr
2014        print >>pFile,sigstr
2015   
2016    PhFrExtPOSig = {}
2017    SizeMuStrSig = {}
2018    ScalExtSig = {}
2019    BabSig = {}
2020    wtFrSum = {}
2021    for phase in Phases:
2022        HistoPhase = Phases[phase]['Histograms']
2023        General = Phases[phase]['General']
2024        SGData = General['SGData']
2025        pId = Phases[phase]['pId']
2026        histoList = HistoPhase.keys()
2027        histoList.sort()
2028        for histogram in histoList:
2029            try:
2030                Histogram = Histograms[histogram]
2031            except KeyError:                       
2032                #skip if histogram not included e.g. in a sequential refinement
2033                continue
2034            hapData = HistoPhase[histogram]
2035            hId = Histogram['hId']
2036            pfx = str(pId)+':'+str(hId)+':'
2037            if hId not in wtFrSum:
2038                wtFrSum[hId] = 0.
2039            if 'PWDR' in histogram:
2040                for item in ['Scale','Extinction']:
2041                    hapData[item][0] = parmDict[pfx+item]
2042                    if pfx+item in sigDict:
2043                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2044                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2045                if hapData['Pref.Ori.'][0] == 'MD':
2046                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2047                    if pfx+'MD' in sigDict:
2048                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2049                else:                           #'SH' spherical harmonics
2050                    for item in hapData['Pref.Ori.'][5]:
2051                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2052                        if pfx+item in sigDict:
2053                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2054                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2055                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2056                    pfx+'HStrain':{}})                 
2057                for item in ['Mustrain','Size']:
2058                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2059                    hapData[item][1][2] = min(1.,max(0.1,hapData[item][1][2]))
2060                    if pfx+item+';mx' in sigDict:
2061                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2062                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2063                        hapData[item][1][0] = parmDict[pfx+item+';i']
2064                        if item == 'Size':
2065                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2066                        if pfx+item+';i' in sigDict: 
2067                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2068                        if hapData[item][0] == 'uniaxial':
2069                            hapData[item][1][1] = parmDict[pfx+item+';a']
2070                            if item == 'Size':
2071                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2072                            if pfx+item+';a' in sigDict:
2073                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2074                    else:       #generalized for mustrain or ellipsoidal for size
2075                        Nterms = len(hapData[item][4])
2076                        for i in range(Nterms):
2077                            sfx = ':'+str(i)
2078                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2079                            if pfx+item+sfx in sigDict:
2080                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2081                names = G2spc.HStrainNames(SGData)
2082                for i,name in enumerate(names):
2083                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2084                    if pfx+name in sigDict:
2085                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2086                for name in ['BabA','BabU']:
2087                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2088                    if pfx+name in sigDict:
2089                        BabSig[pfx+name] = sigDict[pfx+name]               
2090               
2091            elif 'HKLF' in histogram:
2092                for item in ['Scale',]:
2093                    if parmDict.get(pfx+item):
2094                        hapData[item][0] = parmDict[pfx+item]
2095                        if pfx+item in sigDict:
2096                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2097                for item in ['Ep','Eg','Es']:
2098                    if parmDict.get(pfx+item):
2099                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2100                        if pfx+item in sigDict:
2101                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2102                for name in ['BabA','BabU']:
2103                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2104                    if pfx+name in sigDict:
2105                        BabSig[pfx+name] = sigDict[pfx+name]               
2106
2107    if Print:
2108        for phase in Phases:
2109            HistoPhase = Phases[phase]['Histograms']
2110            General = Phases[phase]['General']
2111            SGData = General['SGData']
2112            pId = Phases[phase]['pId']
2113            histoList = HistoPhase.keys()
2114            histoList.sort()
2115            for histogram in histoList:
2116                try:
2117                    Histogram = Histograms[histogram]
2118                except KeyError:                       
2119                    #skip if histogram not included e.g. in a sequential refinement
2120                    continue
2121                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2122                print >>pFile,130*'-'
2123                hapData = HistoPhase[histogram]
2124                hId = Histogram['hId']
2125                Histogram['Residuals'][str(pId)+'::Name'] = phase
2126                pfx = str(pId)+':'+str(hId)+':'
2127                if 'PWDR' in histogram:
2128                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2129                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2130               
2131                    if pfx+'Scale' in PhFrExtPOSig:
2132                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2133                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2134                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2135                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2136                    if pfx+'Extinction' in PhFrExtPOSig:
2137                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2138                    if hapData['Pref.Ori.'][0] == 'MD':
2139                        if pfx+'MD' in PhFrExtPOSig:
2140                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2141                    else:
2142                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2143                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2144                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2145                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2146                    if len(BabSig):
2147                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2148                   
2149                elif 'HKLF' in histogram:
2150                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2151                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2152                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2153                    if pfx+'Scale' in ScalExtSig:
2154                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2155                    if hapData['Extinction'][0] != 'None':
2156                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2157                    if len(BabSig):
2158                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2159
2160################################################################################
2161##### Histogram data
2162################################################################################       
2163                   
2164def GetHistogramData(Histograms,Print=True,pFile=None):
2165    'needs a doc string'
2166   
2167    def GetBackgroundParms(hId,Background):
2168        Back = Background[0]
2169        DebyePeaks = Background[1]
2170        bakType,bakFlag = Back[:2]
2171        backVals = Back[3:]
2172        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2173        backDict = dict(zip(backNames,backVals))
2174        backVary = []
2175        if bakFlag:
2176            backVary = backNames
2177        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2178        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2179        debyeDict = {}
2180        debyeList = []
2181        for i in range(DebyePeaks['nDebye']):
2182            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2183            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2184            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2185        debyeVary = []
2186        for item in debyeList:
2187            if item[1]:
2188                debyeVary.append(item[0])
2189        backDict.update(debyeDict)
2190        backVary += debyeVary
2191        peakDict = {}
2192        peakList = []
2193        for i in range(DebyePeaks['nPeaks']):
2194            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2195                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2196            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2197            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2198        peakVary = []
2199        for item in peakList:
2200            if item[1]:
2201                peakVary.append(item[0])
2202        backDict.update(peakDict)
2203        backVary += peakVary
2204        return bakType,backDict,backVary           
2205       
2206    def GetInstParms(hId,Inst):     
2207        dataType = Inst['Type'][0]
2208        instDict = {}
2209        insVary = []
2210        pfx = ':'+str(hId)+':'
2211        insKeys = Inst.keys()
2212        insKeys.sort()
2213        for item in insKeys:
2214            insName = pfx+item
2215            instDict[insName] = Inst[item][1]
2216            if Inst[item][2]:
2217                insVary.append(insName)
2218#        instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001)
2219#        instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001)
2220        instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2221        return dataType,instDict,insVary
2222       
2223    def GetSampleParms(hId,Sample):
2224        sampVary = []
2225        hfx = ':'+str(hId)+':'       
2226        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2227            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2228        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2229            if key in Sample:
2230                sampDict[hfx+key] = Sample[key]
2231        Type = Sample['Type']
2232        if 'Bragg' in Type:             #Bragg-Brentano
2233            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2234                sampDict[hfx+item] = Sample[item][0]
2235                if Sample[item][1]:
2236                    sampVary.append(hfx+item)
2237        elif 'Debye' in Type:        #Debye-Scherrer
2238            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2239                sampDict[hfx+item] = Sample[item][0]
2240                if Sample[item][1]:
2241                    sampVary.append(hfx+item)
2242        return Type,sampDict,sampVary
2243       
2244    def PrintBackground(Background):
2245        Back = Background[0]
2246        DebyePeaks = Background[1]
2247        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2248        line = ' Coefficients: '
2249        for i,back in enumerate(Back[3:]):
2250            line += '%10.3f'%(back)
2251            if i and not i%10:
2252                line += '\n'+15*' '
2253        print >>pFile,line
2254        if DebyePeaks['nDebye']:
2255            print >>pFile,'\n Debye diffuse scattering coefficients'
2256            parms = ['DebyeA','DebyeR','DebyeU']
2257            line = ' names :  '
2258            for parm in parms:
2259                line += '%8s refine?'%(parm)
2260            print >>pFile,line
2261            for j,term in enumerate(DebyePeaks['debyeTerms']):
2262                line = ' term'+'%2d'%(j)+':'
2263                for i in range(3):
2264                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2265                print >>pFile,line
2266        if DebyePeaks['nPeaks']:
2267            print >>pFile,'\n Single peak coefficients'
2268            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2269            line = ' names :  '
2270            for parm in parms:
2271                line += '%8s refine?'%(parm)
2272            print >>pFile,line
2273            for j,term in enumerate(DebyePeaks['peaksList']):
2274                line = ' peak'+'%2d'%(j)+':'
2275                for i in range(4):
2276                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2277                print >>pFile,line
2278       
2279    def PrintInstParms(Inst):
2280        print >>pFile,'\n Instrument Parameters:'
2281        ptlbls = ' name  :'
2282        ptstr =  ' value :'
2283        varstr = ' refine:'
2284        insKeys = Inst.keys()
2285        insKeys.sort()
2286        for item in insKeys:
2287            if item != 'Type':
2288                ptlbls += '%12s' % (item)
2289                ptstr += '%12.6f' % (Inst[item][1])
2290                if item in ['Lam1','Lam2','Azimuth']:
2291                    varstr += 12*' '
2292                else:
2293                    varstr += '%12s' % (str(bool(Inst[item][2])))
2294        print >>pFile,ptlbls
2295        print >>pFile,ptstr
2296        print >>pFile,varstr
2297       
2298    def PrintSampleParms(Sample):
2299        print >>pFile,'\n Sample Parameters:'
2300        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2301            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2302        ptlbls = ' name  :'
2303        ptstr =  ' value :'
2304        varstr = ' refine:'
2305        if 'Bragg' in Sample['Type']:
2306            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2307                ptlbls += '%14s'%(item)
2308                ptstr += '%14.4f'%(Sample[item][0])
2309                varstr += '%14s'%(str(bool(Sample[item][1])))
2310           
2311        elif 'Debye' in Type:        #Debye-Scherrer
2312            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2313                ptlbls += '%14s'%(item)
2314                ptstr += '%14.4f'%(Sample[item][0])
2315                varstr += '%14s'%(str(bool(Sample[item][1])))
2316
2317        print >>pFile,ptlbls
2318        print >>pFile,ptstr
2319        print >>pFile,varstr
2320       
2321    histDict = {}
2322    histVary = []
2323    controlDict = {}
2324    histoList = Histograms.keys()
2325    histoList.sort()
2326    for histogram in histoList:
2327        if 'PWDR' in histogram:
2328            Histogram = Histograms[histogram]
2329            hId = Histogram['hId']
2330            pfx = ':'+str(hId)+':'
2331            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2332            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2333            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2334            for excl in controlDict[pfx+'Exclude']:
2335                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2336            if controlDict[pfx+'Exclude']:
2337                ma.mask_rows(Histogram['Data'])
2338            Background = Histogram['Background']
2339            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2340            controlDict[pfx+'bakType'] = Type
2341            histDict.update(bakDict)
2342            histVary += bakVary
2343           
2344            Inst = Histogram['Instrument Parameters'][0]
2345            Type,instDict,insVary = GetInstParms(hId,Inst)
2346            controlDict[pfx+'histType'] = Type
2347            if pfx+'Lam1' in instDict:
2348                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2349            else:
2350                controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2351            histDict.update(instDict)
2352            histVary += insVary
2353           
2354            Sample = Histogram['Sample Parameters']
2355            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2356            controlDict[pfx+'instType'] = Type
2357            histDict.update(sampDict)
2358            histVary += sampVary
2359           
2360   
2361            if Print: 
2362                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2363                print >>pFile,135*'-'
2364                Units = {'C':' deg','T':' msec'}
2365                units = Units[controlDict[pfx+'histType'][2]]
2366                Limits = controlDict[pfx+'Limits']
2367                print >>pFile,' Instrument type: ',Sample['Type']
2368                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2369                if len(controlDict[pfx+'Exclude']):
2370                    excls = controlDict[pfx+'Exclude']
2371                    for excl in excls:
2372                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
2373                PrintSampleParms(Sample)
2374                PrintInstParms(Inst)
2375                PrintBackground(Background)
2376        elif 'HKLF' in histogram:
2377            Histogram = Histograms[histogram]
2378            hId = Histogram['hId']
2379            pfx = ':'+str(hId)+':'
2380            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2381            Inst = Histogram['Instrument Parameters'][0]
2382            controlDict[pfx+'histType'] = Inst['Type'][0]
2383            histDict[pfx+'Lam'] = Inst['Lam'][1]
2384            controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2385    return histVary,histDict,controlDict
2386   
2387def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2388    'needs a doc string'
2389   
2390    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2391        Back = Background[0]
2392        DebyePeaks = Background[1]
2393        lenBack = len(Back[3:])
2394        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2395        for i in range(lenBack):
2396            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2397            if pfx+'Back:'+str(i) in sigDict:
2398                backSig[i] = sigDict[pfx+'Back:'+str(i)]
2399        if DebyePeaks['nDebye']:
2400            for i in range(DebyePeaks['nDebye']):
2401                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2402                for j,name in enumerate(names):
2403                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2404                    if name in sigDict:
2405                        backSig[lenBack+3*i+j] = sigDict[name]           
2406        if DebyePeaks['nPeaks']:
2407            for i in range(DebyePeaks['nPeaks']):
2408                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2409                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2410                for j,name in enumerate(names):
2411                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2412                    if name in sigDict:
2413                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2414        return backSig
2415       
2416    def SetInstParms(pfx,Inst,parmDict,sigDict):
2417        instSig = {}
2418        insKeys = Inst.keys()
2419        insKeys.sort()
2420        for item in insKeys:
2421            insName = pfx+item
2422            Inst[item][1] = parmDict[insName]
2423            if insName in sigDict:
2424                instSig[item] = sigDict[insName]
2425            else:
2426                instSig[item] = 0
2427        return instSig
2428       
2429    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2430        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2431            sampSig = [0 for i in range(5)]
2432            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2433                Sample[item][0] = parmDict[pfx+item]
2434                if pfx+item in sigDict:
2435                    sampSig[i] = sigDict[pfx+item]
2436        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2437            sampSig = [0 for i in range(4)]
2438            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2439                Sample[item][0] = parmDict[pfx+item]
2440                if pfx+item in sigDict:
2441                    sampSig[i] = sigDict[pfx+item]
2442        return sampSig
2443       
2444    def PrintBackgroundSig(Background,backSig):
2445        Back = Background[0]
2446        DebyePeaks = Background[1]
2447        lenBack = len(Back[3:])
2448        valstr = ' value : '
2449        sigstr = ' sig   : '
2450        refine = False
2451        for i,back in enumerate(Back[3:]):
2452            valstr += '%10.4g'%(back)
2453            if Back[1]:
2454                refine = True
2455                sigstr += '%10.4g'%(backSig[i])
2456            else:
2457                sigstr += 10*' '
2458        if refine:
2459            print >>pFile,'\n Background function: ',Back[0]
2460            print >>pFile,valstr
2461            print >>pFile,sigstr
2462        if DebyePeaks['nDebye']:
2463            ifAny = False
2464            ptfmt = "%12.3f"
2465            names =  ' names :'
2466            ptstr =  ' values:'
2467            sigstr = ' esds  :'
2468            for item in sigDict:
2469                if 'Debye' in item:
2470                    ifAny = True
2471                    names += '%12s'%(item)
2472                    ptstr += ptfmt%(parmDict[item])
2473                    sigstr += ptfmt%(sigDict[item])
2474            if ifAny:
2475                print >>pFile,'\n Debye diffuse scattering coefficients'
2476                print >>pFile,names
2477                print >>pFile,ptstr
2478                print >>pFile,sigstr
2479        if DebyePeaks['nPeaks']:
2480            ifAny = False
2481            ptfmt = "%14.3f"
2482            names =  ' names :'
2483            ptstr =  ' values:'
2484            sigstr = ' esds  :'
2485            for item in sigDict:
2486                if 'BkPk' in item:
2487                    ifAny = True
2488                    names += '%14s'%(item)
2489                    ptstr += ptfmt%(parmDict[item])
2490                    sigstr += ptfmt%(sigDict[item])
2491            if ifAny:
2492                print >>pFile,'\n Single peak coefficients'
2493                print >>pFile,names
2494                print >>pFile,ptstr
2495                print >>pFile,sigstr
2496       
2497    def PrintInstParmsSig(Inst,instSig):
2498        ptlbls = ' names :'
2499        ptstr =  ' value :'
2500        sigstr = ' sig   :'
2501        refine = False
2502        insKeys = instSig.keys()
2503        insKeys.sort()
2504        for name in insKeys:
2505            if name not in  ['Type','Lam1','Lam2','Azimuth']:
2506                ptlbls += '%12s' % (name)
2507                ptstr += '%12.6f' % (Inst[name][1])
2508                if instSig[name]:
2509                    refine = True
2510                    sigstr += '%12.6f' % (instSig[name])
2511                else:
2512                    sigstr += 12*' '
2513        if refine:
2514            print >>pFile,'\n Instrument Parameters:'
2515            print >>pFile,ptlbls
2516            print >>pFile,ptstr
2517            print >>pFile,sigstr
2518       
2519    def PrintSampleParmsSig(Sample,sampleSig):
2520        ptlbls = ' names :'
2521        ptstr =  ' values:'
2522        sigstr = ' sig   :'
2523        refine = False
2524        if 'Bragg' in Sample['Type']:
2525            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2526                ptlbls += '%14s'%(item)
2527                ptstr += '%14.4f'%(Sample[item][0])
2528                if sampleSig[i]:
2529                    refine = True
2530                    sigstr += '%14.4f'%(sampleSig[i])
2531                else:
2532                    sigstr += 14*' '
2533           
2534        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2535            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2536                ptlbls += '%14s'%(item)
2537                ptstr += '%14.4f'%(Sample[item][0])
2538                if sampleSig[i]:
2539                    refine = True
2540                    sigstr += '%14.4f'%(sampleSig[i])
2541                else:
2542                    sigstr += 14*' '
2543
2544        if refine:
2545            print >>pFile,'\n Sample Parameters:'
2546            print >>pFile,ptlbls
2547            print >>pFile,ptstr
2548            print >>pFile,sigstr
2549       
2550    histoList = Histograms.keys()
2551    histoList.sort()
2552    for histogram in histoList:
2553        if 'PWDR' in histogram:
2554            Histogram = Histograms[histogram]
2555            hId = Histogram['hId']
2556            pfx = ':'+str(hId)+':'
2557            Background = Histogram['Background']
2558            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2559           
2560            Inst = Histogram['Instrument Parameters'][0]
2561            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2562       
2563            Sample = Histogram['Sample Parameters']
2564            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2565
2566            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2567            print >>pFile,135*'-'
2568            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2569            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2570            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2571                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2572            if Print:
2573                print >>pFile,' Instrument type: ',Sample['Type']
2574                PrintSampleParmsSig(Sample,sampSig)
2575                PrintInstParmsSig(Inst,instSig)
2576                PrintBackgroundSig(Background,backSig)
2577               
Note: See TracBrowser for help on using the repository browser.