1 | #GSASIIstructure - structure computation routines |
---|
2 | ########### SVN repository information ################### |
---|
3 | # $Date: 2012-01-19 22:56:20 +0000 (Thu, 19 Jan 2012) $ |
---|
4 | # $Author: vondreele $ |
---|
5 | # $Revision: 453 $ |
---|
6 | # $URL: trunk/GSASIIstruct.py $ |
---|
7 | # $Id: GSASIIstruct.py 453 2012-01-19 22:56:20Z vondreele $ |
---|
8 | ########### SVN repository information ################### |
---|
9 | import sys |
---|
10 | import os |
---|
11 | import os.path as ospath |
---|
12 | import numpy as np |
---|
13 | import numpy.linalg as nl |
---|
14 | import cPickle |
---|
15 | import time |
---|
16 | import math |
---|
17 | import GSASIIpath |
---|
18 | import GSASIIElem as G2el |
---|
19 | import GSASIIlattice as G2lat |
---|
20 | import GSASIIspc as G2spc |
---|
21 | import GSASIIpwd as G2pwd |
---|
22 | import GSASIImapvars as G2mv |
---|
23 | import GSASIImath as G2mth |
---|
24 | import scipy.optimize as so |
---|
25 | |
---|
26 | sind = lambda x: np.sin(x*np.pi/180.) |
---|
27 | cosd = lambda x: np.cos(x*np.pi/180.) |
---|
28 | tand = lambda x: np.tan(x*np.pi/180.) |
---|
29 | asind = lambda x: 180.*np.arcsin(x)/np.pi |
---|
30 | atan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
---|
31 | |
---|
32 | |
---|
33 | def ShowBanner(): |
---|
34 | print 80*'*' |
---|
35 | print ' General Structure Analysis System-II Crystal Structure Refinement' |
---|
36 | print ' by Robert B. Von Dreele, Argonne National Laboratory(C), 2010' |
---|
37 | print ' This product includes software developed by the UChicago Argonne, LLC,' |
---|
38 | print ' as Operator of Argonne National Laboratory.' |
---|
39 | print 80*'*','\n' |
---|
40 | |
---|
41 | def GetControls(GPXfile): |
---|
42 | ''' Returns dictionary of control items found in GSASII gpx file |
---|
43 | input: |
---|
44 | GPXfile = .gpx full file name |
---|
45 | return: |
---|
46 | Controls = dictionary of control items |
---|
47 | ''' |
---|
48 | Controls = {'deriv type':'analytical','min dM/M':0.0001,'shift factor':1.} |
---|
49 | file = open(GPXfile,'rb') |
---|
50 | while True: |
---|
51 | try: |
---|
52 | data = cPickle.load(file) |
---|
53 | except EOFError: |
---|
54 | break |
---|
55 | datum = data[0] |
---|
56 | if datum[0] == 'Controls': |
---|
57 | Controls.update(datum[1]) |
---|
58 | file.close() |
---|
59 | return Controls |
---|
60 | |
---|
61 | def ShowControls(Controls): |
---|
62 | print ' Least squares controls:' |
---|
63 | print ' Derivative type: ',Controls['deriv type'] |
---|
64 | print ' Minimum delta-M/M for convergence: ','%.2g'%(Controls['min dM/M']) |
---|
65 | print ' Initial shift factor: ','%.3f'%(Controls['shift factor']) |
---|
66 | |
---|
67 | def GetConstraints(GPXfile): |
---|
68 | constList = [] |
---|
69 | file = open(GPXfile,'rb') |
---|
70 | while True: |
---|
71 | try: |
---|
72 | data = cPickle.load(file) |
---|
73 | except EOFError: |
---|
74 | break |
---|
75 | datum = data[0] |
---|
76 | if datum[0] == 'Constraints': |
---|
77 | constDict = datum[1] |
---|
78 | for item in constDict: |
---|
79 | constList += constDict[item] |
---|
80 | file.close() |
---|
81 | constDict = [] |
---|
82 | constFlag = [] |
---|
83 | fixedList = [] |
---|
84 | for item in constList: |
---|
85 | if item[-2]: |
---|
86 | fixedList.append(str(item[-2])) |
---|
87 | else: |
---|
88 | fixedList.append('0') |
---|
89 | if item[-1]: |
---|
90 | constFlag.append(['VARY']) |
---|
91 | else: |
---|
92 | constFlag.append([]) |
---|
93 | itemDict = {} |
---|
94 | for term in item[:-2]: |
---|
95 | itemDict[term[1]] = term[0] |
---|
96 | constDict.append(itemDict) |
---|
97 | return constDict,constFlag,fixedList |
---|
98 | |
---|
99 | def GetPhaseNames(GPXfile): |
---|
100 | ''' Returns a list of phase names found under 'Phases' in GSASII gpx file |
---|
101 | input: |
---|
102 | GPXfile = gpx full file name |
---|
103 | return: |
---|
104 | PhaseNames = list of phase names |
---|
105 | ''' |
---|
106 | file = open(GPXfile,'rb') |
---|
107 | PhaseNames = [] |
---|
108 | while True: |
---|
109 | try: |
---|
110 | data = cPickle.load(file) |
---|
111 | except EOFError: |
---|
112 | break |
---|
113 | datum = data[0] |
---|
114 | if 'Phases' == datum[0]: |
---|
115 | for datus in data[1:]: |
---|
116 | PhaseNames.append(datus[0]) |
---|
117 | file.close() |
---|
118 | return PhaseNames |
---|
119 | |
---|
120 | def GetAllPhaseData(GPXfile,PhaseName): |
---|
121 | ''' Returns the entire dictionary for PhaseName from GSASII gpx file |
---|
122 | input: |
---|
123 | GPXfile = gpx full file name |
---|
124 | PhaseName = phase name |
---|
125 | return: |
---|
126 | phase dictionary |
---|
127 | ''' |
---|
128 | file = open(GPXfile,'rb') |
---|
129 | General = {} |
---|
130 | Atoms = [] |
---|
131 | while True: |
---|
132 | try: |
---|
133 | data = cPickle.load(file) |
---|
134 | except EOFError: |
---|
135 | break |
---|
136 | datum = data[0] |
---|
137 | if 'Phases' == datum[0]: |
---|
138 | for datus in data[1:]: |
---|
139 | if datus[0] == PhaseName: |
---|
140 | break |
---|
141 | file.close() |
---|
142 | return datus[1] |
---|
143 | |
---|
144 | def GetHistograms(GPXfile,hNames): |
---|
145 | """ Returns a dictionary of histograms found in GSASII gpx file |
---|
146 | input: |
---|
147 | GPXfile = .gpx full file name |
---|
148 | hNames = list of histogram names |
---|
149 | return: |
---|
150 | Histograms = dictionary of histograms (types = PWDR & HKLF) |
---|
151 | """ |
---|
152 | file = open(GPXfile,'rb') |
---|
153 | Histograms = {} |
---|
154 | while True: |
---|
155 | try: |
---|
156 | data = cPickle.load(file) |
---|
157 | except EOFError: |
---|
158 | break |
---|
159 | datum = data[0] |
---|
160 | hist = datum[0] |
---|
161 | if hist in hNames: |
---|
162 | if 'PWDR' in hist[:4]: |
---|
163 | PWDRdata = {} |
---|
164 | PWDRdata['Data'] = datum[1][1] #powder data arrays |
---|
165 | PWDRdata[data[2][0]] = data[2][1] #Limits |
---|
166 | PWDRdata[data[3][0]] = data[3][1] #Background |
---|
167 | PWDRdata[data[4][0]] = data[4][1] #Instrument parameters |
---|
168 | PWDRdata[data[5][0]] = data[5][1] #Sample parameters |
---|
169 | try: |
---|
170 | PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing |
---|
171 | except IndexError: |
---|
172 | PWDRdata['Reflection lists'] = {} |
---|
173 | |
---|
174 | Histograms[hist] = PWDRdata |
---|
175 | elif 'HKLF' in hist[:4]: |
---|
176 | HKLFdata = [] |
---|
177 | datum = data[0] |
---|
178 | HKLFdata = datum[1:][0] |
---|
179 | Histograms[hist] = HKLFdata |
---|
180 | file.close() |
---|
181 | return Histograms |
---|
182 | |
---|
183 | def GetHistogramNames(GPXfile,hType): |
---|
184 | """ Returns a list of histogram names found in GSASII gpx file |
---|
185 | input: |
---|
186 | GPXfile = .gpx full file name |
---|
187 | hType = list ['PWDR','HKLF'] |
---|
188 | return: |
---|
189 | HistogramNames = list of histogram names (types = PWDR & HKLF) |
---|
190 | """ |
---|
191 | file = open(GPXfile,'rb') |
---|
192 | HistogramNames = [] |
---|
193 | while True: |
---|
194 | try: |
---|
195 | data = cPickle.load(file) |
---|
196 | except EOFError: |
---|
197 | break |
---|
198 | datum = data[0] |
---|
199 | if datum[0][:4] in hType: |
---|
200 | HistogramNames.append(datum[0]) |
---|
201 | file.close() |
---|
202 | return HistogramNames |
---|
203 | |
---|
204 | def GetUsedHistogramsAndPhases(GPXfile): |
---|
205 | ''' Returns all histograms that are found in any phase |
---|
206 | and any phase that uses a histogram |
---|
207 | input: |
---|
208 | GPXfile = .gpx full file name |
---|
209 | return: |
---|
210 | Histograms = dictionary of histograms as {name:data,...} |
---|
211 | Phases = dictionary of phases that use histograms |
---|
212 | ''' |
---|
213 | phaseNames = GetPhaseNames(GPXfile) |
---|
214 | histoList = GetHistogramNames(GPXfile,['PWDR','HKLF']) |
---|
215 | allHistograms = GetHistograms(GPXfile,histoList) |
---|
216 | phaseData = {} |
---|
217 | for name in phaseNames: |
---|
218 | phaseData[name] = GetAllPhaseData(GPXfile,name) |
---|
219 | Histograms = {} |
---|
220 | Phases = {} |
---|
221 | for phase in phaseData: |
---|
222 | Phase = phaseData[phase] |
---|
223 | if Phase['Histograms']: |
---|
224 | if phase not in Phases: |
---|
225 | pId = phaseNames.index(phase) |
---|
226 | Phase['pId'] = pId |
---|
227 | Phases[phase] = Phase |
---|
228 | for hist in Phase['Histograms']: |
---|
229 | if hist not in Histograms: |
---|
230 | Histograms[hist] = allHistograms[hist] |
---|
231 | #future restraint, etc. histograms here |
---|
232 | hId = histoList.index(hist) |
---|
233 | Histograms[hist]['hId'] = hId |
---|
234 | return Histograms,Phases |
---|
235 | |
---|
236 | def GPXBackup(GPXfile,makeBack=True): |
---|
237 | import distutils.file_util as dfu |
---|
238 | GPXpath,GPXname = ospath.split(GPXfile) |
---|
239 | if GPXpath == '': GPXpath = '.' |
---|
240 | Name = ospath.splitext(GPXname)[0] |
---|
241 | files = os.listdir(GPXpath) |
---|
242 | last = 0 |
---|
243 | for name in files: |
---|
244 | name = name.split('.') |
---|
245 | if len(name) == 3 and name[0] == Name and 'bak' in name[1]: |
---|
246 | if makeBack: |
---|
247 | last = max(last,int(name[1].strip('bak'))+1) |
---|
248 | else: |
---|
249 | last = max(last,int(name[1].strip('bak'))) |
---|
250 | GPXback = ospath.join(GPXpath,ospath.splitext(GPXname)[0]+'.bak'+str(last)+'.gpx') |
---|
251 | dfu.copy_file(GPXfile,GPXback) |
---|
252 | return GPXback |
---|
253 | |
---|
254 | def SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,CovData,makeBack=True): |
---|
255 | ''' Updates gpxfile from all histograms that are found in any phase |
---|
256 | and any phase that used a histogram |
---|
257 | input: |
---|
258 | GPXfile = .gpx full file name |
---|
259 | Histograms = dictionary of histograms as {name:data,...} |
---|
260 | Phases = dictionary of phases that use histograms |
---|
261 | CovData = dictionary of refined variables, varyList, & covariance matrix |
---|
262 | makeBack = True if new backup of .gpx file is to be made; else use the last one made |
---|
263 | ''' |
---|
264 | |
---|
265 | GPXback = GPXBackup(GPXfile,makeBack) |
---|
266 | print '\n',135*'-' |
---|
267 | print 'Read from file:',GPXback |
---|
268 | print 'Save to file :',GPXfile |
---|
269 | infile = open(GPXback,'rb') |
---|
270 | outfile = open(GPXfile,'wb') |
---|
271 | while True: |
---|
272 | try: |
---|
273 | data = cPickle.load(infile) |
---|
274 | except EOFError: |
---|
275 | break |
---|
276 | datum = data[0] |
---|
277 | # print 'read: ',datum[0] |
---|
278 | if datum[0] == 'Phases': |
---|
279 | for iphase in range(len(data)): |
---|
280 | if data[iphase][0] in Phases: |
---|
281 | phaseName = data[iphase][0] |
---|
282 | data[iphase][1].update(Phases[phaseName]) |
---|
283 | elif datum[0] == 'Covariance': |
---|
284 | data[0][1] = CovData |
---|
285 | try: |
---|
286 | histogram = Histograms[datum[0]] |
---|
287 | # print 'found ',datum[0] |
---|
288 | data[0][1][1] = histogram['Data'] |
---|
289 | for datus in data[1:]: |
---|
290 | # print ' read: ',datus[0] |
---|
291 | if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: |
---|
292 | datus[1] = histogram[datus[0]] |
---|
293 | except KeyError: |
---|
294 | pass |
---|
295 | |
---|
296 | cPickle.dump(data,outfile,1) |
---|
297 | infile.close() |
---|
298 | outfile.close() |
---|
299 | print 'refinement save successful' |
---|
300 | |
---|
301 | def SetSeqResult(GPXfile,Histograms,SeqResult): |
---|
302 | GPXback = GPXBackup(GPXfile) |
---|
303 | print '\n',135*'-' |
---|
304 | print 'Read from file:',GPXback |
---|
305 | print 'Save to file :',GPXfile |
---|
306 | infile = open(GPXback,'rb') |
---|
307 | outfile = open(GPXfile,'wb') |
---|
308 | while True: |
---|
309 | try: |
---|
310 | data = cPickle.load(infile) |
---|
311 | except EOFError: |
---|
312 | break |
---|
313 | datum = data[0] |
---|
314 | if datum[0] == 'Sequental results': |
---|
315 | data[0][1] = SeqResult |
---|
316 | try: |
---|
317 | histogram = Histograms[datum[0]] |
---|
318 | data[0][1][1] = histogram['Data'] |
---|
319 | for datus in data[1:]: |
---|
320 | if datus[0] in ['Background','Instrument Parameters','Sample Parameters','Reflection Lists']: |
---|
321 | datus[1] = histogram[datus[0]] |
---|
322 | except KeyError: |
---|
323 | pass |
---|
324 | |
---|
325 | cPickle.dump(data,outfile,1) |
---|
326 | infile.close() |
---|
327 | outfile.close() |
---|
328 | print 'refinement save successful' |
---|
329 | |
---|
330 | |
---|
331 | def GetPWDRdata(GPXfile,PWDRname): |
---|
332 | ''' Returns powder data from GSASII gpx file |
---|
333 | input: |
---|
334 | GPXfile = .gpx full file name |
---|
335 | PWDRname = powder histogram name as obtained from GetHistogramNames |
---|
336 | return: |
---|
337 | PWDRdata = powder data dictionary with: |
---|
338 | Data - powder data arrays, Limits, Instrument Parameters, Sample Parameters |
---|
339 | |
---|
340 | ''' |
---|
341 | file = open(GPXfile,'rb') |
---|
342 | PWDRdata = {} |
---|
343 | while True: |
---|
344 | try: |
---|
345 | data = cPickle.load(file) |
---|
346 | except EOFError: |
---|
347 | break |
---|
348 | datum = data[0] |
---|
349 | if datum[0] == PWDRname: |
---|
350 | PWDRdata['Data'] = datum[1][1] #powder data arrays |
---|
351 | PWDRdata[data[2][0]] = data[2][1] #Limits |
---|
352 | PWDRdata[data[3][0]] = data[3][1] #Background |
---|
353 | PWDRdata[data[4][0]] = data[4][1] #Instrument parameters |
---|
354 | PWDRdata[data[5][0]] = data[5][1] #Sample parameters |
---|
355 | try: |
---|
356 | PWDRdata[data[9][0]] = data[9][1] #Reflection lists might be missing |
---|
357 | except IndexError: |
---|
358 | PWDRdata['Reflection lists'] = {} |
---|
359 | file.close() |
---|
360 | return PWDRdata |
---|
361 | |
---|
362 | def GetHKLFdata(GPXfile,HKLFname): |
---|
363 | ''' Returns single crystal data from GSASII gpx file |
---|
364 | input: |
---|
365 | GPXfile = .gpx full file name |
---|
366 | HKLFname = single crystal histogram name as obtained from GetHistogramNames |
---|
367 | return: |
---|
368 | HKLFdata = single crystal data list of reflections: for each reflection: |
---|
369 | HKLF = [np.array([h,k,l]),FoSq,sigFoSq,FcSq,Fcp,Fcpp,phase] |
---|
370 | ''' |
---|
371 | file = open(GPXfile,'rb') |
---|
372 | HKLFdata = [] |
---|
373 | while True: |
---|
374 | try: |
---|
375 | data = cPickle.load(file) |
---|
376 | except EOFError: |
---|
377 | break |
---|
378 | datum = data[0] |
---|
379 | if datum[0] == HKLFname: |
---|
380 | HKLFdata = datum[1:][0] |
---|
381 | file.close() |
---|
382 | return HKLFdata |
---|
383 | |
---|
384 | def GetFFtable(General): |
---|
385 | ''' returns a dictionary of form factor data for atom types found in General |
---|
386 | input: |
---|
387 | General = dictionary of phase info.; includes AtomTypes |
---|
388 | return: |
---|
389 | FFtable = dictionary of form factor data; key is atom type |
---|
390 | ''' |
---|
391 | atomTypes = General['AtomTypes'] |
---|
392 | FFtable = {} |
---|
393 | for El in atomTypes: |
---|
394 | FFs = G2el.GetFormFactorCoeff(El.split('+')[0].split('-')[0]) |
---|
395 | for item in FFs: |
---|
396 | if item['Symbol'] == El.upper(): |
---|
397 | FFtable[El] = item |
---|
398 | return FFtable |
---|
399 | |
---|
400 | def GetBLtable(General): |
---|
401 | ''' returns a dictionary of neutron scattering length data for atom types & isotopes found in General |
---|
402 | input: |
---|
403 | General = dictionary of phase info.; includes AtomTypes & Isotopes |
---|
404 | return: |
---|
405 | BLtable = dictionary of scattering length data; key is atom type |
---|
406 | ''' |
---|
407 | atomTypes = General['AtomTypes'] |
---|
408 | BLtable = {} |
---|
409 | isotopes = General['Isotopes'] |
---|
410 | isotope = General['Isotope'] |
---|
411 | for El in atomTypes: |
---|
412 | BLtable[El] = [isotope[El],isotopes[El][isotope[El]]] |
---|
413 | return BLtable |
---|
414 | |
---|
415 | def GetPawleyConstr(SGLaue,PawleyRef,pawleyVary): |
---|
416 | if SGLaue in ['-1','2/m','mmm']: |
---|
417 | return #no Pawley symmetry required constraints |
---|
418 | for i,varyI in enumerate(pawleyVary): |
---|
419 | refI = int(varyI.split(':')[-1]) |
---|
420 | ih,ik,il = PawleyRef[refI][:3] |
---|
421 | for varyJ in pawleyVary[0:i]: |
---|
422 | refJ = int(varyJ.split(':')[-1]) |
---|
423 | jh,jk,jl = PawleyRef[refJ][:3] |
---|
424 | if SGLaue in ['4/m','4/mmm']: |
---|
425 | isum = ih**2+ik**2 |
---|
426 | jsum = jh**2+jk**2 |
---|
427 | if abs(il) == abs(jl) and isum == jsum: |
---|
428 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
---|
429 | elif SGLaue in ['3R','3mR']: |
---|
430 | isum = ih**2+ik**2+il**2 |
---|
431 | jsum = jh**2+jk**2*jl**2 |
---|
432 | isum2 = ih*ik+ih*il+ik*il |
---|
433 | jsum2 = jh*jk+jh*jl+jk*jl |
---|
434 | if isum == jsum and isum2 == jsum2: |
---|
435 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
---|
436 | elif SGLaue in ['3','3m1','31m','6/m','6/mmm']: |
---|
437 | isum = ih**2+ik**2+ih*ik |
---|
438 | jsum = jh**2+jk**2+jh*jk |
---|
439 | if abs(il) == abs(jl) and isum == jsum: |
---|
440 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
---|
441 | elif SGLaue in ['m3','m3m']: |
---|
442 | isum = ih**2+ik**2+il**2 |
---|
443 | jsum = jh**2+jk**2+jl**2 |
---|
444 | if isum == jsum: |
---|
445 | G2mv.StoreEquivalence(varyJ,(varyI,)) |
---|
446 | |
---|
447 | def cellVary(pfx,SGData): |
---|
448 | if SGData['SGLaue'] in ['-1',]: |
---|
449 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
---|
450 | elif SGData['SGLaue'] in ['2/m',]: |
---|
451 | if SGData['SGUniq'] == 'a': |
---|
452 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3'] |
---|
453 | elif SGData['SGUniq'] == 'b': |
---|
454 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A4'] |
---|
455 | else: |
---|
456 | return [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A5'] |
---|
457 | elif SGData['SGLaue'] in ['mmm',]: |
---|
458 | return [pfx+'A0',pfx+'A1',pfx+'A2'] |
---|
459 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
460 | return [pfx+'A0',pfx+'A2'] |
---|
461 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
462 | return [pfx+'A0',pfx+'A2'] |
---|
463 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
464 | return [pfx+'A0',pfx+'A3'] |
---|
465 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
466 | return [pfx+'A0',] |
---|
467 | |
---|
468 | |
---|
469 | def GetPhaseData(PhaseData,Print=True): |
---|
470 | |
---|
471 | def PrintFFtable(FFtable): |
---|
472 | print '\n X-ray scattering factors:' |
---|
473 | print ' Symbol fa fb fc' |
---|
474 | print 99*'-' |
---|
475 | for Ename in FFtable: |
---|
476 | ffdata = FFtable[Ename] |
---|
477 | fa = ffdata['fa'] |
---|
478 | fb = ffdata['fb'] |
---|
479 | print ' %8s %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f %9.5f' % \ |
---|
480 | (Ename.ljust(8),fa[0],fa[1],fa[2],fa[3],fb[0],fb[1],fb[2],fb[3],ffdata['fc']) |
---|
481 | |
---|
482 | def PrintBLtable(BLtable): |
---|
483 | print '\n Neutron scattering factors:' |
---|
484 | print ' Symbol isotope mass b resonant terms' |
---|
485 | print 99*'-' |
---|
486 | for Ename in BLtable: |
---|
487 | bldata = BLtable[Ename] |
---|
488 | isotope = bldata[0] |
---|
489 | mass = bldata[1][0] |
---|
490 | blen = bldata[1][1] |
---|
491 | bres = [] |
---|
492 | if len(bldata[1]) > 2: |
---|
493 | bres = bldata[1][2:] |
---|
494 | line = ' %8s%11s %10.3f %8.3f'%(Ename.ljust(8),isotope.center(11),mass,blen) |
---|
495 | for item in bres: |
---|
496 | line += '%10.5g'%(item) |
---|
497 | print line |
---|
498 | |
---|
499 | def PrintAtoms(General,Atoms): |
---|
500 | print '\n Atoms:' |
---|
501 | line = ' name type refine? x y z '+ \ |
---|
502 | ' frac site sym mult I/A Uiso U11 U22 U33 U12 U13 U23' |
---|
503 | if General['Type'] == 'magnetic': |
---|
504 | line += ' Mx My Mz' |
---|
505 | elif General['Type'] == 'macromolecular': |
---|
506 | line = ' res no residue chain '+line |
---|
507 | print line |
---|
508 | if General['Type'] == 'nuclear': |
---|
509 | print 135*'-' |
---|
510 | for i,at in enumerate(Atoms): |
---|
511 | line = '%7s'%(at[0])+'%7s'%(at[1])+'%7s'%(at[2])+'%10.5f'%(at[3])+'%10.5f'%(at[4])+ \ |
---|
512 | '%10.5f'%(at[5])+'%8.3f'%(at[6])+'%7s'%(at[7])+'%5d'%(at[8])+'%5s'%(at[9]) |
---|
513 | if at[9] == 'I': |
---|
514 | line += '%8.4f'%(at[10])+48*' ' |
---|
515 | else: |
---|
516 | line += 8*' ' |
---|
517 | for j in range(6): |
---|
518 | line += '%8.4f'%(at[11+j]) |
---|
519 | print line |
---|
520 | |
---|
521 | def PrintTexture(textureData): |
---|
522 | topstr = '\n Spherical harmonics texture: Order:' + \ |
---|
523 | str(textureData['Order']) |
---|
524 | if textureData['Order']: |
---|
525 | print topstr+' Refine? '+str(textureData['SH Coeff'][0]) |
---|
526 | else: |
---|
527 | print topstr |
---|
528 | return |
---|
529 | names = ['omega','chi','phi'] |
---|
530 | line = '\n' |
---|
531 | for name in names: |
---|
532 | line += ' SH '+name+':'+'%12.4f'%(textureData['Sample '+name][1])+' Refine? '+str(textureData['Sample '+name][0]) |
---|
533 | print line |
---|
534 | print '\n Texture coefficients:' |
---|
535 | ptlbls = ' names :' |
---|
536 | ptstr = ' values:' |
---|
537 | SHcoeff = textureData['SH Coeff'][1] |
---|
538 | for item in SHcoeff: |
---|
539 | ptlbls += '%12s'%(item) |
---|
540 | ptstr += '%12.4f'%(SHcoeff[item]) |
---|
541 | print ptlbls |
---|
542 | print ptstr |
---|
543 | |
---|
544 | if Print: print ' Phases:' |
---|
545 | phaseVary = [] |
---|
546 | phaseDict = {} |
---|
547 | phaseConstr = {} |
---|
548 | pawleyLookup = {} |
---|
549 | FFtables = {} #scattering factors - xrays |
---|
550 | BLtables = {} # neutrons |
---|
551 | Natoms = {} |
---|
552 | AtMults = {} |
---|
553 | AtIA = {} |
---|
554 | shModels = ['cylindrical','none','shear - 2/m','rolling - mmm'] |
---|
555 | SamSym = dict(zip(shModels,['0','-1','2/m','mmm'])) |
---|
556 | for name in PhaseData: |
---|
557 | General = PhaseData[name]['General'] |
---|
558 | pId = PhaseData[name]['pId'] |
---|
559 | pfx = str(pId)+'::' |
---|
560 | FFtable = GetFFtable(General) |
---|
561 | BLtable = GetBLtable(General) |
---|
562 | FFtables.update(FFtable) |
---|
563 | BLtables.update(BLtable) |
---|
564 | Atoms = PhaseData[name]['Atoms'] |
---|
565 | try: |
---|
566 | PawleyRef = PhaseData[name]['Pawley ref'] |
---|
567 | except KeyError: |
---|
568 | PawleyRef = [] |
---|
569 | SGData = General['SGData'] |
---|
570 | SGtext = G2spc.SGPrint(SGData) |
---|
571 | cell = General['Cell'] |
---|
572 | A = G2lat.cell2A(cell[1:7]) |
---|
573 | phaseDict.update({pfx+'A0':A[0],pfx+'A1':A[1],pfx+'A2':A[2],pfx+'A3':A[3],pfx+'A4':A[4],pfx+'A5':A[5]}) |
---|
574 | if cell[0]: |
---|
575 | phaseVary += cellVary(pfx,SGData) |
---|
576 | Natoms[pfx] = 0 |
---|
577 | if Atoms: |
---|
578 | if General['Type'] == 'nuclear': |
---|
579 | Natoms[pfx] = len(Atoms) |
---|
580 | for i,at in enumerate(Atoms): |
---|
581 | phaseDict.update({pfx+'Atype:'+str(i):at[1],pfx+'Afrac:'+str(i):at[6],pfx+'Amul:'+str(i):at[8], |
---|
582 | pfx+'Ax:'+str(i):at[3],pfx+'Ay:'+str(i):at[4],pfx+'Az:'+str(i):at[5], |
---|
583 | pfx+'dAx:'+str(i):0.,pfx+'dAy:'+str(i):0.,pfx+'dAz:'+str(i):0., #refined shifts for x,y,z |
---|
584 | pfx+'AI/A:'+str(i):at[9],}) |
---|
585 | if at[9] == 'I': |
---|
586 | phaseDict[pfx+'AUiso:'+str(i)] = at[10] |
---|
587 | else: |
---|
588 | phaseDict.update({pfx+'AU11:'+str(i):at[11],pfx+'AU22:'+str(i):at[12],pfx+'AU33:'+str(i):at[13], |
---|
589 | pfx+'AU12:'+str(i):at[14],pfx+'AU13:'+str(i):at[15],pfx+'AU23:'+str(i):at[16]}) |
---|
590 | if 'F' in at[2]: |
---|
591 | phaseVary.append(pfx+'Afrac:'+str(i)) |
---|
592 | if 'X' in at[2]: |
---|
593 | xId,xCoef = G2spc.GetCSxinel(at[7]) |
---|
594 | delnames = [pfx+'dAx:'+str(i),pfx+'dAy:'+str(i),pfx+'dAz:'+str(i)] |
---|
595 | for j in range(3): |
---|
596 | if xId[j] > 0: |
---|
597 | phaseVary.append(delnames[j]) |
---|
598 | for k in range(j): |
---|
599 | if xId[j] == xId[k]: |
---|
600 | G2mv.StoreEquivalence(delnames[k],((delnames[j],xCoef[j]),)) |
---|
601 | if 'U' in at[2]: |
---|
602 | if at[9] == 'I': |
---|
603 | phaseVary.append(pfx+'AUiso:'+str(i)) |
---|
604 | else: |
---|
605 | uId,uCoef = G2spc.GetCSuinel(at[7])[:2] |
---|
606 | names = [pfx+'AU11:'+str(i),pfx+'AU22:'+str(i),pfx+'AU33:'+str(i), |
---|
607 | pfx+'AU12:'+str(i),pfx+'AU13:'+str(i),pfx+'AU23:'+str(i)] |
---|
608 | for j in range(6): |
---|
609 | if uId[j] > 0: |
---|
610 | phaseVary.append(names[j]) |
---|
611 | for k in range(j): |
---|
612 | if uId[j] == uId[k]: |
---|
613 | G2mv.StoreEquivalence(names[k],((names[j],uCoef[j]),)) |
---|
614 | # elif General['Type'] == 'magnetic': |
---|
615 | # elif General['Type'] == 'macromolecular': |
---|
616 | |
---|
617 | |
---|
618 | if 'SH Texture' in General: |
---|
619 | textureData = General['SH Texture'] |
---|
620 | phaseDict[pfx+'SHmodel'] = SamSym[textureData['Model']] |
---|
621 | phaseDict[pfx+'SHorder'] = textureData['Order'] |
---|
622 | for name in ['omega','chi','phi']: |
---|
623 | phaseDict[pfx+'SH '+name] = textureData['Sample '+name][1] |
---|
624 | if textureData['Sample '+name][0]: |
---|
625 | phaseVary.append(pfx+'SH '+name) |
---|
626 | for name in textureData['SH Coeff'][1]: |
---|
627 | phaseDict[pfx+name] = textureData['SH Coeff'][1][name] |
---|
628 | if textureData['SH Coeff'][0]: |
---|
629 | phaseVary.append(pfx+name) |
---|
630 | |
---|
631 | if Print: |
---|
632 | print '\n Phase name: ',General['Name'] |
---|
633 | print 135*'-' |
---|
634 | PrintFFtable(FFtable) |
---|
635 | PrintBLtable(BLtable) |
---|
636 | print '' |
---|
637 | for line in SGtext: print line |
---|
638 | PrintAtoms(General,Atoms) |
---|
639 | print '\n Unit cell: a =','%.5f'%(cell[1]),' b =','%.5f'%(cell[2]),' c =','%.5f'%(cell[3]), \ |
---|
640 | ' alpha =','%.3f'%(cell[4]),' beta =','%.3f'%(cell[5]),' gamma =', \ |
---|
641 | '%.3f'%(cell[6]),' volume =','%.3f'%(cell[7]),' Refine?',cell[0] |
---|
642 | if 'SH Texture' in General: |
---|
643 | PrintTexture(textureData) |
---|
644 | |
---|
645 | elif PawleyRef: |
---|
646 | pawleyVary = [] |
---|
647 | for i,refl in enumerate(PawleyRef): |
---|
648 | phaseDict[pfx+'PWLref:'+str(i)] = refl[6] |
---|
649 | pawleyLookup[pfx+'%d,%d,%d'%(refl[0],refl[1],refl[2])] = i |
---|
650 | if refl[5]: |
---|
651 | pawleyVary.append(pfx+'PWLref:'+str(i)) |
---|
652 | GetPawleyConstr(SGData['SGLaue'],PawleyRef,pawleyVary) #does G2mv.StoreEquivalence |
---|
653 | phaseVary += pawleyVary |
---|
654 | |
---|
655 | return Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables |
---|
656 | |
---|
657 | def getVCov(varyNames,varyList,covMatrix): |
---|
658 | vcov = np.zeros((len(varyNames),len(varyNames))) |
---|
659 | for i1,name1 in enumerate(varyNames): |
---|
660 | for i2,name2 in enumerate(varyNames): |
---|
661 | try: |
---|
662 | vcov[i1][i2] = covMatrix[varyList.index(name1)][varyList.index(name2)] |
---|
663 | except ValueError: |
---|
664 | vcov[i1][i2] = 0.0 |
---|
665 | return vcov |
---|
666 | |
---|
667 | def cellFill(pfx,SGData,parmDict,sigDict): |
---|
668 | if SGData['SGLaue'] in ['-1',]: |
---|
669 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
670 | parmDict[pfx+'A3'],parmDict[pfx+'A4'],parmDict[pfx+'A5']] |
---|
671 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
672 | sigDict[pfx+'A3'],sigDict[pfx+'A4'],sigDict[pfx+'A5']] |
---|
673 | elif SGData['SGLaue'] in ['2/m',]: |
---|
674 | if SGData['SGUniq'] == 'a': |
---|
675 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
676 | parmDict[pfx+'A3'],0,0] |
---|
677 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
678 | sigDict[pfx+'A3'],0,0] |
---|
679 | elif SGData['SGUniq'] == 'b': |
---|
680 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
681 | 0,parmDict[pfx+'A4'],0] |
---|
682 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
683 | 0,sigDict[pfx+'A4'],0] |
---|
684 | else: |
---|
685 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'], |
---|
686 | 0,0,parmDict[pfx+'A5']] |
---|
687 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'], |
---|
688 | 0,0,sigDict[pfx+'A5']] |
---|
689 | elif SGData['SGLaue'] in ['mmm',]: |
---|
690 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A1'],parmDict[pfx+'A2'],0,0,0] |
---|
691 | sigA = [sigDict[pfx+'A0'],sigDict[pfx+'A1'],sigDict[pfx+'A2'],0,0,0] |
---|
692 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
693 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'],0,0,0] |
---|
694 | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
---|
695 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
696 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A2'], |
---|
697 | parmDict[pfx+'A0'],0,0] |
---|
698 | sigA = [sigDict[pfx+'A0'],0,sigDict[pfx+'A2'],0,0,0] |
---|
699 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
700 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'], |
---|
701 | parmDict[pfx+'A3'],parmDict[pfx+'A3'],parmDict[pfx+'A3']] |
---|
702 | sigA = [sigDict[pfx+'A0'],0,0,sigDict[pfx+'A3'],0,0] |
---|
703 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
704 | A = [parmDict[pfx+'A0'],parmDict[pfx+'A0'],parmDict[pfx+'A0'],0,0,0] |
---|
705 | sigA = [sigDict[pfx+'A0'],0,0,0,0,0] |
---|
706 | return A,sigA |
---|
707 | |
---|
708 | def getCellEsd(pfx,SGData,A,covData): |
---|
709 | dpr = 180./np.pi |
---|
710 | rVsq = G2lat.calc_rVsq(A) |
---|
711 | G,g = G2lat.A2Gmat(A) #get recip. & real metric tensors |
---|
712 | cell = np.array(G2lat.Gmat2cell(g)) #real cell |
---|
713 | cellst = np.array(G2lat.Gmat2cell(G)) #recip. cell |
---|
714 | scos = cosd(cellst[3:6]) |
---|
715 | ssin = sind(cellst[3:6]) |
---|
716 | scot = scos/ssin |
---|
717 | rcos = cosd(cell[3:6]) |
---|
718 | rsin = sind(cell[3:6]) |
---|
719 | rcot = rcos/rsin |
---|
720 | RMnames = [pfx+'A0',pfx+'A1',pfx+'A2',pfx+'A3',pfx+'A4',pfx+'A5'] |
---|
721 | varyList = covData['varyList'] |
---|
722 | covMatrix = covData['covMatrix'] |
---|
723 | vcov = getVCov(RMnames,varyList,covMatrix) |
---|
724 | Ax = np.array(A) |
---|
725 | Ax[3:] /= 2. |
---|
726 | drVdA = np.array([Ax[1]*Ax[2]-Ax[5]**2,Ax[0]*Ax[2]-Ax[4]**2,Ax[0]*Ax[1]-Ax[3]**2, |
---|
727 | Ax[4]*Ax[5]-Ax[2]*Ax[3],Ax[3]*Ax[5]-Ax[1]*Ax[4],Ax[3]*Ax[4]-Ax[0]*Ax[5]]) |
---|
728 | srcvlsq = np.inner(drVdA,np.inner(vcov,drVdA.T)) |
---|
729 | Vol = 1/np.sqrt(rVsq) |
---|
730 | sigVol = Vol**3*np.sqrt(srcvlsq)/2. |
---|
731 | R123 = Ax[0]*Ax[1]*Ax[2] |
---|
732 | dsasdg = np.zeros((3,6)) |
---|
733 | dadg = np.zeros((6,6)) |
---|
734 | for i0 in range(3): #0 1 2 |
---|
735 | i1 = (i0+1)%3 #1 2 0 |
---|
736 | i2 = (i1+1)%3 #2 0 1 |
---|
737 | i3 = 5-i2 #3 5 4 |
---|
738 | i4 = 5-i1 #4 3 5 |
---|
739 | i5 = 5-i0 #5 4 3 |
---|
740 | dsasdg[i0][i1] = 0.5*scot[i0]*scos[i0]/Ax[i1] |
---|
741 | dsasdg[i0][i2] = 0.5*scot[i0]*scos[i0]/Ax[i2] |
---|
742 | dsasdg[i0][i5] = -scot[i0]/np.sqrt(Ax[i1]*Ax[i2]) |
---|
743 | denmsq = Ax[i0]*(R123-Ax[i1]*Ax[i4]**2-Ax[i2]*Ax[i3]**2+(Ax[i4]*Ax[i3])**2) |
---|
744 | denom = np.sqrt(denmsq) |
---|
745 | dadg[i5][i0] = -Ax[i5]/denom-rcos[i0]/denmsq*(R123-0.5*Ax[i1]*Ax[i4]**2-0.5*Ax[i2]*Ax[i3]**2) |
---|
746 | dadg[i5][i1] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i2]-Ax[i0]*Ax[i4]**2) |
---|
747 | dadg[i5][i2] = -0.5*rcos[i0]/denmsq*(Ax[i0]**2*Ax[i1]-Ax[i0]*Ax[i3]**2) |
---|
748 | dadg[i5][i3] = Ax[i4]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i2]*Ax[i3]-Ax[i3]*Ax[i4]**2) |
---|
749 | dadg[i5][i4] = Ax[i3]/denom+rcos[i0]/denmsq*(Ax[i0]*Ax[i1]*Ax[i4]-Ax[i3]**2*Ax[i4]) |
---|
750 | dadg[i5][i5] = -Ax[i0]/denom |
---|
751 | for i0 in range(3): |
---|
752 | i1 = (i0+1)%3 |
---|
753 | i2 = (i1+1)%3 |
---|
754 | i3 = 5-i2 |
---|
755 | for ij in range(6): |
---|
756 | dadg[i0][ij] = cell[i0]*(rcot[i2]*dadg[i3][ij]/rsin[i2]-dsasdg[i1][ij]/ssin[i1]) |
---|
757 | if ij == i0: |
---|
758 | dadg[i0][ij] = dadg[i0][ij]-0.5*cell[i0]/Ax[i0] |
---|
759 | dadg[i3][ij] = -dadg[i3][ij]*rsin[2-i0]*dpr |
---|
760 | sigMat = np.inner(dadg,np.inner(vcov,dadg.T)) |
---|
761 | var = np.diag(sigMat) |
---|
762 | CS = np.where(var>0.,np.sqrt(var),0.) |
---|
763 | cellSig = [CS[0],CS[1],CS[2],CS[5],CS[4],CS[3],sigVol] #exchange sig(alp) & sig(gam) to get in right order |
---|
764 | return cellSig |
---|
765 | |
---|
766 | def SetPhaseData(parmDict,sigDict,Phases,covData): |
---|
767 | |
---|
768 | def PrintAtomsAndSig(General,Atoms,atomsSig): |
---|
769 | print '\n Atoms:' |
---|
770 | line = ' name x y z frac Uiso U11 U22 U33 U12 U13 U23' |
---|
771 | if General['Type'] == 'magnetic': |
---|
772 | line += ' Mx My Mz' |
---|
773 | elif General['Type'] == 'macromolecular': |
---|
774 | line = ' res no residue chain '+line |
---|
775 | print line |
---|
776 | if General['Type'] == 'nuclear': |
---|
777 | print 135*'-' |
---|
778 | fmt = {0:'%7s',1:'%7s',3:'%10.5f',4:'%10.5f',5:'%10.5f',6:'%8.3f',10:'%8.5f', |
---|
779 | 11:'%8.5f',12:'%8.5f',13:'%8.5f',14:'%8.5f',15:'%8.5f',16:'%8.5f'} |
---|
780 | noFXsig = {3:[10*' ','%10s'],4:[10*' ','%10s'],5:[10*' ','%10s'],6:[8*' ','%8s']} |
---|
781 | for i,at in enumerate(Atoms): |
---|
782 | name = fmt[0]%(at[0])+fmt[1]%(at[1])+':' |
---|
783 | valstr = ' values:' |
---|
784 | sigstr = ' sig :' |
---|
785 | for ind in [3,4,5,6]: |
---|
786 | sigind = str(i)+':'+str(ind) |
---|
787 | valstr += fmt[ind]%(at[ind]) |
---|
788 | if sigind in atomsSig: |
---|
789 | sigstr += fmt[ind]%(atomsSig[sigind]) |
---|
790 | else: |
---|
791 | sigstr += noFXsig[ind][1]%(noFXsig[ind][0]) |
---|
792 | if at[9] == 'I': |
---|
793 | valstr += fmt[10]%(at[10]) |
---|
794 | if str(i)+':10' in atomsSig: |
---|
795 | sigstr += fmt[10]%(atomsSig[str(i)+':10']) |
---|
796 | else: |
---|
797 | sigstr += 8*' ' |
---|
798 | else: |
---|
799 | valstr += 8*' ' |
---|
800 | sigstr += 8*' ' |
---|
801 | for ind in [11,12,13,14,15,16]: |
---|
802 | sigind = str(i)+':'+str(ind) |
---|
803 | valstr += fmt[ind]%(at[ind]) |
---|
804 | if sigind in atomsSig: |
---|
805 | sigstr += fmt[ind]%(atomsSig[sigind]) |
---|
806 | else: |
---|
807 | sigstr += 8*' ' |
---|
808 | print name |
---|
809 | print valstr |
---|
810 | print sigstr |
---|
811 | |
---|
812 | def PrintSHtextureAndSig(textureData,SHtextureSig): |
---|
813 | print '\n Spherical harmonics texture: Order:' + str(textureData['Order']) |
---|
814 | names = ['omega','chi','phi'] |
---|
815 | namstr = ' names :' |
---|
816 | ptstr = ' values:' |
---|
817 | sigstr = ' esds :' |
---|
818 | for name in names: |
---|
819 | namstr += '%12s'%(name) |
---|
820 | ptstr += '%12.3f'%(textureData['Sample '+name][1]) |
---|
821 | if 'Sample '+name in SHtextureSig: |
---|
822 | sigstr += '%12.3f'%(SHtextureSig['Sample '+name]) |
---|
823 | else: |
---|
824 | sigstr += 12*' ' |
---|
825 | print namstr |
---|
826 | print ptstr |
---|
827 | print sigstr |
---|
828 | print '\n Texture coefficients:' |
---|
829 | namstr = ' names :' |
---|
830 | ptstr = ' values:' |
---|
831 | sigstr = ' esds :' |
---|
832 | SHcoeff = textureData['SH Coeff'][1] |
---|
833 | for name in SHcoeff: |
---|
834 | namstr += '%12s'%(name) |
---|
835 | ptstr += '%12.3f'%(SHcoeff[name]) |
---|
836 | if name in SHtextureSig: |
---|
837 | sigstr += '%12.3f'%(SHtextureSig[name]) |
---|
838 | else: |
---|
839 | sigstr += 12*' ' |
---|
840 | print namstr |
---|
841 | print ptstr |
---|
842 | print sigstr |
---|
843 | |
---|
844 | |
---|
845 | print '\n Phases:' |
---|
846 | for phase in Phases: |
---|
847 | print ' Result for phase: ',phase |
---|
848 | Phase = Phases[phase] |
---|
849 | General = Phase['General'] |
---|
850 | SGData = General['SGData'] |
---|
851 | Atoms = Phase['Atoms'] |
---|
852 | cell = General['Cell'] |
---|
853 | pId = Phase['pId'] |
---|
854 | pfx = str(pId)+'::' |
---|
855 | if cell[0]: |
---|
856 | A,sigA = cellFill(pfx,SGData,parmDict,sigDict) |
---|
857 | cellSig = getCellEsd(pfx,SGData,A,covData) #includes sigVol |
---|
858 | print ' Reciprocal metric tensor: ' |
---|
859 | ptfmt = "%15.9f" |
---|
860 | names = ['A11','A22','A33','A12','A13','A23'] |
---|
861 | namstr = ' names :' |
---|
862 | ptstr = ' values:' |
---|
863 | sigstr = ' esds :' |
---|
864 | for name,a,siga in zip(names,A,sigA): |
---|
865 | namstr += '%15s'%(name) |
---|
866 | ptstr += ptfmt%(a) |
---|
867 | if siga: |
---|
868 | sigstr += ptfmt%(siga) |
---|
869 | else: |
---|
870 | sigstr += 15*' ' |
---|
871 | print namstr |
---|
872 | print ptstr |
---|
873 | print sigstr |
---|
874 | cell[1:7] = G2lat.A2cell(A) |
---|
875 | cell[7] = G2lat.calc_V(A) |
---|
876 | print ' New unit cell:' |
---|
877 | ptfmt = ["%12.6f","%12.6f","%12.6f","%12.4f","%12.4f","%12.4f","%12.3f"] |
---|
878 | names = ['a','b','c','alpha','beta','gamma','Volume'] |
---|
879 | namstr = ' names :' |
---|
880 | ptstr = ' values:' |
---|
881 | sigstr = ' esds :' |
---|
882 | for name,fmt,a,siga in zip(names,ptfmt,cell[1:8],cellSig): |
---|
883 | namstr += '%12s'%(name) |
---|
884 | ptstr += fmt%(a) |
---|
885 | if siga: |
---|
886 | sigstr += fmt%(siga) |
---|
887 | else: |
---|
888 | sigstr += 12*' ' |
---|
889 | print namstr |
---|
890 | print ptstr |
---|
891 | print sigstr |
---|
892 | |
---|
893 | if 'Pawley' in Phase['General']['Type']: |
---|
894 | pawleyRef = Phase['Pawley ref'] |
---|
895 | for i,refl in enumerate(pawleyRef): |
---|
896 | key = pfx+'PWLref:'+str(i) |
---|
897 | refl[6] = abs(parmDict[key]) #suppress negative Fsq |
---|
898 | if key in sigDict: |
---|
899 | refl[7] = sigDict[key] |
---|
900 | else: |
---|
901 | refl[7] = 0 |
---|
902 | else: |
---|
903 | atomsSig = {} |
---|
904 | if General['Type'] == 'nuclear': |
---|
905 | for i,at in enumerate(Atoms): |
---|
906 | names = {3:pfx+'Ax:'+str(i),4:pfx+'Ay:'+str(i),5:pfx+'Az:'+str(i),6:pfx+'Afrac:'+str(i), |
---|
907 | 10:pfx+'AUiso:'+str(i),11:pfx+'AU11:'+str(i),12:pfx+'AU22:'+str(i),13:pfx+'AU33:'+str(i), |
---|
908 | 14:pfx+'AU12:'+str(i),15:pfx+'AU13:'+str(i),16:pfx+'AU23:'+str(i)} |
---|
909 | for ind in [3,4,5,6]: |
---|
910 | at[ind] = parmDict[names[ind]] |
---|
911 | if ind in [3,4,5]: |
---|
912 | name = names[ind].replace('A','dA') |
---|
913 | else: |
---|
914 | name = names[ind] |
---|
915 | if name in sigDict: |
---|
916 | atomsSig[str(i)+':'+str(ind)] = sigDict[name] |
---|
917 | if at[9] == 'I': |
---|
918 | at[10] = parmDict[names[10]] |
---|
919 | if names[10] in sigDict: |
---|
920 | atomsSig[str(i)+':10'] = sigDict[names[10]] |
---|
921 | else: |
---|
922 | for ind in [11,12,13,14,15,16]: |
---|
923 | at[ind] = parmDict[names[ind]] |
---|
924 | if names[ind] in sigDict: |
---|
925 | atomsSig[str(i)+':'+str(ind)] = sigDict[names[ind]] |
---|
926 | PrintAtomsAndSig(General,Atoms,atomsSig) |
---|
927 | |
---|
928 | if 'SH Texture' in General: |
---|
929 | textureData = General['SH Texture'] |
---|
930 | if textureData['Order']: |
---|
931 | SHtextureSig = {} |
---|
932 | for name in ['omega','chi','phi']: |
---|
933 | aname = pfx+'SH '+name |
---|
934 | textureData['Sample '+name][1] = parmDict[aname] |
---|
935 | if aname in sigDict: |
---|
936 | SHtextureSig['Sample '+name] = sigDict[aname] |
---|
937 | for name in textureData['SH Coeff'][1]: |
---|
938 | aname = pfx+name |
---|
939 | textureData['SH Coeff'][1][name] = parmDict[aname] |
---|
940 | if aname in sigDict: |
---|
941 | SHtextureSig[name] = sigDict[aname] |
---|
942 | PrintSHtextureAndSig(textureData,SHtextureSig) |
---|
943 | |
---|
944 | def GetHistogramPhaseData(Phases,Histograms,Print=True): |
---|
945 | |
---|
946 | def PrintSize(hapData): |
---|
947 | if hapData[0] in ['isotropic','uniaxial']: |
---|
948 | line = '\n Size model : %9s'%(hapData[0]) |
---|
949 | line += ' equatorial:'+'%12.3f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
---|
950 | if hapData[0] == 'uniaxial': |
---|
951 | line += ' axial:'+'%12.3f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
---|
952 | print line |
---|
953 | else: |
---|
954 | print '\n Size model : %s'%(hapData[0]) |
---|
955 | Snames = ['S11','S22','S33','S12','S13','S23'] |
---|
956 | ptlbls = ' names :' |
---|
957 | ptstr = ' values:' |
---|
958 | varstr = ' refine:' |
---|
959 | for i,name in enumerate(Snames): |
---|
960 | ptlbls += '%12s' % (name) |
---|
961 | ptstr += '%12.6f' % (hapData[4][i]) |
---|
962 | varstr += '%12s' % (str(hapData[5][i])) |
---|
963 | print ptlbls |
---|
964 | print ptstr |
---|
965 | print varstr |
---|
966 | |
---|
967 | def PrintMuStrain(hapData,SGData): |
---|
968 | if hapData[0] in ['isotropic','uniaxial']: |
---|
969 | line = '\n Mustrain model: %9s'%(hapData[0]) |
---|
970 | line += ' equatorial:'+'%12.1f'%(hapData[1][0])+' Refine? '+str(hapData[2][0]) |
---|
971 | if hapData[0] == 'uniaxial': |
---|
972 | line += ' axial:'+'%12.1f'%(hapData[1][1])+' Refine? '+str(hapData[2][1]) |
---|
973 | print line |
---|
974 | else: |
---|
975 | print '\n Mustrain model: %s'%(hapData[0]) |
---|
976 | Snames = G2spc.MustrainNames(SGData) |
---|
977 | ptlbls = ' names :' |
---|
978 | ptstr = ' values:' |
---|
979 | varstr = ' refine:' |
---|
980 | for i,name in enumerate(Snames): |
---|
981 | ptlbls += '%12s' % (name) |
---|
982 | ptstr += '%12.6f' % (hapData[4][i]) |
---|
983 | varstr += '%12s' % (str(hapData[5][i])) |
---|
984 | print ptlbls |
---|
985 | print ptstr |
---|
986 | print varstr |
---|
987 | |
---|
988 | def PrintHStrain(hapData,SGData): |
---|
989 | print '\n Hydrostatic/elastic strain: ' |
---|
990 | Hsnames = G2spc.HStrainNames(SGData) |
---|
991 | ptlbls = ' names :' |
---|
992 | ptstr = ' values:' |
---|
993 | varstr = ' refine:' |
---|
994 | for i,name in enumerate(Hsnames): |
---|
995 | ptlbls += '%12s' % (name) |
---|
996 | ptstr += '%12.6f' % (hapData[0][i]) |
---|
997 | varstr += '%12s' % (str(hapData[1][i])) |
---|
998 | print ptlbls |
---|
999 | print ptstr |
---|
1000 | print varstr |
---|
1001 | |
---|
1002 | def PrintSHPO(hapData): |
---|
1003 | print '\n Spherical harmonics preferred orientation: Order:' + \ |
---|
1004 | str(hapData[4])+' Refine? '+str(hapData[2]) |
---|
1005 | ptlbls = ' names :' |
---|
1006 | ptstr = ' values:' |
---|
1007 | for item in hapData[5]: |
---|
1008 | ptlbls += '%12s'%(item) |
---|
1009 | ptstr += '%12.3f'%(hapData[5][item]) |
---|
1010 | print ptlbls |
---|
1011 | print ptstr |
---|
1012 | |
---|
1013 | hapDict = {} |
---|
1014 | hapVary = [] |
---|
1015 | controlDict = {} |
---|
1016 | poType = {} |
---|
1017 | poAxes = {} |
---|
1018 | spAxes = {} |
---|
1019 | spType = {} |
---|
1020 | |
---|
1021 | for phase in Phases: |
---|
1022 | HistoPhase = Phases[phase]['Histograms'] |
---|
1023 | SGData = Phases[phase]['General']['SGData'] |
---|
1024 | cell = Phases[phase]['General']['Cell'][1:7] |
---|
1025 | A = G2lat.cell2A(cell) |
---|
1026 | pId = Phases[phase]['pId'] |
---|
1027 | histoList = HistoPhase.keys() |
---|
1028 | histoList.sort() |
---|
1029 | for histogram in histoList: |
---|
1030 | try: |
---|
1031 | Histogram = Histograms[histogram] |
---|
1032 | except KeyError: |
---|
1033 | #skip if histogram not included e.g. in a sequential refinement |
---|
1034 | continue |
---|
1035 | hapData = HistoPhase[histogram] |
---|
1036 | hId = Histogram['hId'] |
---|
1037 | limits = Histogram['Limits'][1] |
---|
1038 | inst = Histogram['Instrument Parameters'] |
---|
1039 | inst = dict(zip(inst[3],inst[1])) |
---|
1040 | Zero = inst['Zero'] |
---|
1041 | if 'C' in inst['Type']: |
---|
1042 | try: |
---|
1043 | wave = inst['Lam'] |
---|
1044 | except KeyError: |
---|
1045 | wave = inst['Lam1'] |
---|
1046 | dmin = wave/(2.0*sind(limits[1]/2.0)) |
---|
1047 | pfx = str(pId)+':'+str(hId)+':' |
---|
1048 | for item in ['Scale','Extinction']: |
---|
1049 | hapDict[pfx+item] = hapData[item][0] |
---|
1050 | if hapData[item][1]: |
---|
1051 | hapVary.append(pfx+item) |
---|
1052 | names = G2spc.HStrainNames(SGData) |
---|
1053 | for i,name in enumerate(names): |
---|
1054 | hapDict[pfx+name] = hapData['HStrain'][0][i] |
---|
1055 | if hapData['HStrain'][1][i]: |
---|
1056 | hapVary.append(pfx+name) |
---|
1057 | controlDict[pfx+'poType'] = hapData['Pref.Ori.'][0] |
---|
1058 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
1059 | hapDict[pfx+'MD'] = hapData['Pref.Ori.'][1] |
---|
1060 | controlDict[pfx+'MDAxis'] = hapData['Pref.Ori.'][3] |
---|
1061 | if hapData['Pref.Ori.'][2]: |
---|
1062 | hapVary.append(pfx+'MD') |
---|
1063 | else: #'SH' spherical harmonics |
---|
1064 | controlDict[pfx+'SHord'] = hapData['Pref.Ori.'][4] |
---|
1065 | controlDict[pfx+'SHncof'] = len(hapData['Pref.Ori.'][5]) |
---|
1066 | for item in hapData['Pref.Ori.'][5]: |
---|
1067 | hapDict[pfx+item] = hapData['Pref.Ori.'][5][item] |
---|
1068 | if hapData['Pref.Ori.'][2]: |
---|
1069 | hapVary.append(pfx+item) |
---|
1070 | for item in ['Mustrain','Size']: |
---|
1071 | controlDict[pfx+item+'Type'] = hapData[item][0] |
---|
1072 | if hapData[item][0] in ['isotropic','uniaxial']: |
---|
1073 | hapDict[pfx+item+':i'] = hapData[item][1][0] |
---|
1074 | if hapData[item][2][0]: |
---|
1075 | hapVary.append(pfx+item+':i') |
---|
1076 | if hapData[item][0] == 'uniaxial': |
---|
1077 | controlDict[pfx+item+'Axis'] = hapData[item][3] |
---|
1078 | hapDict[pfx+item+':a'] = hapData[item][1][1] |
---|
1079 | if hapData[item][2][1]: |
---|
1080 | hapVary.append(pfx+item+':a') |
---|
1081 | else: #generalized for mustrain or ellipsoidal for size |
---|
1082 | if item == 'Mustrain': |
---|
1083 | names = G2spc.MustrainNames(SGData) |
---|
1084 | pwrs = [] |
---|
1085 | for name in names: |
---|
1086 | h,k,l = name[1:] |
---|
1087 | pwrs.append([int(h),int(k),int(l)]) |
---|
1088 | controlDict[pfx+'MuPwrs'] = pwrs |
---|
1089 | for i in range(len(hapData[item][4])): |
---|
1090 | sfx = ':'+str(i) |
---|
1091 | hapDict[pfx+item+sfx] = hapData[item][4][i] |
---|
1092 | if hapData[item][5][i]: |
---|
1093 | hapVary.append(pfx+item+sfx) |
---|
1094 | |
---|
1095 | if Print: |
---|
1096 | print '\n Phase: ',phase,' in histogram: ',histogram |
---|
1097 | print 135*'-' |
---|
1098 | print ' Phase fraction : %10.4f'%(hapData['Scale'][0]),' Refine?',hapData['Scale'][1] |
---|
1099 | print ' Extinction coeff: %10.4f'%(hapData['Extinction'][0]),' Refine?',hapData['Extinction'][1] |
---|
1100 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
1101 | Ax = hapData['Pref.Ori.'][3] |
---|
1102 | print ' March-Dollase PO: %10.4f'%(hapData['Pref.Ori.'][1]),' Refine?',hapData['Pref.Ori.'][2], \ |
---|
1103 | ' Axis: %d %d %d'%(Ax[0],Ax[1],Ax[2]) |
---|
1104 | else: #'SH' for spherical harmonics |
---|
1105 | PrintSHPO(hapData['Pref.Ori.']) |
---|
1106 | PrintSize(hapData['Size']) |
---|
1107 | PrintMuStrain(hapData['Mustrain'],SGData) |
---|
1108 | PrintHStrain(hapData['HStrain'],SGData) |
---|
1109 | HKLd = np.array(G2lat.GenHLaue(dmin,SGData,A)) |
---|
1110 | refList = [] |
---|
1111 | for h,k,l,d in HKLd: |
---|
1112 | ext,mul,Uniq,phi = G2spc.GenHKLf([h,k,l],SGData) |
---|
1113 | if ext: |
---|
1114 | continue |
---|
1115 | if 'C' in inst['Type']: |
---|
1116 | pos = 2.0*asind(wave/(2.0*d))+Zero |
---|
1117 | if limits[0] < pos < limits[1]: |
---|
1118 | refList.append([h,k,l,mul,d,pos,0.0,0.0,0.0,0.0,0.0,Uniq,phi,0.0]) |
---|
1119 | else: |
---|
1120 | raise ValueError |
---|
1121 | Histogram['Reflection Lists'][phase] = refList |
---|
1122 | return hapVary,hapDict,controlDict |
---|
1123 | |
---|
1124 | def SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms,Print=True): |
---|
1125 | |
---|
1126 | def PrintSizeAndSig(hapData,sizeSig): |
---|
1127 | line = '\n Size model: %9s'%(hapData[0]) |
---|
1128 | refine = False |
---|
1129 | if hapData[0] in ['isotropic','uniaxial']: |
---|
1130 | line += ' equatorial:%12.3f'%(hapData[1][0]) |
---|
1131 | if sizeSig[0][0]: |
---|
1132 | line += ', sig: %8.3f'%(sizeSig[0][0]) |
---|
1133 | refine = True |
---|
1134 | if hapData[0] == 'uniaxial': |
---|
1135 | line += ' axial:%12.3f'%(hapData[1][1]) |
---|
1136 | if sizeSig[0][1]: |
---|
1137 | refine = True |
---|
1138 | line += ', sig: %8.3f'%(sizeSig[0][1]) |
---|
1139 | if refine: |
---|
1140 | print line |
---|
1141 | else: |
---|
1142 | Snames = ['S11','S22','S33','S12','S13','S23'] |
---|
1143 | ptlbls = ' name :' |
---|
1144 | ptstr = ' value :' |
---|
1145 | sigstr = ' sig :' |
---|
1146 | for i,name in enumerate(Snames): |
---|
1147 | ptlbls += '%12s' % (name) |
---|
1148 | ptstr += '%12.6f' % (hapData[4][i]) |
---|
1149 | if sizeSig[1][i]: |
---|
1150 | refine = True |
---|
1151 | sigstr += '%12.6f' % (sizeSig[1][i]) |
---|
1152 | else: |
---|
1153 | sigstr += 12*' ' |
---|
1154 | if refine: |
---|
1155 | print line |
---|
1156 | print ptlbls |
---|
1157 | print ptstr |
---|
1158 | print sigstr |
---|
1159 | |
---|
1160 | def PrintMuStrainAndSig(hapData,mustrainSig,SGData): |
---|
1161 | line = '\n Mustrain model: %9s'%(hapData[0]) |
---|
1162 | refine = False |
---|
1163 | if hapData[0] in ['isotropic','uniaxial']: |
---|
1164 | line += ' equatorial:%12.1f'%(hapData[1][0]) |
---|
1165 | if mustrainSig[0][0]: |
---|
1166 | line += ', sig: %8.1f'%(mustrainSig[0][0]) |
---|
1167 | refine = True |
---|
1168 | if hapData[0] == 'uniaxial': |
---|
1169 | line += ' axial:%12.1f'%(hapData[1][1]) |
---|
1170 | if mustrainSig[0][1]: |
---|
1171 | line += ', sig: %8.1f'%(mustrainSig[0][1]) |
---|
1172 | if refine: |
---|
1173 | print line |
---|
1174 | else: |
---|
1175 | Snames = G2spc.MustrainNames(SGData) |
---|
1176 | ptlbls = ' name :' |
---|
1177 | ptstr = ' value :' |
---|
1178 | sigstr = ' sig :' |
---|
1179 | for i,name in enumerate(Snames): |
---|
1180 | ptlbls += '%12s' % (name) |
---|
1181 | ptstr += '%12.6f' % (hapData[4][i]) |
---|
1182 | if mustrainSig[1][i]: |
---|
1183 | refine = True |
---|
1184 | sigstr += '%12.6f' % (mustrainSig[1][i]) |
---|
1185 | else: |
---|
1186 | sigstr += 12*' ' |
---|
1187 | if refine: |
---|
1188 | print line |
---|
1189 | print ptlbls |
---|
1190 | print ptstr |
---|
1191 | print sigstr |
---|
1192 | |
---|
1193 | def PrintHStrainAndSig(hapData,strainSig,SGData): |
---|
1194 | Hsnames = G2spc.HStrainNames(SGData) |
---|
1195 | ptlbls = ' name :' |
---|
1196 | ptstr = ' value :' |
---|
1197 | sigstr = ' sig :' |
---|
1198 | refine = False |
---|
1199 | for i,name in enumerate(Hsnames): |
---|
1200 | ptlbls += '%12s' % (name) |
---|
1201 | ptstr += '%12.6g' % (hapData[0][i]) |
---|
1202 | if name in strainSig: |
---|
1203 | refine = True |
---|
1204 | sigstr += '%12.6g' % (strainSig[name]) |
---|
1205 | else: |
---|
1206 | sigstr += 12*' ' |
---|
1207 | if refine: |
---|
1208 | print '\n Hydrostatic/elastic strain: ' |
---|
1209 | print ptlbls |
---|
1210 | print ptstr |
---|
1211 | print sigstr |
---|
1212 | |
---|
1213 | def PrintSHPOAndSig(hapData,POsig): |
---|
1214 | print '\n Spherical harmonics preferred orientation: Order:'+str(hapData[4]) |
---|
1215 | ptlbls = ' names :' |
---|
1216 | ptstr = ' values:' |
---|
1217 | sigstr = ' sig :' |
---|
1218 | for item in hapData[5]: |
---|
1219 | ptlbls += '%12s'%(item) |
---|
1220 | ptstr += '%12.3f'%(hapData[5][item]) |
---|
1221 | if item in POsig: |
---|
1222 | sigstr += '%12.3f'%(POsig[item]) |
---|
1223 | else: |
---|
1224 | sigstr += 12*' ' |
---|
1225 | print ptlbls |
---|
1226 | print ptstr |
---|
1227 | print sigstr |
---|
1228 | |
---|
1229 | for phase in Phases: |
---|
1230 | HistoPhase = Phases[phase]['Histograms'] |
---|
1231 | SGData = Phases[phase]['General']['SGData'] |
---|
1232 | pId = Phases[phase]['pId'] |
---|
1233 | histoList = HistoPhase.keys() |
---|
1234 | histoList.sort() |
---|
1235 | for histogram in histoList: |
---|
1236 | try: |
---|
1237 | Histogram = Histograms[histogram] |
---|
1238 | except KeyError: |
---|
1239 | #skip if histogram not included e.g. in a sequential refinement |
---|
1240 | continue |
---|
1241 | print '\n Phase: ',phase,' in histogram: ',histogram |
---|
1242 | print 130*'-' |
---|
1243 | hapData = HistoPhase[histogram] |
---|
1244 | hId = Histogram['hId'] |
---|
1245 | pfx = str(pId)+':'+str(hId)+':' |
---|
1246 | print ' Final refinement RF, RF^2 = %.2f%%, %.2f%% on %d reflections' \ |
---|
1247 | %(Histogram[pfx+'Rf'],Histogram[pfx+'Rf^2'],Histogram[pfx+'Nref']) |
---|
1248 | |
---|
1249 | PhFrExtPOSig = {} |
---|
1250 | for item in ['Scale','Extinction']: |
---|
1251 | hapData[item][0] = parmDict[pfx+item] |
---|
1252 | if pfx+item in sigDict: |
---|
1253 | PhFrExtPOSig[item] = sigDict[pfx+item] |
---|
1254 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
1255 | hapData['Pref.Ori.'][1] = parmDict[pfx+'MD'] |
---|
1256 | if pfx+'MD' in sigDict: |
---|
1257 | PhFrExtPOSig['MD'] = sigDict[pfx+'MD'] |
---|
1258 | else: #'SH' spherical harmonics |
---|
1259 | for item in hapData['Pref.Ori.'][5]: |
---|
1260 | hapData['Pref.Ori.'][5][item] = parmDict[pfx+item] |
---|
1261 | if pfx+item in sigDict: |
---|
1262 | PhFrExtPOSig[item] = sigDict[pfx+item] |
---|
1263 | if Print: |
---|
1264 | if 'Scale' in PhFrExtPOSig: |
---|
1265 | print ' Phase fraction : %10.4f, sig %10.4f'%(hapData['Scale'][0],PhFrExtPOSig['Scale']) |
---|
1266 | if 'Extinction' in PhFrExtPOSig: |
---|
1267 | print ' Extinction coeff: %10.4f, sig %10.4f'%(hapData['Extinction'][0],PhFrExtPOSig['Extinction']) |
---|
1268 | if hapData['Pref.Ori.'][0] == 'MD': |
---|
1269 | if 'MD' in PhFrExtPOSig: |
---|
1270 | print ' March-Dollase PO: %10.4f, sig %10.4f'%(hapData['Pref.Ori.'][1],PhFrExtPOSig['MD']) |
---|
1271 | else: |
---|
1272 | PrintSHPOAndSig(hapData['Pref.Ori.'],PhFrExtPOSig) |
---|
1273 | SizeMuStrSig = {'Mustrain':[[0,0],[0 for i in range(len(hapData['Mustrain'][4]))]], |
---|
1274 | 'Size':[[0,0],[0 for i in range(len(hapData['Size'][4]))]], |
---|
1275 | 'HStrain':{}} |
---|
1276 | for item in ['Mustrain','Size']: |
---|
1277 | if hapData[item][0] in ['isotropic','uniaxial']: |
---|
1278 | hapData[item][1][0] = parmDict[pfx+item+':i'] |
---|
1279 | if item == 'Size': |
---|
1280 | hapData[item][1][0] = min(10.,max(0.001,hapData[item][1][0])) |
---|
1281 | if pfx+item+':i' in sigDict: |
---|
1282 | SizeMuStrSig[item][0][0] = sigDict[pfx+item+':i'] |
---|
1283 | if hapData[item][0] == 'uniaxial': |
---|
1284 | hapData[item][1][1] = parmDict[pfx+item+':a'] |
---|
1285 | if item == 'Size': |
---|
1286 | hapData[item][1][1] = min(10.,max(0.001,hapData[item][1][1])) |
---|
1287 | if pfx+item+':a' in sigDict: |
---|
1288 | SizeMuStrSig[item][0][1] = sigDict[pfx+item+':a'] |
---|
1289 | else: #generalized for mustrain or ellipsoidal for size |
---|
1290 | for i in range(len(hapData[item][4])): |
---|
1291 | sfx = ':'+str(i) |
---|
1292 | hapData[item][4][i] = parmDict[pfx+item+sfx] |
---|
1293 | if pfx+item+sfx in sigDict: |
---|
1294 | SizeMuStrSig[item][1][i] = sigDict[pfx+item+sfx] |
---|
1295 | names = G2spc.HStrainNames(SGData) |
---|
1296 | for i,name in enumerate(names): |
---|
1297 | hapData['HStrain'][0][i] = parmDict[pfx+name] |
---|
1298 | if pfx+name in sigDict: |
---|
1299 | SizeMuStrSig['HStrain'][name] = sigDict[pfx+name] |
---|
1300 | if Print: |
---|
1301 | PrintSizeAndSig(hapData['Size'],SizeMuStrSig['Size']) |
---|
1302 | PrintMuStrainAndSig(hapData['Mustrain'],SizeMuStrSig['Mustrain'],SGData) |
---|
1303 | PrintHStrainAndSig(hapData['HStrain'],SizeMuStrSig['HStrain'],SGData) |
---|
1304 | |
---|
1305 | def GetHistogramData(Histograms,Print=True): |
---|
1306 | |
---|
1307 | def GetBackgroundParms(hId,Background): |
---|
1308 | Back = Background[0] |
---|
1309 | Debye = Background[1] |
---|
1310 | bakType,bakFlag = Back[:2] |
---|
1311 | backVals = Back[3:] |
---|
1312 | backNames = [':'+str(hId)+':Back:'+str(i) for i in range(len(backVals))] |
---|
1313 | backDict = dict(zip(backNames,backVals)) |
---|
1314 | backVary = [] |
---|
1315 | if bakFlag: |
---|
1316 | backVary = backNames |
---|
1317 | backDict[':'+str(hId)+':nDebye'] = Debye['nDebye'] |
---|
1318 | debyeDict = {} |
---|
1319 | debyeList = [] |
---|
1320 | for i in range(Debye['nDebye']): |
---|
1321 | debyeNames = [':'+str(hId)+':DebyeA:'+str(i),':'+str(hId)+':DebyeR:'+str(i),':'+str(hId)+':DebyeU:'+str(i)] |
---|
1322 | debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2]))) |
---|
1323 | debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2]) |
---|
1324 | debyeVary = [] |
---|
1325 | for item in debyeList: |
---|
1326 | if item[1]: |
---|
1327 | debyeVary.append(item[0]) |
---|
1328 | backDict.update(debyeDict) |
---|
1329 | backVary += debyeVary |
---|
1330 | return bakType,backDict,backVary |
---|
1331 | |
---|
1332 | def GetInstParms(hId,Inst): |
---|
1333 | insVals,insFlags,insNames = Inst[1:4] |
---|
1334 | dataType = insVals[0] |
---|
1335 | instDict = {} |
---|
1336 | insVary = [] |
---|
1337 | pfx = ':'+str(hId)+':' |
---|
1338 | for i,flag in enumerate(insFlags): |
---|
1339 | insName = pfx+insNames[i] |
---|
1340 | instDict[insName] = insVals[i] |
---|
1341 | if flag: |
---|
1342 | insVary.append(insName) |
---|
1343 | instDict[pfx+'X'] = max(instDict[pfx+'X'],0.001) |
---|
1344 | instDict[pfx+'Y'] = max(instDict[pfx+'Y'],0.001) |
---|
1345 | instDict[pfx+'SH/L'] = max(instDict[pfx+'SH/L'],0.0005) |
---|
1346 | return dataType,instDict,insVary |
---|
1347 | |
---|
1348 | def GetSampleParms(hId,Sample): |
---|
1349 | sampVary = [] |
---|
1350 | hfx = ':'+str(hId)+':' |
---|
1351 | sampDict = {hfx+'Gonio. radius':Sample['Gonio. radius'],hfx+'Omega':Sample['Omega'], |
---|
1352 | hfx+'Chi':Sample['Chi'],hfx+'Phi':Sample['Phi']} |
---|
1353 | Type = Sample['Type'] |
---|
1354 | if 'Bragg' in Type: #Bragg-Brentano |
---|
1355 | for item in ['Scale','Shift','Transparency']: #surface roughness?, diffuse scattering? |
---|
1356 | sampDict[hfx+item] = Sample[item][0] |
---|
1357 | if Sample[item][1]: |
---|
1358 | sampVary.append(hfx+item) |
---|
1359 | elif 'Debye' in Type: #Debye-Scherrer |
---|
1360 | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
---|
1361 | sampDict[hfx+item] = Sample[item][0] |
---|
1362 | if Sample[item][1]: |
---|
1363 | sampVary.append(hfx+item) |
---|
1364 | return Type,sampDict,sampVary |
---|
1365 | |
---|
1366 | def PrintBackground(Background): |
---|
1367 | Back = Background[0] |
---|
1368 | Debye = Background[1] |
---|
1369 | print '\n Background function: ',Back[0],' Refine?',bool(Back[1]) |
---|
1370 | line = ' Coefficients: ' |
---|
1371 | for i,back in enumerate(Back[3:]): |
---|
1372 | line += '%10.3f'%(back) |
---|
1373 | if i and not i%10: |
---|
1374 | line += '\n'+15*' ' |
---|
1375 | print line |
---|
1376 | if Debye['nDebye']: |
---|
1377 | print '\n Debye diffuse scattering coefficients' |
---|
1378 | parms = ['DebyeA','DebyeR','DebyeU'] |
---|
1379 | line = ' names :' |
---|
1380 | for parm in parms: |
---|
1381 | line += '%16s'%(parm) |
---|
1382 | print line |
---|
1383 | for j,term in enumerate(Debye['debyeTerms']): |
---|
1384 | line = ' term'+'%2d'%(j)+':' |
---|
1385 | for i in range(3): |
---|
1386 | line += '%10.4g %5s'%(term[2*i],bool(term[2*i+1])) |
---|
1387 | print line |
---|
1388 | |
---|
1389 | def PrintInstParms(Inst): |
---|
1390 | print '\n Instrument Parameters:' |
---|
1391 | ptlbls = ' name :' |
---|
1392 | ptstr = ' value :' |
---|
1393 | varstr = ' refine:' |
---|
1394 | instNames = Inst[3][1:] |
---|
1395 | for i,name in enumerate(instNames): |
---|
1396 | ptlbls += '%12s' % (name) |
---|
1397 | ptstr += '%12.6f' % (Inst[1][i+1]) |
---|
1398 | if name in ['Lam1','Lam2','Azimuth']: |
---|
1399 | varstr += 12*' ' |
---|
1400 | else: |
---|
1401 | varstr += '%12s' % (str(bool(Inst[2][i+1]))) |
---|
1402 | print ptlbls |
---|
1403 | print ptstr |
---|
1404 | print varstr |
---|
1405 | |
---|
1406 | def PrintSampleParms(Sample): |
---|
1407 | print '\n Sample Parameters:' |
---|
1408 | print ' Goniometer omega = %.2f, chi = %.2f, phi = %.2f'% \ |
---|
1409 | (Sample['Omega'],Sample['Chi'],Sample['Phi']) |
---|
1410 | ptlbls = ' name :' |
---|
1411 | ptstr = ' value :' |
---|
1412 | varstr = ' refine:' |
---|
1413 | if 'Bragg' in Sample['Type']: |
---|
1414 | for item in ['Scale','Shift','Transparency']: |
---|
1415 | ptlbls += '%14s'%(item) |
---|
1416 | ptstr += '%14.4f'%(Sample[item][0]) |
---|
1417 | varstr += '%14s'%(str(bool(Sample[item][1]))) |
---|
1418 | |
---|
1419 | elif 'Debye' in Type: #Debye-Scherrer |
---|
1420 | for item in ['Scale','Absorption','DisplaceX','DisplaceY']: |
---|
1421 | ptlbls += '%14s'%(item) |
---|
1422 | ptstr += '%14.4f'%(Sample[item][0]) |
---|
1423 | varstr += '%14s'%(str(bool(Sample[item][1]))) |
---|
1424 | |
---|
1425 | print ptlbls |
---|
1426 | print ptstr |
---|
1427 | print varstr |
---|
1428 | |
---|
1429 | |
---|
1430 | histDict = {} |
---|
1431 | histVary = [] |
---|
1432 | controlDict = {} |
---|
1433 | histoList = Histograms.keys() |
---|
1434 | histoList.sort() |
---|
1435 | for histogram in histoList: |
---|
1436 | Histogram = Histograms[histogram] |
---|
1437 | hId = Histogram['hId'] |
---|
1438 | pfx = ':'+str(hId)+':' |
---|
1439 | controlDict[pfx+'Limits'] = Histogram['Limits'][1] |
---|
1440 | |
---|
1441 | Background = Histogram['Background'] |
---|
1442 | Type,bakDict,bakVary = GetBackgroundParms(hId,Background) |
---|
1443 | controlDict[pfx+'bakType'] = Type |
---|
1444 | histDict.update(bakDict) |
---|
1445 | histVary += bakVary |
---|
1446 | |
---|
1447 | Inst = Histogram['Instrument Parameters'] |
---|
1448 | Type,instDict,insVary = GetInstParms(hId,Inst) |
---|
1449 | controlDict[pfx+'histType'] = Type |
---|
1450 | if pfx+'Lam1' in instDict: |
---|
1451 | controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam1'] |
---|
1452 | else: |
---|
1453 | controlDict[pfx+'keV'] = 12.397639/instDict[pfx+'Lam'] |
---|
1454 | histDict.update(instDict) |
---|
1455 | histVary += insVary |
---|
1456 | |
---|
1457 | Sample = Histogram['Sample Parameters'] |
---|
1458 | Type,sampDict,sampVary = GetSampleParms(hId,Sample) |
---|
1459 | controlDict[pfx+'instType'] = Type |
---|
1460 | histDict.update(sampDict) |
---|
1461 | histVary += sampVary |
---|
1462 | |
---|
1463 | if Print: |
---|
1464 | print '\n Histogram: ',histogram,' histogram Id: ',hId |
---|
1465 | print 135*'-' |
---|
1466 | Units = {'C':' deg','T':' msec'} |
---|
1467 | units = Units[controlDict[pfx+'histType'][2]] |
---|
1468 | Limits = controlDict[pfx+'Limits'] |
---|
1469 | print ' Instrument type: ',Sample['Type'] |
---|
1470 | print ' Histogram limits: %8.2f%s to %8.2f%s'%(Limits[0],units,Limits[1],units) |
---|
1471 | PrintSampleParms(Sample) |
---|
1472 | PrintInstParms(Inst) |
---|
1473 | PrintBackground(Background) |
---|
1474 | |
---|
1475 | return histVary,histDict,controlDict |
---|
1476 | |
---|
1477 | def SetHistogramData(parmDict,sigDict,Histograms,Print=True): |
---|
1478 | |
---|
1479 | def SetBackgroundParms(pfx,Background,parmDict,sigDict): |
---|
1480 | Back = Background[0] |
---|
1481 | Debye = Background[1] |
---|
1482 | lenBack = len(Back[3:]) |
---|
1483 | backSig = [0 for i in range(lenBack+3*Debye['nDebye'])] |
---|
1484 | for i in range(lenBack): |
---|
1485 | Back[3+i] = parmDict[pfx+'Back:'+str(i)] |
---|
1486 | if pfx+'Back:'+str(i) in sigDict: |
---|
1487 | backSig[i] = sigDict[pfx+'Back:'+str(i)] |
---|
1488 | if Debye['nDebye']: |
---|
1489 | for i in range(Debye['nDebye']): |
---|
1490 | names = [pfx+'DebyeA:'+str(i),pfx+'DebyeR:'+str(i),pfx+'DebyeU:'+str(i)] |
---|
1491 | for j,name in enumerate(names): |
---|
1492 | Debye['debyeTerms'][i][2*j] = parmDict[name] |
---|
1493 | if name in sigDict: |
---|
1494 | backSig[lenBack+3*i+j] = sigDict[name] |
---|
1495 | return backSig |
---|
1496 | |
---|
1497 | def SetInstParms(pfx,Inst,parmDict,sigDict): |
---|
1498 | insVals,insFlags,insNames = Inst[1:4] |
---|
1499 | instSig = [0 for i in range(len(insVals))] |
---|
1500 | for i,flag in enumerate(insFlags): |
---|
1501 | insName = pfx+insNames[i] |
---|
1502 | insVals[i] = parmDict[insName] |
---|
1503 | if insName in sigDict: |
---|
1504 | instSig[i] = sigDict[insName] |
---|
1505 | return instSig |
---|
1506 | |
---|
1507 | def SetSampleParms(pfx,Sample,parmDict,sigDict): |
---|
1508 | if 'Bragg' in Sample['Type']: #Bragg-Brentano |
---|
1509 | sampSig = [0 for i in range(3)] |
---|
1510 | for i,item in enumerate(['Scale','Shift','Transparency']): #surface roughness?, diffuse scattering? |
---|
1511 | Sample[item][0] = parmDict[pfx+item] |
---|
1512 | if pfx+item in sigDict: |
---|
1513 | sampSig[i] = sigDict[pfx+item] |
---|
1514 | elif 'Debye' in Sample['Type']: #Debye-Scherrer |
---|
1515 | sampSig = [0 for i in range(4)] |
---|
1516 | for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']): |
---|
1517 | Sample[item][0] = parmDict[pfx+item] |
---|
1518 | if pfx+item in sigDict: |
---|
1519 | sampSig[i] = sigDict[pfx+item] |
---|
1520 | return sampSig |
---|
1521 | |
---|
1522 | def PrintBackgroundSig(Background,backSig): |
---|
1523 | Back = Background[0] |
---|
1524 | Debye = Background[1] |
---|
1525 | lenBack = len(Back[3:]) |
---|
1526 | valstr = ' value : ' |
---|
1527 | sigstr = ' sig : ' |
---|
1528 | refine = False |
---|
1529 | for i,back in enumerate(Back[3:]): |
---|
1530 | valstr += '%10.4g'%(back) |
---|
1531 | if Back[1]: |
---|
1532 | refine = True |
---|
1533 | sigstr += '%10.4g'%(backSig[i]) |
---|
1534 | else: |
---|
1535 | sigstr += 10*' ' |
---|
1536 | if refine: |
---|
1537 | print '\n Background function: ',Back[0] |
---|
1538 | print valstr |
---|
1539 | print sigstr |
---|
1540 | if Debye['nDebye']: |
---|
1541 | ifAny = False |
---|
1542 | ptfmt = "%12.5f" |
---|
1543 | names = ' names :' |
---|
1544 | ptstr = ' values:' |
---|
1545 | sigstr = ' esds :' |
---|
1546 | for item in sigDict: |
---|
1547 | if 'Debye' in item: |
---|
1548 | ifAny = True |
---|
1549 | names += '%12s'%(item) |
---|
1550 | ptstr += ptfmt%(parmDict[item]) |
---|
1551 | sigstr += ptfmt%(sigDict[item]) |
---|
1552 | if ifAny: |
---|
1553 | print '\n Debye diffuse scattering coefficients' |
---|
1554 | print names |
---|
1555 | print ptstr |
---|
1556 | print sigstr |
---|
1557 | |
---|
1558 | def PrintInstParmsSig(Inst,instSig): |
---|
1559 | ptlbls = ' names :' |
---|
1560 | ptstr = ' value :' |
---|
1561 | sigstr = ' sig :' |
---|
1562 | instNames = Inst[3][1:] |
---|
1563 | refine = False |
---|
1564 | for i,name in enumerate(instNames): |
---|
1565 | ptlbls += '%12s' % (name) |
---|
1566 | ptstr += '%12.6f' % (Inst[1][i+1]) |
---|
1567 | if instSig[i+1]: |
---|
1568 | refine = True |
---|
1569 | sigstr += '%12.6f' % (instSig[i+1]) |
---|
1570 | else: |
---|
1571 | sigstr += 12*' ' |
---|
1572 | if refine: |
---|
1573 | print '\n Instrument Parameters:' |
---|
1574 | print ptlbls |
---|
1575 | print ptstr |
---|
1576 | print sigstr |
---|
1577 | |
---|
1578 | def PrintSampleParmsSig(Sample,sampleSig): |
---|
1579 | ptlbls = ' names :' |
---|
1580 | ptstr = ' values:' |
---|
1581 | sigstr = ' sig :' |
---|
1582 | refine = False |
---|
1583 | if 'Bragg' in Sample['Type']: |
---|
1584 | for i,item in enumerate(['Scale','Shift','Transparency']): |
---|
1585 | ptlbls += '%14s'%(item) |
---|
1586 | ptstr += '%14.4f'%(Sample[item][0]) |
---|
1587 | if sampleSig[i]: |
---|
1588 | refine = True |
---|
1589 | sigstr += '%14.4f'%(sampleSig[i]) |
---|
1590 | else: |
---|
1591 | sigstr += 14*' ' |
---|
1592 | |
---|
1593 | elif 'Debye' in Sample['Type']: #Debye-Scherrer |
---|
1594 | for i,item in enumerate(['Scale','Absorption','DisplaceX','DisplaceY']): |
---|
1595 | ptlbls += '%14s'%(item) |
---|
1596 | ptstr += '%14.4f'%(Sample[item][0]) |
---|
1597 | if sampleSig[i]: |
---|
1598 | refine = True |
---|
1599 | sigstr += '%14.4f'%(sampleSig[i]) |
---|
1600 | else: |
---|
1601 | sigstr += 14*' ' |
---|
1602 | |
---|
1603 | if refine: |
---|
1604 | print '\n Sample Parameters:' |
---|
1605 | print ptlbls |
---|
1606 | print ptstr |
---|
1607 | print sigstr |
---|
1608 | |
---|
1609 | histoList = Histograms.keys() |
---|
1610 | histoList.sort() |
---|
1611 | for histogram in histoList: |
---|
1612 | if 'PWDR' in histogram: |
---|
1613 | Histogram = Histograms[histogram] |
---|
1614 | hId = Histogram['hId'] |
---|
1615 | pfx = ':'+str(hId)+':' |
---|
1616 | Background = Histogram['Background'] |
---|
1617 | backSig = SetBackgroundParms(pfx,Background,parmDict,sigDict) |
---|
1618 | |
---|
1619 | Inst = Histogram['Instrument Parameters'] |
---|
1620 | instSig = SetInstParms(pfx,Inst,parmDict,sigDict) |
---|
1621 | |
---|
1622 | Sample = Histogram['Sample Parameters'] |
---|
1623 | sampSig = SetSampleParms(pfx,Sample,parmDict,sigDict) |
---|
1624 | |
---|
1625 | print '\n Histogram: ',histogram,' histogram Id: ',hId |
---|
1626 | print 135*'-' |
---|
1627 | print ' Final refinement wRp = %.2f%% on %d observations in this histogram'%(Histogram['wRp'],Histogram['Nobs']) |
---|
1628 | if Print: |
---|
1629 | print ' Instrument type: ',Sample['Type'] |
---|
1630 | PrintSampleParmsSig(Sample,sampSig) |
---|
1631 | PrintInstParmsSig(Inst,instSig) |
---|
1632 | PrintBackgroundSig(Background,backSig) |
---|
1633 | |
---|
1634 | def GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict): |
---|
1635 | Natoms = calcControls['Natoms'][pfx] |
---|
1636 | Tdata = Natoms*[' ',] |
---|
1637 | Mdata = np.zeros(Natoms) |
---|
1638 | IAdata = Natoms*[' ',] |
---|
1639 | Fdata = np.zeros(Natoms) |
---|
1640 | FFdata = [] |
---|
1641 | BLdata = [] |
---|
1642 | Xdata = np.zeros((3,Natoms)) |
---|
1643 | dXdata = np.zeros((3,Natoms)) |
---|
1644 | Uisodata = np.zeros(Natoms) |
---|
1645 | Uijdata = np.zeros((6,Natoms)) |
---|
1646 | keys = {'Atype:':Tdata,'Amul:':Mdata,'Afrac:':Fdata,'AI/A:':IAdata, |
---|
1647 | 'dAx:':dXdata[0],'dAy:':dXdata[1],'dAz:':dXdata[2], |
---|
1648 | 'Ax:':Xdata[0],'Ay:':Xdata[1],'Az:':Xdata[2],'AUiso:':Uisodata, |
---|
1649 | 'AU11:':Uijdata[0],'AU22:':Uijdata[1],'AU33:':Uijdata[2], |
---|
1650 | 'AU12:':Uijdata[3],'AU13:':Uijdata[4],'AU23:':Uijdata[5]} |
---|
1651 | for iatm in range(Natoms): |
---|
1652 | for key in keys: |
---|
1653 | parm = pfx+key+str(iatm) |
---|
1654 | if parm in parmDict: |
---|
1655 | keys[key][iatm] = parmDict[parm] |
---|
1656 | FFdata.append(FFtables[Tdata[iatm]]) |
---|
1657 | BLdata.append(BLtables[Tdata[iatm]][1]) |
---|
1658 | return FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata |
---|
1659 | |
---|
1660 | def StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict): |
---|
1661 | ''' Compute structure factors for all h,k,l for phase |
---|
1662 | input: |
---|
1663 | refList: [ref] where each ref = h,k,l,m,d,...,[equiv h,k,l],phase[equiv] |
---|
1664 | G: reciprocal metric tensor |
---|
1665 | pfx: phase id string |
---|
1666 | SGData: space group info. dictionary output from SpcGroup |
---|
1667 | calcControls: |
---|
1668 | ParmDict: |
---|
1669 | puts result F^2 in each ref[8] in refList |
---|
1670 | ''' |
---|
1671 | twopi = 2.0*np.pi |
---|
1672 | twopisq = 2.0*np.pi**2 |
---|
1673 | ast = np.sqrt(np.diag(G)) |
---|
1674 | Mast = twopisq*np.multiply.outer(ast,ast) |
---|
1675 | FFtables = calcControls['FFtables'] |
---|
1676 | BLtables = calcControls['BLtables'] |
---|
1677 | FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict) |
---|
1678 | if 'N' in parmDict[hfx+'Type']: |
---|
1679 | FP,FPP = G2el.BlenRes(BLdata,parmDict[hfx+'Lam']) |
---|
1680 | else: |
---|
1681 | FP = np.array([El[hfx+'FP'] for El in FFdata]) |
---|
1682 | FPP = np.array([El[hfx+'FPP'] for El in FFdata]) |
---|
1683 | maxPos = len(SGData['SGOps']) |
---|
1684 | Uij = np.array(G2lat.U6toUij(Uijdata)) |
---|
1685 | bij = Mast*Uij.T |
---|
1686 | for refl in refList: |
---|
1687 | fbs = np.array([0,0]) |
---|
1688 | H = refl[:3] |
---|
1689 | SQ = 1./(2.*refl[4])**2 |
---|
1690 | if 'N' in parmDict[hfx+'Type']: |
---|
1691 | FF = np.array([El[1] for El in BLdata]) |
---|
1692 | else: #'X' |
---|
1693 | FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata]) |
---|
1694 | SQfactor = 4.0*SQ*twopisq |
---|
1695 | Uniq = refl[11] |
---|
1696 | phi = refl[12] |
---|
1697 | phase = twopi*(np.inner(Uniq,(dXdata.T+Xdata.T))+phi[:,np.newaxis]) |
---|
1698 | sinp = np.sin(phase) |
---|
1699 | cosp = np.cos(phase) |
---|
1700 | occ = Mdata*Fdata/len(Uniq) |
---|
1701 | biso = -SQfactor*Uisodata |
---|
1702 | Tiso = np.where(biso<1.,np.exp(biso),1.0) |
---|
1703 | HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq]) |
---|
1704 | Tuij = np.where(HbH<1.,np.exp(HbH),1.0) |
---|
1705 | Tcorr = Tiso*Tuij |
---|
1706 | fa = np.array([(FF+FP)*occ*cosp*Tcorr,-FPP*occ*sinp*Tcorr]) |
---|
1707 | fas = np.sum(np.sum(fa,axis=1),axis=1) #real |
---|
1708 | if not SGData['SGInv']: |
---|
1709 | fb = np.array([(FF+FP)*occ*sinp*Tcorr,FPP*occ*cosp*Tcorr]) |
---|
1710 | fbs = np.sum(np.sum(fb,axis=1),axis=1) |
---|
1711 | fasq = fas**2 |
---|
1712 | fbsq = fbs**2 #imaginary |
---|
1713 | refl[9] = np.sum(fasq)+np.sum(fbsq) |
---|
1714 | refl[10] = atan2d(fbs[0],fas[0]) |
---|
1715 | return refList |
---|
1716 | |
---|
1717 | def StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict): |
---|
1718 | twopi = 2.0*np.pi |
---|
1719 | twopisq = 2.0*np.pi**2 |
---|
1720 | ast = np.sqrt(np.diag(G)) |
---|
1721 | Mast = twopisq*np.multiply.outer(ast,ast) |
---|
1722 | FFtables = calcControls['FFtables'] |
---|
1723 | BLtables = calcControls['BLtables'] |
---|
1724 | FFdata,BLdata,Mdata,Fdata,Xdata,dXdata,IAdata,Uisodata,Uijdata = GetAtomFXU(pfx,FFtables,BLtables,calcControls,parmDict) |
---|
1725 | if 'N' in parmDict[hfx+'Type']: |
---|
1726 | FP = 0. |
---|
1727 | FPP = 0. |
---|
1728 | else: |
---|
1729 | FP = np.array([El[hfx+'FP'] for El in FFdata]) |
---|
1730 | FPP = np.array([El[hfx+'FPP'] for El in FFdata]) |
---|
1731 | maxPos = len(SGData['SGOps']) |
---|
1732 | Uij = np.array(G2lat.U6toUij(Uijdata)) |
---|
1733 | bij = Mast*Uij.T |
---|
1734 | dFdvDict = {} |
---|
1735 | dFdfr = np.zeros((len(refList),len(Mdata))) |
---|
1736 | dFdx = np.zeros((len(refList),len(Mdata),3)) |
---|
1737 | dFdui = np.zeros((len(refList),len(Mdata))) |
---|
1738 | dFdua = np.zeros((len(refList),len(Mdata),6)) |
---|
1739 | for iref,refl in enumerate(refList): |
---|
1740 | H = np.array(refl[:3]) |
---|
1741 | SQ = 1./(2.*refl[4])**2 # or (sin(theta)/lambda)**2 |
---|
1742 | if 'N' in parmDict[hfx+'Type']: |
---|
1743 | FF = np.array([El[1] for El in BLdata]) |
---|
1744 | else: #'X' |
---|
1745 | FF = np.array([G2el.ScatFac(El,SQ)[0] for El in FFdata]) |
---|
1746 | SQfactor = 8.0*SQ*np.pi**2 |
---|
1747 | Uniq = refl[11] |
---|
1748 | phi = refl[12] |
---|
1749 | phase = twopi*(np.inner((dXdata.T+Xdata.T),Uniq)+phi[np.newaxis,:]) |
---|
1750 | sinp = np.sin(phase) |
---|
1751 | cosp = np.cos(phase) |
---|
1752 | occ = Mdata*Fdata/len(Uniq) |
---|
1753 | biso = -SQfactor*Uisodata |
---|
1754 | Tiso = np.where(biso<1.,np.exp(biso),1.0) |
---|
1755 | # HbH = np.array([-np.inner(h,np.inner(bij,h)) for h in Uniq]) |
---|
1756 | HbH = -np.inner(H,np.inner(bij,H)) |
---|
1757 | Hij = np.array([Mast*np.multiply.outer(U,U) for U in Uniq]) |
---|
1758 | Hij = np.array([G2lat.UijtoU6(Uij) for Uij in Hij]) |
---|
1759 | Tuij = np.where(HbH<1.,np.exp(HbH),1.0) |
---|
1760 | Tcorr = Tiso*Tuij |
---|
1761 | fot = (FF+FP)*occ*Tcorr |
---|
1762 | fotp = FPP*occ*Tcorr |
---|
1763 | fa = np.array([fot[:,np.newaxis]*cosp,fotp[:,np.newaxis]*cosp]) #non positions |
---|
1764 | fb = np.array([fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp]) |
---|
1765 | |
---|
1766 | fas = np.sum(np.sum(fa,axis=1),axis=1) |
---|
1767 | fbs = np.sum(np.sum(fb,axis=1),axis=1) |
---|
1768 | fax = np.array([-fot[:,np.newaxis]*sinp,-fotp[:,np.newaxis]*sinp]) #positions |
---|
1769 | fbx = np.array([fot[:,np.newaxis]*cosp,-fot[:,np.newaxis]*cosp]) |
---|
1770 | #sum below is over Uniq |
---|
1771 | dfadfr = np.sum(fa/occ[:,np.newaxis],axis=2) |
---|
1772 | dfadx = np.sum(twopi*Uniq*fax[:,:,:,np.newaxis],axis=2) |
---|
1773 | dfadui = np.sum(-SQfactor*fa,axis=2) |
---|
1774 | dfadua = np.sum(-Hij*fa[:,:,:,np.newaxis],axis=2) |
---|
1775 | #NB: the above have been checked against PA(1:10,1:2) in strfctr.for |
---|
1776 | dFdfr[iref] = 2.*(fas[0]*dfadfr[0]+fas[1]*dfadfr[1])*Mdata/len(Uniq) |
---|
1777 | dFdx[iref] = 2.*(fas[0]*dfadx[0]+fas[1]*dfadx[1]) |
---|
1778 | dFdui[iref] = 2.*(fas[0]*dfadui[0]+fas[1]*dfadui[1]) |
---|
1779 | dFdua[iref] = 2.*(fas[0]*dfadua[0]+fas[1]*dfadua[1]) |
---|
1780 | if not SGData['SGInv']: |
---|
1781 | dfbdfr = np.sum(fb/occ[:,np.newaxis],axis=2) #problem here if occ=0 for some atom |
---|
1782 | dfbdx = np.sum(twopi*Uniq*fbx[:,:,:,np.newaxis],axis=2) |
---|
1783 | dfbdui = np.sum(-SQfactor*fb,axis=2) |
---|
1784 | dfbdua = np.sum(-Hij*fb[:,:,:,np.newaxis],axis=2) |
---|
1785 | dFdfr[iref] += 2.*(fbs[0]*dfbdfr[0]-fbs[1]*dfbdfr[1])*Mdata/len(Uniq) |
---|
1786 | dFdx[iref] += 2.*(fbs[0]*dfbdx[0]+fbs[1]*dfbdx[1]) |
---|
1787 | dFdui[iref] += 2.*(fbs[0]*dfbdui[0]-fbs[1]*dfbdui[1]) |
---|
1788 | dFdua[iref] += 2.*(fbs[0]*dfbdua[0]+fbs[1]*dfbdua[1]) |
---|
1789 | #loop over atoms - each dict entry is list of derivatives for all the reflections |
---|
1790 | for i in range(len(Mdata)): |
---|
1791 | dFdvDict[pfx+'Afrac:'+str(i)] = dFdfr.T[i] |
---|
1792 | dFdvDict[pfx+'dAx:'+str(i)] = dFdx.T[0][i] |
---|
1793 | dFdvDict[pfx+'dAy:'+str(i)] = dFdx.T[1][i] |
---|
1794 | dFdvDict[pfx+'dAz:'+str(i)] = dFdx.T[2][i] |
---|
1795 | dFdvDict[pfx+'AUiso:'+str(i)] = dFdui.T[i] |
---|
1796 | dFdvDict[pfx+'AU11:'+str(i)] = dFdua.T[0][i] |
---|
1797 | dFdvDict[pfx+'AU22:'+str(i)] = dFdua.T[1][i] |
---|
1798 | dFdvDict[pfx+'AU33:'+str(i)] = dFdua.T[2][i] |
---|
1799 | dFdvDict[pfx+'AU12:'+str(i)] = 2.*dFdua.T[3][i] |
---|
1800 | dFdvDict[pfx+'AU13:'+str(i)] = 2.*dFdua.T[4][i] |
---|
1801 | dFdvDict[pfx+'AU23:'+str(i)] = 2.*dFdua.T[5][i] |
---|
1802 | return dFdvDict |
---|
1803 | |
---|
1804 | def Dict2Values(parmdict, varylist): |
---|
1805 | '''Use before call to leastsq to setup list of values for the parameters |
---|
1806 | in parmdict, as selected by key in varylist''' |
---|
1807 | return [parmdict[key] for key in varylist] |
---|
1808 | |
---|
1809 | def Values2Dict(parmdict, varylist, values): |
---|
1810 | ''' Use after call to leastsq to update the parameter dictionary with |
---|
1811 | values corresponding to keys in varylist''' |
---|
1812 | parmdict.update(zip(varylist,values)) |
---|
1813 | |
---|
1814 | def GetNewCellParms(parmDict,varyList): |
---|
1815 | newCellDict = {} |
---|
1816 | Ddict = dict(zip(['D11','D22','D33','D12','D13','D23'],['A'+str(i) for i in range(6)])) |
---|
1817 | for item in varyList: |
---|
1818 | keys = item.split(':') |
---|
1819 | if keys[2] in Ddict: |
---|
1820 | key = keys[0]+'::'+Ddict[keys[2]] |
---|
1821 | parm = keys[0]+'::'+keys[2] |
---|
1822 | newCellDict[parm] = [key,parmDict[key]+parmDict[item]] |
---|
1823 | return newCellDict |
---|
1824 | |
---|
1825 | def ApplyXYZshifts(parmDict,varyList): |
---|
1826 | ''' takes atom x,y,z shift and applies it to corresponding atom x,y,z value |
---|
1827 | input: |
---|
1828 | parmDict - parameter dictionary |
---|
1829 | varyList - list of variables |
---|
1830 | returns: |
---|
1831 | newAtomDict - dictitemionary of new atomic coordinate names & values; |
---|
1832 | key is parameter shift name |
---|
1833 | ''' |
---|
1834 | newAtomDict = {} |
---|
1835 | for item in parmDict: |
---|
1836 | if 'dA' in item: |
---|
1837 | parm = ''.join(item.split('d')) |
---|
1838 | parmDict[parm] += parmDict[item] |
---|
1839 | newAtomDict[item] = [parm,parmDict[parm]] |
---|
1840 | return newAtomDict |
---|
1841 | |
---|
1842 | def SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict): |
---|
1843 | IFCoup = 'Bragg' in calcControls[hfx+'instType'] |
---|
1844 | odfCor = 1.0 |
---|
1845 | H = refl[:3] |
---|
1846 | cell = G2lat.Gmat2cell(g) |
---|
1847 | Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']] |
---|
1848 | Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] |
---|
1849 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
---|
1850 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) #ignore 2 sets of angle derivs. |
---|
1851 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) |
---|
1852 | for item in SHnames: |
---|
1853 | L,M,N = eval(item.strip('C')) |
---|
1854 | Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta) |
---|
1855 | Ksl,x,x = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam) |
---|
1856 | Lnorm = G2lat.Lnorm(L) |
---|
1857 | odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl |
---|
1858 | return odfCor |
---|
1859 | |
---|
1860 | def SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict): |
---|
1861 | FORPI = 12.5663706143592 |
---|
1862 | IFCoup = 'Bragg' in calcControls[hfx+'instType'] |
---|
1863 | odfCor = 1.0 |
---|
1864 | dFdODF = {} |
---|
1865 | dFdSA = [0,0,0] |
---|
1866 | H = refl[:3] |
---|
1867 | cell = G2lat.Gmat2cell(g) |
---|
1868 | Sangls = [parmDict[pfx+'SH omega'],parmDict[pfx+'SH chi'],parmDict[pfx+'SH phi']] |
---|
1869 | Gangls = [parmDict[hfx+'Omega'],parmDict[hfx+'Chi'],parmDict[hfx+'Phi'],parmDict[hfx+'Azimuth']] |
---|
1870 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
---|
1871 | psi,gam,dPSdA,dGMdA = G2lat.SamAng(refl[5]/2.,Gangls,Sangls,IFCoup) |
---|
1872 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],parmDict[pfx+'SHmodel'],parmDict[pfx+'SHorder']) |
---|
1873 | for item in SHnames: |
---|
1874 | L,M,N = eval(item.strip('C')) |
---|
1875 | Kcl = G2lat.GetKcl(L,N,SGData['SGLaue'],phi,beta) |
---|
1876 | Ksl,dKsdp,dKsdg = G2lat.GetKsl(L,M,parmDict[pfx+'SHmodel'],psi,gam) |
---|
1877 | Lnorm = G2lat.Lnorm(L) |
---|
1878 | odfCor += parmDict[pfx+item]*Lnorm*Kcl*Ksl |
---|
1879 | dFdODF[pfx+item] = Lnorm*Kcl*Ksl |
---|
1880 | for i in range(3): |
---|
1881 | dFdSA[i] += parmDict[pfx+item]*Lnorm*Kcl*(dKsdp*dPSdA[i]+dKsdg*dGMdA[i]) |
---|
1882 | return odfCor,dFdODF,dFdSA |
---|
1883 | |
---|
1884 | def SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict): |
---|
1885 | odfCor = 1.0 |
---|
1886 | H = refl[:3] |
---|
1887 | cell = G2lat.Gmat2cell(g) |
---|
1888 | Sangl = [0.,0.,0.] |
---|
1889 | if 'Bragg' in calcControls[hfx+'instType']: |
---|
1890 | Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']] |
---|
1891 | IFCoup = True |
---|
1892 | else: |
---|
1893 | Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']] |
---|
1894 | IFCoup = False |
---|
1895 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
---|
1896 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. |
---|
1897 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) |
---|
1898 | for item in SHnames: |
---|
1899 | L,N = eval(item.strip('C')) |
---|
1900 | Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) |
---|
1901 | odfCor += parmDict[phfx+item]*Lnorm*Kcsl |
---|
1902 | return odfCor |
---|
1903 | |
---|
1904 | def SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict): |
---|
1905 | FORPI = 12.5663706143592 |
---|
1906 | odfCor = 1.0 |
---|
1907 | dFdODF = {} |
---|
1908 | H = refl[:3] |
---|
1909 | cell = G2lat.Gmat2cell(g) |
---|
1910 | Sangl = [0.,0.,0.] |
---|
1911 | if 'Bragg' in calcControls[hfx+'instType']: |
---|
1912 | Gangls = [0.,90.,0.,parmDict[hfx+'Azimuth']] |
---|
1913 | IFCoup = True |
---|
1914 | else: |
---|
1915 | Gangls = [0.,0.,0.,parmDict[hfx+'Azimuth']] |
---|
1916 | IFCoup = False |
---|
1917 | phi,beta = G2lat.CrsAng(H,cell,SGData) |
---|
1918 | psi,gam,x,x = G2lat.SamAng(refl[5]/2.,Gangls,Sangl,IFCoup) #ignore 2 sets of angle derivs. |
---|
1919 | SHnames = G2lat.GenSHCoeff(SGData['SGLaue'],'0',calcControls[phfx+'SHord'],False) |
---|
1920 | for item in SHnames: |
---|
1921 | L,N = eval(item.strip('C')) |
---|
1922 | Kcsl,Lnorm = G2lat.GetKclKsl(L,N,SGData['SGLaue'],psi,phi,beta) |
---|
1923 | odfCor += parmDict[phfx+item]*Lnorm*Kcsl |
---|
1924 | dFdODF[phfx+item] = Kcsl*Lnorm |
---|
1925 | return odfCor,dFdODF |
---|
1926 | |
---|
1927 | def GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): |
---|
1928 | if calcControls[phfx+'poType'] == 'MD': |
---|
1929 | MD = parmDict[phfx+'MD'] |
---|
1930 | MDAxis = calcControls[phfx+'MDAxis'] |
---|
1931 | sumMD = 0 |
---|
1932 | for H in refl[11]: |
---|
1933 | cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) |
---|
1934 | A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) |
---|
1935 | sumMD += A**3 |
---|
1936 | POcorr = sumMD/len(refl[11]) |
---|
1937 | else: #spherical harmonics |
---|
1938 | POcorr = SHPOcal(refl,g,phfx,hfx,SGData,calcControls,parmDict) |
---|
1939 | return POcorr |
---|
1940 | |
---|
1941 | def GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict): |
---|
1942 | POderv = {} |
---|
1943 | if calcControls[phfx+'poType'] == 'MD': |
---|
1944 | MD = parmDict[phfx+'MD'] |
---|
1945 | MDAxis = calcControls[phfx+'MDAxis'] |
---|
1946 | sumMD = 0 |
---|
1947 | sumdMD = 0 |
---|
1948 | for H in refl[11]: |
---|
1949 | cosP,sinP = G2lat.CosSinAngle(H,MDAxis,G) |
---|
1950 | A = 1.0/np.sqrt((MD*cosP)**2+sinP**2/MD) |
---|
1951 | sumMD += A**3 |
---|
1952 | sumdMD -= (1.5*A**5)*(2.0*MD*cosP**2-(sinP/MD)**2) |
---|
1953 | POcorr = sumMD/len(refl[11]) |
---|
1954 | POderv[phfx+'MD'] = sumdMD/len(refl[11]) |
---|
1955 | else: #spherical harmonics |
---|
1956 | POcorr,POderv = SHPOcalDerv(refl,g,phfx,hfx,SGData,calcControls,parmDict) |
---|
1957 | return POcorr,POderv |
---|
1958 | |
---|
1959 | def GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): |
---|
1960 | Icorr = parmDict[phfx+'Scale']*parmDict[hfx+'Scale']*refl[3] #scale*multiplicity |
---|
1961 | if 'X' in parmDict[hfx+'Type']: |
---|
1962 | Icorr *= G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth'])[0] |
---|
1963 | Icorr *= GetPrefOri(refl,G,g,phfx,hfx,SGData,calcControls,parmDict) |
---|
1964 | if pfx+'SHorder' in parmDict: |
---|
1965 | Icorr *= SHTXcal(refl,g,pfx,hfx,SGData,calcControls,parmDict) |
---|
1966 | refl[13] = Icorr |
---|
1967 | |
---|
1968 | def GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict): |
---|
1969 | dIdsh = 1./parmDict[hfx+'Scale'] |
---|
1970 | dIdsp = 1./parmDict[phfx+'Scale'] |
---|
1971 | if 'X' in parmDict[hfx+'Type']: |
---|
1972 | pola,dIdPola = G2pwd.Polarization(parmDict[hfx+'Polariz.'],refl[5],parmDict[hfx+'Azimuth']) |
---|
1973 | dIdPola /= pola |
---|
1974 | else: #'N' |
---|
1975 | dIdPola = 0.0 |
---|
1976 | POcorr,dIdPO = GetPrefOriDerv(refl,G,g,phfx,hfx,SGData,calcControls,parmDict) |
---|
1977 | for iPO in dIdPO: |
---|
1978 | dIdPO[iPO] /= POcorr |
---|
1979 | dFdODF = {} |
---|
1980 | dFdSA = [0,0,0] |
---|
1981 | if pfx+'SHorder' in parmDict: |
---|
1982 | odfCor,dFdODF,dFdSA = SHTXcalDerv(refl,g,pfx,hfx,SGData,calcControls,parmDict) |
---|
1983 | for iSH in dFdODF: |
---|
1984 | dFdODF[iSH] /= odfCor |
---|
1985 | for i in range(3): |
---|
1986 | dFdSA[i] /= odfCor |
---|
1987 | return dIdsh,dIdsp,dIdPola,dIdPO,dFdODF,dFdSA |
---|
1988 | |
---|
1989 | def GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict): |
---|
1990 | costh = cosd(refl[5]/2.) |
---|
1991 | #crystallite size |
---|
1992 | if calcControls[phfx+'SizeType'] == 'isotropic': |
---|
1993 | gam = 1.8*wave/(np.pi*parmDict[phfx+'Size:i']*costh) |
---|
1994 | elif calcControls[phfx+'SizeType'] == 'uniaxial': |
---|
1995 | H = np.array(refl[:3]) |
---|
1996 | P = np.array(calcControls[phfx+'SizeAxis']) |
---|
1997 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
---|
1998 | gam = (1.8*wave/np.pi)/(parmDict[phfx+'Size:i']*parmDict[phfx+'Size:a']*costh) |
---|
1999 | gam *= np.sqrt((sinP*parmDict[phfx+'Size:a'])**2+(cosP*parmDict[phfx+'Size:i'])**2) |
---|
2000 | else: #ellipsoidal crystallites |
---|
2001 | Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] |
---|
2002 | H = np.array(refl[:3]) |
---|
2003 | lenR = G2pwd.ellipseSize(H,Sij,GB) |
---|
2004 | gam = 1.8*wave/(np.pi*costh*lenR) |
---|
2005 | #microstrain |
---|
2006 | if calcControls[phfx+'MustrainType'] == 'isotropic': |
---|
2007 | gam += 0.018*parmDict[phfx+'Mustrain:i']*tand(refl[5]/2.)/np.pi |
---|
2008 | elif calcControls[phfx+'MustrainType'] == 'uniaxial': |
---|
2009 | H = np.array(refl[:3]) |
---|
2010 | P = np.array(calcControls[phfx+'MustrainAxis']) |
---|
2011 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
---|
2012 | Si = parmDict[phfx+'Mustrain:i'] |
---|
2013 | Sa = parmDict[phfx+'Mustrain:a'] |
---|
2014 | gam += 0.018*Si*Sa*tand(refl[5]/2.)/(np.pi*np.sqrt((Si*cosP)**2+(Sa*sinP)**2)) |
---|
2015 | else: #generalized - P.W. Stephens model |
---|
2016 | pwrs = calcControls[phfx+'MuPwrs'] |
---|
2017 | sum = 0 |
---|
2018 | for i,pwr in enumerate(pwrs): |
---|
2019 | sum += parmDict[phfx+'Mustrain:'+str(i)]*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] |
---|
2020 | gam += 0.018*refl[4]**2*tand(refl[5]/2.)*sum |
---|
2021 | return gam |
---|
2022 | |
---|
2023 | def GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict): |
---|
2024 | gamDict = {} |
---|
2025 | costh = cosd(refl[5]/2.) |
---|
2026 | tanth = tand(refl[5]/2.) |
---|
2027 | #crystallite size derivatives |
---|
2028 | if calcControls[phfx+'SizeType'] == 'isotropic': |
---|
2029 | gamDict[phfx+'Size:i'] = -1.80*wave/(np.pi*costh) |
---|
2030 | elif calcControls[phfx+'SizeType'] == 'uniaxial': |
---|
2031 | H = np.array(refl[:3]) |
---|
2032 | P = np.array(calcControls[phfx+'SizeAxis']) |
---|
2033 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
---|
2034 | Si = parmDict[phfx+'Size:i'] |
---|
2035 | Sa = parmDict[phfx+'Size:a'] |
---|
2036 | gami = (1.8*wave/np.pi)/(Si*Sa) |
---|
2037 | sqtrm = np.sqrt((sinP*Sa)**2+(cosP*Si)**2) |
---|
2038 | gam = gami*sqtrm/costh |
---|
2039 | gamDict[phfx+'Size:i'] = gami*Si*cosP**2/(sqtrm*costh)-gam/Si |
---|
2040 | gamDict[phfx+'Size:a'] = gami*Sa*sinP**2/(sqtrm*costh)-gam/Sa |
---|
2041 | else: #ellipsoidal crystallites |
---|
2042 | const = 1.8*wave/(np.pi*costh) |
---|
2043 | Sij =[parmDict[phfx+'Size:%d'%(i)] for i in range(6)] |
---|
2044 | H = np.array(refl[:3]) |
---|
2045 | R,dRdS = G2pwd.ellipseSizeDerv(H,Sij,GB) |
---|
2046 | for i,item in enumerate([phfx+'Size:%d'%(j) for j in range(6)]): |
---|
2047 | gamDict[item] = -(const/R**2)*dRdS[i] |
---|
2048 | #microstrain derivatives |
---|
2049 | if calcControls[phfx+'MustrainType'] == 'isotropic': |
---|
2050 | gamDict[phfx+'Mustrain:i'] = 0.018*tanth/np.pi |
---|
2051 | elif calcControls[phfx+'MustrainType'] == 'uniaxial': |
---|
2052 | H = np.array(refl[:3]) |
---|
2053 | P = np.array(calcControls[phfx+'MustrainAxis']) |
---|
2054 | cosP,sinP = G2lat.CosSinAngle(H,P,G) |
---|
2055 | Si = parmDict[phfx+'Mustrain:i'] |
---|
2056 | Sa = parmDict[phfx+'Mustrain:a'] |
---|
2057 | gami = 0.018*Si*Sa*tanth/np.pi |
---|
2058 | sqtrm = np.sqrt((Si*cosP)**2+(Sa*sinP)**2) |
---|
2059 | gam = gami/sqtrm |
---|
2060 | gamDict[phfx+'Mustrain:i'] = gam/Si-gami*Si*cosP**2/sqtrm**3 |
---|
2061 | gamDict[phfx+'Mustrain:a'] = gam/Sa-gami*Sa*sinP**2/sqtrm**3 |
---|
2062 | else: #generalized - P.W. Stephens model |
---|
2063 | pwrs = calcControls[phfx+'MuPwrs'] |
---|
2064 | const = 0.018*refl[4]**2*tanth |
---|
2065 | for i,pwr in enumerate(pwrs): |
---|
2066 | gamDict[phfx+'Mustrain:'+str(i)] = const*refl[0]**pwr[0]*refl[1]**pwr[1]*refl[2]**pwr[2] |
---|
2067 | return gamDict |
---|
2068 | |
---|
2069 | def GetReflPos(refl,wave,G,hfx,calcControls,parmDict): |
---|
2070 | h,k,l = refl[:3] |
---|
2071 | dsq = 1./G2lat.calc_rDsq2(np.array([h,k,l]),G) |
---|
2072 | d = np.sqrt(dsq) |
---|
2073 | |
---|
2074 | refl[4] = d |
---|
2075 | pos = 2.0*asind(wave/(2.0*d))+parmDict[hfx+'Zero'] |
---|
2076 | const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns |
---|
2077 | if 'Bragg' in calcControls[hfx+'instType']: |
---|
2078 | pos -= const*(4.*parmDict[hfx+'Shift']*cosd(pos/2.0)+ \ |
---|
2079 | parmDict[hfx+'Transparency']*sind(pos)*100.0) #trans(=1/mueff) in cm |
---|
2080 | else: #Debye-Scherrer - simple but maybe not right |
---|
2081 | pos -= const*(parmDict[hfx+'DisplaceX']*cosd(pos)+parmDict[hfx+'DisplaceY']*sind(pos)) |
---|
2082 | return pos |
---|
2083 | |
---|
2084 | def GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict): |
---|
2085 | dpr = 180./np.pi |
---|
2086 | h,k,l = refl[:3] |
---|
2087 | dstsq = G2lat.calc_rDsq(np.array([h,k,l]),A) |
---|
2088 | dst = np.sqrt(dstsq) |
---|
2089 | pos = refl[5]-parmDict[hfx+'Zero'] |
---|
2090 | const = dpr/np.sqrt(1.0-wave**2*dstsq/4.0) |
---|
2091 | dpdw = const*dst |
---|
2092 | dpdA = np.array([h**2,k**2,l**2,h*k,h*l,k*l]) |
---|
2093 | dpdA *= const*wave/(2.0*dst) |
---|
2094 | dpdZ = 1.0 |
---|
2095 | const = 9.e-2/(np.pi*parmDict[hfx+'Gonio. radius']) #shifts in microns |
---|
2096 | if 'Bragg' in calcControls[hfx+'instType']: |
---|
2097 | dpdSh = -4.*const*cosd(pos/2.0) |
---|
2098 | dpdTr = -const*sind(pos)*100.0 |
---|
2099 | return dpdA,dpdw,dpdZ,dpdSh,dpdTr,0.,0. |
---|
2100 | else: #Debye-Scherrer - simple but maybe not right |
---|
2101 | dpdXd = -const*cosd(pos) |
---|
2102 | dpdYd = -const*sind(pos) |
---|
2103 | return dpdA,dpdw,dpdZ,0.,0.,dpdXd,dpdYd |
---|
2104 | |
---|
2105 | def GetHStrainShift(refl,SGData,phfx,parmDict): |
---|
2106 | laue = SGData['SGLaue'] |
---|
2107 | uniq = SGData['SGUniq'] |
---|
2108 | h,k,l = refl[:3] |
---|
2109 | if laue in ['m3','m3m']: |
---|
2110 | Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+ \ |
---|
2111 | refl[4]**2*parmDict[phfx+'eA']*((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2 |
---|
2112 | elif laue in ['6/m','6/mmm','3m1','31m','3']: |
---|
2113 | Dij = parmDict[phfx+'D11']*(h**2+k**2+h*k)+parmDict[phfx+'D33']*l**2 |
---|
2114 | elif laue in ['3R','3mR']: |
---|
2115 | Dij = parmDict[phfx+'D11']*(h**2+k**2+l**2)+parmDict[phfx+'D12']*(h*k+h*l+k*l) |
---|
2116 | elif laue in ['4/m','4/mmm']: |
---|
2117 | Dij = parmDict[phfx+'D11']*(h**2+k**2)+parmDict[phfx+'D33']*l**2 |
---|
2118 | elif laue in ['mmm']: |
---|
2119 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2 |
---|
2120 | elif laue in ['2/m']: |
---|
2121 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2 |
---|
2122 | if uniq == 'a': |
---|
2123 | Dij += parmDict[phfx+'D23']*k*l |
---|
2124 | elif uniq == 'b': |
---|
2125 | Dij += parmDict[phfx+'D13']*h*l |
---|
2126 | elif uniq == 'c': |
---|
2127 | Dij += parmDict[phfx+'D12']*h*k |
---|
2128 | else: |
---|
2129 | Dij = parmDict[phfx+'D11']*h**2+parmDict[phfx+'D22']*k**2+parmDict[phfx+'D33']*l**2+ \ |
---|
2130 | parmDict[phfx+'D12']*h*k+parmDict[phfx+'D13']*h*l+parmDict[phfx+'D23']*k*l |
---|
2131 | return Dij*refl[4]**2*tand(refl[5]/2.0) |
---|
2132 | |
---|
2133 | def GetHStrainShiftDerv(refl,SGData,phfx): |
---|
2134 | laue = SGData['SGLaue'] |
---|
2135 | uniq = SGData['SGUniq'] |
---|
2136 | h,k,l = refl[:3] |
---|
2137 | if laue in ['m3','m3m']: |
---|
2138 | dDijDict = {phfx+'D11':h**2+k**2+l**2, |
---|
2139 | phfx+'eA':((h*k)**2+(h*l)**2+(k*l)**2)/(h**2+k**2+l**2)**2} |
---|
2140 | elif laue in ['6/m','6/mmm','3m1','31m','3']: |
---|
2141 | dDijDict = {phfx+'D11':h**2+k**2+h*k,phfx+'D33':l**2} |
---|
2142 | elif laue in ['3R','3mR']: |
---|
2143 | dDijDict = {phfx+'D11':h**2+k**2+l**2,phfx+'D12':h*k+h*l+k*l} |
---|
2144 | elif laue in ['4/m','4/mmm']: |
---|
2145 | dDijDict = {phfx+'D11':h**2+k**2,phfx+'D33':l**2} |
---|
2146 | elif laue in ['mmm']: |
---|
2147 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2} |
---|
2148 | elif laue in ['2/m']: |
---|
2149 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2} |
---|
2150 | if uniq == 'a': |
---|
2151 | dDijDict[phfx+'D23'] = k*l |
---|
2152 | elif uniq == 'b': |
---|
2153 | dDijDict[phfx+'D13'] = h*l |
---|
2154 | elif uniq == 'c': |
---|
2155 | dDijDict[phfx+'D12'] = h*k |
---|
2156 | names.append() |
---|
2157 | else: |
---|
2158 | dDijDict = {phfx+'D11':h**2,phfx+'D22':k**2,phfx+'D33':l**2, |
---|
2159 | phfx+'D12':h*k,phfx+'D13':h*l,phfx+'D23':k*l} |
---|
2160 | for item in dDijDict: |
---|
2161 | dDijDict[item] *= refl[4]**2*tand(refl[5]/2.0) |
---|
2162 | return dDijDict |
---|
2163 | |
---|
2164 | def GetFprime(controlDict,Histograms): |
---|
2165 | FFtables = controlDict['FFtables'] |
---|
2166 | if not FFtables: |
---|
2167 | return |
---|
2168 | histoList = Histograms.keys() |
---|
2169 | histoList.sort() |
---|
2170 | for histogram in histoList: |
---|
2171 | if 'PWDR' in histogram[:4]: |
---|
2172 | Histogram = Histograms[histogram] |
---|
2173 | hId = Histogram['hId'] |
---|
2174 | hfx = ':%d:'%(hId) |
---|
2175 | keV = controlDict[hfx+'keV'] |
---|
2176 | for El in FFtables: |
---|
2177 | Orbs = G2el.GetXsectionCoeff(El.split('+')[0].split('-')[0]) |
---|
2178 | FP,FPP,Mu = G2el.FPcalc(Orbs, keV) |
---|
2179 | FFtables[El][hfx+'FP'] = FP |
---|
2180 | FFtables[El][hfx+'FPP'] = FPP |
---|
2181 | |
---|
2182 | def getPowderProfile(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): |
---|
2183 | |
---|
2184 | def GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict): |
---|
2185 | U = parmDict[hfx+'U'] |
---|
2186 | V = parmDict[hfx+'V'] |
---|
2187 | W = parmDict[hfx+'W'] |
---|
2188 | X = parmDict[hfx+'X'] |
---|
2189 | Y = parmDict[hfx+'Y'] |
---|
2190 | tanPos = tand(refl[5]/2.0) |
---|
2191 | sig = U*tanPos**2+V*tanPos+W #save peak sigma |
---|
2192 | sig = max(0.001,sig) |
---|
2193 | gam = X/cosd(refl[5]/2.0)+Y*tanPos+GetSampleGam(refl,wave,G,GB,phfx,calcControls,parmDict) #save peak gamma |
---|
2194 | gam = max(0.001,gam) |
---|
2195 | return sig,gam |
---|
2196 | |
---|
2197 | hId = Histogram['hId'] |
---|
2198 | hfx = ':%d:'%(hId) |
---|
2199 | bakType = calcControls[hfx+'bakType'] |
---|
2200 | yb = G2pwd.getBackground(hfx,parmDict,bakType,x) |
---|
2201 | yc = np.zeros_like(yb) |
---|
2202 | |
---|
2203 | if 'C' in calcControls[hfx+'histType']: |
---|
2204 | shl = max(parmDict[hfx+'SH/L'],0.002) |
---|
2205 | Ka2 = False |
---|
2206 | if hfx+'Lam1' in parmDict.keys(): |
---|
2207 | wave = parmDict[hfx+'Lam1'] |
---|
2208 | Ka2 = True |
---|
2209 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
---|
2210 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
---|
2211 | else: |
---|
2212 | wave = parmDict[hfx+'Lam'] |
---|
2213 | else: |
---|
2214 | print 'TOF Undefined at present' |
---|
2215 | raise ValueError |
---|
2216 | for phase in Histogram['Reflection Lists']: |
---|
2217 | refList = Histogram['Reflection Lists'][phase] |
---|
2218 | Phase = Phases[phase] |
---|
2219 | pId = Phase['pId'] |
---|
2220 | pfx = '%d::'%(pId) |
---|
2221 | phfx = '%d:%d:'%(pId,hId) |
---|
2222 | hfx = ':%d:'%(hId) |
---|
2223 | SGData = Phase['General']['SGData'] |
---|
2224 | A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] |
---|
2225 | G,g = G2lat.A2Gmat(A) #recip & real metric tensors |
---|
2226 | GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies |
---|
2227 | Vst = np.sqrt(nl.det(G)) #V* |
---|
2228 | if 'Pawley' not in Phase['General']['Type']: |
---|
2229 | refList = StructureFactor(refList,G,hfx,pfx,SGData,calcControls,parmDict) |
---|
2230 | for refl in refList: |
---|
2231 | if 'C' in calcControls[hfx+'histType']: |
---|
2232 | h,k,l = refl[:3] |
---|
2233 | refl[5] = GetReflPos(refl,wave,G,hfx,calcControls,parmDict) #corrected reflection position |
---|
2234 | Lorenz = 1./(2.*sind(refl[5]/2.)**2*cosd(refl[5]/2.)) #Lorentz correction |
---|
2235 | refl[5] += GetHStrainShift(refl,SGData,phfx,parmDict) #apply hydrostatic strain shift |
---|
2236 | refl[6:8] = GetReflSIgGam(refl,wave,G,GB,hfx,phfx,calcControls,parmDict) #peak sig & gam |
---|
2237 | GetIntensityCorr(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) #puts corrections in refl[13] |
---|
2238 | refl[13] *= Vst*Lorenz |
---|
2239 | if 'Pawley' in Phase['General']['Type']: |
---|
2240 | try: |
---|
2241 | refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) |
---|
2242 | except KeyError: |
---|
2243 | # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) |
---|
2244 | continue |
---|
2245 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
---|
2246 | iBeg = np.searchsorted(x,refl[5]-fmin) |
---|
2247 | iFin = np.searchsorted(x,refl[5]+fmax) |
---|
2248 | if not iBeg+iFin: #peak below low limit - skip peak |
---|
2249 | continue |
---|
2250 | elif not iBeg-iFin: #peak above high limit - done |
---|
2251 | break |
---|
2252 | yc[iBeg:iFin] += refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here |
---|
2253 | if Ka2: |
---|
2254 | pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) |
---|
2255 | Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) |
---|
2256 | iBeg = np.searchsorted(x,pos2-fmin) |
---|
2257 | iFin = np.searchsorted(x,pos2+fmax) |
---|
2258 | if not iBeg+iFin: #peak below low limit - skip peak |
---|
2259 | continue |
---|
2260 | elif not iBeg-iFin: #peak above high limit - done |
---|
2261 | return yc,yb |
---|
2262 | yc[iBeg:iFin] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg:iFin]) #and here |
---|
2263 | elif 'T' in calcControls[hfx+'histType']: |
---|
2264 | print 'TOF Undefined at present' |
---|
2265 | raise Exception #no TOF yet |
---|
2266 | return yc,yb |
---|
2267 | |
---|
2268 | def GetFobsSq(Histograms,Phases,parmDict,calcControls): |
---|
2269 | histoList = Histograms.keys() |
---|
2270 | histoList.sort() |
---|
2271 | for histogram in histoList: |
---|
2272 | if 'PWDR' in histogram[:4]: |
---|
2273 | Histogram = Histograms[histogram] |
---|
2274 | hId = Histogram['hId'] |
---|
2275 | hfx = ':%d:'%(hId) |
---|
2276 | Limits = calcControls[hfx+'Limits'] |
---|
2277 | shl = max(parmDict[hfx+'SH/L'],0.002) |
---|
2278 | Ka2 = False |
---|
2279 | kRatio = 0.0 |
---|
2280 | if hfx+'Lam1' in parmDict.keys(): |
---|
2281 | Ka2 = True |
---|
2282 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
---|
2283 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
---|
2284 | x,y,w,yc,yb,yd = Histogram['Data'] |
---|
2285 | ymb = np.array(y-yb) |
---|
2286 | ycmb = np.array(yc-yb) |
---|
2287 | ratio = np.where(ycmb!=0.,ymb/ycmb,0.0) |
---|
2288 | xB = np.searchsorted(x,Limits[0]) |
---|
2289 | xF = np.searchsorted(x,Limits[1]) |
---|
2290 | refLists = Histogram['Reflection Lists'] |
---|
2291 | for phase in refLists: |
---|
2292 | Phase = Phases[phase] |
---|
2293 | pId = Phase['pId'] |
---|
2294 | phfx = '%d:%d:'%(pId,hId) |
---|
2295 | refList = refLists[phase] |
---|
2296 | sumFo = 0.0 |
---|
2297 | sumdF = 0.0 |
---|
2298 | sumFosq = 0.0 |
---|
2299 | sumdFsq = 0.0 |
---|
2300 | for refl in refList: |
---|
2301 | if 'C' in calcControls[hfx+'histType']: |
---|
2302 | yp = np.zeros_like(yb) |
---|
2303 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
---|
2304 | iBeg = np.searchsorted(x[xB:xF],refl[5]-fmin) |
---|
2305 | iFin = np.searchsorted(x[xB:xF],refl[5]+fmax) |
---|
2306 | iFin2 = iFin |
---|
2307 | yp[iBeg:iFin] = refl[13]*refl[9]*G2pwd.getFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) #>90% of time spent here |
---|
2308 | if Ka2: |
---|
2309 | pos2 = refl[5]+lamRatio*tand(refl[5]/2.0) # + 360/pi * Dlam/lam * tan(th) |
---|
2310 | Wd,fmin,fmax = G2pwd.getWidths(pos2,refl[6],refl[7],shl) |
---|
2311 | iBeg2 = np.searchsorted(x,pos2-fmin) |
---|
2312 | iFin2 = np.searchsorted(x,pos2+fmax) |
---|
2313 | yp[iBeg2:iFin2] += refl[13]*refl[9]*kRatio*G2pwd.getFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) #and here |
---|
2314 | refl[8] = np.sum(np.where(ratio[iBeg:iFin2]>0.,yp[iBeg:iFin2]*ratio[iBeg:iFin2]/(refl[13]*(1.+kRatio)),0.0)) |
---|
2315 | elif 'T' in calcControls[hfx+'histType']: |
---|
2316 | print 'TOF Undefined at present' |
---|
2317 | raise Exception #no TOF yet |
---|
2318 | Fo = np.sqrt(np.abs(refl[8])) |
---|
2319 | Fc = np.sqrt(np.abs(refl[9])) |
---|
2320 | sumFo += Fo |
---|
2321 | sumFosq += refl[8]**2 |
---|
2322 | sumdF += np.abs(Fo-Fc) |
---|
2323 | sumdFsq += (refl[8]-refl[9])**2 |
---|
2324 | Histogram[phfx+'Rf'] = min(100.,(sumdF/sumFo)*100.) |
---|
2325 | Histogram[phfx+'Rf^2'] = min(100.,np.sqrt(sumdFsq/sumFosq)*100.) |
---|
2326 | Histogram[phfx+'Nref'] = len(refList) |
---|
2327 | |
---|
2328 | def getPowderProfileDerv(parmDict,x,varylist,Histogram,Phases,calcControls,pawleyLookup): |
---|
2329 | |
---|
2330 | def cellVaryDerv(pfx,SGData,dpdA): |
---|
2331 | if SGData['SGLaue'] in ['-1',]: |
---|
2332 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]], |
---|
2333 | [pfx+'A3',dpdA[3]],[pfx+'A4',dpdA[4]],[pfx+'A5',dpdA[5]]] |
---|
2334 | elif SGData['SGLaue'] in ['2/m',]: |
---|
2335 | if SGData['SGUniq'] == 'a': |
---|
2336 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A3',dpdA[3]]] |
---|
2337 | elif SGData['SGUniq'] == 'b': |
---|
2338 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A4',dpdA[4]]] |
---|
2339 | else: |
---|
2340 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]],[pfx+'A5',dpdA[5]]] |
---|
2341 | elif SGData['SGLaue'] in ['mmm',]: |
---|
2342 | return [[pfx+'A0',dpdA[0]],[pfx+'A1',dpdA[1]],[pfx+'A2',dpdA[2]]] |
---|
2343 | elif SGData['SGLaue'] in ['4/m','4/mmm']: |
---|
2344 | # return [[pfx+'A0',dpdA[0]+dpdA[1]],[pfx+'A2',dpdA[2]]] |
---|
2345 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
---|
2346 | elif SGData['SGLaue'] in ['6/m','6/mmm','3m1', '31m', '3']: |
---|
2347 | # return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[3]],[pfx+'A2',dpdA[2]]] |
---|
2348 | return [[pfx+'A0',dpdA[0]],[pfx+'A2',dpdA[2]]] |
---|
2349 | elif SGData['SGLaue'] in ['3R', '3mR']: |
---|
2350 | return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]],[pfx+'A3',dpdA[3]+dpdA[4]+dpdA[5]]] |
---|
2351 | elif SGData['SGLaue'] in ['m3m','m3']: |
---|
2352 | # return [[pfx+'A0',dpdA[0]+dpdA[1]+dpdA[2]]] |
---|
2353 | return [[pfx+'A0',dpdA[0]]] |
---|
2354 | # create a list of dependent variables and set up a dictionary to hold their derivatives |
---|
2355 | dependentVars = G2mv.GetDependentVars() |
---|
2356 | depDerivDict = {} |
---|
2357 | for j in dependentVars: |
---|
2358 | depDerivDict[j] = np.zeros(shape=(len(x))) |
---|
2359 | #print 'dependent vars',dependentVars |
---|
2360 | lenX = len(x) |
---|
2361 | hId = Histogram['hId'] |
---|
2362 | hfx = ':%d:'%(hId) |
---|
2363 | bakType = calcControls[hfx+'bakType'] |
---|
2364 | dMdv = np.zeros(shape=(len(varylist),len(x))) |
---|
2365 | dMdb,dMddb = G2pwd.getBackgroundDerv(hfx,parmDict,bakType,x) |
---|
2366 | if hfx+'Back:0' in varylist: # for now assume that Back:x vars to not appear in constraints |
---|
2367 | bBpos =varylist.index(hfx+'Back:0') |
---|
2368 | dMdv[bBpos:bBpos+len(dMdb)] = dMdb |
---|
2369 | names = [hfx+'DebyeA',hfx+'DebyeR',hfx+'DebyeU'] |
---|
2370 | for name in varylist: |
---|
2371 | if 'Debye' in name: |
---|
2372 | id = int(name.split(':')[-1]) |
---|
2373 | parm = name[:int(name.rindex(':'))] |
---|
2374 | ip = names.index(parm) |
---|
2375 | dMdv[varylist.index(name)] = dMddb[3*id+ip] |
---|
2376 | if 'C' in calcControls[hfx+'histType']: |
---|
2377 | dx = x[1]-x[0] |
---|
2378 | shl = max(parmDict[hfx+'SH/L'],0.002) |
---|
2379 | Ka2 = False |
---|
2380 | if hfx+'Lam1' in parmDict.keys(): |
---|
2381 | wave = parmDict[hfx+'Lam1'] |
---|
2382 | Ka2 = True |
---|
2383 | lamRatio = 360*(parmDict[hfx+'Lam2']-parmDict[hfx+'Lam1'])/(np.pi*parmDict[hfx+'Lam1']) |
---|
2384 | kRatio = parmDict[hfx+'I(L2)/I(L1)'] |
---|
2385 | else: |
---|
2386 | wave = parmDict[hfx+'Lam'] |
---|
2387 | else: |
---|
2388 | print 'TOF Undefined at present' |
---|
2389 | raise ValueError |
---|
2390 | for phase in Histogram['Reflection Lists']: |
---|
2391 | refList = Histogram['Reflection Lists'][phase] |
---|
2392 | Phase = Phases[phase] |
---|
2393 | SGData = Phase['General']['SGData'] |
---|
2394 | pId = Phase['pId'] |
---|
2395 | pfx = '%d::'%(pId) |
---|
2396 | phfx = '%d:%d:'%(pId,hId) |
---|
2397 | A = [parmDict[pfx+'A%d'%(i)] for i in range(6)] |
---|
2398 | G,g = G2lat.A2Gmat(A) #recip & real metric tensors |
---|
2399 | GA,GB = G2lat.Gmat2AB(G) #Orthogonalization matricies |
---|
2400 | if 'Pawley' not in Phase['General']['Type']: |
---|
2401 | dFdvDict = StructureFactorDerv(refList,G,hfx,pfx,SGData,calcControls,parmDict) |
---|
2402 | for iref,refl in enumerate(refList): |
---|
2403 | if 'C' in calcControls[hfx+'histType']: #CW powder |
---|
2404 | h,k,l = refl[:3] |
---|
2405 | dIdsh,dIdsp,dIdpola,dIdPO,dFdODF,dFdSA = GetIntensityDerv(refl,G,g,pfx,phfx,hfx,SGData,calcControls,parmDict) |
---|
2406 | if 'Pawley' in Phase['General']['Type']: |
---|
2407 | try: |
---|
2408 | refl[9] = abs(parmDict[pfx+'PWLref:%d'%(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])]) |
---|
2409 | except KeyError: |
---|
2410 | # print ' ***Error %d,%d,%d missing from Pawley reflection list ***'%(h,k,l) |
---|
2411 | continue |
---|
2412 | Wd,fmin,fmax = G2pwd.getWidths(refl[5],refl[6],refl[7],shl) |
---|
2413 | iBeg = np.searchsorted(x,refl[5]-fmin) |
---|
2414 | iFin = np.searchsorted(x,refl[5]+fmax) |
---|
2415 | if not iBeg+iFin: #peak below low limit - skip peak |
---|
2416 | continue |
---|
2417 | elif not iBeg-iFin: #peak above high limit - done |
---|
2418 | break |
---|
2419 | pos = refl[5] |
---|
2420 | tanth = tand(pos/2.0) |
---|
2421 | costh = cosd(pos/2.0) |
---|
2422 | dMdpk = np.zeros(shape=(6,len(x))) |
---|
2423 | dMdipk = G2pwd.getdFCJVoigt3(refl[5],refl[6],refl[7],shl,x[iBeg:iFin]) |
---|
2424 | for i in range(1,5): |
---|
2425 | dMdpk[i][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[i] |
---|
2426 | dMdpk[0][iBeg:iFin] += 100.*dx*refl[13]*refl[9]*dMdipk[0] |
---|
2427 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} |
---|
2428 | if Ka2: |
---|
2429 | pos2 = refl[5]+lamRatio*tanth # + 360/pi * Dlam/lam * tan(th) |
---|
2430 | kdelt = int((pos2-refl[5])/dx) |
---|
2431 | iBeg2 = min(lenX,iBeg+kdelt) |
---|
2432 | iFin2 = min(lenX,iFin+kdelt) |
---|
2433 | if iBeg2-iFin2: |
---|
2434 | dMdipk2 = G2pwd.getdFCJVoigt3(pos2,refl[6],refl[7],shl,x[iBeg2:iFin2]) |
---|
2435 | for i in range(1,5): |
---|
2436 | dMdpk[i][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[i] |
---|
2437 | dMdpk[0][iBeg2:iFin2] += 100.*dx*refl[13]*refl[9]*kRatio*dMdipk2[0] |
---|
2438 | dMdpk[5][iBeg2:iFin2] += 100.*dx*refl[13]*dMdipk2[0] |
---|
2439 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*refl[9]} |
---|
2440 | if 'Pawley' in Phase['General']['Type']: |
---|
2441 | try: |
---|
2442 | idx = varylist.index(pfx+'PWLref:'+str(pawleyLookup[pfx+'%d,%d,%d'%(h,k,l)])) |
---|
2443 | dMdv[idx] = dervDict['int']/refl[9] |
---|
2444 | # Assuming Pawley variables not in constraints |
---|
2445 | except ValueError: |
---|
2446 | pass |
---|
2447 | dpdA,dpdw,dpdZ,dpdSh,dpdTr,dpdX,dpdY = GetReflPosDerv(refl,wave,A,hfx,calcControls,parmDict) |
---|
2448 | names = {hfx+'Scale':[dIdsh,'int'],hfx+'Polariz.':[dIdpola,'int'],phfx+'Scale':[dIdsp,'int'], |
---|
2449 | hfx+'U':[tanth**2,'sig'],hfx+'V':[tanth,'sig'],hfx+'W':[1.0,'sig'], |
---|
2450 | hfx+'X':[1.0/costh,'gam'],hfx+'Y':[tanth,'gam'],hfx+'SH/L':[1.0,'shl'], |
---|
2451 | hfx+'I(L2)/I(L1)':[1.0,'L1/L2'],hfx+'Zero':[dpdZ,'pos'],hfx+'Lam':[dpdw,'pos'], |
---|
2452 | hfx+'Shift':[dpdSh,'pos'],hfx+'Transparency':[dpdTr,'pos'],hfx+'DisplaceX':[dpdX,'pos'], |
---|
2453 | hfx+'DisplaceY':[dpdY,'pos'],} |
---|
2454 | for name in names: |
---|
2455 | item = names[name] |
---|
2456 | if name in varylist: |
---|
2457 | dMdv[varylist.index(name)] += item[0]*dervDict[item[1]] |
---|
2458 | elif name in dependentVars: |
---|
2459 | depDerivDict[name] += item[0]*dervDict[item[1]] |
---|
2460 | |
---|
2461 | for iPO in dIdPO: |
---|
2462 | if iPO in varylist: |
---|
2463 | dMdv[varylist.index(iPO)] += dIdPO[iPO]*dervDict['int'] |
---|
2464 | elif iPO in dependentVars: |
---|
2465 | depDerivDict[iPO] = dIdPO[iPO]*dervDict['int'] |
---|
2466 | |
---|
2467 | for i,name in enumerate(['omega','chi','phi']): |
---|
2468 | aname = pfx+'SH '+name |
---|
2469 | if aname in varylist: |
---|
2470 | dMdv[varylist.index(aname)] += dFdSA[i]*dervDict['int'] |
---|
2471 | elif aname in dependentVars: |
---|
2472 | depDerivDict[aname] += dFdSA[i]*dervDict['int'] |
---|
2473 | for iSH in dFdODF: |
---|
2474 | if iSH in varylist: |
---|
2475 | dMdv[varylist.index(iSH)] += dFdODF[iSH]*dervDict['int'] |
---|
2476 | elif iSH in dependentVars: |
---|
2477 | depDerivDict[iSH] += dFdODF[iSH]*dervDict['int'] |
---|
2478 | cellDervNames = cellVaryDerv(pfx,SGData,dpdA) |
---|
2479 | for name,dpdA in cellDervNames: |
---|
2480 | if name in varylist: |
---|
2481 | dMdv[varylist.index(name)] += dpdA*dervDict['pos'] |
---|
2482 | elif name in dependentVars: |
---|
2483 | depDerivDict[name] += dpdA*dervDict['pos'] |
---|
2484 | dDijDict = GetHStrainShiftDerv(refl,SGData,phfx) |
---|
2485 | for name in dDijDict: |
---|
2486 | if name in varylist: |
---|
2487 | dMdv[varylist.index(name)] += dDijDict[name]*dervDict['pos'] |
---|
2488 | elif name in dependentVars: |
---|
2489 | depDerivDict[name] += dDijDict[name]*dervDict['pos'] |
---|
2490 | gamDict = GetSampleGamDerv(refl,wave,G,GB,phfx,calcControls,parmDict) |
---|
2491 | for name in gamDict: |
---|
2492 | if name in varylist: |
---|
2493 | dMdv[varylist.index(name)] += gamDict[name]*dervDict['gam'] |
---|
2494 | elif name in dependentVars: |
---|
2495 | depDerivDict[name] += gamDict[name]*dervDict['gam'] |
---|
2496 | |
---|
2497 | elif 'T' in calcControls[hfx+'histType']: |
---|
2498 | print 'TOF Undefined at present' |
---|
2499 | raise Exception #no TOF yet |
---|
2500 | #do atom derivatives - for F,X & U so far |
---|
2501 | corr = dervDict['int']/refl[9] |
---|
2502 | for name in varylist+dependentVars: |
---|
2503 | try: |
---|
2504 | aname = name.split(pfx)[1][:2] |
---|
2505 | if aname not in ['Af','dA','AU']: continue # skip anything not an atom param |
---|
2506 | except IndexError: |
---|
2507 | continue |
---|
2508 | if name in varylist: |
---|
2509 | dMdv[varylist.index(name)] += dFdvDict[name][iref]*corr |
---|
2510 | elif name in dependentVars: |
---|
2511 | depDerivDict[name] += dFdvDict[name][iref]*corr |
---|
2512 | # now process derivatives in constraints |
---|
2513 | G2mv.Dict2Deriv(varylist,depDerivDict,dMdv) |
---|
2514 | return dMdv |
---|
2515 | |
---|
2516 | def dervRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
---|
2517 | parmdict.update(zip(varylist,values)) |
---|
2518 | G2mv.Dict2Map(parmdict,varylist) |
---|
2519 | Histograms,Phases = HistoPhases |
---|
2520 | nvar = len(varylist) |
---|
2521 | dMdv = np.empty(0) |
---|
2522 | histoList = Histograms.keys() |
---|
2523 | histoList.sort() |
---|
2524 | for histogram in histoList: |
---|
2525 | if 'PWDR' in histogram[:4]: |
---|
2526 | Histogram = Histograms[histogram] |
---|
2527 | hId = Histogram['hId'] |
---|
2528 | hfx = ':%d:'%(hId) |
---|
2529 | Limits = calcControls[hfx+'Limits'] |
---|
2530 | x,y,w,yc,yb,yd = Histogram['Data'] |
---|
2531 | xB = np.searchsorted(x,Limits[0]) |
---|
2532 | xF = np.searchsorted(x,Limits[1]) |
---|
2533 | dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF], |
---|
2534 | varylist,Histogram,Phases,calcControls,pawleyLookup) |
---|
2535 | if len(dMdv): |
---|
2536 | dMdv = np.concatenate((dMdv.T,dMdvh.T)).T |
---|
2537 | else: |
---|
2538 | dMdv = dMdvh |
---|
2539 | return dMdv |
---|
2540 | |
---|
2541 | def HessRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
---|
2542 | parmdict.update(zip(varylist,values)) |
---|
2543 | G2mv.Dict2Map(parmdict,varylist) |
---|
2544 | Histograms,Phases = HistoPhases |
---|
2545 | nvar = len(varylist) |
---|
2546 | Hess = np.empty(0) |
---|
2547 | histoList = Histograms.keys() |
---|
2548 | histoList.sort() |
---|
2549 | for histogram in histoList: |
---|
2550 | if 'PWDR' in histogram[:4]: |
---|
2551 | Histogram = Histograms[histogram] |
---|
2552 | hId = Histogram['hId'] |
---|
2553 | hfx = ':%d:'%(hId) |
---|
2554 | Limits = calcControls[hfx+'Limits'] |
---|
2555 | x,y,w,yc,yb,yd = Histogram['Data'] |
---|
2556 | dy = y-yc |
---|
2557 | xB = np.searchsorted(x,Limits[0]) |
---|
2558 | xF = np.searchsorted(x,Limits[1]) |
---|
2559 | dMdvh = np.sqrt(w[xB:xF])*getPowderProfileDerv(parmdict,x[xB:xF], |
---|
2560 | varylist,Histogram,Phases,calcControls,pawleyLookup) |
---|
2561 | if dlg: |
---|
2562 | dlg.Update(Histogram['wRp'],newmsg='Hessian for histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0] |
---|
2563 | if len(Hess): |
---|
2564 | Vec += np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) |
---|
2565 | Hess += np.inner(dMdvh,dMdvh) |
---|
2566 | else: |
---|
2567 | Vec = np.sum(dMdvh*np.sqrt(w[xB:xF])*dy[xB:xF],axis=1) |
---|
2568 | Hess = np.inner(dMdvh,dMdvh) |
---|
2569 | return Vec,Hess |
---|
2570 | |
---|
2571 | def errRefine(values,HistoPhases,parmdict,varylist,calcControls,pawleyLookup,dlg): |
---|
2572 | parmdict.update(zip(varylist,values)) |
---|
2573 | Values2Dict(parmdict, varylist, values) |
---|
2574 | G2mv.Dict2Map(parmdict,varylist) |
---|
2575 | Histograms,Phases = HistoPhases |
---|
2576 | M = np.empty(0) |
---|
2577 | sumwYo = 0 |
---|
2578 | Nobs = 0 |
---|
2579 | histoList = Histograms.keys() |
---|
2580 | histoList.sort() |
---|
2581 | for histogram in histoList: |
---|
2582 | if 'PWDR' in histogram[:4]: |
---|
2583 | Histogram = Histograms[histogram] |
---|
2584 | hId = Histogram['hId'] |
---|
2585 | hfx = ':%d:'%(hId) |
---|
2586 | Limits = calcControls[hfx+'Limits'] |
---|
2587 | x,y,w,yc,yb,yd = Histogram['Data'] |
---|
2588 | yc *= 0.0 #zero full calcd profiles |
---|
2589 | yb *= 0.0 |
---|
2590 | yd *= 0.0 |
---|
2591 | xB = np.searchsorted(x,Limits[0]) |
---|
2592 | xF = np.searchsorted(x,Limits[1]) |
---|
2593 | Histogram['Nobs'] = xF-xB |
---|
2594 | Nobs += Histogram['Nobs'] |
---|
2595 | Histogram['sumwYo'] = np.sum(w[xB:xF]*y[xB:xF]**2) |
---|
2596 | sumwYo += Histogram['sumwYo'] |
---|
2597 | yc[xB:xF],yb[xB:xF] = getPowderProfile(parmdict,x[xB:xF], |
---|
2598 | varylist,Histogram,Phases,calcControls,pawleyLookup) |
---|
2599 | yc[xB:xF] += yb[xB:xF] |
---|
2600 | yd[xB:xF] = y[xB:xF]-yc[xB:xF] |
---|
2601 | Histogram['sumwYd'] = np.sum(np.sqrt(w[xB:xF])*(yd[xB:xF])) |
---|
2602 | wdy = -np.sqrt(w[xB:xF])*(yd[xB:xF]) |
---|
2603 | Histogram['wRp'] = min(100.,np.sqrt(np.sum(wdy**2)/Histogram['sumwYo'])*100.) |
---|
2604 | if dlg: |
---|
2605 | dlg.Update(Histogram['wRp'],newmsg='For histogram %d Rwp=%8.3f%s'%(hId,Histogram['wRp'],'%'))[0] |
---|
2606 | M = np.concatenate((M,wdy)) |
---|
2607 | Histograms['sumwYo'] = sumwYo |
---|
2608 | Histograms['Nobs'] = Nobs |
---|
2609 | Rwp = min(100.,np.sqrt(np.sum(M**2)/sumwYo)*100.) |
---|
2610 | if dlg: |
---|
2611 | GoOn = dlg.Update(Rwp,newmsg='%s%8.3f%s'%('Powder profile Rwp =',Rwp,'%'))[0] |
---|
2612 | if not GoOn: |
---|
2613 | parmDict['saved values'] = values |
---|
2614 | raise Exception #Abort!! |
---|
2615 | return M |
---|
2616 | |
---|
2617 | |
---|
2618 | def Refine(GPXfile,dlg): |
---|
2619 | import cPickle |
---|
2620 | import pytexture as ptx |
---|
2621 | ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics |
---|
2622 | |
---|
2623 | ShowBanner() |
---|
2624 | varyList = [] |
---|
2625 | parmDict = {} |
---|
2626 | calcControls = {} |
---|
2627 | G2mv.InitVars() |
---|
2628 | Controls = GetControls(GPXfile) |
---|
2629 | ShowControls(Controls) |
---|
2630 | constrDict,constrFlag,fixedList = GetConstraints(GPXfile) |
---|
2631 | Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) |
---|
2632 | if not Phases: |
---|
2633 | print ' *** ERROR - you have no histograms to refine! ***' |
---|
2634 | print ' *** Refine aborted ***' |
---|
2635 | raise Exception |
---|
2636 | if not Histograms: |
---|
2637 | print ' *** ERROR - you have no data to refine with! ***' |
---|
2638 | print ' *** Refine aborted ***' |
---|
2639 | raise Exception |
---|
2640 | Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases) |
---|
2641 | calcControls['Natoms'] = Natoms |
---|
2642 | calcControls['FFtables'] = FFtables |
---|
2643 | calcControls['BLtables'] = BLtables |
---|
2644 | hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histograms) |
---|
2645 | calcControls.update(controlDict) |
---|
2646 | histVary,histDict,controlDict = GetHistogramData(Histograms) |
---|
2647 | calcControls.update(controlDict) |
---|
2648 | varyList = phaseVary+hapVary+histVary |
---|
2649 | parmDict.update(phaseDict) |
---|
2650 | parmDict.update(hapDict) |
---|
2651 | parmDict.update(histDict) |
---|
2652 | GetFprime(calcControls,Histograms) |
---|
2653 | # do constraint processing |
---|
2654 | try: |
---|
2655 | groups,parmlist = G2mv.GroupConstraints(constrDict) |
---|
2656 | G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList) |
---|
2657 | except: |
---|
2658 | print ' *** ERROR - your constraints are internally inconsistent ***' |
---|
2659 | raise Exception(' *** Refine aborted ***') |
---|
2660 | # check to see which generated parameters are fully varied |
---|
2661 | msg = G2mv.SetVaryFlags(varyList) |
---|
2662 | if msg: |
---|
2663 | print ' *** ERROR - you have not set the refine flags for constraints consistently! ***' |
---|
2664 | print msg |
---|
2665 | raise Exception(' *** Refine aborted ***') |
---|
2666 | G2mv.Map2Dict(parmDict,varyList) |
---|
2667 | # print G2mv.VarRemapShow(varyList) |
---|
2668 | |
---|
2669 | while True: |
---|
2670 | begin = time.time() |
---|
2671 | values = np.array(Dict2Values(parmDict, varyList)) |
---|
2672 | Ftol = Controls['min dM/M'] |
---|
2673 | Factor = Controls['shift factor'] |
---|
2674 | maxCyc = Controls['max cyc'] |
---|
2675 | if 'Jacobian' in Controls['deriv type']: |
---|
2676 | result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True, |
---|
2677 | ftol=Ftol,col_deriv=True,factor=Factor, |
---|
2678 | args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
---|
2679 | ncyc = int(result[2]['nfev']/2) |
---|
2680 | elif 'Hessian' in Controls['deriv type']: |
---|
2681 | result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc, |
---|
2682 | args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
---|
2683 | ncyc = result[2]['num cyc']+1 |
---|
2684 | else: #'numeric' |
---|
2685 | result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor, |
---|
2686 | args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
---|
2687 | ncyc = int(result[2]['nfev']/len(varyList)) |
---|
2688 | # table = dict(zip(varyList,zip(values,result[0],(result[0]-values)))) |
---|
2689 | # for item in table: print item,table[item] #useful debug - are things shifting? |
---|
2690 | runtime = time.time()-begin |
---|
2691 | chisq = np.sum(result[2]['fvec']**2) |
---|
2692 | Values2Dict(parmDict, varyList, result[0]) |
---|
2693 | G2mv.Dict2Map(parmDict,varyList) |
---|
2694 | |
---|
2695 | Rwp = np.sqrt(chisq/Histograms['sumwYo'])*100. #to % |
---|
2696 | GOF = chisq/(Histograms['Nobs']-len(varyList)) |
---|
2697 | print '\n Refinement results:' |
---|
2698 | print 135*'-' |
---|
2699 | print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histograms['Nobs'],' Number of parameters: ',len(varyList) |
---|
2700 | print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc) |
---|
2701 | print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF) |
---|
2702 | try: |
---|
2703 | covMatrix = result[1]*GOF |
---|
2704 | sig = np.sqrt(np.diag(covMatrix)) |
---|
2705 | if np.any(np.isnan(sig)): |
---|
2706 | print '*** Least squares aborted - some invalid esds possible ***' |
---|
2707 | # table = dict(zip(varyList,zip(values,result[0],(result[0]-values)/sig))) |
---|
2708 | # for item in table: print item,table[item] #useful debug - are things shifting? |
---|
2709 | break #refinement succeeded - finish up! |
---|
2710 | except TypeError: #result[1] is None on singular matrix |
---|
2711 | print '**** Refinement failed - singular matrix ****' |
---|
2712 | if 'Hessian' in Controls['deriv type']: |
---|
2713 | for i in result[2]['psing'].reverse(): |
---|
2714 | print 'Removing parameter: ',varyList[i] |
---|
2715 | del(varyList[i]) |
---|
2716 | else: |
---|
2717 | Ipvt = result[2]['ipvt'] |
---|
2718 | for i,ipvt in enumerate(Ipvt): |
---|
2719 | if not np.sum(result[2]['fjac'],axis=1)[i]: |
---|
2720 | print 'Removing parameter: ',varyList[ipvt-1] |
---|
2721 | del(varyList[ipvt-1]) |
---|
2722 | break |
---|
2723 | |
---|
2724 | # print 'dependentParmList: ',G2mv.dependentParmList |
---|
2725 | # print 'arrayList: ',G2mv.arrayList |
---|
2726 | # print 'invarrayList: ',G2mv.invarrayList |
---|
2727 | # print 'indParmList: ',G2mv.indParmList |
---|
2728 | # print 'fixedDict: ',G2mv.fixedDict |
---|
2729 | # print 'test1' |
---|
2730 | GetFobsSq(Histograms,Phases,parmDict,calcControls) |
---|
2731 | # print 'test2' |
---|
2732 | sigDict = dict(zip(varyList,sig)) |
---|
2733 | newCellDict = GetNewCellParms(parmDict,varyList) |
---|
2734 | newAtomDict = ApplyXYZshifts(parmDict,varyList) |
---|
2735 | covData = {'variables':result[0],'varyList':varyList,'sig':sig, |
---|
2736 | 'covMatrix':covMatrix,'title':GPXfile,'newAtomDict':newAtomDict,'newCellDict':newCellDict} |
---|
2737 | # add the uncertainties into the esd dictionary (sigDict) |
---|
2738 | sigDict.update(G2mv.ComputeDepESD(covMatrix,varyList,parmDict)) |
---|
2739 | SetPhaseData(parmDict,sigDict,Phases,covData) |
---|
2740 | SetHistogramPhaseData(parmDict,sigDict,Phases,Histograms) |
---|
2741 | SetHistogramData(parmDict,sigDict,Histograms) |
---|
2742 | G2mv.PrintIndependentVars(parmDict,varyList,sigDict) |
---|
2743 | SetUsedHistogramsAndPhases(GPXfile,Histograms,Phases,covData) |
---|
2744 | #for testing purposes!!! |
---|
2745 | # file = open('structTestdata.dat','wb') |
---|
2746 | # cPickle.dump(parmDict,file,1) |
---|
2747 | # cPickle.dump(varyList,file,1) |
---|
2748 | # for histogram in Histograms: |
---|
2749 | # if 'PWDR' in histogram[:4]: |
---|
2750 | # Histogram = Histograms[histogram] |
---|
2751 | # cPickle.dump(Histogram,file,1) |
---|
2752 | # cPickle.dump(Phases,file,1) |
---|
2753 | # cPickle.dump(calcControls,file,1) |
---|
2754 | # cPickle.dump(pawleyLookup,file,1) |
---|
2755 | # file.close() |
---|
2756 | if dlg: |
---|
2757 | return Rwp |
---|
2758 | |
---|
2759 | def SeqRefine(GPXfile,dlg): |
---|
2760 | import cPickle |
---|
2761 | import pytexture as ptx |
---|
2762 | ptx.pyqlmninit() #initialize fortran arrays for spherical harmonics |
---|
2763 | |
---|
2764 | ShowBanner() |
---|
2765 | print ' Sequential Refinement' |
---|
2766 | G2mv.InitVars() |
---|
2767 | Controls = GetControls(GPXfile) |
---|
2768 | ShowControls(Controls) |
---|
2769 | constrDict,constrFlag,fixedList = GetConstraints(GPXfile) |
---|
2770 | Histograms,Phases = GetUsedHistogramsAndPhases(GPXfile) |
---|
2771 | if not Phases: |
---|
2772 | print ' *** ERROR - you have no histograms to refine! ***' |
---|
2773 | print ' *** Refine aborted ***' |
---|
2774 | raise Exception |
---|
2775 | if not Histograms: |
---|
2776 | print ' *** ERROR - you have no data to refine with! ***' |
---|
2777 | print ' *** Refine aborted ***' |
---|
2778 | raise Exception |
---|
2779 | Natoms,phaseVary,phaseDict,pawleyLookup,FFtables,BLtables = GetPhaseData(Phases,False) |
---|
2780 | if 'Seq Data' in Controls: |
---|
2781 | histNames = Controls['Seq Data'] |
---|
2782 | else: |
---|
2783 | histNames = GetHistogramNames(GPXfile,['PWDR',]) |
---|
2784 | if 'Reverse Seq' in Controls: |
---|
2785 | if Controls['Reverse Seq']: |
---|
2786 | histNames.reverse() |
---|
2787 | SeqResult = {'histNames':histNames} |
---|
2788 | makeBack = True |
---|
2789 | for ihst,histogram in enumerate(histNames): |
---|
2790 | ifPrint = False |
---|
2791 | if dlg: |
---|
2792 | dlg.SetTitle('Residual for histogram '+str(ihst)) |
---|
2793 | calcControls = {} |
---|
2794 | calcControls['Natoms'] = Natoms |
---|
2795 | calcControls['FFtables'] = FFtables |
---|
2796 | calcControls['BLtables'] = BLtables |
---|
2797 | varyList = [] |
---|
2798 | parmDict = {} |
---|
2799 | Histo = {histogram:Histograms[histogram],} |
---|
2800 | hapVary,hapDict,controlDict = GetHistogramPhaseData(Phases,Histo,False) |
---|
2801 | calcControls.update(controlDict) |
---|
2802 | histVary,histDict,controlDict = GetHistogramData(Histo,False) |
---|
2803 | calcControls.update(controlDict) |
---|
2804 | varyList = phaseVary+hapVary+histVary |
---|
2805 | if not ihst: |
---|
2806 | saveVaryList = varyList[:] |
---|
2807 | for i,item in enumerate(saveVaryList): |
---|
2808 | items = item.split(':') |
---|
2809 | if items[1]: |
---|
2810 | items[1] = '' |
---|
2811 | item = ':'.join(items) |
---|
2812 | saveVaryList[i] = item |
---|
2813 | SeqResult['varyList'] = saveVaryList |
---|
2814 | else: |
---|
2815 | newVaryList = varyList[:] |
---|
2816 | for i,item in enumerate(newVaryList): |
---|
2817 | items = item.split(':') |
---|
2818 | if items[1]: |
---|
2819 | items[1] = '' |
---|
2820 | item = ':'.join(items) |
---|
2821 | newVaryList[i] = item |
---|
2822 | if newVaryList != SeqResult['varyList']: |
---|
2823 | print newVaryList |
---|
2824 | print SeqResult['varyList'] |
---|
2825 | print '**** ERROR - variable list for this histogram does not match previous' |
---|
2826 | raise Exception |
---|
2827 | parmDict.update(phaseDict) |
---|
2828 | parmDict.update(hapDict) |
---|
2829 | parmDict.update(histDict) |
---|
2830 | GetFprime(calcControls,Histo) |
---|
2831 | constrDict,constrFlag,fixedList = G2mv.InputParse([]) #constraints go here? |
---|
2832 | groups,parmlist = G2mv.GroupConstraints(constrDict) |
---|
2833 | G2mv.GenerateConstraints(groups,parmlist,varyList,constrDict,constrFlag,fixedList) |
---|
2834 | G2mv.Map2Dict(parmDict,varyList) |
---|
2835 | |
---|
2836 | while True: |
---|
2837 | begin = time.time() |
---|
2838 | values = np.array(Dict2Values(parmDict, varyList)) |
---|
2839 | Ftol = Controls['min dM/M'] |
---|
2840 | Factor = Controls['shift factor'] |
---|
2841 | maxCyc = Controls['max cyc'] |
---|
2842 | |
---|
2843 | if 'Jacobian' in Controls['deriv type']: |
---|
2844 | result = so.leastsq(errRefine,values,Dfun=dervRefine,full_output=True, |
---|
2845 | ftol=Ftol,col_deriv=True,factor=Factor, |
---|
2846 | args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
---|
2847 | ncyc = int(result[2]['nfev']/2) |
---|
2848 | elif 'Hessian' in Controls['deriv type']: |
---|
2849 | result = G2mth.HessianLSQ(errRefine,values,Hess=HessRefine,ftol=Ftol,maxcyc=maxCyc, |
---|
2850 | args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
---|
2851 | ncyc = result[2]['num cyc']+1 |
---|
2852 | else: #'numeric' |
---|
2853 | result = so.leastsq(errRefine,values,full_output=True,ftol=Ftol,epsfcn=1.e-8,factor=Factor, |
---|
2854 | args=([Histograms,Phases],parmDict,varyList,calcControls,pawleyLookup,dlg)) |
---|
2855 | ncyc = int(result[2]['nfev']/len(varyList)) |
---|
2856 | |
---|
2857 | |
---|
2858 | |
---|
2859 | runtime = time.time()-begin |
---|
2860 | chisq = np.sum(result[2]['fvec']**2) |
---|
2861 | Values2Dict(parmDict, varyList, result[0]) |
---|
2862 | G2mv.Dict2Map(parmDict,varyList) |
---|
2863 | |
---|
2864 | Rwp = np.sqrt(chisq/Histo['sumwYo'])*100. #to % |
---|
2865 | GOF = chisq/(Histo['Nobs']-len(varyList)) |
---|
2866 | print '\n Refinement results for histogram: v'+histogram |
---|
2867 | print 135*'-' |
---|
2868 | print ' Number of function calls:',result[2]['nfev'],' Number of observations: ',Histo['Nobs'],' Number of parameters: ',len(varyList) |
---|
2869 | print ' Refinement time = %8.3fs, %8.3fs/cycle, for %d cycles'%(runtime,runtime/ncyc,ncyc) |
---|
2870 | print ' wRp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rwp,chisq,GOF) |
---|
2871 | try: |
---|
2872 | covMatrix = result[1]*GOF |
---|
2873 | sig = np.sqrt(np.diag(covMatrix)) |
---|
2874 | if np.any(np.isnan(sig)): |
---|
2875 | print '*** Least squares aborted - some invalid esds possible ***' |
---|
2876 | ifPrint = True |
---|
2877 | break #refinement succeeded - finish up! |
---|
2878 | except TypeError: #result[1] is None on singular matrix |
---|
2879 | print '**** Refinement failed - singular matrix ****' |
---|
2880 | if 'Hessian' in Controls['deriv type']: |
---|
2881 | for i in result[2]['psing'].reverse(): |
---|
2882 | print 'Removing parameter: ',varyList[i] |
---|
2883 | del(varyList[i]) |
---|
2884 | else: |
---|
2885 | Ipvt = result[2]['ipvt'] |
---|
2886 | for i,ipvt in enumerate(Ipvt): |
---|
2887 | if not np.sum(result[2]['fjac'],axis=1)[i]: |
---|
2888 | print 'Removing parameter: ',varyList[ipvt-1] |
---|
2889 | del(varyList[ipvt-1]) |
---|
2890 | break |
---|
2891 | |
---|
2892 | GetFobsSq(Histo,Phases,parmDict,calcControls) |
---|
2893 | sigDict = dict(zip(varyList,sig)) |
---|
2894 | newCellDict = GetNewCellParms(parmDict,varyList) |
---|
2895 | newAtomDict = ApplyXYZshifts(parmDict,varyList) |
---|
2896 | covData = {'variables':result[0],'varyList':varyList,'sig':sig, |
---|
2897 | 'covMatrix':covMatrix,'title':histogram,'newAtomDict':newAtomDict,'newCellDict':newCellDict} |
---|
2898 | SetHistogramPhaseData(parmDict,sigDict,Phases,Histo,ifPrint) |
---|
2899 | SetHistogramData(parmDict,sigDict,Histo,ifPrint) |
---|
2900 | SeqResult[histogram] = covData |
---|
2901 | SetUsedHistogramsAndPhases(GPXfile,Histo,Phases,covData,makeBack) |
---|
2902 | makeBack = False |
---|
2903 | SetSeqResult(GPXfile,Histograms,SeqResult) |
---|
2904 | |
---|
2905 | |
---|
2906 | def main(): |
---|
2907 | arg = sys.argv |
---|
2908 | if len(arg) > 1: |
---|
2909 | GPXfile = arg[1] |
---|
2910 | if not ospath.exists(GPXfile): |
---|
2911 | print 'ERROR - ',GPXfile," doesn't exist!" |
---|
2912 | exit() |
---|
2913 | GPXpath = ospath.dirname(arg[1]) |
---|
2914 | Refine(GPXfile,None) |
---|
2915 | else: |
---|
2916 | print 'ERROR - missing filename' |
---|
2917 | exit() |
---|
2918 | print "Done" |
---|
2919 | |
---|
2920 | if __name__ == '__main__': |
---|
2921 | main() |
---|