source: trunk/GSASIIstrIO.py @ 1106

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

modify reflection array to be a dictionary with 4 items: RefList?, Uniq, Phi & FF. Each is list. Use patches to convert old format to new in various places.
Fix bug in simulate PWDR code.

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