source: trunk/GSASIIstrIO.py @ 1138

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

major constraints revision

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