source: trunk/GSASIIstrIO.py @ 1459

Last change on this file since 1459 was 1459, checked in by vondreele, 7 years ago

add instrument parameters (flight path & detector 2-theta) needed for TOF
rework reflection list for TOF
change default diff curve & reflection marker offsets
change weighting to instrument constants calibration to be 1/esd2 from peak fit positions - works a lot better
1st shot at TOF Rietveld refinement with derivatives - need to be checked for correctness (some are wrong)

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