Skip to content

Commit a527d47

Browse files
alexvakimovalexvakimov
committed
Merge branch 'devel' of https://github.com/Quantum-Dynamics-Hub/libra-code into devel
2 parents b36b2e3 + 41080bf commit a527d47

19 files changed

+773
-151
lines changed

‎src/dyn/Dynamics.cpp‎

Lines changed: 170 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -1120,6 +1120,76 @@ void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* Ham, nHamiltonia
11201120

11211121
}
11221122

1123+
// Learn * vs &
1124+
void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& Ham, dyn_control_params& prms, Random& rnd){
1125+
1126+
int num_el = prms.num_electronic_substeps;
1127+
double dt = prms.dt / num_el;
1128+
double beta = (hartree/boltzmann) / prms.Temperature;
1129+
double a = prms.kcrpmd_a;
1130+
double b = prms.kcrpmd_b;
1131+
double c = prms.kcrpmd_c;
1132+
double d = prms.kcrpmd_d;
1133+
double eta = prms.kcrpmd_eta;
1134+
double gam = prms.kcrpmd_gamma;
1135+
double gamKP = prms.kcrpmd_gammaKP;
1136+
double sigma;
1137+
double xi;
1138+
double theta;
1139+
1140+
vector<double>& m_aux_var = dyn_var.m_aux_var;
1141+
vector<double>& y_aux_var = dyn_var.y_aux_var;
1142+
vector<double>& p_aux_var = dyn_var.p_aux_var;
1143+
vector<double>& f_aux_var = dyn_var.f_aux_var;
1144+
1145+
if(y_aux_var[0] >= -0.5 && y_aux_var[0] <= 0.5){
1146+
if(gamKP == 0.0){
1147+
p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
1148+
y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt;
1149+
f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
1150+
p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
1151+
}
1152+
else{
1153+
sigma = sqrt(2.0 * gamKP / (beta * m_aux_var[0]));
1154+
xi = rnd.normal();
1155+
theta = rnd.normal();
1156+
p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
1157+
+ 0.5 * (f_aux_var[0] - gamKP * p_aux_var[0]) * dt
1158+
- 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gamKP * sigma * m_aux_var[0] * sqrt(pow(dt,3))
1159+
- 0.125 * gamKP * (f_aux_var[0] - gamKP * p_aux_var[0]) * pow(dt,2);
1160+
y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt + 0.5 * theta / sqrt(3.0) * sigma * sqrt(pow(dt,3));
1161+
f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
1162+
p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
1163+
+ 0.5 * (f_aux_var[0] - gamKP * p_aux_var[0]) * dt
1164+
- 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gamKP * sigma * m_aux_var[0] * sqrt(pow(dt,3))
1165+
- 0.125 * gamKP * (f_aux_var[0] - gamKP * p_aux_var[0]) * pow(dt,2);
1166+
}
1167+
}
1168+
else{
1169+
if(gam == 0.0){
1170+
p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
1171+
y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt;
1172+
f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
1173+
p_aux_var[0] += 0.5 * f_aux_var[0] * dt;
1174+
}
1175+
else{
1176+
sigma = sqrt(2.0 * gam / (beta * m_aux_var[0]));
1177+
xi = rnd.normal();
1178+
theta = rnd.normal();
1179+
p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
1180+
+ 0.5 * (f_aux_var[0] - gam * p_aux_var[0]) * dt
1181+
- 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gam * sigma * m_aux_var[0] * sqrt(pow(dt,3))
1182+
- 0.125 * gam * (f_aux_var[0] - gam * p_aux_var[0]) * pow(dt,2);
1183+
y_aux_var[0] += p_aux_var[0] / m_aux_var[0] * dt + 0.5 * theta / sqrt(3.0) * sigma * sqrt(pow(dt,3));
1184+
f_aux_var[0] = (Ham.kcrpmd_effective_auxiliary_force(y_aux_var, beta, eta, a, b, c, d))[0];
1185+
p_aux_var[0] += 0.5 * xi * sigma * m_aux_var[0] * sqrt(dt)
1186+
+ 0.5 * (f_aux_var[0] - gam * p_aux_var[0]) * dt
1187+
- 0.25 * (xi / 2.0 + theta / sqrt(3.0)) * gam * sigma * m_aux_var[0] * sqrt(pow(dt,3))
1188+
- 0.125 * gam * (f_aux_var[0] - gam * p_aux_var[0]) * pow(dt,2);
1189+
}
1190+
}
1191+
}
1192+
11231193
void renormalize_hopping_probabilities(vector< vector<double> >& g, vector< vector<vector<double> > >& coherence_factors, vector<int>& act_states){
11241194
/**
11251195
Scale the hopping probabilities by the factors and ensure the sum of all
@@ -1268,17 +1338,49 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
12681338

12691339
//***************************** Coherent dynamics *******************************
12701340
//============== Nuclear propagation ===================
1341+
// // NVT dynamics
1342+
// if(prms.ensemble==1){
1343+
// for(idof=0; idof<n_therm_dofs; idof++){
1344+
// dof = prms.thermostat_dofs[idof];
1345+
// for(traj=0; traj<ntraj; traj++){
1346+
// dyn_var.p->scale(dof, traj, therm[traj].vel_scale(0.5*prms.dt));
1347+
// }// traj
1348+
// }// idof
1349+
// }
1350+
//
1351+
// *dyn_var.p = *dyn_var.p + 0.5 * prms.dt * (*dyn_var.f);
1352+
12711353
// NVT dynamics
1272-
if(prms.ensemble==1){
1273-
for(idof=0; idof<n_therm_dofs; idof++){
1274-
dof = prms.thermostat_dofs[idof];
1275-
for(traj=0; traj<ntraj; traj++){
1276-
dyn_var.p->scale(dof, traj, therm[traj].vel_scale(0.5*prms.dt));
1277-
}// traj
1278-
}// idof
1354+
double beta = (hartree/boltzmann) / prms.Temperature;
1355+
MATRIX sigma_therm(ndof,ntraj);
1356+
MATRIX xi_therm(ndof, ntraj);
1357+
MATRIX theta_therm(ndof, ntraj);
1358+
MATRIX extra_force(ndof, ntraj);
1359+
MATRIX extra_momentum(ndof, ntraj);
1360+
if(prms.ensemble==1){
1361+
if(therm[0].thermostat_type=="Nose-Hoover"){
1362+
for(idof=0; idof<n_therm_dofs; idof++){
1363+
dof = prms.thermostat_dofs[idof];
1364+
for(traj=0; traj<ntraj; traj++){
1365+
dyn_var.p->scale(dof, traj, therm[traj].vel_scale(0.5*prms.dt));
1366+
}// traj
1367+
}// idof
1368+
}// Nose-Hoover
1369+
else if(therm[0].thermostat_type=="Langevin"){
1370+
for(idof=0; idof<n_therm_dofs; idof++){
1371+
dof = prms.thermostat_dofs[idof];
1372+
for(traj=0; traj<ntraj; traj++){
1373+
sigma_therm.set(dof, traj, sqrt(2 * invM.get(dof,traj) * therm[traj].nu_therm / beta));
1374+
xi_therm.set(dof, traj, rnd.normal());
1375+
theta_therm.set(dof, traj, rnd.normal());
1376+
extra_force.set(dof, traj, 0.5 * xi_therm.get(dof,traj) * sigma_therm.get(dof,traj) / invM.get(dof,traj) * sqrt(prms.dt) - 0.5 * therm[traj].nu_therm * dyn_var.p->get(dof,traj) * prms.dt - 0.25 * (xi_therm.get(dof,traj) / 2.0 + theta_therm.get(dof,traj) / sqrt(3.0)) * therm[traj].nu_therm * sigma_therm.get(dof,traj) / invM.get(dof,traj) * sqrt(pow(prms.dt,3)) - 0.125 * therm[traj].nu_therm * (dyn_var.f->get(dof,traj) - therm[traj].nu_therm * dyn_var.p->get(dof,traj)) * pow(prms.dt,2));
1377+
extra_momentum.set(dof, traj, theta_therm.get(dof,traj) * sigma_therm.get(dof,traj) / (2.0 * sqrt(3.0)) * sqrt(pow(prms.dt,3)));
1378+
}// traj
1379+
}// idof
1380+
}// Langevin
12791381
}
12801382

1281-
*dyn_var.p = *dyn_var.p + 0.5 * prms.dt * (*dyn_var.f);
1383+
*dyn_var.p = *dyn_var.p + 0.5 * prms.dt * (*dyn_var.f) + extra_force;
12821384
//if(prms.use_qtsh==1){ *dyn_var.p = *dyn_var.p + 0.5 * compute_dkinemat(dyn_var, ham); }
12831385
if(prms.use_qtsh==1){ *dyn_var.p = *dyn_var.p + 0.5 * prms.dt * (*dyn_var.qtsh_f_nc); }
12841386

@@ -1293,7 +1395,7 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
12931395
// Update coordinates of nuclei for all trajectories
12941396
for(traj=0; traj<ntraj; traj++){
12951397
for(dof=0; dof<ndof; dof++){
1296-
dyn_var.q->add(dof, traj, invM.get(dof,0) * dyn_var.p->get(dof,traj) * prms.dt );
1398+
dyn_var.q->add(dof, traj, invM.get(dof,0) * dyn_var.p->get(dof,traj) * prms.dt + extra_momentum.get(dof,traj) );
12971399

12981400
if(prms.entanglement_opt==22){
12991401
dyn_var.q->add(dof, traj, invM.get(dof,0) * gamma.get(dof,traj) * prms.dt );
@@ -1344,8 +1446,10 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
13441446

13451447
// Recompute density matrices in response to the updated amplitudes
13461448
dyn_var.update_amplitudes(prms);
1347-
dyn_var.update_density_matrix(prms);
1348-
1449+
dyn_var.update_density_matrix(prms);
1450+
1451+
if(prms.use_kcrpmd==1){ propagate_electronic_kcrpmd(dyn_var, ham, prms, rnd); }
1452+
13491453
vector<int> old_states( dyn_var.act_states);
13501454

13511455
// In the interval [t, t + dt], we may have experienced the basis reordering, so we need to
@@ -1366,18 +1470,40 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
13661470
update_forces_xf(dyn_var, ham, ham_aux);
13671471
}
13681472

1369-
// NVT dynamics
1370-
if(prms.ensemble==1){
1371-
for(traj=0; traj<ntraj; traj++){
1372-
t2[0] = traj;
1373-
pop_submatrix(p, p_traj, t1, t2);
1374-
double ekin = compute_kinetic_energy(p_traj, invM, prms.thermostat_dofs);
1375-
therm[traj].propagate_nhc(prms.dt, ekin, 0.0, 0.0);
1376-
}
1473+
// // NVT dynamics
1474+
// if(prms.ensemble==1){
1475+
// for(traj=0; traj<ntraj; traj++){
1476+
// t2[0] = traj;
1477+
// pop_submatrix(p, p_traj, t1, t2);
1478+
// double ekin = compute_kinetic_energy(p_traj, invM, prms.thermostat_dofs);
1479+
// therm[traj].propagate_nhc(prms.dt, ekin, 0.0, 0.0);
1480+
// }
1481+
//
1482+
// }
1483+
//
1484+
// *dyn_var.p = *dyn_var.p + 0.5*prms.dt* (*dyn_var.f);
13771485

1486+
// NVT dynamics
1487+
if(prms.ensemble==1){
1488+
if(therm[0].thermostat_type=="Nose-Hoover"){
1489+
for(traj=0; traj<ntraj; traj++){
1490+
t2[0] = traj;
1491+
pop_submatrix(p, p_traj, t1, t2);
1492+
double ekin = compute_kinetic_energy(p_traj, invM, prms.thermostat_dofs);
1493+
therm[traj].propagate_nhc(prms.dt, ekin, 0.0, 0.0);
1494+
}// traj
1495+
}// Nose-Hoover
1496+
else if(therm[0].thermostat_type=="Langevin"){
1497+
for(idof=0; idof<n_therm_dofs; idof++){
1498+
dof = prms.thermostat_dofs[idof];
1499+
for(traj=0; traj<ntraj; traj++){
1500+
extra_force.set(dof, traj, 0.5 * xi_therm.get(dof,traj) * sigma_therm.get(dof,traj) / invM.get(dof,traj) * sqrt(prms.dt) - 0.5 * therm[traj].nu_therm * dyn_var.p->get(dof,traj) * prms.dt - 0.25 * (xi_therm.get(dof,traj) / 2.0 + theta_therm.get(dof,traj) / sqrt(3.0)) * therm[traj].nu_therm * sigma_therm.get(dof,traj) / invM.get(dof,traj) * sqrt(pow(prms.dt,3)) - 0.125 * therm[traj].nu_therm * (dyn_var.f->get(dof,traj) - therm[traj].nu_therm * dyn_var.p->get(dof,traj)) * pow(prms.dt,2));
1501+
}// traj
1502+
}// idof
1503+
}// Langevin
13781504
}
13791505

1380-
*dyn_var.p = *dyn_var.p + 0.5*prms.dt* (*dyn_var.f);
1506+
*dyn_var.p = *dyn_var.p + 0.5 * prms.dt * (*dyn_var.f) + extra_force;
13811507
//if(prms.use_qtsh==1){ *dyn_var.p = *dyn_var.p + 0.5* compute_dkinemat(dyn_var, ham); }
13821508
if(prms.use_qtsh==1){ *dyn_var.p = *dyn_var.p + 0.5 * prms.dt * (*dyn_var.qtsh_f_nc); }
13831509

@@ -1386,16 +1512,29 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
13861512
dyn_var.p->scale(prms.constrained_dofs[cdof], -1, 0.0);
13871513
}
13881514

1515+
// // NVT dynamics
1516+
// if(prms.ensemble==1){
1517+
// for(idof=0; idof<n_therm_dofs; idof++){
1518+
// dof = prms.thermostat_dofs[idof];
1519+
// for(traj=0; traj<ntraj; traj++){
1520+
// dyn_var.p->scale(dof, traj, therm[traj].vel_scale(0.5*prms.dt));
1521+
// }// traj
1522+
// }// idof
1523+
// }
1524+
13891525
// NVT dynamics
1390-
if(prms.ensemble==1){
1391-
for(idof=0; idof<n_therm_dofs; idof++){
1392-
dof = prms.thermostat_dofs[idof];
1393-
for(traj=0; traj<ntraj; traj++){
1394-
dyn_var.p->scale(dof, traj, therm[traj].vel_scale(0.5*prms.dt));
1395-
}// traj
1396-
}// idof
1526+
if(prms.ensemble==1){
1527+
if(therm[0].thermostat_type=="Nose-Hoover"){
1528+
for(idof=0; idof<n_therm_dofs; idof++){
1529+
dof = prms.thermostat_dofs[idof];
1530+
for(traj=0; traj<ntraj; traj++){
1531+
dyn_var.p->scale(dof, traj, therm[traj].vel_scale(0.5*prms.dt));
1532+
}// traj
1533+
}// idof
1534+
}// Nose-Hoover
13971535
}
13981536

1537+
13991538
//ham_aux.copy_content(ham);
14001539
update_Hamiltonian_variables(prms, dyn_var, ham, ham_aux, py_funct, params, 1);
14011540

@@ -1635,7 +1774,10 @@ void compute_dynamics(dyn_variables& dyn_var, bp::dict dyn_params,
16351774
// Saves the current density matrix into the previous - needed for FSSH2 and GFSH (original)
16361775
dyn_var.save_curr_dm_into_prev();
16371776

1638-
1777+
for(traj=0; traj<ntraj; traj++){
1778+
*dyn_var.ave_decoherence_rates += decoherence_rates[traj];
1779+
}
1780+
*dyn_var.ave_decoherence_rates /= ntraj;
16391781

16401782
}
16411783

‎src/dyn/Dynamics.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -74,6 +74,7 @@ int need_active_states_diff_rep(dyn_control_params& prms);
7474

7575
void propagate_electronic(dyn_variables& dyn_var, nHamiltonian& ham, nHamiltonian& ham_prev, dyn_control_params& prms);
7676
void propagate_electronic(dyn_variables& dyn_var, nHamiltonian* ham, nHamiltonian* ham_prev, dyn_control_params& prms);
77+
void propagate_electronic_kcrpmd(dyn_variables& dyn_var, nHamiltonian& ham, dyn_control_params& prms, Random& rnd);
7778

7879
/*
7980
void compute_dynamics(MATRIX& q, MATRIX& p, MATRIX& invM, CMATRIX& C, vector<CMATRIX>& projectors, vector<int>& act_states,

‎src/dyn/Energy_and_Forces.cpp‎

Lines changed: 38 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -278,6 +278,29 @@ vector<double> potential_energies(dyn_control_params& prms, dyn_variables& dyn_v
278278
}// rep_force == 1
279279

280280
}
281+
else if(prms.force_method==4){ // KC-RPMD force
282+
283+
// Diabatic
284+
if(prms.rep_force==0){
285+
double Veff;
286+
double beta = (hartree/boltzmann) / prms.Temperature;
287+
double eta = prms.kcrpmd_eta;
288+
double a = prms.kcrpmd_a;
289+
double b = prms.kcrpmd_b;
290+
double c = prms.kcrpmd_c;
291+
double d = prms.kcrpmd_d;
292+
Veff = ham.kcrpmd_effective_potential(dyn_vars.y_aux_var, *dyn_vars.q, *dyn_vars.iM, beta, eta, a, b, c, d);
293+
for(itraj=0; itraj<ntraj; itraj++){
294+
id[1] = itraj;
295+
res[itraj] = Veff;
296+
}
297+
}
298+
// Adiabatic
299+
else if(prms.rep_force==1){
300+
cout<<"Error in potential_energies(): KC-RPMD works exclusively within diabatic bases\n"; exit(0);
301+
}
302+
303+
}// KC-RPMD
281304

282305

283306
return res;
@@ -407,6 +430,21 @@ void update_forces(dyn_control_params& prms, dyn_variables& dyn_vars, nHamiltoni
407430
}
408431

409432
}// QTSH
433+
else if(prms.force_method==4){ // KC-RPMD forces
434+
// Diabatic
435+
if(prms.rep_force==0){
436+
double beta = (hartree/boltzmann) / prms.Temperature;
437+
double eta = prms.kcrpmd_eta;
438+
double a = prms.kcrpmd_a;
439+
double b = prms.kcrpmd_b;
440+
double c = prms.kcrpmd_c;
441+
double d = prms.kcrpmd_d;
442+
*dyn_vars.f = ham.kcrpmd_effective_force(dyn_vars.y_aux_var, *dyn_vars.q, *dyn_vars.iM, beta, eta, a, b, c, d);
443+
}
444+
445+
// Adiabatic
446+
else if(prms.rep_force==1){ cout<<"Error in update_forces(): KC-RPMD works exclusively within diabatic bases\n"; exit(0); }
447+
}// KC-RPMD
410448

411449
// return F;
412450
}

‎src/dyn/Energy_and_Forces.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -31,6 +31,7 @@
3131
//#include "ensemble/libensemble.h"
3232
#include "dyn_control_params.h"
3333
#include "dyn_variables.h"
34+
#include "../Units.h"
3435

3536

3637
/// liblibra namespace

‎src/dyn/dyn_control_params.cpp‎

Lines changed: 27 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -86,6 +86,13 @@ dyn_control_params::dyn_control_params(){
8686

8787
///================= KCRPMD specific ====================
8888
use_kcrpmd = 0;
89+
kcrpmd_eta = 6.28;
90+
kcrpmd_a = 0.1;
91+
kcrpmd_b = 1000.;
92+
kcrpmd_c = 0.5;
93+
kcrpmd_d = 3.0;
94+
kcrpmd_gamma = 0.0;
95+
kcrpmd_gammaKP = 0.0;
8996

9097
///================= Decoherence options =========================================
9198
decoherence_algo = -1;
@@ -197,6 +204,16 @@ dyn_control_params::dyn_control_params(const dyn_control_params& x){
197204
use_qtsh = x.use_qtsh;
198205
qtsh_force_option = x.qtsh_force_option;
199206

207+
///================= KC-RPMD specific ====================
208+
use_kcrpmd = x.use_kcrpmd;
209+
kcrpmd_eta = x.kcrpmd_eta;
210+
kcrpmd_a = x.kcrpmd_a;
211+
kcrpmd_b = x.kcrpmd_b;
212+
kcrpmd_c = x.kcrpmd_c;
213+
kcrpmd_d = x.kcrpmd_d;
214+
kcrpmd_gamma = x.kcrpmd_gamma;
215+
kcrpmd_gammaKP = x.kcrpmd_gammaKP;
216+
200217
///================= Decoherence options =========================================
201218
decoherence_algo = x.decoherence_algo;
202219
sdm_norm_tolerance = x.sdm_norm_tolerance;
@@ -399,6 +416,16 @@ void dyn_control_params::set_parameters(bp::dict params){
399416
else if(key=="use_qtsh"){ use_qtsh = bp::extract<int>(params.values()[i]); }
400417
else if(key=="qtsh_force_option"){ qtsh_force_option = bp::extract<int>(params.values()[i]); }
401418

419+
///================= KC-RPMD specific ====================
420+
else if(key=="use_kcrpmd"){ use_kcrpmd = bp::extract<int>(params.values()[i]); }
421+
else if(key=="kcrpmd_eta"){ kcrpmd_eta = bp::extract<double>(params.values()[i]); }
422+
else if(key=="kcrpmd_a"){ kcrpmd_a = bp::extract<double>(params.values()[i]); }
423+
else if(key=="kcrpmd_b"){ kcrpmd_b = bp::extract<double>(params.values()[i]); }
424+
else if(key=="kcrpmd_c"){ kcrpmd_c = bp::extract<double>(params.values()[i]); }
425+
else if(key=="kcrpmd_d"){ kcrpmd_d = bp::extract<double>(params.values()[i]); }
426+
else if(key=="kcrpmd_gamma"){ kcrpmd_gamma = bp::extract<double>(params.values()[i]); }
427+
else if(key=="kcrpmd_gammaKP"){ kcrpmd_gammaKP = bp::extract<double>(params.values()[i]); }
428+
402429
///================= Decoherence options =========================================
403430
else if(key=="decoherence_algo"){ decoherence_algo = bp::extract<int>(params.values()[i]); }
404431
else if(key=="sdm_norm_tolerance"){ sdm_norm_tolerance = bp::extract<double>(params.values()[i]); }

0 commit comments

Comments
 (0)