1 | #/usr/bin/env python |
---|
2 | # -*- coding: utf-8 -*- |
---|
3 | ''' |
---|
4 | *GSASIIpwd: Powder calculations module* |
---|
5 | ============================================== |
---|
6 | |
---|
7 | This version hacked to provide Laue Fringe fitting. |
---|
8 | |
---|
9 | ''' |
---|
10 | ########### SVN repository information ################### |
---|
11 | # $Date: 2023-01-17 13:31:50 +0000 (Tue, 17 Jan 2023) $ |
---|
12 | # $Author: toby $ |
---|
13 | # $Revision: 5478 $ |
---|
14 | # $URL: trunk/GSASIIpwd.py $ |
---|
15 | # $Id: GSASIIpwd.py 5478 2023-01-17 13:31:50Z toby $ |
---|
16 | ########### SVN repository information ################### |
---|
17 | from __future__ import division, print_function |
---|
18 | import sys |
---|
19 | import math |
---|
20 | import time |
---|
21 | import os |
---|
22 | import os.path |
---|
23 | import subprocess as subp |
---|
24 | import datetime as dt |
---|
25 | import copy |
---|
26 | |
---|
27 | import numpy as np |
---|
28 | import numpy.linalg as nl |
---|
29 | import numpy.ma as ma |
---|
30 | import random as rand |
---|
31 | import numpy.fft as fft |
---|
32 | import scipy.interpolate as si |
---|
33 | import scipy.stats as st |
---|
34 | import scipy.optimize as so |
---|
35 | import scipy.special as sp |
---|
36 | import scipy.signal as signal |
---|
37 | |
---|
38 | import GSASIIpath |
---|
39 | filversion = "$Revision: 5478 $" |
---|
40 | GSASIIpath.SetVersionNumber("$Revision: 5478 $") |
---|
41 | import GSASIIlattice as G2lat |
---|
42 | import GSASIIspc as G2spc |
---|
43 | import GSASIIElem as G2elem |
---|
44 | import GSASIImath as G2mth |
---|
45 | try: |
---|
46 | import pypowder as pyd |
---|
47 | except ImportError: |
---|
48 | print ('pypowder is not available - profile calcs. not allowed') |
---|
49 | try: |
---|
50 | import pydiffax as pyx |
---|
51 | except ImportError: |
---|
52 | print ('pydiffax is not available for this platform') |
---|
53 | import GSASIIfiles as G2fil |
---|
54 | |
---|
55 | # trig functions in degrees |
---|
56 | tand = lambda x: math.tan(x*math.pi/180.) |
---|
57 | atand = lambda x: 180.*math.atan(x)/math.pi |
---|
58 | atan2d = lambda y,x: 180.*math.atan2(y,x)/math.pi |
---|
59 | cosd = lambda x: math.cos(x*math.pi/180.) |
---|
60 | acosd = lambda x: 180.*math.acos(x)/math.pi |
---|
61 | rdsq2d = lambda x,p: round(1.0/math.sqrt(x),p) |
---|
62 | #numpy versions |
---|
63 | npsind = lambda x: np.sin(x*np.pi/180.) |
---|
64 | npasind = lambda x: 180.*np.arcsin(x)/math.pi |
---|
65 | npcosd = lambda x: np.cos(x*math.pi/180.) |
---|
66 | npacosd = lambda x: 180.*np.arccos(x)/math.pi |
---|
67 | nptand = lambda x: np.tan(x*math.pi/180.) |
---|
68 | npatand = lambda x: 180.*np.arctan(x)/np.pi |
---|
69 | npatan2d = lambda y,x: 180.*np.arctan2(y,x)/np.pi |
---|
70 | npT2stl = lambda tth, wave: 2.0*npsind(tth/2.0)/wave #=d* |
---|
71 | npT2q = lambda tth,wave: 2.0*np.pi*npT2stl(tth,wave) #=2pi*d* |
---|
72 | npq2T = lambda Q,wave: 2.0*npasind(0.25*Q*wave/np.pi) |
---|
73 | ateln2 = 8.0*math.log(2.0) |
---|
74 | sateln2 = np.sqrt(ateln2) |
---|
75 | nxs = np.newaxis |
---|
76 | is_exe = lambda fpath: os.path.isfile(fpath) and os.access(fpath, os.X_OK) |
---|
77 | |
---|
78 | #### Powder utilities ################################################################################ |
---|
79 | def PhaseWtSum(G2frame,histo): |
---|
80 | ''' |
---|
81 | Calculate sum of phase mass*phase fraction for PWDR data (exclude magnetic phases) |
---|
82 | |
---|
83 | :param G2frame: GSASII main frame structure |
---|
84 | :param str histo: histogram name |
---|
85 | :returns: sum(scale*mass) for phases in histo |
---|
86 | ''' |
---|
87 | Histograms,Phases = G2frame.GetUsedHistogramsAndPhasesfromTree() |
---|
88 | wtSum = 0.0 |
---|
89 | for phase in Phases: |
---|
90 | if Phases[phase]['General']['Type'] != 'magnetic': |
---|
91 | if histo in Phases[phase]['Histograms']: |
---|
92 | if not Phases[phase]['Histograms'][histo]['Use']: continue |
---|
93 | mass = Phases[phase]['General']['Mass'] |
---|
94 | phFr = Phases[phase]['Histograms'][histo]['Scale'][0] |
---|
95 | wtSum += mass*phFr |
---|
96 | return wtSum |
---|
97 | |
---|
98 | #### GSASII pwdr & pdf calculation routines ################################################################################ |
---|
99 | def Transmission(Geometry,Abs,Diam): |
---|
100 | ''' |
---|
101 | Calculate sample transmission |
---|
102 | |
---|
103 | :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting flat plate in transmission','Fixed flat plate' |
---|
104 | :param float Abs: absorption coeff in cm-1 |
---|
105 | :param float Diam: sample thickness/diameter in mm |
---|
106 | ''' |
---|
107 | if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample |
---|
108 | MuR = Abs*Diam/20.0 |
---|
109 | if MuR <= 3.0: |
---|
110 | T0 = 16/(3.*math.pi) |
---|
111 | T1 = -0.045780 |
---|
112 | T2 = -0.02489 |
---|
113 | T3 = 0.003045 |
---|
114 | T = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4 |
---|
115 | if T < -20.: |
---|
116 | return 2.06e-9 |
---|
117 | else: |
---|
118 | return math.exp(T) |
---|
119 | else: |
---|
120 | T1 = 1.433902 |
---|
121 | T2 = 0.013869+0.337894 |
---|
122 | T3 = 1.933433+1.163198 |
---|
123 | T4 = 0.044365-0.04259 |
---|
124 | T = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4 |
---|
125 | return T/100. |
---|
126 | elif 'plate' in Geometry: |
---|
127 | MuR = Abs*Diam/10. |
---|
128 | return math.exp(-MuR) |
---|
129 | elif 'Bragg' in Geometry: |
---|
130 | return 0.0 |
---|
131 | |
---|
132 | def SurfaceRough(SRA,SRB,Tth): |
---|
133 | ''' Suortti (J. Appl. Cryst, 5,325-331, 1972) surface roughness correction |
---|
134 | :param float SRA: Suortti surface roughness parameter |
---|
135 | :param float SRB: Suortti surface roughness parameter |
---|
136 | :param float Tth: 2-theta(deg) - can be numpy array |
---|
137 | |
---|
138 | ''' |
---|
139 | sth = npsind(Tth/2.) |
---|
140 | T1 = np.exp(-SRB/sth) |
---|
141 | T2 = SRA+(1.-SRA)*np.exp(-SRB) |
---|
142 | return (SRA+(1.-SRA)*T1)/T2 |
---|
143 | |
---|
144 | def SurfaceRoughDerv(SRA,SRB,Tth): |
---|
145 | ''' Suortti surface roughness correction derivatives |
---|
146 | :param float SRA: Suortti surface roughness parameter (dimensionless) |
---|
147 | :param float SRB: Suortti surface roughness parameter (dimensionless) |
---|
148 | :param float Tth: 2-theta(deg) - can be numpy array |
---|
149 | :return list: [dydSRA,dydSRB] derivatives to be used for intensity derivative |
---|
150 | ''' |
---|
151 | sth = npsind(Tth/2.) |
---|
152 | T1 = np.exp(-SRB/sth) |
---|
153 | T2 = SRA+(1.-SRA)*np.exp(-SRB) |
---|
154 | Trans = (SRA+(1.-SRA)*T1)/T2 |
---|
155 | dydSRA = ((1.-T1)*T2-(1.-np.exp(-SRB))*Trans)/T2**2 |
---|
156 | dydSRB = ((SRA-1.)*T1*T2/sth-Trans*(SRA-T2))/T2**2 |
---|
157 | return [dydSRA,dydSRB] |
---|
158 | |
---|
159 | def Absorb(Geometry,MuR,Tth,Phi=0,Psi=0): |
---|
160 | '''Calculate sample absorption |
---|
161 | :param str Geometry: one of 'Cylinder','Bragg-Brentano','Tilting Flat Plate in transmission','Fixed flat plate' |
---|
162 | :param float MuR: absorption coeff * sample thickness/2 or radius |
---|
163 | :param Tth: 2-theta scattering angle - can be numpy array |
---|
164 | :param float Phi: flat plate tilt angle - future |
---|
165 | :param float Psi: flat plate tilt axis - future |
---|
166 | ''' |
---|
167 | |
---|
168 | def muRunder3(MuR,Sth2): |
---|
169 | T0 = 16.0/(3.*np.pi) |
---|
170 | T1 = (25.99978-0.01911*Sth2**0.25)*np.exp(-0.024551*Sth2)+ \ |
---|
171 | 0.109561*np.sqrt(Sth2)-26.04556 |
---|
172 | T2 = -0.02489-0.39499*Sth2+1.219077*Sth2**1.5- \ |
---|
173 | 1.31268*Sth2**2+0.871081*Sth2**2.5-0.2327*Sth2**3 |
---|
174 | T3 = 0.003045+0.018167*Sth2-0.03305*Sth2**2 |
---|
175 | Trns = -T0*MuR-T1*MuR**2-T2*MuR**3-T3*MuR**4 |
---|
176 | return np.exp(Trns) |
---|
177 | |
---|
178 | def muRover3(MuR,Sth2): |
---|
179 | T1 = 1.433902+11.07504*Sth2-8.77629*Sth2*Sth2+ \ |
---|
180 | 10.02088*Sth2**3-3.36778*Sth2**4 |
---|
181 | T2 = (0.013869-0.01249*Sth2)*np.exp(3.27094*Sth2)+ \ |
---|
182 | (0.337894+13.77317*Sth2)/(1.0+11.53544*Sth2)**1.555039 |
---|
183 | T3 = 1.933433/(1.0+23.12967*Sth2)**1.686715- \ |
---|
184 | 0.13576*np.sqrt(Sth2)+1.163198 |
---|
185 | T4 = 0.044365-0.04259/(1.0+0.41051*Sth2)**148.4202 |
---|
186 | Trns = (T1-T4)/(1.0+T2*(MuR-3.0))**T3+T4 |
---|
187 | return Trns/100. |
---|
188 | |
---|
189 | Sth2 = npsind(Tth/2.0)**2 |
---|
190 | if 'Cylinder' in Geometry: #Lobanov & Alte da Veiga for 2-theta = 0; beam fully illuminates sample |
---|
191 | if 'array' in str(type(MuR)): |
---|
192 | MuRSTh2 = np.vstack((MuR,Sth2)) |
---|
193 | AbsCr = np.where(MuRSTh2[0]<=3.0,muRunder3(MuRSTh2[0],MuRSTh2[1]),muRover3(MuRSTh2[0],MuRSTh2[1])) |
---|
194 | return AbsCr |
---|
195 | else: |
---|
196 | if MuR <= 3.0: |
---|
197 | return muRunder3(MuR,Sth2) |
---|
198 | else: |
---|
199 | return muRover3(MuR,Sth2) |
---|
200 | elif 'Bragg' in Geometry: |
---|
201 | return 1.0 |
---|
202 | elif 'Fixed' in Geometry: #assumes sample plane is perpendicular to incident beam |
---|
203 | # and only defined for 2theta < 90 |
---|
204 | MuT = 2.*MuR |
---|
205 | T1 = np.exp(-MuT) |
---|
206 | T2 = np.exp(-MuT/npcosd(Tth)) |
---|
207 | Tb = MuT-MuT/npcosd(Tth) |
---|
208 | return (T2-T1)/Tb |
---|
209 | elif 'Tilting' in Geometry: #assumes symmetric tilt so sample plane is parallel to diffraction vector |
---|
210 | MuT = 2.*MuR |
---|
211 | cth = npcosd(Tth/2.0) |
---|
212 | return np.exp(-MuT/cth)/cth |
---|
213 | |
---|
214 | def AbsorbDerv(Geometry,MuR,Tth,Phi=0,Psi=0): |
---|
215 | 'needs a doc string' |
---|
216 | dA = 0.001 |
---|
217 | AbsP = Absorb(Geometry,MuR+dA,Tth,Phi,Psi) |
---|
218 | if MuR: |
---|
219 | AbsM = Absorb(Geometry,MuR-dA,Tth,Phi,Psi) |
---|
220 | return (AbsP-AbsM)/(2.0*dA) |
---|
221 | else: |
---|
222 | return (AbsP-1.)/dA |
---|
223 | |
---|
224 | def Polarization(Pola,Tth,Azm=0.0): |
---|
225 | """ Calculate angle dependent x-ray polarization correction (not scaled correctly!) |
---|
226 | |
---|
227 | :param Pola: polarization coefficient e.g 1.0 fully polarized, 0.5 unpolarized |
---|
228 | :param Azm: azimuthal angle e.g. 0.0 in plane of polarization - can be numpy array |
---|
229 | :param Tth: 2-theta scattering angle - can be numpy array |
---|
230 | which (if either) of these is "right"? |
---|
231 | :return: (pola, dpdPola) - both 2-d arrays |
---|
232 | * pola = ((1-Pola)*npcosd(Azm)**2+Pola*npsind(Azm)**2)*npcosd(Tth)**2+ \ |
---|
233 | (1-Pola)*npsind(Azm)**2+Pola*npcosd(Azm)**2 |
---|
234 | * dpdPola: derivative needed for least squares |
---|
235 | |
---|
236 | """ |
---|
237 | cazm = npcosd(Azm)**2 |
---|
238 | sazm = npsind(Azm)**2 |
---|
239 | pola = ((1.0-Pola)*cazm+Pola*sazm)*npcosd(Tth)**2+(1.0-Pola)*sazm+Pola*cazm |
---|
240 | dpdPola = -npsind(Tth)**2*(sazm-cazm) |
---|
241 | return pola,dpdPola |
---|
242 | |
---|
243 | def Oblique(ObCoeff,Tth): |
---|
244 | 'currently assumes detector is normal to beam' |
---|
245 | if ObCoeff: |
---|
246 | K = (1.-ObCoeff)/(1.0-np.exp(np.log(ObCoeff)/npcosd(Tth))) |
---|
247 | return K |
---|
248 | else: |
---|
249 | return 1.0 |
---|
250 | |
---|
251 | def Ruland(RulCoff,wave,Q,Compton): |
---|
252 | 'needs a doc string' |
---|
253 | C = 2.9978e8 |
---|
254 | D = 1.5e-3 |
---|
255 | hmc = 0.024262734687 #Compton wavelength in A |
---|
256 | sinth2 = (Q*wave/(4.0*np.pi))**2 |
---|
257 | dlam = (wave**2)*Compton*Q/C |
---|
258 | dlam_c = 2.0*hmc*sinth2-D*wave**2 |
---|
259 | return 1.0/((1.0+dlam/RulCoff)*(1.0+(np.pi*dlam_c/(dlam+RulCoff))**2)) |
---|
260 | |
---|
261 | def KleinNishina(wave,Q): |
---|
262 | hmc = 0.024262734687 #Compton wavelength in A |
---|
263 | TTh = npq2T(Q,wave) |
---|
264 | P = 1./(1.+(1.-npcosd(TTh)*(hmc/wave))) |
---|
265 | KN = (P**3-(P*npsind(TTh))**2+P)/(1.+npcosd(TTh)**2) |
---|
266 | return KN |
---|
267 | |
---|
268 | def LorchWeight(Q): |
---|
269 | 'needs a doc string' |
---|
270 | return np.sin(np.pi*(Q[-1]-Q)/(2.0*Q[-1])) |
---|
271 | |
---|
272 | def GetAsfMean(ElList,Sthl2): |
---|
273 | '''Calculate various scattering factor terms for PDF calcs |
---|
274 | |
---|
275 | :param dict ElList: element dictionary contains scattering factor coefficients, etc. |
---|
276 | :param np.array Sthl2: numpy array of sin theta/lambda squared values |
---|
277 | :returns: mean(f^2), mean(f)^2, mean(compton) |
---|
278 | ''' |
---|
279 | sumNoAtoms = 0.0 |
---|
280 | FF = np.zeros_like(Sthl2) |
---|
281 | FF2 = np.zeros_like(Sthl2) |
---|
282 | CF = np.zeros_like(Sthl2) |
---|
283 | for El in ElList: |
---|
284 | sumNoAtoms += ElList[El]['FormulaNo'] |
---|
285 | for El in ElList: |
---|
286 | el = ElList[El] |
---|
287 | ff2 = (G2elem.ScatFac(el,Sthl2)+el['fp'])**2+el['fpp']**2 |
---|
288 | cf = G2elem.ComptonFac(el,Sthl2) |
---|
289 | FF += np.sqrt(ff2)*el['FormulaNo']/sumNoAtoms |
---|
290 | FF2 += ff2*el['FormulaNo']/sumNoAtoms |
---|
291 | CF += cf*el['FormulaNo']/sumNoAtoms |
---|
292 | return FF2,FF**2,CF |
---|
293 | |
---|
294 | def GetNumDensity(ElList,Vol): |
---|
295 | 'needs a doc string' |
---|
296 | sumNoAtoms = 0.0 |
---|
297 | for El in ElList: |
---|
298 | sumNoAtoms += ElList[El]['FormulaNo'] |
---|
299 | return sumNoAtoms/Vol |
---|
300 | |
---|
301 | def CalcPDF(data,inst,limits,xydata): |
---|
302 | '''Computes I(Q), S(Q) & G(r) from Sample, Bkg, etc. diffraction patterns loaded into |
---|
303 | dict xydata; results are placed in xydata. |
---|
304 | Calculation parameters are found in dicts data and inst and list limits. |
---|
305 | The return value is at present an empty list. |
---|
306 | ''' |
---|
307 | auxPlot = [] |
---|
308 | if 'T' in inst['Type'][0]: |
---|
309 | Ibeg = 0 |
---|
310 | Ifin = len(xydata['Sample'][1][0]) |
---|
311 | else: |
---|
312 | Ibeg = np.searchsorted(xydata['Sample'][1][0],limits[0]) |
---|
313 | Ifin = np.searchsorted(xydata['Sample'][1][0],limits[1])+1 |
---|
314 | #subtract backgrounds - if any & use PWDR limits |
---|
315 | IofQ = copy.deepcopy(xydata['Sample']) |
---|
316 | IofQ[1] = np.array([I[Ibeg:Ifin] for I in IofQ[1]]) |
---|
317 | if data['Sample Bkg.']['Name']: |
---|
318 | IofQ[1][1] += xydata['Sample Bkg.'][1][1][Ibeg:Ifin]*data['Sample Bkg.']['Mult'] |
---|
319 | if data['Container']['Name']: |
---|
320 | xycontainer = xydata['Container'][1][1]*data['Container']['Mult'] |
---|
321 | if data['Container Bkg.']['Name']: |
---|
322 | xycontainer += xydata['Container Bkg.'][1][1][Ibeg:Ifin]*data['Container Bkg.']['Mult'] |
---|
323 | IofQ[1][1] += xycontainer[Ibeg:Ifin] |
---|
324 | data['IofQmin'] = IofQ[1][1][-1] |
---|
325 | IofQ[1][1] -= data.get('Flat Bkg',0.) |
---|
326 | #get element data & absorption coeff. |
---|
327 | ElList = data['ElList'] |
---|
328 | Tth = IofQ[1][0] #2-theta or TOF! |
---|
329 | if 'X' in inst['Type'][0]: |
---|
330 | Abs = G2lat.CellAbsorption(ElList,data['Form Vol']) |
---|
331 | #Apply angle dependent corrections |
---|
332 | MuR = Abs*data['Diam']/20.0 |
---|
333 | IofQ[1][1] /= Absorb(data['Geometry'],MuR,Tth) |
---|
334 | IofQ[1][1] /= Polarization(inst['Polariz.'][1],Tth,Azm=inst['Azimuth'][1])[0] |
---|
335 | if data['DetType'] == 'Area detector': |
---|
336 | IofQ[1][1] *= Oblique(data['ObliqCoeff'],Tth) |
---|
337 | elif 'T' in inst['Type'][0]: #neutron TOF normalized data - needs wavelength dependent absorption |
---|
338 | wave = 2.*G2lat.TOF2dsp(inst,IofQ[1][0])*npsind(inst['2-theta'][1]/2.) |
---|
339 | Els = ElList.keys() |
---|
340 | Isotope = {El:'Nat. abund.' for El in Els} |
---|
341 | GD = {'AtomTypes':ElList,'Isotope':Isotope} |
---|
342 | BLtables = G2elem.GetBLtable(GD) |
---|
343 | FP,FPP = G2elem.BlenResTOF(Els,BLtables,wave) |
---|
344 | Abs = np.zeros(len(wave)) |
---|
345 | for iel,El in enumerate(Els): |
---|
346 | BL = BLtables[El][1] |
---|
347 | SA = BL['SA']*wave/1.798197+4.0*np.pi*FPP[iel]**2 #+BL['SL'][1]? |
---|
348 | SA *= ElList[El]['FormulaNo']/data['Form Vol'] |
---|
349 | Abs += SA |
---|
350 | MuR = Abs*data['Diam']/2. |
---|
351 | IofQ[1][1] /= Absorb(data['Geometry'],MuR,inst['2-theta'][1]*np.ones(len(wave))) |
---|
352 | # improves look of F(Q) but no impact on G(R) |
---|
353 | # bBut,aBut = signal.butter(8,.5,"lowpass") |
---|
354 | # IofQ[1][1] = signal.filtfilt(bBut,aBut,IofQ[1][1]) |
---|
355 | XY = IofQ[1] |
---|
356 | #convert to Q |
---|
357 | nQpoints = 5000 |
---|
358 | if 'C' in inst['Type'][0]: |
---|
359 | wave = G2mth.getWave(inst) |
---|
360 | minQ = npT2q(Tth[0],wave) |
---|
361 | maxQ = npT2q(Tth[-1],wave) |
---|
362 | Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True) |
---|
363 | dq = Qpoints[1]-Qpoints[0] |
---|
364 | XY[0] = npT2q(XY[0],wave) |
---|
365 | Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][0]) #interpolate I(Q) |
---|
366 | elif 'T' in inst['Type'][0]: |
---|
367 | difC = inst['difC'][1] |
---|
368 | minQ = 2.*np.pi*difC/Tth[-1] |
---|
369 | maxQ = 2.*np.pi*difC/Tth[0] |
---|
370 | Qpoints = np.linspace(0.,maxQ,nQpoints,endpoint=True) |
---|
371 | dq = Qpoints[1]-Qpoints[0] |
---|
372 | XY[0] = 2.*np.pi*difC/XY[0] |
---|
373 | Qdata = si.griddata(XY[0],XY[1],Qpoints,method='linear',fill_value=XY[1][-1]) #interpolate I(Q) |
---|
374 | Qdata -= np.min(Qdata)*data['BackRatio'] |
---|
375 | |
---|
376 | qLimits = data['QScaleLim'] |
---|
377 | maxQ = np.searchsorted(Qpoints,min(Qpoints[-1],qLimits[1]))+1 |
---|
378 | minQ = np.searchsorted(Qpoints,min(qLimits[0],0.90*Qpoints[-1])) |
---|
379 | qLimits = [Qpoints[minQ],Qpoints[maxQ-1]] |
---|
380 | newdata = [] |
---|
381 | if len(IofQ) < 3: |
---|
382 | xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],''] |
---|
383 | else: |
---|
384 | xydata['IofQ'] = [IofQ[0],[Qpoints,Qdata],IofQ[2]] |
---|
385 | for item in xydata['IofQ'][1]: |
---|
386 | newdata.append(item[:maxQ]) |
---|
387 | xydata['IofQ'][1] = newdata |
---|
388 | |
---|
389 | xydata['SofQ'] = copy.deepcopy(xydata['IofQ']) |
---|
390 | if 'XC' in inst['Type'][0]: |
---|
391 | FFSq,SqFF,CF = GetAsfMean(ElList,(xydata['SofQ'][1][0]/(4.0*np.pi))**2) #these are <f^2>,<f>^2,Cf |
---|
392 | else: #TOF |
---|
393 | CF = np.zeros(len(xydata['SofQ'][1][0])) |
---|
394 | FFSq = np.ones(len(xydata['SofQ'][1][0])) |
---|
395 | SqFF = np.ones(len(xydata['SofQ'][1][0])) |
---|
396 | Q = xydata['SofQ'][1][0] |
---|
397 | # auxPlot.append([Q,np.copy(CF),'CF-unCorr']) |
---|
398 | if 'XC' in inst['Type'][0]: |
---|
399 | # CF *= KleinNishina(wave,Q) |
---|
400 | ruland = Ruland(data['Ruland'],wave,Q,CF) |
---|
401 | # auxPlot.append([Q,ruland,'Ruland']) |
---|
402 | CF *= ruland |
---|
403 | # auxPlot.append([Q,CF,'CF-Corr']) |
---|
404 | scale = np.sum((FFSq+CF)[minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ]) |
---|
405 | xydata['SofQ'][1][1] *= scale |
---|
406 | if 'XC' in inst['Type'][0]: |
---|
407 | xydata['SofQ'][1][1] -= CF |
---|
408 | xydata['SofQ'][1][1] = xydata['SofQ'][1][1]/SqFF |
---|
409 | scale = len(xydata['SofQ'][1][1][minQ:maxQ])/np.sum(xydata['SofQ'][1][1][minQ:maxQ]) |
---|
410 | xydata['SofQ'][1][1] *= scale |
---|
411 | xydata['FofQ'] = copy.deepcopy(xydata['SofQ']) |
---|
412 | xydata['FofQ'][1][1] = xydata['FofQ'][1][0]*(xydata['SofQ'][1][1]-1.0) |
---|
413 | if data['Lorch']: |
---|
414 | xydata['FofQ'][1][1] *= LorchWeight(Q) |
---|
415 | xydata['GofR'] = copy.deepcopy(xydata['FofQ']) |
---|
416 | xydata['gofr'] = copy.deepcopy(xydata['FofQ']) |
---|
417 | nR = len(xydata['GofR'][1][1]) |
---|
418 | Rmax = GSASIIpath.GetConfigValue('PDF_Rmax',100.) |
---|
419 | mul = int(round(2.*np.pi*nR/(Rmax*qLimits[1]))) |
---|
420 | # mul = int(round(2.*np.pi*nR/(data.get('Rmax',100.)*qLimits[1]))) |
---|
421 | R = 2.*np.pi*np.linspace(0,nR,nR,endpoint=True)/(mul*qLimits[1]) |
---|
422 | xydata['GofR'][1][0] = R |
---|
423 | xydata['gofr'][1][0] = R |
---|
424 | GR = -dq*np.imag(fft.fft(xydata['FofQ'][1][1],mul*nR)[:nR])*data.get('GR Scale',1.0) |
---|
425 | xydata['GofR'][1][1] = GR |
---|
426 | gr = GR/(np.pi*R) |
---|
427 | xydata['gofr'][1][1] = gr |
---|
428 | numbDen = 0. |
---|
429 | if 'ElList' in data: |
---|
430 | numbDen = GetNumDensity(data['ElList'],data['Form Vol']) |
---|
431 | if data.get('noRing',True): |
---|
432 | Rmin = data['Rmin'] |
---|
433 | xydata['gofr'][1][1] = np.where(R<Rmin,-4.*numbDen,xydata['gofr'][1][1]) |
---|
434 | xydata['GofR'][1][1] = np.where(R<Rmin,-4.*R*np.pi*numbDen,xydata['GofR'][1][1]) |
---|
435 | return auxPlot |
---|
436 | |
---|
437 | def PDFPeakFit(peaks,data): |
---|
438 | rs2pi = 1./np.sqrt(2*np.pi) |
---|
439 | |
---|
440 | def MakeParms(peaks): |
---|
441 | varyList = [] |
---|
442 | parmDict = {'slope':peaks['Background'][1][1]} |
---|
443 | if peaks['Background'][2]: |
---|
444 | varyList.append('slope') |
---|
445 | for i,peak in enumerate(peaks['Peaks']): |
---|
446 | parmDict['PDFpos;'+str(i)] = peak[0] |
---|
447 | parmDict['PDFmag;'+str(i)] = peak[1] |
---|
448 | parmDict['PDFsig;'+str(i)] = peak[2] |
---|
449 | if 'P' in peak[3]: |
---|
450 | varyList.append('PDFpos;'+str(i)) |
---|
451 | if 'M' in peak[3]: |
---|
452 | varyList.append('PDFmag;'+str(i)) |
---|
453 | if 'S' in peak[3]: |
---|
454 | varyList.append('PDFsig;'+str(i)) |
---|
455 | return parmDict,varyList |
---|
456 | |
---|
457 | def SetParms(peaks,parmDict,varyList): |
---|
458 | if 'slope' in varyList: |
---|
459 | peaks['Background'][1][1] = parmDict['slope'] |
---|
460 | for i,peak in enumerate(peaks['Peaks']): |
---|
461 | if 'PDFpos;'+str(i) in varyList: |
---|
462 | peak[0] = parmDict['PDFpos;'+str(i)] |
---|
463 | if 'PDFmag;'+str(i) in varyList: |
---|
464 | peak[1] = parmDict['PDFmag;'+str(i)] |
---|
465 | if 'PDFsig;'+str(i) in varyList: |
---|
466 | peak[2] = parmDict['PDFsig;'+str(i)] |
---|
467 | |
---|
468 | |
---|
469 | def CalcPDFpeaks(parmdict,Xdata): |
---|
470 | Z = parmDict['slope']*Xdata |
---|
471 | ipeak = 0 |
---|
472 | while True: |
---|
473 | try: |
---|
474 | pos = parmdict['PDFpos;'+str(ipeak)] |
---|
475 | mag = parmdict['PDFmag;'+str(ipeak)] |
---|
476 | wid = parmdict['PDFsig;'+str(ipeak)] |
---|
477 | wid2 = 2.*wid**2 |
---|
478 | Z += mag*rs2pi*np.exp(-(Xdata-pos)**2/wid2)/wid |
---|
479 | ipeak += 1 |
---|
480 | except KeyError: #no more peaks to process |
---|
481 | return Z |
---|
482 | |
---|
483 | def errPDFProfile(values,xdata,ydata,parmdict,varylist): |
---|
484 | parmdict.update(zip(varylist,values)) |
---|
485 | M = CalcPDFpeaks(parmdict,xdata)-ydata |
---|
486 | return M |
---|
487 | |
---|
488 | newpeaks = copy.copy(peaks) |
---|
489 | iBeg = np.searchsorted(data[1][0],newpeaks['Limits'][0]) |
---|
490 | iFin = np.searchsorted(data[1][0],newpeaks['Limits'][1])+1 |
---|
491 | X = data[1][0][iBeg:iFin] |
---|
492 | Y = data[1][1][iBeg:iFin] |
---|
493 | parmDict,varyList = MakeParms(peaks) |
---|
494 | if not len(varyList): |
---|
495 | G2fil.G2Print (' Nothing varied') |
---|
496 | return newpeaks,None,None,None,None,None |
---|
497 | |
---|
498 | Rvals = {} |
---|
499 | values = np.array(Dict2Values(parmDict, varyList)) |
---|
500 | result = so.leastsq(errPDFProfile,values,full_output=True,ftol=0.0001, |
---|
501 | args=(X,Y,parmDict,varyList)) |
---|
502 | chisq = np.sum(result[2]['fvec']**2) |
---|
503 | Values2Dict(parmDict, varyList, result[0]) |
---|
504 | SetParms(peaks,parmDict,varyList) |
---|
505 | Rvals['Rwp'] = np.sqrt(chisq/np.sum(Y**2))*100. #to % |
---|
506 | chisq = np.sum(result[2]['fvec']**2)/(len(X)-len(values)) #reduced chi^2 = M/(Nobs-Nvar) |
---|
507 | sigList = list(np.sqrt(chisq*np.diag(result[1]))) |
---|
508 | Z = CalcPDFpeaks(parmDict,X) |
---|
509 | newpeaks['calc'] = [X,Z] |
---|
510 | return newpeaks,result[0],varyList,sigList,parmDict,Rvals |
---|
511 | |
---|
512 | def MakeRDF(RDFcontrols,background,inst,pwddata): |
---|
513 | auxPlot = [] |
---|
514 | if 'C' in inst['Type'][0] or 'B' in inst['Type'][0]: |
---|
515 | Tth = pwddata[0] |
---|
516 | wave = G2mth.getWave(inst) |
---|
517 | minQ = npT2q(Tth[0],wave) |
---|
518 | maxQ = npT2q(Tth[-1],wave) |
---|
519 | powQ = npT2q(Tth,wave) |
---|
520 | elif 'T' in inst['Type'][0]: |
---|
521 | TOF = pwddata[0] |
---|
522 | difC = inst['difC'][1] |
---|
523 | minQ = 2.*np.pi*difC/TOF[-1] |
---|
524 | maxQ = 2.*np.pi*difC/TOF[0] |
---|
525 | powQ = 2.*np.pi*difC/TOF |
---|
526 | piDQ = np.pi/(maxQ-minQ) |
---|
527 | Qpoints = np.linspace(minQ,maxQ,len(pwddata[0]),endpoint=True) |
---|
528 | if RDFcontrols['UseObsCalc'] == 'obs-calc': |
---|
529 | Qdata = si.griddata(powQ,pwddata[1]-pwddata[3],Qpoints,method=RDFcontrols['Smooth'],fill_value=0.) |
---|
530 | elif RDFcontrols['UseObsCalc'] == 'obs-back': |
---|
531 | Qdata = si.griddata(powQ,pwddata[1]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0]) |
---|
532 | elif RDFcontrols['UseObsCalc'] == 'calc-back': |
---|
533 | Qdata = si.griddata(powQ,pwddata[3]-pwddata[4],Qpoints,method=RDFcontrols['Smooth'],fill_value=pwddata[1][0]) |
---|
534 | Qdata *= np.sin((Qpoints-minQ)*piDQ)/piDQ |
---|
535 | Qdata *= 0.5*np.sqrt(Qpoints) #Qbin normalization |
---|
536 | dq = Qpoints[1]-Qpoints[0] |
---|
537 | nR = len(Qdata) |
---|
538 | R = 0.5*np.pi*np.linspace(0,nR,nR)/(4.*maxQ) |
---|
539 | iFin = np.searchsorted(R,RDFcontrols['maxR'])+1 |
---|
540 | bBut,aBut = signal.butter(4,0.01) |
---|
541 | Qsmooth = signal.filtfilt(bBut,aBut,Qdata) |
---|
542 | # auxPlot.append([Qpoints,Qdata,'interpolate:'+RDFcontrols['Smooth']]) |
---|
543 | # auxPlot.append([Qpoints,Qsmooth,'interpolate:'+RDFcontrols['Smooth']]) |
---|
544 | DofR = dq*np.imag(fft.fft(Qsmooth,16*nR)[:nR]) |
---|
545 | auxPlot.append([R[:iFin],DofR[:iFin],'D(R) for '+RDFcontrols['UseObsCalc']]) |
---|
546 | return auxPlot |
---|
547 | |
---|
548 | # PDF optimization ============================================================= |
---|
549 | def OptimizePDF(data,xydata,limits,inst,showFit=True,maxCycles=25): |
---|
550 | import scipy.optimize as opt |
---|
551 | numbDen = GetNumDensity(data['ElList'],data['Form Vol']) |
---|
552 | Min,Init,Done = SetupPDFEval(data,xydata,limits,inst,numbDen) |
---|
553 | xstart = Init() |
---|
554 | bakMul = data['Sample Bkg.']['Mult'] |
---|
555 | if showFit: |
---|
556 | rms = Min(xstart) |
---|
557 | G2fil.G2Print(' Optimizing corrections to improve G(r) at low r') |
---|
558 | if data['Sample Bkg.'].get('Refine',False): |
---|
559 | # data['Flat Bkg'] = 0. |
---|
560 | G2fil.G2Print(' start: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f})'.format( |
---|
561 | data['Ruland'],data['Sample Bkg.']['Mult'],rms)) |
---|
562 | else: |
---|
563 | G2fil.G2Print(' start: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} (RMS:{:.4f})'.format( |
---|
564 | data['Flat Bkg'],data['BackRatio'],data['Ruland'],rms)) |
---|
565 | if data['Sample Bkg.'].get('Refine',False): |
---|
566 | res = opt.minimize(Min,xstart,bounds=([0.01,1.],[1.2*bakMul,0.8*bakMul]), |
---|
567 | method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001) |
---|
568 | else: |
---|
569 | res = opt.minimize(Min,xstart,bounds=([0.,None],[0,1],[0.01,1.]), |
---|
570 | method='L-BFGS-B',options={'maxiter':maxCycles},tol=0.001) |
---|
571 | Done(res['x']) |
---|
572 | if showFit: |
---|
573 | if res['success']: |
---|
574 | msg = 'Converged' |
---|
575 | else: |
---|
576 | msg = 'Not Converged' |
---|
577 | if data['Sample Bkg.'].get('Refine',False): |
---|
578 | G2fil.G2Print(' end: Ruland={:.3f}, Sample Bkg mult={:.3f} (RMS:{:.4f}) *** {} ***\n'.format( |
---|
579 | data['Ruland'],data['Sample Bkg.']['Mult'],res['fun'],msg)) |
---|
580 | else: |
---|
581 | G2fil.G2Print(' end: Flat Bkg={:.1f}, BackRatio={:.3f}, Ruland={:.3f} RMS:{:.4f}) *** {} ***\n'.format( |
---|
582 | data['Flat Bkg'],data['BackRatio'],data['Ruland'],res['fun'],msg)) |
---|
583 | return res |
---|
584 | |
---|
585 | def SetupPDFEval(data,xydata,limits,inst,numbDen): |
---|
586 | Data = copy.deepcopy(data) |
---|
587 | BkgMax = 1. |
---|
588 | def EvalLowPDF(arg): |
---|
589 | '''Objective routine -- evaluates the RMS deviations in G(r) |
---|
590 | from -4(pi)*#density*r for for r<Rmin |
---|
591 | arguments are ['Flat Bkg','BackRatio','Ruland'] scaled so that |
---|
592 | the min & max values are between 0 and 1. |
---|
593 | ''' |
---|
594 | if Data['Sample Bkg.'].get('Refine',False): |
---|
595 | R,S = arg |
---|
596 | Data['Sample Bkg.']['Mult'] = S |
---|
597 | else: |
---|
598 | F,B,R = arg |
---|
599 | Data['Flat Bkg'] = BkgMax*(2.*F-1.) |
---|
600 | Data['BackRatio'] = B |
---|
601 | Data['Ruland'] = R |
---|
602 | CalcPDF(Data,inst,limits,xydata) |
---|
603 | # test low r computation |
---|
604 | g = xydata['GofR'][1][1] |
---|
605 | r = xydata['GofR'][1][0] |
---|
606 | g0 = g[r < Data['Rmin']] + 4*np.pi*r[r < Data['Rmin']]*numbDen |
---|
607 | M = sum(g0**2)/len(g0) |
---|
608 | return M |
---|
609 | def GetCurrentVals(): |
---|
610 | '''Get the current ['Flat Bkg','BackRatio','Ruland'] with scaling |
---|
611 | ''' |
---|
612 | if data['Sample Bkg.'].get('Refine',False): |
---|
613 | return [max(data['Ruland'],.05),data['Sample']['Mult']] |
---|
614 | try: |
---|
615 | F = 0.5+0.5*data['Flat Bkg']/BkgMax |
---|
616 | except: |
---|
617 | F = 0 |
---|
618 | return [F,data['BackRatio'],max(data['Ruland'],.05)] |
---|
619 | def SetFinalVals(arg): |
---|
620 | '''Set the 'Flat Bkg', 'BackRatio' & 'Ruland' values from the |
---|
621 | scaled, refined values and plot corrected region of G(r) |
---|
622 | ''' |
---|
623 | if data['Sample Bkg.'].get('Refine',False): |
---|
624 | R,S = arg |
---|
625 | data['Sample Bkg.']['Mult'] = S |
---|
626 | else: |
---|
627 | F,B,R = arg |
---|
628 | data['Flat Bkg'] = BkgMax*(2.*F-1.) |
---|
629 | data['BackRatio'] = B |
---|
630 | data['Ruland'] = R |
---|
631 | CalcPDF(data,inst,limits,xydata) |
---|
632 | EvalLowPDF(GetCurrentVals()) |
---|
633 | BkgMax = max(xydata['IofQ'][1][1])/50. |
---|
634 | return EvalLowPDF,GetCurrentVals,SetFinalVals |
---|
635 | |
---|
636 | #### GSASII convolution peak fitting routines: Finger, Cox & Jephcoat model |
---|
637 | def factorize(num): |
---|
638 | ''' Provide prime number factors for integer num |
---|
639 | :returns: dictionary of prime factors (keys) & power for each (data) |
---|
640 | ''' |
---|
641 | factors = {} |
---|
642 | orig = num |
---|
643 | |
---|
644 | # we take advantage of the fact that (i +1)**2 = i**2 + 2*i +1 |
---|
645 | i, sqi = 2, 4 |
---|
646 | while sqi <= num: |
---|
647 | while not num%i: |
---|
648 | num /= i |
---|
649 | factors[i] = factors.get(i, 0) + 1 |
---|
650 | |
---|
651 | sqi += 2*i + 1 |
---|
652 | i += 1 |
---|
653 | |
---|
654 | if num != 1 and num != orig: |
---|
655 | factors[num] = factors.get(num, 0) + 1 |
---|
656 | |
---|
657 | if factors: |
---|
658 | return factors |
---|
659 | else: |
---|
660 | return {num:1} #a prime number! |
---|
661 | |
---|
662 | def makeFFTsizeList(nmin=1,nmax=1023,thresh=15): |
---|
663 | ''' Provide list of optimal data sizes for FFT calculations |
---|
664 | |
---|
665 | :param int nmin: minimum data size >= 1 |
---|
666 | :param int nmax: maximum data size > nmin |
---|
667 | :param int thresh: maximum prime factor allowed |
---|
668 | :Returns: list of data sizes where the maximum prime factor is < thresh |
---|
669 | ''' |
---|
670 | plist = [] |
---|
671 | nmin = max(1,nmin) |
---|
672 | nmax = max(nmin+1,nmax) |
---|
673 | for p in range(nmin,nmax): |
---|
674 | if max(list(factorize(p).keys())) < thresh: |
---|
675 | plist.append(p) |
---|
676 | return plist |
---|
677 | |
---|
678 | np.seterr(divide='ignore') |
---|
679 | |
---|
680 | # Normal distribution |
---|
681 | |
---|
682 | # loc = mu, scale = std |
---|
683 | _norm_pdf_C = 1./math.sqrt(2*math.pi) |
---|
684 | class norm_gen(st.rv_continuous): |
---|
685 | 'needs a doc string' |
---|
686 | |
---|
687 | def pdf(self,x,*args,**kwds): |
---|
688 | loc,scale=kwds['loc'],kwds['scale'] |
---|
689 | x = (x-loc)/scale |
---|
690 | return np.exp(-x**2/2.0) * _norm_pdf_C / scale |
---|
691 | |
---|
692 | norm = norm_gen(name='norm',longname='A normal',extradoc=""" |
---|
693 | |
---|
694 | Normal distribution |
---|
695 | |
---|
696 | The location (loc) keyword specifies the mean. |
---|
697 | The scale (scale) keyword specifies the standard deviation. |
---|
698 | |
---|
699 | normal.pdf(x) = exp(-x**2/2)/sqrt(2*pi) |
---|
700 | """) |
---|
701 | |
---|
702 | ## Cauchy |
---|
703 | |
---|
704 | # median = loc |
---|
705 | |
---|
706 | class cauchy_gen(st.rv_continuous): |
---|
707 | 'needs a doc string' |
---|
708 | |
---|
709 | def pdf(self,x,*args,**kwds): |
---|
710 | loc,scale=kwds['loc'],kwds['scale'] |
---|
711 | x = (x-loc)/scale |
---|
712 | return 1.0/np.pi/(1.0+x*x) / scale |
---|
713 | |
---|
714 | cauchy = cauchy_gen(name='cauchy',longname='Cauchy',extradoc=""" |
---|
715 | |
---|
716 | Cauchy distribution |
---|
717 | |
---|
718 | cauchy.pdf(x) = 1/(pi*(1+x**2)) |
---|
719 | |
---|
720 | This is the t distribution with one degree of freedom. |
---|
721 | """) |
---|
722 | |
---|
723 | |
---|
724 | class fcjde_gen(st.rv_continuous): |
---|
725 | """ |
---|
726 | Finger-Cox-Jephcoat D(2phi,2th) function for S/L = H/L |
---|
727 | Ref: J. Appl. Cryst. (1994) 27, 892-900. |
---|
728 | |
---|
729 | :param x: array -1 to 1 |
---|
730 | :param t: 2-theta position of peak |
---|
731 | :param s: sum(S/L,H/L); S: sample height, H: detector opening, |
---|
732 | L: sample to detector opening distance |
---|
733 | :param dx: 2-theta step size in deg |
---|
734 | |
---|
735 | :returns: for fcj.pdf |
---|
736 | |
---|
737 | * T = x*dx+t |
---|
738 | * s = S/L+H/L |
---|
739 | * if x < 0:: |
---|
740 | |
---|
741 | fcj.pdf = [1/sqrt({cos(T)**2/cos(t)**2}-1) - 1/s]/|cos(T)| |
---|
742 | |
---|
743 | * if x >= 0: fcj.pdf = 0 |
---|
744 | |
---|
745 | """ |
---|
746 | def _pdf(self,x,t,s,dx): |
---|
747 | T = dx*x+t |
---|
748 | ax2 = abs(npcosd(T)) |
---|
749 | ax = ax2**2 |
---|
750 | bx = npcosd(t)**2 |
---|
751 | bx = np.where(ax>bx,bx,ax) |
---|
752 | fx = np.where(ax>bx,(np.sqrt(bx/(ax-bx))-1./s)/ax2,0.0) |
---|
753 | fx = np.where(fx > 0.,fx,0.0) |
---|
754 | return fx |
---|
755 | |
---|
756 | def pdf(self,x,*args,**kwds): |
---|
757 | loc=kwds['loc'] |
---|
758 | return self._pdf(x-loc,*args) |
---|
759 | |
---|
760 | fcjde = fcjde_gen(name='fcjde',shapes='t,s,dx') |
---|
761 | |
---|
762 | def getFCJVoigt(pos,intens,sig,gam,shl,xdata): |
---|
763 | '''Compute the Finger-Cox-Jepcoat modified Voigt function for a |
---|
764 | CW powder peak by direct convolution. This version is not used. |
---|
765 | ''' |
---|
766 | DX = xdata[1]-xdata[0] |
---|
767 | widths,fmin,fmax = getWidthsCW(pos,sig,gam,shl) |
---|
768 | x = np.linspace(pos-fmin,pos+fmin,256) |
---|
769 | dx = x[1]-x[0] |
---|
770 | Norm = norm.pdf(x,loc=pos,scale=widths[0]) |
---|
771 | Cauchy = cauchy.pdf(x,loc=pos,scale=widths[1]) |
---|
772 | arg = [pos,shl/57.2958,dx,] |
---|
773 | FCJ = fcjde.pdf(x,*arg,loc=pos) |
---|
774 | if len(np.nonzero(FCJ)[0])>5: |
---|
775 | z = np.column_stack([Norm,Cauchy,FCJ]).T |
---|
776 | Z = fft.fft(z) |
---|
777 | Df = fft.ifft(Z.prod(axis=0)).real |
---|
778 | else: |
---|
779 | z = np.column_stack([Norm,Cauchy]).T |
---|
780 | Z = fft.fft(z) |
---|
781 | Df = fft.fftshift(fft.ifft(Z.prod(axis=0))).real |
---|
782 | Df /= np.sum(Df) |
---|
783 | Df = si.interp1d(x,Df,bounds_error=False,fill_value=0.0) |
---|
784 | return intens*Df(xdata)*DX/dx |
---|
785 | |
---|
786 | #### GSASII peak fitting routine: Finger, Cox & Jephcoat model |
---|
787 | |
---|
788 | def getWidthsCW(pos,sig,gam,shl): |
---|
789 | '''Compute the peak widths used for computing the range of a peak |
---|
790 | for constant wavelength data. On low-angle side, 50 FWHM are used, |
---|
791 | on high-angle side 75 are used, high angle side extended for axial divergence |
---|
792 | (for peaks above 90 deg, these are reversed.) |
---|
793 | |
---|
794 | :param pos: peak position; 2-theta in degrees |
---|
795 | :param sig: Gaussian peak variance in centideg^2 |
---|
796 | :param gam: Lorentzian peak width in centidegrees |
---|
797 | :param shl: axial divergence parameter (S+H)/L |
---|
798 | |
---|
799 | :returns: widths; [Gaussian sigma, Lorentzian gamma] in degrees, and |
---|
800 | low angle, high angle ends of peak; 20 FWHM & 50 FWHM from position |
---|
801 | reversed for 2-theta > 90 deg. |
---|
802 | ''' |
---|
803 | widths = [np.sqrt(sig)/100.,gam/100.] |
---|
804 | fwhm = 2.355*widths[0]+widths[1] |
---|
805 | fmin = 50.*(fwhm+shl*abs(npcosd(pos))) |
---|
806 | fmax = 75.0*fwhm |
---|
807 | if pos > 90: |
---|
808 | fmin,fmax = [fmax,fmin] |
---|
809 | return widths,fmin,fmax |
---|
810 | |
---|
811 | def getWidthsED(pos,sig): |
---|
812 | '''Compute the peak widths used for computing the range of a peak |
---|
813 | for energy dispersive data. On low-energy side, 20 FWHM are used, |
---|
814 | on high-energy side 20 are used |
---|
815 | |
---|
816 | :param pos: peak position; energy in keV (not used) |
---|
817 | :param sig: Gaussian peak variance in keV^2 |
---|
818 | |
---|
819 | :returns: widths; [Gaussian sigma] in keV, and |
---|
820 | low angle, high angle ends of peak; 20 FWHM & 50 FWHM from position |
---|
821 | ''' |
---|
822 | widths = [np.sqrt(sig),.001] |
---|
823 | fwhm = 2.355*widths[0] |
---|
824 | fmin = 5.*fwhm |
---|
825 | fmax = 5.*fwhm |
---|
826 | return widths,fmin,fmax |
---|
827 | |
---|
828 | def getWidthsTOF(pos,alp,bet,sig,gam): |
---|
829 | '''Compute the peak widths used for computing the range of a peak |
---|
830 | for constant wavelength data. 50 FWHM are used on both sides each |
---|
831 | extended by exponential coeff. |
---|
832 | |
---|
833 | param pos: peak position; TOF in musec (not used) |
---|
834 | param alp,bet: TOF peak exponential rise & decay parameters |
---|
835 | param sig: Gaussian peak variance in musec^2 |
---|
836 | param gam: Lorentzian peak width in musec |
---|
837 | |
---|
838 | returns: widths; [Gaussian sigma, Lornetzian gamma] in musec |
---|
839 | returns: low TOF, high TOF ends of peak; 50FWHM from position |
---|
840 | ''' |
---|
841 | widths = [np.sqrt(sig),gam] |
---|
842 | fwhm = 2.355*widths[0]+2.*widths[1] |
---|
843 | fmin = 50.*fwhm*(1.+1./alp) |
---|
844 | fmax = 50.*fwhm*(1.+1./bet) |
---|
845 | return widths,fmin,fmax |
---|
846 | |
---|
847 | def getFWHM(pos,Inst,N=1): |
---|
848 | '''Compute total FWHM from Thompson, Cox & Hastings (1987) , J. Appl. Cryst. 20, 79-83 |
---|
849 | via getgamFW(g,s). |
---|
850 | |
---|
851 | :param pos: float peak position in deg 2-theta or tof in musec |
---|
852 | :param Inst: dict instrument parameters |
---|
853 | :param N: int Inst index (0 for input, 1 for fitted) |
---|
854 | |
---|
855 | :returns float: total FWHM of pseudoVoigt in deg or musec |
---|
856 | ''' |
---|
857 | |
---|
858 | sig = lambda Th,U,V,W: np.sqrt(max(0.001,U*tand(Th)**2+V*tand(Th)+W)) |
---|
859 | sigED = lambda E,A,B,C: np.sqrt(max(0.001,A*E**2+B*E+C)) |
---|
860 | sigTOF = lambda dsp,S0,S1,S2,Sq: np.sqrt(S0+S1*dsp**2+S2*dsp**4+Sq*dsp) |
---|
861 | gam = lambda Th,X,Y,Z: Z+X/cosd(Th)+Y*tand(Th) |
---|
862 | gamTOF = lambda dsp,X,Y,Z: Z+X*dsp+Y*dsp**2 |
---|
863 | alpTOF = lambda dsp,alp: alp/dsp |
---|
864 | betTOF = lambda dsp,bet0,bet1,betq: bet0+bet1/dsp**4+betq/dsp**2 |
---|
865 | alpPink = lambda pos,alp0,alp1: alp0+alp1*tand(pos/2.) |
---|
866 | betPink = lambda pos,bet0,bet1: bet0+bet1*tand(pos/2.) |
---|
867 | if 'LF' in Inst['Type'][0]: |
---|
868 | return 3 |
---|
869 | elif 'T' in Inst['Type'][0]: |
---|
870 | dsp = pos/Inst['difC'][N] |
---|
871 | alp = alpTOF(dsp,Inst['alpha'][N]) |
---|
872 | bet = betTOF(dsp,Inst['beta-0'][1],Inst['beta-1'][N],Inst['beta-q'][N]) |
---|
873 | s = sigTOF(dsp,Inst['sig-0'][N],Inst['sig-1'][N],Inst['sig-2'][N],Inst['sig-q'][N]) |
---|
874 | g = gamTOF(dsp,Inst['X'][N],Inst['Y'][N],Inst['Z'][N]) |
---|
875 | return getgamFW(g,s)+np.log(2.0)*(alp+bet)/(alp*bet) |
---|
876 | elif 'C' in Inst['Type'][0]: |
---|
877 | s = sig(pos/2.,Inst['U'][N],Inst['V'][N],Inst['W'][N]) |
---|
878 | g = gam(pos/2.,Inst['X'][N],Inst['Y'][N],Inst['Z'][N]) |
---|
879 | return getgamFW(g,s)/100. #returns FWHM in deg |
---|
880 | elif 'E' in Inst['Type'][0]: |
---|
881 | s = sigED(pos,Inst['A'][N],Inst['B'][N],Inst['C'][N]) |
---|
882 | return 2.35482*s |
---|
883 | else: #'B' |
---|
884 | alp = alpPink(pos,Inst['alpha-0'][N],Inst['alpha-1'][N]) |
---|
885 | bet = betPink(pos,Inst['beta-0'][N],Inst['beta-1'][N]) |
---|
886 | s = sig(pos/2.,Inst['U'][N],Inst['V'][N],Inst['W'][N]) |
---|
887 | g = gam(pos/2.,Inst['X'][N],Inst['Y'][N],Inst['Z'][N]) |
---|
888 | return getgamFW(g,s)/100.+np.log(2.0)*(alp+bet)/(alp*bet) #returns FWHM in deg |
---|
889 | |
---|
890 | def getgamFW(g,s): |
---|
891 | '''Compute total FWHM from Thompson, Cox & Hastings (1987), J. Appl. Cryst. 20, 79-83 |
---|
892 | lambda fxn needs FWHM for both Gaussian & Lorentzian components |
---|
893 | |
---|
894 | :param g: float Lorentzian gamma = FWHM(L) |
---|
895 | :param s: float Gaussian sig |
---|
896 | |
---|
897 | :returns float: total FWHM of pseudoVoigt |
---|
898 | ''' |
---|
899 | gamFW = lambda s,g: np.exp(np.log(s**5+2.69269*s**4*g+2.42843*s**3*g**2+4.47163*s**2*g**3+0.07842*s*g**4+g**5)/5.) |
---|
900 | return gamFW(2.35482*s,g) #sqrt(8ln2)*sig = FWHM(G) |
---|
901 | |
---|
902 | def getBackground(pfx,parmDict,bakType,dataType,xdata,fixback=None): |
---|
903 | '''Computes the background from vars pulled from gpx file or tree. |
---|
904 | ''' |
---|
905 | if 'T' in dataType: |
---|
906 | q = 2.*np.pi*parmDict[pfx+'difC']/xdata |
---|
907 | elif 'E' in dataType: |
---|
908 | const = 4.*np.pi*npsind(parmDict[pfx+'2-theta']/2.0) |
---|
909 | q = const*xdata |
---|
910 | else: |
---|
911 | wave = parmDict.get(pfx+'Lam',parmDict.get(pfx+'Lam1',1.0)) |
---|
912 | q = npT2q(xdata,wave) |
---|
913 | yb = np.zeros_like(xdata) |
---|
914 | nBak = 0 |
---|
915 | sumBk = [0.,0.,0] |
---|
916 | while True: |
---|
917 | key = pfx+'Back;'+str(nBak) |
---|
918 | if key in parmDict: |
---|
919 | nBak += 1 |
---|
920 | else: |
---|
921 | break |
---|
922 | #empirical functions |
---|
923 | if bakType in ['chebyschev','cosine','chebyschev-1']: |
---|
924 | dt = xdata[-1]-xdata[0] |
---|
925 | for iBak in range(nBak): |
---|
926 | key = pfx+'Back;'+str(iBak) |
---|
927 | if bakType == 'chebyschev': |
---|
928 | ybi = parmDict[key]*(-1.+2.*(xdata-xdata[0])/dt)**iBak |
---|
929 | elif bakType == 'chebyschev-1': |
---|
930 | xpos = -1.+2.*(xdata-xdata[0])/dt |
---|
931 | ybi = parmDict[key]*np.cos(iBak*np.arccos(xpos)) |
---|
932 | elif bakType == 'cosine': |
---|
933 | ybi = parmDict[key]*npcosd(180.*xdata*iBak/xdata[-1]) |
---|
934 | yb += ybi |
---|
935 | sumBk[0] = np.sum(yb) |
---|
936 | elif bakType in ['Q^2 power series','Q^-2 power series']: |
---|
937 | QT = 1. |
---|
938 | yb += np.ones_like(yb)*parmDict[pfx+'Back;0'] |
---|
939 | for iBak in range(nBak-1): |
---|
940 | key = pfx+'Back;'+str(iBak+1) |
---|
941 | if '-2' in bakType: |
---|
942 | QT *= (iBak+1)*q**-2 |
---|
943 | else: |
---|
944 | QT *= q**2/(iBak+1) |
---|
945 | yb += QT*parmDict[key] |
---|
946 | sumBk[0] = np.sum(yb) |
---|
947 | elif bakType in ['lin interpolate','inv interpolate','log interpolate',]: |
---|
948 | if nBak == 1: |
---|
949 | yb = np.ones_like(xdata)*parmDict[pfx+'Back;0'] |
---|
950 | elif nBak == 2: |
---|
951 | dX = xdata[-1]-xdata[0] |
---|
952 | T2 = (xdata-xdata[0])/dX |
---|
953 | T1 = 1.0-T2 |
---|
954 | yb = parmDict[pfx+'Back;0']*T1+parmDict[pfx+'Back;1']*T2 |
---|
955 | else: |
---|
956 | xnomask = ma.getdata(xdata) |
---|
957 | xmin,xmax = xnomask[0],xnomask[-1] |
---|
958 | if bakType == 'lin interpolate': |
---|
959 | bakPos = np.linspace(xmin,xmax,nBak,True) |
---|
960 | elif bakType == 'inv interpolate': |
---|
961 | bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True) |
---|
962 | elif bakType == 'log interpolate': |
---|
963 | bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True)) |
---|
964 | bakPos[0] = xmin |
---|
965 | bakPos[-1] = xmax |
---|
966 | bakVals = np.zeros(nBak) |
---|
967 | for i in range(nBak): |
---|
968 | bakVals[i] = parmDict[pfx+'Back;'+str(i)] |
---|
969 | bakInt = si.interp1d(bakPos,bakVals,'linear') |
---|
970 | yb = bakInt(ma.getdata(xdata)) |
---|
971 | sumBk[0] = np.sum(yb) |
---|
972 | #Debye function |
---|
973 | if pfx+'difC' in parmDict or 'E' in dataType: |
---|
974 | ff = 1. |
---|
975 | else: |
---|
976 | try: |
---|
977 | wave = parmDict[pfx+'Lam'] |
---|
978 | except KeyError: |
---|
979 | wave = parmDict[pfx+'Lam1'] |
---|
980 | SQ = (q/(4.*np.pi))**2 |
---|
981 | FF = G2elem.GetFormFactorCoeff('Si')[0] |
---|
982 | ff = np.array(G2elem.ScatFac(FF,SQ)[0])**2 |
---|
983 | iD = 0 |
---|
984 | while True: |
---|
985 | try: |
---|
986 | dbA = parmDict[pfx+'DebyeA;'+str(iD)] |
---|
987 | dbR = parmDict[pfx+'DebyeR;'+str(iD)] |
---|
988 | dbU = parmDict[pfx+'DebyeU;'+str(iD)] |
---|
989 | ybi = ff*dbA*np.sin(q*dbR)*np.exp(-dbU*q**2)/(q*dbR) |
---|
990 | yb += ybi |
---|
991 | sumBk[1] += np.sum(ybi) |
---|
992 | iD += 1 |
---|
993 | except KeyError: |
---|
994 | break |
---|
995 | #peaks |
---|
996 | iD = 0 |
---|
997 | while True: |
---|
998 | try: |
---|
999 | pkP = parmDict[pfx+'BkPkpos;'+str(iD)] |
---|
1000 | pkI = max(parmDict[pfx+'BkPkint;'+str(iD)],0.1) |
---|
1001 | pkS = max(parmDict[pfx+'BkPksig;'+str(iD)],0.01) |
---|
1002 | pkG = max(parmDict[pfx+'BkPkgam;'+str(iD)],0.1) |
---|
1003 | if 'C' in dataType: |
---|
1004 | Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002) |
---|
1005 | elif 'E' in dataType: |
---|
1006 | Wd,fmin,fmax = getWidthsED(pkP,pkS) |
---|
1007 | else: #'T'OF |
---|
1008 | Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG) |
---|
1009 | iBeg = np.searchsorted(xdata,pkP-fmin) |
---|
1010 | iFin = np.searchsorted(xdata,pkP+fmax) |
---|
1011 | lenX = len(xdata) |
---|
1012 | if not iBeg: |
---|
1013 | iFin = np.searchsorted(xdata,pkP+fmax) |
---|
1014 | elif iBeg == lenX: |
---|
1015 | iFin = iBeg |
---|
1016 | else: |
---|
1017 | iFin = np.searchsorted(xdata,pkP+fmax) |
---|
1018 | if 'C' in dataType: |
---|
1019 | ybi = pkI*getFCJVoigt3(pkP,pkS,pkG,0.002,xdata[iBeg:iFin])[0] |
---|
1020 | elif 'T' in dataType: |
---|
1021 | ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin])[0] |
---|
1022 | elif 'B' in dataType: |
---|
1023 | ybi = pkI*getEpsVoigt(pkP,1.,1.,pkS/100.,pkG/1.e4,xdata[iBeg:iFin])[0] |
---|
1024 | elif 'E' in dataType: |
---|
1025 | ybi = pkI*getPsVoigt(pkP,pkS*10.**4,pkG*100.,xdata[iBeg:iFin])[0] |
---|
1026 | else: |
---|
1027 | raise Exception('dataType of {:} should not happen!'.format(dataType)) |
---|
1028 | yb[iBeg:iFin] += ybi |
---|
1029 | sumBk[2] += np.sum(ybi) |
---|
1030 | iD += 1 |
---|
1031 | except KeyError: |
---|
1032 | break |
---|
1033 | except ValueError: |
---|
1034 | G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****') |
---|
1035 | break |
---|
1036 | if fixback is not None: |
---|
1037 | yb += parmDict[pfx+'BF mult']*fixback |
---|
1038 | sumBk[0] = sum(yb) |
---|
1039 | return yb,sumBk |
---|
1040 | |
---|
1041 | def getBackgroundDerv(hfx,parmDict,bakType,dataType,xdata,fixback=None): |
---|
1042 | 'needs a doc string' |
---|
1043 | if 'T' in dataType: |
---|
1044 | q = 2.*np.pi*parmDict[hfx+'difC']/xdata |
---|
1045 | elif 'E' in dataType: |
---|
1046 | const = 4.*np.pi*npsind(parmDict[hfx+'2-theta']/2.0) |
---|
1047 | q = const*xdata |
---|
1048 | else: |
---|
1049 | wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0)) |
---|
1050 | q = 2.*np.pi*npsind(xdata/2.)/wave |
---|
1051 | nBak = 0 |
---|
1052 | while True: |
---|
1053 | key = hfx+'Back;'+str(nBak) |
---|
1054 | if key in parmDict: |
---|
1055 | nBak += 1 |
---|
1056 | else: |
---|
1057 | break |
---|
1058 | dydb = np.zeros(shape=(nBak,len(xdata))) |
---|
1059 | dyddb = np.zeros(shape=(3*parmDict[hfx+'nDebye'],len(xdata))) |
---|
1060 | dydpk = np.zeros(shape=(4*parmDict[hfx+'nPeaks'],len(xdata))) |
---|
1061 | dydfb = [] |
---|
1062 | |
---|
1063 | if bakType in ['chebyschev','cosine','chebyschev-1']: |
---|
1064 | dt = xdata[-1]-xdata[0] |
---|
1065 | for iBak in range(nBak): |
---|
1066 | if bakType == 'chebyschev': |
---|
1067 | dydb[iBak] = (-1.+2.*(xdata-xdata[0])/dt)**iBak |
---|
1068 | elif bakType == 'chebyschev-1': |
---|
1069 | xpos = -1.+2.*(xdata-xdata[0])/dt |
---|
1070 | dydb[iBak] = np.cos(iBak*np.arccos(xpos)) |
---|
1071 | elif bakType == 'cosine': |
---|
1072 | dydb[iBak] = npcosd(180.*xdata*iBak/xdata[-1]) |
---|
1073 | elif bakType in ['Q^2 power series','Q^-2 power series']: |
---|
1074 | QT = 1. |
---|
1075 | dydb[0] = np.ones_like(xdata) |
---|
1076 | for iBak in range(nBak-1): |
---|
1077 | if '-2' in bakType: |
---|
1078 | QT *= (iBak+1)*q**-2 |
---|
1079 | else: |
---|
1080 | QT *= q**2/(iBak+1) |
---|
1081 | dydb[iBak+1] = QT |
---|
1082 | elif bakType in ['lin interpolate','inv interpolate','log interpolate',]: |
---|
1083 | if nBak == 1: |
---|
1084 | dydb[0] = np.ones_like(xdata) |
---|
1085 | elif nBak == 2: |
---|
1086 | dX = xdata[-1]-xdata[0] |
---|
1087 | T2 = (xdata-xdata[0])/dX |
---|
1088 | T1 = 1.0-T2 |
---|
1089 | dydb = [T1,T2] |
---|
1090 | else: |
---|
1091 | xnomask = ma.getdata(xdata) |
---|
1092 | xmin,xmax = xnomask[0],xnomask[-1] |
---|
1093 | if bakType == 'lin interpolate': |
---|
1094 | bakPos = np.linspace(xmin,xmax,nBak,True) |
---|
1095 | elif bakType == 'inv interpolate': |
---|
1096 | bakPos = 1./np.linspace(1./xmax,1./xmin,nBak,True) |
---|
1097 | elif bakType == 'log interpolate': |
---|
1098 | bakPos = np.exp(np.linspace(np.log(xmin),np.log(xmax),nBak,True)) |
---|
1099 | bakPos[0] = xmin |
---|
1100 | bakPos[-1] = xmax |
---|
1101 | for i,pos in enumerate(bakPos): |
---|
1102 | if i == 0: |
---|
1103 | dydb[0] = np.where(xdata<bakPos[1],(bakPos[1]-xdata)/(bakPos[1]-bakPos[0]),0.) |
---|
1104 | elif i == len(bakPos)-1: |
---|
1105 | dydb[i] = np.where(xdata>bakPos[-2],(bakPos[-1]-xdata)/(bakPos[-1]-bakPos[-2]),0.) |
---|
1106 | else: |
---|
1107 | dydb[i] = np.where(xdata>bakPos[i], |
---|
1108 | np.where(xdata<bakPos[i+1],(bakPos[i+1]-xdata)/(bakPos[i+1]-bakPos[i]),0.), |
---|
1109 | np.where(xdata>bakPos[i-1],(xdata-bakPos[i-1])/(bakPos[i]-bakPos[i-1]),0.)) |
---|
1110 | if hfx+'difC' in parmDict: |
---|
1111 | ff = 1. |
---|
1112 | else: |
---|
1113 | wave = parmDict.get(hfx+'Lam',parmDict.get(hfx+'Lam1',1.0)) |
---|
1114 | q = npT2q(xdata,wave) |
---|
1115 | SQ = (q/(4*np.pi))**2 |
---|
1116 | FF = G2elem.GetFormFactorCoeff('Si')[0] |
---|
1117 | ff = np.array(G2elem.ScatFac(FF,SQ)[0])*np.pi**2 #needs pi^2~10. for cw data (why?) |
---|
1118 | iD = 0 |
---|
1119 | while True: |
---|
1120 | try: |
---|
1121 | if hfx+'difC' in parmDict: |
---|
1122 | q = 2*np.pi*parmDict[hfx+'difC']/xdata |
---|
1123 | dbA = parmDict[hfx+'DebyeA;'+str(iD)] |
---|
1124 | dbR = parmDict[hfx+'DebyeR;'+str(iD)] |
---|
1125 | dbU = parmDict[hfx+'DebyeU;'+str(iD)] |
---|
1126 | sqr = np.sin(q*dbR)/(q*dbR) |
---|
1127 | cqr = np.cos(q*dbR) |
---|
1128 | temp = np.exp(-dbU*q**2) |
---|
1129 | dyddb[3*iD] = ff*sqr*temp |
---|
1130 | dyddb[3*iD+1] = ff*dbA*temp*(cqr-sqr)/(dbR) |
---|
1131 | dyddb[3*iD+2] = -ff*dbA*sqr*temp*q**2 |
---|
1132 | iD += 1 |
---|
1133 | except KeyError: |
---|
1134 | break |
---|
1135 | iD = 0 |
---|
1136 | while True: |
---|
1137 | try: |
---|
1138 | pkP = parmDict[hfx+'BkPkpos;'+str(iD)] |
---|
1139 | pkI = max(parmDict[hfx+'BkPkint;'+str(iD)],0.1) |
---|
1140 | pkS = max(parmDict[hfx+'BkPksig;'+str(iD)],0.01) |
---|
1141 | pkG = max(parmDict[hfx+'BkPkgam;'+str(iD)],0.1) |
---|
1142 | if 'C' in dataType: |
---|
1143 | Wd,fmin,fmax = getWidthsCW(pkP,pkS,pkG,.002) |
---|
1144 | elif 'E' in dataType: |
---|
1145 | Wd,fmin,fmax = getWidthsED(pkP,pkS) |
---|
1146 | else: #'T' or 'B' |
---|
1147 | Wd,fmin,fmax = getWidthsTOF(pkP,1.,1.,pkS,pkG) |
---|
1148 | iBeg = np.searchsorted(xdata,pkP-fmin) |
---|
1149 | iFin = np.searchsorted(xdata,pkP+fmax) |
---|
1150 | lenX = len(xdata) |
---|
1151 | if not iBeg: |
---|
1152 | iFin = np.searchsorted(xdata,pkP+fmax) |
---|
1153 | elif iBeg == lenX: |
---|
1154 | iFin = iBeg |
---|
1155 | else: |
---|
1156 | iFin = np.searchsorted(xdata,pkP+fmax) |
---|
1157 | if 'C' in dataType: |
---|
1158 | Df,dFdp,dFds,dFdg,x = getdFCJVoigt3(pkP,pkS,pkG,.002,xdata[iBeg:iFin]) |
---|
1159 | elif 'E' in dataType: |
---|
1160 | Df,dFdp,dFds,dFdg = getdPsVoigt(pkP,pkS*10.**4,pkG*100.,xdata[iBeg:iFin]) |
---|
1161 | else: #'T'OF |
---|
1162 | Df,dFdp,x,x,dFds,dFdg = getdEpsVoigt(pkP,1.,1.,pkS,pkG,xdata[iBeg:iFin]) |
---|
1163 | dydpk[4*iD][iBeg:iFin] += pkI*dFdp |
---|
1164 | dydpk[4*iD+1][iBeg:iFin] += Df |
---|
1165 | dydpk[4*iD+2][iBeg:iFin] += pkI*dFds |
---|
1166 | dydpk[4*iD+3][iBeg:iFin] += pkI*dFdg |
---|
1167 | iD += 1 |
---|
1168 | except KeyError: |
---|
1169 | break |
---|
1170 | except ValueError: |
---|
1171 | G2fil.G2Print ('**** WARNING - backround peak '+str(iD)+' sigma is negative; fix & try again ****') |
---|
1172 | break |
---|
1173 | # fixed background from file |
---|
1174 | if fixback is not None: |
---|
1175 | dydfb = fixback |
---|
1176 | return dydb,dyddb,dydpk,dydfb |
---|
1177 | |
---|
1178 | #### Using old gsas fortran routines for powder peak shapes & derivatives |
---|
1179 | def getFCJVoigt3(pos,sig,gam,shl,xdata): |
---|
1180 | '''Compute the Finger-Cox-Jepcoat modified Pseudo-Voigt function for a |
---|
1181 | CW powder peak in external Fortran routine |
---|
1182 | |
---|
1183 | param pos: peak position in degrees |
---|
1184 | param sig: Gaussian variance in centideg^2 |
---|
1185 | param gam: Lorentzian width in centideg |
---|
1186 | param shl: axial divergence parameter (S+H)/L |
---|
1187 | param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv) |
---|
1188 | |
---|
1189 | returns: array: calculated peak function at each xdata |
---|
1190 | returns: integral of peak; nominally = 1.0 |
---|
1191 | ''' |
---|
1192 | if len(xdata): |
---|
1193 | cw = np.diff(xdata) |
---|
1194 | cw = np.append(cw,cw[-1]) |
---|
1195 | Df = pyd.pypsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) |
---|
1196 | return Df,np.sum(100.*Df*cw) |
---|
1197 | else: |
---|
1198 | return 0.,1. |
---|
1199 | |
---|
1200 | def getdFCJVoigt3(pos,sig,gam,shl,xdata): |
---|
1201 | '''Compute analytic derivatives the Finger-Cox-Jepcoat modified Pseudo-Voigt |
---|
1202 | function for a CW powder peak |
---|
1203 | |
---|
1204 | param pos: peak position in degrees |
---|
1205 | param sig: Gaussian variance in centideg^2 |
---|
1206 | param gam: Lorentzian width in centideg |
---|
1207 | param shl: axial divergence parameter (S+H)/L |
---|
1208 | param xdata: array; profile points for peak to be calculated; bounded by 20FWHM to 50FWHM (or vv) |
---|
1209 | |
---|
1210 | returns: arrays: function and derivatives of pos, sig, gam, & shl |
---|
1211 | ''' |
---|
1212 | Df,dFdp,dFds,dFdg,dFdsh = pyd.pydpsvfcj(len(xdata),xdata-pos,pos,sig,gam,shl) |
---|
1213 | return Df,dFdp,dFds,dFdg,dFdsh |
---|
1214 | |
---|
1215 | def getPsVoigt(pos,sig,gam,xdata): |
---|
1216 | '''Compute the simple Pseudo-Voigt function for a |
---|
1217 | small angle Bragg peak in external Fortran routine |
---|
1218 | |
---|
1219 | param pos: peak position in degrees |
---|
1220 | param sig: Gaussian variance in centideg^2 |
---|
1221 | param gam: Lorentzian width in centideg |
---|
1222 | |
---|
1223 | returns: array: calculated peak function at each xdata |
---|
1224 | returns: integral of peak; nominally = 1.0 |
---|
1225 | ''' |
---|
1226 | |
---|
1227 | cw = np.diff(xdata) |
---|
1228 | cw = np.append(cw,cw[-1]) |
---|
1229 | Df = pyd.pypsvoigt(len(xdata),xdata-pos,sig,gam) |
---|
1230 | return Df,np.sum(100.*Df*cw) |
---|
1231 | |
---|
1232 | def getdPsVoigt(pos,sig,gam,xdata): |
---|
1233 | '''Compute the simple Pseudo-Voigt function derivatives for a |
---|
1234 | small angle Bragg peak peak in external Fortran routine |
---|
1235 | |
---|
1236 | param pos: peak position in degrees |
---|
1237 | param sig: Gaussian variance in centideg^2 |
---|
1238 | param gam: Lorentzian width in centideg |
---|
1239 | |
---|
1240 | returns: arrays: function and derivatives of pos, sig & gam |
---|
1241 | NB: the pos derivative has the opposite sign compared to that in other profile functions |
---|
1242 | ''' |
---|
1243 | |
---|
1244 | Df,dFdp,dFds,dFdg = pyd.pydpsvoigt(len(xdata),xdata-pos,sig,gam) |
---|
1245 | return Df,dFdp,dFds,dFdg |
---|
1246 | |
---|
1247 | def getEpsVoigt(pos,alp,bet,sig,gam,xdata): |
---|
1248 | '''Compute the double exponential Pseudo-Voigt convolution function for a |
---|
1249 | neutron TOF & CW pink peak in external Fortran routine |
---|
1250 | ''' |
---|
1251 | |
---|
1252 | cw = np.diff(xdata) |
---|
1253 | cw = np.append(cw,cw[-1]) |
---|
1254 | Df = pyd.pyepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam) |
---|
1255 | return Df,np.sum(Df*cw) |
---|
1256 | |
---|
1257 | def getdEpsVoigt(pos,alp,bet,sig,gam,xdata): |
---|
1258 | '''Compute the double exponential Pseudo-Voigt convolution function derivatives for a |
---|
1259 | neutron TOF & CW pink peak in external Fortran routine |
---|
1260 | ''' |
---|
1261 | |
---|
1262 | Df,dFdp,dFda,dFdb,dFds,dFdg = pyd.pydepsvoigt(len(xdata),xdata-pos,alp,bet,sig,gam) |
---|
1263 | return Df,dFdp,dFda,dFdb,dFds,dFdg |
---|
1264 | |
---|
1265 | def ellipseSize(H,Sij,GB): |
---|
1266 | '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) per note from Alexander Brady |
---|
1267 | ''' |
---|
1268 | |
---|
1269 | HX = np.inner(H.T,GB) |
---|
1270 | lenHX = np.sqrt(np.sum(HX**2)) |
---|
1271 | Esize,Rsize = nl.eigh(G2lat.U6toUij(Sij)) |
---|
1272 | R = np.inner(HX/lenHX,Rsize)**2*Esize #want column length for hkl in crystal |
---|
1273 | lenR = 1./np.sqrt(np.sum(R)) |
---|
1274 | return lenR |
---|
1275 | |
---|
1276 | def ellipseSizeDerv(H,Sij,GB): |
---|
1277 | '''Implements r=1/sqrt(sum((1/S)*(q.v)^2) derivative per note from Alexander Brady |
---|
1278 | ''' |
---|
1279 | |
---|
1280 | lenR = ellipseSize(H,Sij,GB) |
---|
1281 | delt = 0.001 |
---|
1282 | dRdS = np.zeros(6) |
---|
1283 | for i in range(6): |
---|
1284 | Sij[i] -= delt |
---|
1285 | lenM = ellipseSize(H,Sij,GB) |
---|
1286 | Sij[i] += 2.*delt |
---|
1287 | lenP = ellipseSize(H,Sij,GB) |
---|
1288 | Sij[i] -= delt |
---|
1289 | dRdS[i] = (lenP-lenM)/(2.*delt) |
---|
1290 | return lenR,dRdS |
---|
1291 | |
---|
1292 | def getMustrain(HKL,G,SGData,muStrData): |
---|
1293 | if muStrData[0] == 'isotropic': |
---|
1294 | return np.ones(HKL.shape[1])*muStrData[1][0] |
---|
1295 | elif muStrData[0] == 'uniaxial': |
---|
1296 | H = np.array(HKL) |
---|
1297 | P = np.array(muStrData[3]) |
---|
1298 | cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T |
---|
1299 | Si = muStrData[1][0] |
---|
1300 | Sa = muStrData[1][1] |
---|
1301 | return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2)) |
---|
1302 | else: #generalized - P.W. Stephens model |
---|
1303 | H = np.array(HKL) |
---|
1304 | rdsq = np.array([G2lat.calc_rDsq2(h,G) for h in H.T]) |
---|
1305 | Strms = np.array(G2spc.MustrainCoeff(H,SGData)) |
---|
1306 | Sum = np.sum(np.array(muStrData[4])[:,nxs]*Strms,axis=0) |
---|
1307 | return np.sqrt(Sum)/rdsq |
---|
1308 | |
---|
1309 | def getCrSize(HKL,G,GB,sizeData): |
---|
1310 | if sizeData[0] == 'isotropic': |
---|
1311 | return np.ones(HKL.shape[1])*sizeData[1][0] |
---|
1312 | elif sizeData[0] == 'uniaxial': |
---|
1313 | H = np.array(HKL) |
---|
1314 | P = np.array(sizeData[3]) |
---|
1315 | cosP,sinP = np.array([G2lat.CosSinAngle(h,P,G) for h in H.T]).T |
---|
1316 | Si = sizeData[1][0] |
---|
1317 | Sa = sizeData[1][1] |
---|
1318 | return Si*Sa/(np.sqrt((Si*cosP)**2+(Sa*sinP)**2)) |
---|
1319 | else: |
---|
1320 | Sij =[sizeData[4][i] for i in range(6)] |
---|
1321 | H = np.array(HKL) |
---|
1322 | return 1./np.array([ellipseSize(h,Sij,GB) for h in H.T])**2 |
---|
1323 | |
---|
1324 | def getHKLpeak(dmin,SGData,A,Inst=None,nodup=False): |
---|
1325 | ''' |
---|
1326 | Generates allowed by symmetry reflections with d >= dmin |
---|
1327 | NB: GenHKLf & checkMagextc return True for extinct reflections |
---|
1328 | |
---|
1329 | :param dmin: minimum d-spacing |
---|
1330 | :param SGData: space group data obtained from SpcGroup |
---|
1331 | :param A: lattice parameter terms A1-A6 |
---|
1332 | :param Inst: instrument parameter info |
---|
1333 | :returns: HKLs: np.array hkl, etc for allowed reflections |
---|
1334 | |
---|
1335 | ''' |
---|
1336 | HKL = G2lat.GenHLaue(dmin,SGData,A) |
---|
1337 | HKLs = [] |
---|
1338 | ds = [] |
---|
1339 | for h,k,l,d in HKL: |
---|
1340 | ext = G2spc.GenHKLf([h,k,l],SGData)[0] |
---|
1341 | if ext and 'MagSpGrp' in SGData: |
---|
1342 | ext = G2spc.checkMagextc([h,k,l],SGData) |
---|
1343 | if not ext: |
---|
1344 | if nodup and int(10000*d) in ds: |
---|
1345 | continue |
---|
1346 | ds.append(int(10000*d)) |
---|
1347 | if Inst == None: |
---|
1348 | HKLs.append([h,k,l,d,0,-1]) |
---|
1349 | else: |
---|
1350 | HKLs.append([h,k,l,d,G2lat.Dsp2pos(Inst,d),-1]) |
---|
1351 | return np.array(HKLs) |
---|
1352 | |
---|
1353 | def getHKLMpeak(dmin,Inst,SGData,SSGData,Vec,maxH,A): |
---|
1354 | 'needs a doc string' |
---|
1355 | HKLs = [] |
---|
1356 | vec = np.array(Vec) |
---|
1357 | vstar = np.sqrt(G2lat.calc_rDsq(vec,A)) #find extra needed for -n SS reflections |
---|
1358 | dvec = 1./(maxH*vstar+1./dmin) |
---|
1359 | HKL = G2lat.GenHLaue(dvec,SGData,A) |
---|
1360 | SSdH = [vec*h for h in range(-maxH,maxH+1)] |
---|
1361 | SSdH = dict(zip(range(-maxH,maxH+1),SSdH)) |
---|
1362 | ifMag = False |
---|
1363 | if 'MagSpGrp' in SGData: |
---|
1364 | ifMag = True |
---|
1365 | for h,k,l,d in HKL: |
---|
1366 | ext = G2spc.GenHKLf([h,k,l],SGData)[0] |
---|
1367 | if not ext and d >= dmin: |
---|
1368 | HKLs.append([h,k,l,0,d,G2lat.Dsp2pos(Inst,d),-1]) |
---|
1369 | for dH in SSdH: |
---|
1370 | if dH: |
---|
1371 | DH = SSdH[dH] |
---|
1372 | H = [h+DH[0],k+DH[1],l+DH[2]] |
---|
1373 | d = float(1/np.sqrt(G2lat.calc_rDsq(H,A))) |
---|
1374 | if d >= dmin: |
---|
1375 | HKLM = np.array([h,k,l,dH]) |
---|
1376 | if G2spc.checkSSextc(HKLM,SSGData) or ifMag: |
---|
1377 | HKLs.append([h,k,l,dH,d,G2lat.Dsp2pos(Inst,d),-1]) |
---|
1378 | return G2lat.sortHKLd(HKLs,True,True,True) |
---|
1379 | |
---|
1380 | peakInstPrmMode = True |
---|
1381 | '''Determines the mode used for peak fitting. When peakInstPrmMode=True peak |
---|
1382 | width parameters are computed from the instrument parameters (UVW,... or |
---|
1383 | alpha,... etc) unless the individual parameter is refined. This allows the |
---|
1384 | instrument parameters to be refined. When peakInstPrmMode=False, the instrument |
---|
1385 | parameters are not used and cannot be refined. |
---|
1386 | The default is peakFitMode=True. This is changed only in |
---|
1387 | :func:`setPeakInstPrmMode`, which is called from :mod:`GSASIIscriptable` |
---|
1388 | or GSASIIphsGUI.OnSetPeakWidMode ('Gen unvaried widths' menu item). |
---|
1389 | ''' |
---|
1390 | |
---|
1391 | def setPeakInstPrmMode(normal=True): |
---|
1392 | '''Determines the mode used for peak fitting. If normal=True (default) |
---|
1393 | peak width parameters are computed from the instrument parameters |
---|
1394 | unless the individual parameter is refined. If normal=False, |
---|
1395 | peak widths are used as supplied for each peak. |
---|
1396 | |
---|
1397 | Note that normal=True unless this routine is called. Also, |
---|
1398 | instrument parameters can only be refined with normal=True. |
---|
1399 | |
---|
1400 | :param bool normal: setting to apply to global variable |
---|
1401 | :data:`peakInstPrmMode` |
---|
1402 | ''' |
---|
1403 | global peakInstPrmMode |
---|
1404 | peakInstPrmMode = normal |
---|
1405 | |
---|
1406 | def getPeakProfile(dataType,parmDict,xdata,fixback,varyList,bakType): |
---|
1407 | '''Computes the profiles from multiple peaks for individual peak fitting |
---|
1408 | for powder patterns. |
---|
1409 | NB: not used for Rietveld refinement |
---|
1410 | ''' |
---|
1411 | |
---|
1412 | yb = getBackground('',parmDict,bakType,dataType,xdata,fixback)[0] |
---|
1413 | yc = np.zeros_like(yb) |
---|
1414 | if 'LF' in dataType: |
---|
1415 | if 'Lam1' in parmDict.keys(): |
---|
1416 | lam = parmDict['Lam1'] |
---|
1417 | lam2 = parmDict['Lam2'] |
---|
1418 | Ka2 = True |
---|
1419 | lamRatio = 360*(lam2-lam)/(np.pi*lam) |
---|
1420 | kRatio = parmDict['I(L2)/I(L1)'] |
---|
1421 | else: |
---|
1422 | lam = parmDict['Lam'] |
---|
1423 | Ka2 = False |
---|
1424 | shol = 0 |
---|
1425 | # loop over peaks |
---|
1426 | iPeak = -1 |
---|
1427 | try: |
---|
1428 | ncells = parmDict['ncell'] |
---|
1429 | clat = parmDict['clat'] |
---|
1430 | except KeyError: # no Laue info must be bkg fit |
---|
1431 | print('Laue Fit: no params, assuming bkg fit') |
---|
1432 | return yb |
---|
1433 | while True: |
---|
1434 | iPeak += 1 |
---|
1435 | try: |
---|
1436 | #Qcen = 2 * np.pi * lam * (iPeak+1) / parmDict['clat'] |
---|
1437 | l = parmDict['l'+str(iPeak)] |
---|
1438 | pos = 360 * np.arcsin(0.5 * lam * l / parmDict['clat']) / np.pi |
---|
1439 | parmDict['pos'+str(iPeak)] = pos |
---|
1440 | #tth = (pos-parmDict['Zero']) |
---|
1441 | intens = parmDict['int'+str(iPeak)] |
---|
1442 | damp = parmDict['damp'+str(iPeak)] |
---|
1443 | asym = parmDict['asym'+str(iPeak)] |
---|
1444 | sig = parmDict['sig'+str(iPeak)] |
---|
1445 | gam = parmDict['gam'+str(iPeak)] |
---|
1446 | fmin = 8 # for now make peaks 8 degrees wide |
---|
1447 | fmin = min(0.9*abs(xdata[-1] - xdata[0]),fmin) # unless the data range is smaller |
---|
1448 | iBeg = np.searchsorted(xdata,pos-fmin/2) |
---|
1449 | iFin = np.searchsorted(xdata,pos+fmin/2) |
---|
1450 | if not iBeg+iFin: # skip peak below low limit |
---|
1451 | continue |
---|
1452 | elif not iBeg-iFin: # got peak above high limit (peaks sorted, so we can stop) |
---|
1453 | break |
---|
1454 | #LF.plotme(fmin,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym) |
---|
1455 | #LaueFringePeakCalc(xdata,yc,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym,fmin,plot=(iPeak==0)) |
---|
1456 | LaueFringePeakCalc(xdata,yc,lam,pos,intens,sig,gam,shol,ncells,clat,damp,asym,fmin,plot=False) |
---|
1457 | if Ka2: |
---|
1458 | pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) |
---|
1459 | iBeg = np.searchsorted(xdata,pos2-fmin) |
---|
1460 | iFin = np.searchsorted(xdata,pos2+fmin) |
---|
1461 | if iBeg-iFin: |
---|
1462 | LaueFringePeakCalc(xdata,yc,lam2,pos2,intens*kRatio,sig,gam,shol,ncells,clat,damp,asym,fmin) |
---|
1463 | except KeyError: #no more peaks to process |
---|
1464 | return yb+yc |
---|
1465 | elif 'C' in dataType: |
---|
1466 | shl = max(parmDict['SH/L'],0.002) |
---|
1467 | Ka2 = False |
---|
1468 | if 'Lam1' in parmDict.keys(): |
---|
1469 | Ka2 = True |
---|
1470 | lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1']) |
---|
1471 | kRatio = parmDict['I(L2)/I(L1)'] |
---|
1472 | iPeak = 0 |
---|
1473 | while True: |
---|
1474 | try: |
---|
1475 | pos = parmDict['pos'+str(iPeak)] |
---|
1476 | tth = (pos-parmDict['Zero']) |
---|
1477 | intens = parmDict['int'+str(iPeak)] |
---|
1478 | sigName = 'sig'+str(iPeak) |
---|
1479 | if sigName in varyList or not peakInstPrmMode: |
---|
1480 | sig = parmDict[sigName] |
---|
1481 | else: |
---|
1482 | sig = G2mth.getCWsig(parmDict,tth) |
---|
1483 | sig = max(sig,0.001) #avoid neg sigma^2 |
---|
1484 | gamName = 'gam'+str(iPeak) |
---|
1485 | if gamName in varyList or not peakInstPrmMode: |
---|
1486 | gam = parmDict[gamName] |
---|
1487 | else: |
---|
1488 | gam = G2mth.getCWgam(parmDict,tth) |
---|
1489 | gam = max(gam,0.001) #avoid neg gamma |
---|
1490 | Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl) |
---|
1491 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1492 | iFin = np.searchsorted(xdata,pos+fmin) |
---|
1493 | if not iBeg+iFin: #peak below low limit |
---|
1494 | iPeak += 1 |
---|
1495 | continue |
---|
1496 | elif not iBeg-iFin: #peak above high limit |
---|
1497 | return yb+yc |
---|
1498 | fp = getFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin])[0] |
---|
1499 | yc[iBeg:iFin] += intens*fp |
---|
1500 | if Ka2: |
---|
1501 | pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) |
---|
1502 | iBeg = np.searchsorted(xdata,pos2-fmin) |
---|
1503 | iFin = np.searchsorted(xdata,pos2+fmin) |
---|
1504 | if iBeg-iFin: |
---|
1505 | fp2 = getFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin])[0] |
---|
1506 | yc[iBeg:iFin] += intens*kRatio*fp2 |
---|
1507 | iPeak += 1 |
---|
1508 | except KeyError: #no more peaks to process |
---|
1509 | return yb+yc |
---|
1510 | elif 'E' in dataType: |
---|
1511 | iPeak = 0 |
---|
1512 | dsp = 1.0 #for now - fix later |
---|
1513 | while True: |
---|
1514 | try: |
---|
1515 | pos = parmDict['pos'+str(iPeak)] |
---|
1516 | intens = parmDict['int'+str(iPeak)] |
---|
1517 | sigName = 'sig'+str(iPeak) |
---|
1518 | if sigName in varyList or not peakInstPrmMode: |
---|
1519 | sig = parmDict[sigName] |
---|
1520 | else: |
---|
1521 | sig = G2mth.getEDsig(parmDict,pos) |
---|
1522 | sig = max(sig,0.001) #avoid neg sigma^2 |
---|
1523 | Wd,fmin,fmax = getWidthsED(pos,sig) |
---|
1524 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1525 | iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin)) |
---|
1526 | if not iBeg+iFin: #peak below low limit |
---|
1527 | iPeak += 1 |
---|
1528 | continue |
---|
1529 | elif not iBeg-iFin: #peak above high limit |
---|
1530 | return yb+yc |
---|
1531 | yc[iBeg:iFin] += intens*getPsVoigt(pos,sig*10.**4,0.001,xdata[iBeg:iFin])[0] |
---|
1532 | iPeak += 1 |
---|
1533 | except KeyError: #no more peaks to process |
---|
1534 | return yb+yc |
---|
1535 | elif 'B' in dataType: |
---|
1536 | iPeak = 0 |
---|
1537 | dsp = 1.0 #for now - fix later |
---|
1538 | while True: |
---|
1539 | try: |
---|
1540 | pos = parmDict['pos'+str(iPeak)] |
---|
1541 | tth = (pos-parmDict['Zero']) |
---|
1542 | intens = parmDict['int'+str(iPeak)] |
---|
1543 | alpName = 'alp'+str(iPeak) |
---|
1544 | if alpName in varyList or not peakInstPrmMode: |
---|
1545 | alp = parmDict[alpName] |
---|
1546 | else: |
---|
1547 | alp = G2mth.getPinkalpha(parmDict,tth) |
---|
1548 | alp = max(0.1,alp) |
---|
1549 | betName = 'bet'+str(iPeak) |
---|
1550 | if betName in varyList or not peakInstPrmMode: |
---|
1551 | bet = parmDict[betName] |
---|
1552 | else: |
---|
1553 | bet = G2mth.getPinkbeta(parmDict,tth) |
---|
1554 | bet = max(0.1,bet) |
---|
1555 | sigName = 'sig'+str(iPeak) |
---|
1556 | if sigName in varyList or not peakInstPrmMode: |
---|
1557 | sig = parmDict[sigName] |
---|
1558 | else: |
---|
1559 | sig = G2mth.getCWsig(parmDict,tth) |
---|
1560 | sig = max(sig,0.001) #avoid neg sigma^2 |
---|
1561 | gamName = 'gam'+str(iPeak) |
---|
1562 | if gamName in varyList or not peakInstPrmMode: |
---|
1563 | gam = parmDict[gamName] |
---|
1564 | else: |
---|
1565 | gam = G2mth.getCWgam(parmDict,tth) |
---|
1566 | gam = max(gam,0.001) #avoid neg gamma |
---|
1567 | Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam) |
---|
1568 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1569 | iFin = np.searchsorted(xdata,pos+fmin) |
---|
1570 | if not iBeg+iFin: #peak below low limit |
---|
1571 | iPeak += 1 |
---|
1572 | continue |
---|
1573 | elif not iBeg-iFin: #peak above high limit |
---|
1574 | return yb+yc |
---|
1575 | yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin])[0] |
---|
1576 | iPeak += 1 |
---|
1577 | except KeyError: #no more peaks to process |
---|
1578 | return yb+yc |
---|
1579 | else: |
---|
1580 | Pdabc = parmDict['Pdabc'] |
---|
1581 | difC = parmDict['difC'] |
---|
1582 | iPeak = 0 |
---|
1583 | while True: |
---|
1584 | try: |
---|
1585 | pos = parmDict['pos'+str(iPeak)] |
---|
1586 | tof = pos-parmDict['Zero'] |
---|
1587 | dsp = tof/difC |
---|
1588 | intens = parmDict['int'+str(iPeak)] |
---|
1589 | alpName = 'alp'+str(iPeak) |
---|
1590 | if alpName in varyList or not peakInstPrmMode: |
---|
1591 | alp = parmDict[alpName] |
---|
1592 | else: |
---|
1593 | if len(Pdabc): |
---|
1594 | alp = np.interp(dsp,Pdabc[0],Pdabc[1]) |
---|
1595 | else: |
---|
1596 | alp = G2mth.getTOFalpha(parmDict,dsp) |
---|
1597 | alp = max(0.1,alp) |
---|
1598 | betName = 'bet'+str(iPeak) |
---|
1599 | if betName in varyList or not peakInstPrmMode: |
---|
1600 | bet = parmDict[betName] |
---|
1601 | else: |
---|
1602 | if len(Pdabc): |
---|
1603 | bet = np.interp(dsp,Pdabc[0],Pdabc[2]) |
---|
1604 | else: |
---|
1605 | bet = G2mth.getTOFbeta(parmDict,dsp) |
---|
1606 | bet = max(0.0001,bet) |
---|
1607 | sigName = 'sig'+str(iPeak) |
---|
1608 | if sigName in varyList or not peakInstPrmMode: |
---|
1609 | sig = parmDict[sigName] |
---|
1610 | else: |
---|
1611 | sig = G2mth.getTOFsig(parmDict,dsp) |
---|
1612 | gamName = 'gam'+str(iPeak) |
---|
1613 | if gamName in varyList or not peakInstPrmMode: |
---|
1614 | gam = parmDict[gamName] |
---|
1615 | else: |
---|
1616 | gam = G2mth.getTOFgamma(parmDict,dsp) |
---|
1617 | gam = max(gam,0.001) #avoid neg gamma |
---|
1618 | Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam) |
---|
1619 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1620 | iFin = np.searchsorted(xdata,pos+fmax) |
---|
1621 | lenX = len(xdata) |
---|
1622 | if not iBeg: |
---|
1623 | iFin = np.searchsorted(xdata,pos+fmax) |
---|
1624 | elif iBeg == lenX: |
---|
1625 | iFin = iBeg |
---|
1626 | else: |
---|
1627 | iFin = np.searchsorted(xdata,pos+fmax) |
---|
1628 | if not iBeg+iFin: #peak below low limit |
---|
1629 | iPeak += 1 |
---|
1630 | continue |
---|
1631 | elif not iBeg-iFin: #peak above high limit |
---|
1632 | return yb+yc |
---|
1633 | yc[iBeg:iFin] += intens*getEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin])[0] |
---|
1634 | iPeak += 1 |
---|
1635 | except KeyError: #no more peaks to process |
---|
1636 | return yb+yc |
---|
1637 | |
---|
1638 | def getPeakProfileDerv(dataType,parmDict,xdata,fixback,varyList,bakType): |
---|
1639 | '''Computes the profile derivatives for a powder pattern for single peak fitting |
---|
1640 | |
---|
1641 | return: np.array([dMdx1,dMdx2,...]) in same order as varylist = backVary,insVary,peakVary order |
---|
1642 | |
---|
1643 | NB: not used for Rietveld refinement |
---|
1644 | ''' |
---|
1645 | dMdv = np.zeros(shape=(len(varyList),len(xdata))) |
---|
1646 | dMdb,dMddb,dMdpk,dMdfb = getBackgroundDerv('',parmDict,bakType,dataType,xdata,fixback) |
---|
1647 | if 'Back;0' in varyList: #background derivs are in front if present |
---|
1648 | dMdv[0:len(dMdb)] = dMdb |
---|
1649 | names = ['DebyeA','DebyeR','DebyeU'] |
---|
1650 | for name in varyList: |
---|
1651 | if 'Debye' in name: |
---|
1652 | parm,Id = name.split(';') |
---|
1653 | ip = names.index(parm) |
---|
1654 | dMdv[varyList.index(name)] = dMddb[3*int(Id)+ip] |
---|
1655 | names = ['BkPkpos','BkPkint','BkPksig','BkPkgam'] |
---|
1656 | for name in varyList: |
---|
1657 | if 'BkPk' in name: |
---|
1658 | parm,Id = name.split(';') |
---|
1659 | ip = names.index(parm) |
---|
1660 | dMdv[varyList.index(name)] = dMdpk[4*int(Id)+ip] |
---|
1661 | if 'LF' in dataType: |
---|
1662 | for i,name in enumerate(varyList): |
---|
1663 | if not np.all(dMdv[i] == 0): continue |
---|
1664 | deltaParmDict = parmDict.copy() |
---|
1665 | delta = max(parmDict[name]/1e5,0.001) |
---|
1666 | deltaParmDict[name] += delta |
---|
1667 | #print('num. deriv for',name,'val',deltaParmDict[name],'delta',delta) |
---|
1668 | intArrP = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType) |
---|
1669 | deltaParmDict[name] -= 2*delta |
---|
1670 | intArrM = getPeakProfile(dataType,deltaParmDict,xdata,fixback,varyList,bakType) |
---|
1671 | dMdv[i] = 0.5 * (intArrP - intArrM) / delta |
---|
1672 | return dMdv |
---|
1673 | if 'C' in dataType: |
---|
1674 | shl = max(parmDict['SH/L'],0.002) |
---|
1675 | Ka2 = False |
---|
1676 | if 'Lam1' in parmDict.keys(): |
---|
1677 | Ka2 = True |
---|
1678 | lamRatio = 360*(parmDict['Lam2']-parmDict['Lam1'])/(np.pi*parmDict['Lam1']) |
---|
1679 | kRatio = parmDict['I(L2)/I(L1)'] |
---|
1680 | iPeak = 0 |
---|
1681 | while True: |
---|
1682 | try: |
---|
1683 | pos = parmDict['pos'+str(iPeak)] |
---|
1684 | tth = (pos-parmDict['Zero']) |
---|
1685 | intens = parmDict['int'+str(iPeak)] |
---|
1686 | sigName = 'sig'+str(iPeak) |
---|
1687 | if sigName in varyList or not peakInstPrmMode: |
---|
1688 | sig = parmDict[sigName] |
---|
1689 | dsdU = dsdV = dsdW = 0 |
---|
1690 | else: |
---|
1691 | sig = G2mth.getCWsig(parmDict,tth) |
---|
1692 | dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth) |
---|
1693 | sig = max(sig,0.001) #avoid neg sigma |
---|
1694 | gamName = 'gam'+str(iPeak) |
---|
1695 | if gamName in varyList or not peakInstPrmMode: |
---|
1696 | gam = parmDict[gamName] |
---|
1697 | dgdX = dgdY = dgdZ = 0 |
---|
1698 | else: |
---|
1699 | gam = G2mth.getCWgam(parmDict,tth) |
---|
1700 | dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth) |
---|
1701 | gam = max(gam,0.001) #avoid neg gamma |
---|
1702 | Wd,fmin,fmax = getWidthsCW(pos,sig,gam,shl) |
---|
1703 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1704 | iFin = max(iBeg+3,np.searchsorted(xdata,pos+fmin)) |
---|
1705 | if not iBeg+iFin: #peak below low limit |
---|
1706 | iPeak += 1 |
---|
1707 | continue |
---|
1708 | elif not iBeg-iFin: #peak above high limit |
---|
1709 | break |
---|
1710 | dMdpk = np.zeros(shape=(6,len(xdata))) |
---|
1711 | dMdipk = getdFCJVoigt3(pos,sig,gam,shl,xdata[iBeg:iFin]) |
---|
1712 | for i in range(1,5): |
---|
1713 | dMdpk[i][iBeg:iFin] += intens*dMdipk[i] |
---|
1714 | dMdpk[0][iBeg:iFin] += dMdipk[0] |
---|
1715 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4]} |
---|
1716 | if Ka2: |
---|
1717 | pos2 = pos+lamRatio*tand(pos/2.0) # + 360/pi * Dlam/lam * tan(th) |
---|
1718 | iBeg = np.searchsorted(xdata,pos2-fmin) |
---|
1719 | iFin = np.searchsorted(xdata,pos2+fmin) |
---|
1720 | if iBeg-iFin: |
---|
1721 | dMdipk2 = getdFCJVoigt3(pos2,sig,gam,shl,xdata[iBeg:iFin]) |
---|
1722 | for i in range(1,5): |
---|
1723 | dMdpk[i][iBeg:iFin] += intens*kRatio*dMdipk2[i] |
---|
1724 | dMdpk[0][iBeg:iFin] += kRatio*dMdipk2[0] |
---|
1725 | dMdpk[5][iBeg:iFin] += dMdipk2[0] |
---|
1726 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'sig':dMdpk[2],'gam':dMdpk[3],'shl':dMdpk[4],'L1/L2':dMdpk[5]*intens} |
---|
1727 | for parmName in ['pos','int','sig','gam']: |
---|
1728 | try: |
---|
1729 | idx = varyList.index(parmName+str(iPeak)) |
---|
1730 | dMdv[idx] = dervDict[parmName] |
---|
1731 | except ValueError: |
---|
1732 | pass |
---|
1733 | if 'U' in varyList: |
---|
1734 | dMdv[varyList.index('U')] += dsdU*dervDict['sig'] |
---|
1735 | if 'V' in varyList: |
---|
1736 | dMdv[varyList.index('V')] += dsdV*dervDict['sig'] |
---|
1737 | if 'W' in varyList: |
---|
1738 | dMdv[varyList.index('W')] += dsdW*dervDict['sig'] |
---|
1739 | if 'X' in varyList: |
---|
1740 | dMdv[varyList.index('X')] += dgdX*dervDict['gam'] |
---|
1741 | if 'Y' in varyList: |
---|
1742 | dMdv[varyList.index('Y')] += dgdY*dervDict['gam'] |
---|
1743 | if 'Z' in varyList: |
---|
1744 | dMdv[varyList.index('Z')] += dgdZ*dervDict['gam'] |
---|
1745 | if 'SH/L' in varyList: |
---|
1746 | dMdv[varyList.index('SH/L')] += dervDict['shl'] #problem here |
---|
1747 | if 'I(L2)/I(L1)' in varyList: |
---|
1748 | dMdv[varyList.index('I(L2)/I(L1)')] += dervDict['L1/L2'] |
---|
1749 | iPeak += 1 |
---|
1750 | except KeyError: #no more peaks to process |
---|
1751 | break |
---|
1752 | elif 'E' in dataType: |
---|
1753 | iPeak = 0 |
---|
1754 | while True: |
---|
1755 | try: |
---|
1756 | pos = parmDict['pos'+str(iPeak)] |
---|
1757 | intens = parmDict['int'+str(iPeak)] |
---|
1758 | sigName = 'sig'+str(iPeak) |
---|
1759 | if sigName in varyList or not peakInstPrmMode: |
---|
1760 | sig = parmDict[sigName] |
---|
1761 | dsdA = dsdB = dsdC = 0 |
---|
1762 | else: |
---|
1763 | sig = G2mth.getEDsig(parmDict,pos) |
---|
1764 | dsdA,dsdB,dsdC = G2mth.getEDsigDeriv(parmDict,pos) |
---|
1765 | sig = max(sig,0.001) #avoid neg sigma |
---|
1766 | Wd,fmin,fmax = getWidthsED(pos,sig) |
---|
1767 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1768 | iFin = np.searchsorted(xdata,pos+fmin) |
---|
1769 | if not iBeg+iFin: #peak below low limit |
---|
1770 | iPeak += 1 |
---|
1771 | continue |
---|
1772 | elif not iBeg-iFin: #peak above high limit |
---|
1773 | break |
---|
1774 | dMdpk = np.zeros(shape=(4,len(xdata))) |
---|
1775 | dMdipk = getdPsVoigt(pos,sig*10.**4,0.001,xdata[iBeg:iFin]) |
---|
1776 | dMdpk[0][iBeg:iFin] += dMdipk[0] |
---|
1777 | for i in range(1,4): |
---|
1778 | dMdpk[i][iBeg:iFin] += intens*dMdipk[i] |
---|
1779 | dervDict = {'int':dMdpk[0],'pos':-dMdpk[1],'sig':dMdpk[2]*10**4} |
---|
1780 | for parmName in ['pos','int','sig']: |
---|
1781 | try: |
---|
1782 | idx = varyList.index(parmName+str(iPeak)) |
---|
1783 | dMdv[idx] = dervDict[parmName] |
---|
1784 | except ValueError: |
---|
1785 | pass |
---|
1786 | if 'A' in varyList: |
---|
1787 | dMdv[varyList.index('A')] += dsdA*dervDict['sig'] |
---|
1788 | if 'B' in varyList: |
---|
1789 | dMdv[varyList.index('B')] += dsdB*dervDict['sig'] |
---|
1790 | if 'C' in varyList: |
---|
1791 | dMdv[varyList.index('C')] += dsdC*dervDict['sig'] |
---|
1792 | iPeak += 1 |
---|
1793 | except KeyError: #no more peaks to process |
---|
1794 | break |
---|
1795 | |
---|
1796 | elif 'B' in dataType: |
---|
1797 | iPeak = 0 |
---|
1798 | while True: |
---|
1799 | try: |
---|
1800 | pos = parmDict['pos'+str(iPeak)] |
---|
1801 | tth = (pos-parmDict['Zero']) |
---|
1802 | intens = parmDict['int'+str(iPeak)] |
---|
1803 | alpName = 'alp'+str(iPeak) |
---|
1804 | if alpName in varyList or not peakInstPrmMode: |
---|
1805 | alp = parmDict[alpName] |
---|
1806 | dada0 = dada1 = 0.0 |
---|
1807 | else: |
---|
1808 | alp = G2mth.getPinkalpha(parmDict,tth) |
---|
1809 | dada0,dada1 = G2mth.getPinkalphaDeriv(tth) |
---|
1810 | alp = max(0.0001,alp) |
---|
1811 | betName = 'bet'+str(iPeak) |
---|
1812 | if betName in varyList or not peakInstPrmMode: |
---|
1813 | bet = parmDict[betName] |
---|
1814 | dbdb0 = dbdb1 = 0.0 |
---|
1815 | else: |
---|
1816 | bet = G2mth.getPinkbeta(parmDict,tth) |
---|
1817 | dbdb0,dbdb1 = G2mth.getPinkbetaDeriv(tth) |
---|
1818 | bet = max(0.0001,bet) |
---|
1819 | sigName = 'sig'+str(iPeak) |
---|
1820 | if sigName in varyList or not peakInstPrmMode: |
---|
1821 | sig = parmDict[sigName] |
---|
1822 | dsdU = dsdV = dsdW = 0 |
---|
1823 | else: |
---|
1824 | sig = G2mth.getCWsig(parmDict,tth) |
---|
1825 | dsdU,dsdV,dsdW = G2mth.getCWsigDeriv(tth) |
---|
1826 | sig = max(sig,0.001) #avoid neg sigma |
---|
1827 | gamName = 'gam'+str(iPeak) |
---|
1828 | if gamName in varyList or not peakInstPrmMode: |
---|
1829 | gam = parmDict[gamName] |
---|
1830 | dgdX = dgdY = dgdZ = 0 |
---|
1831 | else: |
---|
1832 | gam = G2mth.getCWgam(parmDict,tth) |
---|
1833 | dgdX,dgdY,dgdZ = G2mth.getCWgamDeriv(tth) |
---|
1834 | gam = max(gam,0.001) #avoid neg gamma |
---|
1835 | Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig/1.e4,gam/100.) |
---|
1836 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1837 | iFin = np.searchsorted(xdata,pos+fmin) |
---|
1838 | if not iBeg+iFin: #peak below low limit |
---|
1839 | iPeak += 1 |
---|
1840 | continue |
---|
1841 | elif not iBeg-iFin: #peak above high limit |
---|
1842 | break |
---|
1843 | dMdpk = np.zeros(shape=(7,len(xdata))) |
---|
1844 | dMdipk = getdEpsVoigt(pos,alp,bet,sig/1.e4,gam/100.,xdata[iBeg:iFin]) |
---|
1845 | for i in range(1,6): |
---|
1846 | dMdpk[i][iBeg:iFin] += intens*dMdipk[i] |
---|
1847 | dMdpk[0][iBeg:iFin] += dMdipk[0] |
---|
1848 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4]/1.e4,'gam':dMdpk[5]/100.} |
---|
1849 | for parmName in ['pos','int','alp','bet','sig','gam']: |
---|
1850 | try: |
---|
1851 | idx = varyList.index(parmName+str(iPeak)) |
---|
1852 | dMdv[idx] = dervDict[parmName] |
---|
1853 | except ValueError: |
---|
1854 | pass |
---|
1855 | if 'U' in varyList: |
---|
1856 | dMdv[varyList.index('U')] += dsdU*dervDict['sig'] |
---|
1857 | if 'V' in varyList: |
---|
1858 | dMdv[varyList.index('V')] += dsdV*dervDict['sig'] |
---|
1859 | if 'W' in varyList: |
---|
1860 | dMdv[varyList.index('W')] += dsdW*dervDict['sig'] |
---|
1861 | if 'X' in varyList: |
---|
1862 | dMdv[varyList.index('X')] += dgdX*dervDict['gam'] |
---|
1863 | if 'Y' in varyList: |
---|
1864 | dMdv[varyList.index('Y')] += dgdY*dervDict['gam'] |
---|
1865 | if 'Z' in varyList: |
---|
1866 | dMdv[varyList.index('Z')] += dgdZ*dervDict['gam'] |
---|
1867 | if 'alpha-0' in varyList: |
---|
1868 | dMdv[varyList.index('alpha-0')] += dada0*dervDict['alp'] |
---|
1869 | if 'alpha-1' in varyList: |
---|
1870 | dMdv[varyList.index('alpha-1')] += dada1*dervDict['alp'] |
---|
1871 | if 'beta-0' in varyList: |
---|
1872 | dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet'] |
---|
1873 | if 'beta-1' in varyList: |
---|
1874 | dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet'] |
---|
1875 | iPeak += 1 |
---|
1876 | except KeyError: #no more peaks to process |
---|
1877 | break |
---|
1878 | else: |
---|
1879 | Pdabc = parmDict['Pdabc'] |
---|
1880 | difC = parmDict['difC'] |
---|
1881 | iPeak = 0 |
---|
1882 | while True: |
---|
1883 | try: |
---|
1884 | pos = parmDict['pos'+str(iPeak)] |
---|
1885 | tof = pos-parmDict['Zero'] |
---|
1886 | dsp = tof/difC |
---|
1887 | intens = parmDict['int'+str(iPeak)] |
---|
1888 | alpName = 'alp'+str(iPeak) |
---|
1889 | if alpName in varyList or not peakInstPrmMode: |
---|
1890 | alp = parmDict[alpName] |
---|
1891 | else: |
---|
1892 | if len(Pdabc): |
---|
1893 | alp = np.interp(dsp,Pdabc[0],Pdabc[1]) |
---|
1894 | dada0 = 0 |
---|
1895 | else: |
---|
1896 | alp = G2mth.getTOFalpha(parmDict,dsp) |
---|
1897 | dada0 = G2mth.getTOFalphaDeriv(dsp) |
---|
1898 | betName = 'bet'+str(iPeak) |
---|
1899 | if betName in varyList or not peakInstPrmMode: |
---|
1900 | bet = parmDict[betName] |
---|
1901 | else: |
---|
1902 | if len(Pdabc): |
---|
1903 | bet = np.interp(dsp,Pdabc[0],Pdabc[2]) |
---|
1904 | dbdb0 = dbdb1 = dbdb2 = 0 |
---|
1905 | else: |
---|
1906 | bet = G2mth.getTOFbeta(parmDict,dsp) |
---|
1907 | dbdb0,dbdb1,dbdb2 = G2mth.getTOFbetaDeriv(dsp) |
---|
1908 | sigName = 'sig'+str(iPeak) |
---|
1909 | if sigName in varyList or not peakInstPrmMode: |
---|
1910 | sig = parmDict[sigName] |
---|
1911 | dsds0 = dsds1 = dsds2 = dsds3 = 0 |
---|
1912 | else: |
---|
1913 | sig = G2mth.getTOFsig(parmDict,dsp) |
---|
1914 | dsds0,dsds1,dsds2,dsds3 = G2mth.getTOFsigDeriv(dsp) |
---|
1915 | gamName = 'gam'+str(iPeak) |
---|
1916 | if gamName in varyList or not peakInstPrmMode: |
---|
1917 | gam = parmDict[gamName] |
---|
1918 | dsdX = dsdY = dsdZ = 0 |
---|
1919 | else: |
---|
1920 | gam = G2mth.getTOFgamma(parmDict,dsp) |
---|
1921 | dsdX,dsdY,dsdZ = G2mth.getTOFgammaDeriv(dsp) |
---|
1922 | gam = max(gam,0.001) #avoid neg gamma |
---|
1923 | Wd,fmin,fmax = getWidthsTOF(pos,alp,bet,sig,gam) |
---|
1924 | iBeg = np.searchsorted(xdata,pos-fmin) |
---|
1925 | lenX = len(xdata) |
---|
1926 | if not iBeg: |
---|
1927 | iFin = np.searchsorted(xdata,pos+fmax) |
---|
1928 | elif iBeg == lenX: |
---|
1929 | iFin = iBeg |
---|
1930 | else: |
---|
1931 | iFin = np.searchsorted(xdata,pos+fmax) |
---|
1932 | if not iBeg+iFin: #peak below low limit |
---|
1933 | iPeak += 1 |
---|
1934 | continue |
---|
1935 | elif not iBeg-iFin: #peak above high limit |
---|
1936 | break |
---|
1937 | dMdpk = np.zeros(shape=(7,len(xdata))) |
---|
1938 | dMdipk = getdEpsVoigt(pos,alp,bet,sig,gam,xdata[iBeg:iFin]) |
---|
1939 | for i in range(1,6): |
---|
1940 | dMdpk[i][iBeg:iFin] += intens*dMdipk[i] |
---|
1941 | dMdpk[0][iBeg:iFin] += dMdipk[0] |
---|
1942 | dervDict = {'int':dMdpk[0],'pos':dMdpk[1],'alp':dMdpk[2],'bet':dMdpk[3],'sig':dMdpk[4],'gam':dMdpk[5]} |
---|
1943 | for parmName in ['pos','int','alp','bet','sig','gam']: |
---|
1944 | try: |
---|
1945 | idx = varyList.index(parmName+str(iPeak)) |
---|
1946 | dMdv[idx] = dervDict[parmName] |
---|
1947 | except ValueError: |
---|
1948 | pass |
---|
1949 | if 'alpha' in varyList: |
---|
1950 | dMdv[varyList.index('alpha')] += dada0*dervDict['alp'] |
---|
1951 | if 'beta-0' in varyList: |
---|
1952 | dMdv[varyList.index('beta-0')] += dbdb0*dervDict['bet'] |
---|
1953 | if 'beta-1' in varyList: |
---|
1954 | dMdv[varyList.index('beta-1')] += dbdb1*dervDict['bet'] |
---|
1955 | if 'beta-q' in varyList: |
---|
1956 | dMdv[varyList.index('beta-q')] += dbdb2*dervDict['bet'] |
---|
1957 | if 'sig-0' in varyList: |
---|
1958 | dMdv[varyList.index('sig-0')] += dsds0*dervDict['sig'] |
---|
1959 | if 'sig-1' in varyList: |
---|
1960 | dMdv[varyList.index('sig-1')] += dsds1*dervDict['sig'] |
---|
1961 | if 'sig-2' in varyList: |
---|
1962 | dMdv[varyList.index('sig-2')] += dsds2*dervDict['sig'] |
---|
1963 | if 'sig-q' in varyList: |
---|
1964 | dMdv[varyList.index('sig-q')] += dsds3*dervDict['sig'] |
---|
1965 | if 'X' in varyList: |
---|
1966 | dMdv[varyList.index('X')] += dsdX*dervDict['gam'] |
---|
1967 | if 'Y' in varyList: |
---|
1968 | dMdv[varyList.index('Y')] += dsdY*dervDict['gam'] |
---|
1969 | if 'Z' in varyList: |
---|
1970 | dMdv[varyList.index('Z')] += dsdZ*dervDict['gam'] |
---|
1971 | iPeak += 1 |
---|
1972 | except KeyError: #no more peaks to process |
---|
1973 | break |
---|
1974 | if 'BF mult' in varyList: |
---|
1975 | dMdv[varyList.index('BF mult')] = fixback |
---|
1976 | |
---|
1977 | return dMdv |
---|
1978 | |
---|
1979 | def Dict2Values(parmdict, varylist): |
---|
1980 | '''Use before call to leastsq to setup list of values for the parameters |
---|
1981 | in parmdict, as selected by key in varylist''' |
---|
1982 | return [parmdict[key] for key in varylist] |
---|
1983 | |
---|
1984 | def Values2Dict(parmdict, varylist, values): |
---|
1985 | ''' Use after call to leastsq to update the parameter dictionary with |
---|
1986 | values corresponding to keys in varylist''' |
---|
1987 | parmdict.update(zip(varylist,values)) |
---|
1988 | |
---|
1989 | def SetBackgroundParms(Background): |
---|
1990 | 'Loads background parameters into dicts/lists to create varylist & parmdict' |
---|
1991 | if len(Background) == 1: # fix up old backgrounds |
---|
1992 | Background.append({'nDebye':0,'debyeTerms':[]}) |
---|
1993 | bakType,bakFlag = Background[0][:2] |
---|
1994 | backVals = Background[0][3:] |
---|
1995 | backNames = ['Back;'+str(i) for i in range(len(backVals))] |
---|
1996 | Debye = Background[1] #also has background peaks stuff |
---|
1997 | backDict = dict(zip(backNames,backVals)) |
---|
1998 | backVary = [] |
---|
1999 | if bakFlag: |
---|
2000 | backVary = backNames |
---|
2001 | |
---|
2002 | backDict['nDebye'] = Debye['nDebye'] |
---|
2003 | debyeDict = {} |
---|
2004 | debyeList = [] |
---|
2005 | for i in range(Debye['nDebye']): |
---|
2006 | debyeNames = ['DebyeA;'+str(i),'DebyeR;'+str(i),'DebyeU;'+str(i)] |
---|
2007 | debyeDict.update(dict(zip(debyeNames,Debye['debyeTerms'][i][::2]))) |
---|
2008 | debyeList += zip(debyeNames,Debye['debyeTerms'][i][1::2]) |
---|
2009 | debyeVary = [] |
---|
2010 | for item in debyeList: |
---|
2011 | if item[1]: |
---|
2012 | debyeVary.append(item[0]) |
---|
2013 | backDict.update(debyeDict) |
---|
2014 | backVary += debyeVary |
---|
2015 | |
---|
2016 | backDict['nPeaks'] = Debye['nPeaks'] |
---|
2017 | peaksDict = {} |
---|
2018 | peaksList = [] |
---|
2019 | for i in range(Debye['nPeaks']): |
---|
2020 | peaksNames = ['BkPkpos;'+str(i),'BkPkint;'+str(i),'BkPksig;'+str(i),'BkPkgam;'+str(i)] |
---|
2021 | peaksDict.update(dict(zip(peaksNames,Debye['peaksList'][i][::2]))) |
---|
2022 | peaksList += zip(peaksNames,Debye['peaksList'][i][1::2]) |
---|
2023 | peaksVary = [] |
---|
2024 | for item in peaksList: |
---|
2025 | if item[1]: |
---|
2026 | peaksVary.append(item[0]) |
---|
2027 | backDict.update(peaksDict) |
---|
2028 | backVary += peaksVary |
---|
2029 | if 'background PWDR' in Background[1]: |
---|
2030 | backDict['Back File'] = Background[1]['background PWDR'][0] |
---|
2031 | backDict['BF mult'] = Background[1]['background PWDR'][1] |
---|
2032 | if len(Background[1]['background PWDR']) > 2: |
---|
2033 | if Background[1]['background PWDR'][2]: |
---|
2034 | backVary += ['BF mult',] |
---|
2035 | return bakType,backDict,backVary |
---|
2036 | |
---|
2037 | def DoCalibInst(IndexPeaks,Inst): |
---|
2038 | |
---|
2039 | def SetInstParms(): |
---|
2040 | dataType = Inst['Type'][0] |
---|
2041 | insVary = [] |
---|
2042 | insNames = [] |
---|
2043 | insVals = [] |
---|
2044 | for parm in Inst: |
---|
2045 | insNames.append(parm) |
---|
2046 | insVals.append(Inst[parm][1]) |
---|
2047 | if parm in ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']: |
---|
2048 | if Inst[parm][2]: |
---|
2049 | insVary.append(parm) |
---|
2050 | instDict = dict(zip(insNames,insVals)) |
---|
2051 | return dataType,instDict,insVary |
---|
2052 | |
---|
2053 | def GetInstParms(parmDict,Inst,varyList): |
---|
2054 | for name in Inst: |
---|
2055 | Inst[name][1] = parmDict[name] |
---|
2056 | |
---|
2057 | def InstPrint(Inst,sigDict): |
---|
2058 | print ('Instrument Parameters:') |
---|
2059 | if 'C' in Inst['Type'][0] or 'B' in Inst['Type'][0]: |
---|
2060 | ptfmt = "%12.6f" |
---|
2061 | else: |
---|
2062 | ptfmt = "%12.3f" |
---|
2063 | ptlbls = 'names :' |
---|
2064 | ptstr = 'values:' |
---|
2065 | sigstr = 'esds :' |
---|
2066 | for parm in Inst: |
---|
2067 | if parm in ['Lam','difC','difA','difB','Zero','2-theta','XE','YE','ZE']: |
---|
2068 | ptlbls += "%s" % (parm.center(12)) |
---|
2069 | ptstr += ptfmt % (Inst[parm][1]) |
---|
2070 | if parm in sigDict: |
---|
2071 | sigstr += ptfmt % (sigDict[parm]) |
---|
2072 | else: |
---|
2073 | sigstr += 12*' ' |
---|
2074 | print (ptlbls) |
---|
2075 | print (ptstr) |
---|
2076 | print (sigstr) |
---|
2077 | |
---|
2078 | def errPeakPos(values,peakDsp,peakPos,peakWt,dataType,parmDict,varyList): |
---|
2079 | parmDict.update(zip(varyList,values)) |
---|
2080 | return np.sqrt(peakWt)*(G2lat.getPeakPos(dataType,parmDict,peakDsp)-peakPos) |
---|
2081 | |
---|
2082 | peakPos = [] |
---|
2083 | peakDsp = [] |
---|
2084 | peakWt = [] |
---|
2085 | for peak,sig in zip(IndexPeaks[0],IndexPeaks[1]): |
---|
2086 | if peak[2] and peak[3] and sig > 0.: |
---|
2087 | peakPos.append(peak[0]) |
---|
2088 | peakDsp.append(peak[-1]) #d-calc |
---|
2089 | # peakWt.append(peak[-1]**2/sig**2) #weight by d**2 |
---|
2090 | peakWt.append(1./(sig*peak[-1])) # |
---|
2091 | peakPos = np.array(peakPos) |
---|
2092 | peakDsp = np.array(peakDsp) |
---|
2093 | peakWt = np.array(peakWt) |
---|
2094 | dataType,insDict,insVary = SetInstParms() |
---|
2095 | parmDict = {} |
---|
2096 | parmDict.update(insDict) |
---|
2097 | varyList = insVary |
---|
2098 | if not len(varyList): |
---|
2099 | G2fil.G2Print ('**** ERROR - nothing to refine! ****') |
---|
2100 | return False |
---|
2101 | while True: |
---|
2102 | begin = time.time() |
---|
2103 | values = np.array(Dict2Values(parmDict, varyList)) |
---|
2104 | result = so.leastsq(errPeakPos,values,full_output=True,ftol=0.000001, |
---|
2105 | args=(peakDsp,peakPos,peakWt,dataType,parmDict,varyList)) |
---|
2106 | ncyc = int(result[2]['nfev']/2) |
---|
2107 | runtime = time.time()-begin |
---|
2108 | chisq = np.sum(result[2]['fvec']**2) |
---|
2109 | Values2Dict(parmDict, varyList, result[0]) |
---|
2110 | GOF = chisq/(len(peakPos)-len(varyList)) #reduced chi^2 |
---|
2111 | G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],len(peakPos),len(varyList))) |
---|
2112 | G2fil.G2Print ('calib time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)) |
---|
2113 | G2fil.G2Print ('chi**2 = %12.6g, reduced chi**2 = %6.2f'%(chisq,GOF)) |
---|
2114 | try: |
---|
2115 | sig = np.sqrt(np.diag(result[1])*GOF) |
---|
2116 | if np.any(np.isnan(sig)): |
---|
2117 | G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***') |
---|
2118 | break #refinement succeeded - finish up! |
---|
2119 | except ValueError: #result[1] is None on singular matrix |
---|
2120 | G2fil.G2Print ('**** Refinement failed - singular matrix ****') |
---|
2121 | |
---|
2122 | sigDict = dict(zip(varyList,sig)) |
---|
2123 | GetInstParms(parmDict,Inst,varyList) |
---|
2124 | InstPrint(Inst,sigDict) |
---|
2125 | return True |
---|
2126 | |
---|
2127 | def getHeaderInfo(dataType): |
---|
2128 | '''Provide parameter name, label name and formatting information for the |
---|
2129 | contents of the Peak Table and where used in DoPeakFit |
---|
2130 | ''' |
---|
2131 | names = ['pos','int'] |
---|
2132 | lnames = ['position','intensity'] |
---|
2133 | if 'LF' in dataType: |
---|
2134 | names = ['int','sig','gam','damp','asym','l','ttheta'] |
---|
2135 | lnames = ['intensity','sigma','gamma','satellite\ndamping', |
---|
2136 | 'satellite\nasym','00l', |
---|
2137 | #'2theta ' |
---|
2138 | '2\u03B8' |
---|
2139 | ] |
---|
2140 | fmt = ["%10.2f","%10.3f","%10.3f","%10.3f","%10.3f","%4d","%7.3f"] |
---|
2141 | elif 'C' in dataType: |
---|
2142 | names += ['sig','gam'] |
---|
2143 | lnames += ['sigma','gamma'] |
---|
2144 | fmt = ["%10.5f","%10.1f","%10.3f","%10.3f"] |
---|
2145 | elif 'T' in dataType: |
---|
2146 | names += ['alp','bet','sig','gam'] |
---|
2147 | lnames += ['alpha','beta','sigma','gamma'] |
---|
2148 | fmt = ["%10.2f","%10.4f","%8.3f","%8.5f","%10.3f","%10.3f"] |
---|
2149 | elif 'E' in dataType: |
---|
2150 | names += ['sig'] |
---|
2151 | lnames += ['sigma'] |
---|
2152 | fmt = ["%10.5f","%10.1f","%8.3f"] |
---|
2153 | else: # 'B' |
---|
2154 | names += ['alp','bet','sig','gam'] |
---|
2155 | lnames += ['alpha','beta','sigma','gamma'] |
---|
2156 | fmt = ["%10.5f","%10.1f","%8.2f","%8.4f","%10.3f","%10.3f"] |
---|
2157 | return names, fmt, lnames |
---|
2158 | |
---|
2159 | def DoPeakFit(FitPgm,Peaks,Background,Limits,Inst,Inst2,data,fixback=None,prevVaryList=[], |
---|
2160 | oneCycle=False,controls=None,wtFactor=1.0,dlg=None,noFit=False): |
---|
2161 | '''Called to perform a peak fit, refining the selected items in the peak |
---|
2162 | table as well as selected items in the background. |
---|
2163 | |
---|
2164 | :param str FitPgm: type of fit to perform. At present this is ignored. |
---|
2165 | :param list Peaks: a list of peaks. Each peak entry is a list with paired values: |
---|
2166 | The number of pairs depends on the data type (see :func:`getHeaderInfo`). |
---|
2167 | For CW data there are |
---|
2168 | four values each followed by a refine flag where the values are: position, intensity, |
---|
2169 | sigma (Gaussian width) and gamma (Lorentzian width). From the Histogram/"Peak List" |
---|
2170 | tree entry, dict item "peaks". For some types of fits, overall parameters are placed |
---|
2171 | in a dict entry. |
---|
2172 | :param list Background: describes the background. List with two items. |
---|
2173 | Item 0 specifies a background model and coefficients. Item 1 is a dict. |
---|
2174 | From the Histogram/Background tree entry. |
---|
2175 | :param list Limits: min and max x-value to use |
---|
2176 | :param dict Inst: Instrument parameters |
---|
2177 | :param dict Inst2: more Instrument parameters |
---|
2178 | :param numpy.array data: a 5xn array. data[0] is the x-values, |
---|
2179 | data[1] is the y-values, data[2] are weight values, data[3], [4] and [5] are |
---|
2180 | calc, background and difference intensities, respectively. |
---|
2181 | :param array fixback: fixed background array; same size as data[0-5] |
---|
2182 | :param list prevVaryList: Used in sequential refinements to override the |
---|
2183 | variable list. Defaults as an empty list. |
---|
2184 | :param bool oneCycle: True if only one cycle of fitting should be performed |
---|
2185 | :param dict controls: a dict specifying two values, Ftol = controls['min dM/M'] |
---|
2186 | and derivType = controls['deriv type']. If None default values are used. |
---|
2187 | :param float wtFactor: weight multiplier; = 1.0 by default |
---|
2188 | :param wx.Dialog dlg: A dialog box that is updated with progress from the fit. |
---|
2189 | Defaults to None, which means no updates are done. |
---|
2190 | :param bool noFit: When noFit is True, a refinement is not performed. Default |
---|
2191 | is False. |
---|
2192 | |
---|
2193 | ''' |
---|
2194 | def GetBackgroundParms(parmList,Background): |
---|
2195 | iBak = 0 |
---|
2196 | while True: |
---|
2197 | try: |
---|
2198 | bakName = 'Back;'+str(iBak) |
---|
2199 | Background[0][iBak+3] = parmList[bakName] |
---|
2200 | iBak += 1 |
---|
2201 | except KeyError: |
---|
2202 | break |
---|
2203 | iDb = 0 |
---|
2204 | while True: |
---|
2205 | names = ['DebyeA;','DebyeR;','DebyeU;'] |
---|
2206 | try: |
---|
2207 | for i,name in enumerate(names): |
---|
2208 | val = parmList[name+str(iDb)] |
---|
2209 | Background[1]['debyeTerms'][iDb][2*i] = val |
---|
2210 | iDb += 1 |
---|
2211 | except KeyError: |
---|
2212 | break |
---|
2213 | iDb = 0 |
---|
2214 | while True: |
---|
2215 | names = ['BkPkpos;','BkPkint;','BkPksig;','BkPkgam;'] |
---|
2216 | try: |
---|
2217 | for i,name in enumerate(names): |
---|
2218 | val = parmList[name+str(iDb)] |
---|
2219 | Background[1]['peaksList'][iDb][2*i] = val |
---|
2220 | iDb += 1 |
---|
2221 | except KeyError: |
---|
2222 | break |
---|
2223 | if 'BF mult' in parmList: |
---|
2224 | Background[1]['background PWDR'][1] = parmList['BF mult'] |
---|
2225 | |
---|
2226 | def BackgroundPrint(Background,sigDict): |
---|
2227 | print ('Background coefficients for '+Background[0][0]+' function') |
---|
2228 | ptfmt = "%12.5f" |
---|
2229 | ptstr = 'value: ' |
---|
2230 | sigstr = 'esd : ' |
---|
2231 | for i,back in enumerate(Background[0][3:]): |
---|
2232 | ptstr += ptfmt % (back) |
---|
2233 | if Background[0][1]: |
---|
2234 | prm = 'Back;'+str(i) |
---|
2235 | if prm in sigDict: |
---|
2236 | sigstr += ptfmt % (sigDict[prm]) |
---|
2237 | else: |
---|
2238 | sigstr += " "*12 |
---|
2239 | if len(ptstr) > 75: |
---|
2240 | print (ptstr) |
---|
2241 | if Background[0][1]: print (sigstr) |
---|
2242 | ptstr = 'value: ' |
---|
2243 | sigstr = 'esd : ' |
---|
2244 | if len(ptstr) > 8: |
---|
2245 | print (ptstr) |
---|
2246 | if Background[0][1]: print (sigstr) |
---|
2247 | |
---|
2248 | if Background[1]['nDebye']: |
---|
2249 | parms = ['DebyeA;','DebyeR;','DebyeU;'] |
---|
2250 | print ('Debye diffuse scattering coefficients') |
---|
2251 | ptfmt = "%12.5f" |
---|
2252 | print (' term DebyeA esd DebyeR esd DebyeU esd') |
---|
2253 | for term in range(Background[1]['nDebye']): |
---|
2254 | line = ' term %d'%(term) |
---|
2255 | for ip,name in enumerate(parms): |
---|
2256 | line += ptfmt%(Background[1]['debyeTerms'][term][2*ip]) |
---|
2257 | if name+str(term) in sigDict: |
---|
2258 | line += ptfmt%(sigDict[name+str(term)]) |
---|
2259 | else: |
---|
2260 | line += " "*12 |
---|
2261 | print (line) |
---|
2262 | if Background[1]['nPeaks']: |
---|
2263 | print ('Coefficients for Background Peaks') |
---|
2264 | ptfmt = "%15.3f" |
---|
2265 | for j,pl in enumerate(Background[1]['peaksList']): |
---|
2266 | names = 'peak %3d:'%(j+1) |
---|
2267 | ptstr = 'values :' |
---|
2268 | sigstr = 'esds :' |
---|
2269 | for i,lbl in enumerate(['BkPkpos','BkPkint','BkPksig','BkPkgam']): |
---|
2270 | val = pl[2*i] |
---|
2271 | prm = lbl+";"+str(j) |
---|
2272 | names += '%15s'%(prm) |
---|
2273 | ptstr += ptfmt%(val) |
---|
2274 | if prm in sigDict: |
---|
2275 | sigstr += ptfmt%(sigDict[prm]) |
---|
2276 | else: |
---|
2277 | sigstr += " "*15 |
---|
2278 | print (names) |
---|
2279 | print (ptstr) |
---|
2280 | print (sigstr) |
---|
2281 | if 'BF mult' in sigDict: |
---|
2282 | print('Background file mult: %.3f(%d)'%(Background[1]['background PWDR'][1],int(1000*sigDict['BF mult']))) |
---|
2283 | |
---|
2284 | def SetInstParms(Inst): |
---|
2285 | dataType = Inst['Type'][0] |
---|
2286 | insVary = [] |
---|
2287 | insNames = [] |
---|
2288 | insVals = [] |
---|
2289 | for parm in Inst: |
---|
2290 | insNames.append(parm) |
---|
2291 | insVals.append(Inst[parm][1]) |
---|
2292 | if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C', |
---|
2293 | 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1'] and Inst[parm][2]: |
---|
2294 | insVary.append(parm) |
---|
2295 | instDict = dict(zip(insNames,insVals)) |
---|
2296 | if 'SH/L' in instDict: |
---|
2297 | instDict['SH/L'] = max(instDict['SH/L'],0.002) |
---|
2298 | return dataType,instDict,insVary |
---|
2299 | |
---|
2300 | def GetPkInstParms(parmDict,Inst,varyList): |
---|
2301 | for name in Inst: |
---|
2302 | Inst[name][1] = parmDict[name] |
---|
2303 | iPeak = 0 |
---|
2304 | while True: |
---|
2305 | try: |
---|
2306 | sigName = 'sig'+str(iPeak) |
---|
2307 | pos = parmDict['pos'+str(iPeak)] |
---|
2308 | if sigName not in varyList and peakInstPrmMode: |
---|
2309 | if 'T' in Inst['Type'][0]: |
---|
2310 | dsp = G2lat.Pos2dsp(Inst,pos) |
---|
2311 | parmDict[sigName] = G2mth.getTOFsig(parmDict,dsp) |
---|
2312 | if 'E' in Inst['Type'][0]: |
---|
2313 | parmDict[sigName] = G2mth.getEDsig(parmDict,pos) |
---|
2314 | else: |
---|
2315 | parmDict[sigName] = G2mth.getCWsig(parmDict,pos) |
---|
2316 | gamName = 'gam'+str(iPeak) |
---|
2317 | if gamName not in varyList and peakInstPrmMode: |
---|
2318 | if 'T' in Inst['Type'][0]: |
---|
2319 | dsp = G2lat.Pos2dsp(Inst,pos) |
---|
2320 | parmDict[gamName] = G2mth.getTOFgamma(parmDict,dsp) |
---|
2321 | else: |
---|
2322 | parmDict[gamName] = G2mth.getCWgam(parmDict,pos) |
---|
2323 | iPeak += 1 |
---|
2324 | except KeyError: |
---|
2325 | break |
---|
2326 | |
---|
2327 | def InstPrint(Inst,sigDict): |
---|
2328 | print ('Instrument Parameters:') |
---|
2329 | ptfmt = "%12.6f" |
---|
2330 | ptlbls = 'names :' |
---|
2331 | ptstr = 'values:' |
---|
2332 | sigstr = 'esds :' |
---|
2333 | for parm in Inst: |
---|
2334 | if parm in ['U','V','W','X','Y','Z','SH/L','I(L2)/I(L1)','alpha','A','B','C', |
---|
2335 | 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q','alpha-0','alpha-1']: |
---|
2336 | ptlbls += "%s" % (parm.center(12)) |
---|
2337 | ptstr += ptfmt % (Inst[parm][1]) |
---|
2338 | if parm in sigDict: |
---|
2339 | sigstr += ptfmt % (sigDict[parm]) |
---|
2340 | else: |
---|
2341 | sigstr += 12*' ' |
---|
2342 | print (ptlbls) |
---|
2343 | print (ptstr) |
---|
2344 | print (sigstr) |
---|
2345 | |
---|
2346 | def SetPeaksParms(dataType,Peaks): |
---|
2347 | '''Set the contents of peakDict from list Peaks |
---|
2348 | ''' |
---|
2349 | peakDict = {} |
---|
2350 | peakVary = [] |
---|
2351 | names,_,_ = getHeaderInfo(dataType) |
---|
2352 | if 'LF' in dataType: |
---|
2353 | off = 2 |
---|
2354 | names = names[:-1] # drop 2nd 2theta value |
---|
2355 | else: |
---|
2356 | off = 0 |
---|
2357 | for i,peak in enumerate(Peaks): |
---|
2358 | if type(peak) is dict: |
---|
2359 | peakDict.update(peak) |
---|
2360 | continue |
---|
2361 | for j,name in enumerate(names): |
---|
2362 | parName = name+str(i) |
---|
2363 | peakDict[parName] = peak[off+2*j] |
---|
2364 | if peak[off+2*j+1]: |
---|
2365 | peakVary.append(parName) |
---|
2366 | return peakDict,peakVary |
---|
2367 | |
---|
2368 | def GetPeaksParms(Inst,parmDict,Peaks,varyList): |
---|
2369 | '''Put values into the Peaks list from the refinement results from inside |
---|
2370 | the parmDict array |
---|
2371 | ''' |
---|
2372 | names,_,_ = getHeaderInfo(Inst['Type'][0]) |
---|
2373 | off = 0 |
---|
2374 | if 'LF' in Inst['Type'][0]: |
---|
2375 | off = 2 |
---|
2376 | if 'clat' in varyList: |
---|
2377 | Peaks[-1]['clat'] = parmDict['clat'] |
---|
2378 | names = names[:-1] # drop 2nd 2theta value |
---|
2379 | for i,peak in enumerate(Peaks): |
---|
2380 | if type(peak) is dict: |
---|
2381 | continue |
---|
2382 | for j in range(len(names)): |
---|
2383 | parName = names[j]+str(i) |
---|
2384 | if parName in varyList or not peakInstPrmMode: |
---|
2385 | peak[2*j+off] = parmDict[parName] |
---|
2386 | if 'pos'+str(i) not in parmDict: continue |
---|
2387 | pos = parmDict['pos'+str(i)] |
---|
2388 | if 'LF' in Inst['Type'][0]: |
---|
2389 | peak[0] = pos |
---|
2390 | peak[-1] = pos |
---|
2391 | if 'difC' in Inst: |
---|
2392 | dsp = pos/Inst['difC'][1] |
---|
2393 | for j in range(len(names)): |
---|
2394 | parName = names[j]+str(i) |
---|
2395 | if peak[2*j+off + 1] or not peakInstPrmMode: continue |
---|
2396 | if 'alp' in parName: |
---|
2397 | if 'T' in Inst['Type'][0]: |
---|
2398 | peak[2*j+off] = G2mth.getTOFalpha(parmDict,dsp) |
---|
2399 | else: #'B' |
---|
2400 | peak[2*j+off] = G2mth.getPinkalpha(parmDict,pos) |
---|
2401 | elif 'bet' in parName: |
---|
2402 | if 'T' in Inst['Type'][0]: |
---|
2403 | peak[2*j+off] = G2mth.getTOFbeta(parmDict,dsp) |
---|
2404 | else: #'B' |
---|
2405 | peak[2*j+off] = G2mth.getPinkbeta(parmDict,pos) |
---|
2406 | elif 'sig' in parName: |
---|
2407 | if 'T' in Inst['Type'][0]: |
---|
2408 | peak[2*j+off] = G2mth.getTOFsig(parmDict,dsp) |
---|
2409 | elif 'E' in Inst['Type'][0]: |
---|
2410 | peak[2*j+off] = G2mth.getEDsig(parmDict,pos) |
---|
2411 | else: #'C' & 'B' |
---|
2412 | peak[2*j+off] = G2mth.getCWsig(parmDict,pos) |
---|
2413 | elif 'gam' in parName: |
---|
2414 | if 'T' in Inst['Type'][0]: |
---|
2415 | peak[2*j+off] = G2mth.getTOFgamma(parmDict,dsp) |
---|
2416 | else: #'C' & 'B' |
---|
2417 | peak[2*j+off] = G2mth.getCWgam(parmDict,pos) |
---|
2418 | |
---|
2419 | def PeaksPrint(dataType,parmDict,sigDict,varyList,ptsperFW): |
---|
2420 | if 'clat' in varyList: |
---|
2421 | print('c = {:.6f} esd {:.6f}'.format(parmDict['clat'],sigDict['clat'])) |
---|
2422 | print ('Peak coefficients:') |
---|
2423 | names,fmt,_ = getHeaderInfo(dataType) |
---|
2424 | head = 13*' ' |
---|
2425 | for name in names: |
---|
2426 | if name == 'l': |
---|
2427 | head += name.center(3)+' ' |
---|
2428 | elif name == 'ttheta': |
---|
2429 | continue |
---|
2430 | elif name in ['alp','bet']: |
---|
2431 | head += name.center(8)+'esd'.center(8) |
---|
2432 | else: |
---|
2433 | head += name.center(10)+'esd'.center(10) |
---|
2434 | head += 'bins'.center(8) |
---|
2435 | print (head) |
---|
2436 | ptfmt = dict(zip(names,fmt)) |
---|
2437 | for i,peak in enumerate(Peaks): |
---|
2438 | if type(peak) is dict: |
---|
2439 | continue |
---|
2440 | ptstr = ':' |
---|
2441 | for j in range(len(names)): |
---|
2442 | name = names[j] |
---|
2443 | parName = name+str(i) |
---|
2444 | if parName not in parmDict: continue |
---|
2445 | ptstr += ptfmt[name] % (parmDict[parName]) |
---|
2446 | if parName in varyList: |
---|
2447 | ptstr += ptfmt[name] % (sigDict[parName]) |
---|
2448 | else: |
---|
2449 | if name in ['alp','bet']: |
---|
2450 | ptstr += 8*' ' |
---|
2451 | else: |
---|
2452 | ptstr += 10*' ' |
---|
2453 | ptstr += '%9.2f'%(ptsperFW[i]) |
---|
2454 | print ('%s'%(('Peak'+str(i+1)).center(8)),ptstr) |
---|
2455 | |
---|
2456 | def devPeakProfile(values,xdata,ydata,fixback, weights,dataType,parmdict,varylist,bakType,dlg): |
---|
2457 | '''Computes a matrix where each row is the derivative of the calc-obs |
---|
2458 | values (see :func:`errPeakProfile`) with respect to each parameter |
---|
2459 | in backVary,insVary,peakVary. Used for peak fitting. |
---|
2460 | ''' |
---|
2461 | parmdict.update(zip(varylist,values)) |
---|
2462 | return np.sqrt(weights)*getPeakProfileDerv(dataType,parmdict,xdata,fixback,varylist,bakType) |
---|
2463 | |
---|
2464 | def errPeakProfile(values,xdata,ydata,fixback,weights,dataType,parmdict,varylist,bakType,dlg): |
---|
2465 | '''Computes a vector with the weighted calc-obs values differences |
---|
2466 | for peak fitting |
---|
2467 | ''' |
---|
2468 | parmdict.update(zip(varylist,values)) |
---|
2469 | M = np.sqrt(weights)*(getPeakProfile(dataType,parmdict,xdata,fixback,varylist,bakType)-ydata) |
---|
2470 | Rwp = min(100.,np.sqrt(np.sum(M**2)/np.sum(weights*ydata**2))*100.) |
---|
2471 | if dlg: |
---|
2472 | dlg.Raise() |
---|
2473 | GoOn = dlg.Update(int(Rwp),newmsg='%s%8.3f%s'%('Peak fit Rwp =',Rwp,'%'))[0] |
---|
2474 | if not GoOn: |
---|
2475 | return -M #abort!! |
---|
2476 | return M |
---|
2477 | |
---|
2478 | #---- beginning of DoPeakFit --------------------------------------------- |
---|
2479 | if controls: |
---|
2480 | Ftol = controls['min dM/M'] |
---|
2481 | else: |
---|
2482 | Ftol = 0.0001 |
---|
2483 | if oneCycle: |
---|
2484 | Ftol = 1.0 |
---|
2485 | x,y,w,yc,yb,yd = data #these are numpy arrays - remove masks! |
---|
2486 | if fixback is None: |
---|
2487 | fixback = np.zeros_like(y) |
---|
2488 | yc.fill(0.) #set calcd ones to zero |
---|
2489 | yb.fill(0.) |
---|
2490 | yd.fill(0.) |
---|
2491 | xBeg = np.searchsorted(x,Limits[0]) |
---|
2492 | xFin = np.searchsorted(x,Limits[1])+1 |
---|
2493 | # find out what is varied |
---|
2494 | bakType,bakDict,bakVary = SetBackgroundParms(Background) |
---|
2495 | dataType,insDict,insVary = SetInstParms(Inst) |
---|
2496 | peakDict,peakVary = SetPeaksParms(Inst['Type'][0],Peaks) |
---|
2497 | parmDict = {} |
---|
2498 | parmDict.update(bakDict) |
---|
2499 | parmDict.update(insDict) |
---|
2500 | parmDict.update(peakDict) |
---|
2501 | parmDict['Pdabc'] = [] #dummy Pdabc |
---|
2502 | parmDict.update(Inst2) #put in real one if there |
---|
2503 | if prevVaryList: |
---|
2504 | varyList = prevVaryList[:] |
---|
2505 | else: |
---|
2506 | varyList = bakVary+insVary+peakVary |
---|
2507 | if 'LF' in Inst['Type'][0] and Peaks: |
---|
2508 | if Peaks[-1].get('clat-ref'): varyList += ['clat'] |
---|
2509 | fullvaryList = varyList[:] |
---|
2510 | if not peakInstPrmMode: |
---|
2511 | for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1','A','B','C', |
---|
2512 | 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',): |
---|
2513 | if v in varyList: |
---|
2514 | raise Exception('Instrumental profile terms cannot be varied '+ |
---|
2515 | 'after setPeakInstPrmMode(False) is used') |
---|
2516 | if 'LF' in Inst['Type'][0]: |
---|
2517 | warn = [] |
---|
2518 | for v in ('U','V','W','X','Y','Z','alpha','alpha-0','alpha-1', |
---|
2519 | 'beta-0','beta-1','beta-q','sig-0','sig-1','sig-2','sig-q',): |
---|
2520 | if v in varyList: |
---|
2521 | warn.append(v) |
---|
2522 | del varyList[varyList.index(v)] |
---|
2523 | if warn: |
---|
2524 | print('Instrumental profile terms cannot be varied '+ |
---|
2525 | 'in Laue Fringe fits:',warn) |
---|
2526 | |
---|
2527 | while not noFit: |
---|
2528 | begin = time.time() |
---|
2529 | values = np.array(Dict2Values(parmDict, varyList)) |
---|
2530 | Rvals = {} |
---|
2531 | badVary = [] |
---|
2532 | try: |
---|
2533 | result = so.leastsq(errPeakProfile,values,Dfun=devPeakProfile,full_output=True,ftol=Ftol,col_deriv=True, |
---|
2534 | args=(x[xBeg:xFin],y[xBeg:xFin],fixback[xBeg:xFin],wtFactor*w[xBeg:xFin],dataType,parmDict,varyList,bakType,dlg)) |
---|
2535 | except Exception as msg: |
---|
2536 | if GSASIIpath.GetConfigValue('debug'): |
---|
2537 | print('peak fit failure\n',msg) |
---|
2538 | import traceback |
---|
2539 | print (traceback.format_exc()) |
---|
2540 | else: |
---|
2541 | print('peak fit failure') |
---|
2542 | return |
---|
2543 | ncyc = int(result[2]['nfev']/2) |
---|
2544 | runtime = time.time()-begin |
---|
2545 | chisq = np.sum(result[2]['fvec']**2) |
---|
2546 | Values2Dict(parmDict, varyList, result[0]) |
---|
2547 | Rvals['Rwp'] = np.sqrt(chisq/np.sum(wtFactor*w[xBeg:xFin]*y[xBeg:xFin]**2))*100. #to % |
---|
2548 | Rvals['GOF'] = chisq/(xFin-xBeg-len(varyList)) #reduced chi^2 |
---|
2549 | G2fil.G2Print ('Number of function calls: %d Number of observations: %d Number of parameters: %d'%(result[2]['nfev'],xFin-xBeg,len(varyList))) |
---|
2550 | if ncyc: |
---|
2551 | G2fil.G2Print ('fitpeak time = %8.3fs, %8.3fs/cycle'%(runtime,runtime/ncyc)) |
---|
2552 | G2fil.G2Print ('Rwp = %7.2f%%, chi**2 = %12.6g, reduced chi**2 = %6.2f'%(Rvals['Rwp'],chisq,Rvals['GOF'])) |
---|
2553 | sig = [0]*len(varyList) |
---|
2554 | if len(varyList) == 0: break # if nothing was refined |
---|
2555 | try: |
---|
2556 | sig = np.sqrt(np.diag(result[1])*Rvals['GOF']) |
---|
2557 | if np.any(np.isnan(sig)): |
---|
2558 | G2fil.G2Print ('*** Least squares aborted - some invalid esds possible ***') |
---|
2559 | break #refinement succeeded - finish up! |
---|
2560 | except ValueError: #result[1] is None on singular matrix |
---|
2561 | G2fil.G2Print ('**** Refinement failed - singular matrix ****') |
---|
2562 | Ipvt = result[2]['ipvt'] |
---|
2563 | for i,ipvt in enumerate(Ipvt): |
---|
2564 | if not np.sum(result[2]['fjac'],axis=1)[i]: |
---|
2565 | G2fil.G2Print ('Removing parameter: '+varyList[ipvt-1]) |
---|
2566 | badVary.append(varyList[ipvt-1]) |
---|
2567 | del(varyList[ipvt-1]) |
---|
2568 | break |
---|
2569 | else: # nothing removed |
---|
2570 | break |
---|
2571 | if dlg: dlg.Destroy() |
---|
2572 | yb[xBeg:xFin] = getBackground('',parmDict,bakType,dataType,x[xBeg:xFin],fixback[xBeg:xFin])[0] |
---|
2573 | yc[xBeg:xFin] = getPeakProfile(dataType,parmDict,x[xBeg:xFin],fixback[xBeg:xFin],varyList,bakType) |
---|
2574 | yd[xBeg:xFin] = y[xBeg:xFin]-yc[xBeg:xFin] |
---|
2575 | if noFit: |
---|
2576 | GetPeaksParms(Inst,parmDict,Peaks,varyList) |
---|
2577 | return |
---|
2578 | sigDict = dict(zip(varyList,sig)) |
---|
2579 | GetBackgroundParms(parmDict,Background) |
---|
2580 | if bakVary: BackgroundPrint(Background,sigDict) |
---|
2581 | GetPkInstParms(parmDict,Inst,varyList) |
---|
2582 | if insVary: InstPrint(Inst,sigDict) |
---|
2583 | GetPeaksParms(Inst,parmDict,Peaks,varyList) |
---|
2584 | binsperFWHM = [] |
---|
2585 | for peak in Peaks: |
---|
2586 | if type(peak) is dict: |
---|
2587 | continue |
---|
2588 | FWHM = getFWHM(peak[0],Inst) |
---|
2589 | try: |
---|
2590 | xpk = x.searchsorted(peak[0]) |
---|
2591 | cw = x[xpk]-x[xpk-1] |
---|
2592 | binsperFWHM.append(FWHM/cw) |
---|
2593 | except IndexError: |
---|
2594 | binsperFWHM.append(0.) |
---|
2595 | if peakVary: PeaksPrint(dataType,parmDict,sigDict,varyList,binsperFWHM) |
---|
2596 | if len(binsperFWHM): |
---|
2597 | if min(binsperFWHM) < 1.: |
---|
2598 | G2fil.G2Print ('*** Warning: calculated peak widths are too narrow to refine profile coefficients ***') |
---|
2599 | if 'T' in Inst['Type'][0]: |
---|
2600 | G2fil.G2Print (' Manually increase sig-0, 1, or 2 in Instrument Parameters') |
---|
2601 | else: |
---|
2602 | G2fil.G2Print (' Manually increase W in Instrument Parameters') |
---|
2603 | elif min(binsperFWHM) < 4.: |
---|
2604 | G2fil.G2Print ('*** Warning: data binning yields too few data points across peak FWHM for reliable Rietveld refinement ***') |
---|
2605 | G2fil.G2Print ('*** recommended is 6-10; you have %.2f ***'%(min(binsperFWHM))) |
---|
2606 | return sigDict,result,sig,Rvals,varyList,parmDict,fullvaryList,badVary |
---|
2607 | |
---|
2608 | def calcIncident(Iparm,xdata): |
---|
2609 | 'needs a doc string' |
---|
2610 | |
---|
2611 | def IfunAdv(Iparm,xdata): |
---|
2612 | Itype = Iparm['Itype'] |
---|
2613 | Icoef = Iparm['Icoeff'] |
---|
2614 | DYI = np.ones((12,xdata.shape[0])) |
---|
2615 | YI = np.ones_like(xdata)*Icoef[0] |
---|
2616 | |
---|
2617 | x = xdata/1000. #expressions are in ms |
---|
2618 | if Itype == 'Exponential': |
---|
2619 | for i in [1,3,5,7,9]: |
---|
2620 | Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2)) |
---|
2621 | YI += Icoef[i]*Eterm |
---|
2622 | DYI[i] *= Eterm |
---|
2623 | DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2) |
---|
2624 | elif 'Maxwell'in Itype: |
---|
2625 | Eterm = np.exp(-Icoef[2]/x**2) |
---|
2626 | DYI[1] = Eterm/x**5 |
---|
2627 | DYI[2] = -Icoef[1]*DYI[1]/x**2 |
---|
2628 | YI += (Icoef[1]*Eterm/x**5) |
---|
2629 | if 'Exponential' in Itype: |
---|
2630 | for i in range(3,11,2): |
---|
2631 | Eterm = np.exp(-Icoef[i+1]*x**((i+1)/2)) |
---|
2632 | YI += Icoef[i]*Eterm |
---|
2633 | DYI[i] *= Eterm |
---|
2634 | DYI[i+1] *= -Icoef[i]*Eterm*x**((i+1)/2) |
---|
2635 | else: #Chebyschev |
---|
2636 | T = (2./x)-1. |
---|
2637 | Ccof = np.ones((12,xdata.shape[0])) |
---|
2638 | Ccof[1] = T |
---|
2639 | for i in range(2,12): |
---|
2640 | Ccof[i] = 2*T*Ccof[i-1]-Ccof[i-2] |
---|
2641 | for i in range(1,10): |
---|
2642 | YI += Ccof[i]*Icoef[i+2] |
---|
2643 | DYI[i+2] =Ccof[i] |
---|
2644 | return YI,DYI |
---|
2645 | |
---|
2646 | Iesd = np.array(Iparm['Iesd']) |
---|
2647 | Icovar = Iparm['Icovar'] |
---|
2648 | YI,DYI = IfunAdv(Iparm,xdata) |
---|
2649 | YI = np.where(YI>0,YI,1.) |
---|
2650 | WYI = np.zeros_like(xdata) |
---|
2651 | vcov = np.zeros((12,12)) |
---|
2652 | k = 0 |
---|
2653 | for i in range(12): |
---|
2654 | for j in range(i,12): |
---|
2655 | vcov[i][j] = Icovar[k]*Iesd[i]*Iesd[j] |
---|
2656 | vcov[j][i] = Icovar[k]*Iesd[i]*Iesd[j] |
---|
2657 | k += 1 |
---|
2658 | M = np.inner(vcov,DYI.T) |
---|
2659 | WYI = np.sum(M*DYI,axis=0) |
---|
2660 | WYI = np.where(WYI>0.,WYI,0.) |
---|
2661 | return YI,WYI |
---|
2662 | |
---|
2663 | #### RMCutilities ################################################################################ |
---|
2664 | def MakeInst(PWDdata,Name,Size,Mustrain,useSamBrd): |
---|
2665 | inst = PWDdata['Instrument Parameters'][0] |
---|
2666 | sample = PWDdata['Sample Parameters'] |
---|
2667 | Xsb = 0. |
---|
2668 | Ysb = 0. |
---|
2669 | if 'T' in inst['Type'][1]: |
---|
2670 | difC = inst['difC'][1] |
---|
2671 | if useSamBrd[0]: |
---|
2672 | if 'ellipsoidal' not in Size[0]: #take the isotropic term only |
---|
2673 | Xsb = 1.e-4*difC/Size[1][0] |
---|
2674 | if useSamBrd[1]: |
---|
2675 | if 'generalized' not in Mustrain[0]: #take the isotropic term only |
---|
2676 | Ysb = 1.e-6*difC*Mustrain[1][0] |
---|
2677 | prms = ['Bank', |
---|
2678 | 'difC','difA','Zero','2-theta','difB', |
---|
2679 | 'alpha','beta-0','beta-1','beta-q', |
---|
2680 | 'sig-0','sig-1','sig-2','sig-q', |
---|
2681 | 'Z','X','Y'] |
---|
2682 | fname = Name+'.inst' |
---|
2683 | fl = open(fname,'w') |
---|
2684 | fl.write('1\n') |
---|
2685 | fl.write('%d\n'%int(inst[prms[0]][1])) |
---|
2686 | fl.write('%19.11f%19.11f%19.11f%19.11f%19.11f\n'%(inst[prms[1]][1],inst[prms[2]][1],inst[prms[3]][1],inst[prms[4]][1],inst[prms[5]][1],)) |
---|
2687 | fl.write('%12.6e%14.6e%14.6e%14.6e\n'%(inst[prms[6]][1],inst[prms[7]][1],inst[prms[8]][1],inst[prms[9]][1])) |
---|
2688 | fl.write('%12.6e%14.6e%14.6e%14.6e\n'%(inst[prms[10]][1],inst[prms[11]][1],inst[prms[12]][1],inst[prms[13]][1])) |
---|
2689 | fl.write('%12.6e%14.6e%14.6e%14.6e%14.6e\n'%(inst[prms[14]][1],inst[prms[15]][1]+Ysb,inst[prms[16]][1]+Xsb,0.0,0.0)) |
---|
2690 | fl.write('%12.6e\n\n\n'%(sample['Absorption'][0])) |
---|
2691 | fl.close() |
---|
2692 | else: |
---|
2693 | if useSamBrd[0]: |
---|
2694 | wave = G2mth.getWave(inst) |
---|
2695 | if 'ellipsoidal' not in Size[0]: #take the isotropic term only |
---|
2696 | Xsb = 1.8*wave/(np.pi*Size[1][0]) |
---|
2697 | if useSamBrd[1]: |
---|
2698 | if 'generalized' not in Mustrain[0]: #take the isotropic term only |
---|
2699 | Ysb = 0.0180*Mustrain[1][0]/np.pi |
---|
2700 | prms = ['Bank', |
---|
2701 | 'Lam','Zero','Polariz.', |
---|
2702 | 'U','V','W', |
---|
2703 | 'X','Y'] |
---|
2704 | fname = Name+'.inst' |
---|
2705 | fl = open(fname,'w') |
---|
2706 | fl.write('1\n') |
---|
2707 | fl.write('%d\n'%int(inst[prms[0]][1])) |
---|
2708 | fl.write('%10.5f%10.5f%10.4f%10d\n'%(inst[prms[1]][1],100.*inst[prms[2]][1],inst[prms[3]][1],0)) |
---|
2709 | fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[4]][1],inst[prms[5]][1],inst[prms[6]][1])) |
---|
2710 | fl.write('%10.3f%10.3f%10.3f\n'%(inst[prms[7]][1]+Xsb,inst[prms[8]][1]+Ysb,0.0)) |
---|
2711 | fl.write('%10.3f%10.3f%10.3f\n'%(0.0,0.0,0.0)) |
---|
2712 | fl.write('%12.6e\n\n\n'%(sample['Absorption'][0])) |
---|
2713 | fl.close() |
---|
2714 | return fname |
---|
2715 | |
---|
2716 | def MakeBack(PWDdata,Name): |
---|
2717 | Back = PWDdata['Background'][0] |
---|
2718 | inst = PWDdata['Instrument Parameters'][0] |
---|
2719 | if 'chebyschev-1' != Back[0]: |
---|
2720 | return None |
---|
2721 | Nback = Back[2] |
---|
2722 | BackVals = Back[3:] |
---|
2723 | fname = Name+'.back' |
---|
2724 | fl = open(fname,'w') |
---|
2725 | fl.write('%10d\n'%Nback) |
---|
2726 | for val in BackVals: |
---|
2727 | if 'T' in inst['Type'][1]: |
---|
2728 | fl.write('%12.6g\n'%(float(val))) |
---|
2729 | else: |
---|
2730 | fl.write('%12.6g\n'%val) |
---|
2731 | fl.close() |
---|
2732 | return fname |
---|
2733 | |
---|
2734 | def findDup(Atoms): |
---|
2735 | Dup = [] |
---|
2736 | Fracs = [] |
---|
2737 | for iat1,at1 in enumerate(Atoms): |
---|
2738 | if any([at1[0] in dup for dup in Dup]): |
---|
2739 | continue |
---|
2740 | else: |
---|
2741 | Dup.append([at1[0],]) |
---|
2742 | Fracs.append([at1[6],]) |
---|
2743 | for iat2,at2 in enumerate(Atoms[(iat1+1):]): |
---|
2744 | if np.sum((np.array(at1[3:6])-np.array(at2[3:6]))**2) < 0.00001: |
---|
2745 | Dup[-1] += [at2[0],] |
---|
2746 | Fracs[-1]+= [at2[6],] |
---|
2747 | return Dup,Fracs |
---|
2748 | |
---|
2749 | def MakeRMC6f(PWDdata,Name,Phase,RMCPdict): |
---|
2750 | |
---|
2751 | Meta = RMCPdict['metadata'] |
---|
2752 | Atseq = RMCPdict['atSeq'] |
---|
2753 | Supercell = RMCPdict['SuperCell'] |
---|
2754 | generalData = Phase['General'] |
---|
2755 | Dups,Fracs = findDup(Phase['Atoms']) |
---|
2756 | Sfracs = [np.cumsum(fracs) for fracs in Fracs] |
---|
2757 | ifSfracs = any([np.any(sfracs-1.) for sfracs in Sfracs]) |
---|
2758 | Sample = PWDdata['Sample Parameters'] |
---|
2759 | Meta['temperature'] = Sample['Temperature'] |
---|
2760 | Meta['pressure'] = Sample['Pressure'] |
---|
2761 | Cell = generalData['Cell'][1:7] |
---|
2762 | Trans = np.eye(3)*np.array(Supercell) |
---|
2763 | newPhase = copy.deepcopy(Phase) |
---|
2764 | newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1] |
---|
2765 | newPhase['General']['Cell'][1:] = G2lat.TransformCell(Cell,Trans) |
---|
2766 | GB = G2lat.cell2Gmat( newPhase['General']['Cell'][1:7])[0] |
---|
2767 | RMCPdict['Rmax'] = np.min(np.sqrt(np.array([1./G2lat.calc_rDsq2(H,GB) for H in [[1,0,0],[0,1,0],[0,0,1]]])))/2. |
---|
2768 | newPhase,Atcodes = G2lat.TransformPhase(Phase,newPhase,Trans,np.zeros(3),np.zeros(3),ifMag=False,Force=True) |
---|
2769 | Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell |
---|
2770 | Natm = np.count_nonzero(Natm-1) |
---|
2771 | Atoms = newPhase['Atoms'] |
---|
2772 | reset = False |
---|
2773 | |
---|
2774 | if ifSfracs: |
---|
2775 | Natm = np.core.defchararray.count(np.array(Atcodes),'+') #no. atoms in original unit cell |
---|
2776 | Natm = np.count_nonzero(Natm-1) |
---|
2777 | Satoms = [] |
---|
2778 | for i in range(len(Atoms)//Natm): |
---|
2779 | ind = i*Natm |
---|
2780 | Satoms.append(G2mth.sortArray(G2mth.sortArray(G2mth.sortArray(Atoms[ind:ind+Natm],5),4),3)) |
---|
2781 | Natoms = [] |
---|
2782 | for satoms in Satoms: |
---|
2783 | for idup,dup in enumerate(Dups): |
---|
2784 | ldup = len(dup) |
---|
2785 | natm = len(satoms) |
---|
2786 | i = 0 |
---|
2787 | while i < natm: |
---|
2788 | if satoms[i][0] in dup: |
---|
2789 | atoms = satoms[i:i+ldup] |
---|
2790 | try: |
---|
2791 | atom = atoms[np.searchsorted(Sfracs[idup],rand.random())] |
---|
2792 | Natoms.append(atom) |
---|
2793 | except IndexError: #what about vacancies? |
---|
2794 | if 'Va' not in Atseq: |
---|
2795 | reset = True |
---|
2796 | Atseq.append('Va') |
---|
2797 | RMCPdict['aTypes']['Va'] = 0.0 |
---|
2798 | atom = atoms[0] |
---|
2799 | atom[1] = 'Va' |
---|
2800 | Natoms.append(atom) |
---|
2801 | i += ldup |
---|
2802 | else: |
---|
2803 | i += 1 |
---|
2804 | else: |
---|
2805 | Natoms = Atoms |
---|
2806 | |
---|
2807 | NAtype = np.zeros(len(Atseq)) |
---|
2808 | for atom in Natoms: |
---|
2809 | NAtype[Atseq.index(atom[1])] += 1 |
---|
2810 | NAstr = ['%6d'%i for i in NAtype] |
---|
2811 | Cell = newPhase['General']['Cell'][1:7] |
---|
2812 | if os.path.exists(Name+'.his6f'): |
---|
2813 | os.remove(Name+'.his6f') |
---|
2814 | if os.path.exists(Name+'.neigh'): |
---|
2815 | os.remove(Name+'.neigh') |
---|
2816 | fname = Name+'.rmc6f' |
---|
2817 | fl = open(fname,'w') |
---|
2818 | fl.write('(Version 6f format configuration file)\n') |
---|
2819 | for item in Meta: |
---|
2820 | fl.write('%-20s%s\n'%('Metadata '+item+':',Meta[item])) |
---|
2821 | fl.write('Atom types present: %s\n'%' '.join(Atseq)) |
---|
2822 | fl.write('Number of each atom type: %s\n'%''.join(NAstr)) |
---|
2823 | fl.write('Number of atoms: %d\n'%len(Natoms)) |
---|
2824 | fl.write('%-35s%4d%4d%4d\n'%('Supercell dimensions:',Supercell[0],Supercell[1],Supercell[2])) |
---|
2825 | fl.write('Cell (Ang/deg): %12.6f%12.6f%12.6f%12.6f%12.6f%12.6f\n'%( |
---|
2826 | Cell[0],Cell[1],Cell[2],Cell[3],Cell[4],Cell[5])) |
---|
2827 | A,B = G2lat.cell2AB(Cell,True) |
---|
2828 | fl.write('Lattice vectors (Ang):\n') |
---|
2829 | for i in [0,1,2]: |
---|
2830 | fl.write('%12.6f%12.6f%12.6f\n'%(A[i,0],A[i,1],A[i,2])) |
---|
2831 | fl.write('Atoms (fractional coordinates):\n') |
---|
2832 | nat = 0 |
---|
2833 | for atm in Atseq: |
---|
2834 | for iat,atom in enumerate(Natoms): |
---|
2835 | if atom[1] == atm: |
---|
2836 | nat += 1 |
---|
2837 | atcode = Atcodes[iat].split(':') |
---|
2838 | cell = [0,0,0] |
---|
2839 | if '+' in atcode[1]: |
---|
2840 | cell = eval(atcode[1].split('+')[1]) |
---|
2841 | fl.write('%6d%4s [%s]%19.15f%19.15f%19.15f%6d%4d%4d%4d\n'%( |
---|
2842 | nat,atom[1].strip(),atcode[0],atom[3],atom[4],atom[5],(iat)%Natm+1,cell[0],cell[1],cell[2])) |
---|
2843 | fl.close() |
---|
2844 | return fname,reset |
---|
2845 | |
---|
2846 | def MakeBragg(PWDdata,Name,Phase): |
---|
2847 | generalData = Phase['General'] |
---|
2848 | Vol = generalData['Cell'][7] |
---|
2849 | Data = PWDdata['Data'] |
---|
2850 | Inst = PWDdata['Instrument Parameters'][0] |
---|
2851 | Bank = int(Inst['Bank'][1]) |
---|
2852 | Sample = PWDdata['Sample Parameters'] |
---|
2853 | Scale = Sample['Scale'][0] |
---|
2854 | Limits = PWDdata['Limits'][1] |
---|
2855 | Ibeg = np.searchsorted(Data[0],Limits[0]) |
---|
2856 | Ifin = np.searchsorted(Data[0],Limits[1])+1 |
---|
2857 | fname = Name+'.bragg' |
---|
2858 | fl = open(fname,'w') |
---|
2859 | fl.write('%12d%6d%15.7f%15.4f\n'%(Ifin-Ibeg-2,Bank,Scale,Vol)) |
---|
2860 | if 'T' in Inst['Type'][0]: |
---|
2861 | fl.write('%12s%12s\n'%(' TOF,ms',' I(obs)')) |
---|
2862 | for i in range(Ibeg,Ifin-1): |
---|
2863 | fl.write('%12.8f%12.6f\n'%(Data[0][i]/1000.,Data[1][i])) |
---|
2864 | else: |
---|
2865 | fl.write('%12s%12s\n'%(' 2-theta, deg',' I(obs)')) |
---|
2866 | for i in range(Ibeg,Ifin-1): |
---|
2867 | fl.write('%11.6f%15.2f\n'%(Data[0][i],Data[1][i])) |
---|
2868 | fl.close() |
---|
2869 | return fname |
---|
2870 | |
---|
2871 | def MakeRMCPdat(PWDdata,Name,Phase,RMCPdict): |
---|
2872 | Meta = RMCPdict['metadata'] |
---|
2873 | Times = RMCPdict['runTimes'] |
---|
2874 | Atseq = RMCPdict['atSeq'] |
---|
2875 | Natoms = RMCPdict['NoAtoms'] |
---|
2876 | sumatms = np.sum(np.array([Natoms[iatm] for iatm in Natoms])) |
---|
2877 | Isotope = RMCPdict['Isotope'] |
---|
2878 | Isotopes = RMCPdict['Isotopes'] |
---|
2879 | Atypes = RMCPdict['aTypes'] |
---|
2880 | if 'Va' in Atypes: |
---|
2881 | Isotope['Va'] = 'Nat. Abund.' |
---|
2882 | Isotopes['Va'] = {'Nat. Abund.':{'SL':[0.0,0.0]}} |
---|
2883 | atPairs = RMCPdict['Pairs'] |
---|
2884 | Files = RMCPdict['files'] |
---|
2885 | BraggWt = RMCPdict['histogram'][1] |
---|
2886 | inst = PWDdata['Instrument Parameters'][0] |
---|
2887 | try: |
---|
2888 | refList = PWDdata['Reflection Lists'][Name]['RefList'] |
---|
2889 | except TypeError: |
---|
2890 | return 'Error - missing reflection list; you must do Refine first' |
---|
2891 | dMin = refList[-1][4] |
---|
2892 | gsasType = 'xray2' |
---|
2893 | if 'T' in inst['Type'][1]: |
---|
2894 | gsasType = 'gsas3' |
---|
2895 | elif 'X' in inst['Type'][1]: |
---|
2896 | XFF = G2elem.GetFFtable(Atseq) |
---|
2897 | Xfl = open(Name+'.xray','w') |
---|
2898 | for atm in Atseq: |
---|
2899 | fa = XFF[atm]['fa'] |
---|
2900 | fb = XFF[atm]['fb'] |
---|
2901 | fc = XFF[atm]['fc'] |
---|
2902 | Xfl.write('%2s %8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f%8.4f\n'%( |
---|
2903 | atm.upper(),fa[0],fb[0],fa[1],fb[1],fa[2],fb[2],fa[3],fb[3],fc)) |
---|
2904 | Xfl.close() |
---|
2905 | lenA = len(Atseq) |
---|
2906 | Pairs = [] |
---|
2907 | Ncoeff = [] |
---|
2908 | Nblen = [Isotopes[at][Isotope[at]]['SL'][0] for at in Atypes] |
---|
2909 | for pair in [[' %s-%s'%(Atseq[i],Atseq[j]) for j in range(i,lenA)] for i in range(lenA)]: |
---|
2910 | Pairs += pair |
---|
2911 | for pair in Pairs: |
---|
2912 | pair = pair.replace(' ','') |
---|
2913 | at1,at2 = pair.split('-') |
---|
2914 | if at1 == 'Va' or at2 == 'Va': |
---|
2915 | ncoef = 0.0 |
---|
2916 | else: |
---|
2917 | ncoef = Isotopes[at1][Isotope[at1]]['SL'][0]*Natoms[at1]/sumatms |
---|
2918 | ncoef *= Isotopes[at2][Isotope[at2]]['SL'][0]*Natoms[at2]/sumatms |
---|
2919 | if at1 != at2: |
---|
2920 | ncoef *= 2. |
---|
2921 | Ncoeff += [ncoef,] |
---|
2922 | pairMin = [atPairs[pair] if pair in atPairs else [0.0,0.,0.] for pair in Pairs ] |
---|
2923 | maxMoves = [Atypes[atm] if atm in Atypes else 0.0 for atm in Atseq ] |
---|
2924 | fname = Name+'.dat' |
---|
2925 | fl = open(fname,'w') |
---|
2926 | fl.write(' %% Hand edit the following as needed\n') |
---|
2927 | fl.write('TITLE :: '+Name+'\n') |
---|
2928 | fl.write('MATERIAL :: '+Meta['material']+'\n') |
---|
2929 | fl.write('PHASE :: '+Meta['phase']+'\n') |
---|
2930 | fl.write('TEMPERATURE :: '+str(Meta['temperature'])+'\n') |
---|
2931 | fl.write('INVESTIGATOR :: '+Meta['owner']+'\n') |
---|
2932 | if RMCPdict.get('useGPU',False): |
---|
2933 | fl.write('GPU_ACCELERATOR :: 0\n') |
---|
2934 | minHD = ' '.join(['%6.3f'%dist[0] for dist in pairMin]) |
---|
2935 | minD = ' '.join(['%6.3f'%dist[1] for dist in pairMin]) |
---|
2936 | maxD = ' '.join(['%6.3f'%dist[2] for dist in pairMin]) |
---|
2937 | fl.write('MINIMUM_DISTANCES :: %s Angstrom\n'%minHD) |
---|
2938 | maxMv = ' '.join(['%6.3f'%mov for mov in maxMoves]) |
---|
2939 | fl.write('MAXIMUM_MOVES :: %s Angstrom\n'%maxMv) |
---|
2940 | fl.write('R_SPACING :: 0.0200 Angstrom\n') |
---|
2941 | fl.write('PRINT_PERIOD :: 100\n') |
---|
2942 | fl.write('TIME_LIMIT :: %.2f MINUTES\n'%Times[0]) |
---|
2943 | fl.write('SAVE_PERIOD :: %.2f MINUTES\n'%Times[1]) |
---|
2944 | fl.write('\n') |
---|
2945 | fl.write('ATOMS :: '+' '.join(Atseq)+'\n') |
---|
2946 | fl.write('\n') |
---|
2947 | fl.write('FLAGS ::\n') |
---|
2948 | fl.write(' > NO_MOVEOUT\n') |
---|
2949 | fl.write(' > NO_SAVE_CONFIGURATIONS\n') |
---|
2950 | fl.write(' > NO_RESOLUTION_CONVOLUTION\n') |
---|
2951 | fl.write('\n') |
---|
2952 | fl.write('INPUT_CONFIGURATION_FORMAT :: rmc6f\n') |
---|
2953 | fl.write('SAVE_CONFIGURATION_FORMAT :: rmc6f\n') |
---|
2954 | fl.write('IGNORE_HISTORY_FILE ::\n') |
---|
2955 | fl.write('\n') |
---|
2956 | if 'T' in inst['Type'][1]: |
---|
2957 | fl.write('NEUTRON_COEFFICIENTS :: '+''.join(['%9.5f'%coeff for coeff in Ncoeff])+'\n') |
---|
2958 | fl.write('DISTANCE_WINDOW ::\n') |
---|
2959 | fl.write(' > MNDIST :: %s\n'%minD) |
---|
2960 | fl.write(' > MXDIST :: %s\n'%maxD) |
---|
2961 | if len(RMCPdict['Potentials']['Stretch']) or len(RMCPdict['Potentials']['Stretch']): |
---|
2962 | fl.write('\n') |
---|
2963 | fl.write('POTENTIALS ::\n') |
---|
2964 | fl.write(' > TEMPERATURE :: %.1f K\n'%RMCPdict['Potentials']['Pot. Temp.']) |
---|
2965 | fl.write(' > PLOT :: pixels=400, colour=red, zangle=90, zrotation=45 deg\n') |
---|
2966 | if len(RMCPdict['Potentials']['Stretch']): |
---|
2967 | fl.write(' > STRETCH_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Stretch search']) |
---|
2968 | for bond in RMCPdict['Potentials']['Stretch']: |
---|
2969 | fl.write(' > STRETCH :: %s %s %.2f eV %.2f Ang\n'%(bond[0],bond[1],bond[3],bond[2])) |
---|
2970 | if len(RMCPdict['Potentials']['Angles']): |
---|
2971 | fl.write(' > ANGLE_SEARCH :: %.1f%%\n'%RMCPdict['Potentials']['Angle search']) |
---|
2972 | for angle in RMCPdict['Potentials']['Angles']: |
---|
2973 | fl.write(' > ANGLE :: %s %s %s %.2f eV %.2f deg %.2f %.2f Ang\n'% |
---|
2974 | (angle[1],angle[0],angle[2],angle[6],angle[3],angle[4],angle[5])) |
---|
2975 | if RMCPdict['useBVS']: |
---|
2976 | fl.write('BVS ::\n') |
---|
2977 | fl.write(' > ATOM :: '+' '.join(Atseq)+'\n') |
---|
2978 | fl.write(' > WEIGHTS :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']])) |
---|
2979 | oxid = [] |
---|
2980 | for val in RMCPdict['Oxid']: |
---|
2981 | if len(val) == 3: |
---|
2982 | oxid.append(val[0][1:]) |
---|
2983 | else: |
---|
2984 | oxid.append(val[0][2:]) |
---|
2985 | fl.write(' > OXID :: %s\n'%' '.join(oxid)) |
---|
2986 | fl.write(' > RIJ :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][0] for bvs in RMCPdict['BVS']])) |
---|
2987 | fl.write(' > BVAL :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][1] for bvs in RMCPdict['BVS']])) |
---|
2988 | fl.write(' > CUTOFF :: %s\n'%' '.join(['%6.3f'%RMCPdict['BVS'][bvs][2] for bvs in RMCPdict['BVS']])) |
---|
2989 | fl.write(' > SAVE :: 100000\n') |
---|
2990 | fl.write(' > UPDATE :: 100000\n') |
---|
2991 | if len(RMCPdict['Swaps']): |
---|
2992 | fl.write('\n') |
---|
2993 | fl.write('SWAP_MULTI ::\n') |
---|
2994 | for swap in RMCPdict['Swaps']: |
---|
2995 | try: |
---|
2996 | at1 = Atseq.index(swap[0]) |
---|
2997 | at2 = Atseq.index(swap[1]) |
---|
2998 | except ValueError: |
---|
2999 | break |
---|
3000 | fl.write(' > SWAP_ATOMS :: %d %d %.2f\n'%(at1,at2,swap[2])) |
---|
3001 | |
---|
3002 | if len(RMCPdict['FxCN']): |
---|
3003 | fl.write('FIXED_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['FxCN'])) |
---|
3004 | for ifx,fxcn in enumerate(RMCPdict['FxCN']): |
---|
3005 | try: |
---|
3006 | at1 = Atseq.index(fxcn[0]) |
---|
3007 | at2 = Atseq.index(fxcn[1]) |
---|
3008 | except ValueError: |
---|
3009 | break |
---|
3010 | fl.write(' > CSTR%d :: %d %d %.2f %.2f %.2f %.2f %.6f\n'%(ifx+1,at1+1,at2+1,fxcn[2],fxcn[3],fxcn[4],fxcn[5],fxcn[6])) |
---|
3011 | if len(RMCPdict['AveCN']): |
---|
3012 | fl.write('AVERAGE_COORDINATION_CONSTRAINTS :: %d\n'%len(RMCPdict['AveCN'])) |
---|
3013 | for iav,avcn in enumerate(RMCPdict['AveCN']): |
---|
3014 | try: |
---|
3015 | at1 = Atseq.index(avcn[0]) |
---|
3016 | at2 = Atseq.index(avcn[1]) |
---|
3017 | except ValueError: |
---|
3018 | break |
---|
3019 | fl.write(' > CAVSTR%d :: %d %d %.2f %.2f %.2f %.6f\n'%(iav+1,at1+1,at2+1,avcn[2],avcn[3],avcn[4],avcn[5])) |
---|
3020 | for File in Files: |
---|
3021 | if Files[File][0] and Files[File][0] != 'Select': |
---|
3022 | if 'Xray' in File and 'F(Q)' in File: |
---|
3023 | fqdata = open(Files[File][0],'r') |
---|
3024 | lines = int(fqdata.readline()[:-1]) |
---|
3025 | fqdata.close() |
---|
3026 | fl.write('\n') |
---|
3027 | fl.write('%s ::\n'%File.split(';')[0].upper().replace(' ','_')) |
---|
3028 | fl.write(' > FILENAME :: %s\n'%Files[File][0]) |
---|
3029 | fl.write(' > DATA_TYPE :: %s\n'%Files[File][2]) |
---|
3030 | fl.write(' > FIT_TYPE :: %s\n'%Files[File][2]) |
---|
3031 | if 'Xray' not in File: |
---|
3032 | fl.write(' > START_POINT :: 1\n') |
---|
3033 | fl.write(' > END_POINT :: 3000\n') |
---|
3034 | fl.write(' > WEIGHT :: %.4f\n'%Files[File][1]) |
---|
3035 | fl.write(' > CONSTANT_OFFSET 0.000\n') |
---|
3036 | fl.write(' > NO_FITTED_OFFSET\n') |
---|
3037 | if RMCPdict['FitScale']: |
---|
3038 | fl.write(' > FITTED_SCALE\n') |
---|
3039 | else: |
---|
3040 | fl.write(' > NO_FITTED_SCALE\n') |
---|
3041 | if Files[File][3] !='RMC': |
---|
3042 | fl.write(' > %s\n'%Files[File][3]) |
---|
3043 | if 'reciprocal' in File: |
---|
3044 | fl.write(' > CONVOLVE ::\n') |
---|
3045 | if 'Xray' in File: |
---|
3046 | fl.write(' > RECIPROCAL_SPACE_FIT :: 1 %d 1\n'%lines) |
---|
3047 | fl.write(' > RECIPROCAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(lines,Files[File][1])) |
---|
3048 | fl.write(' > REAL_SPACE_FIT :: 1 %d 1\n'%(3*lines//2)) |
---|
3049 | fl.write(' > REAL_SPACE_PARAMETERS :: 1 %d %.4f\n'%(3*lines//2,1./Files[File][1])) |
---|
3050 | fl.write('\n') |
---|
3051 | fl.write('BRAGG ::\n') |
---|
3052 | fl.write(' > BRAGG_SHAPE :: %s\n'%gsasType) |
---|
3053 | fl.write(' > RECALCUATE\n') |
---|
3054 | fl.write(' > DMIN :: %.2f\n'%(dMin-0.02)) |
---|
3055 | fl.write(' > WEIGHT :: %10.3f\n'%BraggWt) |
---|
3056 | if 'T' in inst['Type'][1]: |
---|
3057 | fl.write(' > SCATTERING LENGTH :: '+''.join(['%8.4f'%blen for blen in Nblen])+'\n') |
---|
3058 | fl.write('\n') |
---|
3059 | fl.write('END ::\n') |
---|
3060 | fl.close() |
---|
3061 | return fname |
---|
3062 | |
---|
3063 | # def FindBonds(Phase,RMCPdict): |
---|
3064 | # generalData = Phase['General'] |
---|
3065 | # cx,ct,cs,cia = generalData['AtomPtrs'] |
---|
3066 | # atomData = Phase['Atoms'] |
---|
3067 | # Res = 'RMC' |
---|
3068 | # if 'macro' in generalData['Type']: |
---|
3069 | # Res = atomData[0][ct-3] |
---|
3070 | # AtDict = {atom[ct-1]:atom[ct] for atom in atomData} |
---|
3071 | # Pairs = RMCPdict['Pairs'] #dict! |
---|
3072 | # BondList = [] |
---|
3073 | # notNames = [] |
---|
3074 | # for FrstName in AtDict: |
---|
3075 | # nbrs = G2mth.FindAllNeighbors(Phase,FrstName,list(AtDict.keys()),notName=notNames,Short=True)[0] |
---|
3076 | # Atyp1 = AtDict[FrstName] |
---|
3077 | # if 'Va' in Atyp1: |
---|
3078 | # continue |
---|
3079 | # for nbr in nbrs: |
---|
3080 | # Atyp2 = AtDict[nbr[0]] |
---|
3081 | # if 'Va' in Atyp2: |
---|
3082 | # continue |
---|
3083 | # try: |
---|
3084 | # bndData = Pairs[' %s-%s'%(Atyp1,Atyp2)][1:] |
---|
3085 | # except KeyError: |
---|
3086 | # bndData = Pairs[' %s-%s'%(Atyp2,Atyp1)][1:] |
---|
3087 | # if any(bndData): |
---|
3088 | # if bndData[0] <= nbr[1] <= bndData[1]: |
---|
3089 | # bondStr = str((FrstName,nbr[0])+tuple(bndData))+',\n' |
---|
3090 | # revbondStr = str((nbr[0],FrstName)+tuple(bndData))+',\n' |
---|
3091 | # if bondStr not in BondList and revbondStr not in BondList: |
---|
3092 | # BondList.append(bondStr) |
---|
3093 | # notNames.append(FrstName) |
---|
3094 | # return Res,BondList |
---|
3095 | |
---|
3096 | # def FindAngles(Phase,RMCPdict): |
---|
3097 | # generalData = Phase['General'] |
---|
3098 | # Cell = generalData['Cell'][1:7] |
---|
3099 | # Amat = G2lat.cell2AB(Cell)[0] |
---|
3100 | # cx,ct,cs,cia = generalData['AtomPtrs'] |
---|
3101 | # atomData = Phase['Atoms'] |
---|
3102 | # AtLookup = G2mth.FillAtomLookUp(atomData,cia+8) |
---|
3103 | # AtDict = {atom[ct-1]:atom[ct] for atom in atomData} |
---|
3104 | # Angles = RMCPdict['Angles'] |
---|
3105 | # AngDict = {'%s-%s-%s'%(angle[0],angle[1],angle[2]):angle[3:] for angle in Angles} |
---|
3106 | # AngleList = [] |
---|
3107 | # for MidName in AtDict: |
---|
3108 | # nbrs,nbrIds = G2mth.FindAllNeighbors(Phase,MidName,list(AtDict.keys()),Short=True) |
---|
3109 | # if len(nbrs) < 2: #need 2 neighbors to make an angle |
---|
3110 | # continue |
---|
3111 | # Atyp2 = AtDict[MidName] |
---|
3112 | # for i,nbr1 in enumerate(nbrs): |
---|
3113 | # Atyp1 = AtDict[nbr1[0]] |
---|
3114 | # for j,nbr3 in enumerate(nbrs[i+1:]): |
---|
3115 | # Atyp3 = AtDict[nbr3[0]] |
---|
3116 | # IdList = [nbrIds[1][i],nbrIds[0],nbrIds[1][i+j+1]] |
---|
3117 | # try: |
---|
3118 | # angData = AngDict['%s-%s-%s'%(Atyp1,Atyp2,Atyp3)] |
---|
3119 | # except KeyError: |
---|
3120 | # try: |
---|
3121 | # angData = AngDict['%s-%s-%s'%(Atyp3,Atyp2,Atyp1)] |
---|
3122 | # except KeyError: |
---|
3123 | # continue |
---|
3124 | # XYZ = np.array(G2mth.GetAtomItemsById(atomData,AtLookup,IdList,cx,numItems=3)) |
---|
3125 | # calAngle = G2mth.getRestAngle(XYZ,Amat) |
---|
3126 | # if angData[0] <= calAngle <= angData[1]: |
---|
3127 | # angStr = str((MidName,nbr1[0],nbr3[0])+tuple(angData))+',\n' |
---|
3128 | # revangStr = str((MidName,nbr3[0],nbr1[0])+tuple(angData))+',\n' |
---|
3129 | # if angStr not in AngleList and revangStr not in AngleList: |
---|
3130 | # AngleList.append(angStr) |
---|
3131 | # return AngleList |
---|
3132 | |
---|
3133 | # def GetSqConvolution(XY,d): |
---|
3134 | |
---|
3135 | # n = XY.shape[1] |
---|
3136 | # snew = np.zeros(n) |
---|
3137 | # dq = np.zeros(n) |
---|
3138 | # sold = XY[1] |
---|
3139 | # q = XY[0] |
---|
3140 | # dq[1:] = np.diff(q) |
---|
3141 | # dq[0] = dq[1] |
---|
3142 | |
---|
3143 | # for j in range(n): |
---|
3144 | # for i in range(n): |
---|
3145 | # b = abs(q[i]-q[j]) |
---|
3146 | # t = q[i]+q[j] |
---|
3147 | # if j == i: |
---|
3148 | # snew[j] += q[i]*sold[i]*(d-np.sin(t*d)/t)*dq[i] |
---|
3149 | # else: |
---|
3150 | # snew[j] += q[i]*sold[i]*(np.sin(b*d)/b-np.sin(t*d)/t)*dq[i] |
---|
3151 | # snew[j] /= np.pi*q[j] |
---|
3152 | |
---|
3153 | # snew[0] = snew[1] |
---|
3154 | # return snew |
---|
3155 | |
---|
3156 | # def GetMaxSphere(pdbName): |
---|
3157 | # try: |
---|
3158 | # pFil = open(pdbName,'r') |
---|
3159 | # except FileNotFoundError: |
---|
3160 | # return None |
---|
3161 | # while True: |
---|
3162 | # line = pFil.readline() |
---|
3163 | # if 'Boundary' in line: |
---|
3164 | # line = line.split()[3:] |
---|
3165 | # G = np.array([float(item) for item in line]) |
---|
3166 | # G = np.reshape(G,(3,3))**2 |
---|
3167 | # G = nl.inv(G) |
---|
3168 | # pFil.close() |
---|
3169 | # break |
---|
3170 | # dspaces = [0.5/np.sqrt(G2lat.calc_rDsq2(H,G)) for H in np.eye(3)] |
---|
3171 | # return min(dspaces) |
---|
3172 | |
---|
3173 | def findfullrmc(): |
---|
3174 | '''Find where fullrmc is installed. Tries the following: |
---|
3175 | |
---|
3176 | 1. Returns the Config var 'fullrmc_exec', if defined. If an executable |
---|
3177 | is found at that location it is assumed to run and supply |
---|
3178 | fullrmc 5.0+ |
---|
3179 | 2. The path is checked for a fullrmc image as named by Bachir |
---|
3180 | |
---|
3181 | :returns: the full path to a python executable that is assumed to |
---|
3182 | have fullrmc installed or None, if it was not found. |
---|
3183 | ''' |
---|
3184 | fullrmc_exe = GSASIIpath.GetConfigValue('fullrmc_exec') |
---|
3185 | if fullrmc_exe is not None and is_exe(fullrmc_exe): |
---|
3186 | return fullrmc_exe |
---|
3187 | pathlist = os.environ["PATH"].split(os.pathsep) |
---|
3188 | for p in (GSASIIpath.path2GSAS2,GSASIIpath.binaryPath,os.getcwd(), |
---|
3189 | os.path.split(sys.executable)[0]): |
---|
3190 | if p not in pathlist: pathlist.append(p) |
---|
3191 | import glob |
---|
3192 | for p in pathlist: |
---|
3193 | if sys.platform == "win32": |
---|
3194 | lookfor = "fullrmc5*.exe" |
---|
3195 | else: |
---|
3196 | lookfor = "fullrmc5*64bit" |
---|
3197 | fl = glob.glob(os.path.join(p,lookfor)) |
---|
3198 | if len(fl) > 0: |
---|
3199 | fullrmc_exe = os.path.abspath(sorted(fl)[0]) |
---|
3200 | if GSASIIpath.GetConfigValue('debug'): |
---|
3201 | print('fullrmc found as',fullrmc_exe) |
---|
3202 | return fullrmc_exe |
---|
3203 | |
---|
3204 | def fullrmcDownload(): |
---|
3205 | '''Downloads the fullrmc executable from Bachir's site to the current |
---|
3206 | GSAS-II binary directory. |
---|
3207 | |
---|
3208 | Does some error checking. |
---|
3209 | ''' |
---|
3210 | import os |
---|
3211 | import requests |
---|
3212 | import platform |
---|
3213 | if platform.architecture()[0] != '64bit': |
---|
3214 | return "fullrmc is only available for 64 bit machines. This is 32 bit" |
---|
3215 | setXbit = True |
---|
3216 | if sys.platform == "darwin": |
---|
3217 | URL = "https://github.com/bachiraoun/fullrmc/raw/master/standalones/fullrmc500_3p8p6_macOS-10p16-x86_64-i386-64bit" |
---|
3218 | elif sys.platform == "win32": |
---|
3219 | setXbit = False |
---|
3220 | URL = "https://github.com/bachiraoun/fullrmc/raw/master/standalones/fullrmc500_3p8p10_Windows-10-10p0p19041-SP0.exe" |
---|
3221 | else: |
---|
3222 | if 'aarch' in platform.machine() or 'arm' in platform.machine(): |
---|
3223 | return "Sorry, fullrmc is only available for Intel-compatible machines." |
---|
3224 | URL = "https://github.com/bachiraoun/fullrmc/raw/master/standalones/fullrmc500_3p8p5_Linux-4p19p121-linuxkit-x86_64-with-glibc2p29" |
---|
3225 | |
---|
3226 | GSASIIpath.SetBinaryPath() |
---|
3227 | fil = os.path.join(GSASIIpath.binaryPath,os.path.split(URL)[1]) |
---|
3228 | print('Starting installation of fullrmc\nDownloading from', |
---|
3229 | 'https://github.com/bachiraoun/fullrmc/tree/master/standalones', |
---|
3230 | '\nCreating '+fil, |
---|
3231 | '\nThis may take a while...') |
---|
3232 | open(fil, "wb").write(requests.get(URL).content) |
---|
3233 | print('...Download completed') |
---|
3234 | if setXbit: |
---|
3235 | import stat |
---|
3236 | os.chmod(fil, os.stat(fil).st_mode | stat.S_IEXEC) |
---|
3237 | return '' |
---|
3238 | |
---|
3239 | def findPDFfit(): |
---|
3240 | '''Find if PDFfit2 is installed (may be local to GSAS-II). Does the following: |
---|
3241 | :returns: two items: (1) the full path to a python executable or None, if |
---|
3242 | it was not found and (2) path(s) to the PDFfit2 location(s) as a list. |
---|
3243 | |
---|
3244 | ''' |
---|
3245 | if GSASIIpath.GetConfigValue('pdffit2_exec') is not None and is_exe( |
---|
3246 | GSASIIpath.GetConfigValue('pdffit2_exec')): |
---|
3247 | return GSASIIpath.GetConfigValue('pdffit2_exec'),None |
---|
3248 | pdffitloc = os.path.join(GSASIIpath.path2GSAS2,'PDFfit2') |
---|
3249 | if not os.path.exists(pdffitloc): |
---|
3250 | print('PDFfit2 not found in GSAS-II \n\t(expected in '+pdffitloc+')') |
---|
3251 | return None,[] |
---|
3252 | if pdffitloc not in sys.path: sys.path.append(pdffitloc) |
---|
3253 | try: |
---|
3254 | from diffpy.pdffit2 import PdfFit |
---|
3255 | import diffpy |
---|
3256 | import inspect |
---|
3257 | pdffitloc = [os.path.dirname(os.path.dirname(inspect.getfile(diffpy)))] |
---|
3258 | # is this the original version of diffpy (w/pdffit2.py) |
---|
3259 | try: |
---|
3260 | from diffpy.pdffit2 import pdffit2 |
---|
3261 | except ImportError: |
---|
3262 | # or the GSAS-II version w/o; for this we need to find the binary's location |
---|
3263 | try: |
---|
3264 | import pdffit2 # added for GSAS-II to relocate binary file |
---|
3265 | except ImportError: |
---|
3266 | print('\nError: pdffit2 failed to load with this python\n') |
---|
3267 | return None,[] |
---|
3268 | except ModuleNotFoundError: |
---|
3269 | print('\nGSAS-II does not have a PDFfit2 module compatible\nwith this Python interpreter\n') |
---|
3270 | return None,[] |
---|
3271 | pdffitloc += [os.path.dirname(inspect.getfile(pdffit2))] |
---|
3272 | return sys.executable,pdffitloc |
---|
3273 | except Exception as msg: |
---|
3274 | print('Error importing PDFfit2:\n',msg) |
---|
3275 | return None,[] |
---|
3276 | |
---|
3277 | def GetPDFfitAtomVar(Phase,RMCPdict): |
---|
3278 | ''' Find dict of independent "@n" variables for PDFfit in atom constraints |
---|
3279 | ''' |
---|
3280 | General = Phase['General'] |
---|
3281 | Atoms = Phase['Atoms'] |
---|
3282 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
3283 | AtomVar = RMCPdict['AtomVar'] |
---|
3284 | varnames = [] |
---|
3285 | for iat,atom in enumerate(RMCPdict['AtomConstr']): |
---|
3286 | for it,item in enumerate(atom): |
---|
3287 | if it > 1 and item: |
---|
3288 | itms = item.split('@') |
---|
3289 | for itm in itms[1:]: |
---|
3290 | itnum = itm[:2] |
---|
3291 | varname = '@%s'%itnum |
---|
3292 | varnames.append(varname) |
---|
3293 | if it < 6: |
---|
3294 | if varname not in AtomVar: |
---|
3295 | AtomVar[varname] = 0.0 #put ISODISTORT mode displ here? |
---|
3296 | else: |
---|
3297 | for i in range(3): |
---|
3298 | if varname not in AtomVar: |
---|
3299 | AtomVar[varname] = Atoms[iat][cia+i+2] |
---|
3300 | varnames = set(varnames) |
---|
3301 | for name in list(AtomVar.keys()): #clear out unused parameters |
---|
3302 | if name not in varnames: |
---|
3303 | del AtomVar[name] |
---|
3304 | |
---|
3305 | def MakePDFfitAtomsFile(Phase,RMCPdict): |
---|
3306 | '''Make the PDFfit atoms file |
---|
3307 | ''' |
---|
3308 | General = Phase['General'] |
---|
3309 | if General['SGData']['SpGrp'] != 'P 1': |
---|
3310 | return 'Space group symmetry must be lowered to P 1 for PDFfit' |
---|
3311 | fName = General['Name']+'-PDFfit.stru' |
---|
3312 | fName = fName.replace(' ','_') |
---|
3313 | if 'sequential' in RMCPdict['refinement']: |
---|
3314 | fName = 'Sequential_PDFfit.stru' |
---|
3315 | fatm = open(fName,'w') |
---|
3316 | fatm.write('title structure of '+General['Name']+'\n') |
---|
3317 | fatm.write('format pdffit\n') |
---|
3318 | fatm.write('scale 1.000000\n') #fixed |
---|
3319 | sharp = '%10.6f,%10.6f,%10.6f,%10.6f\n'%(RMCPdict['delta2'][0],RMCPdict['delta1'][0],RMCPdict['sratio'][0],RMCPdict['rcut']) |
---|
3320 | fatm.write('sharp '+sharp) |
---|
3321 | shape = '' |
---|
3322 | if RMCPdict['shape'] == 'sphere' and RMCPdict['spdiameter'][0] > 0.: |
---|
3323 | shape = ' sphere, %10.6f\n'%RMCPdict['spdiameter'][0] |
---|
3324 | elif RMCPdict['stepcut'] > 0.: |
---|
3325 | shape = 'stepcut, %10.6f\n'%RMCPdict['stepcut'] |
---|
3326 | if shape: |
---|
3327 | fatm.write('shape '+shape) |
---|
3328 | fatm.write('spcgr %s\n'%RMCPdict['SGData']['SpGrp'].replace(' ','')) |
---|
3329 | cell = General['Cell'][1:7] |
---|
3330 | fatm.write('cell %10.6f,%10.6f,%10.6f,%10.6f,%10.6f,%10.6f\n'%( |
---|
3331 | cell[0],cell[1],cell[2],cell[3],cell[4],cell[5])) |
---|
3332 | fatm.write('dcell '+5*' 0.000000,'+' 0.000000\n') |
---|
3333 | Atoms = Phase['Atoms'] |
---|
3334 | fatm.write('ncell %8d,%8d,%8d,%10d\n'%(1,1,1,len(Atoms))) |
---|
3335 | fatm.write('atoms\n') |
---|
3336 | cx,ct,cs,cia = General['AtomPtrs'] |
---|
3337 | for atom in Atoms: |
---|
3338 | fatm.write('%4s%18.8f%18.8f%18.8f%13.4f\n'%(atom[ct][:2].ljust(2),atom[cx],atom[cx+1],atom[cx+2],atom[cx+3])) |
---|
3339 | fatm.write(' '+'%18.8f%18.8f%18.8f%13.4f\n'%(0.,0.,0.,0.)) |
---|
3340 | fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(atom[cia+2],atom[cia+3],atom[cia+4])) |
---|
3341 | fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.,)) |
---|
3342 | fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(atom[cia+5],atom[cia+6],atom[cia+7])) |
---|
3343 | fatm.write(' '+'%18.8f%18.8f%18.8f\n'%(0.,0.,0.)) |
---|
3344 | fatm.close() |
---|
3345 | |
---|
3346 | def MakePDFfitRunFile(Phase,RMCPdict): |
---|
3347 | '''Make the PDFfit python run file |
---|
3348 | ''' |
---|
3349 | |
---|
3350 | def GetCellConstr(SGData): |
---|
3351 | if SGData['SGLaue'] in ['m3', 'm3m']: |
---|
3352 | return [1,1,1,0,0,0] |
---|
3353 | elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']: |
---|
3354 | return [1,1,2,0,0,0] |
---|
3355 | elif SGData['SGLaue'] in ['3R','3mR']: |
---|
3356 | return [1,1,1,2,2,2] |
---|
3357 | elif SGData['SGLaue'] == 'mmm': |
---|
3358 | return [1,2,3,0,0,0] |
---|
3359 | elif SGData['SGLaue'] == '2/m': |
---|
3360 | if SGData['SGUniq'] == 'a': |
---|
3361 | return [1,2,3,4,0,0] |
---|
3362 | elif SGData['SGUniq'] == 'b': |
---|
3363 | return [1,2,3,0,4,0] |
---|
3364 | elif SGData['SGUniq'] == 'c': |
---|
3365 | return [1,2,3,0,0,4] |
---|
3366 | else: |
---|
3367 | return [1,2,3,4,5,6] |
---|
3368 | |
---|
3369 | General = Phase['General'] |
---|
3370 | Cell = General['Cell'][1:7] |
---|
3371 | rundata = '''#!/usr/bin/env python |
---|
3372 | # -*- coding: utf-8 -*- |
---|
3373 | import sys,os |
---|
3374 | datadir = r'{:}' |
---|
3375 | pathWrap = lambda f: os.path.join(datadir,f) |
---|
3376 | '''.format(os.path.abspath(os.getcwd())) |
---|
3377 | PDFfit_exe,PDFfit_path = findPDFfit() # returns python loc and path(s) for pdffit |
---|
3378 | if not PDFfit_exe: |
---|
3379 | print('PDFfit2 is not found. Creating .sh file without paths.') |
---|
3380 | if PDFfit_path: |
---|
3381 | for p in PDFfit_path: |
---|
3382 | rundata += "sys.path.append(r'{:}')\n".format(p) |
---|
3383 | rundata += 'from diffpy.pdffit2 import PdfFit\n' |
---|
3384 | rundata += 'pf = PdfFit()\n' |
---|
3385 | Nd = 0 |
---|
3386 | Np = 0 |
---|
3387 | parms = {} |
---|
3388 | parmNames = {} |
---|
3389 | if 'sequential' in RMCPdict['refinement']: |
---|
3390 | Np = 3 |
---|
3391 | rundata += '#sequential data here\n' |
---|
3392 | else: |
---|
3393 | for fil in RMCPdict['files']: |
---|
3394 | filNam = RMCPdict['files'][fil][0] |
---|
3395 | if 'Select' in filNam: |
---|
3396 | continue |
---|
3397 | if 'Neutron' in fil: |
---|
3398 | Nd += 1 |
---|
3399 | dType = 'Ndata' |
---|
3400 | else: |
---|
3401 | Nd += 1 |
---|
3402 | dType = 'Xdata' |
---|
3403 | rundata += "pf.read_data(pathWrap(r'%s'), '%s', 30.0, %.4f)\n"%(filNam,dType[0],RMCPdict[dType]['qdamp'][0]) |
---|
3404 | rundata += 'pf.setdata(%d)\n'%Nd |
---|
3405 | rundata += 'pf.pdfrange(%d, %6.2f, %6.2f)\n'%(Nd,RMCPdict[dType]['Fitrange'][0],RMCPdict[dType]['Fitrange'][1]) |
---|
3406 | for item in ['dscale','qdamp','qbroad']: |
---|
3407 | if RMCPdict[dType][item][1]: |
---|
3408 | Np += 1 |
---|
3409 | rundata += 'pf.constrain(pf.%s(),"@%d")\n'%(item,Np) |
---|
3410 | parms[Np] = RMCPdict[dType][item][0] |
---|
3411 | parmNames[Np] = item |
---|
3412 | fName = General['Name']+'-PDFfit.stru' |
---|
3413 | fName = fName.replace(' ','_') |
---|
3414 | if 'sequential' in RMCPdict['refinement']: |
---|
3415 | fName = 'Sequential_PDFfit.stru' |
---|
3416 | Np = 9 |
---|
3417 | rundata += "pf.read_struct(pathWrap(r'{:}'))\n".format(fName) |
---|
3418 | for item in ['delta1','delta2','sratio']: |
---|
3419 | if RMCPdict[item][1]: |
---|
3420 | Np += 1 |
---|
3421 | rundata += 'pf.constrain(pf.%s,"@%d")\n'%(item,Np) |
---|
3422 | parms[Np] = RMCPdict[item][0] |
---|
3423 | parmNames[Np] = item |
---|
3424 | if 'sphere' in RMCPdict['shape'] and RMCPdict['spdiameter'][1]: |
---|
3425 | Np += 1 |
---|
3426 | rundata += 'pf.constrain(pf.spdiameter,"@%d")\n'%Np |
---|
3427 | parms[Np] = RMCPdict['spdiameter'][0] |
---|
3428 | parmNames[Np] = 'spdiameter' |
---|
3429 | |
---|
3430 | if RMCPdict['cellref']: |
---|
3431 | cellconst = GetCellConstr(RMCPdict['SGData']) |
---|
3432 | used = [] |
---|
3433 | cellNames = ['a','b','c','alpha','beta','gamma'] |
---|
3434 | for ic in range(6): |
---|
3435 | if cellconst[ic]: |
---|
3436 | rundata += 'pf.constrain(pf.lat(%d), "@%d")\n'%(ic+1,Np+cellconst[ic]) |
---|
3437 | if cellconst[ic] not in used: |
---|
3438 | parms[Np+cellconst[ic]] = Cell[ic] |
---|
3439 | parmNames[Np+cellconst[ic]] = cellNames[ic] |
---|
3440 | used.append(cellconst[ic]) |
---|
3441 | #Atom constraints here ------------------------------------------------------- |
---|
3442 | AtomVar = RMCPdict['AtomVar'] |
---|
3443 | used = [] |
---|
3444 | for iat,atom in enumerate(RMCPdict['AtomConstr']): |
---|
3445 | for it,item in enumerate(atom): |
---|
3446 | names = ['pf.x(%d)'%(iat+1),'pf.y(%d)'%(iat+1),'pf.z(%d)'%(iat+1),'pf.occ(%d)'%(iat+1)] |
---|
3447 | if it > 1 and item: |
---|
3448 | itms = item.split('@') |
---|
3449 | once = False |
---|
3450 | for itm in itms[1:]: |
---|
3451 | try: |
---|
3452 | itnum = int(itm[:2]) |
---|
3453 | except ValueError: |
---|
3454 | print(' *** ERROR - invalid string in atom constraint %s ***'%(item)) |
---|
3455 | return None |
---|
3456 | if it < 6: |
---|
3457 | if not once: |
---|
3458 | rundata += 'pf.constrain(%s,"%s")\n'%(names[it-2],item) |
---|
3459 | once = True |
---|
3460 | if itnum not in used: |
---|
3461 | parms[itnum] = AtomVar['@%d'%itnum] |
---|
3462 | parmNames[itnum] = names[it-2].split('.')[1] |
---|
3463 | used.append(itnum) |
---|
3464 | else: |
---|
3465 | uijs = ['pf.u11(%d)'%(iat+1),'pf.u22(%d)'%(iat+1),'pf.u33(%d)'%(iat+1)] |
---|
3466 | for i in range(3): |
---|
3467 | rundata += 'pf.constrain(%s,"%s")\n'%(uijs[i],item) |
---|
3468 | if itnum not in used: |
---|
3469 | parms[itnum] = AtomVar['@%d'%itnum] |
---|
3470 | parmNames[itnum] = uijs[i].split('.')[1] |
---|
3471 | used.append(itnum) |
---|
3472 | |
---|
3473 | if 'sequential' in RMCPdict['refinement']: |
---|
3474 | rundata += '#parameters here\n' |
---|
3475 | RMCPdict['Parms'] = parms #{'n':val,...} |
---|
3476 | RMCPdict['ParmNames'] = parmNames #{'n':name,...} |
---|
3477 | else: |
---|
3478 | # set parameter values |
---|
3479 | for iprm in parms: |
---|
3480 | rundata += 'pf.setpar(%d,%.6f)\n'%(iprm,parms[iprm]) |
---|
3481 | |
---|
3482 | # Save results --------------------------------------------------------------- |
---|
3483 | rundata += 'pf.refine()\n' |
---|
3484 | if 'sequential' in RMCPdict['refinement']: |
---|
3485 | fName = 'Sequential_PDFfit' |
---|
3486 | rfile = open('Seq_PDFfit_template.py','w') |
---|
3487 | rundata += 'pf.save_pdf(1, pathWrap("%s"))\n'%(fName+'.fgr') |
---|
3488 | else: |
---|
3489 | fName = General['Name'].replace(' ','_')+'-PDFfit' |
---|
3490 | rfile = open(fName+'.py','w') |
---|
3491 | Nd = 0 |
---|
3492 | for file in RMCPdict['files']: |
---|
3493 | if 'Select' in RMCPdict['files'][file][0]: #skip unselected |
---|
3494 | continue |
---|
3495 | Nd += 1 |
---|
3496 | rundata += 'pf.save_pdf(%d, pathWrap("%s"))\n'%(Nd,fName+file[0]+'.fgr') |
---|
3497 | |
---|
3498 | rundata += 'pf.save_struct(1, pathWrap("%s"))\n'%(fName+'.rstr') |
---|
3499 | rundata += 'pf.save_res(pathWrap("%s"))\n'%(fName+'.res') |
---|
3500 | |
---|
3501 | rfile.writelines(rundata) |
---|
3502 | rfile.close() |
---|
3503 | |
---|
3504 | return fName+'.py' |
---|
3505 | |
---|
3506 | def GetSeqCell(SGData,parmDict): |
---|
3507 | ''' For use in processing PDFfit sequential results |
---|
3508 | ''' |
---|
3509 | try: |
---|
3510 | if SGData['SGLaue'] in ['m3', 'm3m']: |
---|
3511 | cell = [parmDict['11'][0],parmDict['11'][0],parmDict['11'][0],90.,90.,90.] |
---|
3512 | elif SGData['SGLaue'] in ['3','3m1','31m','6/m','6/mmm','4/m','4/mmm']: |
---|
3513 | cell = [parmDict['11'][0],parmDict['11'][0],parmDict['12'][0],90.,90.,90.] |
---|
3514 | elif SGData['SGLaue'] in ['3R','3mR']: |
---|
3515 | cell = [parmDict['11'][0],parmDict['11'][0],parmDict['11'][0], |
---|
3516 | parmDict['12'][0],parmDict['12'][0],parmDict['12'][0]] |
---|
3517 | elif SGData['SGLaue'] == 'mmm': |
---|
3518 | cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,90.,90.] |
---|
3519 | elif SGData['SGLaue'] == '2/m': |
---|
3520 | if SGData['SGUniq'] == 'a': |
---|
3521 | cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],parmDict['14'][0],90.,90.] |
---|
3522 | elif SGData['SGUniq'] == 'b': |
---|
3523 | cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,parmDict['14'][0],90.] |
---|
3524 | elif SGData['SGUniq'] == 'c': |
---|
3525 | cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0],90.,90.,parmDict['14'][0]] |
---|
3526 | else: |
---|
3527 | cell = [parmDict['11'][0],parmDict['12'][0],parmDict['13'][0], |
---|
3528 | parmDict['14'][0],parmDict['15'][0],parmDict['16'][0]] |
---|
3529 | return G2lat.cell2A(cell) |
---|
3530 | except KeyError: |
---|
3531 | return None |
---|
3532 | |
---|
3533 | def UpdatePDFfit(Phase,RMCPdict): |
---|
3534 | ''' Updates various PDFfit parameters held in GSAS-II |
---|
3535 | ''' |
---|
3536 | |
---|
3537 | General = Phase['General'] |
---|
3538 | if RMCPdict['refinement'] == 'normal': |
---|
3539 | fName = General['Name']+'-PDFfit.rstr' |
---|
3540 | try: |
---|
3541 | rstr = open(fName.replace(' ','_'),'r') |
---|
3542 | except FileNotFoundError: |
---|
3543 | return [fName,'Not found - PDFfit failed'] |
---|
3544 | lines = rstr.readlines() |
---|
3545 | rstr.close() |
---|
3546 | header = [line[:-1].split(' ',1) for line in lines[:7]] |
---|
3547 | resdict = dict(header) |
---|
3548 | for item in ['sharp','cell']: |
---|
3549 | resdict[item] = [float(val) for val in resdict[item].split(',')] |
---|
3550 | General['Cell'][1:7] = resdict['cell'] |
---|
3551 | for inam,name in enumerate(['delta2','delta1','sratio']): |
---|
3552 | RMCPdict[name][0] = float(resdict['sharp'][inam]) |
---|
3553 | if 'shape' in resdict: |
---|
3554 | if 'sphere' in resdict['shape']: |
---|
3555 | RMCPdict['spdiameter'][0] = float(resdict['shape'].split()[-1]) |
---|
3556 | else: |
---|
3557 | RMCPdict['stepcut'][0] = float(resdict['shape'][-1]) |
---|
3558 | cx,ct,cs,ci = G2mth.getAtomPtrs(Phase) |
---|
3559 | Atoms = Phase['Atoms'] |
---|
3560 | atmBeg = 0 |
---|
3561 | for line in lines: |
---|
3562 | atmBeg += 1 |
---|
3563 | if 'atoms' in line: |
---|
3564 | break |
---|
3565 | for atom in Atoms: |
---|
3566 | atstr = lines[atmBeg][:-1].split() |
---|
3567 | Uiistr = lines[atmBeg+2][:-1].split() |
---|
3568 | Uijstr = lines[atmBeg+4][:-1].split() |
---|
3569 | atom[cx:cx+4] = [float(atstr[1]),float(atstr[2]),float(atstr[3]),float(atstr[4])] |
---|
3570 | atom[ci] = 'A' |
---|
3571 | atom[ci+2:ci+5] = [float(Uiistr[0]),float(Uiistr[1]),float(Uiistr[2])] |
---|
3572 | atom[ci+5:ci+8] = [float(Uijstr[0]),float(Uijstr[1]),float(Uijstr[2])] |
---|
3573 | atmBeg += 6 |
---|
3574 | fName = General['Name']+'-PDFfit.res' |
---|
3575 | else: |
---|
3576 | fName = 'Sequential_PDFfit.res' |
---|
3577 | try: |
---|
3578 | res = open(fName.replace(' ','_'),'r') |
---|
3579 | except FileNotFoundError: |
---|
3580 | return [fName,'Not found - PDFfit failed'] |
---|
3581 | lines = res.readlines() |
---|
3582 | res.close() |
---|
3583 | Ibeg = False |
---|
3584 | resline = '' |
---|
3585 | XNdata = {'Xdata':RMCPdict['Xdata'],'Ndata':RMCPdict['Ndata']} |
---|
3586 | for line in lines: |
---|
3587 | if 'Radiation' in line and 'X-Rays' in line: |
---|
3588 | dkey = 'Xdata' |
---|
3589 | if 'Radiation' in line and'Neutrons' in line: |
---|
3590 | dkey = 'Ndata' |
---|
3591 | if 'Qdamp' in line and '(' in line: |
---|
3592 | XNdata[dkey]['qdamp'][0] = float(line.split()[4]) |
---|
3593 | if 'Qbroad' in line and '(' in line: |
---|
3594 | XNdata[dkey]['qbroad'][0] = float(line.split()[4]) |
---|
3595 | if 'Scale' in line and '(' in line: |
---|
3596 | XNdata[dkey]['dscale'][0] = float(line.split()[3]) |
---|
3597 | |
---|
3598 | for iline,line in enumerate(lines): |
---|
3599 | if 'Refinement parameters' in line: |
---|
3600 | Ibeg = True |
---|
3601 | continue |
---|
3602 | if Ibeg: |
---|
3603 | if '---------' in line: |
---|
3604 | break |
---|
3605 | resline += line[:-1] |
---|
3606 | for iline,line in enumerate(lines): |
---|
3607 | if 'Rw - ' in line: |
---|
3608 | if 'nan' in line: |
---|
3609 | Rwp = 100.0 |
---|
3610 | else: |
---|
3611 | Rwp = float(line.split(':')[1]) |
---|
3612 | results = resline.replace('(','').split(')')[:-1] |
---|
3613 | results = ['@'+result.lstrip() for result in results] |
---|
3614 | results = [item.split() for item in results] |
---|
3615 | RMCPdict['Parms'] = dict([[item[0][1:-1],float(item[1])] for item in results]) #{'n':val,...} |
---|
3616 | if RMCPdict['refinement'] == 'normal': |
---|
3617 | fName = General['Name']+'-PDFfit.py' |
---|
3618 | py = open(fName.replace(' ','_'),'r') |
---|
3619 | pylines = py.readlines() |
---|
3620 | py.close() |
---|
3621 | py = open(fName.replace(' ','_'),'w') |
---|
3622 | newpy = [] |
---|
3623 | for pyline in pylines: |
---|
3624 | if 'setpar' in pyline: |
---|
3625 | parm = pyline.split('(')[1].split(',')[0] |
---|
3626 | newpy.append('pf.setpar(%s,%.5f)\n'%(parm,RMCPdict['Parms'][parm])) |
---|
3627 | else: |
---|
3628 | newpy.append(pyline) |
---|
3629 | py.writelines(newpy) |
---|
3630 | py.close() |
---|
3631 | RMCPdict.update(XNdata) |
---|
3632 | results = dict([[item[0][:-1],float(item[1])] for item in results if item[0][:-1] in RMCPdict['AtomVar']]) |
---|
3633 | RMCPdict['AtomVar'].update(results) |
---|
3634 | return None |
---|
3635 | else: #sequential |
---|
3636 | newParms = dict([[item[0][1:-1],[float(item[1]),float(item[2])]] for item in results]) #{'n':[val,esd],...} |
---|
3637 | return newParms,Rwp |
---|
3638 | |
---|
3639 | def MakefullrmcSupercell(Phase,RMCPdict): |
---|
3640 | '''Create a fullrmc supercell from GSAS-II |
---|
3641 | |
---|
3642 | :param dict Phase: phase information from data tree |
---|
3643 | :param dict RMCPdict: fullrmc parameters from GUI |
---|
3644 | :param list grpDict: a list of lists where the inner list |
---|
3645 | contains the atom numbers contained in each group. e.g. |
---|
3646 | [[0,1,2,3,4],[5,6],[4,6]] creates three groups with |
---|
3647 | atoms 0-4 in the first |
---|
3648 | atoms 5 & 6 in the second and |
---|
3649 | atoms 4 & 6 in the third. Note that it is fine that |
---|
3650 | atom 4 appears in two groups. |
---|
3651 | ''' |
---|
3652 | #for i in (0,1): grpDict[i].append(1) # debug: 1st & 2nd atoms in 2nd group |
---|
3653 | cell = Phase['General']['Cell'][1:7] |
---|
3654 | A,B = G2lat.cell2AB(cell) |
---|
3655 | cx,ct,cs,cia = Phase['General']['AtomPtrs'] |
---|
3656 | SGData = Phase['General']['SGData'] |
---|
3657 | atomlist = [] |
---|
3658 | for i,atom in enumerate(Phase['Atoms']): |
---|
3659 | el = ''.join([i for i in atom[ct] if i.isalpha()]) |
---|
3660 | grps = [j for j,g in enumerate(RMCPdict.get('Groups',[])) if i in g] |
---|
3661 | atomlist.append((el, atom[ct-1], grps)) |
---|
3662 | # create a list of coordinates with symmetry & unit cell translation duplicates |
---|
3663 | coordlist = [] |
---|
3664 | cellnum = -1 |
---|
3665 | for a in range(int(0.5-RMCPdict['SuperCell'][0]/2),int(1+RMCPdict['SuperCell'][0]/2)): |
---|
3666 | for b in range(int(0.5-RMCPdict['SuperCell'][1]/2),int(1+RMCPdict['SuperCell'][1]/2)): |
---|
3667 | for c in range(int(0.5-RMCPdict['SuperCell'][2]/2),int(1+RMCPdict['SuperCell'][2]/2)): |
---|
3668 | cellnum += 1 |
---|
3669 | for i,atom in enumerate(Phase['Atoms']): |
---|
3670 | for item in G2spc.GenAtom(atom[cx:cx+3],SGData,Move=False): |
---|
3671 | # if i == 0: print(item[0]+[a,b,c]) |
---|
3672 | xyzOrth = np.inner(A,item[0]+[a,b,c]) |
---|
3673 | #coordlist.append((i,list(xyzOrth),cellnum,list(item[0]+[a,b,c]))) |
---|
3674 | coordlist.append((item[1],cellnum,i,list(xyzOrth))) |
---|
3675 | return atomlist,coordlist |
---|
3676 | |
---|
3677 | def MakefullrmcRun(pName,Phase,RMCPdict): |
---|
3678 | '''Creates a script to run fullrmc. Returns the name of the file that was |
---|
3679 | created. |
---|
3680 | ''' |
---|
3681 | BondList = {} |
---|
3682 | for k in RMCPdict['Pairs']: |
---|
3683 | if RMCPdict['Pairs'][k][1]+RMCPdict['Pairs'][k][2]>0: |
---|
3684 | BondList[k] = (RMCPdict['Pairs'][k][1],RMCPdict['Pairs'][k][2]) |
---|
3685 | AngleList = [] |
---|
3686 | for angle in RMCPdict['Angles']: |
---|
3687 | if angle[3] == angle[4] or angle[5] >= angle[6] or angle[6] <= 0: |
---|
3688 | continue |
---|
3689 | for i in (0,1,2): |
---|
3690 | angle[i] = angle[i].strip() |
---|
3691 | AngleList.append(angle) |
---|
3692 | # rmin = RMCPdict['min Contact'] |
---|
3693 | cell = Phase['General']['Cell'][1:7] |
---|
3694 | SymOpList = G2spc.AllOps(Phase['General']['SGData'])[0] |
---|
3695 | cx,ct,cs,cia = Phase['General']['AtomPtrs'] |
---|
3696 | atomsList = [] |
---|
3697 | for atom in Phase['Atoms']: |
---|
3698 | el = ''.join([i for i in atom[ct] if i.isalpha()]) |
---|
3699 | atomsList.append([el] + atom[cx:cx+4]) |
---|
3700 | projDir,projName = os.path.split(os.path.abspath(pName)) |
---|
3701 | scrname = pName+'-fullrmc.py' |
---|
3702 | restart = '%s_restart.pdb'%pName |
---|
3703 | Files = RMCPdict['files'] |
---|
3704 | rundata = '' |
---|
3705 | rundata += '## fullrmc %s file ##\n## OK to edit this by hand ##\n'%scrname |
---|
3706 | rundata += '# created in '+__file__+" v"+filversion.split()[1] |
---|
3707 | rundata += dt.datetime.strftime(dt.datetime.now()," at %Y-%m-%dT%H:%M\n") |
---|
3708 | rundata += ''' |
---|
3709 | # fullrmc imports (all that are potentially useful) |
---|
3710 | import os,glob |
---|
3711 | import time |
---|
3712 | import pickle |
---|
3713 | import types |
---|
3714 | import copy |
---|
3715 | import numpy as np |
---|
3716 | import matplotlib as mpl |
---|
3717 | import fullrmc |
---|
3718 | from pdbparser import pdbparser |
---|
3719 | from pdbparser.Utilities.Database import __ATOM__ |
---|
3720 | from fullrmc.Core import Collection |
---|
3721 | from fullrmc.Engine import Engine |
---|
3722 | import fullrmc.Constraints.PairDistributionConstraints as fPDF |
---|
3723 | from fullrmc.Constraints.StructureFactorConstraints import ReducedStructureFactorConstraint, StructureFactorConstraint |
---|
3724 | from fullrmc.Constraints.RadialDistributionConstraints import RadialDistributionConstraint |
---|
3725 | from fullrmc.Constraints.StructureFactorConstraints import NormalizedStructureFactorConstraint |
---|
3726 | from fullrmc.Constraints.DistanceConstraints import DistanceConstraint |
---|
3727 | from fullrmc.Constraints.BondConstraints import BondConstraint |
---|
3728 | from fullrmc.Constraints.AngleConstraints import BondsAngleConstraint |
---|
3729 | from fullrmc.Constraints.DihedralAngleConstraints import DihedralAngleConstraint |
---|
3730 | from fullrmc.Generators.Swaps import SwapPositionsGenerator |
---|
3731 | from fullrmc.Core.MoveGenerator import MoveGeneratorCollector |
---|
3732 | from fullrmc.Generators.Translations import TranslationGenerator |
---|
3733 | from fullrmc.Generators.Rotations import RotationGenerator |
---|
3734 | |
---|
3735 | # utility routines |
---|
3736 | def writeHeader(ENGINE,statFP): |
---|
3737 | """header for stats file""" |
---|
3738 | statFP.write('generated-steps, total-error, ') |
---|
3739 | for c in ENGINE.constraints: |
---|
3740 | statFP.write(c.constraintName) |
---|
3741 | statFP.write(', ') |
---|
3742 | statFP.write('\\n') |
---|
3743 | statFP.flush() |
---|
3744 | |
---|
3745 | def writeCurrentStatus(ENGINE,statFP,plotF): |
---|
3746 | """line in stats file & current constraint plots""" |
---|
3747 | statFP.write(str(ENGINE.generated)) |
---|
3748 | statFP.write(', ') |
---|
3749 | statFP.write(str(ENGINE.totalStandardError)) |
---|
3750 | statFP.write(', ') |
---|
3751 | for c in ENGINE.constraints: |
---|
3752 | statFP.write(str(c.standardError)) |
---|
3753 | statFP.write(', ') |
---|
3754 | statFP.write('\\n') |
---|
3755 | statFP.flush() |
---|
3756 | mpl.use('agg') |
---|
3757 | fp = open(plotF,'wb') |
---|
3758 | for c in ENGINE.constraints: |
---|
3759 | p = c.plot(show=False) |
---|
3760 | p[0].canvas.draw() |
---|
3761 | image = p[0].canvas.buffer_rgba() |
---|
3762 | pickle.dump(c.constraintName,fp) |
---|
3763 | pickle.dump(np.array(image),fp) |
---|
3764 | fp.close() |
---|
3765 | |
---|
3766 | def calcRmax(ENGINE): |
---|
3767 | """from Bachir, works for non-orthorhombic cells""" |
---|
3768 | a,b,c = ENGINE.basisVectors |
---|
3769 | lens = [] |
---|
3770 | ts = np.linalg.norm(np.cross(a,b))/2 |
---|
3771 | lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(b)] ) |
---|
3772 | ts = np.linalg.norm(np.cross(b,c))/2 |
---|
3773 | lens.extend( [ts/np.linalg.norm(b), ts/np.linalg.norm(c)] ) |
---|
3774 | ts = np.linalg.norm(np.cross(a,c))/2 |
---|
3775 | lens.extend( [ts/np.linalg.norm(a), ts/np.linalg.norm(c)] ) |
---|
3776 | return min(lens) |
---|
3777 | ''' |
---|
3778 | if RMCPdict.get('Groups',[]): rundata += ''' |
---|
3779 | def makepdb(atoms, coords, bbox=None): |
---|
3780 | """creates a supercell directly from atom info""" |
---|
3781 | # used when ENGINE.build_crystal_set_pdb is not called |
---|
3782 | prevcell = None |
---|
3783 | rec = copy.copy(__ATOM__) |
---|
3784 | rec['residue_name'] = 'MOL' |
---|
3785 | records = [] |
---|
3786 | seqNum = 0 |
---|
3787 | segId = '0' |
---|
3788 | groups = {} |
---|
3789 | for symcell in set([(sym,cell) for sym,cell,atm,xyz in coords]): |
---|
3790 | seqNum += 1 |
---|
3791 | if seqNum == 9999: |
---|
3792 | seqNum = 1 |
---|
3793 | segId = str(int(segId) + 1) |
---|
3794 | for i,(sym,cell,atm,(x,y,z)) in enumerate(coords): |
---|
3795 | if (sym,cell) != symcell: continue |
---|
3796 | rec = copy.copy(rec) |
---|
3797 | for grp in atoms[atm][2]: |
---|
3798 | if (sym,cell) not in groups: |
---|
3799 | groups[(sym,cell)] = {} |
---|
3800 | if grp not in groups[(sym,cell)]: |
---|
3801 | groups[(sym,cell)][grp] = [len(records)] |
---|
3802 | else: |
---|
3803 | groups[(sym,cell)][grp].append(len(records)) |
---|
3804 | rec['coordinates_x'] = x |
---|
3805 | rec['coordinates_y'] = y |
---|
3806 | rec['coordinates_z'] = z |
---|
3807 | rec['element_symbol'] = atoms[atm][0] |
---|
3808 | rec['atom_name'] = atoms[atm][1] |
---|
3809 | rec['sequence_number'] = seqNum |
---|
3810 | rec['segment_identifier'] = segId |
---|
3811 | records.append(rec) |
---|
3812 | # create pdb |
---|
3813 | pdb = pdbparser() |
---|
3814 | pdb.records = records |
---|
3815 | if groups: |
---|
3816 | return pdb,[groups[j][i] for j in groups for i in groups[j]] |
---|
3817 | else: |
---|
3818 | return pdb,[] |
---|
3819 | ''' |
---|
3820 | rundata += ''' |
---|
3821 | ### When True, erases an existing engine to provide a fresh start |
---|
3822 | FRESH_START = {:} |
---|
3823 | dirName = "{:}" |
---|
3824 | prefix = "{:}" |
---|
3825 | project = prefix + "-fullrmc" |
---|
3826 | time0 = time.time() |
---|
3827 | '''.format(RMCPdict['ReStart'][0],projDir,projName) |
---|
3828 | |
---|
3829 | rundata += '# setup structure\n' |
---|
3830 | rundata += 'cell = ' + str(cell) + '\n' |
---|
3831 | rundata += 'supercell = ' + str(RMCPdict['SuperCell']) + '\n' |
---|
3832 | rundata += '\n# define structure info\n' |
---|
3833 | if RMCPdict.get('Groups',[]): |
---|
3834 | # compute bounding box coordinates |
---|
3835 | bbox = [] |
---|
3836 | A,B = G2lat.cell2AB(cell) |
---|
3837 | for i in range(3): |
---|
3838 | for val in int(0.5-RMCPdict['SuperCell'][i]/2),int(1+RMCPdict['SuperCell'][0]/2): |
---|
3839 | fpos = [0,0,0] |
---|
3840 | fpos[i] = val |
---|
3841 | bbox.append(np.inner(A,fpos)) |
---|
3842 | rundata += 'bboxlist = [ # orthogonal coordinate for supercell corners\n' |
---|
3843 | for i in bbox: |
---|
3844 | rundata += ' '+str(list(i))+',\n' |
---|
3845 | rundata += ' ] # bboxlist\n\n' |
---|
3846 | atomlist,coordlist = MakefullrmcSupercell(Phase,RMCPdict) |
---|
3847 | rundata += 'atomlist = [ # [element, label, grouplist]\n' |
---|
3848 | for i in atomlist: |
---|
3849 | rundata += ' '+str(i)+',\n' |
---|
3850 | rundata += ' ] # atomlist\n\n' |
---|
3851 | rundata += 'coordlist = [ # (sym#, cell#, atom#, [ortho coords],)\n' |
---|
3852 | for i in coordlist: |
---|
3853 | rundata += ' '+str(i)+',\n' |
---|
3854 | rundata += ' ] # coordlist\n' |
---|
3855 | else: |
---|
3856 | rundata += "SymOpList = "+str([i.lower() for i in SymOpList]) + '\n' |
---|
3857 | rundata += 'atomList = ' + str(atomsList).replace('],','],\n ') + '\n' |
---|
3858 | |
---|
3859 | rundata += '\n# initialize engine\n' |
---|
3860 | rundata += ''' |
---|
3861 | engineFileName = os.path.join(dirName, project + '.rmc') |
---|
3862 | projectStats = os.path.join(dirName, project + '.stats') |
---|
3863 | projectPlots = os.path.join(dirName, project + '.plots') |
---|
3864 | projectXYZ = os.path.join(dirName, project + '.atoms') |
---|
3865 | pdbFile = os.path.join(dirName, project + '_restart.pdb') |
---|
3866 | # check Engine exists if so (and not FRESH_START) load it otherwise build it |
---|
3867 | ENGINE = Engine(path=None) |
---|
3868 | if not ENGINE.is_engine(engineFileName) or FRESH_START: |
---|
3869 | ENGINE = Engine(path=engineFileName, freshStart=True) |
---|
3870 | ''' |
---|
3871 | if RMCPdict.get('Groups',[]): |
---|
3872 | rundata += ''' |
---|
3873 | # create structure from GSAS-II constructed supercell |
---|
3874 | bbox = (np.array(bboxlist[1::2])-np.array(bboxlist[0::2])).flatten() |
---|
3875 | pdb,grouplist = makepdb(atomlist,coordlist,bbox) |
---|
3876 | ENGINE.set_pdb(pdb) |
---|
3877 | ENGINE.set_boundary_conditions(bbox) |
---|
3878 | if grouplist: ENGINE.set_groups(grouplist) |
---|
3879 | ''' |
---|
3880 | if RMCPdict.get('GroupMode',0) == 0: # 'Rotate & Translate' |
---|
3881 | rundata += ''' |
---|
3882 | for g in ENGINE.groups: |
---|
3883 | TMG = TranslationGenerator(amplitude=0.2) # create translation generator |
---|
3884 | if len(g) > 1: # create rotation generator for groups with more than 1 atom |
---|
3885 | RMG = RotationGenerator(amplitude=2) |
---|
3886 | MG = MoveGeneratorCollector(collection=[TMG,RMG],randomize=True) |
---|
3887 | else: |
---|
3888 | MG = MoveGeneratorCollector(collection=[TMG],randomize=True) |
---|
3889 | g.set_move_generator( MG ) |
---|
3890 | ''' |
---|
3891 | elif RMCPdict.get('GroupMode',0) == 1: # 'Rotate only' |
---|
3892 | rundata += ''' |
---|
3893 | for g in ENGINE.groups: |
---|
3894 | if len(g) > 1: # create rotation generator for groups with more than 1 atom |
---|
3895 | RMG = RotationGenerator(amplitude=2) |
---|
3896 | g.set_move_generator( RMG ) |
---|
3897 | ''' |
---|
3898 | else: # 'Translate only' |
---|
3899 | rundata += ' # translate only set by default' |
---|
3900 | else: |
---|
3901 | rundata += ''' |
---|
3902 | # create structure, let fullrmc construct supercell |
---|
3903 | ENGINE.build_crystal_set_pdb(symOps = SymOpList, |
---|
3904 | atoms = atomList, |
---|
3905 | unitcellBC = cell, |
---|
3906 | supercell = supercell) |
---|
3907 | ENGINE.set_groups_as_atoms() |
---|
3908 | ''' |
---|
3909 | rundata += ' rho0 = len(ENGINE.allNames)/ENGINE.volume\n' |
---|
3910 | rundata += '\n # "Constraints" (includes experimental data) setup\n' |
---|
3911 | # settings that require a new Engine |
---|
3912 | for File in Files: |
---|
3913 | filDat = RMCPdict['files'][File] |
---|
3914 | if not os.path.exists(filDat[0]): continue |
---|
3915 | sfwt = 'neutronCohb' |
---|
3916 | if 'Xray' in File: |
---|
3917 | sfwt = 'atomicNumber' |
---|
3918 | if 'G(r)' in File: |
---|
3919 | rundata += ' GR = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0] |
---|
3920 | if filDat[3] == 0: |
---|
3921 | #rundata += ''' # read and xform G(r) as defined in RMCProfile |
---|
3922 | # see eq. 44 in Keen, J. Appl. Cryst. (2001) 34 172-177\n''' |
---|
3923 | #rundata += ' GR[1] *= 4 * np.pi * GR[0] * rho0 / sumCiBi2\n' |
---|
3924 | #rundata += ' GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt |
---|
3925 | rundata += ' # G(r) as defined in RMCProfile\n' |
---|
3926 | rundata += ' GofR = RadialDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt |
---|
3927 | elif filDat[3] == 1: |
---|
3928 | rundata += ' # This is G(r) as defined in PDFFIT\n' |
---|
3929 | rundata += ' GofR = fPDF.PairDistributionConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt |
---|
3930 | elif filDat[3] == 2: |
---|
3931 | rundata += ' # This is g(r)\n' |
---|
3932 | rundata += ' GofR = fPDF.PairCorrelationConstraint(experimentalData=GR.T, weighting="%s")\n'%sfwt |
---|
3933 | else: |
---|
3934 | raise ValueError('Invalid G(r) type: '+str(filDat[3])) |
---|
3935 | rundata += ' ENGINE.add_constraints([GofR])\n' |
---|
3936 | rundata += ' GofR.set_limits((None, calcRmax(ENGINE)))\n' |
---|
3937 | if RMCPdict['addThermalBroadening']: |
---|
3938 | rundata += " GofR.set_thermal_corrections({'defaultFactor': 0.001})\n" |
---|
3939 | rundata += " GofR.thermalCorrections['factors'] = {\n" |
---|
3940 | RMCPdict['addThermalBroadening'] |
---|
3941 | for atm1 in RMCPdict['aTypes']: |
---|
3942 | for atm2 in RMCPdict['aTypes']: |
---|
3943 | rundata += " ('{}', '{}'): {},\n".format( |
---|
3944 | atm1,atm2, |
---|
3945 | (RMCPdict['ThermalU'].get(atm1,0.005)+ RMCPdict['ThermalU'].get(atm2,0.005))/2) |
---|
3946 | rundata += ' }\n' |
---|
3947 | elif '(Q)' in File: |
---|
3948 | rundata += ' SOQ = np.loadtxt(os.path.join(dirName,"%s")).T\n'%filDat[0] |
---|
3949 | if filDat[3] == 0: |
---|
3950 | rundata += ' # F(Q) as defined in RMCProfile\n' |
---|
3951 | #rundata += ' SOQ[1] *= 1 / sumCiBi2\n' |
---|
3952 | if filDat[4]: |
---|
3953 | rundata += ' SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n' |
---|
3954 | rundata += ' SofQ = NormalizedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt |
---|
3955 | elif filDat[3] == 1: |
---|
3956 | rundata += ' # S(Q) as defined in PDFFIT\n' |
---|
3957 | rundata += ' SOQ[1] -= 1\n' |
---|
3958 | if filDat[4]: |
---|
3959 | rundata += ' SOQ[1] = Collection.sinc_convolution(q=SOQ[0],sq=SOQ[1],rmax=calcRmax(ENGINE))\n' |
---|
3960 | rundata += ' SofQ = ReducedStructureFactorConstraint(experimentalData=SOQ.T, weighting="%s")\n'%sfwt |
---|
3961 | else: |
---|
3962 | raise ValueError('Invalid S(Q) type: '+str(filDat[3])) |
---|
3963 | rundata += ' ENGINE.add_constraints([SofQ])\n' |
---|
3964 | else: |
---|
3965 | print('What is this?') |
---|
3966 | minDists = '' |
---|
3967 | if BondList and RMCPdict.get('useBondConstraints',True): |
---|
3968 | rundata += ''' B_CONSTRAINT = BondConstraint() |
---|
3969 | ENGINE.add_constraints(B_CONSTRAINT) |
---|
3970 | B_CONSTRAINT.create_supercell_bonds(bondsDefinition=[ |
---|
3971 | ''' |
---|
3972 | for pair in BondList: |
---|
3973 | e1,e2 = pair.split('-') |
---|
3974 | d1,d2 = BondList[pair] |
---|
3975 | if d1 == 0: continue |
---|
3976 | if d2 == 0: |
---|
3977 | minDists += '("element","{}","{}",{}),'.format(e1.strip(),e2.strip(),d1) |
---|
3978 | else: |
---|
3979 | rundata += ' ("element","{}","{}",{},{}),\n'.format( |
---|
3980 | e1.strip(),e2.strip(),d1,d2) |
---|
3981 | rundata += ' ])\n' |
---|
3982 | rundata += ' D_CONSTRAINT = DistanceConstraint(defaultLowerDistance={})\n'.format(RMCPdict['min Contact']) |
---|
3983 | if minDists: |
---|
3984 | rundata += " D_CONSTRAINT.set_pairs_definition( {'inter':[" + minDists + "]})\n" |
---|
3985 | rundata += ' ENGINE.add_constraints(D_CONSTRAINT)\n' |
---|
3986 | |
---|
3987 | if AngleList: |
---|
3988 | rundata += ''' A_CONSTRAINT = BondsAngleConstraint() |
---|
3989 | ENGINE.add_constraints(A_CONSTRAINT) |
---|
3990 | A_CONSTRAINT.create_supercell_angles(anglesDefinition=[ |
---|
3991 | ''' |
---|
3992 | for item in AngleList: |
---|
3993 | rundata += (' '+ |
---|
3994 | '("element","{1}","{0}","{2}",{5},{6},{5},{6},{3},{4}),\n'.format(*item)) |
---|
3995 | rundata += ' ])\n' |
---|
3996 | rundata += ''' |
---|
3997 | for f in glob.glob(os.path.join(dirName,prefix+"_*.log")): os.remove(f) |
---|
3998 | ENGINE.save() |
---|
3999 | else: |
---|
4000 | ENGINE = ENGINE.load(path=engineFileName) |
---|
4001 | |
---|
4002 | ENGINE.set_log_file(os.path.join(dirName,prefix)) |
---|
4003 | ''' |
---|
4004 | if RMCPdict['Swaps']: |
---|
4005 | rundata += '\n#set up for site swaps\n' |
---|
4006 | rundata += 'aN = ENGINE.allNames\n' |
---|
4007 | rundata += 'SwapGen = {}\n' |
---|
4008 | for swap in RMCPdict['Swaps']: |
---|
4009 | rundata += 'SwapA = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[0] |
---|
4010 | rundata += 'SwapB = [[idx] for idx in range(len(aN)) if aN[idx]=="%s"]\n'%swap[1] |
---|
4011 | rundata += 'SwapGen["%s-%s"] = [SwapPositionsGenerator(swapList=SwapA),SwapPositionsGenerator(swapList=SwapB),%.2f]\n'%(swap[0],swap[1],swap[2]) |
---|
4012 | rundata += ' for swaps in SwapGen:\n' |
---|
4013 | rundata += ' AB = swaps.split("-")\n' |
---|
4014 | rundata += ' ENGINE.set_groups_as_atoms()\n' |
---|
4015 | rundata += ' for g in ENGINE.groups:\n' |
---|
4016 | rundata += ' if aN[g.indexes[0]]==AB[0]:\n' |
---|
4017 | rundata += ' g.set_move_generator(SwapGen[swaps][0])\n' |
---|
4018 | rundata += ' elif aN[g.indexes[0]]==AB[1]:\n' |
---|
4019 | rundata += ' g.set_move_generator(SwapGen[swaps][1])\n' |
---|
4020 | rundata += ' sProb = SwapGen[swaps][2]\n' |
---|
4021 | rundata += '''for c in ENGINE.constraints: |
---|
4022 | if hasattr(c, '_ExperimentalConstraint__adjustScaleFactor'): |
---|
4023 | def _constraint_copy_needs_lut(self, *args, **kwargs): |
---|
4024 | result = super(self.__class__, self)._constraint_copy_needs_lut(*args, **kwargs) |
---|
4025 | result['_ExperimentalConstraint__adjustScaleFactor'] = '_ExperimentalConstraint__adjustScaleFactor' |
---|
4026 | return result |
---|
4027 | c._constraint_copy_needs_lut = types.MethodType(_constraint_copy_needs_lut, c) |
---|
4028 | ''' |
---|
4029 | # rundata += '\n# set weights -- do this now so values can be changed without a restart\n' |
---|
4030 | # rundata += 'wtDict = {}\n' |
---|
4031 | # for File in Files: |
---|
4032 | # filDat = RMCPdict['files'][File] |
---|
4033 | # if not os.path.exists(filDat[0]): continue |
---|
4034 | # if 'Xray' in File: |
---|
4035 | # sfwt = 'atomicNumber' |
---|
4036 | # else: |
---|
4037 | # sfwt = 'neutronCohb' |
---|
4038 | # if 'G(r)' in File: |
---|
4039 | # typ = 'Pair' |
---|
4040 | # elif '(Q)' in File: |
---|
4041 | # typ = 'Struct' |
---|
4042 | # rundata += 'wtDict["{}-{}"] = {}\n'.format(typ,sfwt,filDat[1]) |
---|
4043 | rundata += '\n# set PDF fitting range\n' |
---|
4044 | rundata += 'for c in ENGINE.constraints: # loop over predefined constraints\n' |
---|
4045 | rundata += ' if type(c) is fPDF.PairDistributionConstraint:\n' |
---|
4046 | # rundata += ' c.set_variance_squared(1./wtDict["Pair-"+c.weighting])\n' |
---|
4047 | rundata += ' c.set_limits((None,calcRmax(ENGINE)))\n' |
---|
4048 | if RMCPdict['FitScale']: |
---|
4049 | rundata += ' c.set_adjust_scale_factor((10, 0.01, 100.))\n' |
---|
4050 | # rundata += ' c.set_variance_squared(1./wtDict["Struct-"+c.weighting])\n' |
---|
4051 | if RMCPdict['FitScale']: |
---|
4052 | rundata += ' elif type(c) is ReducedStructureFactorConstraint:\n' |
---|
4053 | rundata += ' c.set_adjust_scale_factor((10, 0.01, 100.))\n' |
---|
4054 | # torsions difficult to implement, must be internal to cell & named with |
---|
4055 | # fullrmc atom names |
---|
4056 | # if len(RMCPdict['Torsions']): # Torsions currently commented out in GUI |
---|
4057 | # rundata += 'for c in ENGINE.constraints: # look for Dihedral Angle Constraints\n' |
---|
4058 | # rundata += ' if type(c) is DihedralAngleConstraint:\n' |
---|
4059 | # rundata += ' c.set_variance_squared(%f)\n'%RMCPdict['Torsion Weight'] |
---|
4060 | # rundata += ' c.create_angles_by_definition(anglesDefinition={"%s":[\n'%Res |
---|
4061 | # for torsion in RMCPdict['Torsions']: |
---|
4062 | # rundata += ' %s\n'%str(tuple(torsion)) |
---|
4063 | # rundata += ' ]})\n' |
---|
4064 | rundata += ''' |
---|
4065 | if FRESH_START: |
---|
4066 | # initialize engine with one step to get starting config energetics |
---|
4067 | ENGINE.run(restartPdb=pdbFile,numberOfSteps=1, saveFrequency=1) |
---|
4068 | statFP = open(projectStats,'w') |
---|
4069 | writeHeader(ENGINE,statFP) |
---|
4070 | writeCurrentStatus(ENGINE,statFP,projectPlots) |
---|
4071 | else: |
---|
4072 | statFP = open(projectStats,'a') |
---|
4073 | |
---|
4074 | # setup runs for fullrmc |
---|
4075 | ''' |
---|
4076 | rundata += 'steps = {}\n'.format(RMCPdict['Steps/cycle']) |
---|
4077 | rundata += 'for _ in range({}):\n'.format(RMCPdict['Cycles']) |
---|
4078 | rundata += ' expected = ENGINE.generated+steps\n' |
---|
4079 | |
---|
4080 | rundata += ' ENGINE.run(restartPdb=pdbFile,numberOfSteps=steps, saveFrequency=steps)\n' |
---|
4081 | rundata += ' writeCurrentStatus(ENGINE,statFP,projectPlots)\n' |
---|
4082 | rundata += ' if ENGINE.generated != expected: break # run was stopped' |
---|
4083 | rundata += ''' |
---|
4084 | statFP.close() |
---|
4085 | fp = open(projectXYZ,'w') # save final atom positions |
---|
4086 | fp.write('cell: {} {} {} {} {} {}\\n') |
---|
4087 | fp.write('supercell: {} {} {}\\n') |
---|
4088 | '''.format(*cell,*RMCPdict['SuperCell']) |
---|
4089 | rundata += '''# loop over atoms |
---|
4090 | for n,e,(x,y,z) in zip(ENGINE.allNames, |
---|
4091 | ENGINE.allElements,ENGINE.realCoordinates): |
---|
4092 | fp.write('{} {} {:.5f} {:.5f} {:.5f}\\n'.format(n,e,x,y,z)) |
---|
4093 | fp.close() |
---|
4094 | print("ENGINE run time %.2f s"%(time.time()-time0)) |
---|
4095 | ''' |
---|
4096 | rfile = open(scrname,'w') |
---|
4097 | rfile.writelines(rundata) |
---|
4098 | rfile.close() |
---|
4099 | return scrname |
---|
4100 | |
---|
4101 | def GetRMCBonds(general,RMCPdict,Atoms,bondList): |
---|
4102 | bondDist = [] |
---|
4103 | Cell = general['Cell'][1:7] |
---|
4104 | Supercell = RMCPdict['SuperCell'] |
---|
4105 | Trans = np.eye(3)*np.array(Supercell) |
---|
4106 | Cell = G2lat.TransformCell(Cell,Trans)[:6] |
---|
4107 | Amat,Bmat = G2lat.cell2AB(Cell) |
---|
4108 | indices = (-1,0,1) |
---|
4109 | Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) |
---|
4110 | for bonds in bondList: |
---|
4111 | Oxyz = np.array(Atoms[bonds[0]][1:]) |
---|
4112 | Txyz = np.array([Atoms[tgt-1][1:] for tgt in bonds[1]]) |
---|
4113 | Dx = np.array([Txyz-Oxyz+unit for unit in Units]) |
---|
4114 | Dx = np.sqrt(np.sum(np.inner(Dx,Amat)**2,axis=2)) |
---|
4115 | for dx in Dx.T: |
---|
4116 | bondDist.append(np.min(dx)) |
---|
4117 | return np.array(bondDist) |
---|
4118 | |
---|
4119 | def GetRMCAngles(general,RMCPdict,Atoms,angleList): |
---|
4120 | bondAngles = [] |
---|
4121 | Cell = general['Cell'][1:7] |
---|
4122 | Supercell = RMCPdict['SuperCell'] |
---|
4123 | Trans = np.eye(3)*np.array(Supercell) |
---|
4124 | Cell = G2lat.TransformCell(Cell,Trans)[:6] |
---|
4125 | Amat,Bmat = G2lat.cell2AB(Cell) |
---|
4126 | indices = (-1,0,1) |
---|
4127 | Units = np.array([[h,k,l] for h in indices for k in indices for l in indices]) |
---|
4128 | for angle in angleList: |
---|
4129 | Oxyz = np.array(Atoms[angle[0]][1:]) |
---|
4130 | TAxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[0]]) |
---|
4131 | TBxyz = np.array([Atoms[tgt-1][1:] for tgt in angle[1].T[1]]) |
---|
4132 | DAxV = np.inner(np.array([TAxyz-Oxyz+unit for unit in Units]),Amat) |
---|
4133 | DAx = np.sqrt(np.sum(DAxV**2,axis=2)) |
---|
4134 | DBxV = np.inner(np.array([TBxyz-Oxyz+unit for unit in Units]),Amat) |
---|
4135 | DBx = np.sqrt(np.sum(DBxV**2,axis=2)) |
---|
4136 | iDAx = np.argmin(DAx,axis=0) |
---|
4137 | iDBx = np.argmin(DBx,axis=0) |
---|
4138 | for i,[iA,iB] in enumerate(zip(iDAx,iDBx)): |
---|
4139 | DAv = DAxV[iA,i]/DAx[iA,i] |
---|
4140 | DBv = DBxV[iB,i]/DBx[iB,i] |
---|
4141 | bondAngles.append(npacosd(np.sum(DAv*DBv))) |
---|
4142 | return np.array(bondAngles) |
---|
4143 | |
---|
4144 | def ISO2PDFfit(Phase): |
---|
4145 | ''' Creates new phase structure to be used for PDFfit from an ISODISTORT mode displacement phase. |
---|
4146 | It builds the distortion mode parameters to be used as PDFfit variables for atom displacements from |
---|
4147 | the original parent positions as transformed to the child cell wiht symmetry defined from ISODISTORT. |
---|
4148 | |
---|
4149 | :param Phase: dict GSAS-II Phase structure; must contain ISODISTORT dict. NB: not accessed otherwise |
---|
4150 | |
---|
4151 | :returns: dict: GSAS-II Phase structure; will contain ['RMC']['PDFfit'] dict |
---|
4152 | ''' |
---|
4153 | |
---|
4154 | Trans = np.eye(3) |
---|
4155 | Uvec = np.zeros(3) |
---|
4156 | Vvec = np.zeros(3) |
---|
4157 | Phase = copy.deepcopy(Phase) |
---|
4158 | Atoms = Phase['Atoms'] |
---|
4159 | parentXYZ = Phase['ISODISTORT']['G2parentCoords'] #starting point for mode displacements |
---|
4160 | cx,ct,cs,cia = Phase['General']['AtomPtrs'] |
---|
4161 | for iat,atom in enumerate(Atoms): |
---|
4162 | atom[cx:cx+3] = parentXYZ[iat] |
---|
4163 | SGData = copy.deepcopy(Phase['General']['SGData']) |
---|
4164 | SGOps = SGData['SGOps'] |
---|
4165 | newPhase = copy.deepcopy(Phase) |
---|
4166 | newPhase['ranId'] = rand.randint(0,sys.maxsize) |
---|
4167 | newPhase['General']['Name'] += '_PDFfit' |
---|
4168 | newPhase['General']['SGData'] = G2spc.SpcGroup('P 1')[1] #this is for filled unit cell |
---|
4169 | newPhase,atCodes = G2lat.TransformPhase(Phase,newPhase,Trans,Uvec,Vvec,False) |
---|
4170 | newPhase['Histograms'] = {} |
---|
4171 | newPhase['Drawing'] = [] |
---|
4172 | Atoms = newPhase['Atoms'] |
---|
4173 | RMCPdict = newPhase['RMC']['PDFfit'] |
---|
4174 | ISOdict = newPhase['ISODISTORT'] |
---|
4175 | RMCPdict['AtomConstr'] = [] |
---|
4176 | RMCPdict['SGData'] = copy.deepcopy(SGData) #this is from the ISODISTORT child; defines PDFfit constraints |
---|
4177 | Norms = ISOdict['NormList'] |
---|
4178 | ModeMatrix = ISOdict['Mode2VarMatrix'] |
---|
4179 | RMCPdict['AtomVar'] = {'@%d'%(itm+21):val for itm,val in enumerate(ISOdict['modeDispl'])} |
---|
4180 | for iatm,[atom,atcode] in enumerate(zip(Atoms,atCodes)): |
---|
4181 | pid,opid = [int(item) for item in atcode.split(':')] |
---|
4182 | atmConstr = [atom[ct-1],atom[ct],'','','','','',atcode] |
---|
4183 | for ip,pname in enumerate(['%s_d%s'%(atom[ct-1],x) for x in ['x','y','z']]): |
---|
4184 | try: |
---|
4185 | conStr = '' |
---|
4186 | if Atoms[iatm][cx+ip]: |
---|
4187 | conStr += '%.5f'%Atoms[iatm][cx+ip] |
---|
4188 | pid = ISOdict['IsoVarList'].index(pname) |
---|
4189 | consVec = ModeMatrix[pid] |
---|
4190 | for ic,citm in enumerate(consVec): #NB: this assumes orthorhombic or lower symmetry |
---|
4191 | if opid < 0: |
---|
4192 | citm *= -SGOps[100-opid%100-1][0][ip][ip] #remove centering, if any |
---|
4193 | else: |
---|
4194 | citm *= SGOps[opid%100-1][0][ip][ip] |
---|
4195 | if citm > 0.: |
---|
4196 | conStr += '+%.5f*@%d'%(citm*Norms[ic],ic+21) |
---|
4197 | elif citm < 0.: |
---|
4198 | conStr += '%.5f*@%d'%(citm*Norms[ic],ic+21) |
---|
4199 | atmConstr[ip+2] = conStr |
---|
4200 | except ValueError: |
---|
4201 | atmConstr[ip+2] = '' |
---|
4202 | RMCPdict['AtomConstr'].append(atmConstr) |
---|
4203 | return newPhase |
---|
4204 | |
---|
4205 | def GetAtmDispList(ISOdata): |
---|
4206 | atmDispList = [] |
---|
4207 | MT = ISOdata['Mode2VarMatrix'].T |
---|
4208 | DispList = ISOdata['IsoVarList'] |
---|
4209 | N = len(DispList) |
---|
4210 | for I in range(N): |
---|
4211 | vary = [] |
---|
4212 | for i in range(N): |
---|
4213 | if MT[I,i]: |
---|
4214 | vary.append(DispList[i]) |
---|
4215 | atmDispList.append(vary) |
---|
4216 | return atmDispList |
---|
4217 | |
---|
4218 | #### Reflectometry calculations ################################################################################ |
---|
4219 | def REFDRefine(Profile,ProfDict,Inst,Limits,Substances,data): |
---|
4220 | G2fil.G2Print ('fit REFD data by '+data['Minimizer']+' using %.2f%% data resolution'%(data['Resolution'][0])) |
---|
4221 | |
---|
4222 | class RandomDisplacementBounds(object): |
---|
4223 | """random displacement with bounds""" |
---|
4224 | def __init__(self, xmin, xmax, stepsize=0.5): |
---|
4225 | self.xmin = xmin |
---|
4226 | self.xmax = xmax |
---|
4227 | self.stepsize = stepsize |
---|
4228 | |
---|
4229 | def __call__(self, x): |
---|
4230 | """take a random step but ensure the new position is within the bounds""" |
---|
4231 | while True: |
---|
4232 | # this could be done in a much more clever way, but it will work for example purposes |
---|
4233 | steps = self.xmax-self.xmin |
---|
4234 | xnew = x + np.random.uniform(-self.stepsize*steps, self.stepsize*steps, np.shape(x)) |
---|
4235 | if np.all(xnew < self.xmax) and np.all(xnew > self.xmin): |
---|
4236 | break |
---|
4237 | return xnew |
---|
4238 | |
---|
4239 | def GetModelParms(): |
---|
4240 | parmDict = {} |
---|
4241 | varyList = [] |
---|
4242 | values = [] |
---|
4243 | bounds = [] |
---|
4244 | parmDict['dQ type'] = data['dQ type'] |
---|
4245 | parmDict['Res'] = data['Resolution'][0]/(100.*sateln2) #% FWHM-->decimal sig |
---|
4246 | for parm in ['Scale','FltBack']: |
---|
4247 | parmDict[parm] = data[parm][0] |
---|
4248 | if data[parm][1]: |
---|
4249 | varyList.append(parm) |
---|
4250 | values.append(data[parm][0]) |
---|
4251 | bounds.append(Bounds[parm]) |
---|
4252 | parmDict['Layer Seq'] = np.array(['0',]+data['Layer Seq'].split()+[str(len(data['Layers'])-1),],dtype=int) |
---|
4253 | parmDict['nLayers'] = len(parmDict['Layer Seq']) |
---|
4254 | for ilay,layer in enumerate(data['Layers']): |
---|
4255 | name = layer['Name'] |
---|
4256 | cid = str(ilay)+';' |
---|
4257 | parmDict[cid+'Name'] = name |
---|
4258 | for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']: |
---|
4259 | parmDict[cid+parm] = layer.get(parm,[0.,False])[0] |
---|
4260 | if layer.get(parm,[0.,False])[1]: |
---|
4261 | varyList.append(cid+parm) |
---|
4262 | value = layer[parm][0] |
---|
4263 | values.append(value) |
---|
4264 | if value: |
---|
4265 | bound = [value*Bfac,value/Bfac] |
---|
4266 | else: |
---|
4267 | bound = [0.,10.] |
---|
4268 | bounds.append(bound) |
---|
4269 | if name not in ['vacuum','unit scatter']: |
---|
4270 | parmDict[cid+'rho'] = Substances[name]['Scatt density'] |
---|
4271 | parmDict[cid+'irho'] = Substances[name].get('XImag density',0.) |
---|
4272 | return parmDict,varyList,values,bounds |
---|
4273 | |
---|
4274 | def SetModelParms(): |
---|
4275 | line = ' Refined parameters: Histogram scale: %.4g'%(parmDict['Scale']) |
---|
4276 | if 'Scale' in varyList: |
---|
4277 | data['Scale'][0] = parmDict['Scale'] |
---|
4278 | line += ' esd: %.4g'%(sigDict['Scale']) |
---|
4279 | G2fil.G2Print (line) |
---|
4280 | line = ' Flat background: %15.4g'%(parmDict['FltBack']) |
---|
4281 | if 'FltBack' in varyList: |
---|
4282 | data['FltBack'][0] = parmDict['FltBack'] |
---|
4283 | line += ' esd: %15.3g'%(sigDict['FltBack']) |
---|
4284 | G2fil.G2Print (line) |
---|
4285 | for ilay,layer in enumerate(data['Layers']): |
---|
4286 | name = layer['Name'] |
---|
4287 | G2fil.G2Print (' Parameters for layer: %d %s'%(ilay,name)) |
---|
4288 | cid = str(ilay)+';' |
---|
4289 | line = ' ' |
---|
4290 | line2 = ' Scattering density: Real %.5g'%(Substances[name]['Scatt density']*parmDict[cid+'DenMul']) |
---|
4291 | line2 += ' Imag %.5g'%(Substances[name].get('XImag density',0.)*parmDict[cid+'DenMul']) |
---|
4292 | for parm in ['Thick','Rough','DenMul','Mag SLD','iDenMul']: |
---|
4293 | if parm in layer: |
---|
4294 | if parm == 'Rough': |
---|
4295 | layer[parm][0] = abs(parmDict[cid+parm]) #make positive |
---|
4296 | else: |
---|
4297 | layer[parm][0] = parmDict[cid+parm] |
---|
4298 | line += ' %s: %.3f'%(parm,layer[parm][0]) |
---|
4299 | if cid+parm in varyList: |
---|
4300 | line += ' esd: %.3g'%(sigDict[cid+parm]) |
---|
4301 | G2fil.G2Print (line) |
---|
4302 | G2fil.G2Print (line2) |
---|
4303 | |
---|
4304 | def calcREFD(values,Q,Io,wt,Qsig,parmDict,varyList): |
---|
4305 | parmDict.update(zip(varyList,values)) |
---|
4306 | M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io) |
---|
4307 | return M |
---|
4308 | |
---|
4309 | def sumREFD(values,Q,Io,wt,Qsig,parmDict,varyList): |
---|
4310 | parmDict.update(zip(varyList,values)) |
---|
4311 | M = np.sqrt(wt)*(getREFD(Q,Qsig,parmDict)-Io) |
---|
4312 | return np.sum(M**2) |
---|
4313 | |
---|
4314 | def getREFD(Q,Qsig,parmDict): |
---|
4315 | Ic = np.ones_like(Q)*parmDict['FltBack'] |
---|
4316 | Scale = parmDict['Scale'] |
---|
4317 | Nlayers = parmDict['nLayers'] |
---|
4318 | Res = parmDict['Res'] |
---|
4319 | depth = np.zeros(Nlayers) |
---|
4320 | rho = np.zeros(Nlayers) |
---|
4321 | irho = np.zeros(Nlayers) |
---|
4322 | sigma = np.zeros(Nlayers) |
---|
4323 | for ilay,lay in enumerate(parmDict['Layer Seq']): |
---|
4324 | cid = str(lay)+';' |
---|
4325 | depth[ilay] = parmDict[cid+'Thick'] |
---|
4326 | sigma[ilay] = parmDict[cid+'Rough'] |
---|
4327 | if parmDict[cid+'Name'] == u'unit scatter': |
---|
4328 | rho[ilay] = parmDict[cid+'DenMul'] |
---|
4329 | irho[ilay] = parmDict[cid+'iDenMul'] |
---|
4330 | elif 'vacuum' != parmDict[cid+'Name']: |
---|
4331 | rho[ilay] = parmDict[cid+'rho']*parmDict[cid+'DenMul'] |
---|
4332 | irho[ilay] = parmDict[cid+'irho']*parmDict[cid+'DenMul'] |
---|
4333 | if cid+'Mag SLD' in parmDict: |
---|
4334 | rho[ilay] += parmDict[cid+'Mag SLD'] |
---|
4335 | if parmDict['dQ type'] == 'None': |
---|
4336 | AB = abeles(0.5*Q,depth,rho,irho,sigma[1:]) #Q --> k, offset roughness for abeles |
---|
4337 | elif 'const' in parmDict['dQ type']: |
---|
4338 | AB = SmearAbeles(0.5*Q,Q*Res,depth,rho,irho,sigma[1:]) |
---|
4339 | else: #dQ/Q in data |
---|
4340 | AB = SmearAbeles(0.5*Q,Qsig,depth,rho,irho,sigma[1:]) |
---|
4341 | Ic += AB*Scale |
---|
4342 | return Ic |
---|
4343 | |
---|
4344 | def estimateT0(takestep): |
---|
4345 | Mmax = -1.e-10 |
---|
4346 | Mmin = 1.e10 |
---|
4347 | for i in range(100): |
---|
4348 | x0 = takestep(values) |
---|
4349 | M = sumREFD(x0,Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList) |
---|
4350 | Mmin = min(M,Mmin) |
---|
4351 | MMax = max(M,Mmax) |
---|
4352 | return 1.5*(MMax-Mmin) |
---|
4353 | |
---|
4354 | Q,Io,wt,Ic,Ib,Qsig = Profile[:6] |
---|
4355 | if data.get('2% weight'): |
---|
4356 | wt = 1./(0.02*Io)**2 |
---|
4357 | Qmin = Limits[1][0] |
---|
4358 | Qmax = Limits[1][1] |
---|
4359 | wtFactor = ProfDict['wtFactor'] |
---|
4360 | Bfac = data['Toler'] |
---|
4361 | Ibeg = np.searchsorted(Q,Qmin) |
---|
4362 | Ifin = np.searchsorted(Q,Qmax)+1 #include last point |
---|
4363 | Ic[:] = 0 |
---|
4364 | Bounds = {'Scale':[data['Scale'][0]*Bfac,data['Scale'][0]/Bfac],'FltBack':[0.,1.e-6], |
---|
4365 | 'DenMul':[0.,1.],'Thick':[1.,500.],'Rough':[0.,10.],'Mag SLD':[-10.,10.],'iDenMul':[-1.,1.]} |
---|
4366 | parmDict,varyList,values,bounds = GetModelParms() |
---|
4367 | Msg = 'Failed to converge' |
---|
4368 | if varyList: |
---|
4369 | if data['Minimizer'] == 'LMLS': |
---|
4370 | result = so.leastsq(calcREFD,values,full_output=True,epsfcn=1.e-8,ftol=1.e-6, |
---|
4371 | args=(Q[Ibeg:Ifin],Io[Ibeg:Ifin],wtFactor*wt[Ibeg:Ifin],Qsig[Ibeg:Ifin],parmDict,varyList)) |
---|
4372 | parmDict.update(zip(varyList,result[0])) |
---|
4373 | chisq = np.sum(result[2]['fvec']**2) |
---|
4374 | ncalc = result[2]['nfev'] |
---|
4375 | covM = result[1] |
---|
4376 | newVals = result[0] |
---|
4377 | elif data['Minimizer'] == 'Basin Hopping': |
---|
4378 | xyrng = np.array(bounds).T |
---|
4379 | take_step = RandomDisplacementBounds(xyrng[0], xyrng[1 |
---|