diff --git a/lib/colvars/colvarbias_meta.C b/lib/colvars/colvarbias_meta.C
index 71ab221d7..1aac51274 100644
--- a/lib/colvars/colvarbias_meta.C
+++ b/lib/colvars/colvarbias_meta.C
@@ -1,1597 +1,1638 @@
#include <iostream>
#include <sstream>
#include <fstream>
#include <iomanip>
#include <cmath>
#include <algorithm>
// used to set the absolute path of a replica file
#if defined (WIN32) && !defined(__CYGWIN__)
#include <direct.h>
#define CHDIR ::_chdir
#define GETCWD ::_getcwd
#define PATHSEP "\\"
#else
#include <unistd.h>
#define CHDIR ::chdir
#define GETCWD ::getcwd
#define PATHSEP "/"
#endif
#include "colvar.h"
#include "colvarbias_meta.h"
colvarbias_meta::colvarbias_meta()
: colvarbias(),
state_file_step (0),
new_hills_begin (hills.end())
{
}
colvarbias_meta::colvarbias_meta (std::string const &conf, char const *key)
: colvarbias (conf, key),
state_file_step (0),
new_hills_begin (hills.end())
{
if (cvm::n_abf_biases > 0)
cvm::log ("Warning: ABF and metadynamics running at the "
"same time can lead to inconsistent results.\n");
get_keyval (conf, "hillWeight", hill_weight, 0.01);
if (hill_weight == 0.0)
cvm::log ("Warning: hillWeight has been set to zero, "
"this bias will have no effect.\n");
get_keyval (conf, "newHillFrequency", new_hill_freq, 1000);
get_keyval (conf, "hillWidth", hill_width, std::sqrt (2.0 * PI) / 2.0);
get_keyval (conf, "useGrids", use_grids, true);
if (use_grids) {
get_keyval (conf, "gridsUpdateFrequency", grids_freq, new_hill_freq);
get_keyval (conf, "rebinGrids", rebin_grids, false);
expand_grids = false;
for (size_t i = 0; i < colvars.size(); i++) {
if (colvars[i]->expand_boundaries) {
- if (comm == multiple_replicas)
- cvm::fatal_error ("Error: expandBoundaries is not supported when "
- "using metadynamics with multiple replicas.\n");
expand_grids = true;
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": Will expand grids when the colvar \""+
colvars[i]->name+"\" approaches its boundaries.\n");
}
}
get_keyval (conf, "keepHills", keep_hills, false);
get_keyval (conf, "dumpFreeEnergyFile", dump_fes, true);
get_keyval (conf, "saveFreeEnergyFile", dump_fes_save, false);
for (size_t i = 0; i < colvars.size(); i++) {
colvars[i]->enable (colvar::task_grid);
}
hills_energy = new colvar_grid_scalar (colvars);
hills_energy_gradients = new colvar_grid_gradient (colvars);
} else {
rebin_grids = false;
keep_hills = false;
dump_fes = false;
dump_fes_save = false;
+ dump_replica_fes = false;
}
{
bool b_replicas = false;
get_keyval (conf, "multipleReplicas", b_replicas, false);
if (b_replicas)
comm = multiple_replicas;
else
comm = single_replica;
}
if (comm != single_replica) {
+ if (expand_grids)
+ cvm::fatal_error ("Error: expandBoundaries is not supported when "
+ "using more than one replicas; please allocate "
+ "wide enough boundaries for each colvar"
+ "ahead of time.\n");
+
+ if (get_keyval (conf, "dumpPartialFreeEnergyFile", dump_replica_fes, false)) {
+ if (dump_replica_fes && (! dump_fes)) {
+ cvm::log ("Enabling \"dumpFreeEnergyFile\".\n");
+ }
+ }
+
get_keyval (conf, "replicaID", replica_id, std::string (""));
if (!replica_id.size())
cvm::fatal_error ("Error: replicaID must be defined "
"when using more than one replica.\n");
get_keyval (conf, "replicasRegistry",
replicas_registry_file,
(this->name+".replicas.registry.txt"));
get_keyval (conf, "replicaUpdateFrequency",
replica_update_freq, new_hill_freq);
if (keep_hills)
cvm::log ("Warning: in metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": keepHills with more than one replica can lead to a very "
"large amount input/output and slow down your calculations. "
"Please consider disabling it.\n");
{
+ // TODO: one may want to specify the path manually for intricated filesystems?
char *pwd = new char[3001];
if (GETCWD (pwd, 3000) == NULL)
- cvm::fatal_error ("Error: cannot obtain the path to the current working directory.\n");
+ cvm::fatal_error ("Error: cannot get the path of the current working directory.\n");
replica_list_file =
(std::string (pwd)+std::string (PATHSEP)+
this->name+"."+replica_id+".files.txt");
// replica_hills_file and replica_state_file are those written
// by the current replica; within the mirror biases, they are
// those by another replica
replica_hills_file =
(std::string (pwd)+std::string (PATHSEP)+
cvm::output_prefix+".colvars."+this->name+"."+replica_id+".hills");
replica_state_file =
(std::string (pwd)+std::string (PATHSEP)+
cvm::output_prefix+".colvars."+this->name+"."+replica_id+".state");
delete pwd;
+ }
- // now register this replica
-
- // first check that it isn't already there
- bool registered_replica = false;
- std::ifstream reg_is (replicas_registry_file.c_str());
- if (reg_is.good()) { // the file may not be there yet
- std::string existing_replica ("");
- std::string existing_replica_file ("");
- while ((reg_is >> existing_replica) && existing_replica.size() &&
- (reg_is >> existing_replica_file) && existing_replica_file.size()) {
- if (existing_replica == replica_id) {
- // this replica was already registered
- replica_list_file = existing_replica_file;
- reg_is.close();
- registered_replica = true;
- break;
- }
+ // now register this replica
+
+ // first check that it isn't already there
+ bool registered_replica = false;
+ std::ifstream reg_is (replicas_registry_file.c_str());
+ if (reg_is.good()) { // the file may not be there yet
+ std::string existing_replica ("");
+ std::string existing_replica_file ("");
+ while ((reg_is >> existing_replica) && existing_replica.size() &&
+ (reg_is >> existing_replica_file) && existing_replica_file.size()) {
+ if (existing_replica == replica_id) {
+ // this replica was already registered
+ replica_list_file = existing_replica_file;
+ reg_is.close();
+ registered_replica = true;
+ break;
}
- reg_is.close();
}
+ reg_is.close();
+ }
- // if this replica was not included yet, we should generate a
- // new record for it
-
- // but first, write an updated list file; if we're running
- // without grids, use a growing list of hills files, otherwise just one state file and
- std::ofstream list_os (replica_list_file.c_str(),
- (use_grids ? std::ios::trunc : std::ios::app));
- if (! list_os.good())
- cvm::fatal_error ("Error: in opening file \""+
- replica_list_file+"\" for writing.\n");
- list_os << "stateFile " << replica_state_file << "\n";
- list_os << "hillsFile " << replica_hills_file << "\n";
- list_os.close();
-
- // and then add a new record to the registry
- if (! registered_replica) {
- std::ofstream reg_os (replicas_registry_file.c_str(), std::ios::app);
- if (! reg_os.good())
- cvm::fatal_error ("Error: in opening file \""+
- replicas_registry_file+"\" for writing.\n");
- reg_os << replica_id << " " << replica_list_file << "\n";
- reg_os.close();
- }
- }
+ // if this replica was not included yet, we should generate a
+ // new record for it: but first, we write this replica's files,
+ // for the others to read
+ // open the "hills" buffer file
replica_hills_os.open (replica_hills_file.c_str());
if (!replica_hills_os.good())
cvm::fatal_error ("Error: in opening file \""+
replica_hills_file+"\" for writing.\n");
replica_hills_os.setf (std::ios::scientific, std::ios::floatfield);
+
+ // write the state file (so that there is always one available)
+ write_replica_state_file();
+ // schedule to read the state files of the other replicas
+ for (size_t ir = 0; ir < replicas.size(); ir++) {
+ (replicas[ir])->replica_state_file_in_sync = false;
+ }
+
+ // if we're running without grids, use a growing list of "hills" files
+ // otherwise, just one state file and one "hills" file as buffer
+ std::ofstream list_os (replica_list_file.c_str(),
+ (use_grids ? std::ios::trunc : std::ios::app));
+ if (! list_os.good())
+ cvm::fatal_error ("Error: in opening file \""+
+ replica_list_file+"\" for writing.\n");
+ list_os << "stateFile " << replica_state_file << "\n";
+ list_os << "hillsFile " << replica_hills_file << "\n";
+ list_os.close();
+
+ // finally, if add a new record for this replica to the registry
+ if (! registered_replica) {
+ std::ofstream reg_os (replicas_registry_file.c_str(), std::ios::app);
+ if (! reg_os.good())
+ cvm::fatal_error ("Error: in opening file \""+
+ replicas_registry_file+"\" for writing.\n");
+ reg_os << replica_id << " " << replica_list_file << "\n";
+ reg_os.close();
+ }
}
get_keyval (conf, "writeHillsTrajectory", b_hills_traj, false);
if (b_hills_traj) {
std::string const traj_file_name (cvm::output_prefix+
".colvars."+this->name+
( (comm != single_replica) ?
("."+replica_id) :
("") )+
".hills.traj");
hills_traj_os.open (traj_file_name.c_str());
if (!hills_traj_os.good())
cvm::fatal_error ("Error: in opening hills output file \"" +
traj_file_name + "\".\n");
}
if (cvm::debug())
cvm::log ("Done initializing the metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
save_delimiters = false;
}
colvarbias_meta::~colvarbias_meta()
{
if (hills_energy) {
delete hills_energy;
hills_energy = NULL;
}
if (hills_energy_gradients) {
delete hills_energy_gradients;
hills_energy_gradients = NULL;
}
if (replica_hills_os.good())
replica_hills_os.close();
if (hills_traj_os.good())
hills_traj_os.close();
}
// **********************************************************************
// Hill management member functions
// **********************************************************************
std::list<colvarbias_meta::hill>::const_iterator
colvarbias_meta::create_hill (colvarbias_meta::hill const &h)
{
hill_iter const hills_end = hills.end();
hills.push_back (h);
if (new_hills_begin == hills_end) {
// if new_hills_begin is unset, set it for the first time
new_hills_begin = hills.end();
new_hills_begin--;
}
if (use_grids) {
// also add it to the list of hills that are off-grid, which may
// need to be computed analytically when the colvar returns
// off-grid
cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h.centers);
if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
hills_off_grid.push_back (h);
}
}
// output to trajectory (if specified)
if (hills_traj_os.good()) {
hills_traj_os << (hills.back()).output_traj();
if (cvm::debug()) {
hills_traj_os.flush();
}
}
has_data = true;
return hills.end();
}
std::list<colvarbias_meta::hill>::const_iterator
colvarbias_meta::delete_hill (hill_iter &h)
{
if (cvm::debug()) {
cvm::log ("Deleting hill from the metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
", with step number "+
cvm::to_str (h->it)+(h->replica.size() ?
", replica id \""+h->replica :
"")+".\n");
}
if (use_grids && hills_off_grid.size()) {
for (hill_iter hoff = hills_off_grid.begin();
hoff != hills_off_grid.end(); hoff++) {
if (*h == *hoff) {
hills_off_grid.erase (hoff);
break;
}
}
}
if (hills_traj_os.good()) {
// output to the trajectory
hills_traj_os << "# DELETED this hill: "
<< (hills.back()).output_traj()
<< "\n";
if (cvm::debug())
hills_traj_os.flush();
}
return hills.erase (h);
}
cvm::real colvarbias_meta::update()
{
if (cvm::debug())
cvm::log ("Updating the metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
if (use_grids) {
std::vector<int> curr_bin = hills_energy->get_colvars_index();
if (expand_grids) {
// first of all, expand the grids, if specified
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": current coordinates on the grid: "+
cvm::to_str (curr_bin)+".\n");
bool changed_grids = false;
int const min_buffer =
(3 * (size_t) std::floor (hill_width)) + 1;
std::vector<int> new_sizes (hills_energy->sizes());
std::vector<colvarvalue> new_lower_boundaries (hills_energy->lower_boundaries);
std::vector<colvarvalue> new_upper_boundaries (hills_energy->upper_boundaries);
for (size_t i = 0; i < colvars.size(); i++) {
if (! colvars[i]->expand_boundaries)
continue;
cvm::real &new_lb = new_lower_boundaries[i].real_value;
cvm::real &new_ub = new_upper_boundaries[i].real_value;
int &new_size = new_sizes[i];
bool changed_lb = false, changed_ub = false;
if (curr_bin[i] < min_buffer) {
int const extra_points = (min_buffer - curr_bin[i]);
new_lb -= extra_points * colvars[i]->width;
new_size += extra_points;
// changed offset in this direction => the pointer needs to
// be changed, too
curr_bin[i] += extra_points;
changed_lb = true;
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new lower boundary for colvar \""+
colvars[i]->name+"\", at "+
cvm::to_str (new_lower_boundaries[i])+".\n");
}
if (curr_bin[i] > new_size - min_buffer - 1) {
int const extra_points = (curr_bin[i] - (new_size - 1) + min_buffer);
new_ub += extra_points * colvars[i]->width;
new_size += extra_points;
changed_ub = true;
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": new upper boundary for colvar \""+
colvars[i]->name+"\", at "+
cvm::to_str (new_upper_boundaries[i])+".\n");
}
if (changed_lb || changed_ub)
changed_grids = true;
}
if (changed_grids) {
// map everything into new grids
colvar_grid_scalar *new_hills_energy =
new colvar_grid_scalar (*hills_energy);
colvar_grid_gradient *new_hills_energy_gradients =
new colvar_grid_gradient (*hills_energy_gradients);
// supply new boundaries to the new grids
new_hills_energy->lower_boundaries = new_lower_boundaries;
new_hills_energy->upper_boundaries = new_upper_boundaries;
new_hills_energy->create (new_sizes, 0.0, 1);
new_hills_energy_gradients->lower_boundaries = new_lower_boundaries;
new_hills_energy_gradients->upper_boundaries = new_upper_boundaries;
new_hills_energy_gradients->create (new_sizes, 0.0, colvars.size());
new_hills_energy->map_grid (*hills_energy);
new_hills_energy_gradients->map_grid (*hills_energy_gradients);
delete hills_energy;
delete hills_energy_gradients;
hills_energy = new_hills_energy;
hills_energy_gradients = new_hills_energy_gradients;
curr_bin = hills_energy->get_colvars_index();
if (cvm::debug())
cvm::log ("Coordinates on the new grid: "+
cvm::to_str (curr_bin)+".\n");
}
}
}
// add a new hill if the required time interval has passed
if ((cvm::step_absolute() % new_hill_freq) == 0) {
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding a new hill at step "+cvm::to_str (cvm::step_absolute())+".\n");
switch (comm) {
case single_replica:
create_hill (hill (hill_weight, colvars, hill_width));
break;
case multiple_replicas:
create_hill (hill (hill_weight, colvars, hill_width, replica_id));
if (replica_hills_os.good()) {
replica_hills_os << hills.back();
} else {
cvm::fatal_error ("Error: in metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
" while writing hills for the other replicas.\n");
}
break;
}
}
// sync with the other replicas (if needed)
if (comm != single_replica) {
// reread the replicas registry
if ((cvm::step_absolute() % replica_update_freq) == 0) {
update_replicas_registry();
// empty the output buffer
replica_hills_os.flush();
read_replica_files();
}
}
// calculate the biasing energy and forces
bias_energy = 0.0;
for (size_t ir = 0; ir < colvars.size(); ir++) {
colvar_forces[ir].reset();
}
if (comm == multiple_replicas)
for (size_t ir = 0; ir < replicas.size(); ir++) {
replicas[ir]->bias_energy = 0.0;
for (size_t ic = 0; ic < colvars.size(); ic++) {
replicas[ir]->colvar_forces[ic].reset();
}
}
if (use_grids) {
// get the forces from the grid
if ((cvm::step_absolute() % grids_freq) == 0) {
// map the most recent gaussians to the grids
project_hills (new_hills_begin, hills.end(),
hills_energy, hills_energy_gradients);
new_hills_begin = hills.end();
// TODO: we may want to condense all into one replicas array,
// including "this" as the first element
if (comm == multiple_replicas) {
for (size_t ir = 0; ir < replicas.size(); ir++) {
replicas[ir]->project_hills (replicas[ir]->new_hills_begin,
replicas[ir]->hills.end(),
replicas[ir]->hills_energy,
replicas[ir]->hills_energy_gradients);
replicas[ir]->new_hills_begin = replicas[ir]->hills.end();
}
}
}
std::vector<int> curr_bin = hills_energy->get_colvars_index();
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": current coordinates on the grid: "+
cvm::to_str (curr_bin)+".\n");
if (hills_energy->index_ok (curr_bin)) {
// within the grid: add the energy and the forces from there
bias_energy += hills_energy->value (curr_bin);
for (size_t ic = 0; ic < colvars.size(); ic++) {
cvm::real const *f = &(hills_energy_gradients->value (curr_bin));
colvar_forces[ic].real_value += -1.0 * f[ic];
// the gradients are stored, not the forces
}
if (comm == multiple_replicas)
for (size_t ir = 0; ir < replicas.size(); ir++) {
bias_energy += replicas[ir]->hills_energy->value (curr_bin);
cvm::real const *f = &(replicas[ir]->hills_energy_gradients->value (curr_bin));
for (size_t ic = 0; ic < colvars.size(); ic++) {
colvar_forces[ic].real_value += -1.0 * f[ic];
}
}
} else {
// off the grid: compute analytically only the hills at the grid's edges
calc_hills (hills_off_grid.begin(), hills_off_grid.end(), bias_energy);
for (size_t ic = 0; ic < colvars.size(); ic++) {
calc_hills_force (ic, hills_off_grid.begin(), hills_off_grid.end(), colvar_forces);
}
if (comm == multiple_replicas)
for (size_t ir = 0; ir < replicas.size(); ir++) {
calc_hills (replicas[ir]->hills_off_grid.begin(),
replicas[ir]->hills_off_grid.end(),
bias_energy);
for (size_t ic = 0; ic < colvars.size(); ic++) {
calc_hills_force (ic,
replicas[ir]->hills_off_grid.begin(),
replicas[ir]->hills_off_grid.end(),
colvar_forces);
}
}
}
}
// now include the hills that have not been binned yet (starting
// from new_hills_begin)
calc_hills (new_hills_begin, hills.end(), bias_energy);
for (size_t ic = 0; ic < colvars.size(); ic++) {
calc_hills_force (ic, new_hills_begin, hills.end(), colvar_forces);
}
if (cvm::debug())
cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
", hills forces = "+cvm::to_str (colvar_forces)+".\n");
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": adding the forces from the other replicas.\n");
if (comm == multiple_replicas)
for (size_t ir = 0; ir < replicas.size(); ir++) {
calc_hills (replicas[ir]->new_hills_begin,
replicas[ir]->hills.end(),
bias_energy);
for (size_t ic = 0; ic < colvars.size(); ic++) {
calc_hills_force (ic,
replicas[ir]->new_hills_begin,
replicas[ir]->hills.end(),
colvar_forces);
}
if (cvm::debug())
cvm::log ("Hills energy = "+cvm::to_str (bias_energy)+
", hills forces = "+cvm::to_str (colvar_forces)+".\n");
}
return bias_energy;
}
void colvarbias_meta::calc_hills (colvarbias_meta::hill_iter h_first,
colvarbias_meta::hill_iter h_last,
cvm::real &energy,
std::vector<colvarvalue> const &colvar_values)
{
std::vector<colvarvalue> curr_values (colvars.size());
for (size_t i = 0; i < colvars.size(); i++) {
curr_values[i].type (colvars[i]->type());
}
if (colvar_values.size()) {
for (size_t i = 0; i < colvars.size(); i++) {
curr_values[i] = colvar_values[i];
}
} else {
for (size_t i = 0; i < colvars.size(); i++) {
curr_values[i] = colvars[i]->value();
}
}
for (hill_iter h = h_first; h != h_last; h++) {
// compute the gaussian exponent
cvm::real cv_sqdev = 0.0;
for (size_t i = 0; i < colvars.size(); i++) {
colvarvalue const &x = curr_values[i];
colvarvalue const ¢er = h->centers[i];
cvm::real const half_width = 0.5 * h->widths[i];
cv_sqdev += (colvars[i]->dist2 (x, center)) / (half_width*half_width);
}
// compute the gaussian
if (cv_sqdev > 23.0) {
// set it to zero if the exponent is more negative than log(1.0E-05)
h->value (0.0);
} else {
h->value (std::exp (-0.5*cv_sqdev));
}
energy += h->energy();
}
}
void colvarbias_meta::calc_hills_force (size_t const &i,
colvarbias_meta::hill_iter h_first,
colvarbias_meta::hill_iter h_last,
std::vector<colvarvalue> &forces,
std::vector<colvarvalue> const &values)
{
// Retrieve the value of the colvar
colvarvalue x (values.size() ? values[i].type() : colvars[i]->type());
x = (values.size() ? values[i] : colvars[i]->value());
// do the type check only once (all colvarvalues in the hills series
// were already saved with their types matching those in the
// colvars)
switch (colvars[i]->type()) {
case colvarvalue::type_scalar:
for (hill_iter h = h_first; h != h_last; h++) {
if (h->value() == 0.0) continue;
colvarvalue const ¢er = h->centers[i];
cvm::real const half_width = 0.5 * h->widths[i];
forces[i].real_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
(colvars[i]->dist2_lgrad (x, center)).real_value );
}
break;
case colvarvalue::type_vector:
case colvarvalue::type_unitvector:
for (hill_iter h = h_first; h != h_last; h++) {
if (h->value() == 0.0) continue;
colvarvalue const ¢er = h->centers[i];
cvm::real const half_width = 0.5 * h->widths[i];
forces[i].rvector_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
(colvars[i]->dist2_lgrad (x, center)).rvector_value );
}
break;
case colvarvalue::type_quaternion:
for (hill_iter h = h_first; h != h_last; h++) {
if (h->value() == 0.0) continue;
colvarvalue const ¢er = h->centers[i];
cvm::real const half_width = 0.5 * h->widths[i];
forces[i].quaternion_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
(colvars[i]->dist2_lgrad (x, center)).quaternion_value );
}
break;
case colvarvalue::type_notset:
break;
}
}
// **********************************************************************
// grid management functions
// **********************************************************************
void colvarbias_meta::project_hills (colvarbias_meta::hill_iter h_first,
colvarbias_meta::hill_iter h_last,
colvar_grid_scalar *he,
colvar_grid_gradient *hg)
{
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": projecting hills.\n");
// TODO: improve it by looping over a small subgrid instead of the whole grid
std::vector<colvarvalue> colvar_values (colvars.size());
std::vector<cvm::real> colvar_forces_scalar (colvars.size());
std::vector<int> he_ix = he->new_index();
std::vector<int> hg_ix = (hg != NULL) ? hg->new_index() : std::vector<int> (0);
cvm::real hills_energy_here = 0.0;
std::vector<colvarvalue> hills_forces_here (colvars.size(), 0.0);
if (hg != NULL) {
// loop over the points of the grid
for ( ;
(he->index_ok (he_ix)) && (hg->index_ok (hg_ix));
) {
for (size_t i = 0; i < colvars.size(); i++) {
colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
}
// loop over the hills and increment the energy grid locally
hills_energy_here = 0.0;
calc_hills (h_first, h_last, hills_energy_here, colvar_values);
he->acc_value (he_ix, hills_energy_here);
for (size_t i = 0; i < colvars.size(); i++) {
hills_forces_here[i].reset();
calc_hills_force (i, h_first, h_last, hills_forces_here, colvar_values);
colvar_forces_scalar[i] = hills_forces_here[i].real_value;
}
hg->acc_force (hg_ix, &(colvar_forces_scalar.front()));
he->incr (he_ix);
hg->incr (hg_ix);
}
} else {
// simpler version, with just the energy
for ( ; (he->index_ok (he_ix)); ) {
for (size_t i = 0; i < colvars.size(); i++) {
colvar_values[i] = hills_energy->bin_to_value_scalar (he_ix[i], i);
}
hills_energy_here = 0.0;
calc_hills (h_first, h_last, hills_energy_here, colvar_values);
he->acc_value (he_ix, hills_energy_here);
he->incr (he_ix);
}
}
if (! keep_hills) {
hills.erase (hills.begin(), hills.end());
}
}
void colvarbias_meta::recount_hills_off_grid (colvarbias_meta::hill_iter h_first,
colvarbias_meta::hill_iter h_last,
colvar_grid_scalar *he)
{
hills_off_grid.clear();
for (hill_iter h = h_first; h != h_last; h++) {
cvm::real const min_dist = hills_energy->bin_distance_from_boundaries (h->centers);
if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
hills_off_grid.push_back (*h);
}
}
}
// **********************************************************************
// multiple replicas functions
// **********************************************************************
void colvarbias_meta::update_replicas_registry()
{
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
": updating the list of replicas, currently containing "+
cvm::to_str (replicas.size())+" elements.\n");
{
// copy the whole file into a string for convenience
std::string line ("");
std::ifstream reg_file (replicas_registry_file.c_str());
if (reg_file.good()) {
replicas_registry.clear();
while (colvarparse::getline_nocomments (reg_file, line))
replicas_registry.append (line+"\n");
} else {
cvm::fatal_error ("Error: failed to open file \""+replicas_registry_file+
"\" for reading.\n");
}
}
// now parse it
std::istringstream reg_is (replicas_registry);
if (reg_is.good()) {
std::string new_replica ("");
std::string new_replica_file ("");
while ((reg_is >> new_replica) && new_replica.size() &&
(reg_is >> new_replica_file) && new_replica_file.size()) {
if (new_replica == this->replica_id) {
// this is the record for this same replica, skip it
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": skipping this replica's own record: \""+
new_replica+"\", \""+new_replica_file+"\"\n");
new_replica_file.clear();
new_replica.clear();
continue;
}
bool already_loaded = false;
for (size_t ir = 0; ir < replicas.size(); ir++) {
if (new_replica == (replicas[ir])->replica_id) {
// this replica was already added
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": skipping a replica already loaded, \""+
(replicas[ir])->replica_id+"\".\n");
already_loaded = true;
break;
}
}
if (!already_loaded) {
// add this replica to the registry
cvm::log ("Metadynamics bias \""+this->name+"\""+
- ": found a new replica, \""+new_replica+"\".\n");
+ ": accessing replica \""+new_replica+"\".\n");
replicas.push_back (new colvarbias_meta());
(replicas.back())->replica_id = new_replica;
(replicas.back())->replica_list_file = new_replica_file;
(replicas.back())->replica_state_file = "";
(replicas.back())->replica_state_file_in_sync = false;
// Note: the following could become a copy constructor?
(replicas.back())->name = this->name;
(replicas.back())->colvars = colvars;
(replicas.back())->use_grids = use_grids;
(replicas.back())->dump_fes = false;
(replicas.back())->expand_grids = false;
(replicas.back())->rebin_grids = false;
(replicas.back())->keep_hills = false;
(replicas.back())->colvar_forces = colvar_forces;
(replicas.back())->comm = multiple_replicas;
if (use_grids) {
(replicas.back())->hills_energy = new colvar_grid_scalar (colvars);
(replicas.back())->hills_energy_gradients = new colvar_grid_gradient (colvars);
}
}
}
// continue the cycle
new_replica_file = "";
new_replica = "";
} else {
cvm::fatal_error ("Error: cannot read the replicas registry file \""+
replicas_registry+"\".\n");
}
// now (re)read the list file of each replica
for (size_t ir = 0; ir < replicas.size(); ir++) {
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
": reading the list file for replica \""+
(replicas[ir])->replica_id+"\".\n");
std::ifstream list_is ((replicas[ir])->replica_list_file.c_str());
std::string key;
std::string new_state_file, new_hills_file;
if (!(list_is >> key) ||
!(key == std::string ("stateFile")) ||
!(list_is >> new_state_file) ||
!(list_is >> key) ||
!(key == std::string ("hillsFile")) ||
!(list_is >> new_hills_file)) {
cvm::log ("Metadynamics bias \""+this->name+"\""+
": failed to read the file \""+
(replicas[ir])->replica_list_file+"\": will try again after "+
cvm::to_str (replica_update_freq)+" steps.\n");
(replicas[ir])->update_status++;
} else {
(replicas[ir])->update_status = 0;
if (new_state_file != (replicas[ir])->replica_state_file) {
cvm::log ("Metadynamics bias \""+this->name+"\""+
": replica \""+(replicas[ir])->replica_id+
"\" has supplied a new state file, \""+new_state_file+
"\".\n");
(replicas[ir])->replica_state_file_in_sync = false;
(replicas[ir])->replica_state_file = new_state_file;
(replicas[ir])->replica_hills_file = new_hills_file;
}
}
}
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\": the list of replicas contains "+
cvm::to_str (replicas.size())+" elements.\n");
}
void colvarbias_meta::read_replica_files()
{
for (size_t ir = 0; ir < replicas.size(); ir++) {
if (! (replicas[ir])->replica_state_file_in_sync) {
// if a new state file is being read, the hills file is also new
(replicas[ir])->replica_hills_file_pos = 0;
}
// (re)read the state file if necessary
if ( (! (replicas[ir])->has_data) ||
(! (replicas[ir])->replica_state_file_in_sync) ) {
cvm::log ("Metadynamics bias \""+this->name+"\""+
": reading the state of replica \""+
(replicas[ir])->replica_id+"\" from file \""+
(replicas[ir])->replica_state_file+"\".\n");
std::ifstream is ((replicas[ir])->replica_state_file.c_str());
if (! (replicas[ir])->read_restart (is)) {
cvm::log ("Reading from file \""+(replicas[ir])->replica_state_file+
"\" failed or incomplete: will try again in "+
cvm::to_str (replica_update_freq)+" steps.\n");
} else {
// state file has been read successfully
(replicas[ir])->replica_state_file_in_sync = true;
(replicas[ir])->update_status = 0;
}
is.close();
}
// now read the hills added after writing the state file
if ((replicas[ir])->replica_hills_file.size()) {
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
": checking for new hills from replica \""+
(replicas[ir])->replica_id+"\" in the file \""+
(replicas[ir])->replica_hills_file+"\".\n");
// read hills from the other replicas' files; for each file, resume
// the position recorded previously
std::ifstream is ((replicas[ir])->replica_hills_file.c_str());
if (is.good()) {
// try to resume the previous position
is.seekg ((replicas[ir])->replica_hills_file_pos, std::ios::beg);
if (!is.good()){
// if fail (the file may have been overwritten), reset this
// position
is.clear();
is.seekg (0, std::ios::beg);
// reset the counter
(replicas[ir])->replica_hills_file_pos = 0;
// schedule to reread the state file
(replicas[ir])->replica_state_file_in_sync = false;
// and record the failure
(replicas[ir])->update_status++;
cvm::log ("Failed to read the file \""+(replicas[ir])->replica_hills_file+
"\" at the previous position: will try again in "+
cvm::to_str (replica_update_freq)+" steps.\n");
} else {
while ((replicas[ir])->read_hill (is)) {
- if (cvm::debug())
+ // if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
- ": found a new hill, created by replica \""+
+ ": received a hill from replica \""+
(replicas[ir])->replica_id+
"\" at step "+
cvm::to_str (((replicas[ir])->hills.back()).it)+".\n");
}
is.clear();
// store the position for the next read
(replicas[ir])->replica_hills_file_pos = is.tellg();
if (cvm::debug())
cvm::log ("Metadynamics bias \""+this->name+"\""+
": stopped reading file \""+(replicas[ir])->replica_hills_file+
"\" at position "+
cvm::to_str ((replicas[ir])->replica_hills_file_pos)+".\n");
// test whether this is the end of the file
is.seekg (0, std::ios::end);
if (is.tellg() > (replicas[ir])->replica_hills_file_pos+1) {
(replicas[ir])->update_status++;
} else {
(replicas[ir])->update_status = 0;
}
}
} else {
cvm::fatal_error ("Error: cannot read from file \""+
(replicas[ir])->replica_hills_file+"\".\n");
}
is.close();
}
- // size_t const n_flush = (replica_update_freq/new_hill_freq + 1);
- // if ((replicas[ir])->update_status > 3*n_flush) {
- // // TODO: suspend the calculation?
- // cvm::log ("Warning: in metadynamics bias \""+this->name+"\""+
- // ": output files from replica \""+
- // (replicas[ir])->replica_id+
- // "\" have not been updated in about "+
- // cvm::to_str (n_flush*new_hill_freq)+
- // " steps. Please check that its simulation is still running.\n");
- // }
-
+ size_t const n_flush = (replica_update_freq/new_hill_freq + 1);
+ if ((replicas[ir])->update_status > 3*n_flush) {
+ // TODO: suspend the calculation?
+ cvm::log ("Warning: in metadynamics bias \""+this->name+"\""+
+ ": failet do read completely output files from replica \""+
+ (replicas[ir])->replica_id+
+ "\" after more than "+
+ cvm::to_str (n_flush*new_hill_freq)+
+ " steps. Please check that that simulation is still running.\n");
+ }
}
}
// **********************************************************************
// input functions
// **********************************************************************
std::istream & colvarbias_meta::read_restart (std::istream& is)
{
size_t const start_pos = is.tellg();
- if (! has_data)
+ if (comm == single_replica) {
+ // if using a multiple replicas scheme, output messages
+ // are printed before and after calling this function
cvm::log ("Restarting metadynamics bias \""+this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
".\n");
- else
- cvm::log ("Rereading the state of metadynamics bias \""+this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
- ".\n");
-
+ }
std::string key, brace, conf;
if ( !(is >> key) || !(key == "metadynamics") ||
!(is >> brace) || !(brace == "{") ||
!(is >> colvarparse::read_block ("configuration", conf)) ) {
- if (! has_data)
+ if (comm == single_replica)
cvm::log ("Error: in reading restart configuration for metadynamics bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
(replica_state_file_in_sync ? ("at position "+
cvm::to_str (start_pos)+
" in the state file") : "")+".\n");
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
std::string name = "";
if ( colvarparse::get_keyval (conf, "name", name,
std::string (""), colvarparse::parse_silent) &&
(name != this->name) )
cvm::fatal_error ("Error: in the restart file, the "
"\"metadynamics\" block has a different name: different system?\n");
if (name.size() == 0) {
cvm::fatal_error ("Error: \"metadynamics\" block within the restart file "
"has no identifiers.\n");
}
if (comm != single_replica) {
std::string replica = "";
if (colvarparse::get_keyval (conf, "replicaID", replica,
std::string (""), colvarparse::parse_silent) &&
(replica != this->replica_id))
cvm::fatal_error ("Error: in the restart file, the "
"\"metadynamics\" block has a different replicaID: different system?\n");
colvarparse::get_keyval (conf, "step", state_file_step,
cvm::step_absolute(), colvarparse::parse_silent);
}
bool grids_from_restart_file = use_grids;
if (use_grids) {
if (expand_grids) {
// the boundaries of the colvars may have been changed; TODO:
// this reallocation is only for backward-compatibility, and may
// be deleted when grid_parameters (i.e. colvargrid's own
// internal reallocation) has kicked in
delete hills_energy;
delete hills_energy_gradients;
hills_energy = new colvar_grid_scalar (colvars);
hills_energy_gradients = new colvar_grid_gradient (colvars);
}
colvar_grid_scalar *hills_energy_backup = NULL;
colvar_grid_gradient *hills_energy_gradients_backup = NULL;
if (has_data) {
if (cvm::debug())
cvm::log ("Backupping grids for metadynamics bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+".\n");
hills_energy_backup = hills_energy;
hills_energy_gradients_backup = hills_energy_gradients;
hills_energy = new colvar_grid_scalar (colvars);
hills_energy_gradients = new colvar_grid_gradient (colvars);
}
size_t const hills_energy_pos = is.tellg();
if (!(is >> key)) {
if (hills_energy_backup != NULL) {
delete hills_energy;
delete hills_energy_gradients;
hills_energy = hills_energy_backup;
hills_energy_gradients = hills_energy_gradients_backup;
}
is.clear();
is.seekg (hills_energy_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
} else if (!(key == std::string ("hills_energy")) ||
!(hills_energy->read_restart (is))) {
is.clear();
is.seekg (hills_energy_pos, std::ios::beg);
grids_from_restart_file = false;
if (!rebin_grids) {
if (hills_energy_backup == NULL)
cvm::fatal_error ("Error: couldn't read the free energy grid for metadynamics bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
- "; if grids were not enabled in the previous runs, "
- "enable rebinGrids now to regenerate them.\n");
+ "; if useGrids was off when the state file was written, "
+ "enable rebinGrids now to regenerate the grids.\n");
else {
- cvm::log ("Error: couldn't read the free energy grid for metadynamics bias \""+
- this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
+ if (comm == single_replica)
+ cvm::log ("Error: couldn't read the free energy grid for metadynamics bias \""+
+ this->name+"\".\n");
delete hills_energy;
delete hills_energy_gradients;
hills_energy = hills_energy_backup;
hills_energy_gradients = hills_energy_gradients_backup;
is.setstate (std::ios::failbit);
return is;
}
}
}
size_t const hills_energy_gradients_pos = is.tellg();
if (!(is >> key)) {
if (hills_energy_backup != NULL) {
delete hills_energy;
delete hills_energy_gradients;
hills_energy = hills_energy_backup;
hills_energy_gradients = hills_energy_gradients_backup;
}
is.clear();
is.seekg (hills_energy_gradients_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
} else if (!(key == std::string ("hills_energy_gradients")) ||
!(hills_energy_gradients->read_restart (is))) {
is.clear();
is.seekg (hills_energy_gradients_pos, std::ios::beg);
grids_from_restart_file = false;
if (hills_energy_backup == NULL) {
if (!rebin_grids)
cvm::fatal_error ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
- "; if grids were not enabled in the previous runs, "
- "enable rebinGrids now to regenerate them.\n");
+ "; if useGrids was off when the state file was written, "
+ "enable rebinGrids now to regenerate the grids.\n");
else {
- cvm::log ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
- this->name+"\""+
- ((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
+ if (comm == single_replica)
+ cvm::log ("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
+ this->name+"\".\n");
delete hills_energy;
delete hills_energy_gradients;
hills_energy = hills_energy_backup;
hills_energy_gradients = hills_energy_gradients_backup;
is.setstate (std::ios::failbit);
return is;
}
}
}
if (cvm::debug())
cvm::log ("Successfully read new grids for bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
if (hills_energy_backup != NULL) {
// now that we have successfully updated the grids, delete the
// backup copies
if (cvm::debug())
cvm::log ("Deallocating the older grids.\n");
delete hills_energy_backup;
delete hills_energy_gradients_backup;
}
}
bool const existing_hills = (hills.size() > 0);
size_t const old_hills_size = hills.size();
hill_iter old_hills_end = hills.end();
hill_iter old_hills_off_grid_end = hills_off_grid.end();
// read the hills explicitly written (if there are any)
while (read_hill (is)) {
if (cvm::debug())
cvm::log ("Read a previously saved hill under the "
"metadynamics bias \""+
this->name+"\", created at step "+
cvm::to_str ((hills.back()).it)+".\n");
}
is.clear();
new_hills_begin = hills.end();
if (grids_from_restart_file) {
if (hills.size() > old_hills_size)
cvm::log ("Read "+cvm::to_str (hills.size())+
" hills in addition to the grids.\n");
} else {
if (hills.size())
cvm::log ("Read "+cvm::to_str (hills.size())+" hills.\n");
}
if (rebin_grids) {
// allocate new grids (based on the new boundaries and widths just
// read from the configuration file), and project onto them the
// grids just read from the restart file
colvar_grid_scalar *new_hills_energy =
new colvar_grid_scalar (colvars);
colvar_grid_gradient *new_hills_energy_gradients =
new colvar_grid_gradient (colvars);
if (!grids_from_restart_file || (keep_hills && hills.size())) {
// if there are hills, recompute the new grids from them
cvm::log ("Rebinning the energy and forces grids from "+
cvm::to_str (hills.size())+" hills (this may take a while)...\n");
project_hills (hills.begin(), hills.end(),
new_hills_energy, new_hills_energy_gradients);
cvm::log ("rebinning done.\n");
} else {
// otherwise, use the grids in the restart file
cvm::log ("Rebinning the energy and forces grids "
"from the grids in the restart file.\n");
new_hills_energy->map_grid (*hills_energy);
new_hills_energy_gradients->map_grid (*hills_energy_gradients);
}
delete hills_energy;
delete hills_energy_gradients;
hills_energy = new_hills_energy;
hills_energy_gradients = new_hills_energy_gradients;
// assuming that some boundaries have expanded, eliminate those
// off-grid hills that aren't necessary any more
if (hills.size())
recount_hills_off_grid (hills.begin(), hills.end(), hills_energy);
}
if (use_grids) {
if (hills_off_grid.size()) {
cvm::log (cvm::to_str (hills_off_grid.size())+" hills are near the "
"grid boundaries: they will be computed analytically "
"and saved to the state files.\n");
}
}
is >> brace;
if (brace != "}") {
cvm::log ("Incomplete restart information for metadynamics bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
": no closing brace at position "+
cvm::to_str (is.tellg())+" in the file.\n");
is.setstate (std::ios::failbit);
return is;
}
if (cvm::debug())
cvm::log ("colvarbias_meta::read_restart() done\n");
if (existing_hills) {
hills.erase (hills.begin(), old_hills_end);
hills_off_grid.erase (hills_off_grid.begin(), old_hills_off_grid_end);
}
has_data = true;
if (comm != single_replica) {
read_replica_files();
}
return is;
}
std::istream & colvarbias_meta::read_hill (std::istream &is)
{
if (!is) return is; // do nothing if failbit is set
size_t const start_pos = is.tellg();
std::string data;
if ( !(is >> read_block ("hill", data)) ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
size_t h_it;
get_keyval (data, "step", h_it, 0, parse_silent);
if (h_it <= state_file_step) {
if (cvm::debug())
cvm::log ("Skipping a hill older than the state file for metadynamics bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+"\n");
return is;
}
cvm::real h_weight;
get_keyval (data, "weight", h_weight, hill_weight, parse_silent);
std::vector<colvarvalue> h_centers (colvars.size());
for (size_t i = 0; i < colvars.size(); i++) {
h_centers[i].type ((colvars[i]->value()).type());
}
{
// it is safer to read colvarvalue objects one at a time;
// TODO: change this it later
std::string centers_input;
key_lookup (data, "centers", centers_input);
std::istringstream centers_is (centers_input);
for (size_t i = 0; i < colvars.size(); i++) {
centers_is >> h_centers[i];
}
}
std::vector<cvm::real> h_widths (colvars.size());
get_keyval (data, "widths", h_widths,
std::vector<cvm::real> (colvars.size(), (std::sqrt (2.0 * PI) / 2.0)),
parse_silent);
std::string h_replica = "";
if (comm != single_replica) {
get_keyval (data, "replicaID", h_replica, replica_id, parse_silent);
if (h_replica != replica_id)
cvm::fatal_error ("Error: trying to read a hill created by replica \""+h_replica+
"\" for replica \""+replica_id+
"\"; did you swap output files?\n");
}
hill_iter const hills_end = hills.end();
hills.push_back (hill (h_it, h_weight, h_centers, h_widths, h_replica));
if (new_hills_begin == hills_end) {
// if new_hills_begin is unset, set it for the first time
new_hills_begin = hills.end();
new_hills_begin--;
}
if (use_grids) {
// add this also to the list of hills that are off-grid, which will
// be computed analytically
cvm::real const min_dist =
hills_energy->bin_distance_from_boundaries ((hills.back()).centers);
if (min_dist < (3.0 * std::floor (hill_width)) + 1.0) {
hills_off_grid.push_back (hills.back());
}
}
has_data = true;
return is;
}
// **********************************************************************
// output functions
// **********************************************************************
std::ostream & colvarbias_meta::write_restart (std::ostream& os)
{
os << "metadynamics {\n"
<< " configuration {\n"
<< " step " << cvm::step_absolute() << "\n"
<< " name " << this->name << "\n";
if (this->comm != single_replica)
os << " replicaID " << this->replica_id << "\n";
os << " }\n\n";
if (use_grids) {
// this is a very good time to project hills, if you haven't done
// it already!
project_hills (new_hills_begin, hills.end(),
hills_energy, hills_energy_gradients);
new_hills_begin = hills.end();
// write down the grids to the restart file
os << " hills_energy\n";
hills_energy->write_restart (os);
os << " hills_energy_gradients\n";
hills_energy_gradients->write_restart (os);
}
if ( (!use_grids) || keep_hills ) {
// write all hills currently in memory
for (std::list<hill>::const_iterator h = this->hills.begin();
h != this->hills.end();
h++) {
os << *h;
}
} else {
// write just those that are near the grid boundaries
for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
h != this->hills_off_grid.end();
h++) {
os << *h;
}
}
os << "}\n\n";
if (comm != single_replica) {
write_replica_state_file();
// schedule to reread the state files of the other replicas (they
// have also rewritten them)
for (size_t ir = 0; ir < replicas.size(); ir++) {
(replicas[ir])->replica_state_file_in_sync = false;
}
}
if (dump_fes) {
write_pmf();
}
return os;
}
void colvarbias_meta::write_pmf()
{
// allocate a new grid to store the pmf
colvar_grid_scalar *pmf = new colvar_grid_scalar (colvars);
- pmf->add_grid (*hills_energy);
- if (comm != single_replica)
+
+ std::string fes_file_name_prefix (cvm::output_prefix);
+
+ if ((cvm::n_meta_biases > 1) || (cvm::n_abf_biases > 0)) {
+ // if this is not the only free energy integrator, append
+ // this bias's name, to distinguish it from the output of the other
+ // biases producing a .pmf file
+ // TODO: fix for ABF with updateBias == no
+ fes_file_name_prefix += this->name;
+ }
+
+ if ((comm == single_replica) || (dump_replica_fes)) {
+ // output the PMF from this instance or replica
+ pmf->reset();
+ pmf->add_grid (*hills_energy);
+ cvm::real const max = pmf->maximum_value();
+ pmf->add_constant (-1.0 * max);
+ pmf->multiply_constant (-1.0);
+ if (comm != single_replica)
+ fes_file_name_prefix += (".partial");
+ std::string const fes_file_name (fes_file_name_prefix +
+ (dump_fes_save ?
+ "."+cvm::to_str (cvm::step_absolute()) : ""));
+ cvm::backup_file (fes_file_name.c_str());
+ std::ofstream fes_os (fes_file_name.c_str());
+ pmf->write_multicol (fes_os);
+ fes_os.close();
+ }
+
+ if (comm != single_replica) {
+ // output the combined PMF from all replicas
+ pmf->reset();
+ pmf->add_grid (*hills_energy);
for (size_t ir = 0; ir < replicas.size(); ir++) {
pmf->add_grid (*(replicas[ir]->hills_energy));
}
- cvm::real const max = pmf->maximum_value();
- pmf->add_constant (-1.0 * max);
- pmf->multiply_constant (-1.0);
- // if this is the only free energy integrator, the pmf file
- // name is general, otherwise there is a label with the bias
- // name TODO: fix for ABF with updateBias == no
- std::string const fes_file_name =
- ((cvm::n_meta_biases == 1) && (cvm::n_abf_biases == 0)) ?
- std::string (cvm::output_prefix+
- (dump_fes_save ? "."+cvm::to_str (cvm::step_absolute()) : "")+
- ".pmf") :
- std::string (cvm::output_prefix+"."+this->name+
- (dump_fes_save ? "."+cvm::to_str (cvm::step_absolute()) : "")+
- ".pmf");
- cvm::backup_file (fes_file_name.c_str());
- std::ofstream fes_os (fes_file_name.c_str());
- pmf->write_multicol (fes_os);
- fes_os.close();
+ cvm::real const max = pmf->maximum_value();
+ pmf->add_constant (-1.0 * max);
+ pmf->multiply_constant (-1.0);
+ std::string const fes_file_name (fes_file_name_prefix +
+ (dump_fes_save ?
+ "."+cvm::to_str (cvm::step_absolute()) : ""));
+ cvm::backup_file (fes_file_name.c_str());
+ std::ofstream fes_os (fes_file_name.c_str());
+ pmf->write_multicol (fes_os);
+ fes_os.close();
+ }
+
delete pmf;
}
void colvarbias_meta::write_replica_state_file()
{
// write down also the restart for the other replicas: TODO: this
// is duplicated code, that could be cleaned up later
cvm::backup_file (replica_state_file.c_str());
std::ofstream rep_state_os (replica_state_file.c_str());
if (!rep_state_os.good())
cvm::fatal_error ("Error: in opening file \""+
replica_state_file+"\" for writing.\n");
rep_state_os.setf (std::ios::scientific, std::ios::floatfield);
rep_state_os << "\n"
<< "metadynamics {\n"
<< " configuration {\n"
<< " name " << this->name << "\n"
<< " step " << cvm::step_absolute() << "\n";
if (this->comm != single_replica) {
rep_state_os << " replicaID " << this->replica_id << "\n";
}
rep_state_os << " }\n\n";
rep_state_os << " hills_energy\n";
rep_state_os << std::setprecision (cvm::cv_prec)
<< std::setw (cvm::cv_width);
hills_energy->write_restart (rep_state_os);
rep_state_os << " hills_energy_gradients\n";
rep_state_os << std::setprecision (cvm::cv_prec)
<< std::setw (cvm::cv_width);
hills_energy_gradients->write_restart (rep_state_os);
if ( (!use_grids) || keep_hills ) {
// write all hills currently in memory
for (std::list<hill>::const_iterator h = this->hills.begin();
h != this->hills.end();
h++) {
rep_state_os << *h;
}
} else {
// write just those that are near the grid boundaries
for (std::list<hill>::const_iterator h = this->hills_off_grid.begin();
h != this->hills_off_grid.end();
h++) {
rep_state_os << *h;
}
}
rep_state_os << "}\n\n";
rep_state_os.close();
// reopen the hills file
replica_hills_os.close();
replica_hills_os.open (replica_hills_file.c_str());
if (!replica_hills_os.good())
cvm::fatal_error ("Error: in opening file \""+
replica_hills_file+"\" for writing.\n");
replica_hills_os.setf (std::ios::scientific, std::ios::floatfield);
}
std::string colvarbias_meta::hill::output_traj()
{
std::ostringstream os;
os.setf (std::ios::fixed, std::ios::floatfield);
os << std::setw (cvm::it_width) << it << " ";
os.setf (std::ios::scientific, std::ios::floatfield);
os << " ";
for (size_t i = 0; i < centers.size(); i++) {
os << " ";
os << std::setprecision (cvm::cv_prec)
<< std::setw (cvm::cv_width) << centers[i];
}
os << " ";
for (size_t i = 0; i < widths.size(); i++) {
os << " ";
os << std::setprecision (cvm::cv_prec)
<< std::setw (cvm::cv_width) << widths[i];
}
os << " ";
os << std::setprecision (cvm::en_prec)
<< std::setw (cvm::en_width) << W << "\n";
return os.str();
}
std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h)
{
os.setf (std::ios::scientific, std::ios::floatfield);
os << "hill {\n";
os << " step " << std::setw (cvm::it_width) << h.it << "\n";
os << " weight "
<< std::setprecision (cvm::en_prec)
<< std::setw (cvm::en_width)
<< h.W << "\n";
if (h.replica.size())
os << " replicaID " << h.replica << "\n";
os << " centers ";
for (size_t i = 0; i < (h.centers).size(); i++) {
os << " "
<< std::setprecision (cvm::cv_prec)
<< std::setw (cvm::cv_width)
<< h.centers[i];
}
os << "\n";
os << " widths ";
for (size_t i = 0; i < (h.widths).size(); i++) {
os << " "
<< std::setprecision (cvm::cv_prec)
<< std::setw (cvm::cv_width)
<< h.widths[i];
}
os << "\n";
os << "}\n";
return os;
}
diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h
index f19382b00..6e3adb1c5 100644
--- a/lib/colvars/colvarbias_meta.h
+++ b/lib/colvars/colvarbias_meta.h
@@ -1,414 +1,418 @@
#ifndef COLVARBIAS_META_H
#define COLVARBIAS_META_H
#include <vector>
#include <list>
#include <sstream>
#include <fstream>
#include "colvarbias.h"
#include "colvargrid.h"
/// Metadynamics bias (implementation of \link colvarbias \endlink)
class colvarbias_meta : public colvarbias {
public:
/// Communication between different replicas
enum Communication {
/// One replica (default)
single_replica,
/// Hills added concurrently by several replicas
multiple_replicas
};
/// Communication between different replicas
Communication comm;
/// Constructor
colvarbias_meta (std::string const &conf, char const *key);
/// Default constructor
colvarbias_meta();
/// Destructor
virtual ~colvarbias_meta();
virtual cvm::real update();
virtual std::istream & read_restart (std::istream &is);
virtual std::ostream & write_restart (std::ostream &os);
virtual void write_pmf();
class hill;
typedef std::list<hill>::iterator hill_iter;
protected:
/// \brief width of a hill
///
/// The local width of each collective variable, multiplied by this
/// number, provides the hill width along that direction
cvm::real hill_width;
/// \brief Number of simulation steps between two hills
size_t new_hill_freq;
/// Write the hill logfile
bool b_hills_traj;
/// Logfile of hill management (creation and deletion)
std::ofstream hills_traj_os;
/// \brief List of hills used on this bias (total); if a grid is
/// employed, these don't need to be updated at every time step
std::list<hill> hills;
/// \brief Iterator to the first of the "newest" hills (when using
/// grids, those who haven't been mapped yet)
hill_iter new_hills_begin;
/// \brief List of hills used on this bias that are on the boundary
/// edges; these are updated regardless of whether hills are used
std::list<hill> hills_off_grid;
/// \brief Same as new_hills_begin, but for the off-grid ones
hill_iter new_hills_off_grid_begin;
/// Regenerate the hills_off_grid list
void recount_hills_off_grid (hill_iter h_first, hill_iter h_last,
colvar_grid_scalar *ge);
/// Read a hill from a file
std::istream & read_hill (std::istream &is);
/// \brief step present in a state file
///
/// When using grids and reading state files containing them
/// (multiple replicas), this is used to check whether a hill is
/// newer or older than the grids
size_t state_file_step;
/// \brief Add a new hill; if a .hills trajectory is written,
/// write it there; if there is more than one replica, communicate
/// it to the others
virtual std::list<hill>::const_iterator create_hill (hill const &h);
/// \brief Remove a previously saved hill (returns an iterator for
/// the next hill in the list)
virtual std::list<hill>::const_iterator delete_hill (hill_iter &h);
/// \brief Calculate the values of the hills, incrementing
/// bias_energy
virtual void calc_hills (hill_iter h_first,
hill_iter h_last,
cvm::real &energy,
std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0));
/// \brief Calculate the forces acting on the i-th colvar,
/// incrementing colvar_forces[i]; must be called after calc_hills
/// each time the values of the colvars are changed
virtual void calc_hills_force (size_t const &i,
hill_iter h_first,
hill_iter h_last,
std::vector<colvarvalue> &forces,
std::vector<colvarvalue> const &values = std::vector<colvarvalue> (0));
/// Height of new hills
cvm::real hill_weight;
/// \brief Bin the hills on grids of energy and forces, and use them
/// to force the colvars (as opposed to deriving the hills analytically)
bool use_grids;
/// \brief Rebin the hills upon restarting
bool rebin_grids;
/// \brief Should the grids be expanded if necessary?
bool expand_grids;
/// \brief How often the hills should be projected onto the grids
size_t grids_freq;
/// \brief Whether to keep the hills in the restart file (e.g. to do
/// meaningful accurate rebinning afterwards)
bool keep_hills;
/// \brief Dump the free energy surface (.pmf file) every restartFrequency
bool dump_fes;
- /// \brief Keep the free energy surface files at different
- /// iterations, appending a step number to each
+ /// \brief Dump the free energy surface (.pmf file) every restartFrequency
+ /// using only the hills from this replica (only applicable to more than one replica)
+ bool dump_replica_fes;
+
+ /// \brief Dump the free energy surface files at different
+ /// time steps, appending the step number to each file
bool dump_fes_save;
/// \brief Try to read the restart information by allocating new
/// grids before replacing the current ones (used e.g. in
/// multiple_replicas)
bool safely_read_restart;
/// Hill energy, cached on a grid
colvar_grid_scalar *hills_energy;
/// Hill forces, cached on a grid
colvar_grid_gradient *hills_energy_gradients;
/// \brief Project the selected hills onto grids
void project_hills (hill_iter h_first, hill_iter h_last,
colvar_grid_scalar *ge, colvar_grid_gradient *gf);
// Multiple Replicas variables and functions
/// \brief Identifier for this replica
std::string replica_id;
/// \brief File containing the paths to the output files from this replica
std::string replica_file_name;
/// \brief Read the existing replicas on registry
virtual void update_replicas_registry();
/// \brief Read new data from replicas' files
virtual void read_replica_files();
/// \brief Write data to other replicas
virtual void write_replica_state_file();
/// \brief Additional, "mirror" metadynamics biases, to collect info
/// from the other replicas
///
/// These are supposed to be synchronized by reading data from the
/// other replicas, and not be modified by the "local" replica
std::vector<colvarbias_meta *> replicas;
/// \brief Frequency at which data the "mirror" biases are updated
size_t replica_update_freq;
/// List of replicas (and their output list files): contents are
/// copied into replicas_registry for convenience
std::string replicas_registry_file;
/// List of replicas (and their output list files)
std::string replicas_registry;
/// List of files written by this replica
std::string replica_list_file;
/// Hills energy and gradients written specifically for other
/// replica (in addition to its own restart file)
std::string replica_state_file;
/// Whether a mirror bias has read the latest version of its state file
bool replica_state_file_in_sync;
/// If there was a failure reading one of the files (because they
/// are not complete), this counter is incremented
size_t update_status;
/// Explicit hills communicated between replicas
///
/// This file becomes empty after replica_state_file is rewritten
std::string replica_hills_file;
/// \brief Output stream corresponding to replica_hills_file
std::ofstream replica_hills_os;
/// Position within replica_hills_file (when reading it)
size_t replica_hills_file_pos;
};
/// \brief A hill for the metadynamics bias
class colvarbias_meta::hill {
protected:
/// Value of the hill function (ranges between 0 and 1)
cvm::real hill_value;
/// Scale factor, which could be modified at runtime (default: 1)
cvm::real sW;
/// Maximum height in energy of the hill
cvm::real W;
/// Center of the hill in the collective variable space
std::vector<colvarvalue> centers;
/// Widths of the hill in the collective variable space
std::vector<cvm::real> widths;
public:
friend class colvarbias_meta;
/// Time step at which this hill was added
size_t it;
/// Identity of the replica who added this hill (only in multiple replica simulations)
std::string replica;
/// \brief Runtime constructor: data are read directly from
/// collective variables \param weight Weight of the hill \param
/// cv Pointer to the array of collective variables involved \param
/// replica (optional) Identity of the replica which creates the
/// hill
inline hill (cvm::real const &W_in,
std::vector<colvar *> &cv,
cvm::real const &hill_width,
std::string const &replica_in = "")
: sW (1.0),
W (W_in),
centers (cv.size()),
widths (cv.size()),
it (cvm::it),
replica (replica_in)
{
for (size_t i = 0; i < cv.size(); i++) {
centers[i].type (cv[i]->type());
centers[i] = cv[i]->value();
widths[i] = cv[i]->width * hill_width;
}
if (cvm::debug())
cvm::log ("New hill, applied to "+cvm::to_str (cv.size())+
" collective variables, with centers "+
cvm::to_str (centers)+", widths "+
cvm::to_str (widths)+" and weight "+
cvm::to_str (W)+".\n");
}
/// \brief General constructor: all data are explicitly passed as
/// arguments (used for instance when reading hills saved on a
/// file) \param it Time step of creation of the hill \param
/// weight Weight of the hill \param centers Center of the hill
/// \param widths Width of the hill around centers \param replica
/// (optional) Identity of the replica which creates the hill
inline hill (size_t const &it_in,
cvm::real const &W_in,
std::vector<colvarvalue> const ¢ers_in,
std::vector<cvm::real> const &widths_in,
std::string const &replica_in = "")
: sW (1.0),
W (W_in),
centers (centers_in),
widths (widths_in),
it (it_in),
replica (replica_in)
{}
/// Copy constructor
inline hill (colvarbias_meta::hill const &h)
: sW (1.0),
W (h.W),
centers (h.centers),
widths (h.widths),
it (h.it),
replica (h.replica)
{}
/// Destructor
inline ~hill()
{}
/// Get the energy
inline cvm::real energy()
{
return W * sW * hill_value;
}
/// Get the energy using another hill weight
inline cvm::real energy (cvm::real const &new_weight)
{
return new_weight * sW * hill_value;
}
/// Get the current hill value
inline cvm::real const &value()
{
return hill_value;
}
/// Set the hill value as specified
inline void value (cvm::real const &new_value)
{
hill_value = new_value;
}
/// Get the weight
inline cvm::real weight()
{
return W * sW;
}
/// Scale the weight with this factor (by default 1.0 is used)
inline void scale (cvm::real const &new_scale_fac)
{
sW = new_scale_fac;
}
/// Get the center of the hill
inline std::vector<colvarvalue> & center()
{
return centers;
}
/// Get the i-th component of the center
inline colvarvalue & center (size_t const &i)
{
return centers[i];
}
/// Comparison operator
inline friend bool operator < (hill const &h1, hill const &h2)
{
if (h1.it < h2.it) return true;
else return false;
}
/// Comparison operator
inline friend bool operator <= (hill const &h1, hill const &h2)
{
if (h1.it <= h2.it) return true;
else return false;
}
/// Comparison operator
inline friend bool operator > (hill const &h1, hill const &h2)
{
if (h1.it > h2.it) return true;
else return false;
}
/// Comparison operator
inline friend bool operator >= (hill const &h1, hill const &h2)
{
if (h1.it >= h2.it) return true;
else return false;
}
/// Comparison operator
inline friend bool operator == (hill const &h1, hill const &h2)
{
if ( (h1.it >= h2.it) && (h1.replica == h2.replica) ) return true;
else return false;
}
/// Represent the hill ina string suitable for a trajectory file
std::string output_traj();
/// Write the hill to an output stream
inline friend std::ostream & operator << (std::ostream &os,
hill const &h);
};
#endif
// Emacs
// Local Variables:
// mode: C++
// End:
diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h
index 7bbf45d7f..f5171fa01 100644
--- a/lib/colvars/colvargrid.h
+++ b/lib/colvars/colvargrid.h
@@ -1,1157 +1,1163 @@
// -*- c++ -*-
#ifndef COLVARGRID_H
#define COLVARGRID_H
#include <iostream>
#include <iomanip>
#include <cmath>
#include "colvar.h"
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
/// \brief Grid of values of a function of several collective
/// variables \param T The data type
///
/// Only scalar colvars supported so far
template <class T> class colvar_grid : public colvarparse {
protected:
/// Number of dimensions
size_t nd;
/// Number of points along each dimension
std::vector<int> nx;
/// Cumulative number of points along each dimension
std::vector<int> nxc;
/// \brief Multiplicity of each datum (allow the binning of
/// non-scalar types)
size_t mult;
/// Total number of grid points
size_t nt;
/// Low-level array of values
std::vector<T> data;
/// Newly read data (used for count grids, when adding several grids read from disk)
std::vector<size_t> new_data;
/// Colvars collected in this grid
std::vector<colvar *> cv;
/// Get the low-level index corresponding to an index
inline size_t address (std::vector<int> const &ix) const
{
size_t addr = 0;
for (size_t i = 0; i < nd; i++) {
addr += ix[i]*nxc[i];
if (cvm::debug()) {
if (ix[i] >= nx[i])
cvm::fatal_error ("Error: exceeding bounds in colvar_grid.\n");
}
}
return addr;
}
public:
/// Lower boundaries of the colvars in this grid
std::vector<colvarvalue> lower_boundaries;
/// Upper boundaries of the colvars in this grid
std::vector<colvarvalue> upper_boundaries;
/// Whether some colvars are periodic
std::vector<bool> periodic;
/// Widths of the colvars in this grid
std::vector<cvm::real> widths;
/// True if this is a count grid related to another grid of data
bool has_parent_data;
/// Whether this grid has been filled with data or is still empty
bool has_data;
/// Return the number of colvars
inline size_t number_of_colvars() const
{
return nd;
}
/// Return the number of points in the i-th direction, if provided, or
/// the total number
inline size_t number_of_points (int const icv = -1) const
{
if (icv < 0) {
return nt;
} else {
return nx[icv];
}
}
/// Get the sizes in each direction
inline std::vector<int> const & sizes() const
{
return nx;
}
/// Set the sizes in each direction
inline void set_sizes (std::vector<int> const &new_sizes)
{
nx = new_sizes;
}
/// Return the multiplicity of the type used
inline size_t multiplicity() const
{
return mult;
}
/// \brief Allocate data (allow initialization also after construction)
void create (std::vector<int> const &nx_i,
T const &t = T(),
size_t const &mult_i = 1)
{
mult = mult_i;
nd = nx_i.size();
nxc.resize (nd);
nx = nx_i;
nt = mult;
for (int i = nd-1; i >= 0; i--) {
if (nx_i[i] <= 0)
cvm::fatal_error ("Error: providing an invalid number of points, "+
cvm::to_str (nx_i[i])+".\n");
nxc[i] = nt;
nt *= nx[i];
}
data.reserve (nt);
data.assign (nt, t);
}
/// \brief Allocate data (allow initialization also after construction)
void create()
{
create (this->nx, T(), this->mult);
}
+ /// \brief Reset data (in case the grid is being reused)
+ void reset (T const &t = T())
+ {
+ data.assign (nt, t);
+ }
+
/// Default constructor
colvar_grid() : has_data (false)
{
save_delimiters = false;
nd = nt = 0;
}
/// Destructor
virtual ~colvar_grid()
{}
/// \brief "Almost copy-constructor": only copies configuration
/// parameters from another grid, but doesn't reallocate stuff;
/// create() must be called after that;
colvar_grid (colvar_grid<T> const &g) : has_data (false),
nd (g.nd),
mult (g.mult),
cv (g.cv),
lower_boundaries (g.lower_boundaries),
upper_boundaries (g.upper_boundaries),
periodic (g.periodic),
widths (g.widths),
data()
{
save_delimiters = false;
}
/// \brief Constructor from explicit grid sizes \param nx_i Number
/// of grid points along each dimension \param t Initial value for
/// the function at each point (optional) \param mult_i Multiplicity
/// of each value
colvar_grid (std::vector<int> const &nx_i,
T const &t = T(),
size_t const &mult_i = 1) : has_data (false)
{
save_delimiters = false;
this->create (nx_i, t, mult_i);
}
/// \brief Constructor from a vector of colvars
colvar_grid (std::vector<colvar *> const &colvars,
T const &t = T(),
size_t const &mult_i = 1,
bool margin = false)
: cv (colvars), has_data (false)
{
save_delimiters = false;
std::vector<int> nx_i;
if (cvm::debug())
cvm::log ("Allocating a grid for "+cvm::to_str (colvars.size())+
" collective variables.\n");
for (size_t i = 0; i < cv.size(); i++) {
if (cv[i]->type() != colvarvalue::type_scalar) {
cvm::fatal_error ("Colvar grids can only be automatically "
"constructed for scalar variables. "
"ABF and histogram can not be used; metadynamics "
"can be used with useGrids disabled.\n");
}
if (cv[i]->width <= 0.0) {
cvm::fatal_error ("Tried to initialize a grid on a "
"variable with negative or zero width.\n");
}
if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) {
cvm::fatal_error ("Tried to initialize a grid on a "
"variable with undefined boundaries.\n");
}
widths.push_back (cv[i]->width);
periodic.push_back (cv[i]->periodic_boundaries());
if (margin) {
if (periodic[i]) {
// Shift the grid by half the bin width (values at edges instead of center of bins)
lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
upper_boundaries.push_back (cv[i]->upper_boundary.real_value - 0.5 * widths[i]);
} else {
// Make this grid larger by one bin width
lower_boundaries.push_back (cv[i]->lower_boundary.real_value - 0.5 * widths[i]);
upper_boundaries.push_back (cv[i]->upper_boundary.real_value + 0.5 * widths[i]);
}
} else {
lower_boundaries.push_back (cv[i]->lower_boundary);
upper_boundaries.push_back (cv[i]->upper_boundary);
}
{
cvm::real nbins = ( upper_boundaries[i].real_value -
lower_boundaries[i].real_value ) / widths[i];
int nbins_round = (int)(nbins+0.5);
if (std::fabs (nbins - cvm::real (nbins_round)) > 1.0E-10) {
cvm::log ("Warning: grid interval ("+
cvm::to_str (lower_boundaries[i], cvm::cv_width, cvm::cv_prec)+" - "+
cvm::to_str (upper_boundaries[i], cvm::cv_width, cvm::cv_prec)+
") is not commensurate to its bin width ("+
cvm::to_str (widths[i], cvm::cv_width, cvm::cv_prec)+").\n");
upper_boundaries[i].real_value = lower_boundaries[i].real_value +
(nbins_round * widths[i]);
}
if (cvm::debug())
cvm::log ("Number of points is "+cvm::to_str ((int) nbins_round)+
" for the colvar no. "+cvm::to_str (i+1)+".\n");
nx_i.push_back (nbins_round);
}
}
create (nx_i, t, mult_i);
}
/// Wrap an index vector around periodic boundary conditions
/// also checks validity of non-periodic indices
inline void wrap (std::vector<int> & ix) const
{
for (size_t i = 0; i < nd; i++) {
if (periodic[i]) {
ix[i] = (ix[i] + nx[i]) % nx[i]; //to ensure non-negative result
} else {
if (ix[i] < 0 || ix[i] >= nx[i])
cvm::fatal_error ("Trying to wrap illegal index vector (non-PBC): "
+ cvm::to_str (ix));
}
}
}
/// \brief Report the bin corresponding to the current value of variable i
inline int current_bin_scalar(int const i) const
{
return value_to_bin_scalar (cv[i]->value(), i);
}
/// \brief Use the lower boundary and the width to report which bin
/// the provided value is in
inline int value_to_bin_scalar (colvarvalue const &value, const int i) const
{
return (int) std::floor ( (value.real_value - lower_boundaries[i].real_value) / widths[i] );
}
/// \brief Same as the standard version, but uses another grid definition
inline int value_to_bin_scalar (colvarvalue const &value,
colvarvalue const &new_offset,
cvm::real const &new_width) const
{
return (int) std::floor ( (value.real_value - new_offset.real_value) / new_width );
}
/// \brief Use the two boundaries and the width to report the
/// central value corresponding to a bin index
inline colvarvalue bin_to_value_scalar (int const &i_bin, int const i) const
{
return lower_boundaries[i].real_value + widths[i] * (0.5 + i_bin);
}
/// \brief Same as the standard version, but uses different parameters
inline colvarvalue bin_to_value_scalar (int const &i_bin,
colvarvalue const &new_offset,
cvm::real const &new_width) const
{
return new_offset.real_value + new_width * (0.5 + i_bin);
}
/// Set the value at the point with index ix
inline void set_value (std::vector<int> const &ix,
T const &t,
size_t const &imult = 0)
{
data[this->address (ix)+imult] = t;
has_data = true;
}
/// \brief Get the binned value indexed by ix, or the first of them
/// if the multiplicity is larger than 1
inline T const & value (std::vector<int> const &ix,
size_t const &imult = 0) const
{
return data[this->address (ix) + imult];
}
/// \brief Add a constant to all elements (fast loop)
inline void add_constant (T const &t)
{
for (size_t i = 0; i < nt; i++)
data[i] += t;
has_data = true;
}
/// \brief Multiply all elements by a scalar constant (fast loop)
inline void multiply_constant (cvm::real const &a)
{
for (size_t i = 0; i < nt; i++)
data[i] *= a;
}
/// \brief Get the bin indices corresponding to the provided values of
/// the colvars
inline std::vector<int> const get_colvars_index (std::vector<colvarvalue> const &values) const
{
std::vector<int> index = new_index();
for (size_t i = 0; i < nd; i++) {
index[i] = value_to_bin_scalar (values[i], i);
}
return index;
}
/// \brief Get the bin indices corresponding to the current values
/// of the colvars
inline std::vector<int> const get_colvars_index() const
{
std::vector<int> index = new_index();
for (size_t i = 0; i < nd; i++) {
index[i] = current_bin_scalar (i);
}
return index;
}
/// \brief Get the minimal distance (in number of bins) from the
/// boundaries; a negative number is returned if the given point is
/// off-grid
inline cvm::real bin_distance_from_boundaries (std::vector<colvarvalue> const &values)
{
cvm::real minimum = 1.0E+16;
for (size_t i = 0; i < nd; i++) {
if (periodic[i]) continue;
cvm::real dl = std::sqrt (cv[i]->dist2 (values[i], lower_boundaries[i])) / widths[i];
cvm::real du = std::sqrt (cv[i]->dist2 (values[i], upper_boundaries[i])) / widths[i];
if (values[i].real_value < lower_boundaries[i])
dl *= -1.0;
if (values[i].real_value > upper_boundaries[i])
du *= -1.0;
if (dl < minimum)
minimum = dl;
if (du < minimum)
minimum = du;
}
return minimum;
}
/// \brief Add data from another grid of the same type
///
/// Note: this function maps other_grid inside this one regardless
/// of whether it fits or not.
void map_grid (colvar_grid<T> const &other_grid)
{
if (other_grid.multiplicity() != this->multiplicity())
cvm::fatal_error ("Error: trying to merge two grids with values of "
"different multiplicity.\n");
std::vector<colvarvalue> const &gb = this->lower_boundaries;
std::vector<cvm::real> const &gw = this->widths;
std::vector<colvarvalue> const &ogb = other_grid.lower_boundaries;
std::vector<cvm::real> const &ogw = other_grid.widths;
std::vector<int> ix = this->new_index();
std::vector<int> oix = other_grid.new_index();
if (cvm::debug())
cvm::log ("Remapping grid...\n");
for ( ; this->index_ok (ix); this->incr (ix)) {
for (size_t i = 0; i < nd; i++) {
oix[i] =
value_to_bin_scalar (bin_to_value_scalar (ix[i], gb[i], gw[i]),
ogb[i],
ogw[i]);
}
if (! other_grid.index_ok (oix)) {
continue;
}
for (size_t im = 0; im < mult; im++) {
this->set_value (ix, other_grid.value (oix, im), im);
}
}
has_data = true;
if (cvm::debug())
cvm::log ("Remapping done.\n");
}
/// \brief Add data from another grid of the same type, AND
/// identical definition (boundaries, widths)
void add_grid (colvar_grid<T> const &other_grid,
cvm::real scale_factor = 1.0)
{
if (other_grid.multiplicity() != this->multiplicity())
cvm::fatal_error ("Error: trying to sum togetehr two grids with values of "
"different multiplicity.\n");
if (scale_factor != 1.0)
for (size_t i = 0; i < data.size(); i++) {
data[i] += scale_factor * other_grid.data[i];
}
else
// skip multiplication if possible
for (size_t i = 0; i < data.size(); i++) {
data[i] += other_grid.data[i];
}
has_data = true;
}
/// \brief Return the value suitable for output purposes (so that it
/// may be rescaled or manipulated without changing it permanently)
virtual inline T value_output (std::vector<int> const &ix,
size_t const &imult = 0)
{
return value (ix, imult);
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (the two may be different,
/// e.g. when using colvar_grid_count)
virtual inline void value_input (std::vector<int> const &ix,
T const &t,
size_t const &imult = 0,
bool add = false)
{
if ( add )
data[address (ix) + imult] += t;
else
data[address (ix) + imult] = t;
has_data = true;
}
// /// Get the pointer to the binned value indexed by ix
// inline T const *value_p (std::vector<int> const &ix)
// {
// return &(data[address (ix)]);
// }
/// \brief Get the index corresponding to the "first" bin, to be
/// used as the initial value for an index in looping
inline std::vector<int> const new_index() const
{
return std::vector<int> (nd, 0);
}
/// \brief Check that the index is within range in each of the
/// dimensions
inline bool index_ok (std::vector<int> const &ix) const
{
for (size_t i = 0; i < nd; i++) {
if ( (ix[i] < 0) || (ix[i] >= int (nx[i])) )
return false;
}
return true;
}
/// \brief Increment the index, in a way that will make it loop over
/// the whole nd-dimensional array
inline void incr (std::vector<int> &ix) const
{
for (int i = ix.size()-1; i >= 0; i--) {
ix[i]++;
if (ix[i] >= nx[i]) {
if (i > 0) {
ix[i] = 0;
continue;
} else {
// this is the last iteration, a non-valid index is being
// set for the outer index, which will be caught by
// index_ok()
ix[0] = nx[0];
return;
}
} else {
return;
}
}
}
/// \brief Write the grid parameters (number of colvars, boundaries, width and number of points)
std::ostream & write_params (std::ostream &os)
{
os << "grid_parameters {\n n_colvars " << nd << "\n";
os << " lower_boundaries ";
for (size_t i = 0; i < nd; i++)
os << " " << lower_boundaries[i];
os << "\n";
os << " upper_boundaries ";
for (size_t i = 0; i < nd; i++)
os << " " << upper_boundaries[i];
os << "\n";
os << " widths ";
for (size_t i = 0; i < nd; i++)
os << " " << widths[i];
os << "\n";
os << " sizes ";
for (size_t i = 0; i < nd; i++)
os << " " << nx[i];
os << "\n";
os << "}\n";
return os;
}
bool parse_params (std::string const &conf)
{
std::vector<int> old_nx = nx;
std::vector<colvarvalue> old_lb = lower_boundaries;
{
size_t nd_in = 0;
colvarparse::get_keyval (conf, "n_colvars", nd_in, nd, colvarparse::parse_silent);
if (nd_in != nd)
cvm::fatal_error ("Error: trying to read data for a grid "
"that contains a different number of colvars ("+
cvm::to_str (nd_in)+") than the grid defined "
"in the configuration file ("+cvm::to_str (nd)+
").\n");
}
colvarparse::get_keyval (conf, "lower_boundaries",
lower_boundaries, lower_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval (conf, "upper_boundaries",
upper_boundaries, upper_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval (conf, "widths", widths, widths, colvarparse::parse_silent);
colvarparse::get_keyval (conf, "sizes", nx, nx, colvarparse::parse_silent);
bool new_params = false;
for (size_t i = 0; i < nd; i++) {
if ( (old_nx[i] != nx[i]) ||
(std::sqrt (cv[i]->dist2 (old_lb[i],
lower_boundaries[i])) > 1.0E-10) ) {
new_params = true;
}
}
// reallocate the array in case the grid params have just changed
if (new_params) {
data.resize (0);
this->create (nx, T(), mult);
}
return true;
}
/// \brief Check that the grid information inside (boundaries,
/// widths, ...) is consistent with the current setting of the
/// colvars
void check_consistency()
{
for (size_t i = 0; i < nd; i++) {
if ( (std::sqrt (cv[i]->dist2 (cv[i]->lower_boundary,
lower_boundaries[i])) > 1.0E-10) ||
(std::sqrt (cv[i]->dist2 (cv[i]->upper_boundary,
upper_boundaries[i])) > 1.0E-10) ||
(std::sqrt (cv[i]->dist2 (cv[i]->width,
widths[i])) > 1.0E-10) ) {
cvm::fatal_error ("Error: restart information for a grid is "
"inconsistent with that of its colvars.\n");
}
}
}
/// \brief Check that the grid information inside (boundaries,
/// widths, ...) is consistent with the one of another grid
void check_consistency (colvar_grid<T> const &other_grid)
{
for (size_t i = 0; i < nd; i++) {
// we skip dist2(), because periodicities and the like should
// matter: boundaries should be EXACTLY the same (otherwise,
// map_grid() should be used)
if ( (std::fabs (other_grid.lower_boundaries[i] -
lower_boundaries[i]) > 1.0E-10) ||
(std::fabs (other_grid.upper_boundaries[i] -
upper_boundaries[i]) > 1.0E-10) ||
(std::fabs (other_grid.widths[i] -
widths[i]) > 1.0E-10) ||
(data.size() != other_grid.data.size()) ) {
cvm::fatal_error ("Error: inconsistency between "
"two grids that are supposed to be equal, "
"aside from the data stored.\n");
}
}
}
/// \brief Write the grid data without labels, as they are
/// represented in memory
/// \param buf_size Number of values per line
std::ostream & write_raw (std::ostream &os,
size_t const buf_size = 3)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
std::vector<int> ix = new_index();
size_t count = 0;
for ( ; index_ok (ix); incr (ix)) {
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw (w) << std::setprecision (p)
<< value_output (ix, imult);
if (((++count) % buf_size) == 0)
os << "\n";
}
}
// write a final newline only if buffer is not empty
if ((count % buf_size) != 0)
os << "\n";
return os;
}
/// \brief Read data written by colvar_grid::write_raw()
std::istream & read_raw (std::istream &is)
{
size_t const start_pos = is.tellg();
for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix)) {
for (size_t imult = 0; imult < mult; imult++) {
T new_value;
if (is >> new_value) {
value_input (ix, new_value, imult);
} else {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
}
}
has_data = true;
return is;
}
/// \brief To be called after colvar_grid::read_raw() returns an error
void read_raw_error()
{
cvm::fatal_error ("Error: failed to read all of the grid points from file. Possible explanations: grid parameters in the configuration (lowerBoundary, upperBoundary, width) are different from those in the file, or the file is corrupt/incomplete.\n");
}
/// \brief Write the grid in a format which is both human readable
/// and suitable for visualization e.g. with gnuplot
void write_multicol (std::ostream &os)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
// Data in the header: nColvars, then for each
// xiMin, dXi, nPoints, periodic
os << std::setw (2) << "# " << nd << "\n";
for (size_t i = 0; i < nd; i++) {
os << "# "
<< std::setw (10) << lower_boundaries[i]
<< std::setw (10) << widths[i]
<< std::setw (10) << nx[i] << " "
<< periodic[i] << "\n";
}
for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
if (ix.back() == 0) {
// if the last index is 0, add a new line to mark the new record
os << "\n";
}
for (size_t i = 0; i < nd; i++) {
os << " "
<< std::setw (w) << std::setprecision (p)
<< bin_to_value_scalar (ix[i], i);
}
os << " ";
for (size_t imult = 0; imult < mult; imult++) {
os << " "
<< std::setw (w) << std::setprecision (p)
<< value_output (ix, imult);
}
os << "\n";
}
}
/// \brief Read a grid written by colvar_grid::write_multicol()
/// Adding data if add is true, replacing if false
std::istream & read_multicol (std::istream &is, bool add = false)
{
// Data in the header: nColvars, then for each
// xiMin, dXi, nPoints, periodic
std::string hash;
cvm::real lower, width, x;
size_t n, periodic;
bool remap;
std::vector<T> new_value;
std::vector<int> nx_read;
std::vector<int> bin;
if ( cv.size() != nd ) {
cvm::fatal_error ("Cannot read grid file: missing reference to colvars.");
}
if ( !(is >> hash) || (hash != "#") ) {
cvm::fatal_error ("Error reading grid at position "+
cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
}
is >> n;
if ( n != nd ) {
cvm::fatal_error ("Error reading grid: wrong number of collective variables.\n");
}
nx_read.resize (n);
bin.resize (n);
new_value.resize (mult);
if (this->has_parent_data && add) {
new_data.resize (data.size());
}
remap = false;
for (size_t i = 0; i < nd; i++ ) {
if ( !(is >> hash) || (hash != "#") ) {
cvm::fatal_error ("Error reading grid at position "+
cvm::to_str (is.tellg())+" in stream (read \"" + hash + "\")\n");
}
is >> lower >> width >> nx_read[i] >> periodic;
if ( (std::fabs (lower - lower_boundaries[i].real_value) > 1.0e-10) ||
(std::fabs (width - widths[i] ) > 1.0e-10) ||
(nx_read[i] != nx[i]) ) {
cvm::log ("Warning: reading from different grid definition (colvar "
+ cvm::to_str (i+1) + "); remapping data on new grid.\n");
remap = true;
}
}
if ( remap ) {
// re-grid data
while (is.good()) {
bool end_of_file = false;
for (size_t i = 0; i < nd; i++ ) {
if ( !(is >> x) ) end_of_file = true;
bin[i] = value_to_bin_scalar (x, i);
}
if (end_of_file) break;
for (size_t imult = 0; imult < mult; imult++) {
is >> new_value[imult];
}
if ( index_ok(bin) ) {
for (size_t imult = 0; imult < mult; imult++) {
value_input (bin, new_value[imult], imult, add);
}
}
}
} else {
// do not re-grid the data but assume the same grid is used
for (std::vector<int> ix = new_index(); index_ok (ix); incr (ix) ) {
for (size_t i = 0; i < nd; i++ ) {
is >> x;
}
for (size_t imult = 0; imult < mult; imult++) {
is >> new_value[imult];
value_input (ix, new_value[imult], imult, add);
}
}
}
has_data = true;
return is;
}
};
/// \brief Colvar_grid derived class to hold counters in discrete
/// n-dim colvar space
class colvar_grid_count : public colvar_grid<size_t>
{
public:
/// Default constructor
colvar_grid_count();
/// Destructor
virtual inline ~colvar_grid_count()
{}
/// Constructor
colvar_grid_count (std::vector<int> const &nx_i,
size_t const &def_count = 0);
/// Constructor from a vector of colvars
colvar_grid_count (std::vector<colvar *> &colvars,
size_t const &def_count = 0);
/// Increment the counter at given position
inline void incr_count (std::vector<int> const &ix)
{
++(data[this->address (ix)]);
}
/// \brief Get the binned count indexed by ix from the newly read data
inline size_t const & new_count (std::vector<int> const &ix,
size_t const &imult = 0)
{
return new_data[address (ix) + imult];
}
/// \brief Read the grid from a restart
std::istream & read_restart (std::istream &is);
/// \brief Write the grid to a restart
std::ostream & write_restart (std::ostream &os);
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual inline void value_input (std::vector<int> const &ix,
size_t const &t,
size_t const &imult = 0,
bool add = false)
{
if (add) {
data[address (ix)] += t;
if (this->has_parent_data) {
// save newly read data for inputting parent grid
new_data[address (ix)] = t;
}
} else {
data[address (ix)] = t;
}
has_data = true;
}
};
/// Class for accumulating a scalar function on a grid
class colvar_grid_scalar : public colvar_grid<cvm::real>
{
public:
/// \brief Provide the associated sample count by which each binned value
/// should be divided
colvar_grid_count *samples;
/// Default constructor
colvar_grid_scalar();
/// Copy constructor (needed because of the grad pointer)
colvar_grid_scalar (colvar_grid_scalar const &g);
/// Destructor
~colvar_grid_scalar();
/// Constructor from specific sizes arrays
colvar_grid_scalar (std::vector<int> const &nx_i);
/// Constructor from a vector of colvars
colvar_grid_scalar (std::vector<colvar *> &colvars,
bool margin = 0);
/// Accumulate the value
inline void acc_value (std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0)
{
// only legal value of imult here is 0
data[address (ix)] += new_value;
if (samples)
samples->incr_count (ix);
has_data = true;
}
/// Return the gradient of the scalar field from finite differences
inline const cvm::real * gradient_finite_diff ( const std::vector<int> &ix0 )
{
cvm::real A0, A1;
std::vector<int> ix;
if (nd != 2) cvm::fatal_error ("Finite differences available in dimension 2 only.");
for (int n = 0; n < nd; n++) {
ix = ix0;
A0 = data[address (ix)];
ix[n]++; wrap (ix);
A1 = data[address (ix)];
ix[1-n]++; wrap (ix);
A1 += data[address (ix)];
ix[n]--; wrap (ix);
A0 += data[address (ix)];
grad[n] = 0.5 * (A1 - A0) / widths[n];
}
return grad;
}
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual cvm::real value_output (std::vector<int> const &ix,
size_t const &imult = 0)
{
if (imult > 0)
cvm::fatal_error ("Error: trying to access a component "
"larger than 1 in a scalar data grid.\n");
if (samples)
return (samples->value (ix) > 0) ?
(data[address (ix)] / cvm::real (samples->value (ix))) :
0.0;
else
return data[address (ix)];
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual void value_input (std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0,
bool add = false)
{
if (imult > 0)
cvm::fatal_error ("Error: trying to access a component "
"larger than 1 in a scalar data grid.\n");
if (add) {
if (samples)
data[address (ix)] += new_value * samples->new_count (ix);
else
data[address (ix)] += new_value;
} else {
if (samples)
data[address (ix)] = new_value * samples->value (ix);
else
data[address (ix)] = new_value;
}
has_data = true;
}
/// \brief Read the grid from a restart
std::istream & read_restart (std::istream &is);
/// \brief Write the grid to a restart
std::ostream & write_restart (std::ostream &os);
/// \brief Return the highest value
inline cvm::real maximum_value()
{
cvm::real max = data[0];
for (size_t i = 0; i < nt; i++) {
if (data[i] > max) max = data[i];
}
return max;
}
/// \brief Return the lowest value
inline cvm::real minimum_value()
{
cvm::real min = data[0];
for (size_t i = 0; i < nt; i++) {
if (data[i] < min) min = data[i];
}
return min;
}
private:
// gradient
cvm::real * grad;
};
/// Class for accumulating the gradient of a scalar function on a grid
class colvar_grid_gradient : public colvar_grid<cvm::real>
{
public:
/// \brief Provide the sample count by which each binned value
/// should be divided
colvar_grid_count *samples;
/// Default constructor
colvar_grid_gradient();
/// Destructor
virtual inline ~colvar_grid_gradient()
{}
/// Constructor from specific sizes arrays
colvar_grid_gradient (std::vector<int> const &nx_i);
/// Constructor from a vector of colvars
colvar_grid_gradient (std::vector<colvar *> &colvars);
/// \brief Accumulate the gradient
inline void acc_grad (std::vector<int> const &ix, cvm::real const *grads) {
for (size_t imult = 0; imult < mult; imult++) {
data[address (ix) + imult] += grads[imult];
}
if (samples)
samples->incr_count (ix);
}
/// \brief Accumulate the gradient based on the force (i.e. sums the
/// opposite of the force)
inline void acc_force (std::vector<int> const &ix, cvm::real const *forces) {
for (size_t imult = 0; imult < mult; imult++) {
data[address (ix) + imult] -= forces[imult];
}
if (samples)
samples->incr_count (ix);
}
/// \brief Return the value of the function at ix divided by its
/// number of samples (if the count grid is defined)
virtual inline cvm::real value_output (std::vector<int> const &ix,
size_t const &imult = 0)
{
if (samples)
return (samples->value (ix) > 0) ?
(data[address (ix) + imult] / cvm::real (samples->value (ix))) :
0.0;
else
return data[address (ix) + imult];
}
/// \brief Get the value from a formatted output and transform it
/// into the internal representation (it may have been rescaled or
/// manipulated)
virtual inline void value_input (std::vector<int> const &ix,
cvm::real const &new_value,
size_t const &imult = 0,
bool add = false)
{
if (add) {
if (samples)
data[address (ix) + imult] += new_value * samples->new_count (ix);
else
data[address (ix) + imult] += new_value;
} else {
if (samples)
data[address (ix) + imult] = new_value * samples->value (ix);
else
data[address (ix) + imult] = new_value;
}
has_data = true;
}
/// \brief Read the grid from a restart
std::istream & read_restart (std::istream &is);
/// \brief Write the grid to a restart
std::ostream & write_restart (std::ostream &os);
/// Compute and return average value for a 1D gradient grid
inline cvm::real average()
{
size_t n = 0;
if (nd != 1 || nx[0] == 0) {
return 0.0;
}
cvm::real sum = 0.0;
std::vector<int> ix = new_index();
if (samples) {
for ( ; index_ok (ix); incr (ix)) {
if ( (n = samples->value (ix)) )
sum += value (ix) / n;
}
} else {
for ( ; index_ok (ix); incr (ix)) {
sum += value (ix);
}
}
return (sum / cvm::real (nx[0]));
}
/// \brief If the grid is 1-dimensional, integrate it and write the
/// integral to a file
void write_1D_integral (std::ostream &os);
};
#endif
diff --git a/lib/colvars/colvarmodule.C b/lib/colvars/colvarmodule.C
index 43d459dd9..6a9f5a3e8 100644
--- a/lib/colvars/colvarmodule.C
+++ b/lib/colvars/colvarmodule.C
@@ -1,1357 +1,1360 @@
#include "colvarmodule.h"
#include "colvarparse.h"
#include "colvarproxy.h"
#include "colvar.h"
#include "colvarbias.h"
#include "colvarbias_meta.h"
#include "colvarbias_abf.h"
colvarmodule::colvarmodule (char const *config_filename,
colvarproxy *proxy_in)
{
// pointer to the proxy object
if (proxy == NULL) {
proxy = proxy_in;
parse = new colvarparse();
} else {
cvm::fatal_error ("Error: trying to allocate twice the collective "
"variable module.\n");
}
cvm::log (cvm::line_marker);
cvm::log ("Initializing the collective variables module, version "+
cvm::to_str(COLVARS_VERSION)+".\n");
- // "it_restart" will be set by the input restart file, if any;
+ // "it_restart" will be set by the input state file, if any;
// "it" should be updated by the proxy
it = it_restart = 0;
+ it_restart_from_state_file = true;
// open the configfile
config_s.open (config_filename);
if (!config_s)
cvm::fatal_error ("Error: in opening configuration file \""+
std::string (config_filename)+"\".\n");
// read the config file into a string
std::string conf = "";
{
std::string line;
while (colvarparse::getline_nocomments (config_s, line))
conf.append (line+"\n");
// don't need the stream any more
config_s.close();
}
parse->get_keyval (conf, "analysis", b_analysis, false);
if (cvm::debug())
parse->get_keyval (conf, "debugGradientsStepSize", debug_gradients_step_size, 1.0e-03,
colvarparse::parse_silent);
parse->get_keyval (conf, "eigenvalueCrossingThreshold",
colvarmodule::rotation::crossing_threshold, 1.0e-04,
colvarparse::parse_silent);
parse->get_keyval (conf, "colvarsTrajFrequency", cv_traj_freq, 100);
parse->get_keyval (conf, "colvarsRestartFrequency", restart_out_freq,
proxy->restart_frequency());
// by default overwrite the existing trajectory file
parse->get_keyval (conf, "colvarsTrajAppend", cv_traj_append, false);
// input restart file
restart_in_name = proxy->input_prefix().size() ?
std::string (proxy->input_prefix()+".colvars.state") :
std::string ("") ;
// output restart file
restart_out_name = proxy->restart_output_prefix().size() ?
std::string (proxy->restart_output_prefix()+".colvars.state") :
std::string ("");
if (restart_out_name.size())
cvm::log ("The restart output state file will be \""+restart_out_name+"\".\n");
output_prefix = proxy->output_prefix();
cvm::log ("The final output state file will be \""+
(output_prefix.size() ?
std::string (output_prefix+".colvars.state") :
std::string ("colvars.state"))+"\".\n");
cv_traj_name =
(output_prefix.size() ?
std::string (output_prefix+".colvars.traj") :
std::string ("colvars.traj"));
cvm::log ("The trajectory file will be \""+
cv_traj_name+"\".\n");
// open trajectory file
if (cv_traj_append) {
cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+
"\".\n");
cv_traj_os.open (cv_traj_name.c_str(), std::ios::app);
} else {
proxy->backup_file (cv_traj_name.c_str());
cv_traj_os.open (cv_traj_name.c_str(), std::ios::out);
}
cv_traj_os.setf (std::ios::scientific, std::ios::floatfield);
// parse the options for collective variables
init_colvars (conf);
// parse the options for biases
init_biases (conf);
// done with the parsing, check that all keywords are valid
parse->check_keywords (conf, "colvarmodule");
cvm::log (cvm::line_marker);
// read the restart configuration, if available
if (restart_in_name.size()) {
// read the restart file
std::ifstream input_is (restart_in_name.c_str());
if (!input_is.good())
fatal_error ("Error: in opening restart file \""+
std::string (restart_in_name)+"\".\n");
else {
cvm::log ("Restarting from file \""+restart_in_name+"\".\n");
read_restart (input_is);
cvm::log (cvm::line_marker);
}
}
// check if it is possible to save output configuration
if ((!output_prefix.size()) && (!restart_out_name.size())) {
cvm::fatal_error ("Error: neither the final output state file or "
"the output restart file could be defined, exiting.\n");
}
cvm::log ("Collective variables module initialized.\n");
cvm::log (cvm::line_marker);
}
std::istream & colvarmodule::read_restart (std::istream &is)
{
{
// read global restart information
std::string restart_conf;
if (is >> colvarparse::read_block ("configuration", restart_conf)) {
- parse->get_keyval (restart_conf, "step",
- it_restart, (size_t) 0,
- colvarparse::parse_silent);
- it = it_restart;
+ if (it_restart_from_state_file) {
+ parse->get_keyval (restart_conf, "step",
+ it_restart, (size_t) 0,
+ colvarparse::parse_silent);
+ it = it_restart;
+ }
}
is.clear();
}
// colvars restart
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
if ( !((*cvi)->read_restart (is)) )
cvm::fatal_error ("Error: in reading restart configuration for collective variable \""+
(*cvi)->name+"\".\n");
}
// biases restart
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if (!((*bi)->read_restart (is)))
fatal_error ("Error: in reading restart configuration for bias \""+
(*bi)->name+"\".\n");
}
cvm::decrease_depth();
return is;
}
void colvarmodule::init_colvars (std::string const &conf)
{
if (cvm::debug())
cvm::log ("Initializing the collective variables.\n");
std::string colvar_conf = "";
size_t pos = 0;
while (parse->key_lookup (conf, "colvar", colvar_conf, pos)) {
if (colvar_conf.size()) {
cvm::log (cvm::line_marker);
cvm::increase_depth();
colvars.push_back (new colvar (colvar_conf));
(colvars.back())->check_keywords (colvar_conf, "colvar");
cvm::decrease_depth();
} else {
cvm::log ("Warning: \"colvar\" keyword found without any configuration.\n");
}
colvar_conf = "";
}
if (!colvars.size())
cvm::fatal_error ("Error: no collective variables defined.\n");
if (colvars.size())
cvm::log (cvm::line_marker);
cvm::log ("Collective variables initialized, "+
cvm::to_str (colvars.size())+
" in total.\n");
}
void colvarmodule::init_biases (std::string const &conf)
{
if (cvm::debug())
cvm::log ("Initializing the collective variables biases.\n");
{
/// initialize ABF instances
std::string abf_conf = "";
size_t abf_pos = 0;
while (parse->key_lookup (conf, "abf", abf_conf, abf_pos)) {
if (abf_conf.size()) {
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_abf (abf_conf, "abf"));
(biases.back())->check_keywords (abf_conf, "abf");
cvm::decrease_depth();
n_abf_biases++;
} else {
cvm::log ("Warning: \"abf\" keyword found without configuration.\n");
}
abf_conf = "";
}
}
{
/// initialize harmonic restraints
std::string harm_conf = "";
size_t harm_pos = 0;
while (parse->key_lookup (conf, "harmonic", harm_conf, harm_pos)) {
if (harm_conf.size()) {
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_harmonic (harm_conf, "harmonic"));
(biases.back())->check_keywords (harm_conf, "harmonic");
cvm::decrease_depth();
n_harm_biases++;
} else {
cvm::log ("Warning: \"harmonic\" keyword found without configuration.\n");
}
harm_conf = "";
}
}
{
/// initialize histograms
std::string histo_conf = "";
size_t histo_pos = 0;
while (parse->key_lookup (conf, "histogram", histo_conf, histo_pos)) {
if (histo_conf.size()) {
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_histogram (histo_conf, "histogram"));
(biases.back())->check_keywords (histo_conf, "histogram");
cvm::decrease_depth();
n_histo_biases++;
} else {
cvm::log ("Warning: \"histogram\" keyword found without configuration.\n");
}
histo_conf = "";
}
}
{
/// initialize metadynamics instances
std::string meta_conf = "";
size_t meta_pos = 0;
while (parse->key_lookup (conf, "metadynamics", meta_conf, meta_pos)) {
if (meta_conf.size()) {
cvm::log (cvm::line_marker);
cvm::increase_depth();
biases.push_back (new colvarbias_meta (meta_conf, "metadynamics"));
(biases.back())->check_keywords (meta_conf, "metadynamics");
cvm::decrease_depth();
n_meta_biases++;
} else {
cvm::log ("Warning: \"metadynamics\" keyword found without configuration.\n");
}
meta_conf = "";
}
}
if (biases.size())
cvm::log (cvm::line_marker);
cvm::log ("Collective variables biases initialized, "+
cvm::to_str (biases.size())+" in total.\n");
}
void colvarmodule::change_configuration(
std::string const &name, std::string const &conf)
{
cvm::increase_depth();
int found = 0;
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ( (*bi)->name == name ) {
++found;
(*bi)->change_configuration(conf);
}
}
if ( found < 1 ) cvm::fatal_error ("Error: bias not found");
if ( found > 1 ) cvm::fatal_error ("Error: duplicate bias name");
cvm::decrease_depth();
}
cvm::real colvarmodule::energy_difference(
std::string const &name, std::string const &conf)
{
cvm::increase_depth();
cvm::real energy_diff = 0.;
int found = 0;
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ( (*bi)->name == name ) {
++found;
energy_diff = (*bi)->energy_difference(conf);
}
}
if ( found < 1 ) cvm::fatal_error ("Error: bias not found");
if ( found > 1 ) cvm::fatal_error ("Error: duplicate bias name");
cvm::decrease_depth();
return energy_diff;
}
void colvarmodule::calc() {
cvm::real total_bias_energy = 0.0;
cvm::real total_colvar_energy = 0.0;
if (cvm::debug()) {
cvm::log (cvm::line_marker);
cvm::log ("Collective variables module, step no. "+
cvm::to_str (cvm::step_absolute())+"\n");
}
// calculate collective variables and their gradients
if (cvm::debug())
cvm::log ("Calculating collective variables.\n");
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->calc();
}
cvm::decrease_depth();
// update the biases and communicate their forces to the collective
// variables
if (cvm::debug() && biases.size())
cvm::log ("Updating collective variable biases.\n");
cvm::increase_depth();
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
total_bias_energy += (*bi)->update();
}
cvm::decrease_depth();
// sum the forces from all biases for each collective variable
if (cvm::debug() && biases.size())
cvm::log ("Collecting forces from all biases.\n");
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->reset_bias_force();
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->communicate_forces();
}
cvm::decrease_depth();
if (cvm::b_analysis) {
// perform runtime analysis of colvars and biases
if (cvm::debug() && biases.size())
cvm::log ("Perform runtime analyses.\n");
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->analyse();
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->analyse();
}
cvm::decrease_depth();
}
// sum up the forces for each colvar and integrate any internal
// equation of motion
if (cvm::debug())
cvm::log ("Updating the internal degrees of freedom "
"of colvars (if they have any).\n");
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
total_colvar_energy += (*cvi)->update();
}
cvm::decrease_depth();
proxy->add_energy (total_bias_energy + total_colvar_energy);
// make collective variables communicate their forces to their
// coupled degrees of freedom (i.e. atoms)
if (cvm::debug())
cvm::log ("Communicating forces from the colvars to the atoms.\n");
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
if ((*cvi)->tasks[colvar::task_gradients])
(*cvi)->communicate_forces();
}
cvm::decrease_depth();
// write restart file, if needed
if (restart_out_freq && !cvm::b_analysis) {
if ( (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log ("Writing the state file \""+
restart_out_name+"\".\n");
proxy->backup_file (restart_out_name.c_str());
restart_out_os.open (restart_out_name.c_str());
restart_out_os.setf (std::ios::scientific, std::ios::floatfield);
if (!write_restart (restart_out_os))
cvm::fatal_error ("Error: in writing restart file.\n");
restart_out_os.close();
}
}
// write trajectory file, if needed
if (cv_traj_freq) {
if (cvm::debug())
cvm::log ("Writing trajectory file.\n");
// (re)open trajectory file
if (!cv_traj_os.good()) {
if (cv_traj_append) {
cvm::log ("Appending to colvar trajectory file \""+cv_traj_name+
"\".\n");
cv_traj_os.open (cv_traj_name.c_str(), std::ios::app);
} else {
cvm::log ("Overwriting colvar trajectory file \""+cv_traj_name+
"\".\n");
proxy->backup_file (cv_traj_name.c_str());
cv_traj_os.open (cv_traj_name.c_str(), std::ios::out);
}
cv_traj_os.setf (std::ios::scientific, std::ios::floatfield);
}
// write labels in the traj file every 1000 lines and at first ts
cvm::increase_depth();
if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
cv_traj_os << "# " << cvm::wrap_string ("step", cvm::it_width-2)
<< " ";
if (cvm::debug())
cv_traj_os.flush();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_traj_label (cv_traj_os);
}
cv_traj_os << "\n";
if (cvm::debug())
cv_traj_os.flush();
}
cvm::decrease_depth();
// write collective variable values to the traj file
cvm::increase_depth();
if ((cvm::step_absolute() % cv_traj_freq) == 0) {
cv_traj_os << std::setw (cvm::it_width) << it
<< " ";
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_traj (cv_traj_os);
}
cv_traj_os << "\n";
if (cvm::debug())
cv_traj_os.flush();
}
cvm::decrease_depth();
if (restart_out_freq) {
// flush the trajectory file if we are at the restart frequency
if ( (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log ("Synchronizing (emptying the buffer of) trajectory file \""+
cv_traj_name+"\".\n");
cv_traj_os.flush();
}
}
} // end if (cv_traj_freq)
}
void colvarmodule::analyze()
{
if (cvm::debug()) {
cvm::log ("colvarmodule::analyze(), step = "+cvm::to_str (it)+".\n");
}
if (cvm::step_relative() == 0)
cvm::log ("Performing analysis.\n");
// perform colvar-specific analysis
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
cvm::increase_depth();
(*cvi)->analyse();
cvm::decrease_depth();
}
// perform bias-specific analysis
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
cvm::increase_depth();
(*bi)->analyse();
cvm::decrease_depth();
}
}
colvarmodule::~colvarmodule()
{
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
delete *cvi;
}
colvars.clear();
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
delete *bi;
}
biases.clear();
if (cv_traj_os.good()) {
cv_traj_os.close();
}
delete parse;
proxy = NULL;
}
void colvarmodule::write_output_files()
{
// if this is a simulation run (i.e. not a postprocessing), output data
// must be written to be able to restart the simulation
std::string const out_name =
(output_prefix.size() ?
std::string (output_prefix+".colvars.state") :
std::string ("colvars.state"));
cvm::log ("Saving collective variables state to \""+out_name+"\".\n");
proxy->backup_file (out_name.c_str());
std::ofstream out (out_name.c_str());
out.setf (std::ios::scientific, std::ios::floatfield);
this->write_restart (out);
out.close();
// do not close to avoid problems with multiple NAMD runs
cv_traj_os.flush();
}
bool colvarmodule::read_traj (char const *traj_filename,
size_t traj_read_begin,
size_t traj_read_end)
{
cvm::log ("Opening trajectory file \""+
std::string (traj_filename)+"\".\n");
std::ifstream traj_is (traj_filename);
while (true) {
while (true) {
std::string line ("");
do {
if (!colvarparse::getline_nocomments (traj_is, line)) {
cvm::log ("End of file \""+std::string (traj_filename)+
"\" reached, or corrupted file.\n");
traj_is.close();
return false;
}
} while (line.find_first_not_of (colvarparse::white_space) == std::string::npos);
std::istringstream is (line);
if (!(is >> it)) return false;
if ( (it < traj_read_begin) ) {
if ((it % 1000) == 0)
std::cerr << "Skipping trajectory step " << it
<< " \r";
continue;
} else {
if ((it % 1000) == 0)
std::cerr << "Reading from trajectory, step = " << it
<< " \r";
if ( (traj_read_end > traj_read_begin) &&
(it > traj_read_end) ) {
std::cerr << "\n";
cvm::log ("Reached the end of the trajectory, "
"read_end = "+cvm::to_str (traj_read_end)+"\n");
return false;
}
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
if (!(*cvi)->read_traj (is)) {
cvm::log ("Error: in reading colvar \""+(*cvi)->name+
"\" from trajectory file \""+
std::string (traj_filename)+"\".\n");
return false;
}
}
break;
}
}
}
return true;
}
std::ostream & colvarmodule::write_restart (std::ostream &os)
{
os << "configuration {\n"
<< " step " << std::setw (it_width)
<< it << "\n"
<< " dt " << dt() << "\n"
<< "}\n\n";
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
(*cvi)->write_restart (os);
}
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
(*bi)->write_restart (os);
}
cvm::decrease_depth();
return os;
}
void cvm::log (std::string const &message)
{
if (depth > 0)
proxy->log ((std::string (2*depth, ' '))+message);
else
proxy->log (message);
}
void cvm::increase_depth()
{
depth++;
}
void cvm::decrease_depth()
{
if (depth) depth--;
}
void cvm::fatal_error (std::string const &message)
{
proxy->fatal_error (message);
}
void cvm::exit (std::string const &message)
{
proxy->exit (message);
}
// static pointers
std::vector<colvar *> colvarmodule::colvars;
std::vector<colvarbias *> colvarmodule::biases;
size_t colvarmodule::n_abf_biases = 0;
size_t colvarmodule::n_harm_biases = 0;
size_t colvarmodule::n_histo_biases = 0;
size_t colvarmodule::n_meta_biases = 0;
colvarproxy *colvarmodule::proxy = NULL;
// static runtime data
cvm::real colvarmodule::debug_gradients_step_size = 1.0e-03;
size_t colvarmodule::it = 0;
size_t colvarmodule::it_restart = 0;
size_t colvarmodule::restart_out_freq = 0;
size_t colvarmodule::cv_traj_freq = 0;
size_t colvarmodule::depth = 0;
bool colvarmodule::b_analysis = false;
cvm::real colvarmodule::rotation::crossing_threshold = 1.0E-04;
// file name prefixes
std::string colvarmodule::output_prefix = "";
std::string colvarmodule::input_prefix = "";
std::string colvarmodule::restart_in_name = "";
// i/o constants
size_t const colvarmodule::it_width = 12;
size_t const colvarmodule::cv_prec = 14;
size_t const colvarmodule::cv_width = 21;
size_t const colvarmodule::en_prec = 14;
size_t const colvarmodule::en_width = 21;
std::string const colvarmodule::line_marker =
"----------------------------------------------------------------------\n";
std::ostream & operator << (std::ostream &os, colvarmodule::rvector const &v)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
os.width (2);
os << "( ";
os.width (w); os.precision (p);
os << v.x << " , ";
os.width (w); os.precision (p);
os << v.y << " , ";
os.width (w); os.precision (p);
os << v.z << " )";
return os;
}
std::istream & operator >> (std::istream &is, colvarmodule::rvector &v)
{
size_t const start_pos = is.tellg();
char sep;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> v.x) || !(is >> sep) || !(sep == ',') ||
!(is >> v.y) || !(is >> sep) || !(sep == ',') ||
!(is >> v.z) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
return is;
}
std::ostream & operator << (std::ostream &os, colvarmodule::quaternion const &q)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
os.width (2);
os << "( ";
os.width (w); os.precision (p);
os << q.q0 << " , ";
os.width (w); os.precision (p);
os << q.q1 << " , ";
os.width (w); os.precision (p);
os << q.q2 << " , ";
os.width (w); os.precision (p);
os << q.q3 << " )";
return os;
}
std::istream & operator >> (std::istream &is, colvarmodule::quaternion &q)
{
size_t const start_pos = is.tellg();
std::string euler ("");
if ( (is >> euler) && (colvarparse::to_lower_cppstr (euler) ==
std::string ("euler")) ) {
// parse the Euler angles
char sep;
cvm::real phi, theta, psi;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> phi) || !(is >> sep) || !(sep == ',') ||
!(is >> theta) || !(is >> sep) || !(sep == ',') ||
!(is >> psi) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
q = colvarmodule::quaternion (phi, theta, psi);
} else {
// parse the quaternion components
is.seekg (start_pos, std::ios::beg);
char sep;
if ( !(is >> sep) || !(sep == '(') ||
!(is >> q.q0) || !(is >> sep) || !(sep == ',') ||
!(is >> q.q1) || !(is >> sep) || !(sep == ',') ||
!(is >> q.q2) || !(is >> sep) || !(sep == ',') ||
!(is >> q.q3) || !(is >> sep) || !(sep == ')') ) {
is.clear();
is.seekg (start_pos, std::ios::beg);
is.setstate (std::ios::failbit);
return is;
}
}
return is;
}
cvm::quaternion
cvm::quaternion::position_derivative_inner (cvm::rvector const &pos,
cvm::rvector const &vec) const
{
cvm::quaternion result (0.0, 0.0, 0.0, 0.0);
result.q0 = 2.0 * pos.x * q0 * vec.x
+2.0 * pos.y * q0 * vec.y
+2.0 * pos.z * q0 * vec.z
-2.0 * pos.y * q3 * vec.x
+2.0 * pos.z * q2 * vec.x
+2.0 * pos.x * q3 * vec.y
-2.0 * pos.z * q1 * vec.y
-2.0 * pos.x * q2 * vec.z
+2.0 * pos.y * q1 * vec.z;
result.q1 = +2.0 * pos.x * q1 * vec.x
-2.0 * pos.y * q1 * vec.y
-2.0 * pos.z * q1 * vec.z
+2.0 * pos.y * q2 * vec.x
+2.0 * pos.z * q3 * vec.x
+2.0 * pos.x * q2 * vec.y
-2.0 * pos.z * q0 * vec.y
+2.0 * pos.x * q3 * vec.z
+2.0 * pos.y * q0 * vec.z;
result.q2 = -2.0 * pos.x * q2 * vec.x
+2.0 * pos.y * q2 * vec.y
-2.0 * pos.z * q2 * vec.z
+2.0 * pos.y * q1 * vec.x
+2.0 * pos.z * q0 * vec.x
+2.0 * pos.x * q1 * vec.y
+2.0 * pos.z * q3 * vec.y
-2.0 * pos.x * q0 * vec.z
+2.0 * pos.y * q3 * vec.z;
result.q3 = -2.0 * pos.x * q3 * vec.x
-2.0 * pos.y * q3 * vec.y
+2.0 * pos.z * q3 * vec.z
-2.0 * pos.y * q0 * vec.x
+2.0 * pos.z * q1 * vec.x
+2.0 * pos.x * q0 * vec.y
+2.0 * pos.z * q2 * vec.y
+2.0 * pos.x * q1 * vec.z
+2.0 * pos.y * q2 * vec.z;
return result;
}
// Calculate the optimal rotation between two groups, and implement it
// as a quaternion. The method is the one documented in: Coutsias EA,
// Seok C, Dill KA. Using quaternions to calculate RMSD. J Comput
// Chem. 25(15):1849-57 (2004) DOI: 10.1002/jcc.20110 PubMed: 15376254
void colvarmodule::rotation::build_matrix (std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2,
matrix2d<cvm::real, 4, 4> &S)
{
cvm::rmatrix C;
// build the correlation matrix
C.reset();
for (size_t i = 0; i < pos1.size(); i++) {
C.xx() += pos1[i].x * pos2[i].x;
C.xy() += pos1[i].x * pos2[i].y;
C.xz() += pos1[i].x * pos2[i].z;
C.yx() += pos1[i].y * pos2[i].x;
C.yy() += pos1[i].y * pos2[i].y;
C.yz() += pos1[i].y * pos2[i].z;
C.zx() += pos1[i].z * pos2[i].x;
C.zy() += pos1[i].z * pos2[i].y;
C.zz() += pos1[i].z * pos2[i].z;
}
// build the "overlap" matrix, whose eigenvectors are stationary
// points of the RMSD in the space of rotations
S[0][0] = C.xx() + C.yy() + C.zz();
S[1][0] = C.yz() - C.zy();
S[0][1] = S[1][0];
S[2][0] = - C.xz() + C.zx() ;
S[0][2] = S[2][0];
S[3][0] = C.xy() - C.yx();
S[0][3] = S[3][0];
S[1][1] = C.xx() - C.yy() - C.zz();
S[2][1] = C.xy() + C.yx();
S[1][2] = S[2][1];
S[3][1] = C.xz() + C.zx();
S[1][3] = S[3][1];
S[2][2] = - C.xx() + C.yy() - C.zz();
S[3][2] = C.yz() + C.zy();
S[2][3] = S[3][2];
S[3][3] = - C.xx() - C.yy() + C.zz();
// if (cvm::debug()) {
// for (size_t i = 0; i < 4; i++) {
// std::string line ("");
// for (size_t j = 0; j < 4; j++) {
// line += std::string (" S["+cvm::to_str (i)+
// "]["+cvm::to_str (j)+"] ="+cvm::to_str (S[i][j]));
// }
// cvm::log (line+"\n");
// }
// }
}
void colvarmodule::rotation::diagonalize_matrix (matrix2d<cvm::real, 4, 4> &S,
cvm::real S_eigval[4],
matrix2d<cvm::real, 4, 4> &S_eigvec)
{
// diagonalize
int jac_nrot = 0;
jacobi (S, 4, S_eigval, S_eigvec, &jac_nrot);
eigsrt (S_eigval, S_eigvec, 4);
// jacobi saves eigenvectors by columns
transpose (S_eigvec, 4);
// normalize eigenvectors
for (size_t ie = 0; ie < 4; ie++) {
cvm::real norm2 = 0.0;
for (size_t i = 0; i < 4; i++) norm2 += std::pow (S_eigvec[ie][i], int (2));
cvm::real const norm = std::sqrt (norm2);
for (size_t i = 0; i < 4; i++) S_eigvec[ie][i] /= norm;
}
}
// Calculate the rotation, plus its derivatives
void colvarmodule::rotation::calc_optimal_rotation
(std::vector<cvm::atom_pos> const &pos1,
std::vector<cvm::atom_pos> const &pos2)
{
matrix2d<cvm::real, 4, 4> S;
matrix2d<cvm::real, 4, 4> S_backup;
cvm::real S_eigval[4];
matrix2d<cvm::real, 4, 4> S_eigvec;
// if (cvm::debug()) {
// cvm::atom_pos cog1 (0.0, 0.0, 0.0);
// for (size_t i = 0; i < pos1.size(); i++) {
// cog1 += pos1[i];
// }
// cog1 /= cvm::real (pos1.size());
// cvm::atom_pos cog2 (0.0, 0.0, 0.0);
// for (size_t i = 0; i < pos2.size(); i++) {
// cog2 += pos2[i];
// }
// cog2 /= cvm::real (pos1.size());
// cvm::log ("calc_optimal_rotation: centers of geometry are: "+
// cvm::to_str (cog1, cvm::cv_width, cvm::cv_prec)+
// " and "+cvm::to_str (cog2, cvm::cv_width, cvm::cv_prec)+".\n");
// }
build_matrix (pos1, pos2, S);
S_backup = S;
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log ("S = "+cvm::to_str (cvm::to_str (S_backup), cvm::cv_width, cvm::cv_prec)+"\n");
}
}
diagonalize_matrix (S, S_eigval, S_eigvec);
// eigenvalues and eigenvectors
cvm::real const &L0 = S_eigval[0];
cvm::real const &L1 = S_eigval[1];
cvm::real const &L2 = S_eigval[2];
cvm::real const &L3 = S_eigval[3];
cvm::real const *Q0 = S_eigvec[0];
cvm::real const *Q1 = S_eigvec[1];
cvm::real const *Q2 = S_eigvec[2];
cvm::real const *Q3 = S_eigvec[3];
lambda = L0;
q = cvm::quaternion (Q0);
if (q_old.norm2() > 0.0) {
q.match (q_old);
if (q_old.inner (q) < (1.0 - crossing_threshold)) {
cvm::log ("Warning: discontinuous rotation!\n");
}
}
q_old = q;
if (cvm::debug()) {
if (b_debug_gradients) {
cvm::log ("L0 = "+cvm::to_str (L0, cvm::cv_width, cvm::cv_prec)+
", Q0 = "+cvm::to_str (cvm::quaternion (Q0), cvm::cv_width, cvm::cv_prec)+
", Q0*Q0 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q0)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L1 = "+cvm::to_str (L1, cvm::cv_width, cvm::cv_prec)+
", Q1 = "+cvm::to_str (cvm::quaternion (Q1), cvm::cv_width, cvm::cv_prec)+
", Q0*Q1 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q1)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L2 = "+cvm::to_str (L2, cvm::cv_width, cvm::cv_prec)+
", Q2 = "+cvm::to_str (cvm::quaternion (Q2), cvm::cv_width, cvm::cv_prec)+
", Q0*Q2 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q2)), cvm::cv_width, cvm::cv_prec)+
"\n");
cvm::log ("L3 = "+cvm::to_str (L3, cvm::cv_width, cvm::cv_prec)+
", Q3 = "+cvm::to_str (cvm::quaternion (Q3), cvm::cv_width, cvm::cv_prec)+
", Q0*Q3 = "+cvm::to_str (cvm::quaternion (Q0).inner (cvm::quaternion (Q3)), cvm::cv_width, cvm::cv_prec)+
"\n");
}
}
// calculate derivatives of L0 and Q0 with respect to each atom in
// either group; note: if dS_1 is a null vector, nothing will be
// calculated
for (size_t ia = 0; ia < dS_1.size(); ia++) {
cvm::real const &a2x = pos2[ia].x;
cvm::real const &a2y = pos2[ia].y;
cvm::real const &a2z = pos2[ia].z;
matrix2d<cvm::rvector, 4, 4> &ds_1 = dS_1[ia];
// derivative of the S matrix
ds_1.reset();
ds_1[0][0] = cvm::rvector ( a2x, a2y, a2z);
ds_1[1][0] = cvm::rvector ( 0.0, a2z, -a2y);
ds_1[0][1] = ds_1[1][0];
ds_1[2][0] = cvm::rvector (-a2z, 0.0, a2x);
ds_1[0][2] = ds_1[2][0];
ds_1[3][0] = cvm::rvector ( a2y, -a2x, 0.0);
ds_1[0][3] = ds_1[3][0];
ds_1[1][1] = cvm::rvector ( a2x, -a2y, -a2z);
ds_1[2][1] = cvm::rvector ( a2y, a2x, 0.0);
ds_1[1][2] = ds_1[2][1];
ds_1[3][1] = cvm::rvector ( a2z, 0.0, a2x);
ds_1[1][3] = ds_1[3][1];
ds_1[2][2] = cvm::rvector (-a2x, a2y, -a2z);
ds_1[3][2] = cvm::rvector ( 0.0, a2z, a2y);
ds_1[2][3] = ds_1[3][2];
ds_1[3][3] = cvm::rvector (-a2x, -a2y, a2z);
cvm::rvector &dl0_1 = dL0_1[ia];
vector1d<cvm::rvector, 4> &dq0_1 = dQ0_1[ia];
// matrix multiplications; derivatives of L_0 and Q_0 are
// calculated using Hellmann-Feynman theorem (i.e. exploiting the
// fact that the eigenvectors Q_i form an orthonormal basis)
dl0_1.reset();
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dl0_1 += Q0[i] * ds_1[i][j] * Q0[j];
}
}
dq0_1.reset();
for (size_t p = 0; p < 4; p++) {
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dq0_1[p] +=
(Q1[i] * ds_1[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
(Q2[i] * ds_1[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
(Q3[i] * ds_1[i][j] * Q0[j]) / (L0-L3) * Q3[p];
}
}
}
}
// do the same for the second group
for (size_t ia = 0; ia < dS_2.size(); ia++) {
cvm::real const &a1x = pos1[ia].x;
cvm::real const &a1y = pos1[ia].y;
cvm::real const &a1z = pos1[ia].z;
matrix2d<cvm::rvector, 4, 4> &ds_2 = dS_2[ia];
ds_2.reset();
ds_2[0][0] = cvm::rvector ( a1x, a1y, a1z);
ds_2[1][0] = cvm::rvector ( 0.0, -a1z, a1y);
ds_2[0][1] = ds_2[1][0];
ds_2[2][0] = cvm::rvector ( a1z, 0.0, -a1x);
ds_2[0][2] = ds_2[2][0];
ds_2[3][0] = cvm::rvector (-a1y, a1x, 0.0);
ds_2[0][3] = ds_2[3][0];
ds_2[1][1] = cvm::rvector ( a1x, -a1y, -a1z);
ds_2[2][1] = cvm::rvector ( a1y, a1x, 0.0);
ds_2[1][2] = ds_2[2][1];
ds_2[3][1] = cvm::rvector ( a1z, 0.0, a1x);
ds_2[1][3] = ds_2[3][1];
ds_2[2][2] = cvm::rvector (-a1x, a1y, -a1z);
ds_2[3][2] = cvm::rvector ( 0.0, a1z, a1y);
ds_2[2][3] = ds_2[3][2];
ds_2[3][3] = cvm::rvector (-a1x, -a1y, a1z);
cvm::rvector &dl0_2 = dL0_2[ia];
vector1d<cvm::rvector, 4> &dq0_2 = dQ0_2[ia];
dl0_2.reset();
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dl0_2 += Q0[i] * ds_2[i][j] * Q0[j];
}
}
dq0_2.reset();
for (size_t p = 0; p < 4; p++) {
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
dq0_2[p] +=
(Q1[i] * ds_2[i][j] * Q0[j]) / (L0-L1) * Q1[p] +
(Q2[i] * ds_2[i][j] * Q0[j]) / (L0-L2) * Q2[p] +
(Q3[i] * ds_2[i][j] * Q0[j]) / (L0-L3) * Q3[p];
}
}
}
if (cvm::debug()) {
if (b_debug_gradients) {
matrix2d<cvm::real, 4, 4> S_new;
cvm::real S_new_eigval[4];
matrix2d<cvm::real, 4, 4> S_new_eigvec;
// make an infitesimal move along each cartesian coordinate of
// this atom, and solve again the eigenvector problem
for (size_t comp = 0; comp < 3; comp++) {
S_new = S_backup;
// diagonalize the new overlap matrix
for (size_t i = 0; i < 4; i++) {
for (size_t j = 0; j < 4; j++) {
S_new[i][j] +=
colvarmodule::debug_gradients_step_size * ds_2[i][j][comp];
}
}
// cvm::log ("S_new = "+cvm::to_str (cvm::to_str (S_new), cvm::cv_width, cvm::cv_prec)+"\n");
diagonalize_matrix (S_new, S_new_eigval, S_new_eigvec);
cvm::real const &L0_new = S_new_eigval[0];
cvm::real const *Q0_new = S_new_eigvec[0];
cvm::real const DL0 = (dl0_2[comp]) * colvarmodule::debug_gradients_step_size;
cvm::quaternion const q0 (Q0);
cvm::quaternion const DQ0 (dq0_2[0][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[1][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[2][comp] * colvarmodule::debug_gradients_step_size,
dq0_2[3][comp] * colvarmodule::debug_gradients_step_size);
cvm::log ( "|(l_0+dl_0) - l_0^new|/l_0 = "+
cvm::to_str (std::fabs (L0+DL0 - L0_new)/L0, cvm::cv_width, cvm::cv_prec)+
", |(q_0+dq_0) - q_0^new| = "+
cvm::to_str ((Q0+DQ0 - Q0_new).norm(), cvm::cv_width, cvm::cv_prec)+
"\n");
}
}
}
}
}
// Numerical Recipes routine for diagonalization
#define ROTATE(a,i,j,k,l) g=a[i][j];h=a[k][l];a[i][j]=g-s*(h+g*tau); \
a[k][l]=h+s*(g-h*tau);
void jacobi(cvm::real **a, int n, cvm::real d[], cvm::real **v, int *nrot)
{
int j,iq,ip,i;
cvm::real tresh,theta,tau,t,sm,s,h,g,c;
std::vector<cvm::real> b (n, 0.0);
std::vector<cvm::real> z (n, 0.0);
for (ip=0;ip<n;ip++) {
for (iq=0;iq<n;iq++) v[ip][iq]=0.0;
v[ip][ip]=1.0;
}
for (ip=0;ip<n;ip++) {
b[ip]=d[ip]=a[ip][ip];
z[ip]=0.0;
}
*nrot=0;
for (i=0;i<=50;i++) {
sm=0.0;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++)
sm += std::fabs(a[ip][iq]);
}
if (sm == 0.0) {
return;
}
if (i < 4)
tresh=0.2*sm/(n*n);
else
tresh=0.0;
for (ip=0;ip<n-1;ip++) {
for (iq=ip+1;iq<n;iq++) {
g=100.0*std::fabs(a[ip][iq]);
if (i > 4 && (cvm::real)(std::fabs(d[ip])+g) == (cvm::real)std::fabs(d[ip])
&& (cvm::real)(std::fabs(d[iq])+g) == (cvm::real)std::fabs(d[iq]))
a[ip][iq]=0.0;
else if (std::fabs(a[ip][iq]) > tresh) {
h=d[iq]-d[ip];
if ((cvm::real)(std::fabs(h)+g) == (cvm::real)std::fabs(h))
t=(a[ip][iq])/h;
else {
theta=0.5*h/(a[ip][iq]);
t=1.0/(std::fabs(theta)+std::sqrt(1.0+theta*theta));
if (theta < 0.0) t = -t;
}
c=1.0/std::sqrt(1+t*t);
s=t*c;
tau=s/(1.0+c);
h=t*a[ip][iq];
z[ip] -= h;
z[iq] += h;
d[ip] -= h;
d[iq] += h;
a[ip][iq]=0.0;
for (j=0;j<=ip-1;j++) {
ROTATE(a,j,ip,j,iq)
}
for (j=ip+1;j<=iq-1;j++) {
ROTATE(a,ip,j,j,iq)
}
for (j=iq+1;j<n;j++) {
ROTATE(a,ip,j,iq,j)
}
for (j=0;j<n;j++) {
ROTATE(v,j,ip,j,iq)
}
++(*nrot);
}
}
}
for (ip=0;ip<n;ip++) {
b[ip] += z[ip];
d[ip]=b[ip];
z[ip]=0.0;
}
}
cvm::fatal_error ("Too many iterations in routine jacobi.\n");
}
#undef ROTATE
void eigsrt(cvm::real d[], cvm::real **v, int n)
{
int k,j,i;
cvm::real p;
for (i=0;i<n;i++) {
p=d[k=i];
for (j=i+1;j<n;j++)
if (d[j] >= p) p=d[k=j];
if (k != i) {
d[k]=d[i];
d[i]=p;
for (j=0;j<n;j++) {
p=v[j][i];
v[j][i]=v[j][k];
v[j][k]=p;
}
}
}
}
void transpose(cvm::real **v, int n)
{
cvm::real p;
for (int i=0;i<n;i++) {
for (int j=i+1;j<n;j++) {
p=v[i][j];
v[i][j]=v[j][i];
v[j][i]=p;
}
}
}
diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h
index 5febf2f8e..8ff727df2 100644
--- a/lib/colvars/colvarmodule.h
+++ b/lib/colvars/colvarmodule.h
@@ -1,489 +1,493 @@
#ifndef COLVARMODULE_H
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
-#define COLVARS_VERSION "2012-03-21"
+#define COLVARS_VERSION "2012-03-23"
#endif
#ifndef COLVARS_DEBUG
#define COLVARS_DEBUG false
#endif
/// \file colvarmodule.h
/// \brief Collective variables main module
///
/// This file declares the main class for defining and manipulating
/// collective variables: there should be only one instance of this
/// class, because several variables are made static (i.e. they are
/// shared between all object instances) to be accessed from other
/// objects.
#include <iostream>
#include <iomanip>
#include <string>
#include <sstream>
#include <fstream>
#include <cmath>
#include <vector>
#include <list>
class colvarparse;
class colvar;
class colvarbias;
class colvarproxy;
/// \brief Collective variables module (main class)
///
/// Class to control the collective variables calculation. An object
/// (usually one) of this class is spawned from the MD program,
/// containing all i/o routines and general interface.
///
/// At initialization, the colvarmodule object creates a proxy object
/// to provide a transparent interface between the MD program and the
/// child objects
class colvarmodule {
private:
/// Impossible to initialize the main object without arguments
colvarmodule();
public:
friend class colvarproxy;
/// Defining an abstract real number allows to switch precision
typedef double real;
/// Residue identifier
typedef int residue_id;
class rvector;
template <class T,
size_t const length> class vector1d;
template <class T,
size_t const outer_length,
size_t const inner_length> class matrix2d;
class quaternion;
class rotation;
/// \brief Atom position (different type name from rvector, to make
/// possible future PBC-transparent implementations)
typedef rvector atom_pos;
/// \brief 3x3 matrix of real numbers
class rmatrix;
// allow these classes to access protected data
class atom;
class atom_group;
friend class atom;
friend class atom_group;
typedef std::vector<atom>::iterator atom_iter;
typedef std::vector<atom>::const_iterator atom_const_iter;
/// Current step number
static size_t it;
/// Starting step number for this run
static size_t it_restart;
/// Return the current step number from the beginning of this run
static inline size_t step_relative()
{
return it - it_restart;
}
/// Return the current step number from the beginning of the whole
/// calculation
static inline size_t step_absolute()
{
return it;
}
+ /// If true, get it_restart from the state file; if set to false,
+ /// the MD program is providing it
+ bool it_restart_from_state_file;
+
/// \brief Finite difference step size (if there is no dynamics, or
/// if gradients need to be tested independently from the size of
/// dt)
static real debug_gradients_step_size;
/// Prefix for all output files for this run
static std::string output_prefix;
/// Prefix for files from a previous run (including restart/output)
static std::string input_prefix;
/// input restart file name (determined from input_prefix)
static std::string restart_in_name;
/// Array of collective variables
static std::vector<colvar *> colvars;
/* TODO: implement named CVCs
/// Array of named (reusable) collective variable components
static std::vector<cvc *> cvcs;
/// Named cvcs register themselves at initialization time
inline void register_cvc (cvc *p) {
cvcs.push_back(p);
}
*/
/// Array of collective variable biases
static std::vector<colvarbias *> biases;
/// \brief Number of ABF biases initialized (in normal conditions
/// should be 1)
static size_t n_abf_biases;
/// \brief Number of metadynamics biases initialized (in normal
/// conditions should be 1)
static size_t n_meta_biases;
/// \brief Number of harmonic biases initialized (no limit on the
/// number)
static size_t n_harm_biases;
/// \brief Number of histograms initialized (no limit on the
/// number)
static size_t n_histo_biases;
/// \brief Whether debug output should be enabled (compile-time option)
static inline bool debug()
{
return COLVARS_DEBUG;
}
/// \brief Constructor \param config_name Configuration file name
/// \param restart_name (optional) Restart file name
colvarmodule (char const *config_name,
colvarproxy *proxy_in);
/// Destructor
~colvarmodule();
/// Initialize collective variables
void init_colvars (std::string const &conf);
/// Initialize collective variable biases
void init_biases (std::string const &conf);
/// Load new configuration - force constant and/or centers only
void change_configuration(std::string const &name, std::string const &conf);
/// Calculate change in energy from using alternate configuration
real energy_difference(std::string const &name, std::string const &conf);
/// Calculate collective variables and biases
void calc();
/// Read the input restart file
std::istream & read_restart (std::istream &is);
/// Write the output restart file
std::ostream & write_restart (std::ostream &os);
/// Write all output files (called by the proxy)
void write_output_files();
/// \brief Call colvarproxy::backup_file()
static void backup_file (char const *filename);
/// Perform analysis
void analyze();
/// \brief Read a collective variable trajectory (post-processing
/// only, not called at runtime)
bool read_traj (char const *traj_filename,
size_t traj_read_begin,
size_t traj_read_end);
/// Get the pointer of a colvar from its name (returns NULL if not found)
static colvar * colvar_p (std::string const &name);
/// Quick conversion of an object to a string
template<typename T> static std::string to_str (T const &x,
size_t const &width = 0,
size_t const &prec = 0);
/// Quick conversion of a vector of objects to a string
template<typename T> static std::string to_str (std::vector<T> const &x,
size_t const &width = 0,
size_t const &prec = 0);
/// Reduce the number of characters in a string
static inline std::string wrap_string (std::string const &s,
size_t const &nchars)
{
if (!s.size())
return std::string (nchars, ' ');
else
return ( (s.size() <= size_t (nchars)) ?
(s+std::string (nchars-s.size(), ' ')) :
(std::string (s, 0, nchars)) );
}
/// Number of characters to represent a time step
static size_t const it_width;
/// Number of digits to represent a collective variables value(s)
static size_t const cv_prec;
/// Number of characters to represent a collective variables value(s)
static size_t const cv_width;
/// Number of digits to represent the collective variables energy
static size_t const en_prec;
/// Number of characters to represent the collective variables energy
static size_t const en_width;
/// Line separator in the log output
static std::string const line_marker;
// proxy functions
/// \brief Value of the unit for atomic coordinates with respect to
/// angstroms (used by some variables for hard-coded default values)
static real unit_angstrom();
/// \brief Boltmann constant
static real boltzmann();
/// \brief Temperature of the simulation (K)
static real temperature();
/// \brief Time step of MD integrator (fs)
static real dt();
/// Request calculation of system force from MD engine
static void request_system_force();
/// Print a message to the main log
static void log (std::string const &message);
/// Print a message to the main log and exit with error code
static void fatal_error (std::string const &message);
/// Print a message to the main log and exit normally
static void exit (std::string const &message);
/// \brief Get the distance between two atomic positions with pbcs handled
/// correctly
static rvector position_distance (atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Get the square distance between two positions (with
/// periodic boundary conditions handled transparently)
///
/// Note: in the case of periodic boundary conditions, this provides
/// an analytical square distance (while taking the square of
/// position_distance() would produce leads to a cusp)
static real position_dist2 (atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Get the closest periodic image to a reference position
/// \param pos The position to look for the closest periodic image
/// \param ref_pos (optional) The reference position
static void select_closest_image (atom_pos &pos,
atom_pos const &ref_pos);
/// \brief Perform select_closest_image() on a set of atomic positions
///
/// After that, distance vectors can then be calculated directly,
/// without using position_distance()
static void select_closest_images (std::vector<atom_pos> &pos,
atom_pos const &ref_pos);
/// \brief Create atoms from a file \param filename name of the file
/// (usually a PDB) \param atoms array of the atoms to be allocated
/// \param pdb_field (optiona) if "filename" is a PDB file, use this
/// field to determine which are the atoms to be set
static void load_atoms (char const *filename,
std::vector<atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from a file
/// (usually a PDB); the number of atoms in "filename" must match
/// the number of elements in "pos"
static void load_coords (char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
/// Frequency for collective variables trajectory output
static size_t cv_traj_freq;
/// \brief True if only analysis is performed and not a run
static bool b_analysis;
/// Frequency for saving output restarts
static size_t restart_out_freq;
/// Output restart file name
std::string restart_out_name;
/// Pseudo-random number with Gaussian distribution
static real rand_gaussian (void);
protected:
/// Configuration file
std::ifstream config_s;
/// Configuration file parser object
colvarparse *parse;
/// Name of the trajectory file
std::string cv_traj_name;
/// Collective variables output trajectory file
std::ofstream cv_traj_os;
/// Appending to the existing trajectory file?
bool cv_traj_append;
/// Output restart file
std::ofstream restart_out_os;
/// \brief Pointer to the proxy object, used to retrieve atomic data
/// from the hosting program; it is static in order to be accessible
/// from static functions in the colvarmodule class
static colvarproxy *proxy;
/// \brief Counter for the current depth in the object hierarchy (useg e.g. in outpu
static size_t depth;
public:
/// Increase the depth (number of indentations in the output)
static void increase_depth();
/// Decrease the depth (number of indentations in the output)
static void decrease_depth();
};
/// Shorthand for the frequently used type prefix
typedef colvarmodule cvm;
#include "colvartypes.h"
std::ostream & operator << (std::ostream &os, cvm::rvector const &v);
std::istream & operator >> (std::istream &is, cvm::rvector &v);
template<typename T> std::string cvm::to_str (T const &x,
size_t const &width,
size_t const &prec) {
std::ostringstream os;
if (width) os.width (width);
if (prec) {
os.setf (std::ios::scientific, std::ios::floatfield);
os.precision (prec);
}
os << x;
return os.str();
}
template<typename T> std::string cvm::to_str (std::vector<T> const &x,
size_t const &width,
size_t const &prec) {
if (!x.size()) return std::string ("");
std::ostringstream os;
if (prec) {
os.setf (std::ios::scientific, std::ios::floatfield);
}
os << "{ ";
if (width) os.width (width);
if (prec) os.precision (prec);
os << x[0];
for (size_t i = 1; i < x.size(); i++) {
os << ", ";
if (width) os.width (width);
if (prec) os.precision (prec);
os << x[i];
}
os << " }";
return os.str();
}
#include "colvarproxy.h"
inline cvm::real cvm::unit_angstrom()
{
return proxy->unit_angstrom();
}
inline cvm::real cvm::boltzmann()
{
return proxy->boltzmann();
}
inline cvm::real cvm::temperature()
{
return proxy->temperature();
}
inline cvm::real cvm::dt()
{
return proxy->dt();
}
inline void cvm::request_system_force()
{
proxy->request_system_force (true);
}
inline void cvm::select_closest_image (atom_pos &pos,
atom_pos const &ref_pos)
{
proxy->select_closest_image (pos, ref_pos);
}
inline void cvm::select_closest_images (std::vector<atom_pos> &pos,
atom_pos const &ref_pos)
{
proxy->select_closest_images (pos, ref_pos);
}
inline cvm::rvector cvm::position_distance (atom_pos const &pos1,
atom_pos const &pos2)
{
return proxy->position_distance (pos1, pos2);
}
inline cvm::real cvm::position_dist2 (cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
return proxy->position_dist2 (pos1, pos2);
}
inline void cvm::load_atoms (char const *file_name,
std::vector<cvm::atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value)
{
proxy->load_atoms (file_name, atoms, pdb_field, pdb_field_value);
}
inline void cvm::load_coords (char const *file_name,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value)
{
proxy->load_coords (file_name, pos, indices, pdb_field, pdb_field_value);
}
inline void cvm::backup_file (char const *filename)
{
proxy->backup_file (filename);
}
inline cvm::real cvm::rand_gaussian (void)
{
return proxy->rand_gaussian();
}
#endif
// Emacs
// Local Variables:
// mode: C++
// End:
diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp
index 45e687e66..ad0ecf3f0 100644
--- a/src/USER-COLVARS/fix_colvars.cpp
+++ b/src/USER-COLVARS/fix_colvars.cpp
@@ -1,859 +1,859 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Axel Kohlmeyer (TempleU)
------------------------------------------------------------------------- */
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <errno.h>
#include "fix_colvars.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "group.h"
#include "memory.h"
#include "modify.h"
#include "respa.h"
#include "update.h"
#include "colvarproxy_lammps.h"
/* re-usable integer hash table code with static linkage. */
/** hash table top level data structure */
typedef struct inthash_t {
struct inthash_node_t **bucket; /* array of hash nodes */
int size; /* size of the array */
int entries; /* number of entries in table */
int downshift; /* shift cound, used in hash function */
int mask; /* used to select bits for hashing */
} inthash_t;
/** hash table node data structure */
typedef struct inthash_node_t {
int data; /* data in hash node */
int key; /* key for hash lookup */
struct inthash_node_t *next; /* next node in hash chain */
} inthash_node_t;
#define HASH_FAIL -1
#define HASH_LIMIT 0.5
/* initialize new hash table */
static void inthash_init(inthash_t *tptr, int buckets);
/* lookup entry in hash table */
static int inthash_lookup(void *tptr, int key);
/* insert an entry into hash table. */
static int inthash_insert(inthash_t *tptr, int key, int data);
/* delete the hash table */
static void inthash_destroy(inthash_t *tptr);
/* adapted sort for in-place sorting of map indices. */
static void id_sort(int *idmap, int left, int right);
/************************************************************************
* integer hash code:
************************************************************************/
/* inthash() - Hash function returns a hash number for a given key.
* tptr: Pointer to a hash table, key: The key to create a hash number for */
static int inthash(const inthash_t *tptr, int key) {
int hashvalue;
hashvalue = (((key*1103515249)>>tptr->downshift) & tptr->mask);
if (hashvalue < 0) {
hashvalue = 0;
}
return hashvalue;
}
/*
* rebuild_table_int() - Create new hash table when old one fills up.
*
* tptr: Pointer to a hash table
*/
static void rebuild_table_int(inthash_t *tptr) {
inthash_node_t **old_bucket, *old_hash, *tmp;
int old_size, h, i;
old_bucket=tptr->bucket;
old_size=tptr->size;
/* create a new table and rehash old buckets */
inthash_init(tptr, old_size<<1);
for (i=0; i<old_size; i++) {
old_hash=old_bucket[i];
while(old_hash) {
tmp=old_hash;
old_hash=old_hash->next;
h=inthash(tptr, tmp->key);
tmp->next=tptr->bucket[h];
tptr->bucket[h]=tmp;
tptr->entries++;
} /* while */
} /* for */
/* free memory used by old table */
free(old_bucket);
return;
}
/*
* inthash_init() - Initialize a new hash table.
*
* tptr: Pointer to the hash table to initialize
* buckets: The number of initial buckets to create
*/
void inthash_init(inthash_t *tptr, int buckets) {
/* make sure we allocate something */
if (buckets==0)
buckets=16;
/* initialize the table */
tptr->entries=0;
tptr->size=2;
tptr->mask=1;
tptr->downshift=29;
/* ensure buckets is a power of 2 */
while (tptr->size<buckets) {
tptr->size<<=1;
tptr->mask=(tptr->mask<<1)+1;
tptr->downshift--;
} /* while */
/* allocate memory for table */
tptr->bucket=(inthash_node_t **) calloc(tptr->size, sizeof(inthash_node_t *));
return;
}
/*
* inthash_lookup() - Lookup an entry in the hash table and return a pointer to
* it or HASH_FAIL if it wasn't found.
*
* tptr: Pointer to the hash table
* key: The key to lookup
*/
int inthash_lookup(void *ptr, int key) {
const inthash_t *tptr = (const inthash_t *) ptr;
int h;
inthash_node_t *node;
/* find the entry in the hash table */
h=inthash(tptr, key);
for (node=tptr->bucket[h]; node!=NULL; node=node->next) {
if (node->key == key)
break;
}
/* return the entry if it exists, or HASH_FAIL */
return(node ? node->data : HASH_FAIL);
}
/*
* inthash_insert() - Insert an entry into the hash table. If the entry already
* exists return a pointer to it, otherwise return HASH_FAIL.
*
* tptr: A pointer to the hash table
* key: The key to insert into the hash table
* data: A pointer to the data to insert into the hash table
*/
int inthash_insert(inthash_t *tptr, int key, int data) {
int tmp;
inthash_node_t *node;
int h;
/* check to see if the entry exists */
if ((tmp=inthash_lookup(tptr, key)) != HASH_FAIL)
return(tmp);
/* expand the table if needed */
while (tptr->entries>=HASH_LIMIT*tptr->size)
rebuild_table_int(tptr);
/* insert the new entry */
h=inthash(tptr, key);
node=(struct inthash_node_t *) malloc(sizeof(inthash_node_t));
node->data=data;
node->key=key;
node->next=tptr->bucket[h];
tptr->bucket[h]=node;
tptr->entries++;
return HASH_FAIL;
}
/*
* inthash_destroy() - Delete the entire table, and all remaining entries.
*
*/
void inthash_destroy(inthash_t *tptr) {
inthash_node_t *node, *last;
int i;
for (i=0; i<tptr->size; i++) {
node = tptr->bucket[i];
while (node != NULL) {
last = node;
node = node->next;
free(last);
}
}
/* free the entire array of buckets */
if (tptr->bucket != NULL) {
free(tptr->bucket);
memset(tptr, 0, sizeof(inthash_t));
}
}
/************************************************************************
* integer list sort code:
************************************************************************/
/* sort for integer map. initial call id_sort(idmap, 0, natoms - 1); */
static void id_sort(int *idmap, int left, int right)
{
int pivot, l_hold, r_hold;
l_hold = left;
r_hold = right;
pivot = idmap[left];
while (left < right) {
while ((idmap[right] >= pivot) && (left < right))
right--;
if (left != right) {
idmap[left] = idmap[right];
left++;
}
while ((idmap[left] <= pivot) && (left < right))
left++;
if (left != right) {
idmap[right] = idmap[left];
right--;
}
}
idmap[left] = pivot;
pivot = left;
left = l_hold;
right = r_hold;
if (left < pivot)
id_sort(idmap, left, pivot-1);
if (right > pivot)
id_sort(idmap, pivot+1, right);
}
/***************************************************************/
using namespace LAMMPS_NS;
using namespace FixConst;
// initialize static class members
int FixColvars::instances=0;
/***************************************************************
create class and parse arguments in LAMMPS script. Syntax:
fix ID group-ID colvars <config_file> [optional flags...]
optional keyword value pairs:
input <input prefix> (for restarting/continuing, defaults to
NULL, but set to <output prefix> at end)
output <output prefix> (defaults to 'out')
seed <integer> (seed for RNG, defaults to '1966')
tstat <fix label> (label of thermostatting fix)
TODO: add (optional) arguments for RNG seed, temperature compute
***************************************************************/
FixColvars::FixColvars(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg)
{
if (narg < 4)
error->all(FLERR,"Illegal fix colvars command: too few arguments");
if (atom->tag_enable == 0)
error->all(FLERR,"Cannot use fix colvars without atom IDs defined");
if (atom->rmass_flag)
error->all(FLERR,"Cannot use fix colvars for atoms with rmass attribute");
if (instances)
error->all(FLERR,"Only one fix colvars can be active at a time");
++instances;
scalar_flag = 1;
global_freq = 1;
nevery = 1;
extscalar = 1;
me = comm->me;
conf_file = strdup(arg[3]);
rng_seed = 1966;
inp_name = NULL;
out_name = NULL;
tmp_name = NULL;
/* parse optional arguments */
int argsdone = 4;
while (argsdone+1 < narg) {
if (0 == strcmp(arg[argsdone], "input")) {
inp_name = strdup(arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "output")) {
out_name = strdup(arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "seed")) {
rng_seed = atoi(arg[argsdone+1]);
} else if (0 == strcmp(arg[argsdone], "tstat")) {
tmp_name = strdup(arg[argsdone+1]);
} else {
- error->all(FLERR,"Unknown fix imd parameter");
+ error->all(FLERR,"Unknown fix colvars parameter");
}
++argsdone; ++argsdone;
}
if (!out_name) out_name = strdup("out");
/* initialize various state variables. */
tstat_id = -1;
energy = 0.0;
nlevels_respa = 0;
num_coords = 0;
coords = forces = oforce = NULL;
comm_buf = NULL;
force_buf = NULL;
proxy = NULL;
idmap = NULL;
/* storage required to communicate a single coordinate or force. */
size_one = sizeof(struct commdata);
}
/*********************************
* Clean up on deleting the fix. *
*********************************/
FixColvars::~FixColvars()
{
memory->sfree(conf_file);
memory->sfree(inp_name);
memory->sfree(out_name);
memory->sfree(tmp_name);
deallocate();
--instances;
}
/* ---------------------------------------------------------------------- */
void FixColvars::deallocate()
{
memory->destroy(comm_buf);
if (proxy) {
delete proxy;
inthash_t *hashtable = (inthash_t *)idmap;
inthash_destroy(hashtable);
delete hashtable;
}
proxy = NULL;
idmap = NULL;
coords = NULL;
forces = NULL;
oforce = NULL;
comm_buf = NULL;
}
/* ---------------------------------------------------------------------- */
void FixColvars::post_run()
{
deallocate();
memory->sfree(inp_name);
inp_name = strdup(out_name);
}
/* ---------------------------------------------------------------------- */
int FixColvars::setmask()
{
int mask = 0;
mask |= THERMO_ENERGY;
mask |= MIN_POST_FORCE;
mask |= POST_FORCE;
mask |= POST_FORCE_RESPA;
mask |= POST_RUN;
mask |= END_OF_STEP;
return mask;
}
/* ---------------------------------------------------------------------- */
// initial setup of colvars run.
void FixColvars::init()
{
if (strstr(update->integrate_style,"respa"))
nlevels_respa = ((Respa *) update->integrate)->nlevels;
int i,nme,tmp,ndata;
MPI_Status status;
MPI_Request request;
// collect a list of atom type by atom id for the entire system.
// the colvar module requires this information to set masses. :-(
int *typemap,*type_buf;
int nlocal_max,tag_max,max;
const int * const tag = atom->tag;
const int * const type = atom->type;
int nlocal = atom->nlocal;
max=0;
for (i = 0; i < nlocal; i++) max = MAX(max,tag[i]);
MPI_Allreduce(&max,&tag_max,1,MPI_INT,MPI_MAX,world);
MPI_Allreduce(&nlocal,&nlocal_max,1,MPI_INT,MPI_MAX,world);
if (me == 0) {
typemap = new int[tag_max+1];
memset(typemap,0,sizeof(int)*tag_max);
}
type_buf = new int[2*nlocal_max];
if (me == 0) {
for (i=0; i<nlocal; ++i)
typemap[tag[i]] = type[i];
// loop over procs to receive and apply remote data
for (i=1; i < comm->nprocs; ++i) {
MPI_Irecv(type_buf, 2*nlocal_max, MPI_INT, i, 0, world, &request);
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
MPI_Wait(&request, &status);
MPI_Get_count(&status, MPI_INT, &ndata);
for (int k=0; k<ndata; k+=2)
typemap[type_buf[k]] = type_buf[k+1];
}
} else { // me != 0
// copy tag/type data into communication buffer
nme = 0;
for (i=0; i<nlocal; ++i) {
type_buf[nme] = tag[i];
type_buf[nme+1] = type[i];
nme +=2;
}
/* blocking receive to wait until it is our turn to send data. */
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
MPI_Rsend(type_buf, nme, MPI_INT, 0, 0, world);
}
// now create and initialize the colvar proxy
if (me == 0) {
if (inp_name) {
if (strcmp(inp_name,"NULL") == 0)
memory->sfree(inp_name);
inp_name = NULL;
}
double t_target = 0.0;
if (tmp_name) {
if (strcmp(tmp_name,"NULL") == 0)
tstat_id = -1;
else {
tstat_id = modify->find_fix(tmp_name);
if (tstat_id < 0) error->one(FLERR,"Could not find tstat fix ID");
double *tt = (double*)modify->fix[tstat_id]->extract("t_target",tmp);
if (tt) t_target = *tt;
}
}
proxy = new colvarproxy_lammps(lmp,conf_file,inp_name,out_name,
rng_seed,typemap);
coords = proxy->get_coords();
forces = proxy->get_forces();
oforce = proxy->get_oforce();
num_coords = coords->size();
}
// send the list of all colvar atom IDs to all nodes.
// also initialize and build hashtable on master.
MPI_Bcast(&num_coords, 1, MPI_INT, 0, world);
memory->create(taglist,num_coords,"colvars:taglist");
memory->create(force_buf,3*num_coords,"colvars:force_buf");
if (me == 0) {
std::vector<int> *tags_list = proxy->get_tags();
std::vector<int> &tl = *tags_list;
inthash_t *hashtable=new inthash_t;
inthash_init(hashtable, num_coords);
idmap = (void *)hashtable;
for (i=0; i < num_coords; ++i) {
taglist[i] = tl[i];
inthash_insert(hashtable, tl[i], i);
}
}
MPI_Bcast(taglist, num_coords, MPI_INT, 0, world);
// determine size of comm buffer
nme=0;
for (i=0; i < num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal))
++nme;
}
MPI_Allreduce(&nme,&nmax,1,MPI_INT,MPI_MAX,world);
memory->create(comm_buf,nmax,"colvars:comm_buf");
const double * const * const x = atom->x;
if (me == 0) {
std::vector<struct commdata> &cd = *coords;
std::vector<struct commdata> &of = *oforce;
// store coordinate data in holding array, clear old forces
for (i=0; i<num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal)) {
of[i].tag = cd[i].tag = tag[k];
of[i].type = cd[i].type = type[k];
cd[i].x = x[k][0];
cd[i].y = x[k][1];
cd[i].z = x[k][2];
of[i].x = of[i].y = of[i].z = 0.0;
}
}
// loop over procs to receive and apply remote data
for (i=1; i < comm->nprocs; ++i) {
int maxbuf = nmax*size_one;
MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request);
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
MPI_Wait(&request, &status);
MPI_Get_count(&status, MPI_BYTE, &ndata);
ndata /= size_one;
for (int k=0; k<ndata; ++k) {
const int j = inthash_lookup(idmap, comm_buf[k].tag);
if (j != HASH_FAIL) {
of[j].tag = cd[j].tag = comm_buf[k].tag;
of[j].type = cd[j].type = comm_buf[k].type;
cd[j].x = comm_buf[k].x;
cd[j].y = comm_buf[k].y;
cd[j].z = comm_buf[k].z;
of[j].x = of[j].y = of[j].z = 0.0;
}
}
}
} else { // me != 0
// copy coordinate data into communication buffer
nme = 0;
for (i=0; i<num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal)) {
comm_buf[nme].tag = tag[k];
comm_buf[nme].type = type[k];
comm_buf[nme].x = x[k][0];
comm_buf[nme].y = x[k][1];
comm_buf[nme].z = x[k][2];
++nme;
}
}
/* blocking receive to wait until it is our turn to send data. */
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
}
// clear temporary storage
if (me == 0) delete typemap;
delete type_buf;
}
/* ---------------------------------------------------------------------- */
void FixColvars::setup(int vflag)
{
if (strstr(update->integrate_style,"verlet"))
post_force(vflag);
else
post_force_respa(vflag,0,0);
}
/* ---------------------------------------------------------------------- */
/* Main colvars handler:
* Send coodinates and add colvar forces to atoms. */
void FixColvars::post_force(int vflag)
{
// some housekeeping: update status of the proxy as needed.
if (me == 0) {
if (tstat_id < 0) {
proxy->set_temperature(0.0);
} else {
int tmp;
// get thermostat target temperature from corresponding fix,
// if the fix supports extraction.
double *tt = (double *) modify->fix[tstat_id]->extract("t_target",tmp);
if (tt)
proxy->set_temperature(*tt);
else
proxy->set_temperature(0.0);
}
}
const int * const tag = atom->tag;
const double * const * const x = atom->x;
double * const * const f = atom->f;
const int nlocal = atom->nlocal;
/* check and potentially grow local communication buffers. */
int i,nmax_new,nme=0;
for (i=0; i < num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal))
++nme;
}
MPI_Allreduce(&nme,&nmax_new,1,MPI_INT,MPI_MAX,world);
if (nmax_new > nmax) {
nmax = nmax_new;
memory->grow(comm_buf,nmax,"colvars:comm_buf");
}
MPI_Status status;
MPI_Request request;
int tmp, ndata;
if (me == 0) {
std::vector<struct commdata> &cd = *coords;
// store coordinate data
for (i=0; i<num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal)) {
cd[i].x = x[k][0];
cd[i].y = x[k][1];
cd[i].z = x[k][2];
}
}
/* loop over procs to receive remote data */
for (i=1; i < comm->nprocs; ++i) {
int maxbuf = nmax*size_one;
MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request);
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
MPI_Wait(&request, &status);
MPI_Get_count(&status, MPI_BYTE, &ndata);
ndata /= size_one;
for (int k=0; k<ndata; ++k) {
const int j = inthash_lookup(idmap, comm_buf[k].tag);
if (j != HASH_FAIL) {
cd[j].x = comm_buf[k].x;
cd[j].y = comm_buf[k].y;
cd[j].z = comm_buf[k].z;
}
}
}
} else { // me != 0
/* copy coordinate data into communication buffer */
nme = 0;
for (i=0; i<num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal)) {
comm_buf[nme].tag = tag[k];
comm_buf[nme].x = x[k][0];
comm_buf[nme].y = x[k][1];
comm_buf[nme].z = x[k][2];
++nme;
}
}
/* blocking receive to wait until it is our turn to send data. */
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
}
////////////////////////////////////////////////////////////////////////
// call our workhorse and retrieve additional information.
if (me == 0) {
energy = proxy->compute();
store_forces = proxy->need_system_forces();
}
////////////////////////////////////////////////////////////////////////
// broadcast store_forces flag and energy data to all processors
MPI_Bcast(&energy, 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&store_forces, 1, MPI_INT, 0, world);
// broadcast and apply biasing forces
if (me == 0) {
std::vector<struct commdata> &fo = *forces;
double *fbuf = force_buf;
for (int j=0; j < num_coords; ++j) {
*fbuf++ = fo[j].x;
*fbuf++ = fo[j].y;
*fbuf++ = fo[j].z;
}
}
MPI_Bcast(force_buf, 3*num_coords, MPI_DOUBLE, 0, world);
for (int i=0; i < num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal)) {
f[k][0] += force_buf[3*i+0];
f[k][1] += force_buf[3*i+1];
f[k][2] += force_buf[3*i+2];
}
}
}
/* ---------------------------------------------------------------------- */
void FixColvars::min_post_force(int vflag)
{
post_force(vflag);
}
/* ---------------------------------------------------------------------- */
void FixColvars::post_force_respa(int vflag, int ilevel, int iloop)
{
/* only process colvar forces on the outmost RESPA level. */
if (ilevel == nlevels_respa-1) post_force(vflag);
return;
}
/* ---------------------------------------------------------------------- */
void FixColvars::end_of_step()
{
if (store_forces) {
const int * const tag = atom->tag;
double * const * const f = atom->f;
const int nlocal = atom->nlocal;
/* check and potentially grow local communication buffers. */
int i,nmax_new,nme=0;
for (i=0; i < num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal))
++nme;
}
MPI_Allreduce(&nme,&nmax_new,1,MPI_INT,MPI_MAX,world);
if (nmax_new > nmax) {
nmax = nmax_new;
memory->grow(comm_buf,nmax,"colvars:comm_buf");
}
MPI_Status status;
MPI_Request request;
int tmp, ndata;
if (me == 0) {
// store old force data
std::vector<struct commdata> &of = *oforce;
for (i=0; i<num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal)) {
const int j = inthash_lookup(idmap, tag[i]);
if (j != HASH_FAIL) {
of[j].x = f[k][0];
of[j].y = f[k][1];
of[j].z = f[k][2];
}
}
}
/* loop over procs to receive remote data */
for (i=1; i < comm->nprocs; ++i) {
int maxbuf = nmax*size_one;
MPI_Irecv(comm_buf, maxbuf, MPI_BYTE, i, 0, world, &request);
MPI_Send(&tmp, 0, MPI_INT, i, 0, world);
MPI_Wait(&request, &status);
MPI_Get_count(&status, MPI_BYTE, &ndata);
ndata /= size_one;
for (int k=0; k<ndata; ++k) {
const int j = inthash_lookup(idmap, comm_buf[k].tag);
if (j != HASH_FAIL) {
of[j].x = comm_buf[k].x;
of[j].y = comm_buf[k].y;
of[j].z = comm_buf[k].z;
}
}
}
} else { // me != 0
/* copy total force data into communication buffer */
nme = 0;
for (i=0; i<num_coords; ++i) {
const int k = atom->map(taglist[i]);
if ((k >= 0) && (k < nlocal)) {
comm_buf[nme].tag = tag[k];
comm_buf[nme].x = f[k][0];
comm_buf[nme].y = f[k][1];
comm_buf[nme].z = f[k][2];
++nme;
}
}
/* blocking receive to wait until it is our turn to send data. */
MPI_Recv(&tmp, 0, MPI_INT, 0, 0, world, &status);
MPI_Rsend(comm_buf, nme*size_one, MPI_BYTE, 0, 0, world);
}
}
}
/* ---------------------------------------------------------------------- */
double FixColvars::compute_scalar()
{
return energy;
}
/* ---------------------------------------------------------------------- */
/* local memory usage. approximately. */
double FixColvars::memory_usage(void)
{
double bytes = (double) (num_coords * (2*sizeof(int)+3*sizeof(double)));
bytes += (double) (nmax*size_one) + sizeof(this);
return bytes;
}