source: trunk/GSASIIstrIO.py @ 1107

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

make reflist np.array
fix bug - hfx+'Type' needs to be in parmDict

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