source: trunk/GSASIIstrIO.py @ 934

Last change on this file since 934 was 934, checked in by vondreele, 9 years ago

changes for MC/SA

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