source: trunk/GSASIIstrIO.py @ 1604

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

make modulation wave maps
fix atom index bugs
begin modulated structure imput to least squares

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