Skip to content
Permalink
Browse files

core: thermostat: Removed head_up/cool_down mechanism for Langevin ty…

…pe thermostats
  • Loading branch information
fweik committed Nov 26, 2019
1 parent 9415194 commit 2e4f50d048721251c3563bc07f9949b49085a129
@@ -106,16 +106,6 @@ int thermalized_bond_set_params(int bond_type, double temp_com,
return ES_OK;
}

void thermalized_bond_heat_up() {
double pref_scale = sqrt(3);
thermalized_bond_update_params(pref_scale);
}

void thermalized_bond_cool_down() {
double pref_scale = 1.0 / sqrt(3);
thermalized_bond_update_params(pref_scale);
}

void thermalized_bond_init() {

for (auto &bonded_ia_param : bonded_ia_params) {
@@ -59,9 +59,6 @@ bool thermalized_bond_is_seed_required();
uint64_t thermalized_bond_get_rng_state();
void thermalized_bond_set_rng_state(uint64_t counter);

/* Setup */
void thermalized_bond_heat_up();
void thermalized_bond_cool_down();
void thermalized_bond_update_params(double pref_scale);
void thermalized_bond_init();

@@ -72,7 +69,6 @@ void thermalized_bond_init();
2. Salt (decorrelates different counter)
3. Particle ID, Partner ID
*/

inline Utils::Vector3d v_noise(int particle_id, int partner_id) {

using rng_type = r123::Philox4x64;
@@ -103,16 +103,6 @@ void dpd_set_rng_state(const uint64_t counter) {

uint64_t dpd_get_rng_state() { return dpd_rng_counter->value(); }

void dpd_heat_up() {
double pref_scale = sqrt(3);
dpd_update_params(pref_scale);
}

void dpd_cool_down() {
double pref_scale = 1.0 / sqrt(3);
dpd_update_params(pref_scale);
}

int dpd_set_params(int part_type_a, int part_type_b, double gamma, double r_c,
int wf, double tgamma, double tr_c, int twf) {
IA_parameters &ia_params = *get_ia_param_safe(part_type_a, part_type_b);
@@ -43,8 +43,6 @@ struct DPDParameters {
double pref = 0.0;
};

void dpd_heat_up();
void dpd_cool_down();
int dpd_set_params(int part_type_a, int part_type_b, double gamma, double r_c,
int wf, double tgamma, double tr_c, int twf);
void dpd_init();
@@ -178,8 +178,6 @@ void integrate_vv(int n_steps, int reuse_forces) {
*/
if (reuse_forces == -1 || (recalc_forces && reuse_forces != 1)) {
ESPRESSO_PROFILER_MARK_BEGIN("Initial Force Calculation");
thermo_heat_up();

lb_lbcoupling_deactivate();

#ifdef VIRTUAL_SITES
@@ -197,8 +195,6 @@ void integrate_vv(int n_steps, int reuse_forces) {
#endif
}

thermo_cool_down();

ESPRESSO_PROFILER_MARK_END("Initial Force Calculation");
}

@@ -161,44 +161,3 @@ void thermo_init() {
thermo_init_npt_isotropic();
#endif
}

void langevin_heat_up() {
langevin_pref2_buffer = langevin_pref2;
langevin_pref2 *= sqrt(3);

langevin_pref2_rotation_buffer = langevin_pref2_rotation;
langevin_pref2_rotation *= sqrt(3);
}

void thermo_heat_up() {
if (thermo_switch & THERMO_LANGEVIN) {
langevin_heat_up();
}
#ifdef DPD
if (thermo_switch & THERMO_DPD) {
dpd_heat_up();
}
#endif
if (n_thermalized_bonds) {
thermalized_bond_heat_up();
}
}

void langevin_cool_down() {
langevin_pref2 = langevin_pref2_buffer;
langevin_pref2_rotation = langevin_pref2_rotation_buffer;
}

void thermo_cool_down() {
if (thermo_switch & THERMO_LANGEVIN) {
langevin_cool_down();
}
#ifdef DPD
if (thermo_switch & THERMO_DPD) {
dpd_cool_down();
}
#endif
if (n_thermalized_bonds) {
thermalized_bond_cool_down();
}
}
@@ -108,18 +108,6 @@ uint64_t langevin_get_rng_state();
start of integration */
void thermo_init();

/** very nasty: if we recalculate force when leaving/reentering the integrator,
a(t) and a((t-dt)+dt) are NOT equal in the vv algorithm. The random numbers
are drawn twice, resulting in a different variance of the random force.
This is corrected by additional heat when restarting the integrator here.
Currently only works for the Langevin thermostat, although probably also
others are affected.
*/
void thermo_heat_up();

/** pendant to \ref thermo_heat_up */
void thermo_cool_down();

#ifdef NPT
/** add velocity-dependent noise and friction for NpT-sims to the particle's
velocity

0 comments on commit 2e4f50d

Please sign in to comment.
You can’t perform that action at this time.