source: trunk/GSASIIstruct.py @ 230

Last change on this file since 230 was 230, checked in by vondreele, 11 years ago

add svn keywords & add log intensity plotting

  • Property svn:keywords set to Date Author Revision URL Id
File size: 11.4 KB
Line 
1#GSASIIstructure - structure computation routines
2########### SVN repository information ###################
3# $Date: 2011-01-07 19:27:24 +0000 (Fri, 07 Jan 2011) $
4# $Author: vondreele $
5# $Revision: 230 $
6# $URL: trunk/GSASIIstruct.py $
7# $Id: GSASIIstruct.py 230 2011-01-07 19:27:24Z vondreele $
8########### SVN repository information ###################
9import sys
10import os.path as ospath
11import numpy as np
12import cPickle
13import time
14import math
15import GSASIIpath
16import GSASIIElem as G2el
17import GSASIIlattice as G2lat
18import GSASIIspc as G2spc
19
20def ShowBanner():
21    print 80*'*'
22    print '   General Structure Analysis System-II Crystal Structure Refinement'
23    print '     by Robert B. Von Dreele, Argonne National Laboratory(C), 2010'
24    print ' This product includes software developed by the UChicago Argonne, LLC,' 
25    print '            as Operator of Argonne National Laboratory.'
26    print 80*'*','\n'
27
28def GetControls(GPXfile):
29    ''' Returns dictionary of control items found in GSASII gpx file
30    input:
31        GPXfile = .gpx full file name
32    return:
33        Controls = dictionary of control items
34    '''
35    file = open(GPXfile,'rb')
36    HistogramNames = []
37    while True:
38        try:
39            data = cPickle.load(file)
40        except EOFError:
41            break
42        datum = data[0]
43        if datum[0] == 'Controls':
44            Controls = datum[1]
45    file.close()
46    return Controls
47   
48def ShowControls(Controls):
49    print ' Controls:'
50    if Controls['bandWidth']:
51        print ' Band approximated least squares matrix refinement, width: ',Controls['bandWidth']
52    else:
53        print ' Full matrix least squares refinement'
54    print ' Maximum number of refinement cycles: ',Controls['Ncycles']
55    print ' Minimum sum(shift/esd)^2 for convergence: ','%.2f'%(Controls['minSumShftESD'])
56    print ' The Marquardt damping factor: ','%.2f'%(Controls['Marquardt'])
57    print ' Maximum allowed atom shift: ','%.2f'%(Controls['maxShift']),'A'
58    if Controls['restraintWeight']:
59        print ' The weights of the restraint histograms will be modified to normalize their contribution','\n'
60    else:
61        print ' The restraint histogram weights are fixed','\n'
62   
63def GetPhaseNames(GPXfile):
64    ''' Returns a list of phase names found under 'Phases' in GSASII gpx file
65    input:
66        GPXfile = gpx full file name
67    return:
68        PhaseNames = list of phase names
69    '''
70    file = open(GPXfile,'rb')
71    PhaseNames = []
72    while True:
73        try:
74            data = cPickle.load(file)
75        except EOFError:
76            break
77        datum = data[0]
78        if 'Phases' == datum[0]:
79            for datus in data[1:]:
80                PhaseNames.append(datus[0])
81    file.close()
82    return PhaseNames
83
84def GetPhaseData(GPXfile,PhaseName):
85    ''' Returns the 'General' and 'Atoms' objects for PhaseName from GSASII gpx file
86    input:
87        GPXfile = gpx full file name
88        PhaseName = phase name
89    return:
90        General = dictionary of general phase info.
91        Atoms = list of atom parameters
92        these are returned empty if PhaseName not found
93    '''       
94    file = open(GPXfile,'rb')
95    General = {}
96    Atoms = []
97    while True:
98        try:
99            data = cPickle.load(file)
100        except EOFError:
101            break
102        datum = data[0]
103        if 'Phases' == datum[0]:
104            for datus in data[1:]:
105                if datus[0] == PhaseName:
106                    General = datus[1]['General']
107                    Atoms = datus[1]['Atoms']                                           
108    file.close()
109    return {'General':General,'Atoms':Atoms}
110   
111def ShowPhaseData(phaseNames,PhaseData):
112    print ' Phases:'
113    for name in phaseNames:
114        General = PhaseData[name]['General']
115        Atoms = PhaseData[name]['Atoms']
116        print '\n Phase name: ',General['Name']
117        SGtext = G2spc.SGPrint(General['SGData'])
118        for line in SGtext: print line
119        cell = General['Cell']
120        print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \
121            ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =','%.3f'%(cell[6]),' volume =','%.3f'%(cell[7])
122        print ' Refine?',cell[0]
123        print '\n Atoms:'
124        line = '   name    type  refine?   x         y         z    '+ \
125            '  frac site sym  mult I/A   Uiso     U11     U22     U33     U12     U13     U23'
126        if General['Type'] == 'magnetic':
127            line += '   Mx     My     Mz'
128        elif General['Type'] == 'macromolecular':
129            line = ' res no  residue  chain '+line
130        print line
131        if General['Type'] == 'nuclear':
132            print 135*'-'
133            for at in Atoms:
134                line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \
135                    '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9])
136                if at[9] == 'I':
137                    line += '%8.4f'%(at[10])+48*' '
138                else:
139                    line += 8*' '
140                    for i in range(6):
141                        line += '%8.4f'%(at[11+i])
142                print line
143#        elif General['Type'] == 'magnetic':
144#        elif General['Type'] == 'macromolecular':
145       
146
147def GetHistogramNames(GPXfile):
148    ''' Returns a list of histogram names found in GSASII gpx file
149    input:
150        GPXfile = .gpx full file name
151    return:
152        HistogramNames = list of histogram names (types = PWDR & HKLF)
153    '''
154    file = open(GPXfile,'rb')
155    HistogramNames = []
156    while True:
157        try:
158            data = cPickle.load(file)
159        except EOFError:
160            break
161        datum = data[0]
162        if datum[0][:4] in ['PWDR','HKLF']:
163            HistogramNames.append(datum[0])
164    file.close()
165    return HistogramNames
166   
167def GetPWDRdata(GPXfile,PWDRname):
168    ''' Returns powder data from GSASII gpx file
169    input:
170        GPXfile = .gpx full file name
171        PWDRname = powder histogram name as obtained from GetHistogramNames
172    return:
173        PWDRdata = powder data list:
174       
175    '''
176    file = open(GPXfile,'rb')
177    PWDRdata = []
178    while True:
179        try:
180            data = cPickle.load(file)
181        except EOFError:
182            break
183        datum = data[0]
184        if datum[0] == PWDRname:
185            PWDRdata = datum[1:][0]
186                                               
187           
188    file.close()
189    return PWDRdata#,Limits,InstrumentParms
190   
191def GetHKLFdata(GPXfile,HKLFname):
192    ''' Returns single crystal data from GSASII gpx file
193    input:
194        GPXfile = .gpx full file name
195        HKLFname = single crystal histogram name as obtained from GetHistogramNames
196    return:
197        HKLFdata = single crystal data list of reflections: for each reflection:
198            HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase]
199    '''
200    file = open(GPXfile,'rb')
201    HKLFdata = []
202    while True:
203        try:
204            data = cPickle.load(file)
205        except EOFError:
206            break
207        datum = data[0]
208        if datum[0] == HKLFname:
209            HKLFdata = datum[1:][0]
210    file.close()
211    return HKLFdata
212   
213def GetFFtable(General):
214    ''' returns a dictionary of form factor data for atom types found in General
215    input:
216        General = dictionary of phase info.; includes AtomTypes
217    return:
218        FFtable = dictionary of form factor data; key is atom type
219    '''
220    atomTypes = General['AtomTypes']
221    FFtable = {}
222    for El in atomTypes:
223        FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0])
224        for item in FFs:
225            if item['Symbol'] == El:
226                FFtable[El] = item
227    return FFtable
228       
229def SetupSFcalc(General,Atoms):
230    ''' setup data for use in StructureFactor; mostly rearranging arrays
231    input:
232        General = dictionary of general phase info.
233        Atoms = list of atom parameters
234    returns:
235        G = reciprocal metric tensor
236        StrData = list of arrays; one entry per atom:
237            T = atom types
238            R = refinement flags, e.g. 'FXU'
239            F = site fractions
240            X = atom coordinates as numpy array
241            M = atom multiplicities
242            IA = isotropic/anisotropic thermal flags
243            Uiso = isotropic thermal parameters if IA = 'I'; else = 0
244            Uij = numpy array of 6 anisotropic thermal parameters if IA='A'; else = zeros
245    '''           
246    SGData = General['SGData']
247    Cell = General['Cell']
248    G,g = G2lat.cell2Gmat(Cell[1:7])        #skip refine & volume; get recip & real metric tensors
249    cx,ct,cs,cia = General['AtomPtrs']
250    X = [];F = [];T = [];IA = [];Uiso = [];Uij = [];R = [];M = []
251    for atom in Atoms:
252        T.append(atom[ct])
253        R.append(atom[ct+1])
254        F.append(atom[cx+3])
255        X.append(np.array(atom[cx:cx+3]))
256        M.append(atom[cia-1])
257        IA.append(atom[cia])
258        Uiso.append(atom[cia+1])
259        Uij.append(np.array(atom[cia+2:cia+8]))
260    return G,[T,R,np.array(F),np.array(X),np.array(M),IA,np.array(Uiso),np.array(Uij)]
261           
262def StructureFactor(H,G,SGData,StrData,FFtable):
263    ''' Compute structure factor for a single h,k,l
264    H = np.array(h,k,l)
265    G = reciprocal metric tensor
266    SGData = space group info. dictionary output from SpcGroup
267    StrData = [
268        [atom types],
269        [refinement flags],
270        [site fractions],
271        np.array(coordinates),
272        [multiplicities],
273        [I/A flag],
274        [Uiso],
275        np.array(Uij)]
276    FFtable = dictionary of form factor coeff. for atom types used in StrData
277    '''
278    twopi = 2.0*math.pi
279    twopisq = 2.0*math.pi**2
280    SQ = G2lat.calc_rDsq2(H,G)/4.0          # SQ = (sin(theta)/lambda)**2
281    SQfactor = 8.0*SQ*math.pi**2
282    print 'SQ',SQfactor
283    FF = {}
284    for type in FFtable.keys():
285        FF[type] = G2el.ScatFac(FFtable[type],SQ)
286    print 'FF',FF
287    iabsnt,mulp,Uniq,Phs = G2spc.GenHKL(H,SGData,False)
288    fa = [0,0,0]        #real
289    fb = [0,0,0]        #imaginary
290    if not iabsnt:
291        phase = twopi*np.inner(Uniq,StrData[3])       #2pi[hx+ky+lz] for each atom in each equiv. hkl
292        sinp = np.sin(phase)
293        cosp = np.cos(phase)
294        occ = StrData[2]*StrData[4]/mulp
295        Tiso = np.multiply(StrData[6],SQfactor)
296        Tuij = np.multiply(-SQfactor,np.inner(H,np.inner(G2spc.Uij2U(StrData[7]),H)))
297        print 'sinp',sinp
298        print 'cosp',cosp
299        print 'occ',occ
300        print 'Tiso',Tiso
301        print 'Tuij',Tuij
302    else:
303        print 'Absent'
304       
305def Refine(GPXfile):
306    ShowBanner()
307    Controls = GetControls(GPXfile)
308    ShowControls(Controls)
309   
310    phaseNames = GetPhaseNames(GPXfile)
311    phaseData = {}
312    for name in phaseNames: 
313        phaseData[name] =  GetPhaseData(GPXfile,name)
314    ShowPhaseData(phaseNames,phaseData)
315       
316    histograms = GetHistogramNames(GPXfile)
317    for hist in histograms:
318        if 'PWDR' in hist[:4]: 
319            print hist
320            PWDRdata = GetPWDRdata(GPXfile,hist)
321            print PWDRdata
322
323def main():
324    arg = sys.argv
325    if len(arg) > 1:
326        GPXfile = arg[1]
327        if not ospath.exists(GPXfile):
328            print 'ERROR - ',GPXfile," doesn't exist!"
329            exit()
330        GPXpath = ospath.dirname(arg[1])
331        Refine(GPXfile)
332    else:
333        print 'ERROR - missing filename'
334        exit()
335    print "Done"
336         
337if __name__ == '__main__':
338    main()
Note: See TracBrowser for help on using the repository browser.