Line data Source code
1 : /*----------------------------------------------------------------------------*/
2 : /* CP2K: A general program to perform molecular dynamics simulations */
3 : /* Copyright 2000-2025 CP2K developers group <https://cp2k.org> */
4 : /* */
5 : /* SPDX-License-Identifier: MIT */
6 : /*----------------------------------------------------------------------------*/
7 :
8 : /*
9 : * libgrpp - a library for the evaluation of integrals over
10 : * generalized relativistic pseudopotentials.
11 : *
12 : * Copyright (C) 2021-2023 Alexander Oleynichenko
13 : */
14 :
15 : #include "grpp_nuclear_models.h"
16 : #include "grpp_specfunc.h"
17 : #include <math.h>
18 :
19 : #ifndef M_PI
20 : #define M_PI 3.14159265358979323846
21 : #endif
22 :
23 : extern double libgrpp_fermi_model_Sk(int k, double x);
24 :
25 : extern double libgrpp_fermi_model_norm_factor(double c, double a);
26 :
27 : extern double libgrpp_fermi_bubble_model_norm_factor(double c, double a,
28 : double k);
29 :
30 : /*
31 : * estimates for root mean square radius of nuclear charge distribution
32 : */
33 :
34 : /**
35 : * estimate for R(rms) from:
36 : * W. R. Johnson, G. Soff. The Lamb shift in hydrogen-like atoms, 1 <= Z <=>
37 : * 110. At. Data Nucl. Data Tables, 33(3), 405 (1985).
38 : *
39 : * A = mass number
40 : * returns R(rms) in [fm] units
41 : */
42 0 : double libgrpp_estimate_nuclear_rms_radius_johnson_1985(int A) {
43 0 : return 0.836 * pow(A, 1.0 / 3.0) + 0.570;
44 : }
45 :
46 : /**
47 : * estimate for R(rms) from:
48 : * O. A. Golovko, I. A. Goidenko, I. I. Tupitsyn.
49 : * Quantum electrodynamic corrections for valence electrons in Eka-Hg.
50 : * Opt. Spectrosc. 104(5), 662 (2008).
51 : *
52 : * A = mass number
53 : * returns R(rms) in [fm] units
54 : */
55 0 : double libgrpp_estimate_nuclear_rms_radius_golovko_2008(int A) {
56 0 : return 0.77 * pow(A, 1.0 / 3.0) + 0.98;
57 : }
58 :
59 : /**
60 : * estimate parameters of the Fermi model using default formulas for 'c' and
61 : * 'a'. return -1 if such estimate cannot be performed, 1 otherwise.
62 : *
63 : * units: Fermi for R_rms, c, a.
64 : */
65 0 : int libgrpp_estimate_fermi_model_parameters(double R_rms, double *c,
66 : double *a) {
67 0 : const double t = 2.3;
68 0 : *a = t / (4 * log(3));
69 :
70 0 : const double c2 =
71 0 : 5.0 / 3.0 * R_rms * R_rms - 7.0 / 3.0 * M_PI * M_PI * (*a) * (*a);
72 0 : if (c2 < 0) {
73 : return -1;
74 : }
75 :
76 0 : *c = sqrt(c2);
77 0 : return 1;
78 : }
79 :
80 : /*
81 : * radially-local electrostatic potentials and density functions
82 : * for different nuclear finite charge distribution models
83 : */
84 :
85 : /**
86 : * point charge
87 : */
88 0 : double libgrpp_coulomb_potential_point(double r, double Z) { return -Z / r; }
89 :
90 : /**
91 : * uniformly charged ball: density
92 : */
93 0 : double libgrpp_charge_density_ball(double r, double Z, double R_rms) {
94 0 : const double R0 = sqrt(5.0 / 3.0) * R_rms;
95 :
96 0 : if (r <= R0) {
97 0 : return 3.0 * Z / (4.0 * M_PI * R0 * R0 * R0);
98 : } else {
99 : return 0.0;
100 : }
101 : }
102 :
103 : /**
104 : * uniformly charged ball: potential
105 : *
106 : * formula from:
107 : * L. Visscher, K. G. Dyall,
108 : * Dirac-Fock atomic electronic structure calculations using different nuclear
109 : * charge distributions. Atomic Data and Nuclear Data Tables, 67, 207–224 (1997)
110 : */
111 0 : double libgrpp_coulomb_potential_ball(double r, double Z, double R_rms) {
112 0 : const double R0 = sqrt(5.0 / 3.0) * R_rms;
113 :
114 0 : if (r <= R0) {
115 0 : return -Z / (2.0 * R0) * (3.0 - r * r / (R0 * R0));
116 : } else {
117 0 : return -Z / r;
118 : }
119 : }
120 :
121 : /**
122 : * Gaussian charge distribution: density
123 : */
124 0 : double libgrpp_charge_density_gaussian(double r, double Z, double R_rms) {
125 0 : const double xi = 3.0 / (2.0 * R_rms * R_rms);
126 0 : const double rho0 = Z * pow(xi / M_PI, 1.5);
127 :
128 0 : return rho0 * exp(-xi * r * r);
129 : }
130 :
131 : /**
132 : * Gaussian charge distribution: potential
133 : *
134 : * formula from:
135 : * L. Visscher, K. G. Dyall,
136 : * Dirac-Fock atomic electronic structure calculations using different nuclear
137 : * charge distributions. Atomic Data and Nuclear Data Tables, 67, 207–224 (1997)
138 : */
139 0 : double libgrpp_coulomb_potential_gaussian(double r, double Z, double R_rms) {
140 0 : const double xi = 3.0 / (2.0 * R_rms * R_rms);
141 :
142 0 : return -Z / r * erf(sqrt(xi) * r);
143 : }
144 :
145 : /**
146 : * Fermi charge distribution: density
147 : */
148 0 : double libgrpp_charge_density_fermi(double r, double Z, double c, double a) {
149 0 : const double N = libgrpp_fermi_model_norm_factor(c, a);
150 0 : const double c3 = c * c * c;
151 0 : const double rho0 = 3.0 * Z / (4 * M_PI * c3 * N);
152 :
153 0 : return rho0 / (1.0 + exp((r - c) / a));
154 : }
155 :
156 : /**
157 : * Fermi charge distribution: rms radius
158 : */
159 0 : double libgrpp_rms_radius_fermi(int Z, double c, double a) {
160 0 : const double N = libgrpp_fermi_model_norm_factor(c, a);
161 0 : const double c3 = c * c * c;
162 0 : const double rho0 = 3.0 * Z / (4 * M_PI * c3 * N);
163 :
164 0 : const double r2 =
165 0 : 4 * M_PI * rho0 / Z * pow(c, 5) / 5.0 *
166 0 : (1.0 + 10.0 / 3.0 * a * a * M_PI * M_PI / (c * c) +
167 0 : 7.0 / 3.0 * pow(M_PI, 4) * pow(a, 4) / pow(c, 4) -
168 0 : 120.0 * pow(a, 5) / pow(c, 5) * libgrpp_specfunc_fermi_sk(5, -c / a));
169 :
170 0 : return sqrt(r2);
171 : }
172 :
173 : /**
174 : * Fermi charge distribution: potential
175 : *
176 : * formula from:
177 : * N. S. Mosyagin, A. V. Zaitsevskii, A. V. Titov
178 : * Generalized relativistic effective core potentials for superheavy elements
179 : * Int. J. Quantum Chem. e26076 (2019)
180 : * doi: https://doi.org/10.1002/qua.26076
181 : */
182 0 : double libgrpp_coulomb_potential_fermi(double r, double Z, double c, double a) {
183 0 : const double a2 = a * a;
184 0 : const double a3 = a * a * a;
185 0 : const double c3 = c * c * c;
186 0 : const double N = libgrpp_fermi_model_norm_factor(c, a);
187 :
188 0 : if (r > c) {
189 0 : const double S2 = libgrpp_specfunc_fermi_sk(2, (c - r) / a);
190 0 : const double S3 = libgrpp_specfunc_fermi_sk(3, (c - r) / a);
191 :
192 0 : return -Z / (N * r) * (N + 3 * a2 * r / c3 * S2 + 6 * a3 / c3 * S3);
193 : } else {
194 0 : const double P2 = libgrpp_specfunc_fermi_sk(2, (r - c) / a);
195 0 : const double P3 = libgrpp_specfunc_fermi_sk(3, (r - c) / a);
196 0 : const double S3 = libgrpp_specfunc_fermi_sk(3, -c / a);
197 0 : const double r3 = r * r * r;
198 :
199 0 : return -Z / (N * r) *
200 0 : (1.5 * r / c - r3 / (2 * c3) + M_PI * M_PI * a2 * r / (2 * c3) -
201 0 : 3 * a2 * r / c3 * P2 + 6 * a3 / c3 * (P3 - S3));
202 : }
203 : }
204 :
205 : /**
206 : * normalization factor for the Fermi nuclear charge distribution
207 : */
208 0 : double libgrpp_fermi_model_norm_factor(double c, double a) {
209 0 : const double a2 = a * a;
210 0 : const double a3 = a * a * a;
211 0 : const double c2 = c * c;
212 0 : const double c3 = c * c * c;
213 :
214 0 : return 1.0 + M_PI * M_PI * a2 / c2 -
215 0 : 6.0 * a3 / c3 * libgrpp_specfunc_fermi_sk(3, -c / a);
216 : }
217 :
218 : /**
219 : * "Fermi bubble" charge distribution: density
220 : */
221 0 : double libgrpp_charge_density_fermi_bubble(double r, double Z, double c,
222 : double a, double k) {
223 0 : const double Nk = libgrpp_fermi_bubble_model_norm_factor(c, a, k);
224 0 : const double c3 = c * c * c;
225 0 : const double rho0 = 3.0 * Z / (4 * M_PI * c3 * Nk);
226 :
227 0 : return rho0 * (1 + k * pow(r / c, 2)) / (1.0 + exp((r - c) / a));
228 : }
229 :
230 : /**
231 : * "Fermi bubble" charge distribution: rms radius
232 : */
233 0 : double libgrpp_rms_radius_fermi_bubble(int Z, double c, double a, double k) {
234 0 : const double Nk = libgrpp_fermi_bubble_model_norm_factor(c, a, k);
235 0 : const double c3 = c * c * c;
236 0 : const double rho0 = 3.0 * Z / (4 * M_PI * c3 * Nk);
237 :
238 0 : const double part_r4 =
239 0 : pow(c, 5) / 5.0 *
240 0 : (1.0 + 10.0 / 3.0 * a * a * M_PI * M_PI / (c * c) +
241 0 : 7.0 / 3.0 * pow(M_PI, 4) * pow(a, 4) / pow(c, 4) -
242 0 : 120.0 * pow(a, 5) / pow(c, 5) * libgrpp_specfunc_fermi_sk(5, -c / a));
243 :
244 0 : const double part_r6 =
245 0 : pow(c, 7) / 7.0 *
246 0 : (1.0 + 7.0 * a * a * M_PI * M_PI / (c * c) +
247 0 : 49.0 / 3.0 * pow(M_PI, 4) * pow(a, 4) / pow(c, 4) +
248 0 : 31.0 / 3.0 * pow(M_PI, 6) * pow(a, 6) / pow(c, 6) -
249 0 : 5040.0 * pow(a, 7) / pow(c, 7) * libgrpp_specfunc_fermi_sk(7, -c / a));
250 :
251 0 : const double r2 = 4 * M_PI * rho0 / Z * (part_r4 + k / (c * c) * part_r6);
252 :
253 0 : return sqrt(r2);
254 : }
255 :
256 : /**
257 : * "Fermi bubble" charge distribution: potential.
258 : *
259 : * derivation of the formula is based on:
260 : *
261 : * F. A. Parpia and A. K. Mohanty,
262 : * Relativistic basis-set calculations for atoms with Fermi nuclei.
263 : * Phys. Rev. A 46, 3735 (1992)
264 : * doi: 10.1103/PhysRevA.46.3735
265 : *
266 : */
267 0 : double libgrpp_coulomb_potential_fermi_bubble(double r, double Z, double c,
268 : double a, double k) {
269 : // const double a2 = a * a;
270 : // const double a3 = a * a * a;
271 : // const double c2 = c * c;
272 : // const double c3 = c * c * c;
273 :
274 0 : const double Nk = libgrpp_fermi_bubble_model_norm_factor(c, a, k);
275 : // const double rho0 = 3 * Z / (4 * M_PI * M_PI * c * c * c * Nk);
276 :
277 0 : double F0 = 0.0;
278 0 : double F2 = 0.0;
279 :
280 0 : if (r < c) {
281 0 : const double S2 = libgrpp_specfunc_fermi_sk(2, (r - c) / a);
282 0 : const double S3 = libgrpp_specfunc_fermi_sk(3, (r - c) / a);
283 0 : const double S4 = libgrpp_specfunc_fermi_sk(4, (r - c) / a);
284 0 : const double S5 = libgrpp_specfunc_fermi_sk(5, (r - c) / a);
285 :
286 : // contribution from the "classical" Fermi term
287 0 : F0 = -pow(r, 3) / 6.0 - r * a * a * S2 + 2.0 * pow(a, 3) * S3 +
288 0 : r * c * c / 2.0 + M_PI * M_PI / 6.0 * r * a * a -
289 0 : 2.0 * pow(a, 3) * libgrpp_specfunc_fermi_sk(3, -c / a);
290 :
291 : // contribution from the quadratic, "hole" term
292 0 : F2 = -pow(r, 5) / 20.0 - pow(r, 3) * pow(a, 2) * S2 +
293 0 : 6.0 * pow(r, 2) * pow(a, 3) * S3 - 18.0 * r * pow(a, 4) * S4 +
294 0 : 24.0 * pow(a, 5) * S5 + r * pow(c, 4) / 4.0 +
295 0 : r * M_PI * M_PI * c * c * a * a / 2.0 +
296 0 : r * pow(a, 4) * pow(M_PI, 4) * 7.0 / 60.0 -
297 0 : 24.0 * pow(a, 5) * libgrpp_specfunc_fermi_sk(5, -c / a);
298 : } else {
299 0 : const double S2 = libgrpp_specfunc_fermi_sk(2, (c - r) / a);
300 0 : const double S3 = libgrpp_specfunc_fermi_sk(3, (c - r) / a);
301 0 : const double S4 = libgrpp_specfunc_fermi_sk(4, (c - r) / a);
302 0 : const double S5 = libgrpp_specfunc_fermi_sk(5, (c - r) / a);
303 :
304 : // contribution from the "classical" Fermi term
305 0 : F0 = pow(c, 3) / 3.0 + M_PI * M_PI / 3.0 * c * a * a -
306 0 : 2.0 * pow(a, 3) * libgrpp_specfunc_fermi_sk(3, -c / a) +
307 0 : r * a * a * S2 + 2.0 * pow(a, 3) * S3;
308 :
309 : // contribution from the quadratic, "hole" term
310 0 : F2 = pow(c, 5) / 5.0 + 2.0 * pow(c, 3) * a * a * M_PI * M_PI / 3.0 +
311 0 : 7.0 * pow(a, 4) * c * pow(M_PI, 4) / 15.0 -
312 0 : 24.0 * pow(a, 5) * libgrpp_specfunc_fermi_sk(5, -c / a) +
313 0 : pow(a, 2) * pow(r, 3) * S2 + 6.0 * pow(a, 3) * pow(r, 2) * S3 +
314 0 : 18.0 * r * pow(a, 4) * S4 + 24.0 * pow(a, 5) * S5;
315 : }
316 :
317 0 : return -Z / (Nk * r) * 3.0 / pow(c, 3) * (F0 + k / (c * c) * F2);
318 : }
319 :
320 : /**
321 : * normalization factor for the "Fermi bubble" nuclear charge distribution
322 : */
323 0 : double libgrpp_fermi_bubble_model_norm_factor(double c, double a, double k) {
324 0 : const double a2 = a * a;
325 0 : const double a3 = a * a2;
326 0 : const double a4 = a * a3;
327 0 : const double a5 = a * a4;
328 0 : const double c2 = c * c;
329 0 : const double c3 = c * c2;
330 0 : const double c4 = c * c3;
331 0 : const double c5 = c * c4;
332 :
333 0 : return libgrpp_fermi_model_norm_factor(c, a) + 3.0 / 5.0 * k +
334 0 : 2.0 * M_PI * M_PI * a2 * k / c2 +
335 0 : 7.0 * M_PI * M_PI * M_PI * M_PI * a4 * k / (5.0 * c4) -
336 0 : 72.0 * a5 * k / c5 * libgrpp_specfunc_fermi_sk(5, -c / a);
337 : }
|
