source: trunk/GSASIIstrIO.py @ 1529

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

revise display of space group info so that operators are in lined up columns
in bot the LS output & the popup window. - new class SGMessageBox (from wx.Dialog!)
G2spc.SGPrint now returns 2 items: Text & Tables
Work on supersymmetry - now better display of operators
operators more complete
use ErrorDialog? for more error messages instead of putting the message on the console

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