source: trunk/GSASIIstrIO.py @ 1496

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

Add sequential peak fitting.
Change Back:n to Back;n for background parameters
Also DebyeA, etc.

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