version 1.4 | version 1.5 |
---|
| |
// -*- c++ -*- | // -*- c++ -*- |
| |
| // This file is part of the Collective Variables module (Colvars). |
| // The original version of Colvars and its updates are located at: |
| // https://github.com/colvars/colvars |
| // Please update all Colvars source files before making any changes. |
| // If you wish to distribute your changes, please submit them to the |
| // Colvars repository at GitHub. |
| |
#include "colvarmodule.h" | #include "colvarmodule.h" |
#include "colvarvalue.h" | #include "colvarvalue.h" |
#include "colvarbias_restraint.h" | #include "colvarbias_restraint.h" |
| |
| |
| |
colvarbias_restraint::colvarbias_restraint(char const *key) | colvarbias_restraint::colvarbias_restraint(char const *key) |
: colvarbias(key) | : colvarbias(key) |
{ | { |
| |
int colvarbias_restraint::init(std::string const &conf) | int colvarbias_restraint::init(std::string const &conf) |
{ | { |
colvarbias::init(conf); | colvarbias::init(conf); |
| enable(f_cvb_apply_force); |
| |
if (cvm::debug()) | if (cvm::debug()) |
cvm::log("Initializing a new restraint bias.\n"); | cvm::log("Initializing a new restraint bias.\n"); |
| |
// TODO move these initializations to constructor and let get_keyval | return COLVARS_OK; |
// only override existing values | } |
target_nstages = 0; | |
target_nsteps = 0; | |
force_k = 0.0; | |
| |
get_keyval(conf, "forceConstant", force_k, 1.0); | |
| |
| int colvarbias_restraint::update() |
{ | { |
// get the initial restraint centers | // Update base class (bias_energy and colvar_forces are zeroed there) |
colvar_centers.resize(colvars.size()); | colvarbias::update(); |
colvar_centers_raw.resize(colvars.size()); | |
size_t i; | |
| |
enable(f_cvb_apply_force); | // Force and energy calculation |
| for (size_t i = 0; i < num_variables(); i++) { |
| bias_energy += restraint_potential(i); |
| colvar_forces[i].type(variables(i)->value()); |
| colvar_forces[i].is_derivative(); |
| colvar_forces[i] = restraint_force(i); |
| } |
| |
for (i = 0; i < colvars.size(); i++) { | if (cvm::debug()) |
colvar_centers[i].type(colvars[i]->value()); | cvm::log("Done updating the restraint bias \""+this->name+"\".\n"); |
colvar_centers_raw[i].type(colvars[i]->value()); | |
if (cvm::debug()) { | if (cvm::debug()) |
cvm::log("colvarbias_restraint: center size = "+ | cvm::log("Current forces for the restraint bias \""+ |
cvm::to_str(colvar_centers[i].vector1d_value.size())+"\n"); | this->name+"\": "+cvm::to_str(colvar_forces)+".\n"); |
| |
| return COLVARS_OK; |
| } |
| |
| |
| colvarbias_restraint::~colvarbias_restraint() |
| { |
| } |
| |
| |
| std::string const colvarbias_restraint::get_state_params() const |
| { |
| return colvarbias::get_state_params(); |
| } |
| |
| |
| int colvarbias_restraint::set_state_params(std::string const &conf) |
| { |
| return colvarbias::set_state_params(conf); |
| } |
| |
| |
| std::ostream & colvarbias_restraint::write_traj_label(std::ostream &os) |
| { |
| return colvarbias::write_traj_label(os); |
| } |
| |
| |
| std::ostream & colvarbias_restraint::write_traj(std::ostream &os) |
| { |
| return colvarbias::write_traj(os); |
} | } |
| |
| |
| |
| colvarbias_restraint_centers::colvarbias_restraint_centers(char const *key) |
| : colvarbias(key), colvarbias_restraint(key) |
| { |
} | } |
| |
| |
| int colvarbias_restraint_centers::init(std::string const &conf) |
| { |
| size_t i; |
| |
| bool null_centers = (colvar_centers.size() == 0); |
| if (null_centers) { |
| // try to initialize the restraint centers for the first time |
| colvar_centers.resize(num_variables()); |
| colvar_centers_raw.resize(num_variables()); |
| for (i = 0; i < num_variables(); i++) { |
| colvar_centers[i].type(variables(i)->value()); |
| colvar_centers[i].reset(); |
| colvar_centers_raw[i].type(variables(i)->value()); |
| colvar_centers_raw[i].reset(); |
| } |
| } |
| |
if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { | if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { |
for (i = 0; i < colvars.size(); i++) { | for (i = 0; i < num_variables(); i++) { |
if (cvm::debug()) { | if (cvm::debug()) { |
cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n"); | cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n"); |
} | } |
| |
colvar_centers[i].apply_constraints(); | |
colvar_centers_raw[i] = colvar_centers[i]; | colvar_centers_raw[i] = colvar_centers[i]; |
| colvar_centers[i].apply_constraints(); |
} | } |
} else { | null_centers = false; |
| } |
| |
| if (null_centers) { |
colvar_centers.clear(); | colvar_centers.clear(); |
cvm::error("Error: must define the initial centers of the restraints.\n"); | cvm::error("Error: must define the initial centers of the restraints.\n", INPUT_ERROR); |
| return INPUT_ERROR; |
} | } |
| |
if (colvar_centers.size() != colvars.size()) { | if (colvar_centers.size() != num_variables()) { |
cvm::error("Error: number of centers does not match " | cvm::error("Error: number of centers does not match " |
"that of collective variables.\n"); | "that of collective variables.\n", INPUT_ERROR); |
} | return INPUT_ERROR; |
} | } |
| |
{ | return COLVARS_OK; |
if (cvm::debug()) { | |
cvm::log("colvarbias_restraint: parsing target centers.\n"); | |
} | } |
| |
size_t i; | |
if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) { | int colvarbias_restraint_centers::change_configuration(std::string const &conf) |
if (colvar_centers.size() != colvars.size()) { | { |
cvm::error("Error: number of target centers does not match " | if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { |
"that of collective variables.\n"); | for (size_t i = 0; i < num_variables(); i++) { |
} | colvar_centers[i].type(variables(i)->value()); |
b_chg_centers = true; | colvar_centers[i].apply_constraints(); |
for (i = 0; i < target_centers.size(); i++) { | colvar_centers_raw[i].type(variables(i)->value()); |
target_centers[i].apply_constraints(); | colvar_centers_raw[i] = colvar_centers[i]; |
} | } |
} else { | |
b_chg_centers = false; | |
target_centers.clear(); | |
} | } |
| return COLVARS_OK; |
} | } |
| |
if (get_keyval(conf, "targetForceConstant", target_force_k, 0.0)) { | |
if (b_chg_centers) | |
cvm::error("Error: cannot specify both targetCenters and targetForceConstant.\n"); | |
| |
starting_force_k = force_k; | |
b_chg_force_k = true; | |
| |
get_keyval(conf, "targetEquilSteps", target_equil_steps, 0); | colvarbias_restraint_k::colvarbias_restraint_k(char const *key) |
| : colvarbias(key), colvarbias_restraint(key) |
| { |
| force_k = -1.0; |
| } |
| |
get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule); | |
if (lambda_schedule.size()) { | int colvarbias_restraint_k::init(std::string const &conf) |
// There is one more lambda-point than stages | { |
target_nstages = lambda_schedule.size() - 1; | get_keyval(conf, "forceConstant", force_k, (force_k > 0.0 ? force_k : 1.0)); |
| if (force_k < 0.0) { |
| cvm::error("Error: undefined or invalid force constant.\n", INPUT_ERROR); |
| return INPUT_ERROR; |
} | } |
} else { | return COLVARS_OK; |
b_chg_force_k = false; | |
} | } |
| |
if (b_chg_centers || b_chg_force_k) { | |
get_keyval(conf, "targetNumSteps", target_nsteps, 0); | |
if (!target_nsteps) | |
cvm::error("Error: targetNumSteps must be non-zero.\n"); | |
| |
if (get_keyval(conf, "targetNumStages", target_nstages, target_nstages) && | int colvarbias_restraint_k::change_configuration(std::string const &conf) |
lambda_schedule.size()) { | { |
cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n"); | get_keyval(conf, "forceConstant", force_k, force_k); |
| return COLVARS_OK; |
} | } |
| |
if (target_nstages) { | |
// This means that either numStages of lambdaSchedule has been provided | |
| colvarbias_restraint_moving::colvarbias_restraint_moving(char const *key) |
| { |
| target_nstages = 0; |
| target_nsteps = 0; |
stage = 0; | stage = 0; |
restraint_FE = 0.0; | b_chg_centers = false; |
| b_chg_force_k = false; |
} | } |
| |
if (get_keyval(conf, "targetForceExponent", force_k_exp, 1.0)) { | |
if (! b_chg_force_k) | int colvarbias_restraint_moving::init(std::string const &conf) |
cvm::log("Warning: not changing force constant: targetForceExponent will be ignored\n"); | { |
if (force_k_exp < 1.0) | if (b_chg_centers && b_chg_force_k) { |
cvm::log("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n"); | cvm::error("Error: cannot specify both targetCenters and targetForceConstant.\n", |
} | INPUT_ERROR); |
| return INPUT_ERROR; |
} | } |
| |
get_keyval(conf, "outputCenters", b_output_centers, false); | if (b_chg_centers || b_chg_force_k) { |
if (b_chg_centers) { | |
get_keyval(conf, "outputAccumulatedWork", b_output_acc_work, false); | get_keyval(conf, "targetNumSteps", target_nsteps, target_nsteps); |
} else { | if (!target_nsteps) { |
b_output_acc_work = false; | cvm::error("Error: targetNumSteps must be non-zero.\n", INPUT_ERROR); |
| return cvm::get_error(); |
} | } |
acc_work = 0.0; | |
| |
if (cvm::debug()) | if (get_keyval(conf, "targetNumStages", target_nstages, target_nstages) && |
cvm::log("Done initializing a new restraint bias.\n"); | lambda_schedule.size()) { |
| cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", INPUT_ERROR); |
| return cvm::get_error(); |
| } |
| } |
| |
return COLVARS_OK; | return COLVARS_OK; |
} | } |
| |
| |
colvarbias_restraint::~colvarbias_restraint() | std::string const colvarbias_restraint_moving::get_state_params() const |
{ | { |
if (cvm::n_rest_biases > 0) | std::ostringstream os; |
cvm::n_rest_biases -= 1; | os.setf(std::ios::scientific, std::ios::floatfield); |
| if (b_chg_centers || b_chg_force_k) { |
| // TODO move this |
| if (target_nstages) { |
| os << "stage " << std::setw(cvm::it_width) |
| << stage << "\n"; |
| } |
| } |
| return os.str(); |
} | } |
| |
| |
void colvarbias_restraint::change_configuration(std::string const &conf) | int colvarbias_restraint_moving::set_state_params(std::string const &conf) |
{ | { |
get_keyval(conf, "forceConstant", force_k, force_k); | if (b_chg_centers || b_chg_force_k) { |
if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { | if (target_nstages) { |
for (size_t i = 0; i < colvars.size(); i++) { | // cvm::log ("Reading current stage from the restart.\n"); |
colvar_centers[i].type(colvars[i]->value()); | if (!get_keyval(conf, "stage", stage)) |
colvar_centers[i].apply_constraints(); | cvm::error("Error: current stage is missing from the restart.\n"); |
colvar_centers_raw[i].type(colvars[i]->value()); | |
colvar_centers_raw[i] = colvar_centers[i]; | |
} | } |
} | } |
| return COLVARS_OK; |
} | } |
| |
| |
cvm::real colvarbias_restraint::energy_difference(std::string const &conf) | |
| colvarbias_restraint_centers_moving::colvarbias_restraint_centers_moving(char const *key) |
| : colvarbias(key), |
| colvarbias_restraint(key), |
| colvarbias_restraint_centers(key), |
| colvarbias_restraint_moving(key) |
{ | { |
std::vector<colvarvalue> alt_colvar_centers; | b_chg_centers = false; |
cvm::real alt_force_k; | b_output_centers = false; |
cvm::real alt_bias_energy = 0.0; | b_output_acc_work = false; |
| acc_work = 0.0; |
| } |
| |
get_keyval(conf, "forceConstant", alt_force_k, force_k); | |
| |
alt_colvar_centers.resize(colvars.size()); | int colvarbias_restraint_centers_moving::init(std::string const &conf) |
| { |
| colvarbias_restraint_centers::init(conf); |
| |
| if (cvm::debug()) { |
| cvm::log("colvarbias_restraint: parsing target centers.\n"); |
| } |
| |
size_t i; | size_t i; |
for (i = 0; i < colvars.size(); i++) { | if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) { |
alt_colvar_centers[i].type(colvars[i]->value()); | if (colvar_centers.size() != num_variables()) { |
| cvm::error("Error: number of target centers does not match " |
| "that of collective variables.\n"); |
} | } |
if (get_keyval(conf, "centers", alt_colvar_centers, colvar_centers)) { | b_chg_centers = true; |
for (i = 0; i < colvars.size(); i++) { | for (i = 0; i < target_centers.size(); i++) { |
alt_colvar_centers[i].apply_constraints(); | target_centers[i].apply_constraints(); |
} | } |
} | } |
| |
for (i = 0; i < colvars.size(); i++) { | if (b_chg_centers) { |
alt_bias_energy += restraint_potential(restraint_convert_k(alt_force_k, colvars[i]->width), | // parse moving restraint options |
colvars[i], | colvarbias_restraint_moving::init(conf); |
alt_colvar_centers[i]); | } else { |
| target_centers.clear(); |
| return COLVARS_OK; |
} | } |
| |
return alt_bias_energy - bias_energy; | get_keyval(conf, "outputCenters", b_output_centers, b_output_centers); |
| get_keyval(conf, "outputAccumulatedWork", b_output_acc_work, b_output_acc_work); |
| |
| return COLVARS_OK; |
} | } |
| |
| |
int colvarbias_restraint::update() | int colvarbias_restraint_centers_moving::update() |
{ | { |
bias_energy = 0.0; | if (b_chg_centers) { |
| |
if (cvm::debug()) | |
cvm::log("Updating the restraint bias \""+this->name+"\".\n"); | |
| |
// Setup first stage of staged variable force constant calculation | if (cvm::debug()) { |
if (b_chg_force_k && target_nstages && cvm::step_absolute() == 0) { | cvm::log("Updating centers for the restraint bias \""+ |
cvm::real lambda; | this->name+"\": "+cvm::to_str(colvar_centers)+".\n"); |
if (lambda_schedule.size()) { | |
lambda = lambda_schedule[0]; | |
} else { | |
lambda = 0.0; | |
} | |
force_k = starting_force_k + (target_force_k - starting_force_k) | |
* std::pow(lambda, force_k_exp); | |
cvm::log("Restraint " + this->name + ", stage " + | |
cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); | |
cvm::log("Setting force constant to " + cvm::to_str(force_k)); | |
} | } |
| |
if (b_chg_centers) { | |
| |
if (!centers_incr.size()) { | if (!centers_incr.size()) { |
// if this is the first calculation, calculate the advancement | // if this is the first calculation, calculate the advancement |
// at each simulation step (or stage, if applicable) | // at each simulation step (or stage, if applicable) |
// (take current stage into account: it can be non-zero | // (take current stage into account: it can be non-zero |
// if we are restarting a staged calculation) | // if we are restarting a staged calculation) |
centers_incr.resize(colvars.size()); | centers_incr.resize(num_variables()); |
for (size_t i = 0; i < colvars.size(); i++) { | for (size_t i = 0; i < num_variables(); i++) { |
centers_incr[i].type(colvars[i]->value()); | centers_incr[i].type(variables(i)->value()); |
centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) / | centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) / |
cvm::real( target_nstages ? (target_nstages - stage) : | cvm::real( target_nstages ? (target_nstages - stage) : |
(target_nsteps - cvm::step_absolute())); | (target_nsteps - cvm::step_absolute())); |
| |
&& (cvm::step_absolute() % target_nsteps) == 0 | && (cvm::step_absolute() % target_nsteps) == 0 |
&& stage < target_nstages) { | && stage < target_nstages) { |
| |
for (size_t i = 0; i < colvars.size(); i++) { | for (size_t i = 0; i < num_variables(); i++) { |
colvar_centers_raw[i] += centers_incr[i]; | colvar_centers_raw[i] += centers_incr[i]; |
colvar_centers[i] = colvar_centers_raw[i]; | colvar_centers[i] = colvar_centers_raw[i]; |
colvars[i]->wrap(colvar_centers[i]); | variables(i)->wrap(colvar_centers[i]); |
colvar_centers[i].apply_constraints(); | colvar_centers[i].apply_constraints(); |
} | } |
stage++; | stage++; |
| |
} else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) { | } else if ((cvm::step_relative() > 0) && (cvm::step_absolute() <= target_nsteps)) { |
// move the restraint centers in the direction of the targets | // move the restraint centers in the direction of the targets |
// (slow growth) | // (slow growth) |
for (size_t i = 0; i < colvars.size(); i++) { | for (size_t i = 0; i < num_variables(); i++) { |
colvar_centers_raw[i] += centers_incr[i]; | colvar_centers_raw[i] += centers_incr[i]; |
colvar_centers[i] = colvar_centers_raw[i]; | colvar_centers[i] = colvar_centers_raw[i]; |
colvars[i]->wrap(colvar_centers[i]); | variables(i)->wrap(colvar_centers[i]); |
colvar_centers[i].apply_constraints(); | colvar_centers[i].apply_constraints(); |
} | } |
} | } |
| |
if (cvm::debug()) | if (cvm::debug()) { |
cvm::log("Current centers for the restraint bias \""+ | cvm::log("New centers for the restraint bias \""+ |
this->name+"\": "+cvm::to_str(colvar_centers)+".\n"); | this->name+"\": "+cvm::to_str(colvar_centers)+".\n"); |
} | } |
| } |
| |
| return COLVARS_OK; |
| } |
| |
| |
| int colvarbias_restraint_centers_moving::update_acc_work() |
| { |
| if (b_output_acc_work) { |
| if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) { |
| for (size_t i = 0; i < num_variables(); i++) { |
| // project forces on the calculated increments at this step |
| acc_work += colvar_forces[i] * centers_incr[i]; |
| } |
| } |
| } |
| return COLVARS_OK; |
| } |
| |
| |
| std::string const colvarbias_restraint_centers_moving::get_state_params() const |
| { |
| std::ostringstream os; |
| os.setf(std::ios::scientific, std::ios::floatfield); |
| |
| if (b_chg_centers) { |
| size_t i; |
| os << "centers "; |
| for (i = 0; i < num_variables(); i++) { |
| os << " " |
| << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) |
| << colvar_centers[i]; |
| } |
| os << "\n"; |
| os << "centers_raw "; |
| for (i = 0; i < num_variables(); i++) { |
| os << " " |
| << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) |
| << colvar_centers_raw[i]; |
| } |
| os << "\n"; |
| |
| if (b_output_acc_work) { |
| os << "accumulatedWork " |
| << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width) |
| << acc_work << "\n"; |
| } |
| } |
| |
| return colvarbias_restraint_moving::get_state_params() + os.str(); |
| } |
| |
| |
| int colvarbias_restraint_centers_moving::set_state_params(std::string const &conf) |
| { |
| colvarbias_restraint::set_state_params(conf); |
| |
| if (b_chg_centers) { |
| // cvm::log ("Reading the updated restraint centers from the restart.\n"); |
| if (!get_keyval(conf, "centers", colvar_centers)) |
| cvm::error("Error: restraint centers are missing from the restart.\n"); |
| if (!get_keyval(conf, "centers_raw", colvar_centers_raw)) |
| cvm::error("Error: \"raw\" restraint centers are missing from the restart.\n"); |
| if (b_output_acc_work) { |
| if (!get_keyval(conf, "accumulatedWork", acc_work)) |
| cvm::error("Error: accumulatedWork is missing from the restart.\n"); |
| } |
| } |
| |
| return COLVARS_OK; |
| } |
| |
| |
| std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostream &os) |
| { |
| if (b_output_centers) { |
| for (size_t i = 0; i < num_variables(); i++) { |
| size_t const this_cv_width = (variables(i)->value()).output_width(cvm::cv_width); |
| os << " x0_" |
| << cvm::wrap_string(variables(i)->name, this_cv_width-3); |
| } |
| } |
| |
| if (b_output_acc_work) { |
| os << " W_" |
| << cvm::wrap_string(this->name, cvm::en_width-2); |
| } |
| |
| return os; |
| } |
| |
| |
| std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os) |
| { |
| if (b_output_centers) { |
| for (size_t i = 0; i < num_variables(); i++) { |
| os << " " |
| << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) |
| << colvar_centers[i]; |
| } |
| } |
| |
| if (b_output_acc_work) { |
| os << " " |
| << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width) |
| << acc_work; |
| } |
| |
| return os; |
| } |
| |
| |
| |
| colvarbias_restraint_k_moving::colvarbias_restraint_k_moving(char const *key) |
| : colvarbias(key), |
| colvarbias_restraint(key), |
| colvarbias_restraint_k(key), |
| colvarbias_restraint_moving(key) |
| { |
| b_chg_force_k = false; |
| target_equil_steps = 0; |
| target_force_k = 0.0; |
| starting_force_k = 0.0; |
| force_k_exp = 1.0; |
| restraint_FE = 0.0; |
| } |
| |
| |
| int colvarbias_restraint_k_moving::init(std::string const &conf) |
| { |
| colvarbias_restraint_k::init(conf); |
| |
| if (get_keyval(conf, "targetForceConstant", target_force_k, target_force_k)) { |
| starting_force_k = force_k; |
| b_chg_force_k = true; |
| } |
| |
if (b_chg_force_k) { | if (b_chg_force_k) { |
// Coupling parameter, between 0 and 1 | // parse moving restraint options |
| colvarbias_restraint_moving::init(conf); |
| } else { |
| return COLVARS_OK; |
| } |
| |
| get_keyval(conf, "targetEquilSteps", target_equil_steps, target_equil_steps); |
| |
| if (get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule) && |
| target_nstages > 0) { |
| cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", INPUT_ERROR); |
| return cvm::get_error(); |
| } |
| |
| if (lambda_schedule.size()) { |
| // There is one more lambda-point than stages |
| target_nstages = lambda_schedule.size() - 1; |
| } |
| |
| if (get_keyval(conf, "targetForceExponent", force_k_exp, force_k_exp)) { |
| if (! b_chg_force_k) |
| cvm::log("Warning: not changing force constant: targetForceExponent will be ignored\n"); |
| } |
| if (force_k_exp < 1.0) { |
| cvm::log("Warning: for all practical purposes, targetForceExponent should be 1.0 or greater.\n"); |
| } |
| |
| return COLVARS_OK; |
| } |
| |
| |
| int colvarbias_restraint_k_moving::update() |
| { |
| if (b_chg_force_k) { |
| |
cvm::real lambda; | cvm::real lambda; |
| |
if (target_nstages) { | if (target_nstages) { |
| |
| if (cvm::step_absolute() == 0) { |
| // Setup first stage of staged variable force constant calculation |
| if (lambda_schedule.size()) { |
| lambda = lambda_schedule[0]; |
| } else { |
| lambda = 0.0; |
| } |
| force_k = starting_force_k + (target_force_k - starting_force_k) |
| * std::pow(lambda, force_k_exp); |
| cvm::log("Restraint " + this->name + ", stage " + |
| cvm::to_str(stage) + " : lambda = " + cvm::to_str(lambda)); |
| cvm::log("Setting force constant to " + cvm::to_str(force_k)); |
| } |
| |
// TI calculation: estimate free energy derivative | // TI calculation: estimate free energy derivative |
// need current lambda | // need current lambda |
if (lambda_schedule.size()) { | if (lambda_schedule.size()) { |
| |
| |
// Square distance normalized by square colvar width | // Square distance normalized by square colvar width |
cvm::real dist_sq = 0.0; | cvm::real dist_sq = 0.0; |
for (size_t i = 0; i < colvars.size(); i++) { | for (size_t i = 0; i < num_variables(); i++) { |
dist_sq += colvars[i]->dist2(colvars[i]->value(), colvar_centers[i]) | dist_sq += d_restraint_potential_dk(i); |
/ (colvars[i]->width * colvars[i]->width); | |
} | } |
| |
restraint_FE += 0.5 * force_k_exp * std::pow(lambda, force_k_exp - 1.0) | restraint_FE += 0.5 * force_k_exp * std::pow(lambda, force_k_exp - 1.0) |
| |
cvm::log("Setting force constant to " + cvm::to_str(force_k)); | cvm::log("Setting force constant to " + cvm::to_str(force_k)); |
} | } |
} | } |
| |
} else if (cvm::step_absolute() <= target_nsteps) { | } else if (cvm::step_absolute() <= target_nsteps) { |
| |
// update force constant (slow growth) | // update force constant (slow growth) |
lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps); | lambda = cvm::real(cvm::step_absolute()) / cvm::real(target_nsteps); |
force_k = starting_force_k + (target_force_k - starting_force_k) | force_k = starting_force_k + (target_force_k - starting_force_k) |
| |
} | } |
} | } |
| |
if (cvm::debug()) | return COLVARS_OK; |
cvm::log("Done updating the restraint bias \""+this->name+"\".\n"); | } |
| |
// Force and energy calculation | |
for (size_t i = 0; i < colvars.size(); i++) { | std::string const colvarbias_restraint_k_moving::get_state_params() const |
colvar_forces[i].type(colvars[i]->value()); | { |
colvar_forces[i] = -1.0 * restraint_force(restraint_convert_k(force_k, colvars[i]->width), | std::ostringstream os; |
colvars[i], | os.setf(std::ios::scientific, std::ios::floatfield); |
colvar_centers[i]); | if (b_chg_force_k) { |
bias_energy += restraint_potential(restraint_convert_k(force_k, colvars[i]->width), | os << "forceConstant " |
colvars[i], | << std::setprecision(cvm::en_prec) |
colvar_centers[i]); | << std::setw(cvm::en_width) << force_k << "\n"; |
if (cvm::debug()) { | |
cvm::log("dist_grad["+cvm::to_str(i)+ | |
"] = "+cvm::to_str(colvars[i]->dist2_lgrad(colvars[i]->value(), | |
colvar_centers[i]))+"\n"); | |
} | } |
| return colvarbias_restraint_moving::get_state_params() + os.str(); |
} | } |
| |
if (b_output_acc_work) { | |
if ((cvm::step_relative() > 0) || (cvm::step_absolute() == 0)) { | int colvarbias_restraint_k_moving::set_state_params(std::string const &conf) |
for (size_t i = 0; i < colvars.size(); i++) { | { |
// project forces on the calculated increments at this step | colvarbias_restraint::set_state_params(conf); |
acc_work += colvar_forces[i] * centers_incr[i]; | |
| if (b_chg_force_k) { |
| // cvm::log ("Reading the updated force constant from the restart.\n"); |
| if (!get_keyval(conf, "forceConstant", force_k, force_k)) |
| cvm::error("Error: force constant is missing from the restart.\n"); |
} | } |
| |
| return COLVARS_OK; |
} | } |
| |
| |
| std::ostream & colvarbias_restraint_k_moving::write_traj_label(std::ostream &os) |
| { |
| return os; |
} | } |
| |
if (cvm::debug()) | |
cvm::log("Current forces for the restraint bias \""+ | |
this->name+"\": "+cvm::to_str(colvar_forces)+".\n"); | |
| |
return COLVARS_OK; | std::ostream & colvarbias_restraint_k_moving::write_traj(std::ostream &os) |
| { |
| return os; |
} | } |
| |
| |
std::istream & colvarbias_restraint::read_restart(std::istream &is) | |
| // redefined due to legacy state file keyword "harmonic" |
| std::istream & colvarbias_restraint::read_state(std::istream &is) |
{ | { |
size_t const start_pos = is.tellg(); | size_t const start_pos = is.tellg(); |
| |
cvm::log("Restarting restraint bias \""+ | |
this->name+"\".\n"); | |
| |
std::string key, brace, conf; | std::string key, brace, conf; |
if ( !(is >> key) || !(key == "restraint" || key == "harmonic") || | if ( !(is >> key) || !(key == "restraint" || key == "harmonic") || |
!(is >> brace) || !(brace == "{") || | !(is >> brace) || !(brace == "{") || |
!(is >> colvarparse::read_block("configuration", conf)) ) { | !(is >> colvarparse::read_block("configuration", conf)) || |
| (set_state_params(conf) != COLVARS_OK) ) { |
cvm::log("Error: in reading restart configuration for restraint bias \""+ | cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+ |
this->name+"\" at position "+ | this->name+"\" at position "+ |
cvm::to_str(is.tellg())+" in stream.\n"); | cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR); |
is.clear(); | is.clear(); |
is.seekg(start_pos, std::ios::beg); | is.seekg(start_pos, std::ios::beg); |
is.setstate(std::ios::failbit); | is.setstate(std::ios::failbit); |
return is; | return is; |
} | } |
| |
// int id = -1; | if (!read_state_data(is)) { |
std::string name = ""; | cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+ |
// if ( ( (colvarparse::get_keyval (conf, "id", id, -1, colvarparse::parse_silent)) && | this->name+"\" at position "+ |
// (id != this->id) ) || | cvm::to_str(is.tellg())+" in stream.\n", INPUT_ERROR); |
if ( (colvarparse::get_keyval(conf, "name", name, std::string(""), colvarparse::parse_silent)) && | is.clear(); |
(name != this->name) ) | is.seekg(start_pos, std::ios::beg); |
cvm::error("Error: in the restart file, the " | is.setstate(std::ios::failbit); |
"\"restraint\" block has a wrong name\n"); | |
// if ( (id == -1) && (name == "") ) { | |
if (name.size() == 0) { | |
cvm::error("Error: \"restraint\" block in the restart file " | |
"has no identifiers.\n"); | |
} | |
| |
if (b_chg_centers) { | |
// cvm::log ("Reading the updated restraint centers from the restart.\n"); | |
if (!get_keyval(conf, "centers", colvar_centers)) | |
cvm::error("Error: restraint centers are missing from the restart.\n"); | |
if (!get_keyval(conf, "centers_raw", colvar_centers_raw)) | |
cvm::error("Error: \"raw\" restraint centers are missing from the restart.\n"); | |
} | |
| |
if (b_chg_force_k) { | |
// cvm::log ("Reading the updated force constant from the restart.\n"); | |
if (!get_keyval(conf, "forceConstant", force_k)) | |
cvm::error("Error: force constant is missing from the restart.\n"); | |
} | |
| |
if (target_nstages) { | |
// cvm::log ("Reading current stage from the restart.\n"); | |
if (!get_keyval(conf, "stage", stage)) | |
cvm::error("Error: current stage is missing from the restart.\n"); | |
} | |
| |
if (b_output_acc_work) { | |
if (!get_keyval(conf, "accumulatedWork", acc_work)) | |
cvm::error("Error: accumulatedWork is missing from the restart.\n"); | |
} | } |
| |
is >> brace; | is >> brace; |
if (brace != "}") { | if (brace != "}") { |
cvm::error("Error: corrupt restart information for restraint bias \""+ | cvm::log("brace = "+brace+"\n"); |
| cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+ |
this->name+"\": no matching brace at position "+ | this->name+"\": no matching brace at position "+ |
cvm::to_str(is.tellg())+" in the restart file.\n"); | cvm::to_str(is.tellg())+" in stream.\n"); |
is.setstate(std::ios::failbit); | is.setstate(std::ios::failbit); |
} | } |
| |
return is; | return is; |
} | } |
| |
| |
std::ostream & colvarbias_restraint::write_restart(std::ostream &os) | std::ostream & colvarbias_restraint::write_state(std::ostream &os) |
{ | { |
| os.setf(std::ios::scientific, std::ios::floatfield); |
os << "restraint {\n" | os << "restraint {\n" |
<< " configuration {\n" | << " configuration {\n"; |
// << " id " << this->id << "\n" | std::istringstream is(get_state_params()); |
<< " name " << this->name << "\n"; | std::string line; |
| while (std::getline(is, line)) { |
| os << " " << line << "\n"; |
| } |
| os << " }\n"; |
| write_state_data(os); |
| os << "}\n\n"; |
| return os; |
| } |
| |
if (b_chg_centers) { | |
size_t i; | |
os << " centers "; | colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key) |
for (i = 0; i < colvars.size(); i++) { | : colvarbias(key), |
os << " " << colvar_centers[i]; | colvarbias_restraint(key), |
| colvarbias_restraint_centers(key), |
| colvarbias_restraint_moving(key), |
| colvarbias_restraint_k(key), |
| colvarbias_restraint_centers_moving(key), |
| colvarbias_restraint_k_moving(key) |
| { |
} | } |
os << "\n"; | |
os << " centers_raw "; | |
for (i = 0; i < colvars.size(); i++) { | int colvarbias_restraint_harmonic::init(std::string const &conf) |
os << " " << colvar_centers_raw[i]; | { |
| colvarbias_restraint::init(conf); |
| colvarbias_restraint_moving::init(conf); |
| colvarbias_restraint_centers_moving::init(conf); |
| colvarbias_restraint_k_moving::init(conf); |
| |
| for (size_t i = 0; i < num_variables(); i++) { |
| if (variables(i)->width != 1.0) |
| cvm::log("The force constant for colvar \""+variables(i)->name+ |
| "\" will be rescaled to "+ |
| cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+ |
| " according to the specified width.\n"); |
} | } |
os << "\n"; | |
| return COLVARS_OK; |
} | } |
| |
if (b_chg_force_k) { | |
os << " forceConstant " | int colvarbias_restraint_harmonic::update() |
<< std::setprecision(cvm::en_prec) | { |
<< std::setw(cvm::en_width) << force_k << "\n"; | // update parameters (centers or force constant) |
| colvarbias_restraint_centers_moving::update(); |
| colvarbias_restraint_k_moving::update(); |
| |
| // update restraint energy and forces |
| colvarbias_restraint::update(); |
| |
| // update accumulated work using the current forces |
| colvarbias_restraint_centers_moving::update_acc_work(); |
| |
| return COLVARS_OK; |
} | } |
| |
if (target_nstages) { | |
os << " stage " << std::setw(cvm::it_width) | cvm::real colvarbias_restraint_harmonic::restraint_potential(size_t i) const |
<< stage << "\n"; | { |
| return 0.5 * force_k / (variables(i)->width * variables(i)->width) * |
| variables(i)->dist2(variables(i)->value(), colvar_centers[i]); |
} | } |
| |
if (b_output_acc_work) { | |
os << " accumulatedWork " << acc_work << "\n"; | colvarvalue const colvarbias_restraint_harmonic::restraint_force(size_t i) const |
| { |
| return -0.5 * force_k / (variables(i)->width * variables(i)->width) * |
| variables(i)->dist2_lgrad(variables(i)->value(), colvar_centers[i]); |
} | } |
| |
os << " }\n" | |
<< "}\n\n"; | |
| |
return os; | cvm::real colvarbias_restraint_harmonic::d_restraint_potential_dk(size_t i) const |
| { |
| return 0.5 / (variables(i)->width * variables(i)->width) * |
| variables(i)->dist2(variables(i)->value(), colvar_centers[i]); |
} | } |
| |
| |
std::ostream & colvarbias_restraint::write_traj_label(std::ostream &os) | std::string const colvarbias_restraint_harmonic::get_state_params() const |
{ | { |
os << " "; | return colvarbias_restraint::get_state_params() + |
| colvarbias_restraint_centers_moving::get_state_params() + |
| colvarbias_restraint_k_moving::get_state_params(); |
| } |
| |
if (b_output_energy) | |
os << " E_" | |
<< cvm::wrap_string(this->name, cvm::en_width-2); | |
| |
if (b_output_centers) | int colvarbias_restraint_harmonic::set_state_params(std::string const &conf) |
for (size_t i = 0; i < colvars.size(); i++) { | { |
size_t const this_cv_width = (colvars[i]->value()).output_width(cvm::cv_width); | int error_code = COLVARS_OK; |
os << " x0_" | error_code |= colvarbias_restraint::set_state_params(conf); |
<< cvm::wrap_string(colvars[i]->name, this_cv_width-3); | error_code |= colvarbias_restraint_centers_moving::set_state_params(conf); |
| error_code |= colvarbias_restraint_k_moving::set_state_params(conf); |
| return error_code; |
} | } |
| |
if (b_output_acc_work) | |
os << " W_" | |
<< cvm::wrap_string(this->name, cvm::en_width-2); | |
| |
| std::ostream & colvarbias_restraint_harmonic::write_traj_label(std::ostream &os) |
| { |
| colvarbias_restraint::write_traj_label(os); |
| colvarbias_restraint_centers_moving::write_traj_label(os); |
| colvarbias_restraint_k_moving::write_traj_label(os); |
return os; | return os; |
} | } |
| |
| |
std::ostream & colvarbias_restraint::write_traj(std::ostream &os) | std::ostream & colvarbias_restraint_harmonic::write_traj(std::ostream &os) |
{ | { |
os << " "; | colvarbias_restraint::write_traj(os); |
| colvarbias_restraint_centers_moving::write_traj(os); |
| colvarbias_restraint_k_moving::write_traj(os); |
| return os; |
| } |
| |
if (b_output_energy) | |
os << " " | |
<< std::setprecision(cvm::en_prec) << std::setw(cvm::en_width) | |
<< bias_energy; | |
| |
if (b_output_centers) | int colvarbias_restraint_harmonic::change_configuration(std::string const &conf) |
for (size_t i = 0; i < colvars.size(); i++) { | { |
os << " " | return colvarbias_restraint_centers::change_configuration(conf) | |
<< std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) | colvarbias_restraint_k::change_configuration(conf); |
<< colvar_centers[i]; | |
} | } |
| |
if (b_output_acc_work) | |
os << " " | |
<< std::setprecision(cvm::en_prec) << std::setw(cvm::en_width) | |
<< acc_work; | |
| |
return os; | cvm::real colvarbias_restraint_harmonic::energy_difference(std::string const &conf) |
| { |
| cvm::real const old_bias_energy = bias_energy; |
| cvm::real const old_force_k = force_k; |
| std::vector<colvarvalue> const old_centers = colvar_centers; |
| |
| change_configuration(conf); |
| update(); |
| |
| cvm::real const result = (bias_energy - old_bias_energy); |
| |
| bias_energy = old_bias_energy; |
| force_k = old_force_k; |
| colvar_centers = old_centers; |
| |
| return result; |
} | } |
| |
| |
colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key) | |
: colvarbias_restraint(key) | colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char const *key) |
| : colvarbias(key), |
| colvarbias_restraint(key), |
| colvarbias_restraint_k(key), |
| colvarbias_restraint_moving(key), |
| colvarbias_restraint_k_moving(key) |
{ | { |
} | } |
| |
| |
int colvarbias_restraint_harmonic::init(std::string const &conf) | int colvarbias_restraint_harmonic_walls::init(std::string const &conf) |
{ | { |
colvarbias_restraint::init(conf); | colvarbias_restraint::init(conf); |
| colvarbias_restraint_moving::init(conf); |
| colvarbias_restraint_k_moving::init(conf); |
| |
| provide(f_cvb_scalar_variables); |
| enable(f_cvb_scalar_variables); |
| |
| size_t i; |
| |
| bool b_null_lower_walls = false; |
| if (lower_walls.size() == 0) { |
| b_null_lower_walls = true; |
| lower_walls.resize(num_variables()); |
| for (i = 0; i < num_variables(); i++) { |
| lower_walls[i].type(variables(i)->value()); |
| lower_walls[i].reset(); |
| } |
| } |
| if (!get_keyval(conf, "lowerWalls", lower_walls, lower_walls) && |
| b_null_lower_walls) { |
| cvm::log("Lower walls were not provided.\n"); |
| lower_walls.resize(0); |
| } |
| |
| bool b_null_upper_walls = false; |
| if (upper_walls.size() == 0) { |
| b_null_upper_walls = true; |
| upper_walls.resize(num_variables()); |
| for (i = 0; i < num_variables(); i++) { |
| upper_walls[i].type(variables(i)->value()); |
| upper_walls[i].reset(); |
| } |
| } |
| if (!get_keyval(conf, "upperWalls", upper_walls, upper_walls) && |
| b_null_upper_walls) { |
| cvm::log("Upper walls were not provided.\n"); |
| upper_walls.resize(0); |
| } |
| |
| if ((lower_walls.size() == 0) && (upper_walls.size() == 0)) { |
| cvm::error("Error: no walls provided.\n", INPUT_ERROR); |
| return INPUT_ERROR; |
| } |
| |
| if ((lower_walls.size() == 0) || (upper_walls.size() == 0)) { |
| for (i = 0; i < num_variables(); i++) { |
| if (variables(i)->is_enabled(f_cv_periodic)) { |
| cvm::error("Error: at least one variable is periodic, " |
| "both walls must be provided .\n", INPUT_ERROR); |
| return INPUT_ERROR; |
| } |
| } |
| } |
| |
| if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) { |
| for (i = 0; i < num_variables(); i++) { |
| if (lower_walls[i] >= upper_walls[i]) { |
| cvm::error("Error: one upper wall, "+ |
| cvm::to_str(upper_walls[i])+ |
| ", is not higher than the lower wall, "+ |
| cvm::to_str(lower_walls[i])+".\n", |
| INPUT_ERROR); |
| } |
| } |
| } |
| |
for (size_t i = 0; i < colvars.size(); i++) { | for (i = 0; i < num_variables(); i++) { |
if (colvars[i]->width != 1.0) | if (variables(i)->width != 1.0) |
cvm::log("The force constant for colvar \""+colvars[i]->name+ | cvm::log("The force constant for colvar \""+variables(i)->name+ |
"\" will be rescaled to "+ | "\" will be rescaled to "+ |
cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+ | cvm::to_str(force_k / (variables(i)->width * variables(i)->width))+ |
" according to the specified width.\n"); | " according to the specified width.\n"); |
} | } |
| |
| return COLVARS_OK; |
| } |
| |
| |
| int colvarbias_restraint_harmonic_walls::update() |
| { |
| colvarbias_restraint_k_moving::update(); |
| |
| colvarbias_restraint::update(); |
| |
return COLVARS_OK; | return COLVARS_OK; |
} | } |
| |
| |
cvm::real colvarbias_restraint_harmonic::restraint_potential(cvm::real k, | void colvarbias_restraint_harmonic_walls::communicate_forces() |
colvar const *x, | { |
colvarvalue const &xcenter) const | for (size_t i = 0; i < num_variables(); i++) { |
| if (cvm::debug()) { |
| cvm::log("Communicating a force to colvar \""+ |
| variables(i)->name+"\".\n"); |
| } |
| variables(i)->add_bias_force_actual_value(colvar_forces[i]); |
| } |
| } |
| |
| |
| cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const |
| { |
| colvar *cv = variables(i); |
| colvarvalue const &cvv = variables(i)->actual_value(); |
| |
| // For a periodic colvar, both walls may be applicable at the same time |
| // in which case we pick the closer one |
| |
| if (cv->is_enabled(f_cv_periodic)) { |
| cvm::real const lower_wall_dist2 = cv->dist2(cvv, lower_walls[i]); |
| cvm::real const upper_wall_dist2 = cv->dist2(cvv, upper_walls[i]); |
| if (lower_wall_dist2 < upper_wall_dist2) { |
| cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]); |
| if (grad < 0.0) { return 0.5 * grad; } |
| } else { |
| cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]); |
| if (grad > 0.0) { return 0.5 * grad; } |
| } |
| return 0.0; |
| } |
| |
| if (lower_walls.size() > 0) { |
| cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]); |
| if (grad < 0.0) { return 0.5 * grad; } |
| } |
| if (upper_walls.size() > 0) { |
| cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]); |
| if (grad > 0.0) { return 0.5 * grad; } |
| } |
| return 0.0; |
| } |
| |
| |
| cvm::real colvarbias_restraint_harmonic_walls::restraint_potential(size_t i) const |
| { |
| cvm::real const dist = colvar_distance(i); |
| return 0.5 * force_k / (variables(i)->width * variables(i)->width) * |
| dist * dist; |
| } |
| |
| |
| colvarvalue const colvarbias_restraint_harmonic_walls::restraint_force(size_t i) const |
| { |
| cvm::real const dist = colvar_distance(i); |
| return - force_k / (variables(i)->width * variables(i)->width) * dist; |
| } |
| |
| |
| cvm::real colvarbias_restraint_harmonic_walls::d_restraint_potential_dk(size_t i) const |
| { |
| cvm::real const dist = colvar_distance(i); |
| return 0.5 / (variables(i)->width * variables(i)->width) * |
| dist * dist; |
| } |
| |
| |
| std::string const colvarbias_restraint_harmonic_walls::get_state_params() const |
| { |
| return colvarbias_restraint::get_state_params() + |
| colvarbias_restraint_k_moving::get_state_params(); |
| } |
| |
| |
| int colvarbias_restraint_harmonic_walls::set_state_params(std::string const &conf) |
{ | { |
return 0.5 * k * x->dist2(x->value(), xcenter); | int error_code = COLVARS_OK; |
| error_code |= colvarbias_restraint::set_state_params(conf); |
| error_code |= colvarbias_restraint_k_moving::set_state_params(conf); |
| return error_code; |
} | } |
| |
| |
colvarvalue colvarbias_restraint_harmonic::restraint_force(cvm::real k, | std::ostream & colvarbias_restraint_harmonic_walls::write_traj_label(std::ostream &os) |
colvar const *x, | |
colvarvalue const &xcenter) const | |
{ | { |
return 0.5 * k * x->dist2_lgrad(x->value(), xcenter); | colvarbias_restraint::write_traj_label(os); |
| colvarbias_restraint_k_moving::write_traj_label(os); |
| return os; |
} | } |
| |
| |
cvm::real colvarbias_restraint_harmonic::restraint_convert_k(cvm::real k, | std::ostream & colvarbias_restraint_harmonic_walls::write_traj(std::ostream &os) |
cvm::real dist_measure) const | |
{ | { |
return k / (dist_measure * dist_measure); | colvarbias_restraint::write_traj(os); |
| colvarbias_restraint_k_moving::write_traj(os); |
| return os; |
} | } |
| |
| |
| |
colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key) | colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key) |
: colvarbias_restraint(key) | : colvarbias(key), |
| colvarbias_restraint(key), |
| colvarbias_restraint_centers(key), |
| colvarbias_restraint_moving(key), |
| colvarbias_restraint_k(key), |
| colvarbias_restraint_centers_moving(key), |
| colvarbias_restraint_k_moving(key) |
{ | { |
} | } |
| |
| |
int colvarbias_restraint_linear::init(std::string const &conf) | int colvarbias_restraint_linear::init(std::string const &conf) |
{ | { |
colvarbias_restraint::init(conf); | colvarbias_restraint::init(conf); |
| colvarbias_restraint_moving::init(conf); |
for (size_t i = 0; i < colvars.size(); i++) { | colvarbias_restraint_centers_moving::init(conf); |
if (colvars[i]->width != 1.0) | colvarbias_restraint_k_moving::init(conf); |
cvm::log("The force constant for colvar \""+colvars[i]->name+ | |
| for (size_t i = 0; i < num_variables(); i++) { |
| if (variables(i)->is_enabled(f_cv_periodic)) { |
| cvm::error("Error: linear biases cannot be applied to periodic variables.\n", |
| INPUT_ERROR); |
| return INPUT_ERROR; |
| } |
| if (variables(i)->width != 1.0) |
| cvm::log("The force constant for colvar \""+variables(i)->name+ |
"\" will be rescaled to "+ | "\" will be rescaled to "+ |
cvm::to_str(restraint_convert_k(force_k, colvars[i]->width))+ | cvm::to_str(force_k / variables(i)->width)+ |
" according to the specified width.\n"); | " according to the specified width.\n"); |
} | } |
| |
| return COLVARS_OK; |
| } |
| |
| |
| int colvarbias_restraint_linear::update() |
| { |
| // update parameters (centers or force constant) |
| colvarbias_restraint_centers_moving::update(); |
| colvarbias_restraint_k_moving::update(); |
| |
| // update restraint energy and forces |
| colvarbias_restraint::update(); |
| |
| // update accumulated work using the current forces |
| colvarbias_restraint_centers_moving::update_acc_work(); |
| |
return COLVARS_OK; | return COLVARS_OK; |
} | } |
| |
| |
cvm::real colvarbias_restraint_linear::restraint_potential(cvm::real k, | int colvarbias_restraint_linear::change_configuration(std::string const &conf) |
colvar const *x, | { |
colvarvalue const &xcenter) const | // Only makes sense to change the force constant |
| return colvarbias_restraint_k::change_configuration(conf); |
| } |
| |
| |
| cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf) |
| { |
| cvm::real const old_bias_energy = bias_energy; |
| cvm::real const old_force_k = force_k; |
| |
| change_configuration(conf); |
| update(); |
| |
| cvm::real const result = (bias_energy - old_bias_energy); |
| |
| bias_energy = old_bias_energy; |
| force_k = old_force_k; |
| |
| return result; |
| } |
| |
| |
| cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const |
| { |
| return force_k / variables(i)->width * (variables(i)->value() - colvar_centers[i]); |
| } |
| |
| |
| colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const |
| { |
| return -1.0 * force_k / variables(i)->width; |
| } |
| |
| |
| cvm::real colvarbias_restraint_linear::d_restraint_potential_dk(size_t i) const |
| { |
| return 1.0 / variables(i)->width * (variables(i)->value() - colvar_centers[i]); |
| } |
| |
| |
| std::string const colvarbias_restraint_linear::get_state_params() const |
| { |
| return colvarbias_restraint::get_state_params() + |
| colvarbias_restraint_centers_moving::get_state_params() + |
| colvarbias_restraint_k_moving::get_state_params(); |
| } |
| |
| |
| int colvarbias_restraint_linear::set_state_params(std::string const &conf) |
{ | { |
return k * (x->value() - xcenter); | int error_code = COLVARS_OK; |
| error_code |= colvarbias_restraint::set_state_params(conf); |
| error_code |= colvarbias_restraint_centers_moving::set_state_params(conf); |
| error_code |= colvarbias_restraint_k_moving::set_state_params(conf); |
| return error_code; |
} | } |
| |
| |
colvarvalue colvarbias_restraint_linear::restraint_force(cvm::real k, | std::ostream & colvarbias_restraint_linear::write_traj_label(std::ostream &os) |
colvar const *x, | |
colvarvalue const &xcenter) const | |
{ | { |
return k; | colvarbias_restraint::write_traj_label(os); |
| colvarbias_restraint_centers_moving::write_traj_label(os); |
| colvarbias_restraint_k_moving::write_traj_label(os); |
| return os; |
} | } |
| |
| |
cvm::real colvarbias_restraint_linear::restraint_convert_k(cvm::real k, | std::ostream & colvarbias_restraint_linear::write_traj(std::ostream &os) |
cvm::real dist_measure) const | |
{ | { |
return k / dist_measure; | colvarbias_restraint::write_traj(os); |
| colvarbias_restraint_centers_moving::write_traj(os); |
| colvarbias_restraint_k_moving::write_traj(os); |
| return os; |
} | } |
| |
| |
| |
| |
size_t vector_size = 0; | size_t vector_size = 0; |
size_t icv; | size_t icv; |
for (icv = 0; icv < colvars.size(); icv++) { | for (icv = 0; icv < num_variables(); icv++) { |
vector_size += colvars[icv]->value().size(); | vector_size += variables(icv)->value().size(); |
} | } |
| |
cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size); | cvm::real const norm = 1.0/(std::sqrt(2.0*PI)*gaussian_width*vector_size); |
| |
// calculate the histogram | // calculate the histogram |
p.reset(); | p.reset(); |
for (icv = 0; icv < colvars.size(); icv++) { | for (icv = 0; icv < num_variables(); icv++) { |
colvarvalue const &cv = colvars[icv]->value(); | colvarvalue const &cv = variables(icv)->value(); |
if (cv.type() == colvarvalue::type_scalar) { | if (cv.type() == colvarvalue::type_scalar) { |
cvm::real const cv_value = cv.real_value; | cvm::real const cv_value = cv.real_value; |
size_t igrid; | size_t igrid; |
| |
} | } |
} | } |
} else { | } else { |
// TODO | cvm::error("Error: unsupported type for variable "+variables(icv)->name+".\n", |
| COLVARS_NOT_IMPLEMENTED); |
| return COLVARS_NOT_IMPLEMENTED; |
} | } |
} | } |
| |
| |
bias_energy = 0.5 * force_k_cv * p_diff * p_diff; | bias_energy = 0.5 * force_k_cv * p_diff * p_diff; |
| |
// calculate the forces | // calculate the forces |
for (icv = 0; icv < colvars.size(); icv++) { | for (icv = 0; icv < num_variables(); icv++) { |
colvarvalue const &cv = colvars[icv]->value(); | colvarvalue const &cv = variables(icv)->value(); |
colvarvalue &cv_force = colvar_forces[icv]; | colvarvalue &cv_force = colvar_forces[icv]; |
cv_force.type(cv); | cv_force.type(cv); |
cv_force.reset(); | cv_force.reset(); |
| |
std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os) | std::ostream & colvarbias_restraint_histogram::write_restart(std::ostream &os) |
{ | { |
if (b_write_histogram) { | if (b_write_histogram) { |
std::string file_name(cvm::output_prefix+"."+this->name+".hist.dat"); | std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat"); |
std::ofstream os(file_name.c_str()); | std::ofstream os(file_name.c_str()); |
os << "# " << cvm::wrap_string(colvars[0]->name, cvm::cv_width) | os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width) |
<< " " << "p(" << cvm::wrap_string(colvars[0]->name, cvm::cv_width-3) | << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3) |
<< ")\n"; | << ")\n"; |
size_t igrid; | size_t igrid; |
for (igrid = 0; igrid < p.size(); igrid++) { | for (igrid = 0; igrid < p.size(); igrid++) { |