source: trunk/GSASIIstrIO.py @ 1594

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

begin Pawley refinement stuff for incommensurate structures
Pawley SS reflections done
start on Refine for SS
only allow (3+1) supersymmetry - remove all (3+2) & (3+3) stuff
fix more bugs...

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