source: trunk/GSASIIstrIO.py @ 1489

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

correct Uij to U6 conversions Uij = Uk/2 for i != j; k > 2 in G2lattice & G2spc
This fixes derivative issues for aniso tropic atoms in Al2O3
change names of sytsyms e.g. 2(100) to 2(x) for non-hex/trig space groups. Removes possible ambiguity with ortho, etc. cases.
Still checking extinction issues

  • 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-07 18:37:26 +0000 (Sun, 07 Sep 2014) $
4# $Author: vondreele $
5# $Revision: 1489 $
6# $URL: trunk/GSASIIstrIO.py $
7# $Id: GSASIIstrIO.py 1489 2014-09-07 18:37:26Z 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: 1489 $")
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                controlDict[pfx+'Tbar'] = extParms['Tbar']
1861                controlDict[pfx+'Cos2TM'] = extParms['Cos2TM']
1862                if 'Primary' in extType:
1863                    Ekey = ['Ep',]
1864                elif 'I & II' in extType:
1865                    Ekey = ['Eg','Es']
1866                elif 'Secondary Type II' == extType:
1867                    Ekey = ['Es',]
1868                elif 'Secondary Type I' == extType:
1869                    Ekey = ['Eg',]
1870                else:   #'None'
1871                    Ekey = []
1872                for eKey in Ekey:
1873                    hapDict[pfx+eKey] = extParms[eKey][0]
1874                    if extParms[eKey][1]:
1875                        hapVary.append(pfx+eKey)
1876                for bab in ['BabA','BabU']:
1877                    hapDict[pfx+bab] = hapData['Babinet'][bab][0]
1878                    if hapData['Babinet'][bab][1]:
1879                        hapVary.append(pfx+bab)
1880                if Print: 
1881                    print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
1882                    print >>pFile,135*'-'
1883                    print >>pFile,' Scale factor     : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1]
1884                    if extType != 'None':
1885                        print >>pFile,' Extinction  Type: %15s'%(extType),' approx: %10s'%(extApprox),' tbar: %6.3f'%(extParms['Tbar'])
1886                        text = ' Parameters       :'
1887                        for eKey in Ekey:
1888                            text += ' %4s : %10.3e Refine? '%(eKey,extParms[eKey][0])+str(extParms[eKey][1])
1889                        print >>pFile,text
1890                    if hapData['Babinet']['BabA'][0]:
1891                        PrintBabinet(hapData['Babinet'])
1892                Histogram['Reflection Lists'] = phase       
1893               
1894    return hapVary,hapDict,controlDict
1895   
1896def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True,pFile=None):
1897    'needs a doc string'
1898   
1899    def PrintSizeAndSig(hapData,sizeSig):
1900        line = '\n Size model:     %9s'%(hapData[0])
1901        refine = False
1902        if hapData[0] in ['isotropic','uniaxial']:
1903            line += ' equatorial:%12.5f'%(hapData[1][0])
1904            if sizeSig[0][0]:
1905                line += ', sig:%8.4f'%(sizeSig[0][0])
1906                refine = True
1907            if hapData[0] == 'uniaxial':
1908                line += ' axial:%12.4f'%(hapData[1][1])
1909                if sizeSig[0][1]:
1910                    refine = True
1911                    line += ', sig:%8.4f'%(sizeSig[0][1])
1912            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1913            if sizeSig[0][2]:
1914                refine = True
1915                line += ', sig:%8.4f'%(sizeSig[0][2])
1916            if refine:
1917                print >>pFile,line
1918        else:
1919            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1920            if sizeSig[0][2]:
1921                refine = True
1922                line += ', sig:%8.4f'%(sizeSig[0][2])
1923            Snames = ['S11','S22','S33','S12','S13','S23']
1924            ptlbls = ' name  :'
1925            ptstr =  ' value :'
1926            sigstr = ' sig   :'
1927            for i,name in enumerate(Snames):
1928                ptlbls += '%12s' % (name)
1929                ptstr += '%12.6f' % (hapData[4][i])
1930                if sizeSig[1][i]:
1931                    refine = True
1932                    sigstr += '%12.6f' % (sizeSig[1][i])
1933                else:
1934                    sigstr += 12*' '
1935            if refine:
1936                print >>pFile,line
1937                print >>pFile,ptlbls
1938                print >>pFile,ptstr
1939                print >>pFile,sigstr
1940       
1941    def PrintMuStrainAndSig(hapData,mustrainSig,SGData):
1942        line = '\n Mustrain model: %9s'%(hapData[0])
1943        refine = False
1944        if hapData[0] in ['isotropic','uniaxial']:
1945            line += ' equatorial:%12.1f'%(hapData[1][0])
1946            if mustrainSig[0][0]:
1947                line += ', sig:%8.1f'%(mustrainSig[0][0])
1948                refine = True
1949            if hapData[0] == 'uniaxial':
1950                line += ' axial:%12.1f'%(hapData[1][1])
1951                if mustrainSig[0][1]:
1952                     line += ', sig:%8.1f'%(mustrainSig[0][1])
1953            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1954            if mustrainSig[0][2]:
1955                refine = True
1956                line += ', sig:%8.3f'%(mustrainSig[0][2])
1957            if refine:
1958                print >>pFile,line
1959        else:
1960            line += ' LG mix coeff.:%12.4f'%(hapData[1][2])
1961            if mustrainSig[0][2]:
1962                refine = True
1963                line += ', sig:%8.3f'%(mustrainSig[0][2])
1964            Snames = G2spc.MustrainNames(SGData)
1965            ptlbls = ' name  :'
1966            ptstr =  ' value :'
1967            sigstr = ' sig   :'
1968            for i,name in enumerate(Snames):
1969                ptlbls += '%12s' % (name)
1970                ptstr += '%12.6f' % (hapData[4][i])
1971                if mustrainSig[1][i]:
1972                    refine = True
1973                    sigstr += '%12.6f' % (mustrainSig[1][i])
1974                else:
1975                    sigstr += 12*' '
1976            if refine:
1977                print >>pFile,line
1978                print >>pFile,ptlbls
1979                print >>pFile,ptstr
1980                print >>pFile,sigstr
1981           
1982    def PrintHStrainAndSig(hapData,strainSig,SGData):
1983        Hsnames = G2spc.HStrainNames(SGData)
1984        ptlbls = ' name  :'
1985        ptstr =  ' value :'
1986        sigstr = ' sig   :'
1987        refine = False
1988        for i,name in enumerate(Hsnames):
1989            ptlbls += '%12s' % (name)
1990            ptstr += '%12.4g' % (hapData[0][i])
1991            if name in strainSig:
1992                refine = True
1993                sigstr += '%12.4g' % (strainSig[name])
1994            else:
1995                sigstr += 12*' '
1996        if refine:
1997            print >>pFile,'\n Hydrostatic/elastic strain: '
1998            print >>pFile,ptlbls
1999            print >>pFile,ptstr
2000            print >>pFile,sigstr
2001       
2002    def PrintSHPOAndSig(pfx,hapData,POsig):
2003        print >>pFile,'\n Spherical harmonics preferred orientation: Order:'+str(hapData[4])
2004        ptlbls = ' names :'
2005        ptstr =  ' values:'
2006        sigstr = ' sig   :'
2007        for item in hapData[5]:
2008            ptlbls += '%12s'%(item)
2009            ptstr += '%12.3f'%(hapData[5][item])
2010            if pfx+item in POsig:
2011                sigstr += '%12.3f'%(POsig[pfx+item])
2012            else:
2013                sigstr += 12*' ' 
2014        print >>pFile,ptlbls
2015        print >>pFile,ptstr
2016        print >>pFile,sigstr
2017       
2018    def PrintExtAndSig(pfx,hapData,ScalExtSig):
2019        print >>pFile,'\n Single crystal extinction: Type: ',hapData[0],' Approx: ',hapData[1]
2020        text = ''
2021        for item in hapData[2]:
2022            if pfx+item in ScalExtSig:
2023                text += '       %s: '%(item)
2024                text += '%12.2e'%(hapData[2][item][0])
2025                if pfx+item in ScalExtSig:
2026                    text += ' sig: %12.2e'%(ScalExtSig[pfx+item])
2027        print >>pFile,text       
2028       
2029    def PrintBabinetAndSig(pfx,hapData,BabSig):
2030        print >>pFile,'\n Babinet form factor modification: '
2031        ptlbls = ' names :'
2032        ptstr =  ' values:'
2033        sigstr = ' sig   :'
2034        for item in hapData:
2035            ptlbls += '%12s'%(item)
2036            ptstr += '%12.3f'%(hapData[item][0])
2037            if pfx+item in BabSig:
2038                sigstr += '%12.3f'%(BabSig[pfx+item])
2039            else:
2040                sigstr += 12*' ' 
2041        print >>pFile,ptlbls
2042        print >>pFile,ptstr
2043        print >>pFile,sigstr
2044   
2045    PhFrExtPOSig = {}
2046    SizeMuStrSig = {}
2047    ScalExtSig = {}
2048    BabSig = {}
2049    wtFrSum = {}
2050    for phase in Phases:
2051        HistoPhase = Phases[phase]['Histograms']
2052        General = Phases[phase]['General']
2053        SGData = General['SGData']
2054        pId = Phases[phase]['pId']
2055        histoList = HistoPhase.keys()
2056        histoList.sort()
2057        for histogram in histoList:
2058            try:
2059                Histogram = Histograms[histogram]
2060            except KeyError:                       
2061                #skip if histogram not included e.g. in a sequential refinement
2062                continue
2063            hapData = HistoPhase[histogram]
2064            hId = Histogram['hId']
2065            pfx = str(pId)+':'+str(hId)+':'
2066            if hId not in wtFrSum:
2067                wtFrSum[hId] = 0.
2068            if 'PWDR' in histogram:
2069                for item in ['Scale','Extinction']:
2070                    hapData[item][0] = parmDict[pfx+item]
2071                    if pfx+item in sigDict:
2072                        PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2073                wtFrSum[hId] += hapData['Scale'][0]*General['Mass']
2074                if hapData['Pref.Ori.'][0] == 'MD':
2075                    hapData['Pref.Ori.'][1] = parmDict[pfx+'MD']
2076                    if pfx+'MD' in sigDict:
2077                        PhFrExtPOSig.update({pfx+'MD':sigDict[pfx+'MD'],})
2078                else:                           #'SH' spherical harmonics
2079                    for item in hapData['Pref.Ori.'][5]:
2080                        hapData['Pref.Ori.'][5][item] = parmDict[pfx+item]
2081                        if pfx+item in sigDict:
2082                            PhFrExtPOSig.update({pfx+item:sigDict[pfx+item],})
2083                SizeMuStrSig.update({pfx+'Mustrain':[[0,0,0],[0 for i in range(len(hapData['Mustrain'][4]))]],
2084                    pfx+'Size':[[0,0,0],[0 for i in range(len(hapData['Size'][4]))]],
2085                    pfx+'HStrain':{}})                 
2086                for item in ['Mustrain','Size']:
2087                    hapData[item][1][2] = parmDict[pfx+item+';mx']
2088                    hapData[item][1][2] = min(1.,max(0.,hapData[item][1][2]))
2089                    if pfx+item+';mx' in sigDict:
2090                        SizeMuStrSig[pfx+item][0][2] = sigDict[pfx+item+';mx']
2091                    if hapData[item][0] in ['isotropic','uniaxial']:                   
2092                        hapData[item][1][0] = parmDict[pfx+item+';i']
2093                        if item == 'Size':
2094                            hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0]))
2095                        if pfx+item+';i' in sigDict: 
2096                            SizeMuStrSig[pfx+item][0][0] = sigDict[pfx+item+';i']
2097                        if hapData[item][0] == 'uniaxial':
2098                            hapData[item][1][1] = parmDict[pfx+item+';a']
2099                            if item == 'Size':
2100                                hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1]))                       
2101                            if pfx+item+';a' in sigDict:
2102                                SizeMuStrSig[pfx+item][0][1] = sigDict[pfx+item+';a']
2103                    else:       #generalized for mustrain or ellipsoidal for size
2104                        Nterms = len(hapData[item][4])
2105                        for i in range(Nterms):
2106                            sfx = ':'+str(i)
2107                            hapData[item][4][i] = parmDict[pfx+item+sfx]
2108                            if pfx+item+sfx in sigDict:
2109                                SizeMuStrSig[pfx+item][1][i] = sigDict[pfx+item+sfx]
2110                names = G2spc.HStrainNames(SGData)
2111                for i,name in enumerate(names):
2112                    hapData['HStrain'][0][i] = parmDict[pfx+name]
2113                    if pfx+name in sigDict:
2114                        SizeMuStrSig[pfx+'HStrain'][name] = sigDict[pfx+name]
2115                for name in ['BabA','BabU']:
2116                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2117                    if pfx+name in sigDict:
2118                        BabSig[pfx+name] = sigDict[pfx+name]               
2119               
2120            elif 'HKLF' in histogram:
2121                for item in ['Scale',]:
2122                    if parmDict.get(pfx+item):
2123                        hapData[item][0] = parmDict[pfx+item]
2124                        if pfx+item in sigDict:
2125                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2126                for item in ['Ep','Eg','Es']:
2127                    if parmDict.get(pfx+item):
2128                        hapData['Extinction'][2][item][0] = parmDict[pfx+item]
2129                        if pfx+item in sigDict:
2130                            ScalExtSig[pfx+item] = sigDict[pfx+item]
2131                for name in ['BabA','BabU']:
2132                    hapData['Babinet'][name][0] = parmDict[pfx+name]
2133                    if pfx+name in sigDict:
2134                        BabSig[pfx+name] = sigDict[pfx+name]               
2135
2136    if Print:
2137        for phase in Phases:
2138            HistoPhase = Phases[phase]['Histograms']
2139            General = Phases[phase]['General']
2140            SGData = General['SGData']
2141            pId = Phases[phase]['pId']
2142            histoList = HistoPhase.keys()
2143            histoList.sort()
2144            for histogram in histoList:
2145                try:
2146                    Histogram = Histograms[histogram]
2147                except KeyError:                       
2148                    #skip if histogram not included e.g. in a sequential refinement
2149                    continue
2150                print >>pFile,'\n Phase: ',phase,' in histogram: ',histogram
2151                print >>pFile,130*'-'
2152                hapData = HistoPhase[histogram]
2153                hId = Histogram['hId']
2154                Histogram['Residuals'][str(pId)+'::Name'] = phase
2155                pfx = str(pId)+':'+str(hId)+':'
2156                if 'PWDR' in histogram:
2157                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2158                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2159               
2160                    if pfx+'Scale' in PhFrExtPOSig:
2161                        wtFr = hapData['Scale'][0]*General['Mass']/wtFrSum[hId]
2162                        sigwtFr = PhFrExtPOSig[pfx+'Scale']*wtFr/hapData['Scale'][0]
2163                        print >>pFile,' Phase fraction  : %10.5f, sig %10.5f Weight fraction  : %8.5f, sig %10.5f' \
2164                            %(hapData['Scale'][0],PhFrExtPOSig[pfx+'Scale'],wtFr,sigwtFr)
2165                    if pfx+'Extinction' in PhFrExtPOSig:
2166                        print >>pFile,' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig[pfx+'Extinction'])
2167                    if hapData['Pref.Ori.'][0] == 'MD':
2168                        if pfx+'MD' in PhFrExtPOSig:
2169                            print >>pFile,' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig[pfx+'MD'])
2170                    else:
2171                        PrintSHPOAndSig(pfx,hapData['Pref.Ori.'],PhFrExtPOSig)
2172                    PrintSizeAndSig(hapData['Size'],SizeMuStrSig[pfx+'Size'])
2173                    PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig[pfx+'Mustrain'],SGData)
2174                    PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig[pfx+'HStrain'],SGData)
2175                    if len(BabSig):
2176                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2177                   
2178                elif 'HKLF' in histogram:
2179                    print >>pFile,' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections'   \
2180                        %(Histogram['Residuals'][pfx+'Rf'],Histogram['Residuals'][pfx+'Rf^2'],Histogram['Residuals'][pfx+'Nref'])
2181                    print >>pFile,' HKLF histogram weight factor = ','%.3f'%(Histogram['wtFactor'])
2182                    if pfx+'Scale' in ScalExtSig:
2183                        print >>pFile,' Scale factor : %10.4f, sig %10.4f'%(hapData['Scale'][0],ScalExtSig[pfx+'Scale'])
2184                    if hapData['Extinction'][0] != 'None':
2185                        PrintExtAndSig(pfx,hapData['Extinction'],ScalExtSig)
2186                    if len(BabSig):
2187                        PrintBabinetAndSig(pfx,hapData['Babinet'],BabSig)
2188
2189################################################################################
2190##### Histogram data
2191################################################################################       
2192                   
2193def GetHistogramData(Histograms,Print=True,pFile=None):
2194    'needs a doc string'
2195   
2196    def GetBackgroundParms(hId,Background):
2197        Back = Background[0]
2198        DebyePeaks = Background[1]
2199        bakType,bakFlag = Back[:2]
2200        backVals = Back[3:]
2201        backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))]
2202        backDict = dict(zip(backNames,backVals))
2203        backVary = []
2204        if bakFlag:
2205            backVary = backNames
2206        backDict[':'+str(hId)+':nDebye'] = DebyePeaks['nDebye']
2207        backDict[':'+str(hId)+':nPeaks'] = DebyePeaks['nPeaks']
2208        debyeDict = {}
2209        debyeList = []
2210        for i in range(DebyePeaks['nDebye']):
2211            debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)]
2212            debyeDict.update(dict(zip(debyeNames,DebyePeaks['debyeTerms'][i][::2])))
2213            debyeList += zip(debyeNames,DebyePeaks['debyeTerms'][i][1::2])
2214        debyeVary = []
2215        for item in debyeList:
2216            if item[1]:
2217                debyeVary.append(item[0])
2218        backDict.update(debyeDict)
2219        backVary += debyeVary
2220        peakDict = {}
2221        peakList = []
2222        for i in range(DebyePeaks['nPeaks']):
2223            peakNames = [':'+str(hId)+':BkPkpos;'+str(i),':'+str(hId)+ \
2224                ':BkPkint;'+str(i),':'+str(hId)+':BkPksig;'+str(i),':'+str(hId)+':BkPkgam;'+str(i)]
2225            peakDict.update(dict(zip(peakNames,DebyePeaks['peaksList'][i][::2])))
2226            peakList += zip(peakNames,DebyePeaks['peaksList'][i][1::2])
2227        peakVary = []
2228        for item in peakList:
2229            if item[1]:
2230                peakVary.append(item[0])
2231        backDict.update(peakDict)
2232        backVary += peakVary
2233        return bakType,backDict,backVary           
2234       
2235    def GetInstParms(hId,Inst):     
2236        dataType = Inst['Type'][0]
2237        instDict = {}
2238        insVary = []
2239        pfx = ':'+str(hId)+':'
2240        insKeys = Inst.keys()
2241        insKeys.sort()
2242        for item in insKeys:
2243            insName = pfx+item
2244            instDict[insName] = Inst[item][1]
2245            if len(Inst[item]) > 2 and Inst[item][2]:
2246                insVary.append(insName)
2247        if 'C' in dataType:
2248            instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005)
2249        return dataType,instDict,insVary
2250       
2251    def GetSampleParms(hId,Sample):
2252        sampVary = []
2253        hfx = ':'+str(hId)+':'       
2254        sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'],
2255            hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']}
2256        for key in ('Temperature','Pressure','FreePrm1','FreePrm2','FreePrm3'):
2257            if key in Sample:
2258                sampDict[hfx+key] = Sample[key]
2259        Type = Sample['Type']
2260        if 'Bragg' in Type:             #Bragg-Brentano
2261            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2262                sampDict[hfx+item] = Sample[item][0]
2263                if Sample[item][1]:
2264                    sampVary.append(hfx+item)
2265        elif 'Debye' in Type:        #Debye-Scherrer
2266            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2267                sampDict[hfx+item] = Sample[item][0]
2268                if Sample[item][1]:
2269                    sampVary.append(hfx+item)
2270        return Type,sampDict,sampVary
2271       
2272    def PrintBackground(Background):
2273        Back = Background[0]
2274        DebyePeaks = Background[1]
2275        print >>pFile,'\n Background function: ',Back[0],' Refine?',bool(Back[1])
2276        line = ' Coefficients: '
2277        for i,back in enumerate(Back[3:]):
2278            line += '%10.3f'%(back)
2279            if i and not i%10:
2280                line += '\n'+15*' '
2281        print >>pFile,line
2282        if DebyePeaks['nDebye']:
2283            print >>pFile,'\n Debye diffuse scattering coefficients'
2284            parms = ['DebyeA','DebyeR','DebyeU']
2285            line = ' names :  '
2286            for parm in parms:
2287                line += '%8s refine?'%(parm)
2288            print >>pFile,line
2289            for j,term in enumerate(DebyePeaks['debyeTerms']):
2290                line = ' term'+'%2d'%(j)+':'
2291                for i in range(3):
2292                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2293                print >>pFile,line
2294        if DebyePeaks['nPeaks']:
2295            print >>pFile,'\n Single peak coefficients'
2296            parms =    ['BkPkpos','BkPkint','BkPksig','BkPkgam']
2297            line = ' names :  '
2298            for parm in parms:
2299                line += '%8s refine?'%(parm)
2300            print >>pFile,line
2301            for j,term in enumerate(DebyePeaks['peaksList']):
2302                line = ' peak'+'%2d'%(j)+':'
2303                for i in range(4):
2304                    line += '%10.3f %5s'%(term[2*i],bool(term[2*i+1]))                   
2305                print >>pFile,line
2306       
2307    def PrintInstParms(Inst):
2308        print >>pFile,'\n Instrument Parameters:'
2309        insKeys = Inst.keys()
2310        insKeys.sort()
2311        iBeg = 0
2312        Ok = True
2313        while Ok:
2314            ptlbls = ' name  :'
2315            ptstr =  ' value :'
2316            varstr = ' refine:'
2317            iFin = min(iBeg+9,len(insKeys))
2318            for item in insKeys[iBeg:iFin]:
2319                if item not in ['Type','Source']:
2320                    ptlbls += '%12s' % (item)
2321                    ptstr += '%12.6f' % (Inst[item][1])
2322                    if item in ['Lam1','Lam2','Azimuth','fltPath','2-theta',]:
2323                        varstr += 12*' '
2324                    else:
2325                        varstr += '%12s' % (str(bool(Inst[item][2])))
2326            print >>pFile,ptlbls
2327            print >>pFile,ptstr
2328            print >>pFile,varstr
2329            iBeg = iFin
2330            if iBeg == len(insKeys):
2331                Ok = False
2332            else:
2333                print >>pFile,'\n'
2334       
2335    def PrintSampleParms(Sample):
2336        print >>pFile,'\n Sample Parameters:'
2337        print >>pFile,' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \
2338            (Sample['Omega'],Sample['Chi'],Sample['Phi'])
2339        ptlbls = ' name  :'
2340        ptstr =  ' value :'
2341        varstr = ' refine:'
2342        if 'Bragg' in Sample['Type']:
2343            for item in ['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']:
2344                ptlbls += '%14s'%(item)
2345                ptstr += '%14.4f'%(Sample[item][0])
2346                varstr += '%14s'%(str(bool(Sample[item][1])))
2347           
2348        elif 'Debye' in Type:        #Debye-Scherrer
2349            for item in ['Scale','Absorption','DisplaceX','DisplaceY']:
2350                ptlbls += '%14s'%(item)
2351                ptstr += '%14.4f'%(Sample[item][0])
2352                varstr += '%14s'%(str(bool(Sample[item][1])))
2353
2354        print >>pFile,ptlbls
2355        print >>pFile,ptstr
2356        print >>pFile,varstr
2357       
2358    histDict = {}
2359    histVary = []
2360    controlDict = {}
2361    histoList = Histograms.keys()
2362    histoList.sort()
2363    for histogram in histoList:
2364        if 'PWDR' in histogram:
2365            Histogram = Histograms[histogram]
2366            hId = Histogram['hId']
2367            pfx = ':'+str(hId)+':'
2368            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2369            controlDict[pfx+'Limits'] = Histogram['Limits'][1]
2370            controlDict[pfx+'Exclude'] = Histogram['Limits'][2:]
2371            for excl in controlDict[pfx+'Exclude']:
2372                Histogram['Data'][0] = ma.masked_inside(Histogram['Data'][0],excl[0],excl[1])
2373            if controlDict[pfx+'Exclude']:
2374                ma.mask_rows(Histogram['Data'])
2375            Background = Histogram['Background']
2376            Type,bakDict,bakVary = GetBackgroundParms(hId,Background)
2377            controlDict[pfx+'bakType'] = Type
2378            histDict.update(bakDict)
2379            histVary += bakVary
2380           
2381            Inst = Histogram['Instrument Parameters'][0]
2382            Type,instDict,insVary = GetInstParms(hId,Inst)
2383            controlDict[pfx+'histType'] = Type
2384            if 'XC' in Type:
2385                if pfx+'Lam1' in instDict:
2386                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1']
2387                else:
2388                    controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam']           
2389            histDict.update(instDict)
2390            histVary += insVary
2391           
2392            Sample = Histogram['Sample Parameters']
2393            Type,sampDict,sampVary = GetSampleParms(hId,Sample)
2394            controlDict[pfx+'instType'] = Type
2395            histDict.update(sampDict)
2396            histVary += sampVary
2397           
2398   
2399            if Print: 
2400                print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2401                print >>pFile,135*'-'
2402                Units = {'C':' deg','T':' msec'}
2403                units = Units[controlDict[pfx+'histType'][2]]
2404                Limits = controlDict[pfx+'Limits']
2405                print >>pFile,' Instrument type: ',Sample['Type']
2406                print >>pFile,' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units)
2407                if len(controlDict[pfx+'Exclude']):
2408                    excls = controlDict[pfx+'Exclude']
2409                    for excl in excls:
2410                        print >>pFile,' Excluded region:  %8.2f%s to %8.2f%s'%(excl[0],units,excl[1],units)   
2411                PrintSampleParms(Sample)
2412                PrintInstParms(Inst)
2413                PrintBackground(Background)
2414        elif 'HKLF' in histogram:
2415            Histogram = Histograms[histogram]
2416            hId = Histogram['hId']
2417            pfx = ':'+str(hId)+':'
2418            controlDict[pfx+'wtFactor'] = Histogram['wtFactor']
2419            Inst = Histogram['Instrument Parameters'][0]
2420            controlDict[pfx+'histType'] = Inst['Type'][0]
2421            if 'X' in Inst['Type'][0]:
2422                histDict[pfx+'Lam'] = Inst['Lam'][1]
2423                controlDict[pfx+'keV'] = 12.397639/histDict[pfx+'Lam']                   
2424    return histVary,histDict,controlDict
2425   
2426def SetHistogramData(parmDict,sigDict,Histograms,Print=True,pFile=None):
2427    'needs a doc string'
2428   
2429    def SetBackgroundParms(pfx,Background,parmDict,sigDict):
2430        Back = Background[0]
2431        DebyePeaks = Background[1]
2432        lenBack = len(Back[3:])
2433        backSig = [0 for i in range(lenBack+3*DebyePeaks['nDebye']+4*DebyePeaks['nPeaks'])]
2434        for i in range(lenBack):
2435            Back[3+i] = parmDict[pfx+'Back:'+str(i)]
2436            if pfx+'Back:'+str(i) in sigDict:
2437                backSig[i] = sigDict[pfx+'Back:'+str(i)]
2438        if DebyePeaks['nDebye']:
2439            for i in range(DebyePeaks['nDebye']):
2440                names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)]
2441                for j,name in enumerate(names):
2442                    DebyePeaks['debyeTerms'][i][2*j] = parmDict[name]
2443                    if name in sigDict:
2444                        backSig[lenBack+3*i+j] = sigDict[name]           
2445        if DebyePeaks['nPeaks']:
2446            for i in range(DebyePeaks['nPeaks']):
2447                names = [pfx+'BkPkpos;'+str(i),pfx+'BkPkint;'+str(i),
2448                    pfx+'BkPksig;'+str(i),pfx+'BkPkgam;'+str(i)]
2449                for j,name in enumerate(names):
2450                    DebyePeaks['peaksList'][i][2*j] = parmDict[name]
2451                    if name in sigDict:
2452                        backSig[lenBack+3*DebyePeaks['nDebye']+4*i+j] = sigDict[name]
2453        return backSig
2454       
2455    def SetInstParms(pfx,Inst,parmDict,sigDict):
2456        instSig = {}
2457        insKeys = Inst.keys()
2458        insKeys.sort()
2459        for item in insKeys:
2460            insName = pfx+item
2461            Inst[item][1] = parmDict[insName]
2462            if insName in sigDict:
2463                instSig[item] = sigDict[insName]
2464            else:
2465                instSig[item] = 0
2466        return instSig
2467       
2468    def SetSampleParms(pfx,Sample,parmDict,sigDict):
2469        if 'Bragg' in Sample['Type']:             #Bragg-Brentano
2470            sampSig = [0 for i in range(5)]
2471            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2472                Sample[item][0] = parmDict[pfx+item]
2473                if pfx+item in sigDict:
2474                    sampSig[i] = sigDict[pfx+item]
2475        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2476            sampSig = [0 for i in range(4)]
2477            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2478                Sample[item][0] = parmDict[pfx+item]
2479                if pfx+item in sigDict:
2480                    sampSig[i] = sigDict[pfx+item]
2481        return sampSig
2482       
2483    def PrintBackgroundSig(Background,backSig):
2484        Back = Background[0]
2485        DebyePeaks = Background[1]
2486        lenBack = len(Back[3:])
2487        valstr = ' value : '
2488        sigstr = ' sig   : '
2489        refine = False
2490        for i,back in enumerate(Back[3:]):
2491            valstr += '%10.4g'%(back)
2492            if Back[1]:
2493                refine = True
2494                sigstr += '%10.4g'%(backSig[i])
2495            else:
2496                sigstr += 10*' '
2497        if refine:
2498            print >>pFile,'\n Background function: ',Back[0]
2499            print >>pFile,valstr
2500            print >>pFile,sigstr
2501        if DebyePeaks['nDebye']:
2502            ifAny = False
2503            ptfmt = "%12.3f"
2504            names =  ' names :'
2505            ptstr =  ' values:'
2506            sigstr = ' esds  :'
2507            for item in sigDict:
2508                if 'Debye' in item:
2509                    ifAny = True
2510                    names += '%12s'%(item)
2511                    ptstr += ptfmt%(parmDict[item])
2512                    sigstr += ptfmt%(sigDict[item])
2513            if ifAny:
2514                print >>pFile,'\n Debye diffuse scattering coefficients'
2515                print >>pFile,names
2516                print >>pFile,ptstr
2517                print >>pFile,sigstr
2518        if DebyePeaks['nPeaks']:
2519            ifAny = False
2520            ptfmt = "%14.3f"
2521            names =  ' names :'
2522            ptstr =  ' values:'
2523            sigstr = ' esds  :'
2524            for item in sigDict:
2525                if 'BkPk' in item:
2526                    ifAny = True
2527                    names += '%14s'%(item)
2528                    ptstr += ptfmt%(parmDict[item])
2529                    sigstr += ptfmt%(sigDict[item])
2530            if ifAny:
2531                print >>pFile,'\n Single peak coefficients'
2532                print >>pFile,names
2533                print >>pFile,ptstr
2534                print >>pFile,sigstr
2535       
2536    def PrintInstParmsSig(Inst,instSig):
2537        refine = False
2538        insKeys = instSig.keys()
2539        insKeys.sort()
2540        iBeg = 0
2541        Ok = True
2542        while Ok:
2543            ptlbls = ' names :'
2544            ptstr =  ' value :'
2545            sigstr = ' sig   :'
2546            iFin = min(iBeg+9,len(insKeys))
2547            for name in insKeys[iBeg:iFin]:
2548                if name not in  ['Type','Lam1','Lam2','Azimuth','Source','fltPath']:
2549                    ptlbls += '%12s' % (name)
2550                    ptstr += '%12.6f' % (Inst[name][1])
2551                    if instSig[name]:
2552                        refine = True
2553                        sigstr += '%12.6f' % (instSig[name])
2554                    else:
2555                        sigstr += 12*' '
2556            if refine:
2557                print >>pFile,'\n Instrument Parameters:'
2558                print >>pFile,ptlbls
2559                print >>pFile,ptstr
2560                print >>pFile,sigstr
2561            iBeg = iFin
2562            if iBeg == len(insKeys):
2563                Ok = False
2564       
2565    def PrintSampleParmsSig(Sample,sampleSig):
2566        ptlbls = ' names :'
2567        ptstr =  ' values:'
2568        sigstr = ' sig   :'
2569        refine = False
2570        if 'Bragg' in Sample['Type']:
2571            for i,item in enumerate(['Scale','Shift','Transparency','SurfRoughA','SurfRoughB']):
2572                ptlbls += '%14s'%(item)
2573                ptstr += '%14.4f'%(Sample[item][0])
2574                if sampleSig[i]:
2575                    refine = True
2576                    sigstr += '%14.4f'%(sampleSig[i])
2577                else:
2578                    sigstr += 14*' '
2579           
2580        elif 'Debye' in Sample['Type']:        #Debye-Scherrer
2581            for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']):
2582                ptlbls += '%14s'%(item)
2583                ptstr += '%14.4f'%(Sample[item][0])
2584                if sampleSig[i]:
2585                    refine = True
2586                    sigstr += '%14.4f'%(sampleSig[i])
2587                else:
2588                    sigstr += 14*' '
2589
2590        if refine:
2591            print >>pFile,'\n Sample Parameters:'
2592            print >>pFile,ptlbls
2593            print >>pFile,ptstr
2594            print >>pFile,sigstr
2595       
2596    histoList = Histograms.keys()
2597    histoList.sort()
2598    for histogram in histoList:
2599        if 'PWDR' in histogram:
2600            Histogram = Histograms[histogram]
2601            hId = Histogram['hId']
2602            pfx = ':'+str(hId)+':'
2603            Background = Histogram['Background']
2604            backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict)
2605           
2606            Inst = Histogram['Instrument Parameters'][0]
2607            instSig = SetInstParms(pfx,Inst,parmDict,sigDict)
2608       
2609            Sample = Histogram['Sample Parameters']
2610            sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict)
2611
2612            print >>pFile,'\n Histogram: ',histogram,' histogram Id: ',hId
2613            print >>pFile,135*'-'
2614            print >>pFile,' PWDR histogram weight factor = '+'%.3f'%(Histogram['wtFactor'])
2615            print >>pFile,' Final refinement wR = %.2f%% on %d observations in this histogram'%(Histogram['Residuals']['wR'],Histogram['Residuals']['Nobs'])
2616            print >>pFile,' Other residuals: R = %.2f%%, Rb = %.2f%%, wRb = %.2f%% wRmin = %.2f%%'% \
2617                (Histogram['Residuals']['R'],Histogram['Residuals']['Rb'],Histogram['Residuals']['wRb'],Histogram['Residuals']['wRmin'])
2618            if Print:
2619                print >>pFile,' Instrument type: ',Sample['Type']
2620                PrintSampleParmsSig(Sample,sampSig)
2621                PrintInstParmsSig(Inst,instSig)
2622                PrintBackgroundSig(Background,backSig)
2623               
Note: See TracBrowser for help on using the repository browser.