source: trunk/GSASIIstrIO.py @ 1570

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

fix to G2strIO for space group output.
finish interpretation of super symmetry symbols

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