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_solid
21 : use math_lib
22 : use math_def
23 : use auto_diff
24 : use const_def, only: dp, pi, fine, eulernum
25 :
26 : implicit none
27 :
28 : contains
29 :
30 : !> Computes the classical and quantum anharmonic non-ideal
31 : !! free energy of a one-component plasma in the solid phase
32 : !! using classical fits due to
33 : !! Farouki & Hamaguchi'93.
34 : !!
35 : !! Quantum corrections to this are due to
36 : !! Potekhin & Chabrier 2010 (DOI 10.1002/ctpp.201010017)
37 : !!
38 : !! For a more recent reference and more mathematical detail see Appendix B.2 of
39 : !! A. Y. Potekhin and G. Chabrier: Equation of state for magnetized Coulomb plasmas
40 : !!
41 : !! Based on the routine ANHARM8 due to Potekhin and Chabrier.
42 : !!
43 : !! @param GAMI Ion coupling parameter (Gamma_i)
44 : !! @param TPT effective T_p/T - ion quantum parameter
45 : !! @param F free energy per kT per ion
46 310914 : function ocp_solid_anharmonic_free_energy(GAMI,TPT) result(F)
47 : type(auto_diff_real_2var_order3), intent(in) :: GAMI,TPT
48 : type(auto_diff_real_2var_order3) :: F
49 :
50 : real(dp), parameter :: A1 = 10.9d0
51 : real(dp), parameter :: A2 = 247d0
52 : real(dp), parameter :: A3 = 1.765d5
53 : real(dp), parameter :: B1 = 0.12d0 ! coefficient of \eta^2/\Gamma at T=0
54 :
55 : type(auto_diff_real_2var_order3) :: TPT2
56 :
57 310914 : TPT2=TPT*TPT
58 :
59 310914 : F = -(A1 / GAMI + A2 / (2d0 * pow2(GAMI)) + A3 / (3d0 * pow3(GAMI)))
60 310914 : F = F * exp(-(B1 / A1) * TPT2) ! suppress.factor of classical anharmonicity
61 310914 : F = F - B1 * TPT2 / GAMI ! Quantum correction
62 310914 : end function ocp_solid_anharmonic_free_energy
63 :
64 : !> Computes the harmonic non-ideal free energy of a
65 : !! one-component plasma in the solid phase using fits due to
66 : !! Baiko, Potekhin, & Yakovlev (2001). Includes both classical
67 : !! and quantum corrections.
68 : !!
69 : !! For a more recent reference and more mathematical detail see Appendix B.2 of
70 : !! A. Y. Potekhin and G. Chabrier: Equation of state for magnetized Coulomb plasmas
71 : !!
72 : !! Based on the routines HLfit12 and FHARM12 due to Potekhin and Chabrier.
73 : !!
74 : !! @param GAMI Ion coupling parameter (Gamma_i)
75 : !! @param TPT effective T_p/T - ion quantum parameter
76 : !! @param F free energy per kT per ion
77 310914 : function ocp_solid_harmonic_free_energy(GAMI,TPT_in) result(F)
78 : ! Inputs
79 : type(auto_diff_real_2var_order3), intent(in) :: GAMI,TPT_in
80 :
81 : ! Intermediates
82 : type(auto_diff_real_2var_order3) :: TPT, UP, DN, EA, EB, EG, E0
83 : type(auto_diff_real_2var_order3) :: Fth, U0
84 :
85 : ! Output
86 : type(auto_diff_real_2var_order3) :: F
87 :
88 : real(dp), parameter :: CM = .895929256d0 ! Madelung
89 : real(dp), parameter :: EPS=1d-5
90 :
91 : ! BCC lattice.
92 : ! Empirically the Potekhin & Chabrier fit to the BCC lattice produces
93 : ! a lower free energy than their fit to the FCC lattice in all cases of
94 : ! interest (and no cases of which I am aware), so we can specialize
95 : ! to the BCC case.
96 : real(dp), parameter :: CLM=-2.49389d0 ! 3*ln<\omega/\omega_p>
97 : real(dp), parameter :: U1=0.5113875d0
98 : real(dp), parameter :: ALPHA=0.265764d0
99 : real(dp), parameter :: BETA=0.334547d0
100 : real(dp), parameter :: GAMMA=0.932446d0
101 : real(dp), parameter :: A1=0.1839d0
102 : real(dp), parameter :: A2=0.593586d0
103 : real(dp), parameter :: A3=0.0054814d0
104 : real(dp), parameter :: A4=5.01813d-4
105 : real(dp), parameter :: A6=3.9247d-7
106 : real(dp), parameter :: A8=5.8356d-11
107 : real(dp), parameter :: B0=261.66d0
108 : real(dp), parameter :: B2=7.07997d0
109 : real(dp), parameter :: B4=0.0409484d0
110 : real(dp), parameter :: B5=0.000397355d0
111 : real(dp), parameter :: B6=5.11148d-5
112 : real(dp), parameter :: B7=2.19749d-6
113 : real(dp), parameter :: C9=0.004757014d0
114 : real(dp), parameter :: C11=0.0047770935d0
115 : real(dp), parameter :: B9=A6*C9
116 : real(dp), parameter :: B11=A8*C11
117 :
118 :
119 310914 : TPT = TPT_in
120 :
121 310914 : if (TPT > 1d0/EPS) then ! asymptote of Eq.(13) of BPY'01
122 0 : F=-1d0 / (C11*TPT*TPT*TPT)
123 310914 : else if (TPT < EPS) then ! Eq.(17) of BPY'01
124 0 : F=3d0*log(TPT)+CLM-1.5d0*U1*TPT+TPT*TPT/24.d0
125 : else
126 310914 : UP=1d0+TPT*(A1+TPT*(A2+TPT*(A3+TPT*(A4+TPT*TPT*(A6+TPT*TPT*A8)))))
127 310914 : DN=B0+TPT*TPT*(B2+TPT*TPT*(B4+TPT*(B5+TPT*(B6+TPT*(B7+TPT*TPT*(B9+TPT*TPT*B11))))))
128 :
129 310914 : EA=exp(-ALPHA*TPT)
130 310914 : EB=exp(-BETA*TPT)
131 310914 : EG=exp(-GAMMA*TPT)
132 :
133 310914 : F=log(1.d0-EA)+log(1.d0-EB)+log(1.d0-EG)-UP/DN ! F_{thermal}/NT
134 : end if
135 :
136 310914 : U0=-CM*GAMI ! perfect lattice
137 310914 : E0=1.5d0*U1*TPT ! zero-point energy
138 310914 : Fth=F+E0 ! Thermal component
139 310914 : F=U0+Fth ! Total
140 :
141 310914 : F = F -3d0 * log(TPT) + 1.5d0*log(GAMI) + 1.323515d0 ! Subtract out ideal free energy
142 : ! We use this form because it's what PC fit against
143 : ! in producing FHARM
144 :
145 310914 : end function ocp_solid_harmonic_free_energy
146 :
147 : !> Calculates an exponential with cutoffs to prevent over/underflow.
148 : !!
149 : !! @param x Input to take the exponential of.
150 : !! @param ex Output (exp(x) clipped to avoid over/underflow).
151 310914 : type(auto_diff_real_2var_order3) function safe_exp(x) result(ex)
152 : type(auto_diff_real_2var_order3), intent(in) :: x
153 310914 : ex = exp(max(-1d2,min(1d2,x)))
154 310914 : end function safe_exp
155 :
156 : !> Calculates a tanh with cutoffs to prevent over/underflow.
157 : !!
158 : !! @param x Input to take the exponential of.
159 : !! @param ex Output (exp(x) clipped to avoid over/underflow).
160 310914 : type(auto_diff_real_2var_order3) function safe_tanh(x) result(th)
161 : type(auto_diff_real_2var_order3), intent(in) :: x
162 310914 : th = tanh(max(-1d2,min(1d2,x)))
163 310914 : end function safe_tanh
164 :
165 : !> Calculates the electron-ion screening corrections to the free energy
166 : !! of a one-component plasma in the solid phase using the fits of Potekhin & Chabrier 2013.
167 : !! In the non-degenerate limit we smoothly transition to the liquid-phase screening
168 : !! free energy, which is just an easy way to ensure we reproduce the Debye-Huckel limit.
169 : !!
170 : !! @param Z ion charge
171 : !! @param mi ion mass in amu
172 : !! @param ge electron interaction parameter
173 : !! @param rs non-dimensionalized electron radius
174 : !! @param F non-ideal free energy
175 310914 : function ocp_solid_screening_free_energy_correction(Z, mi, ge, rs) result(F)
176 : use skye_coulomb_liquid, only: me_in_amu, ocp_liquid_screening_free_energy_correction
177 : real(dp), intent(in) :: Z, mi
178 : type(auto_diff_real_2var_order3), intent(in) :: ge, rs
179 :
180 310914 : real(dp) :: s, b1, b2, b3, b4, COTPT
181 : real(dp), parameter :: aTF = 0.00352d0
182 :
183 : type(auto_diff_real_2var_order3) :: TPT, f_inf, A, Q, xr, supp, g, alpha, Fliq, gr, switch
184 : type(auto_diff_real_2var_order3) :: F
185 :
186 310914 : s = 1d0 / (1d0 + 1d-2 * pre_z(int(Z))% logz_3_2 + 0.097d0 / pre_z(int(Z))% z2)
187 310914 : b1 = 1d0 - 1.1866d0 * pre_z(int(Z))% zm0p267 + 0.27d0 / Z
188 : b2 = 1d0 + (2.25d0 * pre_z(int(Z))% zm1_3) * (1d0 + 0.684d0 * pre_z(int(Z))% z5 + 0.222d0 * pre_z(int(Z))% z6) &
189 310914 : / (1d0 + 0.222d0 * pre_z(int(Z))% z6)
190 310914 : b3 = 41.5d0 / (1d0 + pre_z(int(Z))% logz)
191 310914 : b4 = 0.395d0 * pre_z(int(Z))% logz + 0.347d0 * pre_z(int(Z))% zm3_2
192 :
193 310914 : g = ge * pre_z(int(Z))% z5_3
194 :
195 310914 : COTPT = sqrt(3d0 * me_in_amu / mi) / pre_z(int(Z))% z7_6
196 310914 : TPT = g * COTPT / sqrt(RS)
197 310914 : supp = safe_exp(-pow2(0.205d0 * TPT))
198 310914 : Q = sqrt((pow2(0.205d0 * TPT) + log(1d0 + supp)) / log(eulernum - (eulernum - 2d0) * supp))
199 :
200 310914 : xr = pow(9d0 * pi / 4d0, 1d0/3d0) * fine / rs
201 310914 : A = (b3 + 17.9d0 * pow2(xr)) / (1d0 + b4 * pow2(xr))
202 310914 : f_inf = aTF * pre_z(int(Z))% z2_3 * b1 * sqrt(1d0 + b2 / pow2(xr))
203 :
204 310914 : F = -f_inf * g * (1d0 + A * pow(Q / g, s))
205 :
206 310914 : gr = sqrt(1d0 + pow2(xr))
207 310914 : alpha = 3d0 * pow(4d0 / (9d0 * pi), 2d0/3d0) * (rs / ge) * gr
208 :
209 310914 : Fliq = ocp_liquid_screening_free_energy_correction(Z, mi, ge, rs)
210 :
211 310914 : switch = pow3(safe_tanh(2d0*alpha))
212 310914 : F = switch * Fliq + (1d0 - switch) * F
213 :
214 310914 : end function ocp_solid_screening_free_energy_correction
215 :
216 : !> Computes the correction deltaG to the linear mixing rule for a two-component Coulomb solid mixture.
217 : !! From Shuji Ogata, Hiroshi Iyetomi, and Setsuo Ichimaru 1993
218 : !!
219 : !! Based on simulations done for charge ratio 4/3 <= Rz <= 4 and 163 <= Gamma <= 383, where
220 : !! Gamma here is the ionic interaction strength measured for the lower-charge species.
221 : !! deltaG is defined as deltaF / (x1*x2*Gamma), where deltaF is a free energy per kT per ion
222 : !! and x1 is the abundance of species 1.
223 : !!
224 : !! @param x2 Abundance of the higher-charge species
225 : !! @param Rz Charge ratio of species (> 1 by definition).
226 0 : real(dp) function deltaG_Ogata93(x2, Rz) result(dG)
227 : ! Inputs
228 : real(dp), intent(in) :: x2
229 : real(dp), intent(in) :: Rz
230 :
231 : ! Intermediates
232 0 : real(dp) :: CR
233 :
234 0 : CR = 0.05d0 * pow2(Rz - 1d0) / ((1d0 + 0.64d0 * (Rz - 1d0)) * (1d0 + 0.5d0 * pow2(Rz - 1d0)))
235 : dG = CR / (1 + &
236 : (sqrt(x2) * (sqrt(x2) - 0.3d0) * (sqrt(x2) - 0.7d0) * (sqrt(x2) - 1d0)) &
237 0 : * 27d0 * (Rz - 1d0) / (1d0 + 0.1d0 * (Rz - 1d0)))
238 :
239 310914 : end function deltaG_Ogata93
240 :
241 : !> Computes the correction deltaG to the linear mixing rule for a two-component Coulomb solid mixture.
242 : !! Originally from PhysRevE.79.016411 (Equation of state of classical Coulomb plasma mixtures)
243 : !! by Potekhin, Alexander Y. and Chabrier, Gilles and Rogers, Forrest J.
244 : !! https://link.aps.org/doi/10.1103/PhysRevE.79.016411
245 : !!
246 : !! @param x Abundance of the higher-charge species
247 : !! @param Rz Charge ratio of species (> 1 by definition).
248 1300910 : real(dp) function deltaG_PC13(x2, Rz) result(dG)
249 : ! Inputs
250 : real(dp), intent(in) :: x2
251 : real(dp), intent(in) :: Rz
252 :
253 : ! Intermediates
254 1300910 : real(dp) :: x
255 :
256 1300910 : x = x2 / Rz + (1d0 - 1d0 / Rz) * pow(x2, Rz)
257 1300910 : dG = 0.012d0 * ((x*(1d0-x)) / (x2*(1d0-x2))) * (1d0 - 1d0/pow2(Rz)) * (1d0 - x2 + x2 * pow(Rz,5d0/3d0))
258 :
259 1300910 : end function deltaG_PC13
260 :
261 : !> Calculates the correction to the linear mixing rule for a Coulomb solid mixture
262 : !! by extending a two-component deltaG prescription to the multi-component case, using the
263 : !! prescription of Medin & Cumming 2010.
264 : !!
265 : !! @param n Number of species
266 : !! @param AY Array of length NMIX holding the masses of species
267 : !! @param AZion Array of length NMIX holding the charges of species
268 : !! @param GAME election interaction parameter
269 : !! @param F mixing free energy correction per ion per kT.
270 45520 : function solid_mixing_rule_correction(Skye_solid_mixing_rule, n, AY, AZion, GAME) result(F)
271 : ! Inputs
272 : character(len=128), intent(in) :: Skye_solid_mixing_rule
273 : integer, intent(in) :: n
274 : real(dp), intent(in) :: AZion(:), AY(:)
275 : type(auto_diff_real_2var_order3), intent(in) :: GAME
276 :
277 : ! Intermediates
278 : integer :: i,j, num_unique_charges
279 712868 : real(dp) :: unique_charges(n), charge_abundances(n)
280 : logical :: found
281 : integer :: found_index
282 45520 : real(dp) :: RZ, aj, dG
283 : type(auto_diff_real_2var_order3) :: GAMI
284 :
285 : ! Output
286 : type(auto_diff_real_2var_order3) :: F
287 :
288 : ! Parameters
289 : real(dp), parameter :: C = 0.012d0
290 : real(dp), parameter :: eps = 1d-40
291 :
292 : ! Identify and group unique charges
293 45520 : num_unique_charges = 0
294 356434 : do i=1,n
295 1300910 : found = .false.
296 :
297 1300910 : do j=1,num_unique_charges
298 1300910 : if (unique_charges(j) == AZion(i)) then
299 : found = .true.
300 : found_index = j
301 : exit
302 : end if
303 : end do
304 :
305 356434 : if (.not. found) then
306 310914 : num_unique_charges = num_unique_charges + 1
307 310914 : unique_charges(num_unique_charges) = AZion(i)
308 310914 : charge_abundances(num_unique_charges) = AY(i)
309 : else
310 0 : charge_abundances(found_index) = charge_abundances(found_index) + AY(i)
311 : end if
312 : end do
313 :
314 45520 : F = 0d0
315 356434 : do i=1,num_unique_charges
316 310914 : if (unique_charges(i) == 0d0) cycle
317 2647340 : do j=1,num_unique_charges
318 2290906 : if (unique_charges(j) == 0d0) cycle
319 :
320 : ! Expression needs R > 1.
321 : ! From PC2013: 'dF_sol = Sum_i Sum_{j>i} ... where the indices are arranged so that Z_j < Z_{j+1}'
322 : ! So if j > i then Z_j > Z_i, which means R = Z_j / Z_i > 1.
323 : ! We extend to the case of equality by grouping equal-charge species together, as above.
324 2290906 : if (unique_charges(j) < unique_charges(i)) cycle
325 :
326 1300910 : RZ = unique_charges(j)/unique_charges(i) ! Charge ratio
327 :
328 : ! max avoids divergence.
329 : ! The contribution to F scales as abundance_sum^2, so in cases where the max returns eps
330 : ! we don't care much about the error this incurs.
331 1300910 : aj = charge_abundances(j) / max(eps, charge_abundances(i) + charge_abundances(j)) ! = x2 / (x1 + x2) in MC10's language
332 1300910 : if (Skye_solid_mixing_rule == 'Ogata') then
333 0 : dG = deltaG_Ogata93(aj, RZ)
334 1300910 : else if (Skye_solid_mixing_rule == 'PC') then
335 1300910 : dG = deltaG_PC13(aj, RZ)
336 : else
337 0 : write(*,*) 'Error: Invalid choice for Skye_solid_mixing_rule.'
338 0 : stop
339 : end if
340 :
341 1300910 : GAMI=pow(unique_charges(i),5d0/3d0)*GAME
342 2601820 : F = F + GAMI * (charge_abundances(i) * charge_abundances(j) * dG)
343 : end do
344 : end do
345 45520 : end function solid_mixing_rule_correction
346 :
347 :
348 : end module skye_coulomb_solid
|
