12int mspes_SAC(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
13 const double V0 = 0.01e0, V1 = 0.005e0, E0 = 0.0e0, a = 1.6e0, b = 1.0e0;
15 V[0] = (R[0] > 0) ? V0 * (1.0 - exp(-a * R[0])) : -V0 * (1.0 - exp(a * R[0]));
17 V[1] = V1 * exp(-b * R[0] * R[0]);
19 if (flag < 1)
return 0;
21 dV[0] = (R[0] > 0) ? V0 * a * exp(-a * R[0]) : V0 * a * exp(a * R[0]);
23 dV[1] = -2 * b * R[0] * V[1];
25 if (flag < 2)
return 0;
32int mspes_SAC2(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
33 const double V0 = 0.01e0, V1 = 0.005e0, E0 = 0.0e0, a = 1.6e0, b = 1.0e0;
35 V[0] = V0 * tanh(a * R[0]);
37 V[1] = V1 * exp(-b * R[0] * R[0]);
39 if (flag < 1)
return 0;
41 double tmp = cosh(a * R[0]);
42 dV[0] = V0 * a / (tmp * tmp);
44 dV[1] = -2 * b * R[0] * V[1];
46 if (flag < 2)
return 0;
53int mspes_SAC3(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
54 const double V1 = 0.04e0, V2 = 0.01e0, Vc = 0.005e0, a = 1.0e0, b = 1.0e0, Rc = 0.7;
56 V[0] = V1 * (1.0e0 + tanh(a * R[0]));
57 V[3] = V2 * (1.0e0 - tanh(a * R[0]));
58 V[1] = Vc * exp(-b * (R[0] + Rc) * (R[0] + Rc));
60 if (flag < 1)
return 0;
62 double tmp = cosh(a * R[0]);
63 dV[0] = +V1 * a / (tmp * tmp);
64 dV[3] = -V2 * a / (tmp * tmp);
65 dV[1] = -2 * b * (R[0] + Rc) * V[1];
67 if (flag < 2)
return 0;
75int mspes_DAC(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
76 const double V0 = 0.10e0, V1 = 0.015e0, E0 = 0.05e0, a = 0.28e0, b = 0.06e0;
79 V[3] = -V0 * exp(-a * R[0] * R[0]) + E0;
80 V[1] = V1 * exp(-b * R[0] * R[0]);
82 if (flag < 1)
return 0;
85 dV[3] = 2 * a * R[0] * V0 * exp(-a * R[0] * R[0]);
86 dV[1] = -2 * b * R[0] * V[1];
88 if (flag < 2)
return 0;
95int mspes_ECR(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
96 const double V0 = 6.0e-4, V1 = 0.1e0, E0 = 0.0e0, a = 0.0e0, b = 0.9e0;
98 V[0] = -V0, V[3] = V0;
99 V[1] = (R[0] < 0) ? V1 * exp(b * R[0]) : V1 * (2.0 - exp(-b * R[0]));
101 if (flag < 1)
return 0;
103 dV[0] = 0.0e0, dV[3] = 0.0e0;
104 dV[1] = (R[0] < 0) ? V1 * b * exp(b * R[0]) : V1 * b * exp(-b * R[0]);
106 if (flag < 2)
return 0;
113int mspes_DBG(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
114 const double V0 = 6.0e-4, V1 = 0.1e0, E0 = 0.0e0, a = 0.0e0, b = 0.9e0, Z = 10.0e0;
116 V[0] = -V0, V[3] = V0;
118 V[1] = V1 * exp(b * (R[0] - Z)) + V1 * (2.0 - exp(b * (R[0] - Z)));
119 }
else if (R[0] < Z) {
120 V[1] = V1 * exp(b * (R[0] - Z)) + V1 * exp(-b * (R[0] + Z));
122 V[1] = V1 * exp(-b * (R[0] + Z)) + V1 * (2.0 - exp(-b * (R[0] - Z)));
125 if (flag < 1)
return 0;
127 dV[0] = 0.0e0, dV[3] = 0.0e0;
129 dV[1] = V1 * exp(b * (R[0] - Z)) * b + V1 * (-exp(b * (R[0] - Z))) * b;
130 }
else if (R[0] < Z) {
131 dV[1] = V1 * exp(b * (R[0] - Z)) * b + V1 * (exp(-b * (R[0] + Z))) * (-b);
133 dV[1] = V1 * exp(-b * (R[0] + Z)) * (-b) + V1 * (-exp(-b * (R[0] - Z))) * (-b);
136 if (flag < 2)
return 0;
143int mspes_DAG(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
144 const double V0 = 6.0e-4, V1 = 0.1e0, E0 = 0.0e0, a = 0.0e0, b = 0.9e0, Z = 10.0e0;
146 V[0] = -V0, V[3] = V0;
148 V[1] = -V1 * exp(b * (R[0] - Z)) + V1 * exp(b * (R[0] + Z));
149 }
else if (R[0] < Z) {
150 V[1] = -V1 * exp(b * (R[0] - Z)) - V1 * exp(-b * (R[0] + Z)) + 2 * V1;
152 V[1] = V1 * exp(-b * (R[0] - Z)) - V1 * exp(-b * (R[0] + Z));
155 if (flag < 1)
return 0;
157 dV[0] = 0.0e0, dV[3] = 0.0e0;
159 dV[1] = -V1 * exp(b * (R[0] - Z)) * b + V1 * exp(b * (R[0] + Z)) * b;
160 }
else if (R[0] < Z) {
161 dV[1] = -V1 * exp(b * (R[0] - Z)) * b - V1 * exp(-b * (R[0] + Z)) * (-b);
163 dV[1] = V1 * exp(-b * (R[0] - Z)) * (-b) - V1 * exp(-b * (R[0] + Z)) * (-b);
166 if (flag < 2)
return 0;
173int mspes_DRN(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
174 const double E0 = 0.01e0, V1 = 0.03e0, b = 3.2e0, Z = 2.0e0;
178 V[1] = V1 * (exp(-b * (R[0] - Z) * (R[0] - Z)) + exp(-b * (R[0] + Z) * (R[0] + Z)));
180 if (flag < 1)
return 0;
184 dV[1] = V1 * (exp(-b * (R[0] - Z) * (R[0] - Z)) * (-2 * b * (R[0] - Z)) +
185 exp(-b * (R[0] + Z) * (R[0] + Z)) * (-2 * b * (R[0] + Z)));
187 if (flag < 2)
return 0;
195int mspes_MORSE3A(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
197 const double De[3] = {0.003e0, 0.004e0, 0.003e0};
198 const double Be[3] = {0.65e0, 0.60e0, 0.65e0};
199 const double Re[3] = {5.0e0, 4.0e0, 6.0e0};
200 const double C[3] = {0.0e0, 0.01e0, 0.006e0};
201 const double Aijx[3] = {0.002e0, 0.00e0, 0.002e0};
202 const double Aeijx[3] = {16.0e0, 0.00e0, 16.0e0};
203 const double Rijx[3] = {4.80e0, 0.00e0, 3.40e0};
205 V[0] = De[0] * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) + C[0];
206 V[4] = De[1] * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) + C[1];
207 V[8] = De[2] * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) + C[2];
209 V[1] = Aijx[2] * exp(-Aeijx[2] * (R[0] - Rijx[2]) * (R[0] - Rijx[2]));
210 V[2] = Aijx[1] * exp(-Aeijx[1] * (R[0] - Rijx[1]) * (R[0] - Rijx[1]));
211 V[5] = Aijx[0] * exp(-Aeijx[0] * (R[0] - Rijx[0]) * (R[0] - Rijx[0]));
216 if (flag < 1)
return 0;
218 dV[0] = De[0] * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) * exp(-Be[0] * (R[0] - Re[0])) * 2 * Be[0];
219 dV[4] = De[1] * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) * exp(-Be[1] * (R[0] - Re[1])) * 2 * Be[1];
220 dV[8] = De[2] * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) * exp(-Be[2] * (R[0] - Re[2])) * 2 * Be[2];
222 dV[1] = V[1] * (-2 * Aeijx[2] * (R[0] - Rijx[2]));
223 dV[2] = V[2] * (-2 * Aeijx[1] * (R[0] - Rijx[1]));
224 dV[5] = V[5] * (-2 * Aeijx[0] * (R[0] - Rijx[0]));
229 if (flag < 2)
return 0;
234int mspes_MORSE3B(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
236 const double De[3] = {0.020e0, 0.010e0, 0.003e0};
237 const double Be[3] = {0.65e0, 0.40e0, 0.65e0};
238 const double Re[3] = {4.5e0, 4.0e0, 4.4e0};
239 const double C[3] = {0.0e0, 0.01e0, 0.02e0};
240 const double Aijx[3] = {0.00e0, 0.005e0, 0.005e0};
241 const double Aeijx[3] = {0.00e0, 32.0e0, 32.0e0};
242 const double Rijx[3] = {0.00e0, 3.34e0, 3.66e0};
244 V[0] = De[0] * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) + C[0];
245 V[4] = De[1] * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) + C[1];
246 V[8] = De[2] * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) + C[2];
248 V[1] = Aijx[2] * exp(-Aeijx[2] * (R[0] - Rijx[2]) * (R[0] - Rijx[2]));
249 V[2] = Aijx[1] * exp(-Aeijx[1] * (R[0] - Rijx[1]) * (R[0] - Rijx[1]));
250 V[5] = Aijx[0] * exp(-Aeijx[0] * (R[0] - Rijx[0]) * (R[0] - Rijx[0]));
255 if (flag < 1)
return 0;
257 dV[0] = De[0] * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) * exp(-Be[0] * (R[0] - Re[0])) * 2 * Be[0];
258 dV[4] = De[1] * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) * exp(-Be[1] * (R[0] - Re[1])) * 2 * Be[1];
259 dV[8] = De[2] * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) * exp(-Be[2] * (R[0] - Re[2])) * 2 * Be[2];
261 dV[1] = V[1] * (-2 * Aeijx[2] * (R[0] - Rijx[2]));
262 dV[2] = V[2] * (-2 * Aeijx[1] * (R[0] - Rijx[1]));
263 dV[5] = V[5] * (-2 * Aeijx[0] * (R[0] - Rijx[0]));
268 if (flag < 2)
return 0;
273int mspes_MORSE3C(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
275 const double De[3] = {0.020e0, 0.020e0, 0.003e0};
276 const double Be[3] = {0.40e0, 0.65e0, 0.65e0};
277 const double Re[3] = {4.0e0, 4.5e0, 6.0e0};
278 const double C[3] = {0.02e0, 0.00e0, 0.02e0};
279 const double Aijx[3] = {0.00e0, 0.005e0, 0.005e0};
280 const double Aeijx[3] = {0.00e0, 32.0e0, 32.0e0};
281 const double Rijx[3] = {0.00e0, 4.97e0, 3.40e0};
283 V[0] = De[0] * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) + C[0];
284 V[4] = De[1] * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) + C[1];
285 V[8] = De[2] * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) + C[2];
287 V[1] = Aijx[2] * exp(-Aeijx[2] * (R[0] - Rijx[2]) * (R[0] - Rijx[2]));
288 V[2] = Aijx[1] * exp(-Aeijx[1] * (R[0] - Rijx[1]) * (R[0] - Rijx[1]));
289 V[5] = Aijx[0] * exp(-Aeijx[0] * (R[0] - Rijx[0]) * (R[0] - Rijx[0]));
296 if (flag < 1)
return 0;
298 dV[0] = De[0] * (1.0e0 - exp(-Be[0] * (R[0] - Re[0]))) * exp(-Be[0] * (R[0] - Re[0])) * 2 * Be[0];
299 dV[4] = De[1] * (1.0e0 - exp(-Be[1] * (R[0] - Re[1]))) * exp(-Be[1] * (R[0] - Re[1])) * 2 * Be[1];
300 dV[8] = De[2] * (1.0e0 - exp(-Be[2] * (R[0] - Re[2]))) * exp(-Be[2] * (R[0] - Re[2])) * 2 * Be[2];
302 dV[1] = V[1] * (-2 * Aeijx[2] * (R[0] - Rijx[2]));
303 dV[2] = V[2] * (-2 * Aeijx[1] * (R[0] - Rijx[1]));
304 dV[5] = V[5] * (-2 * Aeijx[0] * (R[0] - Rijx[0]));
309 if (flag < 2)
return 0;
315int mspes_MORSE15(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
317 const double Dg = 0.2e0, De = 0.05e0, Dc = 0.05e0, alpha = 0.4e0, lambda = 1.0e0, eta = 0.004e0;
318 const double Re[15] = {3.0000000000000e0, 6.8335793401926175e0, 6.909940519903887e0, 6.988372712202933e0,
319 7.0690037103879675e0, 7.151973244334301e0, 7.237434522736795e0, 7.325556032471846e0,
320 7.416523648311893e0, 7.5105431195440095e0, 7.607843017311827e0, 7.708678249086644e0,
321 7.813334276495011e0, 7.922132212503321e0, 8.035435027584366};
322 for (
int i = 0; i < 225; ++i) V[i] = 0.0e0;
323 V[0] = Dg * (1 - exp(-alpha * (R[0] - Re[0]))) * (1 - exp(-alpha * (R[0] - Re[0])));
324 for (
int i = 1; i < 15; ++i) {
325 V[i * 16] = exp(-alpha * R[0]) + eta * (i + 1) + De;
326 V[i] = Dc * exp(-lambda * (R[0] - Re[i]) * (R[0] - Re[i]));
329 if (flag < 1)
return 0;
331 for (
int i = 0; i < 225; ++i) dV[i] = 0.0e0;
332 dV[0] = Dg * (1 - exp(-alpha * (R[0] - Re[0]))) * exp(-alpha * (R[0] - Re[0])) * 2 * alpha;
333 for (
int i = 1; i < 15; ++i) {
334 dV[i * 16] = -alpha * exp(-alpha * R[0]);
335 dV[i] = V[i] * (-2 * lambda * (R[0] - Re[i]));
339 if (flag < 2)
return 0;
344int mspes_MORSE15C(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
346 const double Dg = 0.2e0, De = 0.05e0, Dc = 0.05e0, alpha = 0.4e0, lambda = 1.0e0, eta = 0.004e0;
347 const double Re[15] = {3.0000000000000e0, 6.8335793401926175e0, 6.909940519903887e0, 6.988372712202933e0,
348 7.0690037103879675e0, 7.151973244334301e0, 7.237434522736795e0, 7.325556032471846e0,
349 7.416523648311893e0, 7.5105431195440095e0, 7.607843017311827e0, 7.708678249086644e0,
350 7.813334276495011e0, 7.922132212503321e0, 8.035435027584366};
351 for (
int i = 0; i < 225; ++i) V[i] = 0.0e0;
352 V[0] = Dg * (1 - exp(-alpha * (R[0] - Re[0]))) * (1 - exp(-alpha * (R[0] - Re[0])));
353 for (
int i = 1; i < 15; ++i) {
354 V[i * 16] = exp(-alpha * R[0]) + eta * (i + 1) + De;
358 if (flag < 1)
return 0;
360 for (
int i = 0; i < 225; ++i) dV[i] = 0.0e0;
361 dV[0] = Dg * (1 - exp(-alpha * (R[0] - Re[0]))) * exp(-alpha * (R[0] - Re[0])) * 2 * alpha;
362 for (
int i = 1; i < 15; ++i) {
363 dV[i * 16] = -alpha * exp(-alpha * R[0]);
368 if (flag < 2)
return 0;
373int mspes_MORSE15E(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
375 const double Dg = 0.2e0, De = 0.05e0, Dc = 0.05e0, alpha = 0.4e0, lambda = 1.0e0, eta = 0.004e0;
376 const double Re[15] = {3.0000000000000e0, 6.8335793401926175e0, 6.909940519903887e0, 6.988372712202933e0,
377 7.0690037103879675e0, 7.151973244334301e0, 7.237434522736795e0, 7.325556032471846e0,
378 7.416523648311893e0, 7.5105431195440095e0, 7.607843017311827e0, 7.708678249086644e0,
379 7.813334276495011e0, 7.922132212503321e0, 8.035435027584366};
380 for (
int i = 0; i < 225; ++i) V[i] = 0.0e0;
381 V[0] = Dg * (1 - exp(-alpha * (R[0] - Re[0]))) * (1 - exp(-alpha * (R[0] - Re[0])));
382 for (
int i = 1; i < 15; ++i) {
383 V[i * 16] = exp(-alpha * R[0]) + eta * (i + 1) + De;
384 V[i] = Dc * exp(-0.2 * R[0]);
387 if (flag < 1)
return 0;
389 for (
int i = 0; i < 225; ++i) dV[i] = 0.0e0;
390 dV[0] = Dg * (1 - exp(-alpha * (R[0] - Re[0]))) * exp(-alpha * (R[0] - Re[0])) * 2 * alpha;
391 for (
int i = 1; i < 15; ++i) {
392 dV[i * 16] = -alpha * exp(-alpha * R[0]);
397 if (flag < 2)
return 0;
402int mspes_CL1D(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
409 if (flag < 1)
return 0;
415 if (flag < 2)
return 0;
418 for (
int i = 0; i < 4; ++i) ddV[i] = 0;
422int mspes_JC1D(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
429 if (flag < 1)
return 0;
435 if (flag < 2)
return 0;
442int mspes_NA_I(
double* V,
double* dV,
double* ddV,
double* R,
int flag,
int rdim,
int fdim) {
443 const double au_2_ev = 27.21138602e0;
444 const double au_2_ang = 0.529177249e0;
445 const double Acov = 3150.0e0 / au_2_ev;
446 const double Bcov12 = pow(2.647e0 / au_2_ang, 12) / au_2_ev;
447 const double Ccov = 1000.0e0 / au_2_ev / pow(au_2_ang, 6);
448 const double Rcov = 0.435e0 / au_2_ang;
449 const double Aion = 2760.0e0 / au_2_ev;
450 const double Bion8 = pow(2.398e0 / au_2_ang, 8) / au_2_ev;
451 const double Cion = 11.3e0 / au_2_ev / pow(au_2_ang, 6);
452 const double Rion = 0.3489e0 / au_2_ang;
453 const double al_Mp = 0.408e0 / pow(au_2_ang, 3);
454 const double al_Xm = 6.431e0 / pow(au_2_ang, 3);
455 const double al_M = 27.0e0 / pow(au_2_ang, 3);
456 const double al_X = 7.0e0 / pow(au_2_ang, 3);
457 const double Eth = 2.075e0 / au_2_ev;
458 const double Aoff = 17.08e0 / au_2_ev;
459 const double Roff = 1.239e0 / au_2_ang;
460 const double mred = (22.989769282e0 * 126.904473e0) / (22.989769282e0 + 126.904473e0) /
phys::au_2_amu;
470 double r12 = r6 * r6;
471 double r13 = r12 * r;
473 double Vcent =
mspes_parm[1] / (2 * mred * r2);
476 + (Acov + Bcov12 / r12) * exp(-r / Rcov)
479 + (Aion + Bion8 / r8) * exp(-r / Rion)
481 - 0.5 * (al_Mp + al_Xm) / r4
483 - 2 * al_Mp * al_Xm / r7
485 V[1] = Aoff * exp(-r / Roff);
488 if (flag < 1)
return 0;
490 double dVcent = -2 *
mspes_parm[1] / (2 * mred * r2 * r);
493 + (-12 * Bcov12 / r13) * exp(-r / Rcov)
494 + (Acov + Bcov12 / r12) * (-1 / Rcov) * exp(-r / Rcov)
497 + (-8 * Bion8 / r9) * exp(-r / Rion)
498 + (Aion + Bion8 / r8) * (-1 / Rion) * exp(-r / Rion)
500 + 2 * (al_Mp + al_Xm) / r5
502 + 14 * al_Mp * al_Xm / r8;
503 dV[1] = (-1 / Roff) * Aoff * exp(-r / Roff);
506 if (flag < 2)
return 0;
522 case NAD1DPolicy::SAC:
523 case NAD1DPolicy::SAC2:
524 case NAD1DPolicy::SAC3:
525 case NAD1DPolicy::DAC:
526 case NAD1DPolicy::ECR:
527 case NAD1DPolicy::DBG:
528 case NAD1DPolicy::DAG:
529 case NAD1DPolicy::DRN: {
530 double x0_read =
_param->get_double(
"x0",
LOC(), -10.0e0);
531 double p0_read =
_param->get_double(
"p0",
LOC(), 10.0e0);
532 double mass_read =
_param->get_double(
"m0",
LOC(), 2000.0e0);
533 double tend_read =
_param->get_double(
"tend",
LOC(), -1);
534 double dt_read =
_param->get_double(
"dt",
LOC(), -1);
536 double tend_rev = std::abs(5 * x0_read * mass_read / p0_read);
537 double dt_rev = std::abs(x0_read * mass_read / p0_read) / 5000;
538 if (tend_read < 0) (*(
_param->pjson()))[
"tend"] = tend_rev;
539 if (dt_read < 0) (*(
_param->pjson()))[
"dt"] = dt_rev;
542 case NAD1DPolicy::NA_I: {
547 double mred = (mass_Na * mass_I) / (mass_Na + mass_I);
549 double parm_E =
_param->get_double(
"parm_E",
LOC(), 1.0e0);
550 double tend_read =
_param->get_double(
"tend",
LOC(), -1);
551 double dt_read =
_param->get_double(
"dt",
LOC(), -1);
554 double vel = sqrt(2 * parm_E / mred);
555 double tend_rev = 2 * x0_max / vel;
556 double dt_rev = tend_rev / 50000;
558 if (tend_read < 0) (*(
_param->pjson()))[
"tend"] = tend_rev;
559 if (dt_read < 0) (*(
_param->pjson()))[
"dt"] = dt_rev;
570 std::ifstream ifs(
"Hsys.dat");
571 std::string H_unit_str;
572 std::string firstline;
573 getline(ifs, firstline);
574 std::stringstream sstr(firstline);
580 if (ifs >> val)
Hsys[i] = val / H_unit;
603 double x0_read =
_param->get_double(
"x0grid",
LOC(), -10.0e0);
604 int Nxgird =
_param->get_int(
"Nxgrid",
LOC(), 101);
605 double dx = (2 * fabs(x0_read)) / (Nxgird - 1);
606 double* xgrid = DS->def_real(
"integrator.xgrid", Nxgird);
607 for (
int i = 0; i < Nxgird; ++i) xgrid[i] = -abs(x0_read) + i * dx;
615 double varx =
_param->get_double(
"varx",
LOC(), 0.5e0);
616 double varp =
_param->get_double(
"varp",
LOC(), 0.5e0);
621 case NAD1DPolicy::SAC:
622 case NAD1DPolicy::SAC2:
623 case NAD1DPolicy::DAC:
624 case NAD1DPolicy::ECR:
625 case NAD1DPolicy::DBG:
626 case NAD1DPolicy::DAG:
627 case NAD1DPolicy::DRN: {
631 case NAD1DPolicy::SAC3: {
633 double gammawidth = 0.25e0;
634 x_sigma[0] = sqrt(0.5e0 / gammawidth);
635 p_sigma[0] = sqrt(0.5e0 * gammawidth);
638 case NAD1DPolicy::MORSE3A:
639 case NAD1DPolicy::MORSE3B:
640 case NAD1DPolicy::MORSE3C: {
644 case NAD1DPolicy::MORSE3A:
647 case NAD1DPolicy::MORSE3B:
650 case NAD1DPolicy::MORSE3C:
655 double wground = 5.0e-03;
660 case NAD1DPolicy::MORSE15:
661 case NAD1DPolicy::MORSE15C:
662 case NAD1DPolicy::MORSE15E: {
667 double Dg = 0.2e0, alpha = 0.4e0;
668 x_sigma[0] = sqrt(0.5e0 / std::sqrt(
mass[0] * 2 * alpha * alpha * Dg));
672 case NAD1DPolicy::CL1D:
673 case NAD1DPolicy::JC1D: {
682 case NAD1DPolicy::NA_I: {
687 mass[0] = (mass_Na * mass_I) / (mass_Na + mass_I);
699 case NAD1DPolicy::NA_I: {
703 double b = sqrt(randu) * bmax;
735 case NAD1DPolicy::SAC:
738 case NAD1DPolicy::SAC2:
741 case NAD1DPolicy::SAC3:
744 case NAD1DPolicy::DAC:
747 case NAD1DPolicy::ECR:
750 case NAD1DPolicy::DBG:
753 case NAD1DPolicy::DAG:
756 case NAD1DPolicy::DRN:
759 case NAD1DPolicy::MORSE3A:
762 case NAD1DPolicy::MORSE3B:
765 case NAD1DPolicy::MORSE3C:
768 case NAD1DPolicy::MORSE15:
771 case NAD1DPolicy::MORSE15C:
774 case NAD1DPolicy::MORSE15E:
777 case NAD1DPolicy::JC1D:
780 case NAD1DPolicy::CL1D:
783 case NAD1DPolicy::NA_I:
786 case NAD1DPolicy::PURE: {
static int rand_uniform(kids_real *res_arr, int N=1, kids_real sigma=1.0)
static int rand_gaussian(kids_real *res_arr, int N=1, kids_real sigma=1.0, kids_real mu=0.0)
std::shared_ptr< Param > _param
Shared pointer to the Param object associated with this kernel.
Status & executeKernel(Status &stat)
Execute the kernel's algorithm and those of its children.
std::shared_ptr< DataSet > _dataset
Shared pointer to the DataSet object associated with this kernel.
virtual void setInputParam_impl(std::shared_ptr< Param > &PM)
Virtual function to set input parameters for the kernel implementation.
virtual const std::string getName()
Get the name of the kernel.
virtual Status & initializeKernel_impl(Status &stat)
Virtual function to initialize the kernel implementation.
virtual int getType() const
Get the type of the kernel.
virtual Status & executeKernel_impl(Status &stat)
Virtual function to execute the kernel implementation.
NAD1DPolicy::_type nad1d_type
virtual void setInputDataSet_impl(std::shared_ptr< DataSet > &DS)
Virtual function to set input data set for the kernel implementation.
static uval parse(const std::string &str)
CONSTEXPR_DECOR value_type conv(const uval &u)
#define LOC()
show the location information for debug
VARIABLE< kids_complex > p_sign
VARIABLE< kids_real > hess
VARIABLE< kids_real > p_sigma
VARIABLE< kids_real > vpes
VARIABLE< kids_real > x_sigma
VARIABLE< kids_real > grad
VARIABLE< kids_real > mass
std::size_t N
Number of nuclear degrees of freedom.
std::size_t F
Number of electronic degrees of freedom.
std::size_t FF
Product of F and F (F * F).
< http://warp.povusers.org/FunctionParser/fparser.html
int mspes_MORSE15C(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_ECR(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_CL1D(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_MORSE15(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_DAC(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_SAC3(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_MORSE3C(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
double kids_real
Alias for real number type.
int mspes_DRN(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_MORSE3B(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_DBG(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_MORSE3A(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_JC1D(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_DAG(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_SAC2(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_NA_I(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_SAC(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
int mspes_MORSE15E(double *V, double *dV, double *ddV, double *R, int flag, int rdim, int fdim)
constexpr std::complex< real_precision > iu(1.0L, 0.0L)
constexpr std::complex< real_precision > iz(0.0L, 0.0L)
static CONSTTYPE real_precision au_2_amu
1mea means we measure a quantity at 1*N level.
static CONSTTYPE real_precision au_2_ang
constexpr uint32_t hash(const char *str)
declaration of variables used in the program.