KIDS  ver-0.0.1
KIDS : Kernel Integrated Dynamics Simulator
Loading...
Searching...
No Matches
Model_LVCM.cpp
Go to the documentation of this file.
1#include "kids/Model_LVCM.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
10const std::string Model_LVCM::getName() { return "Model_LVCM"; }
11
13
14void Model_LVCM::setInputParam_impl(std::shared_ptr<Param> &PM) {
15 lvcm_type = LVCMPolicy::_from(_param->get_string("lvcm_flag", LOC(), "PYR3"));
16 classical_bath = _param->get_bool("classical_bath", LOC(), false);
17}
18
19void Model_LVCM::setInputDataSet_impl(std::shared_ptr<DataSet> &DS) {
20 Hsys = DS->def(DATA::model::Hsys);
21 w = DS->def(DATA::model::w);
22 memset(Hsys, 0, Dimension::FF * sizeof(kids_real));
23 switch (lvcm_type) {
24 case LVCMPolicy::PYR3:
25 case LVCMPolicy::PYR24: {
26 double H_unit = phys::au_2_ev;
27
28 // parameter for PYR3
29 double E_data_PYR3[2] = {3.94f, 4.84f};
30 double w_data_PYR3[3] = {0.074f, 0.126f, 0.118f};
31 double kcoeff_data_PYR3[4] = {-0.105f, 0.149f, 0.037f, -0.254f};
32 double lcoeff_data_PYR3[4] = {0.000f, 0.262f, 0.262f, 0.000f};
33
34 // parameter for PYR24
35 double E_data_PYR24[2] = {-0.4617f, 0.4617f};
36 double w_data_PYR24[24] = {0.0740f, 0.1273f, 0.1568f, 0.1347f, 0.3431f, 0.1157f, 0.3242f, 0.3621f,
37 0.2673f, 0.3052f, 0.0968f, 0.0589f, 0.0400f, 0.1726f, 0.2863f, 0.2484f,
38 0.1536f, 0.2105f, 0.0778f, 0.2294f, 0.1915f, 0.4000f, 0.3810f, 0.0936f};
39 double kcoeff_data_PYR24[46] = {-0.0964f, 0.1194f, 0.0470f, 0.2012f, 0.1594f, 0.0484f, 0.0308f, -0.0308f,
40 0.0782f, -0.0782f, 0.0261f, -0.0261f, 0.0717f, -0.0717f, 0.0780f, -0.0780f,
41 0.0560f, -0.0560f, 0.0625f, -0.0625f, 0.0188f, -0.0188f, 0.0112f, -0.0112f,
42 0.0069f, -0.0069f, 0.0265f, -0.0265f, 0.0433f, -0.0433f, 0.0361f, -0.0361f,
43 0.0210f, -0.0210f, 0.0281f, -0.0281f, 0.0102f, -0.0102f, 0.0284f, -0.0284f,
44 0.0196f, -0.0196f, 0.0306f, -0.0306f, 0.0269f, -0.0269f};
45 double lcoeff_data_PYR24[4] = {0.000f, 0.1825f, 0.1825f, 0.000f};
46
47
48 double *E_data;
49 double *w_data;
50 double *kcoeff_data;
51 double *lcoeff_data;
52 switch (lvcm_type) {
53 case LVCMPolicy::PYR3:
54 N_mode = 2;
55 N_coup = 1;
56 E_data = E_data_PYR3;
57 w_data = w_data_PYR3;
58 kcoeff_data = kcoeff_data_PYR3;
59 lcoeff_data = lcoeff_data_PYR3;
60 break;
61 case LVCMPolicy::PYR24:
62 N_mode = 23;
63 N_coup = 1;
64 E_data = E_data_PYR24;
65 w_data = w_data_PYR24;
66 kcoeff_data = kcoeff_data_PYR24;
67 lcoeff_data = lcoeff_data_PYR24;
68 break;
69 }
70
71 kcoeff = DS->def(DATA::model::kcoeff);
72 lcoeff = DS->def(DATA::model::lcoeff);
73
74 for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) { Hsys[ii] = E_data[i] / H_unit; }
75 for (int j = 0; j < Dimension::N; ++j) w[j] = w_data[j] / H_unit;
76 for (int j = 0, ji = 0; j < N_mode; ++j) {
77 for (int i = 0; i < Dimension::F; ++i, ++ji) { kcoeff[ji] = (kcoeff_data[ji] / H_unit) * sqrt(w[j]); }
78 }
79 for (int j = N_mode, j0ik = 0; j < Dimension::N; ++j) {
80 for (int ik = 0; ik < Dimension::FF; ++ik, ++j0ik) {
81 lcoeff[j0ik] = (lcoeff_data[j0ik] / H_unit) * sqrt(w[j]);
82 }
83 }
84 break;
85 }
86 case LVCMPolicy::BUTA5: {
87 double H_unit = phys::au_2_ev;
88
89 // parameter for BUTA5
90 double E_data_BUTA5[2] = {9.41165f, 9.95575f};
91 double w_data_BUTA5[5] = {0.1089f, 0.1773f, 0.2578f, 0.3713f, 0.0912f};
92 double kcoeff_data_BUTA5[8] = {-0.0531f, -0.0594f, 0.0115f, 0.0100f, //
93 -0.1628f, 0.3422f, -0.0403f, 0.0321f};
94 double lcoeff_data_BUTA5[4] = {0.000f, 0.2880f, 0.2880f, 0.000f};
95
96 double *E_data;
97 double *w_data;
98 double *kcoeff_data;
99 double *lcoeff_data;
100 N_mode = 4;
101 N_coup = 1;
102 E_data = E_data_BUTA5;
103 w_data = w_data_BUTA5;
104 kcoeff_data = kcoeff_data_BUTA5;
105 lcoeff_data = lcoeff_data_BUTA5;
106
107 kcoeff = DS->def(DATA::model::kcoeff);
108 lcoeff = DS->def(DATA::model::lcoeff);
109
110 for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) { Hsys[ii] = E_data[i] / H_unit; }
111 for (int j = 0; j < Dimension::N; ++j) w[j] = w_data[j] / H_unit;
112 for (int j = 0, ji = 0; j < N_mode; ++j) {
113 for (int i = 0; i < Dimension::F; ++i, ++ji) { kcoeff[ji] = (kcoeff_data[ji] / H_unit) * sqrt(w[j]); }
114 }
115 for (int j = N_mode, j0ik = 0; j < Dimension::N; ++j) {
116 for (int ik = 0; ik < Dimension::FF; ++ik, ++j0ik) {
117 lcoeff[j0ik] = (lcoeff_data[j0ik] / H_unit) * sqrt(w[j]);
118 }
119 }
120 break;
121 }
122 case LVCMPolicy::CRC2:
123 case LVCMPolicy::CRC5: {
124 double H_unit = phys::au_2_ev;
125
126 N_mode = 0;
128
129 const double kappa1 = -0.0169f;
130 const double kappa3 = -0.0272f;
131 const double lambda1a = 0.0328f;
132 const double lambda1b = -0.0978f;
133 const double lambda2a = 0.0095f;
134 const double lambda2b = 0.1014f;
135
136 double E_data_CRC5[3] = {0.0424f, 0.0424f, 0.4344f};
137 double w_data_CRC5[5] = {0.0129f, 0.0129f, 0.0342f, 0.0561f, 0.0561f};
138 double kcoeff_data_CRC5[1] = {0.0f};
139 double lcoeff_data_CRC5[45] = {
140 // data
141 0.0f, lambda1a, 0.0f, // x0
142 lambda1a, 0.0f, lambda1b, //
143 0.0f, lambda1b, 0.0f, //
144 -lambda1a, 0.0f, lambda1b, // x1
145 0.0f, lambda1a, 0.0f, //
146 lambda1b, 0.0f, 0.0f, //
147 kappa1, 0.0f, 0.0f, // x2
148 0.0f, kappa1, 0.0f, //
149 0.0f, 0.0f, kappa3, //
150 -lambda2a, 0.0f, lambda2b, // x3
151 0.0f, lambda2a, 0.0f, //
152 lambda2b, 0.0f, 0.0f, //
153 0.0f, lambda2a, 0.0f, // x4
154 lambda2a, 0.0f, lambda2b, //
155 0.0f, lambda2b, 0.0f //
156 };
157
158 double *E_data = E_data_CRC5;
159 double *w_data = w_data_CRC5;
160 double *kcoeff_data = kcoeff_data_CRC5;
161 double *lcoeff_data = lcoeff_data_CRC5;
162
163 kcoeff = DS->def(DATA::model::kcoeff);
164 lcoeff = DS->def(DATA::model::lcoeff);
165
166 for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) { Hsys[ii] = E_data[i] / H_unit; }
167 for (int j = 0; j < Dimension::N; ++j) w[j] = w_data[j] / H_unit;
168 for (int j = 0, ji = 0; j < N_mode; ++j) {
169 for (int i = 0; i < Dimension::F; ++i, ++ji) { kcoeff[ji] = (kcoeff_data[ji] / H_unit) * sqrt(w[j]); }
170 }
171 for (int j = N_mode, j0ik = 0; j < Dimension::N; ++j) {
172 for (int ik = 0; ik < Dimension::FF; ++ik, ++j0ik) {
173 lcoeff[j0ik] = (lcoeff_data[j0ik] / H_unit) * sqrt(w[j]);
174 }
175 }
176 break;
177 }
178 case LVCMPolicy::CED2:
179 case LVCMPolicy::CED3: {
180 const double lightspeed = 137.03599907444f;
181 const double epsilon0 = 0.25f / phys::math::pi;
182
183 double E_data_CED2[2] = {-0.6738f, -0.2798f};
184 double mu_data_CED2[4] = {0.000f, +1.034f, //
185 +1.034f, 0.000f};
186
187 double E_data_CED3[3] = {-0.6738f, -0.2798f, -0.1547f};
188 double mu_data_CED3[9] = {0.000f, +1.034f, 0.000f, //
189 +1.034f, 0.000f, -2.536f, //
190 0.000f, -2.536f, 0.000f};
191
192 double Lcav = 2.362e5;
193 double Rcav = Lcav / 2;
194
195 double *E_data;
196 double *mu_data;
197 switch (lvcm_type) {
198 case LVCMPolicy::CED2:
199 E_data = E_data_CED2;
200 mu_data = mu_data_CED2;
201 break;
202 case LVCMPolicy::CED3:
203 E_data = E_data_CED3;
204 mu_data = mu_data_CED3;
205 break;
206 }
207
208 N_mode = 0;
210
211 lcoeff = DS->def(DATA::model::lcoeff);
212
213 for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) { Hsys[ii] = E_data[i]; }
214 for (int j = 0, jik = 0; j < Dimension::N; ++j) {
215 w[j] = (2 * j + 1) * lightspeed * phys::math::pi / Lcav;
216 for (int ik = 0; ik < Dimension::FF; ++ik, ++jik) {
217 lcoeff[jik] = sqrt(2.0f / (epsilon0 * Lcav)) * sin((2 * j + 1) * phys::math::pi * Rcav / Lcav) *
218 w[j] * mu_data[ik];
219 }
220 }
221 break;
222 }
223 case LVCMPolicy::PYR2CED: {
224 N_mode = 1;
225 N_coup = 1;
226
227 double H_unit = phys::au_2_ev;
228
229 double gcoup = _param->get_double("gcoup", LOC(), phys::energy_d, 0.24f / H_unit);
230 double wcav = _param->get_double("wcav", LOC(), phys::energy_d, 0.62f / H_unit);
231
232 // parameter for PYR2
233 double E_data_PYR2[2] = {3.94f, 4.84f};
234 double w_data_PYR2[3] = {0.074f, 0.118f};
235 double kcoeff_data_PYR2[4] = {-0.105f, 0.149f};
236 double lcoeff_data_PYR2[4] = {0.000f, 0.262f, 0.262f, 0.000f};
237
238 double *E_data = E_data_PYR2;
239 double *w_data = w_data_PYR2;
240 double *kcoeff_data = kcoeff_data_PYR2;
241 double *lcoeff_data = lcoeff_data_PYR2;
242
243 kcoeff = DS->def(DATA::model::kcoeff);
244 lcoeff = DS->def(DATA::model::lcoeff);
245
246 for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1) {
247 Hsys[ii] = i / (Dimension::F / 2) * wcav + E_data[i % 2] / H_unit;
248 }
249 Hsys[0 * 4 + 3] = gcoup;
250 Hsys[1 * 4 + 2] = gcoup;
251 Hsys[2 * 4 + 1] = gcoup;
252 Hsys[3 * 4 + 0] = gcoup;
253 for (int j = 0; j < Dimension::N; ++j) w[j] = w_data[j] / H_unit;
254 for (int j = 0, ji = 0; j < N_mode; ++j) {
255 for (int i = 0; i < Dimension::F; ++i, ++ji) {
256 kcoeff[ji] = (kcoeff_data[j * 2 + i % 2] / H_unit) * sqrt(w[j]);
257 }
258 }
259 for (int j = N_mode, j0 = 0, j0ik = 0; j < Dimension::N; ++j, ++j0) {
260 for (int i = 0, ik = 0; i < Dimension::F; ++i) {
261 for (int k = 0; k < Dimension::F; ++k, ++ik, ++j0ik) {
262 if (i / (Dimension::F / 2) == k / (Dimension::F / 2))
263 lcoeff[j0ik] = (lcoeff_data[j0 * 4 + (i % 2) * 2 + k % 2] / H_unit) * sqrt(w[j]);
264 }
265 }
266 }
267 // ARRAY_SHOW(Hsys, Dimension::F, Dimension::F);
268 // ARRAY_SHOW(kcoeff, N_mode, Dimension::F);
269 // ARRAY_SHOW(lcoeff, N_coup, Dimension::FF);
270 break;
271 }
272 case LVCMPolicy::BEN5: {
273 break;
274 }
275 case LVCMPolicy::Read: {
276 std::string lvcm_readfile = _param->get_string("lvcm_readfile", LOC(), "lvcm.dat");
277 std::ifstream ifs(lvcm_readfile);
278 std::string H_unit_str;
279 std::string firstline;
280 getline(ifs, firstline);
281 std::stringstream sstr(firstline);
282 sstr >> H_unit_str;
283 double H_unit = phys::us::conv(phys::au::unit, phys::us::parse(H_unit_str));
284
285 // read E
286 int dsize;
287 std::string flag;
288 double val;
289 ifs >> flag >> dsize;
290 for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1)
291 if (ifs >> val) Hsys[ii] = val / H_unit;
292
293 // read kcoeff & lcoeff
294 for (int j = 0, ji = 0, j0ik = 0; j < Dimension::N; ++j) {
295 ifs >> flag >> dsize;
296 if (dsize == Dimension::F) {
297 if (ifs >> val) w[j] = val / H_unit; // freq
298 for (int i = 0; i < Dimension::F; ++i, ++ji)
299 if (ifs >> val) kcoeff[ji] = val / H_unit;
300 } else if (dsize = Dimension::FF) {
301 if (ifs >> val) w[j] = val / H_unit; // freq
302 for (int ik = 0; ik < Dimension::FF; ++ik, ++j0ik)
303 if (ifs >> val) lcoeff[j0ik] = val / H_unit;
304 } else {
305 exit(-1);
306 }
307 }
308 ifs >> flag >> dsize;
309 for (int i = 0, ii = 0; i < Dimension::F; ++i, ii += Dimension::Fadd1)
310 if (ifs >> val) Hsys[ii] = val / H_unit;
311
312 for (int i = 0; i < Dimension::FF; ++i)
313 if (ifs >> val) Hsys[i] = val / H_unit;
314 ifs.close();
315 }
316 }
317
319 x0 = DS->def(DATA::model::x0);
320 p0 = DS->def(DATA::model::p0);
321 x_sigma = DS->def(DATA::model::x_sigma);
322 p_sigma = DS->def(DATA::model::p_sigma);
323 switch (lvcm_type) {
324 case LVCMPolicy::CRC2:
325 case LVCMPolicy::CRC5: {
326 double reqb[5] = {0.0f, 14.3514f, -9.9699f, -7.0189f, 0.0f};
327 double alpw[5] = {0.4501f, 0.4286f, 0.6204f, 0.4535f, 0.5539f};
328 for (int j = 0; j < Dimension::N; ++j) {
329 x0[j] = reqb[j] / sqrt(w[j]);
330 p0[j] = 0.0f;
331 x_sigma[j] = alpw[j] / sqrt(w[j]);
332 p_sigma[j] = 0.5f * sqrt(w[j]) / alpw[j];
333 }
334 break;
335 }
336 default: {
337 for (int j = 0; j < Dimension::N; ++j) {
338 x0[j] = 0.0f;
339 p0[j] = 0.0f;
340 x_sigma[j] = sqrt(0.5f / w[j]);
341 p_sigma[j] = sqrt(0.5f * w[j]);
342 }
343 break;
344 }
345 }
346
347 // model field
348 mass = DS->def(DATA::model::mass);
349 for (int j = 0; j < Dimension::N; ++j) mass[j] = 1.0f;
350 vpes = DS->def(DATA::model::vpes);
351 grad = DS->def(DATA::model::grad);
352 hess = DS->def(DATA::model::hess);
353 V = DS->def(DATA::model::V);
354 dV = DS->def(DATA::model::dV);
355 // ddV = DS->def(DATA::model::ddV);
356
357 // init & integrator
358 x = DS->def(DATA::integrator::x);
359 p = DS->def(DATA::integrator::p);
360}
361
363 for (int iP = 0; iP < Dimension::P; ++iP) {
364 kids_real *x = this->x + iP * Dimension::N;
365 kids_real *p = this->p + iP * Dimension::N;
366
369 for (int j = 0; j < Dimension::N; ++j) {
370 if (classical_bath) {
371 x[j] = x0[j];
372 p[j] = p0[j];
373 } else {
374 x[j] = x0[j] + x[j] * x_sigma[j];
375 p[j] = p0[j] + p[j] * p_sigma[j];
376 }
377 }
378 }
379
380 _dataset->def_real("init.x", x, Dimension::PN);
381 _dataset->def_real("init.p", p, Dimension::PN);
382 executeKernel(stat);
383 return stat;
384}
385
387 for (int iP = 0; iP < Dimension::P; ++iP) {
388 kids_real *x = this->x + iP * Dimension::N;
389 kids_real *vpes = this->vpes + iP;
390 kids_real *grad = this->grad + iP * Dimension::N;
391 kids_real *hess = this->hess + iP * Dimension::NN;
392 kids_real *V = this->V + iP * Dimension::FF;
393 kids_real *dV = this->dV + iP * Dimension::NFF;
394 // ARRAY_SHOW(x, 1, Dimension::N);
395
396 // calculate nuclear vpes and grad
397 double term = 0.0f;
398 for (int j = 0; j < Dimension::N; ++j) {
399 term += w[j] * w[j] * x[j] * x[j];
400 grad[j] = w[j] * w[j] * x[j];
401 }
402 vpes[0] = 0.5 * term;
403
404 // electronic pes
405 memset(V, 0, Dimension::FF * sizeof(kids_real));
406 for (int ik = 0; ik < Dimension::FF; ++ik) V[ik] = Hsys[ik];
407 // ARRAY_SHOW(V, Dimension::F, Dimension::F);
408
409 for (int j = 0, ji = 0; j < N_mode; ++j) {
410 for (int i = 0, ii = 0; i < Dimension::F; ++i, ++ji, ii += Dimension::Fadd1) {
411 V[ii] += kcoeff[ji] * x[j]; //
412 };
413 }
414 for (int j = N_mode, j0ik = 0; j < Dimension::N; ++j) {
415 for (int ik = 0; ik < Dimension::FF; ++ik, ++j0ik) {
416 V[ik] += lcoeff[j0ik] * x[j]; // slow
417 };
418 }
419
420 // ARRAY_SHOW(V, Dimension::F, Dimension::F);
421
422 if (count_exec == 0) {
423 // N_mode
424 for (int j = 0, ji = 0, jFF = 0; j < N_mode; ++j, jFF += Dimension::FF) {
425 for (int i = 0, jii = jFF; i < Dimension::F; ++i, ++ji, jii += Dimension::Fadd1) {
426 dV[jii] = kcoeff[ji]; //
427 }
428 }
429 // N_coup
430 for (int j = N_mode, jFF = N_mode * Dimension::FF, j0FF = 0; j < Dimension::N;
431 ++j, jFF += Dimension::FF, j0FF += Dimension::FF) {
432 for (int ik = 0, jik = jFF, j0ik = j0FF; ik < Dimension::FF; ++ik, ++jik, ++j0ik) {
433 dV[jik] = lcoeff[j0ik];
434 }
435 }
436 }
437 }
438 return stat;
439}
440
441
442}; // namespace PROJECT_NS
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
Status & executeKernel_impl(Status &stat)
Virtual function to execute the kernel implementation.
void setInputParam_impl(std::shared_ptr< Param > &PM)
Virtual function to set input parameters for the kernel implementation.
Status & initializeKernel_impl(Status &stat)
Virtual function to initialize the kernel implementation.
virtual int getType() const
Get the type of the kernel.
LVCMPolicy::_type lvcm_type
Definition Model_LVCM.h:34
void setInputDataSet_impl(std::shared_ptr< DataSet > &DS)
Virtual function to set input data set for the kernel implementation.
virtual const std::string getName()
Get the name of the kernel.
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_real > x
VARIABLE< kids_real > hess
VARIABLE< kids_real > p_sigma
VARIABLE< kids_real > vpes
VARIABLE< kids_real > lcoeff
VARIABLE< kids_real > V
VARIABLE< kids_real > p0
VARIABLE< kids_real > dV
VARIABLE< kids_real > x_sigma
VARIABLE< kids_real > Hsys
VARIABLE< kids_real > x0
VARIABLE< kids_real > w
VARIABLE< kids_real > kcoeff
VARIABLE< kids_real > grad
VARIABLE< kids_real > mass
std::size_t PN
Product of P and N (P * N).
Definition vars_list.cpp:14
std::size_t NFF
Product of N, F, and F (N * F * F).
Definition vars_list.cpp:22
std::size_t N
Number of nuclear degrees of freedom.
Definition vars_list.cpp:10
std::size_t P
Number of parallel trajectories (swarms of trajectories) in each Monte Carlo run.
Definition vars_list.cpp:9
std::size_t NN
Product of N and N (N * N).
Definition vars_list.cpp:20
std::size_t F
Number of electronic degrees of freedom.
Definition vars_list.cpp:11
std::size_t Fadd1
F plus 1 (F + 1).
Definition vars_list.cpp:24
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
double kids_real
Alias for real number type.
Definition Types.h:59
constexpr real_precision pi
pi
Definition phys.h:30
constexpr dimension7 energy_d
[M*L^2/T^2] force integrate length
Definition phys.h:251
static CONSTTYPE real_precision au_2_ev
Definition phys.h:994
constexpr uint32_t hash(const char *str)
Definition hash_fnv1a.h:12
declaration of variables used in the program.