source: trunk/GSASIIstrIO.py @ 1078

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

fix a single crystal blunder

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