Line data Source code
1 : ! ***********************************************************************
2 : !
3 : ! Copyright (C) 2020 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 create_EXCOR7_table
21 :
22 : use const_def, only: dp
23 : use chem_def
24 : use utils_lib, only: is_bad
25 : use math_lib
26 :
27 : implicit none
28 :
29 : public :: do_create_EXCOR7_table
30 : private
31 :
32 :
33 : contains
34 :
35 0 : subroutine do_create_EXCOR7_table(fname)
36 : character (len=*), intent(in) :: fname
37 : real(dp) :: logRS, logRS_min, logRS_max, dlogRS
38 0 : real(dp) :: logGAME, logGAME_min, logGAME_max, dlogGAME
39 : integer :: nlogRS, nlogGAME, io_unit, i, j
40 0 : real(dp) :: RS, GAME, FXC, UXC, PXC, CVXC, SXC, PDTXC, PDRXC
41 : include 'formats'
42 :
43 0 : logRS_min = -3.5d0
44 0 : logRS_max = 0.0d0
45 0 : dlogRS = 1d-2
46 0 : nlogRS = (logRS_max - logRS_min)/dlogRS + 1
47 :
48 0 : logGAME_min = -2d0
49 0 : logGAME_max = 4d0
50 0 : dlogGAME = 1d-2
51 0 : nlogGAME = (logGAME_max - logGAME_min)/dlogGAME + 1
52 :
53 : !write(*,'(a)') 'create ' // trim(fname)
54 :
55 0 : open(newunit=io_unit,file=trim(fname))
56 :
57 0 : write(io_unit,'(99(a14))') 'num logRS', 'logRS min', 'logRS max', 'del logRS', &
58 0 : 'num logGAME', 'logGAME min', 'logGAME max', 'del logGAME'
59 :
60 : write(io_unit,'(2(i10,4x,3(f14.4)),i10)') &
61 0 : nlogRS, logRS_min, logRS_max, dlogRS, &
62 0 : nlogGAME, logGAME_min, logGAME_max, dlogGAME, 0
63 :
64 0 : do i = 1, nlogRS
65 0 : logRS = logRS_min + (i-1) * dlogRS
66 0 : RS = exp10(logRS)
67 0 : write(io_unit,'(/,7x,a)') 'logRS'
68 0 : write(io_unit,'(2x,f14.6/)') logRS
69 : write(io_unit,'(99(a26))') &
70 0 : 'logGAME', 'FXC', 'UXC', 'PXC', 'CVXC', 'SXC', 'PDTXC', 'PDRXC', 'RS', 'GAME'
71 0 : do j = 1, nlogGAME
72 0 : logGAME = logGAME_min + (j-1) * dlogGAME
73 0 : GAME = exp10(logGAME)
74 0 : if (is_bad(GAME)) then
75 0 : write(*,3) 'GAME', i, j, GAME, logGAME, logGAME_min, dlogGAME
76 : end if
77 0 : call EXCOR7(RS, GAME, FXC, UXC, PXC, CVXC, SXC, PDTXC, PDRXC)
78 : write(io_unit,'(99(1pd26.16))') &
79 0 : logGAME, FXC, UXC, PXC, CVXC, SXC, PDTXC, PDRXC, RS, GAME
80 : end do
81 0 : write(io_unit,*)
82 : end do
83 :
84 0 : write(io_unit,*)
85 0 : write(io_unit,*)
86 0 : close(io_unit)
87 :
88 0 : end subroutine do_create_EXCOR7_table
89 :
90 :
91 : ! ============== ELECTRON EXCHANGE AND CORRELATION ================ *
92 0 : subroutine EXCOR7(RS,GAME,FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC)
93 : ! Version 09.06.07
94 : ! Exchange-correlation contribution for the electron gas
95 : ! Stems from TANAKA1 v.03.03.96. Added derivatives.
96 : ! Input: RS - electron density parameter =electron-sphere radius in a.u.
97 : ! GAME - electron Coulomb coupling parameter
98 : ! Output: FXC - excess free energy of e-liquid per kT per one electron
99 : ! according to Tanaka & Ichimaru 85-87 and Ichimaru 93
100 : ! UXC - internal energy contr.[per 1 electron, kT]
101 : ! PXC - pressure contribution divided by (n_e kT)
102 : ! CVXC - heat capacity divided by N_e k
103 : ! SXC - entropy divided by N_e k
104 : ! PDTXC,PDRXC = PXC + d ln PXC / d ln(T,\rho)
105 : real(dp), intent(in) :: RS, GAME
106 : real(dp), intent(out) :: FXC,UXC,PXC,CVXC,SXC,PDTXC,PDRXC
107 :
108 0 : real(dp) :: THETA, SQTH, THETA2, THETA3, THETA4, EXP1TH
109 0 : real(dp) :: CHT1, SHT1, CHT2, SHT2
110 0 : real(dp) :: T1, T2, T1DH, T1DHH, T2DH, T2DHH
111 0 : real(dp) :: A0, A0DH, A0DHH, A1, A1DH, A1DHH, A, AH, ADH, ADHH
112 0 : real(dp) :: B0, B0DH, B0DHH, B1, B1DH, B1DHH, B, BH, BDH, BDHH
113 0 : real(dp) :: C, CDH, CDHH, C3, C3DH, C3DHH
114 0 : real(dp) :: D0, D0DH, D0DHH, D1, D1DH, D1DHH, D, DH, DDH, DDHH
115 0 : real(dp) :: E0, E0DH, E0DHH, E1, E1DH, E1DHH, E, EH, EDH, EDHH
116 0 : real(dp) :: DISCR, DIDH, DIDHH
117 0 : real(dp) :: S1, S1H, S1DH, S1DHH, S1DG, S1DHG
118 0 : real(dp) :: B2, B2DH, B2DHH, SQGE, B3, B3DH, B3DHH
119 0 : real(dp) :: S2, S2H, S2DH, S2DHH, S2DG, S2DGG, S2DHG
120 0 : real(dp) :: R3, R3DH, R3DHH, R3DG, R3DGG, R3DHG
121 0 : real(dp) :: S3, S3DH, S3DHH, S3DG, S3DGG, S3DHG
122 0 : real(dp) :: B4, B4DH, B4DHH
123 0 : real(dp) :: C4, C4DH, C4DHH, C4DG, C4DGG, C4DHG
124 0 : real(dp) :: S4A, S4AH, S4ADH, S4ADHH
125 0 : real(dp) :: S4B, S4BDH, S4BDHH, UP1, DN1, UP2, DN2
126 0 : real(dp) :: S4C, S4CDH, S4CDHH, S4CDG, S4CDGG, S4CDHG
127 0 : real(dp) :: S4, S4DH, S4DHH, S4DG, S4DGG, S4DHG
128 0 : real(dp) :: FXCDH, FXCDHH, FXCDG, FXCDGG, FXCDHG
129 0 : real(dp) :: PDLH, PDLG
130 : include 'formats'
131 0 : THETA=0.543d0*RS/GAME ! non-relativistic degeneracy parameter
132 0 : SQTH=sqrt(THETA)
133 0 : THETA2=THETA*THETA
134 0 : THETA3=THETA2*THETA
135 0 : THETA4=THETA3*THETA
136 0 : if (THETA>.005d0) then
137 0 : CHT1=cosh(1.d0/THETA)
138 0 : SHT1=sinh(1.d0/THETA)
139 0 : CHT2=cosh(1.d0/SQTH)
140 0 : SHT2=sinh(1.d0/SQTH)
141 0 : T1=SHT1/CHT1 ! dtanh(1.d0/THETA)
142 0 : T2=SHT2/CHT2 ! dtanh(1./sqrt(THETA))
143 0 : T1DH=-1.d0/((THETA*CHT1)*(THETA*CHT1)) ! d T1 / d\theta
144 0 : T1DHH=2.d0/pow3(THETA*CHT1)*(CHT1-SHT1/THETA)
145 0 : T2DH=-0.5d0*SQTH/((THETA*CHT2)*(THETA*CHT2))
146 0 : T2DHH=(0.75d0*SQTH*CHT2-0.5d0*SHT2)/pow3(THETA*CHT2)
147 : else
148 : T1=1.d0
149 : T2=1.d0
150 : T1DH=0.d0
151 : T2DH=0.d0
152 : T1DHH=0.d0
153 : T2DHH=0.d0
154 : end if
155 0 : A0=0.75d0+3.04363d0*THETA2-0.09227d0*THETA3+1.7035d0*THETA4
156 0 : A0DH=6.08726d0*THETA-0.27681d0*THETA2+6.814d0*THETA3
157 0 : A0DHH=6.08726d0-0.55362d0*THETA+20.442d0*THETA2
158 0 : A1=1d0+8.31051d0*THETA2+5.1105d0*THETA4
159 0 : A1DH=16.62102d0*THETA+20.442d0*THETA3
160 0 : A1DHH=16.62102d0+61.326d0*THETA2
161 0 : A=0.610887d0*A0/A1*T1 ! HF fit of Perrot and Dharma-wardana
162 0 : AH=A0DH/A0-A1DH/A1+T1DH/T1
163 0 : ADH=A*AH
164 : ADHH=ADH*AH+A*(A0DHH/A0-pow2(A0DH/A0)-A1DHH/A1+pow2(A1DH/A1) &
165 0 : + T1DHH/T1-pow2(T1DH/T1))
166 0 : B0=0.341308d0+12.070873d0*THETA2+1.148889d0*THETA4
167 0 : B0DH=24.141746d0*THETA+4.595556d0*THETA3
168 0 : B0DHH=24.141746d0+13.786668d0*THETA2
169 0 : B1=1d0+10.495346d0*THETA2+1.326623d0*THETA4
170 0 : B1DH=20.990692d0*THETA+5.306492d0*THETA3
171 0 : B1DHH=20.990692d0+15.919476d0*THETA2
172 0 : B=SQTH*T2*B0/B1
173 0 : BH=0.5d0/THETA+T2DH/T2+B0DH/B0-B1DH/B1
174 0 : BDH=B*BH
175 : BDHH=BDH*BH+B*(-0.5d0/THETA2+T2DHH/T2-pow2(T2DH/T2) &
176 0 : + B0DHH/B0-pow2(B0DH/B0)-B1DHH/B1+pow2(B1DH/B1))
177 0 : D0=0.614925d0+16.996055d0*THETA2+1.489056d0*THETA4
178 0 : D0DH=33.99211d0*THETA+5.956224d0*THETA3
179 0 : D0DHH=33.99211d0+17.868672d0*THETA2
180 0 : D1=1d0+10.10935d0*THETA2+1.22184d0*THETA4
181 0 : D1DH=20.2187d0*THETA+4.88736d0*THETA3
182 0 : D1DHH=20.2187d0+14.66208d0*THETA2
183 0 : D=SQTH*T2*D0/D1
184 0 : DH=0.5d0/THETA+T2DH/T2+D0DH/D0-D1DH/D1
185 0 : DDH=D*DH
186 : DDHH=DDH*DH+D*(-0.5d0/THETA2+T2DHH/T2-pow2(T2DH/T2) &
187 0 : + D0DHH/D0-pow2(D0DH/D0)-D1DHH/D1+pow2(D1DH/D1))
188 0 : E0=0.539409d0+2.522206d0*THETA2+0.178484d0*THETA4
189 0 : E0DH=5.044412d0*THETA+0.713936d0*THETA3
190 0 : E0DHH=5.044412d0+2.141808d0*THETA2
191 0 : E1=1d0+2.555501d0*THETA2+0.146319d0*THETA4
192 0 : E1DH=5.111002d0*THETA+0.585276d0*THETA3
193 0 : E1DHH=5.111002d0+1.755828d0*THETA2
194 0 : E=THETA*T1*E0/E1
195 0 : EH=1.d0/THETA+T1DH/T1+E0DH/E0-E1DH/E1
196 0 : EDH=E*EH
197 : EDHH=EDH*EH+E*(T1DHH/T1-pow2(T1DH/T1)+E0DHH/E0-pow2(E0DH/E0) &
198 0 : - E1DHH/E1+pow2(E1DH/E1)-1.0d0/THETA2)
199 0 : EXP1TH=exp(-1.d0/THETA)
200 0 : C=(0.872496d0+0.025248d0*EXP1TH)*E
201 0 : CDH=0.025248d0*EXP1TH/THETA2*E+C*EDH/E
202 : CDHH=0.025248d0*EXP1TH/THETA2*(EDH+(1.0d0-2.0d0*THETA)/THETA2*E) &
203 0 : + CDH*EDH/E+C*EDHH/E-C*pow2(EDH/E)
204 0 : DISCR=SQRT(4.0d0*E-D*D)
205 0 : DIDH=0.5d0/DISCR*(4.0d0*EDH-2.0d0*D*DDH)
206 0 : DIDHH=(-pow2((2.0d0*EDH-D*DDH)/DISCR)+2.0d0*EDHH-DDH*DDH-D*DDHH)/DISCR
207 0 : S1=-C/E*GAME
208 0 : S1H=CDH/C-EDH/E
209 0 : S1DH=S1*S1H
210 0 : S1DHH=S1DH*S1H+S1*(CDHH/C-pow2(CDH/C)-EDHH/E+pow2(EDH/E))
211 0 : S1DG=-C/E ! => S1DGG=0
212 0 : S1DHG=S1DG*(CDH/C-EDH/E)
213 0 : B2=B-C*D/E
214 0 : B2DH=BDH-(CDH*D+C*DDH)/E+C*D*EDH/(E*E)
215 0 : B2DHH=BDHH-(CDHH*D+2d0*CDH*DDH+C*DDHH)/E+(2d0*(CDH*D+C*DDH-C*D*EDH/E)*EDH+C*D*EDHH)/(E*E)
216 0 : SQGE=SQRT(GAME)
217 0 : S2=-2.d0/E*B2*SQGE
218 0 : S2H=B2DH/B2-EDH/E
219 0 : S2DH=S2*S2H
220 0 : S2DHH=S2DH*S2H+S2*(B2DHH/B2-pow2(B2DH/B2)-EDHH/E+pow2(EDH/E))
221 0 : S2DG=0.5d0*S2/GAME
222 0 : S2DGG=-0.5d0*S2DG/GAME
223 0 : S2DHG=0.5d0*S2DH/GAME
224 0 : R3=E*GAME+D*SQGE+1.0d0
225 0 : R3DH=EDH*GAME+DDH*SQGE
226 0 : R3DHH=EDHH*GAME+DDHH*SQGE
227 0 : R3DG=E+0.5d0*D/SQGE
228 0 : R3DGG=-0.25d0*D/(GAME*SQGE)
229 0 : R3DHG=EDH+0.5d0*DDH/SQGE
230 0 : B3=A-C/E
231 0 : B3DH=ADH-CDH/E+C*EDH/(E*E)
232 0 : B3DHH=ADHH-CDHH/E+(2.0d0*CDH*EDH+C*EDHH)/(E*E)-2d0*C*EDH*EDH/pow3(E)
233 0 : C3=(D/E*B2-B3)/E
234 0 : C3DH=(DDH*B2+D*B2DH+B3*EDH)/(E*E)-2d0*D*B2*EDH/pow3(E)-B3DH/E
235 : C3DHH=(-B3DHH &
236 : + (DDHH*B2+2d0*DDH*B2DH+D*B2DHH+B3DH*EDH+B3*EDHH+B3DH*EDH)/E &
237 : - 2.0d0*((DDH*B2+D*B2DH+B3*EDH+DDH*B2+D*B2DH)*EDH+D*B2*EDHH)/(E*E) &
238 0 : + 6.0d0*D*B2*EDH*EDH/pow3(E))/E
239 0 : S3=C3*log(R3)
240 0 : S3DH=S3*C3DH/C3+C3*R3DH/R3
241 : S3DHH=(S3DH*C3DH+S3*C3DHH)/C3-S3*pow2(C3DH/C3) &
242 0 : + (C3DH*R3DH+C3*R3DHH)/R3-C3*pow2(R3DH/R3)
243 0 : S3DG=C3*R3DG/R3
244 0 : S3DGG=C3*(R3DGG/R3-pow2(R3DG/R3))
245 0 : S3DHG=(C3DH*R3DG+C3*R3DHG)/R3-C3*R3DG*R3DH/(R3*R3)
246 0 : B4=2.d0-D*D/E
247 0 : B4DH=EDH*(D/E)*(D/E)-2d0*D*DDH/E
248 : B4DHH=EDHH*(D/E)*(D/E)+2d0*EDH*(D/E)*(D/E)*(DDH/D-EDH/E) &
249 0 : - 2d0*(DDH*DDH+D*DDHH)/E+2d0*D*DDH*EDH/(E*E)
250 0 : C4=2d0*E*SQGE+D
251 0 : C4DH=2d0*EDH*SQGE+DDH
252 0 : C4DHH=2d0*EDHH*SQGE+DDHH
253 0 : C4DG=E/SQGE
254 0 : C4DGG=-0.5d0*E/(GAME*SQGE)
255 0 : C4DHG=EDH/SQGE
256 0 : S4A=2.0d0/E/DISCR
257 0 : S4AH=EDH/E+DIDH/DISCR
258 0 : S4ADH=-S4A*S4AH
259 0 : S4ADHH=-S4ADH*S4AH - S4A*(EDHH/E-(EDH/E)*(EDH/E)+DIDHH/DISCR-pow2(DIDH/DISCR))
260 0 : S4B=D*B3+B4*B2
261 0 : S4BDH=DDH*B3+D*B3DH+B4DH*B2+B4*B2DH
262 0 : S4BDHH=DDHH*B3+2d0*DDH*B3DH+D*B3DHH+B4DHH*B2+2d0*B4DH*B2DH+B4*B2DHH
263 0 : S4C=atan(C4/DISCR)-atan(D/DISCR)
264 0 : UP1=C4DH*DISCR-C4*DIDH
265 0 : DN1=DISCR*DISCR+C4*C4
266 0 : UP2=DDH*DISCR-D*DIDH
267 0 : DN2=DISCR*DISCR+D*D
268 0 : S4CDH=UP1/DN1-UP2/DN2
269 : S4CDHH=(C4DHH*DISCR-C4*DIDHH)/DN1 &
270 : - UP1*2d0*(DISCR*DIDH+C4*C4DH)/(DN1*DN1) &
271 0 : - (DDHH*DISCR-D*DIDHH)/DN2+UP2*2d0*(DISCR*DIDH+D*DDH)/(DN2*DN2)
272 0 : S4CDG=C4DG*DISCR/DN1
273 0 : S4CDGG=C4DGG*DISCR/DN1-2d0*C4*DISCR*pow2(C4DG/DN1)
274 0 : S4CDHG=(C4DHG*DISCR+C4DG*DIDH-C4DG*DISCR/DN1*2d0*(DISCR*DIDH+C4*C4DH))/DN1
275 0 : S4=S4A*S4B*S4C
276 0 : S4DH=S4ADH*S4B*S4C+S4A*S4BDH*S4C+S4A*S4B*S4CDH
277 : S4DHH=S4ADHH*S4B*S4C+S4A*S4BDHH*S4C+S4A*S4B*S4CDHH &
278 0 : + 2d0*(S4ADH*S4BDH*S4C+S4ADH*S4B*S4CDH+S4A*S4BDH*S4CDH)
279 0 : S4DG=S4A*S4B*S4CDG
280 0 : S4DGG=S4A*S4B*S4CDGG
281 0 : S4DHG=S4A*S4B*S4CDHG+S4CDG*(S4ADH*S4B+S4A*S4BDH)
282 0 : FXC=S1+S2+S3+S4
283 0 : FXCDH=S1DH+S2DH+S3DH+S4DH
284 0 : FXCDG=S1DG+S2DG+S3DG+S4DG
285 0 : FXCDHH=S1DHH+S2DHH+S3DHH+S4DHH
286 0 : FXCDGG=S2DGG+S3DGG+S4DGG
287 0 : FXCDHG=S1DHG+S2DHG+S3DHG+S4DHG
288 0 : PXC=(GAME*FXCDG-2d0*THETA*FXCDH)/3.d0
289 0 : UXC=GAME*FXCDG-THETA*FXCDH
290 0 : SXC=(GAME*S2DG-S2+GAME*S3DG-S3+S4A*S4B*(GAME*S4CDG-S4C))-THETA*FXCDH
291 0 : if (abs(SXC)<1.d-9*abs(THETA*FXCDH)) SXC=0.d0 ! accuracy loss
292 0 : CVXC=2d0*THETA*(GAME*FXCDHG-FXCDH)-THETA*THETA*FXCDHH-GAME*GAME*FXCDGG
293 0 : if (abs(CVXC)<1.d-9*abs(GAME*GAME*FXCDGG)) CVXC=0.d0 ! accuracy
294 0 : PDLH=THETA*(GAME*FXCDHG-2d0*FXCDH-2d0*THETA*FXCDHH)/3.d0
295 0 : PDLG=GAME*(FXCDG+GAME*FXCDGG-2d0*THETA*FXCDHG)/3.d0
296 0 : PDRXC=PXC+(PDLG-2d0*PDLH)/3.d0
297 0 : PDTXC=GAME*(THETA*FXCDHG-GAME*FXCDGG/3.d0)-THETA*(FXCDH/0.75d0+THETA*FXCDHH/1.5d0)
298 0 : return
299 : end subroutine EXCOR7
300 :
301 :
302 : end module create_EXCOR7_table
|
