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_liquid
21 : use math_lib
22 : use math_def
23 : use auto_diff
24 : use const_def, only: dp, me, amu, fine, pi, ev2erg
25 :
26 : implicit none
27 :
28 : real(dp), parameter :: me_in_amu = me / amu
29 :
30 : contains
31 :
32 : !> Calculates the free energy of a classical one-component
33 : !! plasma in the liquid phase using the fitting form of
34 : !! A.Y.Potekhin and G.Chabrier, Phys.Rev.E62, 8554 (2000)
35 : !! fit to the data of DeWitt & Slattery (1999).
36 : !! This choice ensures consistency with the quantum_ocp_liquid_free_energy_correction
37 : !! routine, which is based on the fits of Baiko & Yakolev 2019 (who were correcting relative
38 : !! to the DeWitt & Slattery fits).
39 : !!
40 : !! @param g ion interaction parameter
41 : !! @param F non-ideal free energy
42 310914 : function classical_ocp_liquid_free_energy(g) result(F)
43 : type(auto_diff_real_2var_order3), intent(in) :: g
44 : type(auto_diff_real_2var_order3) :: FA, FB
45 : type(auto_diff_real_2var_order3) :: F
46 :
47 : real(dp), parameter :: SQ32=sqrt(3d0) / 2d0
48 :
49 : real(dp), parameter :: A1=-0.9070d0
50 : real(dp), parameter :: A2=0.62954d0
51 : real(dp), parameter :: A3 = -SQ32 - A1 / sqrt(A2)
52 :
53 : real(dp), parameter :: B1 = 4.56d-3
54 : real(dp), parameter :: B2 = 211.6d0
55 : real(dp), parameter :: B3 = -1d-4
56 : real(dp), parameter :: B4 = 4.62d-3
57 :
58 : FA = A1 * (sqrt(g * (A2 + g)) - A2 * log(sqrt(g / A2) + sqrt(1 + g / A2))) &
59 310914 : + 2d0 * A3 * (sqrt(g) - atan(sqrt(g)))
60 310914 : FB = B1 * (g - B2 * log(1d0 + g / B2)) + 0.5d0 * B3 * log(1d0 + pow2(g) / B4)
61 310914 : F = FA + FB
62 :
63 310914 : end function classical_ocp_liquid_free_energy
64 :
65 :
66 : !> Calculates the quantum corrections to the free energy of a
67 : !! one-component plasma in the liquid phase using the fits due to
68 : !! Baiko & Yakovlev 2019.
69 : !!
70 : !! @param TPT effective T_p/T - ion quantum parameter
71 : !! @param F non-ideal free energy
72 310914 : function quantum_ocp_liquid_free_energy_correction(TPT) result(F)
73 : type(auto_diff_real_2var_order3), intent(in) :: TPT
74 : type(auto_diff_real_2var_order3) :: eta
75 : type(auto_diff_real_2var_order3) :: F
76 :
77 : real(dp), parameter :: Q1 = 5.994d0
78 : real(dp), parameter :: Q2 = 70.3d0
79 : real(dp), parameter :: Q4 = 22.7d0
80 : real(dp), parameter :: Q3 = (0.25d0 - Q1 / Q2) * Q4
81 :
82 310914 : eta = TPT / sqrt(3d0) ! Note that this is the BK19 definition of eta. PC typically use eta = TPT.
83 :
84 310914 : F = Q1 * eta - Q1 * Q2 * log(1d0 + eta / Q2) + 0.5d0 * Q3 * log(1d0 + pow2(eta) / Q4)
85 310914 : end function quantum_ocp_liquid_free_energy_correction
86 :
87 : !> Calculates the electron-ion screening corrections to the free energy
88 : !! of a one-component plasma in the liquid phase using the fits of Potekhin & Chabrier 2013.
89 : !!
90 : !! @param Z ion charge
91 : !! @param mi ion mass in amu
92 : !! @param ge electron interaction parameter
93 : !! @param rs non-dimensionalized electron radius
94 : !! @param F non-ideal free energy
95 621828 : function ocp_liquid_screening_free_energy_correction(Z, mi, ge, rs) result(F)
96 : real(dp), intent(in) :: Z, mi
97 : type(auto_diff_real_2var_order3), intent(in) :: ge, rs
98 :
99 621828 : real(dp) :: cDH, cTF, a, b, nu, COTPT
100 :
101 : type(auto_diff_real_2var_order3) :: TPT, g, g1, g2, h, gr, xr
102 : type(auto_diff_real_2var_order3) :: F
103 :
104 621828 : a = 1.11d0 * pre_z(int(Z))% z0p475
105 621828 : b = 0.2d0 + 0.078d0 * pre_z(int(Z))% sqlogz
106 621828 : nu = 1.16d0 + 0.08d0 * pre_z(int(Z))% logz
107 621828 : cDH = (Z / sqrt(3d0)) * (pre_z(int(Z))% zp1_3_2 - 1d0 - pre_z(int(Z))% z3_2)
108 : cTF = (18d0 / 175d0) * pow(12d0 / pi, 2d0/3d0) * pre_z(int(Z))% z7_3 &
109 621828 : * (1d0 - pre_z(int(Z))% zm1_3 + 0.2d0 * pre_z(int(Z))% zm1_2)
110 :
111 621828 : g = ge * pre_z(int(Z))% z5_3
112 621828 : COTPT = sqrt(3d0 * me_in_amu / mi) / pre_z(int(Z))% z7_6
113 621828 : TPT = g * COTPT / sqrt(rs)
114 :
115 621828 : xr = pow(9d0 * pi / 4d0, 1d0/3d0) * fine / rs
116 621828 : gr = sqrt(1d0 + pow2(xr))
117 621828 : g1 = 1d0 + 0.78d0 * sqrt(ge / z) / (21d0 + ge * pow3(Z / rs))
118 621828 : g2 = 1d0 + ((Z - 1d0) / 9d0) * (1d0 + 1d0 / (0.001d0 * pre_z(int(Z))% z2 + 2d0 * ge)) * (pow3(rs) / (1d0 + 6d0 * pow2(rs)))
119 : h = (1d0 + 0.2d0 * pow2(xr)) / &
120 621828 : (1d0 + 0.18d0 * xr * pre_z(int(Z))% zm1_4 + 0.37d0 * pre_z(int(Z))% zm1_2 * pow2(xr) + 0.2d0 * pow2(xr))
121 :
122 621828 : F = -ge * (cDH * sqrt(ge) + cTF * a * pow(ge, nu) * g1 * h) / (1d0 + (b * sqrt(ge) + a * g2 * pow(ge, nu) / rs) / gr)
123 :
124 621828 : end function ocp_liquid_screening_free_energy_correction
125 :
126 :
127 : !> Calculates the correction to the linear mixing rule for a Coulomb liquid mixture.
128 : !! Based on CORMIX Version 02.07.09 by Potekhin and Chabrier
129 : !!
130 : !! @param RS Electron density parameter for component species
131 : !! @param GAME Electron coupling parameter (Gamma_i)
132 : !! @param Zmean Mean charge of ions.
133 : !! @param Z2mean Mean of squared charge of ions.
134 : !! @param Z52 Mean of Z^(5/2) for ions.
135 : !! @param Z53 Mean of Z^(5/3) for ions.
136 : !! @param Z321 Mean of Z*(1+Z)^(3/2) for ions.
137 : !! @param FMIX mixing free energy correction per ion per kT.
138 45520 : function liquid_mixing_rule_correction(RS,GAME,Zmean,Z2mean,Z52,Z53,Z321) result(FMIX)
139 : type(auto_diff_real_2var_order3), intent(in) :: RS,GAME
140 : real(dp), intent(in) :: Zmean,Z2mean,Z52,Z53,Z321
141 :
142 : type(auto_diff_real_2var_order3) :: GAMImean, Dif0, DifR, DifFDH, D
143 : type(auto_diff_real_2var_order3) :: P3, D0, GP, FMIX0, Q, R, GQ
144 :
145 : type(auto_diff_real_2var_order3) :: FMIX
146 :
147 : real(dp), parameter :: TINY = 1d-9
148 :
149 :
150 : ! Skip balls of neutrons
151 45520 : if(zmean==0) then
152 0 : FMIX=0d0
153 0 : return
154 : end if
155 :
156 45520 : GAMImean=GAME*Z53
157 45520 : if (RS<TINY) then ! OCP
158 0 : Dif0=Z52-sqrt(Z2mean*Z2mean*Z2mean/Zmean)
159 : else
160 45520 : Dif0=Z321-sqrt((Z2mean+Zmean)*(Z2mean+Zmean)*(Z2mean+Zmean)/Zmean)
161 : end if
162 45520 : DifR=Dif0/Z52
163 45520 : DifFDH=Dif0*GAME*sqrt(GAME/3d0) ! F_DH - F_LM(DH)
164 45520 : D=Z2mean/(Zmean*Zmean)
165 45520 : if (abs(D-1.d0)<TINY) then ! no correction
166 408 : FMIX=0d0
167 408 : return
168 : end if
169 45112 : P3=pow(D,-0.2d0)
170 45112 : D0=(2.6d0*DifR+14d0*DifR*DifR*DifR)/(1.d0-P3)
171 45112 : GP=D0*pow(GAMImean,P3)
172 45112 : FMIX0=DifFDH/(1d0+GP)
173 :
174 45112 : Q=D*D*0.0117d0
175 45112 : R=1.5d0/P3-1d0
176 45112 : GQ=Q*GP
177 45112 : FMIX=FMIX0/pow(1d0+GQ,R)
178 45112 : end function liquid_mixing_rule_correction
179 :
180 : end module skye_coulomb_liquid
181 :
|
