version 1.29 | version 1.30 |
---|
| |
get_keyval(conf, "width", width, 1.0); | get_keyval(conf, "width", width, 1.0); |
if (width <= 0.0) { | if (width <= 0.0) { |
cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); | cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); |
| return; |
} | } |
| |
// NOTE: not porting wall stuff to new deps, as this will change to a separate bias | // NOTE: not porting wall stuff to new deps, as this will change to a separate bias |
| |
return error_code; | return error_code; |
} | } |
| |
| if (cvm::step_relative() > 0) { |
| // Total force depends on Jacobian derivative from previous timestep |
| error_code |= calc_cvc_total_force(first_cvc, num_cvcs); |
| } |
| // atom coordinates are updated by the next line |
error_code |= calc_cvc_values(first_cvc, num_cvcs); | error_code |= calc_cvc_values(first_cvc, num_cvcs); |
error_code |= calc_cvc_gradients(first_cvc, num_cvcs); | error_code |= calc_cvc_gradients(first_cvc, num_cvcs); |
error_code |= calc_cvc_total_force(first_cvc, num_cvcs); | |
error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs); | error_code |= calc_cvc_Jacobians(first_cvc, num_cvcs); |
| |
if (cvm::debug()) | if (cvm::debug()) |
| |
| |
int error_code = COLVARS_OK; | int error_code = COLVARS_OK; |
| |
| if (cvm::step_relative() > 0) { |
| // Total force depends on Jacobian derivative from previous timestep |
| error_code |= collect_cvc_total_forces(); |
| } |
error_code |= collect_cvc_values(); | error_code |= collect_cvc_values(); |
error_code |= collect_cvc_gradients(); | error_code |= collect_cvc_gradients(); |
error_code |= collect_cvc_total_forces(); | |
error_code |= collect_cvc_Jacobians(); | error_code |= collect_cvc_Jacobians(); |
error_code |= calc_colvar_properties(); | error_code |= calc_colvar_properties(); |
| |
| |
// if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) { | // if (!tasks[task_extended_lagrangian] && (cvm::step_relative() > 0)) { |
// Disabled check to allow for explicit total force calculation | // Disabled check to allow for explicit total force calculation |
// even with extended Lagrangian | // even with extended Lagrangian |
| |
if (cvm::step_relative() > 0) { | |
cvm::increase_depth(); | cvm::increase_depth(); |
// get from the cvcs the total forces from the PREVIOUS step | |
for (i = first_cvc, cvc_count = 0; | for (i = first_cvc, cvc_count = 0; |
(i < cvcs.size()) && (cvc_count < cvc_max_count); | (i < cvcs.size()) && (cvc_count < cvc_max_count); |
i++) { | i++) { |
| |
(cvcs[i])->calc_force_invgrads(); | (cvcs[i])->calc_force_invgrads(); |
} | } |
cvm::decrease_depth(); | cvm::decrease_depth(); |
} | |
| |
if (cvm::debug()) | if (cvm::debug()) |
cvm::log("Done calculating total force of colvar \""+this->name+"\".\n"); | cvm::log("Done calculating total force of colvar \""+this->name+"\".\n"); |
| |
if (is_enabled(f_cv_Jacobian)) { | if (is_enabled(f_cv_Jacobian)) { |
// the instantaneous Jacobian force was not included in the reported total force; | // the instantaneous Jacobian force was not included in the reported total force; |
// instead, it is subtracted from the applied force (silent Jacobian correction) | // instead, it is subtracted from the applied force (silent Jacobian correction) |
| // This requires the Jacobian term for the *current* timestep |
if (is_enabled(f_cv_hide_Jacobian)) | if (is_enabled(f_cv_hide_Jacobian)) |
f -= fj; | f -= fj; |
} | } |
| |
if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) { | |
| |
// Wall force | // Wall force |
colvarvalue fw(x); | colvarvalue fw(x); |
fw.reset(); | fw.reset(); |
| |
| if (is_enabled(f_cv_lower_wall) || is_enabled(f_cv_upper_wall)) { |
| |
if (cvm::debug()) | if (cvm::debug()) |
cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n"); | cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n"); |
| |
// For a periodic colvar, both walls may be applicable at the same time | // For a periodic colvar, both walls may be applicable at the same time |
// in which case we pick the closer one | // in which case we pick the closer one |
if ( (!is_enabled(f_cv_upper_wall)) || | if ( (!is_enabled(f_cv_upper_wall)) || |
(this->dist2(x_reported, lower_wall) < this->dist2(x_reported, upper_wall)) ) { | (this->dist2(x, lower_wall) < this->dist2(x, upper_wall)) ) { |
| |
cvm::real const grad = this->dist2_lgrad(x_reported, lower_wall); | cvm::real const grad = this->dist2_lgrad(x, lower_wall); |
if (grad < 0.0) { | if (grad < 0.0) { |
fw = -0.5 * lower_wall_k * grad; | fw = -0.5 * lower_wall_k * grad; |
f += fw; | |
if (cvm::debug()) | if (cvm::debug()) |
cvm::log("Applying a lower wall force("+ | cvm::log("Applying a lower wall force("+ |
cvm::to_str(fw)+") to \""+this->name+"\".\n"); | cvm::to_str(fw)+") to \""+this->name+"\".\n"); |
| |
| |
} else { | } else { |
| |
cvm::real const grad = this->dist2_lgrad(x_reported, upper_wall); | cvm::real const grad = this->dist2_lgrad(x, upper_wall); |
if (grad > 0.0) { | if (grad > 0.0) { |
fw = -0.5 * upper_wall_k * grad; | fw = -0.5 * upper_wall_k * grad; |
f += fw; | |
if (cvm::debug()) | if (cvm::debug()) |
cvm::log("Applying an upper wall force("+ | cvm::log("Applying an upper wall force("+ |
cvm::to_str(fw)+") to \""+this->name+"\".\n"); | cvm::to_str(fw)+") to \""+this->name+"\".\n"); |
| |
} | } |
} | } |
| |
| // At this point f is the force f from external biases that will be applied to the |
| // extended variable if there is one |
| |
if (is_enabled(f_cv_extended_Lagrangian)) { | if (is_enabled(f_cv_extended_Lagrangian)) { |
| |
cvm::real dt = cvm::dt(); | cvm::real dt = cvm::dt(); |
| |
f_ext.reset(); | f_ext.reset(); |
| |
// the total force is applied to the fictitious mass, while the | // the total force is applied to the fictitious mass, while the |
// atoms only feel the harmonic force | // atoms only feel the harmonic force + wall force |
// fr: bias force on extended variable (without harmonic spring), for output in trajectory | // fr: bias force on extended variable (without harmonic spring), for output in trajectory |
// f_ext: total force on extended variable (including harmonic spring) | // f_ext: total force on extended variable (including harmonic spring) |
// f: - initially, external biasing force (including wall forces) | // f: - initially, external biasing force |
// - after this code block, colvar force to be applied to atomic coordinates, ie. spring force | // - after this code block, colvar force to be applied to atomic coordinates |
| // ie. spring force + wall force |
fr = f; | fr = f; |
f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); | f_ext = f + (-0.5 * ext_force_k) * this->dist2_lgrad(xr, x); |
f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x); | f = (-0.5 * ext_force_k) * this->dist2_rgrad(xr, x); |
| |
vr += (0.5 * dt) * f_ext / ext_mass; | vr += (0.5 * dt) * f_ext / ext_mass; |
xr += dt * vr; | xr += dt * vr; |
xr.apply_constraints(); | xr.apply_constraints(); |
if (this->b_periodic) this->wrap(xr); | if (this->is_enabled(f_cv_periodic)) this->wrap(xr); |
} | } |
| |
| f += fw; |
| |
| // Store force to be applied, possibly summed over several timesteps |
f_accumulated += f; | f_accumulated += f; |
| |
if (is_enabled(f_cv_fdiff_velocity)) { | if (is_enabled(f_cv_fdiff_velocity)) { |