source: trunk/GSASIIstrIO.py @ 1596

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

some bug fixes

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