source: trunk/GSASIIstrIO.py @ 1523

Last change on this file since 1523 was 1523, checked in by toby, 7 years ago

change isinstance in patch code

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