Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2022 The MESA Team
4 : !
5 : ! This program is free software: you can redistribute it and/or modify
6 : ! it under the terms of the GNU Lesser General Public License
7 : ! as published by the Free Software Foundation,
8 : ! either version 3 of the License, or (at your option) any later version.
9 : !
10 : ! This program is distributed in the hope that it will be useful,
11 : ! but WITHOUT ANY WARRANTY; without even the implied warranty of
12 : ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
13 : ! See the GNU Lesser General Public License for more details.
14 : !
15 : ! You should have received a copy of the GNU Lesser General Public License
16 : ! along with this program. If not, see <https://www.gnu.org/licenses/>.
17 : !
18 : ! ***********************************************************************
19 :
20 : module skye_coulomb
21 : use math_lib
22 : use math_def
23 : use auto_diff
24 : use const_def, only: dp, pi, rbohr, qe, amu, me, boltzm, five_thirds, kerg
25 : use skye_coulomb_solid
26 : use skye_coulomb_liquid
27 :
28 : implicit none
29 :
30 : logical, parameter :: dbg = .false.
31 :
32 : private
33 :
34 : public :: nonideal_corrections
35 :
36 : contains
37 :
38 : !! Computes the non-ideal free energy correction for a Coulomb system.
39 : !! This is done for both the liquid phase and the solid phase, and the resulting
40 : !! free energies are then combined in a way that blends from the solid phase when
41 : !! F_solid < F_liquid to the liquid phase when F_liquid < F_solid, with a blend width
42 : !! of order kT/abar. The blend is done in a way that propagates derivatives, so that the
43 : !! result is thermodynamically consistent.
44 : !!
45 : !! @param NMIX The number of different elements.
46 : !! @param XA The mass fractions of the different species.
47 : !! @param AZion The charges of the different species.
48 : !! @param ACMI The weight of the different species in amu.
49 : !! @param min_gamma_for_solid The minimum Gamma_i at which to use the solid free energy fit (below this, extrapolate).
50 : !! @param max_gamma_for_liquid The maximum Gamma_i at which to use the liquid free energy fit (above this, extrapolate).
51 : !! @param RHO The density of the system in g/cm^3.
52 : !! @param temp The temperature of the system in K.
53 : !! @param xnefer The electron density in 1/cm^3.
54 : !! @param abar The mean atomic mass number.
55 : !! @param dF The non-ideal correction to the free energy in erg/g.
56 : !! @param phase The blended phase. 0 for liquid, 1 for solid, smoothly interpolates in between.
57 : !! @param latent_ddlnT The latent heat of the smoothed phase transition in lnT (T dS/dlnT)
58 : !! @param latent_ddlnRho The latent heat of the smoothed phase transition in lnRho (T dS/dlnRho)
59 45520 : subroutine nonideal_corrections(NMIX,AY,AZion,ACMI, min_gamma_for_solid, max_gamma_for_liquid,&
60 : Skye_solid_mixing_rule, RHO,temp, xnefer, abar, &
61 : dF, latent_ddlnT, latent_ddlnRho,phase)
62 : integer, intent(in) :: NMIX
63 : real(dp), intent(in) :: AZion(:), ACMI(:), abar, AY(:), min_gamma_for_solid, max_gamma_for_liquid
64 : type(auto_diff_real_2var_order3), intent(in) :: RHO, temp, xnefer
65 : type(auto_diff_real_2var_order3), intent(out) :: dF, phase, latent_ddlnT, latent_ddlnRho
66 : character(len=128), intent(in) :: Skye_solid_mixing_rule
67 :
68 : integer :: IX
69 : integer :: LIQSOL
70 : real(dp) :: Zmean, Z2mean, Z52, Z53, Z321
71 : type(auto_diff_real_2var_order3) :: GAME, RS, DENS, Smix, F_phase_independent
72 : type(auto_diff_real_2var_order3) :: kT, dF_sol, dF_liq
73 :
74 : ! Compute various mean charge quantities
75 45520 : Zmean=0d0
76 45520 : Z2mean=0d0
77 45520 : Z52=0d0
78 45520 : Z53=0d0
79 45520 : Z321=0d0
80 356434 : do IX=1,NMIX
81 310914 : Zmean=Zmean+AY(IX)*AZion(IX)
82 310914 : Z2mean=Z2mean+AY(IX)*AZion(IX)*AZion(IX)
83 310914 : Z53=Z53+AY(IX)*pow(AZION(IX), five_thirds)
84 310914 : Z52=Z52+AY(IX)*pow5(sqrt(AZion(IX)))
85 356434 : Z321=Z321+AY(IX)*AZion(IX)*pow3(sqrt(AZion(IX)+1.d0))
86 : end do
87 :
88 45520 : DENS = xnefer * pow3(rbohr) ! DENS = (electrons per cubic bohr)
89 45520 : RS = pow(3d0 / (4d0 * PI * DENS),1d0/3d0) ! r_s - electron density parameter
90 45520 : GAME = qe * qe / (rbohr * boltzm * temp * RS) ! electron Coulomb parameter Gamma_e
91 :
92 : ! Electron exchange-correlation free energy
93 45520 : F_phase_independent = Zmean * EXCOR7(RS,GAME)
94 :
95 : ! Linear mixing entropy
96 45520 : Smix = linear_mixing_entropy(NMIX, Azion, AY)
97 :
98 : ! Incorporate mixing entropy into phase-independent free energy.
99 : ! Usually this would be incorporated by subtracting T*Smix (because F = E - TS),
100 : ! but here F is per ion per kB T and dSmix is per ion per kB, so
101 : ! F -> F - Smix.
102 45520 : F_phase_independent = F_phase_independent - Smix
103 :
104 :
105 : ! Compute free energy correction for liquid and solid phase.
106 45520 : LIQSOL = 0
107 : dF_liq = nonideal_corrections_phase(NMIX,AY,AZion,ACMI,min_gamma_for_solid, max_gamma_for_liquid,&
108 45520 : Skye_solid_mixing_rule, temp,abar,GAME,RS,LIQSOL,Zmean, Z2mean, Z52, Z53, Z321)
109 :
110 45520 : LIQSOL = 1
111 : dF_sol = nonideal_corrections_phase(NMIX,AY,AZion,ACMI,min_gamma_for_solid, max_gamma_for_liquid,&
112 45520 : Skye_solid_mixing_rule, temp,abar,GAME,RS,LIQSOL,Zmean, Z2mean, Z52, Z53, Z321)
113 :
114 : ! Add electron exchange-correlation energy
115 45520 : dF_liq = dF_liq + F_phase_independent
116 45520 : dF_sol = dF_sol + F_phase_independent
117 :
118 : ! Change the units from (free energy per kT per ion) to (erg/g).
119 45520 : kT = temp * kerg / (abar * amu)
120 45520 : dF_liq = dF_liq * kT
121 45520 : dF_sol = dF_sol * kT
122 :
123 : ! Produce a smoothed version of the phase transition to extract the latent heat
124 45520 : call decide_phase(dF_liq, dF_sol, kT, temp, rho, dF, phase, latent_ddlnT, latent_ddlnRho)
125 :
126 : if (dbg) then
127 : write(*,*) 'GAME',GAME%val,'Phase', phase%val
128 : end if
129 :
130 45520 : end subroutine nonideal_corrections
131 :
132 :
133 : !> Computes the free energy, phase, and latent heat across the phase transition
134 : !! between liquid and solid. The latent heat is blurred / smeared out over a finite
135 : !! width in lnT and lnRho so that the solver doesn't need to deal with a Dirac Delta.
136 : !! @param dF_liq The free energy of the liquid phase in erg/g.
137 : !! @param dF_sol The free energy of the solid phase in erg/g.
138 : !! @param kT The factor kB*T/(abar * amu), the thermal energy per mean ion mass. This is in erg/g.
139 : !! @param temp The temperature in K.
140 : !! @param rho The density in g/cm^3.
141 : !! @param phase The phase (0 for liquid, 1 for solid, in-between to represent the blur).
142 : !! @param latent_ddlnT Equals T^2 dS/dT for the latent heat. Equivalently, d(latent heat)/dlnT.
143 : !! @param latent_ddlnRho Equals T Rho dS/dRho for the latent heat. Equivalently, d(latent heat)/dlnRho.
144 45520 : subroutine decide_phase(dF_liq, dF_sol, kT, temp, rho, dF, phase, latent_ddlnT, latent_ddlnRho)
145 : type(auto_diff_real_2var_order3), intent(in) :: dF_liq, dF_sol, kT, temp, rho
146 : type(auto_diff_real_2var_order3), intent(out) :: dF, phase, latent_ddlnT, latent_ddlnRho
147 :
148 : type(auto_diff_real_2var_order3) :: blur, dF_blur, latent_S
149 : real(dp), parameter :: blur_width = 1d2
150 :
151 : ! Pick the phase with the minimum free energy
152 45520 : dF = min(dF_liq, dF_sol)
153 :
154 : ! Next, blur the free energy so it's twice differentiable across the phase
155 : ! transition, and use that to extract the latent heat.
156 45520 : blur = kT / blur_width
157 :
158 : ! Avoid overflow or underflow
159 45520 : if (dF_liq < dF_sol - 20d0*blur) then
160 45520 : phase = 0d0
161 0 : else if (dF_sol < dF_liq - 20d0*blur) then
162 0 : phase = 1d0
163 : else
164 : ! Mix phases with a blur (softmin)
165 0 : phase = exp((dF_liq-dF_sol)/blur) / (exp((dF_liq-dF_sol)/blur) + 1d0)
166 : end if
167 45520 : dF_blur = dF_liq * (1d0 - phase) + dF_sol * phase
168 :
169 : ! Now subtract off the regular (un-blurred) free energy so we don't
170 : ! double-count the off-transition dS/dT and dS/dRho in the latent heat.
171 45520 : dF_blur = dF_blur - dF
172 :
173 : ! Latent entropy
174 45520 : latent_S = -(differentiate_1(dF_blur)) ! S = -dF/dT
175 :
176 : ! T dS/dlnT = T^2 dS/dT
177 45520 : latent_ddlnT = differentiate_1(latent_S) * pow2(temp)
178 :
179 : ! T dS/dlnRho = T Rho dS/dRho
180 45520 : latent_ddlnRho = temp * rho * differentiate_2(latent_S)
181 :
182 45520 : end subroutine decide_phase
183 :
184 : !> Computes the linear mixing entropy per ion per kB.
185 : !! @param Nmix The number of species in the mixture.
186 : !! @param Azion An array of the charges of those species.
187 : !! @param AY An array of the number fractions of those species.
188 : !! @param Smix The linear mixing entropy per ion per kB.
189 45520 : type(auto_diff_real_2var_order3) function linear_mixing_entropy(Nmix, AZion, AY) result(Smix)
190 : ! Inputs
191 : integer, intent(in) :: Nmix
192 : real(dp), intent(in) :: AZion(:), AY(:)
193 :
194 : ! Intermediates and constants
195 : integer :: i
196 :
197 45520 : Smix = 0d0
198 356434 : do i=1,Nmix
199 356434 : Smix = Smix - AY(i)*log(AY(i))
200 : end do
201 45520 : end function linear_mixing_entropy
202 :
203 : !> Computes the non-ideal one-component plasma corrections to the free energy.
204 : !! This involves a loop over species to compute the free energy of a
205 : !! one-component plasma, which is added to the non-ideal mixing free energy.
206 : !!
207 : !! @param Nmix The number of species in the mixture.
208 : !! @param AY An array of the number fractions of those species.
209 : !! @param Azion An array of the charges in electron charges of those species.
210 : !! @param ACMI An array of the masses in AMU of those species.
211 : !! @param min_gamma_for_solid The minimum Gamma_i at which to use the solid free energy fit (below this, extrapolate).
212 : !! @param max_gamma_for_liquid The maximum Gamma_i at which to use the liquid free energy fit (above this, extrapolate).
213 : !! @param temp The temperature in K.
214 : !! @param abar The mean atomic mass number.
215 : !! @param RS Electron density parameter for component species
216 : !! @param GAME Electron coupling parameter (Gamma_i)
217 : !! @param LIQSOL Integer specifying the phase: 0 for liquid, 1 for solid
218 : !! @param Zmean The mean ion charge (mass fraction weighted)
219 : !! @param Z2mean The mean squared ion charge (mass fraction weighted)
220 : !! @param Z53mean The mean of ion charge to the 5/3 power (mass fraction weighted)
221 : !! @param Z321mean The mean of Z(Z+1)^(3/2), where Z is the ion charge (mass fraction weighted)
222 91040 : function nonideal_corrections_phase(NMIX,AY,AZion,ACMI,min_gamma_for_solid, max_gamma_for_liquid,Skye_solid_mixing_rule,&
223 : temp,abar,GAME,RS,LIQSOL,Zmean, Z2mean, Z52, Z53, Z321) result(dF)
224 : ! Inputs
225 : integer, intent(in) :: NMIX
226 : integer, intent(in) :: LIQSOL
227 : real(dp), intent(in) :: AZion(:), ACMI(:), abar, AY(:), Zmean, Z2mean, Z52, Z53, Z321, &
228 : min_gamma_for_solid, max_gamma_for_liquid
229 : type(auto_diff_real_2var_order3), intent(in) :: temp, GAME, RS
230 : character(len=128), intent(in) :: Skye_solid_mixing_rule
231 :
232 : ! Intermediates and constants
233 : integer :: i
234 : type(auto_diff_real_2var_order3) :: FMIX, f
235 : real(dp), parameter :: TINY=1.d-7
236 :
237 : ! Output
238 : type(auto_diff_real_2var_order3) :: dF
239 :
240 91040 : dF = 0d0
241 :
242 : ! Composition loop
243 712868 : do i=1,nmix
244 712868 : if (AY(i) > TINY .and. AZion(i) /= 0d0) then ! skip low-abundance species and neutrons
245 :
246 : ! Add up non-ideal corrections
247 621828 : f = extrapolate_free_energy(LIQSOL, temp, RS, AZion(i), ACMI(i), min_gamma_for_solid, max_gamma_for_liquid)
248 621828 : if (LIQSOL == 0) then
249 310914 : f = f + ocp_liquid_screening_free_energy_correction(AZion(i), ACMI(i), GAME, RS) ! screening corrections
250 : else
251 310914 : f = f + ocp_solid_screening_free_energy_correction(AZion(i), ACMI(i), GAME, RS) ! screening corrections
252 : end if
253 621828 : dF = dF + AY(i) * f
254 :
255 : end if
256 : end do
257 :
258 : ! Corrections to the linear mixing rule:
259 91040 : if (LIQSOL == 0) then ! liquid phase
260 45520 : FMIX = liquid_mixing_rule_correction(RS,GAME,Zmean,Z2mean,Z52,Z53,Z321)
261 : else ! solid phase (only Madelung contribution) [22.12.12]
262 45520 : FMIX = solid_mixing_rule_correction(Skye_solid_mixing_rule, NMIX, AY, AZion, GAME)
263 : end if
264 91040 : dF = dF + FMIX
265 :
266 91040 : end function nonideal_corrections_phase
267 :
268 : !> Extrapolates the free energy of a one-component plasma beyond the boundaries of the data used
269 : !! to construct the fitting functions. This is done by fixing the probability distribution over states
270 : !! to its value on the boundary, which then fixes the internal energy and entropy to their boundary values.
271 : !! The free energy is then just linear in temperature (F = F(boundary) - S(boundary) * (T - T(boundary))),
272 : !!
273 : !! @param LIQSOL Integer specifying the phase: 0 for liquid, 1 for solid
274 : !! @param temp Temperature (K)
275 : !! @param RS Electron density parameter
276 : !! @param Zion Charge of the species of interest in electron charges.
277 : !! @param CMI Mass of the species of interest in AMU.
278 : !! @param F non-ideal free energy per ion per kT
279 621828 : function extrapolate_free_energy(LIQSOL, temp, RS, Zion, CMI, min_gamma_for_solid, max_gamma_for_liquid) result(F)
280 : ! Inputs
281 : integer, intent(in) :: LIQSOL
282 : type(auto_diff_real_2var_order3), intent(in) :: temp, RS
283 : real(dp), intent(in) :: Zion, CMI, min_gamma_for_solid, max_gamma_for_liquid
284 :
285 : ! Intermediates
286 : real(dp) :: COTPT, gamma_boundary
287 : type(auto_diff_real_2var_order3) :: temp_boundary, fake_dens, GAMI, TPT, g, tp, dF_dlnT
288 :
289 : ! Output
290 : type(auto_diff_real_2var_order3) :: F
291 :
292 621828 : GAMI = pre_z(int(Zion))% z5_3 * qe * qe / (rbohr * boltzm * temp * RS) ! ion Coulomb parameter Gamma_i
293 621828 : COTPT=sqrt(3d0 * me_in_amu / CMI)/pre_z(int(Zion))%z7_6 ! auxiliary coefficient
294 621828 : TPT=GAMI*COTPT/sqrt(RS) ! T_p/T
295 :
296 621828 : if ((LIQSOL == 0 .and. GAMI < max_gamma_for_liquid) .or. (LIQSOL == 1 .and. GAMI > min_gamma_for_solid)) then
297 : ! No extrapolation needed
298 310914 : F = ocp_free_energy(LIQSOL, Zion, CMI, GAMI, TPT)
299 310914 : if (dbg) then
300 : write(*,*) 'Species:', Zion, 'LIQSOL', LIQSOL, 'Normal, GAMI:', GAMI%val, 'F', F%val
301 : end if
302 : else
303 : ! Extrapolate past the boundary
304 : if (dbg) then
305 : write(*,*) 'Species:', Zion, 'LIQSOL', LIQSOL, 'Extrapolated, GAMI:', GAMI%val
306 : end if
307 :
308 : ! Identify the boundary
309 310914 : if (LIQSOL == 0) then
310 0 : gamma_boundary = max_gamma_for_liquid
311 : else
312 310914 : gamma_boundary = min_gamma_for_solid
313 : end if
314 :
315 : ! Find the temperature where Gamma = Gamma_boundary
316 310914 : temp_boundary = temp%val * GAMI%val / gamma_boundary
317 :
318 : ! Make d(temp_boundary)/dT = 1 so we can extract dF/dT at the boundary.
319 310914 : temp_boundary%d1val1 = 1d0
320 :
321 : ! Compute new (differentiable) Gamma and TPT at the boundary
322 310914 : g = pre_z(int(Zion))% z5_3 * qe * qe / (rbohr * boltzm * temp_boundary * RS) ! ion Coulomb parameter Gamma_i
323 310914 : tp=g*COTPT/sqrt(RS) ! T_p/T
324 :
325 : ! Compute boundary free energy
326 310914 : F = ocp_free_energy(LIQSOL, Zion, CMI, g, tp)
327 :
328 : ! Extract derivative at boundary
329 310914 : dF_dlnT = differentiate_1(F) * temp_boundary
330 :
331 : ! Now make the substitution T -> T_boundary(rho).
332 : ! We do this by treating F and dF_dlnT as binary operators of (T, rho),
333 : ! and letting T = T_boundary(rho).
334 :
335 : ! First, find the temperature where Gamma = Gamma_boundary
336 : ! Note that because GAMI \propto 1/temp, this ends up being
337 : ! independent of temperature but dependent on density.
338 310914 : temp_boundary = temp * GAMI / gamma_boundary
339 :
340 : ! Next, construct a dummy operator
341 : ! with d(fake_dens)/d(dens) = 1.
342 : ! The value doesn't matter because only derivatives
343 : ! of this will be used in the chain rule for the
344 : ! substitution.
345 310914 : fake_dens = 0d0
346 310914 : fake_dens%d1val2 = 1d0
347 :
348 : ! Then make the substitution. This uses the chain rule to replace all dependences on temperature with dependences
349 : ! on temp_boundary, which in turn depends on rho. Note that the order of the arguments matters here: above we've
350 : ! made temperature the first independent variable, so the thing we're substituting for it (temp_boundary) has to be
351 : ! the first argument here.
352 : dF_dlnT = make_binop(temp_boundary, fake_dens, &
353 : dF_dlnT%val, dF_dlnT%d1val1, dF_dlnT%d1val2, dF_dlnT%d2val1, dF_dlnT%d1val1_d1val2, &
354 310914 : dF_dlnT%d2val2, dF_dlnT%d3val1, dF_dlnT%d2val1_d1val2, dF_dlnT%d1val1_d2val2, dF_dlnT%d3val2)
355 :
356 : F = make_binop(temp_boundary, fake_dens, F%val, F%d1val1, F%d1val2, F%d2val1, &
357 : F%d1val1_d1val2, F%d2val2, F%d3val1, &
358 310914 : F%d2val1_d1val2, F%d1val1_d2val2, F%d3val2)
359 :
360 : ! Extrapolate
361 310914 : F = F + (1d0 - temp_boundary / temp) * dF_dlnT
362 : end if
363 :
364 621828 : end function extrapolate_free_energy
365 :
366 : !> Calculates the free energy of a one-component
367 : !! plasma in the specified phase. This routine is
368 : !! just responsible for assembling different terms together.
369 : !! Based on EOSFI8 by Potekhin and Chabrier.
370 : !!
371 : !! @param LIQSOL Integer specifying the phase: 0 for liquid, 1 for solid.
372 : !! @param Zion Charge of the species of interest in electron charges.
373 : !! @param CMI Mass of the species of interest in AMU.
374 : !! @param GAMI Ion coupling parameter (Gamma_i)
375 : !! @param TPT effective T_p/T - ion quantum parameter
376 : !! @param F non-ideal free energy per ion per kT
377 621828 : function ocp_free_energy(LIQSOL, Zion, CMI, GAMI, TPT) result(F)
378 : ! Inputs
379 : real(dp), intent(in) :: Zion, CMI
380 : integer, intent(in) :: LIQSOL
381 : type(auto_diff_real_2var_order3), intent(in) :: GAMI, TPT
382 :
383 : ! Output
384 : type(auto_diff_real_2var_order3) :: F
385 :
386 621828 : if (LIQSOL == 0) then
387 310914 : F = classical_ocp_liquid_free_energy(GAMI) ! classical ion-ion interaction
388 310914 : F = F + quantum_ocp_liquid_free_energy_correction(TPT) ! quantum ion-ion corrections
389 : else
390 310914 : F = ocp_solid_harmonic_free_energy(GAMI,TPT) ! harmonic classical and quantum ion-ion corrections
391 310914 : F = F + ocp_solid_anharmonic_free_energy(GAMI,TPT) ! anharmonic classical and quantum ion-ion corrections
392 : end if
393 :
394 621828 : end function ocp_free_energy
395 :
396 :
397 : !> Calculates the electron exchange-correlation non-ideal free energy correction.
398 : !! Based on the results of Tanaka & Ichimaru 85-87.
399 : !! and the routine EXCOR7 Version 09.06.07 by Potekhin and Chabrier
400 : !!
401 : !! @param RS Electron density parameter for component species
402 : !! @param GAME Electron coupling parameter (Gamma_i)
403 : !! @param FXC Electron exchange-correlation non-ideal free energy correction per electron per kT
404 45520 : function EXCOR7(RS,GAME) result(FXC)
405 : ! Inputs
406 : type(auto_diff_real_2var_order3), intent(in) :: RS, GAME
407 :
408 : ! Intermediates
409 : type(auto_diff_real_2var_order3) :: THETA, SQTH, THETA2, THETA3, THETA4, EXP1TH
410 : type(auto_diff_real_2var_order3) :: CHT1, SHT1, CHT2, SHT2
411 : type(auto_diff_real_2var_order3) :: T1, T2
412 : type(auto_diff_real_2var_order3) :: A0, A1, A, B0, B1, B
413 : type(auto_diff_real_2var_order3) :: C, C3
414 : type(auto_diff_real_2var_order3) :: D0, D1, D, E0, E1, E
415 : type(auto_diff_real_2var_order3) :: DISCR, SQGE
416 : type(auto_diff_real_2var_order3) :: B2, B3, R3, S1, S2, S3, B4, C4
417 : type(auto_diff_real_2var_order3) :: S4A, S4B, S4C, S4
418 :
419 : ! Output
420 : type(auto_diff_real_2var_order3) :: FXC
421 :
422 45520 : THETA=0.543d0*RS/GAME ! non-relativistic degeneracy parameter
423 45520 : SQTH=sqrt(THETA)
424 45520 : THETA2=THETA*THETA
425 45520 : THETA3=THETA2*THETA
426 45520 : THETA4=THETA3*THETA
427 45520 : if (THETA>.007d0) then
428 45298 : CHT1=cosh(1.d0/THETA)
429 45298 : SHT1=sinh(1.d0/THETA)
430 45298 : CHT2=cosh(1.d0/SQTH)
431 45298 : SHT2=sinh(1.d0/SQTH)
432 45298 : T1=SHT1/CHT1 ! dtanh(1.d0/THETA)
433 45298 : T2=SHT2/CHT2 ! dtanh(1./sqrt(THETA))
434 : else
435 222 : T1=1.d0
436 222 : T2=1.d0
437 : end if
438 :
439 45520 : A0=0.75d0+3.04363d0*THETA2-0.09227d0*THETA3+1.7035d0*THETA4
440 45520 : A1=1d0+8.31051d0*THETA2+5.1105d0*THETA4
441 45520 : A=0.610887d0*A0/A1*T1 ! HF fit of Perrot and Dharma-wardana
442 :
443 45520 : B0=0.341308d0+12.070873d0*THETA2+1.148889d0*THETA4
444 45520 : B1=1d0+10.495346d0*THETA2+1.326623d0*THETA4
445 45520 : B=SQTH*T2*B0/B1
446 :
447 45520 : D0=0.614925d0+16.996055d0*THETA2+1.489056d0*THETA4
448 45520 : D1=1d0+10.10935d0*THETA2+1.22184d0*THETA4
449 45520 : D=SQTH*T2*D0/D1
450 :
451 45520 : E0=0.539409d0+2.522206d0*THETA2+0.178484d0*THETA4
452 45520 : E1=1d0+2.555501d0*THETA2+0.146319d0*THETA4
453 45520 : E=THETA*T1*E0/E1
454 :
455 45520 : EXP1TH=exp(-1.d0/THETA)
456 45520 : C=(0.872496d0+0.025248d0*EXP1TH)*E
457 :
458 45520 : DISCR=SQRT(4.0d0*E-D*D)
459 :
460 45520 : S1=-C/E*GAME
461 :
462 45520 : B2=B-C*D/E
463 :
464 45520 : SQGE=SQRT(GAME)
465 45520 : S2=-2.d0/E*B2*SQGE
466 :
467 45520 : R3=E*GAME+D*SQGE+1.0d0
468 :
469 45520 : B3=A-C/E
470 45520 : C3=(D/E*B2-B3)/E
471 45520 : S3=C3*log(R3)
472 :
473 45520 : B4=2.d0-D*D/E
474 :
475 45520 : C4=2d0*E*SQGE+D
476 :
477 45520 : S4A=2.0d0/E/DISCR
478 45520 : S4B=D*B3+B4*B2
479 45520 : S4C=atan(C4/DISCR)-atan(D/DISCR)
480 45520 : S4=S4A*S4B*S4C
481 :
482 45520 : FXC=S1+S2+S3+S4
483 :
484 45520 : end function EXCOR7
485 :
486 : end module skye_coulomb
|
