colvarbias_restraint.C

00001 // -*- c++ -*-
00002 
00003 // This file is part of the Collective Variables module (Colvars).
00004 // The original version of Colvars and its updates are located at:
00005 // https://github.com/Colvars/colvars
00006 // Please update all Colvars source files before making any changes.
00007 // If you wish to distribute your changes, please submit them to the
00008 // Colvars repository at GitHub.
00009 
00010 #include <fstream>
00011 #include <iostream>
00012 #include <iomanip>
00013 
00014 #include "colvarmodule.h"
00015 #include "colvarproxy.h"
00016 #include "colvarvalue.h"
00017 #include "colvarbias_restraint.h"
00018 
00019 
00020 
00021 colvarbias_restraint::colvarbias_restraint(char const *key)
00022 : colvarbias(key), colvarbias_ti(key)
00023 {
00024 state_keyword = "restraint";
00025 }
00026 
00027 
00028 int colvarbias_restraint::init(std::string const &conf)
00029 {
00030 colvarbias::init(conf);
00031 enable(f_cvb_apply_force);
00032 
00033 colvarbias_ti::init(conf);
00034 
00035 if (cvm::debug())
00036 cvm::log("Initializing a new restraint bias.\n");
00037 
00038 return COLVARS_OK;
00039 }
00040 
00041 
00042 int colvarbias_restraint::update()
00043 {
00044 // Update base class (bias_energy and colvar_forces are zeroed there)
00045 colvarbias::update();
00046 
00047 // Force and energy calculation
00048 for (size_t i = 0; i < num_variables(); i++) {
00049 bias_energy += restraint_potential(i);
00050 colvar_forces[i].type(variables(i)->value());
00051 colvar_forces[i].is_derivative();
00052 colvar_forces[i] = restraint_force(i);
00053 }
00054 
00055 if (cvm::debug())
00056 cvm::log("Done updating the restraint bias \""+this->name+"\".\n");
00057 
00058 if (cvm::debug())
00059 cvm::log("Current forces for the restraint bias \""+
00060 this->name+"\": "+cvm::to_str(colvar_forces)+".\n");
00061 
00062 return COLVARS_OK;
00063 }
00064 
00065 
00066 colvarbias_restraint::~colvarbias_restraint()
00067 {
00068 }
00069 
00070 
00071 std::string const colvarbias_restraint::get_state_params() const
00072 {
00073 return colvarbias::get_state_params();
00074 }
00075 
00076 
00077 int colvarbias_restraint::set_state_params(std::string const &conf)
00078 {
00079 return colvarbias::set_state_params(conf);
00080 }
00081 
00082 
00083 std::ostream & colvarbias_restraint::write_traj_label(std::ostream &os)
00084 {
00085 return colvarbias::write_traj_label(os);
00086 }
00087 
00088 
00089 std::ostream & colvarbias_restraint::write_traj(std::ostream &os)
00090 {
00091 return colvarbias::write_traj(os);
00092 }
00093 
00094 
00095 
00096 colvarbias_restraint_centers::colvarbias_restraint_centers(char const *key)
00097 : colvarbias(key), colvarbias_ti(key), colvarbias_restraint(key)
00098 {
00099 }
00100 
00101 
00102 int colvarbias_restraint_centers::init(std::string const &conf)
00103 {
00104 size_t i;
00105 
00106 bool null_centers = (colvar_centers.size() == 0);
00107 if (null_centers) {
00108 // try to initialize the restraint centers for the first time
00109 colvar_centers.resize(num_variables());
00110 for (i = 0; i < num_variables(); i++) {
00111 colvar_centers[i].type(variables(i)->value());
00112 colvar_centers[i].reset();
00113 }
00114 }
00115 
00116 if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
00117 for (i = 0; i < num_variables(); i++) {
00118 if (cvm::debug()) {
00119 cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n");
00120 }
00121 colvar_centers[i].apply_constraints();
00122 }
00123 null_centers = false;
00124 }
00125 
00126 if (null_centers) {
00127 colvar_centers.clear();
00128 cvm::error("Error: must define the initial centers of the restraints.\n", COLVARS_INPUT_ERROR);
00129 return COLVARS_INPUT_ERROR;
00130 }
00131 
00132 if (colvar_centers.size() != num_variables()) {
00133 cvm::error("Error: number of centers does not match "
00134 "that of collective variables.\n", COLVARS_INPUT_ERROR);
00135 return COLVARS_INPUT_ERROR;
00136 }
00137 
00138 return COLVARS_OK;
00139 }
00140 
00141 
00142 int colvarbias_restraint_centers::change_configuration(std::string const &conf)
00143 {
00144 if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) {
00145 for (size_t i = 0; i < num_variables(); i++) {
00146 colvar_centers[i].type(variables(i)->value());
00147 colvar_centers[i].apply_constraints();
00148 }
00149 }
00150 return COLVARS_OK;
00151 }
00152 
00153 
00154 
00155 colvarbias_restraint_k::colvarbias_restraint_k(char const *key)
00156 : colvarbias(key), colvarbias_ti(key), colvarbias_restraint(key)
00157 {
00158 force_k = -1.0;
00159 check_positive_k = true;
00160 }
00161 
00162 
00163 int colvarbias_restraint_k::init(std::string const &conf)
00164 {
00165 get_keyval(conf, "forceConstant", force_k, (force_k > 0.0 ? force_k : 1.0));
00166 if (check_positive_k && (force_k < 0.0)) {
00167 cvm::error("Error: undefined or invalid force constant.\n", COLVARS_INPUT_ERROR);
00168 return COLVARS_INPUT_ERROR;
00169 }
00170 return COLVARS_OK;
00171 }
00172 
00173 
00174 int colvarbias_restraint_k::change_configuration(std::string const &conf)
00175 {
00176 get_keyval(conf, "forceConstant", force_k, force_k);
00177 return COLVARS_OK;
00178 }
00179 
00180 
00181 
00182 colvarbias_restraint_moving::colvarbias_restraint_moving(char const * /* key */)
00183 {
00184 target_nstages = 0;
00185 target_nsteps = 0L;
00186 stage = 0;
00187 acc_work = 0.0;
00188 b_chg_centers = false;
00189 b_chg_force_k = false;
00190 }
00191 
00192 
00193 int colvarbias_restraint_moving::init(std::string const &conf)
00194 {
00195 if (b_chg_centers && b_chg_force_k) {
00196 cvm::error("Error: cannot specify both targetCenters and targetForceConstant.\n",
00197 COLVARS_INPUT_ERROR);
00198 return COLVARS_INPUT_ERROR;
00199 }
00200 
00201 if (b_chg_centers || b_chg_force_k) {
00202 
00203 first_step = cvm::step_absolute();
00204 
00205 get_keyval(conf, "targetNumSteps", target_nsteps, target_nsteps);
00206 if (!target_nsteps) {
00207 cvm::error("Error: targetNumSteps must be non-zero.\n", COLVARS_INPUT_ERROR);
00208 return cvm::get_error();
00209 }
00210 
00211 if (get_keyval(conf, "targetNumStages", target_nstages, target_nstages) &&
00212 lambda_schedule.size()) {
00213 cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", COLVARS_INPUT_ERROR);
00214 return cvm::get_error();
00215 }
00216 
00217 get_keyval_feature(this, conf, "outputAccumulatedWork",
00218 f_cvb_output_acc_work,
00219 is_enabled(f_cvb_output_acc_work));
00220 if (is_enabled(f_cvb_output_acc_work) && (target_nstages > 0)) {
00221 return cvm::error("Error: outputAccumulatedWork and targetNumStages "
00222 "are incompatible.\n", COLVARS_INPUT_ERROR);
00223 }
00224 }
00225 
00226 return COLVARS_OK;
00227 }
00228 
00229 
00230 std::string const colvarbias_restraint_moving::get_state_params() const
00231 {
00232 std::ostringstream os;
00233 os.setf(std::ios::scientific, std::ios::floatfield);
00234 if (b_chg_centers || b_chg_force_k) {
00235 // TODO move this
00236 if (target_nstages) {
00237 os << "stage " << std::setw(cvm::it_width)
00238 << stage << "\n";
00239 }
00240 }
00241 return os.str();
00242 }
00243 
00244 
00245 int colvarbias_restraint_moving::set_state_params(std::string const &conf)
00246 {
00247 if (b_chg_centers || b_chg_force_k) {
00248 if (target_nstages) {
00249 get_keyval(conf, "stage", stage, stage,
00250 colvarparse::parse_restart | colvarparse::parse_required);
00251 }
00252 }
00253 return COLVARS_OK;
00254 }
00255 
00256 
00257 
00258 colvarbias_restraint_centers_moving::colvarbias_restraint_centers_moving(char const *key)
00259 : colvarbias(key),
00260 colvarbias_ti(key),
00261 colvarbias_restraint(key),
00262 colvarbias_restraint_centers(key),
00263 colvarbias_restraint_moving(key)
00264 {
00265 b_chg_centers = false;
00266 b_output_centers = false;
00267 }
00268 
00269 
00270 int colvarbias_restraint_centers_moving::init(std::string const &conf)
00271 {
00272 colvarbias_restraint_centers::init(conf);
00273 
00274 if (cvm::debug()) {
00275 cvm::log("colvarbias_restraint: parsing target centers.\n");
00276 }
00277 
00278 size_t i;
00279 if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) {
00280 if (target_centers.size() != num_variables()) {
00281 cvm::error("Error: number of target centers does not match "
00282 "that of collective variables.\n", COLVARS_INPUT_ERROR);
00283 }
00284 b_chg_centers = true;
00285 for (i = 0; i < target_centers.size(); i++) {
00286 target_centers[i].apply_constraints();
00287 centers_incr.push_back(colvar_centers[i]);
00288 centers_incr[i].reset();
00289 }
00290 }
00291 
00292 if (b_chg_centers) {
00293 // parse moving schedule options
00294 colvarbias_restraint_moving::init(conf);
00295 if (initial_centers.size() == 0) {
00296 // One-time init
00297 initial_centers = colvar_centers;
00298 }
00299 // Call to check that the definition is correct
00300 for (i = 0; i < num_variables(); i++) {
00301 colvarvalue const midpoint =
00302 colvarvalue::interpolate(initial_centers[i],
00303 target_centers[i],
00304 0.5);
00305 }
00306 
00307 } else {
00308 target_centers.clear();
00309 }
00310 
00311 // Output restraint centers even when they do not change; some NAMD REUS
00312 // scripts expect this behavior
00313 get_keyval(conf, "outputCenters", b_output_centers, b_output_centers);
00314 
00315 return COLVARS_OK;
00316 }
00317 
00318 
00319 int colvarbias_restraint_centers_moving::update_centers(cvm::real lambda)
00320 {
00321 if (cvm::debug()) {
00322 cvm::log("Updating centers for the restraint bias \""+
00323 this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
00324 }
00325 size_t i;
00326 for (i = 0; i < num_variables(); i++) {
00327 colvarvalue const c_new = colvarvalue::interpolate(initial_centers[i],
00328 target_centers[i],
00329 lambda);
00330 centers_incr[i] = 0.5 * c_new.dist2_grad(colvar_centers[i]);
00331 colvar_centers[i] = c_new;
00332 variables(i)->wrap(colvar_centers[i]);
00333 }
00334 if (cvm::debug()) {
00335 cvm::log("New centers for the restraint bias \""+
00336 this->name+"\": "+cvm::to_str(colvar_centers)+".\n");
00337 }
00338 return cvm::get_error();
00339 }
00340 
00341 
00342 int colvarbias_restraint_centers_moving::update()
00343 {
00344 if (!cvm::main()->proxy->simulation_running()) {
00345 return COLVARS_OK;
00346 }
00347 
00348 if (b_chg_centers) {
00349 
00350 if (target_nstages) {
00351 // Staged update
00352 if (stage <= target_nstages) {
00353 if ((cvm::step_relative() > 0) &&
00354 (((cvm::step_absolute() - first_step) % target_nsteps) == 1)) {
00355 cvm::real const lambda =
00356 cvm::real(stage)/cvm::real(target_nstages);
00357 update_centers(lambda);
00358 stage++;
00359 cvm::log("Moving restraint \"" + this->name +
00360 "\" stage " + cvm::to_str(stage) +
00361 " : setting centers to " + cvm::to_str(colvar_centers) +
00362 " at step " + cvm::to_str(cvm::step_absolute()));
00363 } else {
00364 for (size_t i = 0; i < num_variables(); i++) {
00365 centers_incr[i].reset();
00366 }
00367 }
00368 }
00369 } else {
00370 // Continuous update
00371 if (cvm::step_absolute() - first_step <= target_nsteps) {
00372 cvm::real const lambda =
00373 cvm::real(cvm::step_absolute() - first_step)/cvm::real(target_nsteps);
00374 update_centers(lambda);
00375 } else {
00376 for (size_t i = 0; i < num_variables(); i++) {
00377 centers_incr[i].reset();
00378 }
00379 }
00380 }
00381 
00382 if (cvm::step_relative() == 0) {
00383 for (size_t i = 0; i < num_variables(); i++) {
00384 // finite differences are undefined when restarting
00385 centers_incr[i].reset();
00386 }
00387 }
00388 
00389 if (cvm::debug()) {
00390 cvm::log("Center increment for the restraint bias \""+
00391 this->name+"\": "+cvm::to_str(centers_incr)+
00392 " at stage "+cvm::to_str(stage)+ ".\n");
00393 }
00394 }
00395 
00396 return cvm::get_error();
00397 }
00398 
00399 
00400 int colvarbias_restraint_centers_moving::update_acc_work()
00401 {
00402 if (!cvm::main()->proxy->simulation_running()) {
00403 return COLVARS_OK;
00404 }
00405 if (b_chg_centers) {
00406 if (is_enabled(f_cvb_output_acc_work)) {
00407 if ((cvm::step_relative() > 0) &&
00408 (cvm::step_absolute() - first_step <= target_nsteps)) {
00409 for (size_t i = 0; i < num_variables(); i++) {
00410 // project forces on the calculated increments at this step
00411 acc_work += colvar_forces[i] * centers_incr[i];
00412 }
00413 }
00414 }
00415 }
00416 return COLVARS_OK;
00417 }
00418 
00419 
00420 std::string const colvarbias_restraint_centers_moving::get_state_params() const
00421 {
00422 std::ostringstream os;
00423 os.setf(std::ios::scientific, std::ios::floatfield);
00424 
00425 if (b_chg_centers) {
00426 size_t i;
00427 os << "centers ";
00428 for (i = 0; i < num_variables(); i++) {
00429 os << " "
00430 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
00431 << colvar_centers[i];
00432 }
00433 os << "\n";
00434 
00435 if (is_enabled(f_cvb_output_acc_work)) {
00436 os << "accumulatedWork "
00437 << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
00438 << acc_work << "\n";
00439 }
00440 }
00441 
00442 return os.str();
00443 }
00444 
00445 
00446 int colvarbias_restraint_centers_moving::set_state_params(std::string const &conf)
00447 {
00448 colvarbias_restraint::set_state_params(conf);
00449 
00450 if (b_chg_centers) {
00451 get_keyval(conf, "centers", colvar_centers, colvar_centers,
00452 colvarparse::parse_restart | colvarparse::parse_required);
00453 }
00454 
00455 if (is_enabled(f_cvb_output_acc_work)) {
00456 get_keyval(conf, "accumulatedWork", acc_work, acc_work,
00457 colvarparse::parse_restart | colvarparse::parse_required);
00458 }
00459 
00460 return COLVARS_OK;
00461 }
00462 
00463 
00464 std::ostream & colvarbias_restraint_centers_moving::write_traj_label(std::ostream &os)
00465 {
00466 if (b_output_centers) {
00467 for (size_t i = 0; i < num_variables(); i++) {
00468 size_t const this_cv_width = (variables(i)->value()).output_width(cvm::cv_width);
00469 os << " x0_"
00470 << cvm::wrap_string(variables(i)->name, this_cv_width-3);
00471 }
00472 }
00473 
00474 if (b_chg_centers && is_enabled(f_cvb_output_acc_work)) {
00475 os << " W_"
00476 << cvm::wrap_string(this->name, cvm::en_width-2);
00477 }
00478 
00479 return os;
00480 }
00481 
00482 
00483 std::ostream & colvarbias_restraint_centers_moving::write_traj(std::ostream &os)
00484 {
00485 if (b_output_centers) {
00486 for (size_t i = 0; i < num_variables(); i++) {
00487 os << " "
00488 << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width)
00489 << colvar_centers[i];
00490 }
00491 }
00492 
00493 if (b_chg_centers && is_enabled(f_cvb_output_acc_work)) {
00494 os << " "
00495 << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
00496 << acc_work;
00497 }
00498 
00499 return os;
00500 }
00501 
00502 
00503 
00504 colvarbias_restraint_k_moving::colvarbias_restraint_k_moving(char const *key)
00505 : colvarbias(key),
00506 colvarbias_ti(key),
00507 colvarbias_restraint(key),
00508 colvarbias_restraint_k(key),
00509 colvarbias_restraint_moving(key)
00510 {
00511 b_chg_force_k = false;
00512 b_decoupling = false;
00513 target_equil_steps = 0;
00514 target_force_k = -1.0;
00515 starting_force_k = -1.0;
00516 lambda_exp = 1.0;
00517 restraint_FE = 0.0;
00518 force_k_incr = 0.0;
00519 }
00520 
00521 
00522 int colvarbias_restraint_k_moving::init(std::string const &conf)
00523 {
00524 colvarbias_restraint_k::init(conf);
00525 
00526 get_keyval(conf, "decoupling", b_decoupling, b_decoupling);
00527 if (b_decoupling) {
00528 starting_force_k = 0.0;
00529 target_force_k = force_k;
00530 b_chg_force_k = true;
00531 }
00532 
00533 if (get_keyval(conf, "targetForceConstant", target_force_k, target_force_k)) {
00534 if (b_decoupling) {
00535 cvm::error("Error: targetForceConstant may not be specified together with decoupling.\n", COLVARS_INPUT_ERROR);
00536 return COLVARS_ERROR;
00537 }
00538 starting_force_k = force_k;
00539 b_chg_force_k = true;
00540 }
00541 
00542 if (!b_chg_force_k) {
00543 return COLVARS_OK;
00544 }
00545 
00546 // parse moving restraint options
00547 colvarbias_restraint_moving::init(conf);
00548 
00549 get_keyval(conf, "targetEquilSteps", target_equil_steps, target_equil_steps);
00550 
00551 if (get_keyval(conf, "lambdaSchedule", lambda_schedule, lambda_schedule) &&
00552 target_nstages > 0) {
00553 cvm::error("Error: targetNumStages and lambdaSchedule are incompatible.\n", COLVARS_INPUT_ERROR);
00554 return cvm::get_error();
00555 }
00556 
00557 if (lambda_schedule.size()) {
00558 // There is one more lambda-point than stages
00559 target_nstages = lambda_schedule.size() - 1;
00560 }
00561 
00562 if ((get_keyval(conf, "targetForceExponent", lambda_exp, lambda_exp, parse_deprecated)
00563 || get_keyval(conf, "lambdaExponent", lambda_exp, lambda_exp))
00564 && !b_chg_force_k) {
00565 cvm::error("Error: cannot set lambdaExponent unless a changing force constant is active.\n", COLVARS_INPUT_ERROR);
00566 }
00567 if (lambda_exp < 1.0) {
00568 cvm::log("Warning: for all practical purposes, lambdaExponent should be 1.0 or greater.\n");
00569 }
00570 
00571 return COLVARS_OK;
00572 }
00573 
00574 
00575 int colvarbias_restraint_k_moving::update()
00576 {
00577 if (!cvm::main()->proxy->simulation_running()) {
00578 return COLVARS_OK;
00579 }
00580 if (b_chg_force_k) {
00581 
00582 cvm::real lambda;
00583 
00584 if (target_nstages) {
00585 
00586 if (cvm::step_absolute() == first_step) {
00587 // Setup first stage of staged variable force constant calculation
00588 if (lambda_schedule.size()) {
00589 lambda = lambda_schedule[0];
00590 } else {
00591 lambda = (b_decoupling ? 1.0 : 0.0);
00592 }
00593 force_k = starting_force_k + (target_force_k - starting_force_k)
00594 * cvm::pow(lambda, lambda_exp);
00595 cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
00596 + " : lambda = " + cvm::to_str(lambda)
00597 + ", k = " + cvm::to_str(force_k)+"\n");
00598 }
00599 
00600 // TI calculation: estimate free energy derivative
00601 // need current lambda
00602 if (lambda_schedule.size()) {
00603 lambda = lambda_schedule[stage];
00604 } else {
00605 lambda = cvm::real(stage) / cvm::real(target_nstages);
00606 if (b_decoupling) lambda = 1.0 - lambda;
00607 }
00608 
00609 if (target_equil_steps == 0 || (cvm::step_absolute() - first_step) % target_nsteps >= target_equil_steps) {
00610 // Start averaging after equilibration period, if requested
00611 
00612 // Derivative of energy with respect to force_k
00613 cvm::real dU_dk = 0.0;
00614 for (size_t i = 0; i < num_variables(); i++) {
00615 dU_dk += d_restraint_potential_dk(i);
00616 }
00617 restraint_FE += lambda_exp * cvm::pow(lambda, lambda_exp - 1.0)
00618 * (target_force_k - starting_force_k) * dU_dk;
00619 }
00620 
00621 // Finish current stage...
00622 if ((cvm::step_absolute() - first_step) % target_nsteps == 0 &&
00623 cvm::step_absolute() > first_step) {
00624 
00625 cvm::log("Restraint " + this->name + " Lambda= "
00626 + cvm::to_str(lambda) + " dA/dLambda= "
00627 + cvm::to_str(restraint_FE / cvm::real(target_nsteps - target_equil_steps))+"\n");
00628 
00629 // ...and move on to the next one
00630 if (stage < target_nstages) {
00631 
00632 restraint_FE = 0.0;
00633 stage++;
00634 if (lambda_schedule.size()) {
00635 lambda = lambda_schedule[stage];
00636 } else {
00637 lambda = cvm::real(stage) / cvm::real(target_nstages);
00638 if (b_decoupling) lambda = 1.0 - lambda;
00639 }
00640 force_k = starting_force_k + (target_force_k - starting_force_k)
00641 * cvm::pow(lambda, lambda_exp);
00642 cvm::log("Restraint " + this->name + ", stage " + cvm::to_str(stage)
00643 + " : lambda = " + cvm::to_str(lambda)
00644 + ", k = " + cvm::to_str(force_k)+"\n");
00645 }
00646 }
00647 
00648 } else if (cvm::step_absolute() - first_step <= target_nsteps) {
00649 
00650 // update force constant (slow growth)
00651 lambda = cvm::real(cvm::step_absolute() - first_step) / cvm::real(target_nsteps);
00652 if (b_decoupling) lambda = 1.0 - lambda;
00653 cvm::real const force_k_old = force_k;
00654 force_k = starting_force_k + (target_force_k - starting_force_k)
00655 * cvm::pow(lambda, lambda_exp);
00656 force_k_incr = force_k - force_k_old;
00657 }
00658 }
00659 
00660 return COLVARS_OK;
00661 }
00662 
00663 
00664 int colvarbias_restraint_k_moving::update_acc_work()
00665 {
00666 if (!cvm::main()->proxy->simulation_running()) {
00667 return COLVARS_OK;
00668 }
00669 if (b_chg_force_k) {
00670 if (is_enabled(f_cvb_output_acc_work)) {
00671 if (cvm::step_relative() > 0) {
00672 cvm::real dU_dk = 0.0;
00673 for (size_t i = 0; i < num_variables(); i++) {
00674 dU_dk += d_restraint_potential_dk(i);
00675 }
00676 acc_work += dU_dk * force_k_incr;
00677 }
00678 }
00679 }
00680 return COLVARS_OK;
00681 }
00682 
00683 
00684 std::string const colvarbias_restraint_k_moving::get_state_params() const
00685 {
00686 std::ostringstream os;
00687 os.setf(std::ios::scientific, std::ios::floatfield);
00688 if (b_chg_force_k) {
00689 os << "forceConstant "
00690 << std::setprecision(cvm::en_prec)
00691 << std::setw(cvm::en_width) << force_k << "\n";
00692 
00693 if (is_enabled(f_cvb_output_acc_work)) {
00694 os << "accumulatedWork "
00695 << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
00696 << acc_work << "\n";
00697 }
00698 }
00699 return os.str();
00700 }
00701 
00702 
00703 int colvarbias_restraint_k_moving::set_state_params(std::string const &conf)
00704 {
00705 colvarbias_restraint::set_state_params(conf);
00706 
00707 if (b_chg_force_k) {
00708 get_keyval(conf, "forceConstant", force_k, force_k,
00709 colvarparse::parse_restart | colvarparse::parse_required);
00710 }
00711 
00712 if (is_enabled(f_cvb_output_acc_work)) {
00713 get_keyval(conf, "accumulatedWork", acc_work, acc_work,
00714 colvarparse::parse_restart | colvarparse::parse_required);
00715 }
00716 
00717 return COLVARS_OK;
00718 }
00719 
00720 
00721 std::ostream & colvarbias_restraint_k_moving::write_traj_label(std::ostream &os)
00722 {
00723 if (b_chg_force_k && is_enabled(f_cvb_output_acc_work)) {
00724 os << " W_"
00725 << cvm::wrap_string(this->name, cvm::en_width-2);
00726 }
00727 return os;
00728 }
00729 
00730 
00731 std::ostream & colvarbias_restraint_k_moving::write_traj(std::ostream &os)
00732 {
00733 if (b_chg_force_k && is_enabled(f_cvb_output_acc_work)) {
00734 os << " "
00735 << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
00736 << acc_work;
00737 }
00738 return os;
00739 }
00740 
00741 
00742 
00743 colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key)
00744 : colvarbias(key),
00745 colvarbias_ti(key),
00746 colvarbias_restraint(key),
00747 colvarbias_restraint_centers(key),
00748 colvarbias_restraint_moving(key),
00749 colvarbias_restraint_k(key),
00750 colvarbias_restraint_centers_moving(key),
00751 colvarbias_restraint_k_moving(key)
00752 {
00753 }
00754 
00755 
00756 int colvarbias_restraint_harmonic::init(std::string const &conf)
00757 {
00758 colvarbias_restraint::init(conf);
00759 colvarbias_restraint_moving::init(conf);
00760 colvarbias_restraint_centers_moving::init(conf);
00761 colvarbias_restraint_k_moving::init(conf);
00762 
00763 cvm::main()->cite_feature("Harmonic colvar bias implementation");
00764 
00765 for (size_t i = 0; i < num_variables(); i++) {
00766 cvm::real const w = variables(i)->width;
00767 cvm::log("The force constant for colvar \""+variables(i)->name+
00768 "\" will be rescaled to "+
00769 cvm::to_str(force_k/(w*w))+
00770 " according to the specified width ("+cvm::to_str(w)+").\n");
00771 }
00772 
00773 return COLVARS_OK;
00774 }
00775 
00776 
00777 int colvarbias_restraint_harmonic::update()
00778 {
00779 int error_code = COLVARS_OK;
00780 
00781 // update the TI estimator (if defined)
00782 error_code |= colvarbias_ti::update();
00783 
00784 // update parameters (centers or force constant)
00785 error_code |= colvarbias_restraint_centers_moving::update();
00786 error_code |= colvarbias_restraint_k_moving::update();
00787 
00788 // update restraint energy and forces
00789 error_code |= colvarbias_restraint::update();
00790 
00791 // update accumulated work using the current forces
00792 error_code |= colvarbias_restraint_centers_moving::update_acc_work();
00793 error_code |= colvarbias_restraint_k_moving::update_acc_work();
00794 
00795 return error_code;
00796 }
00797 
00798 
00799 cvm::real colvarbias_restraint_harmonic::restraint_potential(size_t i) const
00800 {
00801 return 0.5 * force_k / (variables(i)->width * variables(i)->width) *
00802 variables(i)->dist2(variables(i)->value(), colvar_centers[i]);
00803 }
00804 
00805 
00806 colvarvalue const colvarbias_restraint_harmonic::restraint_force(size_t i) const
00807 {
00808 return -0.5 * force_k / (variables(i)->width * variables(i)->width) *
00809 variables(i)->dist2_lgrad(variables(i)->value(), colvar_centers[i]);
00810 }
00811 
00812 
00813 cvm::real colvarbias_restraint_harmonic::d_restraint_potential_dk(size_t i) const
00814 {
00815 return 0.5 / (variables(i)->width * variables(i)->width) *
00816 variables(i)->dist2(variables(i)->value(), colvar_centers[i]);
00817 }
00818 
00819 
00820 std::string const colvarbias_restraint_harmonic::get_state_params() const
00821 {
00822 return colvarbias_restraint::get_state_params() +
00823 colvarbias_restraint_moving::get_state_params() +
00824 colvarbias_restraint_centers_moving::get_state_params() +
00825 colvarbias_restraint_k_moving::get_state_params();
00826 }
00827 
00828 
00829 int colvarbias_restraint_harmonic::set_state_params(std::string const &conf)
00830 {
00831 int error_code = COLVARS_OK;
00832 error_code |= colvarbias_restraint::set_state_params(conf);
00833 error_code |= colvarbias_restraint_moving::set_state_params(conf);
00834 error_code |= colvarbias_restraint_centers_moving::set_state_params(conf);
00835 error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
00836 return error_code;
00837 }
00838 
00839 
00840 std::ostream & colvarbias_restraint_harmonic::write_state_data(std::ostream &os)
00841 {
00842 return colvarbias_ti::write_state_data(os);
00843 }
00844 
00845 
00846 std::istream & colvarbias_restraint_harmonic::read_state_data(std::istream &is)
00847 {
00848 return colvarbias_ti::read_state_data(is);
00849 }
00850 
00851 
00852 std::ostream & colvarbias_restraint_harmonic::write_traj_label(std::ostream &os)
00853 {
00854 colvarbias_restraint::write_traj_label(os);
00855 colvarbias_restraint_centers_moving::write_traj_label(os);
00856 colvarbias_restraint_k_moving::write_traj_label(os);
00857 return os;
00858 }
00859 
00860 
00861 std::ostream & colvarbias_restraint_harmonic::write_traj(std::ostream &os)
00862 {
00863 colvarbias_restraint::write_traj(os);
00864 colvarbias_restraint_centers_moving::write_traj(os);
00865 colvarbias_restraint_k_moving::write_traj(os);
00866 return os;
00867 }
00868 
00869 
00870 int colvarbias_restraint_harmonic::change_configuration(std::string const &conf)
00871 {
00872 return colvarbias_restraint_centers::change_configuration(conf) |
00873 colvarbias_restraint_k::change_configuration(conf);
00874 }
00875 
00876 
00877 cvm::real colvarbias_restraint_harmonic::energy_difference(std::string const &conf)
00878 {
00879 cvm::real const old_bias_energy = bias_energy;
00880 cvm::real const old_force_k = force_k;
00881 std::vector<colvarvalue> const old_centers = colvar_centers;
00882 
00883 change_configuration(conf);
00884 update();
00885 
00886 cvm::real const result = (bias_energy - old_bias_energy);
00887 
00888 bias_energy = old_bias_energy;
00889 force_k = old_force_k;
00890 colvar_centers = old_centers;
00891 
00892 return result;
00893 }
00894 
00895 
00896 
00897 colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char const *key)
00898 : colvarbias(key),
00899 colvarbias_ti(key),
00900 colvarbias_restraint(key),
00901 colvarbias_restraint_k(key),
00902 colvarbias_restraint_moving(key),
00903 colvarbias_restraint_k_moving(key)
00904 {
00905 lower_wall_k = -1.0;
00906 upper_wall_k = -1.0;
00907 // This bias implements the bias_actual_colvars feature (most others do not)
00908 provide(f_cvb_bypass_ext_lagrangian);
00909 set_enabled(f_cvb_bypass_ext_lagrangian); // Defaults to enabled
00910 }
00911 
00912 
00913 int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
00914 {
00915 colvarbias_restraint::init(conf);
00916 colvarbias_restraint_moving::init(conf);
00917 colvarbias_restraint_k_moving::init(conf);
00918 
00919 cvm::main()->cite_feature("harmonicWalls colvar bias implementation");
00920 
00921 enable(f_cvb_scalar_variables);
00922 
00923 size_t i;
00924 
00925 bool b_null_lower_walls = false;
00926 if (lower_walls.size() == 0) {
00927 b_null_lower_walls = true;
00928 lower_walls.resize(num_variables());
00929 for (i = 0; i < num_variables(); i++) {
00930 lower_walls[i].type(variables(i)->value());
00931 lower_walls[i].reset();
00932 }
00933 }
00934 if (!get_keyval(conf, "lowerWalls", lower_walls, lower_walls) &&
00935 b_null_lower_walls) {
00936 cvm::log("Lower walls were not provided.\n");
00937 lower_walls.clear();
00938 }
00939 
00940 bool b_null_upper_walls = false;
00941 if (upper_walls.size() == 0) {
00942 b_null_upper_walls = true;
00943 upper_walls.resize(num_variables());
00944 for (i = 0; i < num_variables(); i++) {
00945 upper_walls[i].type(variables(i)->value());
00946 upper_walls[i].reset();
00947 }
00948 }
00949 if (!get_keyval(conf, "upperWalls", upper_walls, upper_walls) &&
00950 b_null_upper_walls) {
00951 cvm::log("Upper walls were not provided.\n");
00952 upper_walls.clear();
00953 }
00954 
00955 if ((lower_walls.size() == 0) && (upper_walls.size() == 0)) {
00956 return cvm::error("Error: no walls provided.\n", COLVARS_INPUT_ERROR);
00957 }
00958 
00959 if (lower_walls.size() > 0) {
00960 get_keyval(conf, "lowerWallConstant", lower_wall_k,
00961 (lower_wall_k > 0.0) ? lower_wall_k : force_k);
00962 }
00963 if (upper_walls.size() > 0) {
00964 get_keyval(conf, "upperWallConstant", upper_wall_k,
00965 (upper_wall_k > 0.0) ? upper_wall_k : force_k);
00966 }
00967 
00968 if ((lower_walls.size() == 0) || (upper_walls.size() == 0)) {
00969 for (i = 0; i < num_variables(); i++) {
00970 if (variables(i)->is_enabled(f_cv_periodic)) {
00971 return cvm::error("Error: at least one variable is periodic, "
00972 "both walls must be provided.\n", COLVARS_INPUT_ERROR);
00973 }
00974 }
00975 }
00976 
00977 if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) {
00978 for (i = 0; i < num_variables(); i++) {
00979 if (lower_walls[i] >= upper_walls[i]) {
00980 return cvm::error("Error: one upper wall, "+
00981 cvm::to_str(upper_walls[i])+
00982 ", is not higher than the lower wall, "+
00983 cvm::to_str(lower_walls[i])+".\n",
00984 COLVARS_INPUT_ERROR);
00985 }
00986 if (variables(i)->dist2(lower_walls[i], upper_walls[i]) < 1.0e-12) {
00987 return cvm::error("Error: lower wall and upper wall are equal "
00988 "in the domain of the variable \""+
00989 variables(i)->name+"\".\n", COLVARS_INPUT_ERROR);
00990 }
00991 }
00992 if (lower_wall_k * upper_wall_k == 0.0) {
00993 cvm::error("Error: lowerWallConstant and upperWallConstant, "
00994 "when defined, must both be positive.\n",
00995 COLVARS_INPUT_ERROR);
00996 return COLVARS_INPUT_ERROR;
00997 }
00998 force_k = cvm::sqrt(lower_wall_k * upper_wall_k);
00999 // transform the two constants to relative values using gemetric mean as ref
01000 // to preserve force_k if provided as single parameter
01001 // (allow changing both via force_k)
01002 lower_wall_k /= force_k;
01003 upper_wall_k /= force_k;
01004 } else {
01005 // If only one wall is defined, need to rescale as well
01006 if (lower_walls.size() > 0) {
01007 force_k = lower_wall_k;
01008 lower_wall_k = 1.0;
01009 }
01010 if (upper_walls.size() > 0) {
01011 force_k = upper_wall_k;
01012 upper_wall_k = 1.0;
01013 }
01014 }
01015 
01016 // Initialize starting value of the force constant (in case it's changing)
01017 starting_force_k = (b_decoupling ? 0.0 : force_k);
01018 
01019 if (lower_walls.size() > 0) {
01020 for (i = 0; i < num_variables(); i++) {
01021 cvm::real const w = variables(i)->width;
01022 cvm::log("The lower wall force constant for colvar \""+
01023 variables(i)->name+"\" will be rescaled to "+
01024 cvm::to_str(lower_wall_k * force_k / (w*w))+
01025 " according to the specified width ("+cvm::to_str(w)+").\n");
01026 }
01027 }
01028 
01029 if (upper_walls.size() > 0) {
01030 for (i = 0; i < num_variables(); i++) {
01031 cvm::real const w = variables(i)->width;
01032 cvm::log("The upper wall force constant for colvar \""+
01033 variables(i)->name+"\" will be rescaled to "+
01034 cvm::to_str(upper_wall_k * force_k / (w*w))+
01035 " according to the specified width ("+cvm::to_str(w)+").\n");
01036 }
01037 }
01038 
01039 return COLVARS_OK;
01040 }
01041 
01042 
01043 int colvarbias_restraint_harmonic_walls::update()
01044 {
01045 int error_code = COLVARS_OK;
01046 
01047 error_code |= colvarbias_ti::update();
01048 
01049 error_code |= colvarbias_restraint_k_moving::update();
01050 
01051 error_code |= colvarbias_restraint::update();
01052 
01053 error_code |= colvarbias_restraint_k_moving::update_acc_work();
01054 
01055 return error_code;
01056 }
01057 
01058 
01059 cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const
01060 {
01061 colvar *cv = variables(i);
01062 
01063 colvarvalue const &cvv = is_enabled(f_cvb_bypass_ext_lagrangian) ?
01064 variables(i)->actual_value() :
01065 variables(i)->value();
01066 
01067 // For a periodic colvar, both walls may be applicable at the same time
01068 // in which case we pick the closer one
01069 
01070 if (cv->is_enabled(f_cv_periodic)) {
01071 cvm::real const lower_wall_dist2 = cv->dist2(cvv, lower_walls[i]);
01072 cvm::real const upper_wall_dist2 = cv->dist2(cvv, upper_walls[i]);
01073 if (lower_wall_dist2 < upper_wall_dist2) {
01074 cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]);
01075 if (grad < 0.0) { return 0.5 * grad; }
01076 } else {
01077 cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]);
01078 if (grad > 0.0) { return 0.5 * grad; }
01079 }
01080 return 0.0;
01081 }
01082 
01083 if (lower_walls.size() > 0) {
01084 cvm::real const grad = cv->dist2_lgrad(cvv, lower_walls[i]);
01085 if (grad < 0.0) { return 0.5 * grad; }
01086 }
01087 if (upper_walls.size() > 0) {
01088 cvm::real const grad = cv->dist2_lgrad(cvv, upper_walls[i]);
01089 if (grad > 0.0) { return 0.5 * grad; }
01090 }
01091 return 0.0;
01092 }
01093 
01094 
01095 cvm::real colvarbias_restraint_harmonic_walls::restraint_potential(size_t i) const
01096 {
01097 cvm::real const dist = colvar_distance(i);
01098 cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k;
01099 return 0.5 * force_k * scale / (variables(i)->width * variables(i)->width) *
01100 dist * dist;
01101 }
01102 
01103 
01104 colvarvalue const colvarbias_restraint_harmonic_walls::restraint_force(size_t i) const
01105 {
01106 cvm::real const dist = colvar_distance(i);
01107 cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k;
01108 return - force_k * scale / (variables(i)->width * variables(i)->width) * dist;
01109 }
01110 
01111 
01112 cvm::real colvarbias_restraint_harmonic_walls::d_restraint_potential_dk(size_t i) const
01113 {
01114 cvm::real const dist = colvar_distance(i);
01115 cvm::real const scale = dist > 0.0 ? upper_wall_k : lower_wall_k;
01116 return 0.5 * scale / (variables(i)->width * variables(i)->width) *
01117 dist * dist;
01118 }
01119 
01120 
01121 std::string const colvarbias_restraint_harmonic_walls::get_state_params() const
01122 {
01123 return colvarbias_restraint::get_state_params() +
01124 colvarbias_restraint_moving::get_state_params() +
01125 colvarbias_restraint_k_moving::get_state_params();
01126 }
01127 
01128 
01129 int colvarbias_restraint_harmonic_walls::set_state_params(std::string const &conf)
01130 {
01131 int error_code = COLVARS_OK;
01132 error_code |= colvarbias_restraint::set_state_params(conf);
01133 error_code |= colvarbias_restraint_moving::set_state_params(conf);
01134 error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
01135 return error_code;
01136 }
01137 
01138 
01139 std::ostream & colvarbias_restraint_harmonic_walls::write_state_data(std::ostream &os)
01140 {
01141 return colvarbias_ti::write_state_data(os);
01142 }
01143 
01144 
01145 std::istream & colvarbias_restraint_harmonic_walls::read_state_data(std::istream &is)
01146 {
01147 return colvarbias_ti::read_state_data(is);
01148 }
01149 
01150 
01151 std::ostream & colvarbias_restraint_harmonic_walls::write_traj_label(std::ostream &os)
01152 {
01153 colvarbias_restraint::write_traj_label(os);
01154 colvarbias_restraint_k_moving::write_traj_label(os);
01155 return os;
01156 }
01157 
01158 
01159 std::ostream & colvarbias_restraint_harmonic_walls::write_traj(std::ostream &os)
01160 {
01161 colvarbias_restraint::write_traj(os);
01162 colvarbias_restraint_k_moving::write_traj(os);
01163 return os;
01164 }
01165 
01166 
01167 
01168 colvarbias_restraint_linear::colvarbias_restraint_linear(char const *key)
01169 : colvarbias(key),
01170 colvarbias_ti(key),
01171 colvarbias_restraint(key),
01172 colvarbias_restraint_centers(key),
01173 colvarbias_restraint_moving(key),
01174 colvarbias_restraint_k(key),
01175 colvarbias_restraint_centers_moving(key),
01176 colvarbias_restraint_k_moving(key)
01177 {
01178 check_positive_k = false;
01179 }
01180 
01181 
01182 int colvarbias_restraint_linear::init(std::string const &conf)
01183 {
01184 colvarbias_restraint::init(conf);
01185 colvarbias_restraint_moving::init(conf);
01186 colvarbias_restraint_centers_moving::init(conf);
01187 colvarbias_restraint_k_moving::init(conf);
01188 
01189 cvm::main()->cite_feature("harmonicWalls colvar bias implementation");
01190 
01191 for (size_t i = 0; i < num_variables(); i++) {
01192 if (variables(i)->is_enabled(f_cv_periodic)) {
01193 cvm::error("Error: linear biases cannot be applied to periodic variables.\n",
01194 COLVARS_INPUT_ERROR);
01195 return COLVARS_INPUT_ERROR;
01196 }
01197 cvm::real const w = variables(i)->width;
01198 cvm::log("The force constant for colvar \""+variables(i)->name+
01199 "\" will be rescaled to "+
01200 cvm::to_str(force_k / w)+
01201 " according to the specified width ("+cvm::to_str(w)+").\n");
01202 }
01203 
01204 return COLVARS_OK;
01205 }
01206 
01207 
01208 int colvarbias_restraint_linear::update()
01209 {
01210 int error_code = COLVARS_OK;
01211 
01212 // update the TI estimator (if defined)
01213 error_code |= colvarbias_ti::update();
01214 
01215 // update parameters (centers or force constant)
01216 error_code |= colvarbias_restraint_centers_moving::update();
01217 error_code |= colvarbias_restraint_k_moving::update();
01218 
01219 // update restraint energy and forces
01220 error_code |= colvarbias_restraint::update();
01221 
01222 // update accumulated work using the current forces
01223 error_code |= colvarbias_restraint_centers_moving::update_acc_work();
01224 error_code |= colvarbias_restraint_k_moving::update_acc_work();
01225 
01226 return error_code;
01227 }
01228 
01229 
01230 int colvarbias_restraint_linear::change_configuration(std::string const &conf)
01231 {
01232 // Only makes sense to change the force constant
01233 return colvarbias_restraint_k::change_configuration(conf);
01234 }
01235 
01236 
01237 cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf)
01238 {
01239 cvm::real const old_bias_energy = bias_energy;
01240 cvm::real const old_force_k = force_k;
01241 
01242 change_configuration(conf);
01243 update();
01244 
01245 cvm::real const result = (bias_energy - old_bias_energy);
01246 
01247 bias_energy = old_bias_energy;
01248 force_k = old_force_k;
01249 
01250 return result;
01251 }
01252 
01253 
01254 cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const
01255 {
01256 return force_k / variables(i)->width * (variables(i)->value() -
01257 colvar_centers[i]).sum();
01258 }
01259 
01260 
01261 colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const
01262 {
01263 colvarvalue dummy(variables(i)->value());
01264 dummy.set_ones();
01265 return -1.0 * force_k / variables(i)->width * dummy;
01266 }
01267 
01268 
01269 cvm::real colvarbias_restraint_linear::d_restraint_potential_dk(size_t i) const
01270 {
01271 return 1.0 / variables(i)->width * (variables(i)->value() -
01272 colvar_centers[i]).sum();
01273 }
01274 
01275 
01276 std::string const colvarbias_restraint_linear::get_state_params() const
01277 {
01278 return colvarbias_restraint::get_state_params() +
01279 colvarbias_restraint_moving::get_state_params() +
01280 colvarbias_restraint_centers_moving::get_state_params() +
01281 colvarbias_restraint_k_moving::get_state_params();
01282 }
01283 
01284 
01285 int colvarbias_restraint_linear::set_state_params(std::string const &conf)
01286 {
01287 int error_code = COLVARS_OK;
01288 error_code |= colvarbias_restraint::set_state_params(conf);
01289 error_code |= colvarbias_restraint_moving::set_state_params(conf);
01290 error_code |= colvarbias_restraint_centers_moving::set_state_params(conf);
01291 error_code |= colvarbias_restraint_k_moving::set_state_params(conf);
01292 return error_code;
01293 }
01294 
01295 
01296 std::ostream & colvarbias_restraint_linear::write_state_data(std::ostream &os)
01297 {
01298 return colvarbias_ti::write_state_data(os);
01299 }
01300 
01301 
01302 std::istream & colvarbias_restraint_linear::read_state_data(std::istream &is)
01303 {
01304 return colvarbias_ti::read_state_data(is);
01305 }
01306 
01307 
01308 std::ostream & colvarbias_restraint_linear::write_traj_label(std::ostream &os)
01309 {
01310 colvarbias_restraint::write_traj_label(os);
01311 colvarbias_restraint_centers_moving::write_traj_label(os);
01312 colvarbias_restraint_k_moving::write_traj_label(os);
01313 return os;
01314 }
01315 
01316 
01317 std::ostream & colvarbias_restraint_linear::write_traj(std::ostream &os)
01318 {
01319 colvarbias_restraint::write_traj(os);
01320 colvarbias_restraint_centers_moving::write_traj(os);
01321 colvarbias_restraint_k_moving::write_traj(os);
01322 return os;
01323 }
01324 
01325 
01326 
01327 colvarbias_restraint_histogram::colvarbias_restraint_histogram(char const *key)
01328 : colvarbias(key)
01329 {
01330 lower_boundary = 0.0;
01331 upper_boundary = 0.0;
01332 width = 0.0;
01333 gaussian_width = 0.0;
01334 }
01335 
01336 
01337 int colvarbias_restraint_histogram::init(std::string const &conf)
01338 {
01339 int error_code = COLVARS_OK;
01340 
01341 colvarbias::init(conf);
01342 enable(f_cvb_apply_force);
01343 
01344 cvm::main()->cite_feature("histogramRestraint colvar bias implementation");
01345 
01346 get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary);
01347 get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary);
01348 get_keyval(conf, "width", width, width);
01349 
01350 if (width <= 0.0) {
01351 error_code |= cvm::error("Error: \"width\" must be positive.\n",
01352 COLVARS_INPUT_ERROR);
01353 }
01354 
01355 get_keyval(conf, "gaussianWidth", gaussian_width, 2.0 * width, colvarparse::parse_silent);
01356 get_keyval(conf, "gaussianSigma", gaussian_width, 2.0 * width);
01357 
01358 if (lower_boundary >= upper_boundary) {
01359 error_code |= cvm::error("Error: the upper boundary, "+
01360 cvm::to_str(upper_boundary)+
01361 ", is not higher than the lower boundary, "+
01362 cvm::to_str(lower_boundary)+".\n",
01363 COLVARS_INPUT_ERROR);
01364 }
01365 
01366 cvm::real const nbins = (upper_boundary - lower_boundary) / width;
01367 int const nbins_round = (int)(nbins);
01368 
01369 if (cvm::fabs(nbins - cvm::real(nbins_round)) > 1.0E-10) {
01370 cvm::log("Warning: grid interval ("+
01371 cvm::to_str(lower_boundary, cvm::cv_width, cvm::cv_prec)+" - "+
01372 cvm::to_str(upper_boundary, cvm::cv_width, cvm::cv_prec)+
01373 ") is not commensurate to its bin width ("+
01374 cvm::to_str(width, cvm::cv_width, cvm::cv_prec)+").\n");
01375 }
01376 
01377 p.resize(nbins_round);
01378 ref_p.resize(nbins_round);
01379 p_diff.resize(nbins_round);
01380 
01381 bool const inline_ref_p =
01382 get_keyval(conf, "refHistogram", ref_p.data_array(), ref_p.data_array());
01383 std::string ref_p_file;
01384 get_keyval(conf, "refHistogramFile", ref_p_file, std::string(""));
01385 if (ref_p_file.size()) {
01386 if (inline_ref_p) {
01387 error_code |= cvm::error("Error: cannot specify both refHistogram and refHistogramFile at the same time.\n",
01388 COLVARS_INPUT_ERROR);
01389 } else {
01390 
01391 std::istream &is =
01392 cvm::main()->proxy->input_stream(ref_p_file,
01393 "reference histogram file");
01394 
01395 std::string data_s = "";
01396 std::string line;
01397 while (getline_nocomments(is, line)) {
01398 data_s.append(line+"\n");
01399 }
01400 if (data_s.size() == 0) {
01401 error_code |= cvm::error("Error: file \""+ref_p_file+
01402 "\" empty or unreadable.\n",
01403 COLVARS_FILE_ERROR);
01404 }
01405 error_code |= cvm::main()->proxy->close_input_stream(ref_p_file);
01406 
01407 cvm::vector1d<cvm::real> data;
01408 if (data.from_simple_string(data_s) != 0) {
01409 error_code |= cvm::error("Error: could not read histogram from file \""+
01410 ref_p_file+"\".\n");
01411 }
01412 if (data.size() == 2*ref_p.size()) {
01413 // file contains both x and p(x)
01414 size_t i;
01415 for (i = 0; i < ref_p.size(); i++) {
01416 ref_p[i] = data[2*i+1];
01417 }
01418 } else if (data.size() == ref_p.size()) {
01419 ref_p = data;
01420 } else {
01421 error_code |= cvm::error("Error: file \""+ref_p_file+
01422 "\" contains a histogram of different length.\n",
01423 COLVARS_INPUT_ERROR);
01424 }
01425 }
01426 }
01427 
01428 cvm::real const ref_integral = ref_p.sum() * width;
01429 if (cvm::fabs(ref_integral - 1.0) > 1.0e-03) {
01430 cvm::log("Reference distribution not normalized, normalizing to unity.\n");
01431 ref_p /= ref_integral;
01432 }
01433 
01434 get_keyval(conf, "writeHistogram", b_write_histogram, false);
01435 get_keyval(conf, "forceConstant", force_k, 1.0);
01436 
01437 return error_code;
01438 }
01439 
01440 
01441 colvarbias_restraint_histogram::~colvarbias_restraint_histogram()
01442 {
01443 p.clear();
01444 ref_p.clear();
01445 p_diff.clear();
01446 }
01447 
01448 
01449 int colvarbias_restraint_histogram::update()
01450 {
01451 if (cvm::debug())
01452 cvm::log("Updating the histogram restraint bias \""+this->name+"\".\n");
01453 
01454 size_t vector_size = 0;
01455 size_t icv;
01456 for (icv = 0; icv < num_variables(); icv++) {
01457 vector_size += variables(icv)->value().size();
01458 }
01459 
01460 cvm::real const norm = 1.0/(cvm::sqrt(2.0*PI)*gaussian_width*vector_size);
01461 
01462 // calculate the histogram
01463 p.reset();
01464 for (icv = 0; icv < num_variables(); icv++) {
01465 colvarvalue const &cv = variables(icv)->value();
01466 if (cv.type() == colvarvalue::type_scalar) {
01467 cvm::real const cv_value = cv.real_value;
01468 size_t igrid;
01469 for (igrid = 0; igrid < p.size(); igrid++) {
01470 cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
01471 p[igrid] += norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
01472 (2.0 * gaussian_width * gaussian_width));
01473 }
01474 } else if (cv.type() == colvarvalue::type_vector) {
01475 size_t idim;
01476 for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
01477 cvm::real const cv_value = cv.vector1d_value[idim];
01478 size_t igrid;
01479 for (igrid = 0; igrid < p.size(); igrid++) {
01480 cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
01481 p[igrid] += norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
01482 (2.0 * gaussian_width * gaussian_width));
01483 }
01484 }
01485 } else {
01486 cvm::error("Error: unsupported type for variable "+variables(icv)->name+".\n",
01487 COLVARS_NOT_IMPLEMENTED);
01488 return COLVARS_NOT_IMPLEMENTED;
01489 }
01490 }
01491 
01492 cvm::real const force_k_cv = force_k * vector_size;
01493 
01494 // calculate the difference between current and reference
01495 p_diff = p - ref_p;
01496 bias_energy = 0.5 * force_k_cv * p_diff * p_diff;
01497 
01498 // calculate the forces
01499 for (icv = 0; icv < num_variables(); icv++) {
01500 colvarvalue const &cv = variables(icv)->value();
01501 colvarvalue &cv_force = colvar_forces[icv];
01502 cv_force.type(cv);
01503 cv_force.reset();
01504 
01505 if (cv.type() == colvarvalue::type_scalar) {
01506 cvm::real const cv_value = cv.real_value;
01507 cvm::real &force = cv_force.real_value;
01508 size_t igrid;
01509 for (igrid = 0; igrid < p.size(); igrid++) {
01510 cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
01511 force += force_k_cv * p_diff[igrid] *
01512 norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
01513 (2.0 * gaussian_width * gaussian_width)) *
01514 (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
01515 }
01516 } else if (cv.type() == colvarvalue::type_vector) {
01517 size_t idim;
01518 for (idim = 0; idim < cv.vector1d_value.size(); idim++) {
01519 cvm::real const cv_value = cv.vector1d_value[idim];
01520 cvm::real &force = cv_force.vector1d_value[idim];
01521 size_t igrid;
01522 for (igrid = 0; igrid < p.size(); igrid++) {
01523 cvm::real const x_grid = (lower_boundary + (igrid+0.5)*width);
01524 force += force_k_cv * p_diff[igrid] *
01525 norm * cvm::exp(-1.0 * (x_grid - cv_value) * (x_grid - cv_value) /
01526 (2.0 * gaussian_width * gaussian_width)) *
01527 (-1.0 * (x_grid - cv_value) / (gaussian_width * gaussian_width));
01528 }
01529 }
01530 } else {
01531 // TODO
01532 }
01533 }
01534 
01535 return COLVARS_OK;
01536 }
01537 
01538 
01539 int colvarbias_restraint_histogram::write_output_files()
01540 {
01541 if (b_write_histogram) {
01542 colvarproxy *proxy = cvm::main()->proxy;
01543 std::string file_name(cvm::output_prefix()+"."+this->name+".hist.dat");
01544 std::ostream &os = proxy->output_stream(file_name,
01545 "histogram output file");
01546 os << "# " << cvm::wrap_string(variables(0)->name, cvm::cv_width)
01547 << " " << "p(" << cvm::wrap_string(variables(0)->name, cvm::cv_width-3)
01548 << ")\n";
01549 
01550 os.setf(std::ios::fixed, std::ios::floatfield);
01551 
01552 size_t igrid;
01553 for (igrid = 0; igrid < p.size(); igrid++) {
01554 cvm::real const x_grid = (lower_boundary + (igrid+1)*width);
01555 os << " "
01556 << std::setprecision(cvm::cv_prec)
01557 << std::setw(cvm::cv_width)
01558 << x_grid
01559 << " "
01560 << std::setprecision(cvm::cv_prec)
01561 << std::setw(cvm::cv_width)
01562 << p[igrid] << "\n";
01563 }
01564 proxy->close_output_stream(file_name);
01565 }
01566 return COLVARS_OK;
01567 }
01568 
01569 
01570 std::ostream & colvarbias_restraint_histogram::write_traj_label(std::ostream &os)
01571 {
01572 os << " ";
01573 if (b_output_energy) {
01574 os << " E_"
01575 << cvm::wrap_string(this->name, cvm::en_width-2);
01576 }
01577 return os;
01578 }
01579 
01580 
01581 std::ostream & colvarbias_restraint_histogram::write_traj(std::ostream &os)
01582 {
01583 os << " ";
01584 if (b_output_energy) {
01585 os << " "
01586 << std::setprecision(cvm::en_prec) << std::setw(cvm::en_width)
01587 << bias_energy;
01588 }
01589 return os;
01590 }