KIDS  ver-0.0.1
KIDS : Kernel Integrated Dynamics Simulator
Loading...
Searching...
No Matches
Model_NAD1D.cpp
Go to the documentation of this file.
1#include "kids/Model_NAD1D.h"
2
4#include "kids/hash_fnv1a.h"
5#include "kids/macro_utils.h"
6#include "kids/vars_list.h"
7
8namespace PROJECT_NS {
9
10double mspes_parm[100];
11
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;
14
15 V[0] = (R[0] > 0) ? V0 * (1.0 - exp(-a * R[0])) : -V0 * (1.0 - exp(a * R[0]));
16 V[3] = -V[0];
17 V[1] = V1 * exp(-b * R[0] * R[0]);
18 V[2] = V[1];
19 if (flag < 1) return 0;
20
21 dV[0] = (R[0] > 0) ? V0 * a * exp(-a * R[0]) : V0 * a * exp(a * R[0]);
22 dV[3] = -dV[0];
23 dV[1] = -2 * b * R[0] * V[1];
24 dV[2] = dV[1];
25 if (flag < 2) return 0;
26
27 // ddV
28 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
29 return 0;
30}
31
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;
34
35 V[0] = V0 * tanh(a * R[0]);
36 V[3] = -V[0];
37 V[1] = V1 * exp(-b * R[0] * R[0]);
38 V[2] = V[1];
39 if (flag < 1) return 0;
40
41 double tmp = cosh(a * R[0]);
42 dV[0] = V0 * a / (tmp * tmp);
43 dV[3] = -dV[0];
44 dV[1] = -2 * b * R[0] * V[1];
45 dV[2] = dV[1];
46 if (flag < 2) return 0;
47
48 // ddV
49 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
50 return 0;
51}
52
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;
55
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));
59 V[2] = V[1];
60 if (flag < 1) return 0;
61
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];
66 dV[2] = dV[1];
67 if (flag < 2) return 0;
68
69 // ddV
70 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
71 return 0;
72}
73
74
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;
77
78 V[0] = 0.0e0;
79 V[3] = -V0 * exp(-a * R[0] * R[0]) + E0;
80 V[1] = V1 * exp(-b * R[0] * R[0]);
81 V[2] = V[1];
82 if (flag < 1) return 0;
83
84 dV[0] = 0.0e0;
85 dV[3] = 2 * a * R[0] * V0 * exp(-a * R[0] * R[0]);
86 dV[1] = -2 * b * R[0] * V[1];
87 dV[2] = dV[1];
88 if (flag < 2) return 0;
89
90 // ddV
91 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
92 return 0;
93}
94
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;
97
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]));
100 V[2] = V[1];
101 if (flag < 1) return 0;
102
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]);
105 dV[2] = dV[1];
106 if (flag < 2) return 0;
107
108 // ddV
109 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
110 return 0;
111}
112
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;
115
116 V[0] = -V0, V[3] = V0;
117 if (R[0] < -Z) {
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));
121 } else {
122 V[1] = V1 * exp(-b * (R[0] + Z)) + V1 * (2.0 - exp(-b * (R[0] - Z)));
123 }
124 V[2] = V[1];
125 if (flag < 1) return 0;
126
127 dV[0] = 0.0e0, dV[3] = 0.0e0;
128 if (R[0] < -Z) {
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);
132 } else {
133 dV[1] = V1 * exp(-b * (R[0] + Z)) * (-b) + V1 * (-exp(-b * (R[0] - Z))) * (-b);
134 }
135 dV[2] = dV[1];
136 if (flag < 2) return 0;
137
138 // ddV
139 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
140 return 0;
141}
142
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;
145
146 V[0] = -V0, V[3] = V0;
147 if (R[0] < -Z) {
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;
151 } else {
152 V[1] = V1 * exp(-b * (R[0] - Z)) - V1 * exp(-b * (R[0] + Z));
153 }
154 V[2] = V[1];
155 if (flag < 1) return 0;
156
157 dV[0] = 0.0e0, dV[3] = 0.0e0;
158 if (R[0] < -Z) {
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);
162 } else {
163 dV[1] = V1 * exp(-b * (R[0] - Z)) * (-b) - V1 * exp(-b * (R[0] + Z)) * (-b);
164 }
165 dV[2] = dV[1];
166 if (flag < 2) return 0;
167
168 // ddV
169 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
170 return 0;
171}
172
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;
175
176 V[0] = 0.0e0;
177 V[3] = E0;
178 V[1] = V1 * (exp(-b * (R[0] - Z) * (R[0] - Z)) + exp(-b * (R[0] + Z) * (R[0] + Z)));
179 V[2] = V[1];
180 if (flag < 1) return 0;
181
182 dV[0] = 0.0e0;
183 dV[3] = 0.0e0;
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)));
186 dV[2] = dV[1];
187 if (flag < 2) return 0;
188
189 // ddV
190 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
191 return 0;
192}
193
194
195int mspes_MORSE3A(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
196 // V
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};
204
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];
208
209 V[1] = Aijx[2] * exp(-Aeijx[2] * (R[0] - Rijx[2]) * (R[0] - Rijx[2])); // V(0,1)
210 V[2] = Aijx[1] * exp(-Aeijx[1] * (R[0] - Rijx[1]) * (R[0] - Rijx[1])); // V(0,2)
211 V[5] = Aijx[0] * exp(-Aeijx[0] * (R[0] - Rijx[0]) * (R[0] - Rijx[0])); // V(1,2)
212 V[3] = V[1]; // V(1,0)
213 V[6] = V[2]; // V(2,0)
214 V[7] = V[5]; // V(2,1)
215
216 if (flag < 1) return 0;
217
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];
221
222 dV[1] = V[1] * (-2 * Aeijx[2] * (R[0] - Rijx[2])); // dV(0,1)
223 dV[2] = V[2] * (-2 * Aeijx[1] * (R[0] - Rijx[1])); // dV(0,2)
224 dV[5] = V[5] * (-2 * Aeijx[0] * (R[0] - Rijx[0])); // dV(1,2)
225 dV[3] = dV[1]; // dV(1,0)
226 dV[6] = dV[2]; // dV(2,0)
227 dV[7] = dV[5]; // dV(2,1)
228
229 if (flag < 2) return 0;
230 // for (int i = 0; i < 9; ++i) ddV[i] = 0;
231 return 0;
232}
233
234int mspes_MORSE3B(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
235 // V
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};
243
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];
247
248 V[1] = Aijx[2] * exp(-Aeijx[2] * (R[0] - Rijx[2]) * (R[0] - Rijx[2])); // V(0,1)
249 V[2] = Aijx[1] * exp(-Aeijx[1] * (R[0] - Rijx[1]) * (R[0] - Rijx[1])); // V(0,2)
250 V[5] = Aijx[0] * exp(-Aeijx[0] * (R[0] - Rijx[0]) * (R[0] - Rijx[0])); // V(1,2)
251 V[3] = V[1]; // V(1,0)
252 V[6] = V[2]; // V(2,0)
253 V[7] = V[5]; // V(2,1)
254
255 if (flag < 1) return 0;
256
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];
260
261 dV[1] = V[1] * (-2 * Aeijx[2] * (R[0] - Rijx[2])); // dV(0,1)
262 dV[2] = V[2] * (-2 * Aeijx[1] * (R[0] - Rijx[1])); // dV(0,2)
263 dV[5] = V[5] * (-2 * Aeijx[0] * (R[0] - Rijx[0])); // dV(1,2)
264 dV[3] = dV[1]; // dV(1,0)
265 dV[6] = dV[2]; // dV(2,0)
266 dV[7] = dV[5]; // dV(2,1)
267
268 if (flag < 2) return 0;
269 // for (int i = 0; i < 9; ++i) ddV[i] = 0;
270 return 0;
271}
272
273int mspes_MORSE3C(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
274 // V
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};
282
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];
286
287 V[1] = Aijx[2] * exp(-Aeijx[2] * (R[0] - Rijx[2]) * (R[0] - Rijx[2])); // V(0,1)
288 V[2] = Aijx[1] * exp(-Aeijx[1] * (R[0] - Rijx[1]) * (R[0] - Rijx[1])); // V(0,2)
289 V[5] = Aijx[0] * exp(-Aeijx[0] * (R[0] - Rijx[0]) * (R[0] - Rijx[0])); // V(1,2)
290 V[3] = V[1]; // V(1,0)
291 V[6] = V[2]; // V(2,0)
292 V[7] = V[5]; // V(2,1)
293
294 // if (R[0] < 0) P[0] = fabs(P[0]);
295
296 if (flag < 1) return 0;
297
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];
301
302 dV[1] = V[1] * (-2 * Aeijx[2] * (R[0] - Rijx[2])); // dV(0,1)
303 dV[2] = V[2] * (-2 * Aeijx[1] * (R[0] - Rijx[1])); // dV(0,2)
304 dV[5] = V[5] * (-2 * Aeijx[0] * (R[0] - Rijx[0])); // dV(1,2)
305 dV[3] = dV[1]; // dV(1,0)
306 dV[6] = dV[2]; // dV(2,0)
307 dV[7] = dV[5]; // dV(2,1)
308
309 if (flag < 2) return 0;
310 // for (int i = 0; i < 9; ++i) ddV[i] = 0;
311 return 0;
312}
313
314
315int mspes_MORSE15(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
316 // V
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]));
327 V[i * 15] = V[i];
328 }
329 if (flag < 1) return 0;
330
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]));
336 dV[i * 15] = dV[i];
337 }
338
339 if (flag < 2) return 0;
340 // for (int i = 0; i < 225; ++i) ddV[i] = 0;
341 return 0;
342}
343
344int mspes_MORSE15C(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
345 // V
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;
355 V[i] = 0.25e0 * Dc;
356 V[i * 15] = V[i];
357 }
358 if (flag < 1) return 0;
359
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]);
364 dV[i] = 0.0e0;
365 dV[i * 15] = dV[i];
366 }
367
368 if (flag < 2) return 0;
369 // for (int i = 0; i < 225; ++i) ddV[i] = 0;
370 return 0;
371}
372
373int mspes_MORSE15E(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
374 // V
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]);
385 V[i * 15] = V[i];
386 }
387 if (flag < 1) return 0;
388
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]);
393 dV[i] = -0.2 * V[i];
394 dV[i * 15] = dV[i];
395 }
396
397 if (flag < 2) return 0;
398 // for (int i = 0; i < 225; ++i) ddV[i] = 0;
399 return 0;
400}
401
402int mspes_CL1D(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
403 // @deps. on mass and mod_W
404 double meanV = 0.5e0 * mspes_parm[3] * mspes_parm[3] * R[0] * R[0]; // mass = 1
405 V[0] = meanV + mspes_parm[0] + mspes_parm[2] * R[0];
406 V[3] = meanV - mspes_parm[0] - mspes_parm[2] * R[0];
407 V[1] = mspes_parm[1];
408 V[2] = mspes_parm[1];
409 if (flag < 1) return 0;
410
411 dV[0] = mspes_parm[3] * mspes_parm[3] * R[0] + mspes_parm[2];
412 dV[3] = -dV[0];
413 dV[1] = 0;
414 dV[2] = 0;
415 if (flag < 2) return 0;
416
417 // ddV
418 for (int i = 0; i < 4; ++i) ddV[i] = 0;
419 return 0;
420}
421
422int mspes_JC1D(double* V, double* dV, double* ddV, double* R, int flag, int rdim, int fdim) {
423 // @deps. on mass and mod_W
424 double meanV = 0.5e0 * mspes_parm[3] * mspes_parm[3] * R[0] * R[0]; // mass = 1
425 V[0] = meanV + mspes_parm[0];
426 V[3] = meanV - mspes_parm[0];
427 V[1] = mspes_parm[1] + mspes_parm[2] * R[0];
428 V[2] = mspes_parm[1] + mspes_parm[2] * R[0];
429 if (flag < 1) return 0;
430
431 dV[0] = mspes_parm[3] * mspes_parm[3] * R[0];
432 dV[3] = -dV[0];
433 dV[1] = mspes_parm[2];
434 dV[2] = mspes_parm[2];
435 if (flag < 2) return 0;
436
437 // ddV
438 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
439 return 0;
440}
441
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;
461
462 double r = R[0];
463 double r2 = r * r;
464 double r4 = r2 * r2;
465 double r5 = r4 * r;
466 double r6 = r2 * r4;
467 double r7 = r6 * r;
468 double r8 = r4 * r4;
469 double r9 = r8 * r;
470 double r12 = r6 * r6;
471 double r13 = r12 * r;
472
473 double Vcent = mspes_parm[1] / (2 * mred * r2);
474
475 V[0] = Vcent //
476 + (Acov + Bcov12 / r12) * exp(-r / Rcov) //
477 - Ccov / r6; //
478 V[3] = Vcent //
479 + (Aion + Bion8 / r8) * exp(-r / Rion) //
480 - 1 / r //
481 - 0.5 * (al_Mp + al_Xm) / r4 //
482 - Cion / r6 //
483 - 2 * al_Mp * al_Xm / r7 //
484 + Eth;
485 V[1] = Aoff * exp(-r / Roff);
486 V[2] = V[1];
487
488 if (flag < 1) return 0;
489
490 double dVcent = -2 * mspes_parm[1] / (2 * mred * r2 * r);
491
492 dV[0] = dVcent //
493 + (-12 * Bcov12 / r13) * exp(-r / Rcov) //
494 + (Acov + Bcov12 / r12) * (-1 / Rcov) * exp(-r / Rcov) //
495 + 6 * Ccov / r7; //
496 dV[3] = dVcent //
497 + (-8 * Bion8 / r9) * exp(-r / Rion) //
498 + (Aion + Bion8 / r8) * (-1 / Rion) * exp(-r / Rion) //
499 + 1 / r2 //
500 + 2 * (al_Mp + al_Xm) / r5 //
501 + 6 * Cion / r7 //
502 + 14 * al_Mp * al_Xm / r8; //
503 dV[1] = (-1 / Roff) * Aoff * exp(-r / Roff);
504 dV[2] = dV[1];
505
506 if (flag < 2) return 0;
507
508 // ddV
509 // for (int i = 0; i < 4; ++i) ddV[i] = 0;
510 return 0;
511}
512
513const std::string Model_NAD1D::getName() { return "Model_NAD1D"; }
514
516
517void Model_NAD1D::setInputParam_impl(std::shared_ptr<Param>& PM) {
518 nad1d_type = NAD1DPolicy::_from(_param->get_string("nad1d_flag", LOC(), "SAC"));
519
520 // revise tend and dt
521 switch (nad1d_type) {
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);
535
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;
540 break;
541 }
542 case NAD1DPolicy::NA_I: {
543 // CHECK_EQ(F, 2);
544 // initializetion in au
545 double mass_Na = 22.989769282e0 / phys::au_2_amu;
546 double mass_I = 126.904473e0 / phys::au_2_amu;
547 double mred = (mass_Na * mass_I) / (mass_Na + mass_I);
548
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);
552
553 double x0_max = 160.0e0 / phys::au_2_ang;
554 double vel = sqrt(2 * parm_E / mred);
555 double tend_rev = 2 * x0_max / vel;
556 double dt_rev = tend_rev / 50000;
557
558 if (tend_read < 0) (*(_param->pjson()))["tend"] = tend_rev;
559 if (dt_read < 0) (*(_param->pjson()))["dt"] = dt_rev;
560 break;
561 }
562 }
563};
564
565void Model_NAD1D::setInputDataSet_impl(std::shared_ptr<DataSet>& DS) {
566 Hsys = DS->def_real("model.Hsys", Dimension::FF);
567 memset(Hsys, 0, Dimension::FF * sizeof(kids_real));
568
569 if (nad1d_type == NAD1DPolicy::PURE) {
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);
575 sstr >> H_unit_str;
576 double H_unit = phys::us::conv(phys::au::unit, phys::us::parse(H_unit_str));
577
578 kids_real val;
579 for (int i = 0; i < Dimension::FF; ++i)
580 if (ifs >> val) Hsys[i] = val / H_unit;
581 ifs.close();
582 }
583
584 // model field
585 mass = DS->def(DATA::model::mass);
586 vpes = DS->def(DATA::model::vpes);
587 grad = DS->def(DATA::model::grad);
588 hess = DS->def(DATA::model::hess);
589 V = DS->def(DATA::model::V);
590 dV = DS->def(DATA::model::dV);
591 // ddV = DS->def(DATA::model::ddV);
592
593 x0 = DS->def(DATA::model::x0);
594 p0 = DS->def(DATA::model::p0);
595 x_sigma = DS->def(DATA::model::x_sigma);
596 p_sigma = DS->def(DATA::model::p_sigma);
597
598 // init & integrator
599 x = DS->def(DATA::integrator::x);
600 p = DS->def(DATA::integrator::p);
602
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;
608
609 DS->def_real("init.x", Dimension::N);
610 DS->def_real("init.p", Dimension::N);
611
612 mass[0] = _param->get_double("m0", LOC(), 2000.0e0);
613 x0[0] = _param->get_double("x0", LOC(), 100.0e0);
614 p0[0] = _param->get_double("p0", LOC(), 100.0e0);
615 double varx = _param->get_double("varx", LOC(), 0.5e0);
616 double varp = _param->get_double("varp", LOC(), 0.5e0);
617 x_sigma[0] = sqrt(varx);
618 p_sigma[0] = sqrt(varp);
619
620 switch (nad1d_type) {
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: {
628 mass[0] = 2000.0e0;
629 break;
630 }
631 case NAD1DPolicy::SAC3: { // asymmetrical SAC
632 mass[0] = 1980.0e0;
633 double gammawidth = 0.25e0;
634 x_sigma[0] = sqrt(0.5e0 / gammawidth);
635 p_sigma[0] = sqrt(0.5e0 * gammawidth);
636 break;
637 }
638 case NAD1DPolicy::MORSE3A:
639 case NAD1DPolicy::MORSE3B:
640 case NAD1DPolicy::MORSE3C: {
641 // CHECK_EQ(F, 3);
642 mass[0] = 20000.0e0;
643 switch (nad1d_type) {
644 case NAD1DPolicy::MORSE3A:
645 x0[0] = 2.9e0;
646 break;
647 case NAD1DPolicy::MORSE3B:
648 x0[0] = 3.3e0;
649 break;
650 case NAD1DPolicy::MORSE3C:
651 x0[0] = 2.1e0;
652 break;
653 }
654 p0[0] = 0.0e0;
655 double wground = 5.0e-03;
656 x_sigma[0] = sqrt(0.5e0 / (mass[0] * wground));
657 p_sigma[0] = 0.5e0 / x_sigma[0];
658 break;
659 }
660 case NAD1DPolicy::MORSE15:
661 case NAD1DPolicy::MORSE15C:
662 case NAD1DPolicy::MORSE15E: {
663 // CHECK_EQ(F, 15);
664 mass[0] = 2000.0e0;
665 x0[0] = 13.0e0;
666 p0[0] = -30.0e0;
667 double Dg = 0.2e0, alpha = 0.4e0; // 2*alpha^2 * D = m*w^2
668 x_sigma[0] = sqrt(0.5e0 / std::sqrt(mass[0] * 2 * alpha * alpha * Dg));
669 p_sigma[0] = 0.5e0 / x_sigma[0];
670 break;
671 }
672 case NAD1DPolicy::CL1D:
673 case NAD1DPolicy::JC1D: {
674 // CHECK_EQ(F, 2);
675 mass[0] = 1.0e0;
676 mspes_parm[0] = _param->get_double("nad1d_e", LOC(), 1.0e0);
677 mspes_parm[1] = _param->get_double("nad1d_d", LOC(), 1.0e0);
678 mspes_parm[2] = _param->get_double("nad1d_c", LOC(), 1.0e0);
679 mspes_parm[3] = _param->get_double("nad1d_w", LOC(), 1.0e0);
680 break;
681 }
682 case NAD1DPolicy::NA_I: {
683 // CHECK_EQ(F, 2);
684 // initializetion in au
685 double mass_Na = 22.989769282e0 / phys::au_2_amu;
686 double mass_I = 126.904473e0 / phys::au_2_amu;
687 mass[0] = (mass_Na * mass_I) / (mass_Na + mass_I);
688
689 mspes_parm[0] = _param->get_double("parm_E", LOC(), 1.0e0);
690 x0[0] = 160.0e0 / phys::au_2_ang;
691 p0[0] = -sqrt(2 * mass[0] * mspes_parm[0]);
692 break;
693 }
694 }
695};
696
698 switch (nad1d_type) {
699 case NAD1DPolicy::NA_I: {
700 double bmax = 9.0e0 / phys::au_2_ang;
701 double randu;
703 double b = sqrt(randu) * bmax;
704 mspes_parm[1] = 2 * mass[0] * mspes_parm[0] * b * b; // l^2 = 2 m E b^2
705 for (int j = 0; j < Dimension::N; ++j) {
706 x[j] = x0[j];
707 p[j] = p0[j];
708 }
709 break;
710 }
711 default: {
714 for (int j = 0; j < Dimension::N; ++j) {
715 x[j] = x0[j] + 0 * x[j] * x_sigma[j];
716 p[j] = p0[j] + 0 * p[j] * p_sigma[j];
717 }
718 break;
719 }
720 }
721 if (p[0] >= 0) {
723 } else {
725 }
726 _dataset->def_real("init.x", x, Dimension::N);
727 _dataset->def_real("init.p", p, Dimension::N);
728
729 executeKernel(stat);
730 return stat;
731}
732
734 switch (nad1d_type) {
735 case NAD1DPolicy::SAC:
736 mspes_SAC(V, dV, ddV, x, 1, 1, Dimension::F);
737 break;
738 case NAD1DPolicy::SAC2:
739 mspes_SAC2(V, dV, ddV, x, 1, 1, Dimension::F);
740 break;
741 case NAD1DPolicy::SAC3:
742 mspes_SAC3(V, dV, ddV, x, 1, 1, Dimension::F);
743 break;
744 case NAD1DPolicy::DAC:
745 mspes_DAC(V, dV, ddV, x, 1, 1, Dimension::F);
746 break;
747 case NAD1DPolicy::ECR:
748 mspes_ECR(V, dV, ddV, x, 1, 1, Dimension::F);
749 break;
750 case NAD1DPolicy::DBG:
751 mspes_DBG(V, dV, ddV, x, 1, 1, Dimension::F);
752 break;
753 case NAD1DPolicy::DAG:
754 mspes_DAG(V, dV, ddV, x, 1, 1, Dimension::F);
755 break;
756 case NAD1DPolicy::DRN:
757 mspes_DRN(V, dV, ddV, x, 1, 1, Dimension::F);
758 break;
759 case NAD1DPolicy::MORSE3A:
760 mspes_MORSE3A(V, dV, ddV, x, 1, 1, Dimension::F);
761 break;
762 case NAD1DPolicy::MORSE3B:
763 mspes_MORSE3B(V, dV, ddV, x, 1, 1, Dimension::F);
764 break;
765 case NAD1DPolicy::MORSE3C:
766 mspes_MORSE3C(V, dV, ddV, x, 1, 1, Dimension::F);
767 break;
768 case NAD1DPolicy::MORSE15:
769 mspes_MORSE15(V, dV, ddV, x, 1, 1, Dimension::F);
770 break;
771 case NAD1DPolicy::MORSE15C:
772 mspes_MORSE15C(V, dV, ddV, x, 1, 1, Dimension::F);
773 break;
774 case NAD1DPolicy::MORSE15E:
775 mspes_MORSE15E(V, dV, ddV, x, 1, 1, Dimension::F);
776 break;
777 case NAD1DPolicy::JC1D:
778 mspes_JC1D(V, dV, ddV, x, 1, 1, Dimension::F);
779 break;
780 case NAD1DPolicy::CL1D:
781 mspes_CL1D(V, dV, ddV, x, 1, 1, Dimension::F);
782 break;
783 case NAD1DPolicy::NA_I:
784 mspes_NA_I(V, dV, ddV, x, 1, 1, Dimension::F);
785 break;
786 case NAD1DPolicy::PURE: {
787 if (count_exec == 0) {
788 for (int i = 0; i < Dimension::FF; ++i) { V[i] = Hsys[i], dV[i] = 0.0e0; }
789 }
790 break;
791 }
792 }
793 if (p[0] >= 0) {
795 } else {
797 }
798 return stat;
799}
800
801
802}; // namespace PROJECT_NS
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.
Definition Kernel.h:273
Status & executeKernel(Status &stat)
Execute the kernel's algorithm and those of its children.
Definition Kernel.cpp:48
std::shared_ptr< DataSet > _dataset
Shared pointer to the DataSet object associated with this kernel.
Definition Kernel.h:278
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.
kids_complex * p_sign
Definition Model_NAD1D.h:51
NAD1DPolicy::_type nad1d_type
Definition Model_NAD1D.h:40
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)
Definition phys.h:790
CONSTEXPR_DECOR value_type conv(const uval &u)
Definition phys.h:772
#define LOC()
show the location information for debug
Definition fmt.h:49
int count_exec
Counter for the number of executions performed by this kernel.
Definition Kernel.h:258
#define FUNCTION_NAME
Definition macro_utils.h:9
VARIABLE< kids_real > p
VARIABLE< kids_complex > p_sign
VARIABLE< kids_real > x
VARIABLE< kids_real > hess
VARIABLE< kids_real > p_sigma
VARIABLE< kids_real > vpes
VARIABLE< kids_real > V
VARIABLE< kids_real > p0
VARIABLE< kids_real > dV
VARIABLE< kids_real > x_sigma
VARIABLE< kids_real > x0
VARIABLE< kids_real > grad
VARIABLE< kids_real > mass
std::size_t N
Number of nuclear degrees of freedom.
Definition vars_list.cpp:10
std::size_t F
Number of electronic degrees of freedom.
Definition vars_list.cpp:11
std::size_t FF
Product of F and F (F * F).
Definition vars_list.cpp:21
< http://warp.povusers.org/FunctionParser/fparser.html
Definition Context.h:39
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)
double mspes_parm[100]
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.
Definition Types.h:59
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.
Definition phys.h:992
static CONSTTYPE real_precision au_2_ang
Definition phys.h:993
constexpr uint32_t hash(const char *str)
Definition hash_fnv1a.h:12
declaration of variables used in the program.