Page Menu
Home
c4science
Search
Configure Global Search
Log In
Files
F73813433
colvarbias_meta.cpp
No One
Temporary
Actions
Download File
Edit File
Delete File
View Transforms
Subscribe
Mute Notifications
Award Token
Subscribers
None
File Metadata
Details
File Info
Storage
Attached
Created
Wed, Jul 24, 15:27
Size
59 KB
Mime Type
text/x-c
Expires
Fri, Jul 26, 15:27 (2 d)
Engine
blob
Format
Raw Data
Handle
19271807
Attached To
rLAMMPS lammps
colvarbias_meta.cpp
View Options
/// -*- c++ -*-
#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(),
new_hills_begin(hills.end()),
state_file_step(0)
{
}
colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key)
: colvarbias(conf, key),
new_hills_begin(hills.end()),
state_file_step(0)
{
if (cvm::n_abf_biases > 0)
cvm::log("Warning: running ABF and metadynamics together is not recommended unless applyBias is off for ABF.\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);
{
bool b_replicas = false;
get_keyval(conf, "multipleReplicas", b_replicas, false);
if (b_replicas)
comm = multiple_replicas;
else
comm = single_replica;
}
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;
size_t i;
for (i = 0; i < colvars.size(); i++) {
if (colvars[i]->expand_boundaries) {
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);
if (! get_keyval(conf, "writeFreeEnergyFile", dump_fes, true))
get_keyval(conf, "dumpFreeEnergyFile", dump_fes, true, colvarparse::parse_silent);
get_keyval(conf, "saveFreeEnergyFile", dump_fes_save, false);
for (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;
hills_energy = NULL;
hills_energy_gradients = NULL;
}
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::error("Error: replicaID must be defined "
"when using more than one replica.\n", INPUT_ERROR);
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 of input/output and slow down your calculations. "
"Please consider disabling it.\n");
}
get_keyval(conf, "writeHillsTrajectory", b_hills_traj, false);
// for well-tempered metadynamics
get_keyval(conf, "wellTempered", well_tempered, false);
get_keyval(conf, "biasTemperature", bias_temperature, -1.0);
if ((bias_temperature == -1.0) && well_tempered) {
cvm::fatal_error("Error: biasTemperature is not set.\n");
}
if (well_tempered) {
cvm::log("Well-tempered metadynamics is used.\n");
cvm::log("The bias temperature is "+cvm::to_str(bias_temperature)+".\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.is_open())
replica_hills_os.close();
if (hills_traj_os.is_open())
hills_traj_os.close();
if (cvm::n_meta_biases > 0)
cvm::n_meta_biases -= 1;
}
// **********************************************************************
// 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, true);
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.is_open()) {
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.empty()) {
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.is_open()) {
// 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 (!colvars[i]->hard_lower_boundary)
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 (!colvars[i]->hard_upper_boundary)
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->setup(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->setup(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:
if (well_tempered) {
cvm::real hills_energy_sum_here = 0.0;
if (use_grids) {
std::vector<int> curr_bin = hills_energy->get_colvars_index();
hills_energy_sum_here = hills_energy->value(curr_bin);
} else {
calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
}
cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
create_hill(hill((hill_weight*exp_weight), colvars, hill_width));
} else {
create_hill(hill(hill_weight, colvars, hill_width));
}
break;
case multiple_replicas:
if (well_tempered) {
cvm::real hills_energy_sum_here = 0.0;
if (use_grids) {
std::vector<int> curr_bin = hills_energy->get_colvars_index();
hills_energy_sum_here = hills_energy->value(curr_bin);
} else {
calc_hills(new_hills_begin, hills.end(), hills_energy_sum_here);
}
cvm::real const exp_weight = std::exp(-1.0*hills_energy_sum_here/(bias_temperature*cvm::boltzmann()));
create_hill(hill((hill_weight*exp_weight), colvars, hill_width, replica_id));
} else {
create_hill(hill(hill_weight, colvars, hill_width, replica_id));
}
if (replica_hills_os.is_open()) {
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
if (replica_hills_os.is_open())
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]->value());
}
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 const 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)
hill_iter h;
switch (colvars[i]->value().type()) {
case colvarvalue::type_scalar:
for (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_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
for (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:
case colvarvalue::type_quaternionderiv:
for (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_vector:
for (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].vector1d_value +=
( h->weight() * h->value() * (0.5 / (half_width*half_width)) *
(colvars[i]->dist2_lgrad(x, center)).vector1d_value );
}
break;
case colvarvalue::type_notset:
case colvarvalue::type_all:
default:
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,
bool print_progress)
{
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);
size_t count = 0;
size_t const print_frequency = ((hills.size() >= 1000000) ? 1 : (1000000/(hills.size()+1)));
if (hg != NULL) {
// loop over the points of the grid
for ( ;
(he->index_ok(he_ix)) && (hg->index_ok(hg_ix));
count++) {
size_t i;
for (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 (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);
if ((count % print_frequency) == 0) {
if (print_progress) {
cvm::real const progress = cvm::real(count) / cvm::real(hg->number_of_points());
std::ostringstream os;
os.setf(std::ios::fixed, std::ios::floatfield);
os << std::setw(6) << std::setprecision(2)
<< 100.0 * progress
<< "% done.";
cvm::log(os.str());
}
}
}
} 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);
count++;
if ((count % print_frequency) == 0) {
if (print_progress) {
cvm::real const progress = cvm::real(count) / cvm::real(he->number_of_points());
std::ostringstream os;
os.setf(std::ios::fixed, std::ios::floatfield);
os << std::setw(6) << std::setprecision(2)
<< 100.0 * progress
<< "% done.";
cvm::log(os.str());
}
}
}
}
if (print_progress) {
cvm::log("100.00% done.");
}
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, true);
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.is_open()) {
replicas_registry.clear();
while (colvarparse::getline_nocomments(reg_file, line))
replicas_registry.append(line+"\n");
} else {
cvm::error("Error: failed to open file \""+replicas_registry_file+
"\" for reading.\n", FILE_ERROR);
}
}
// 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+"\""+
": 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);
}
}
}
} 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.is_open()) {
// try to resume the previous position
is.seekg((replicas[ir])->replica_hills_file_pos, std::ios::beg);
if (!is.is_open()){
// 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())
cvm::log("Metadynamics bias \""+this->name+"\""+
": 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::log("Failed to read the file \""+(replicas[ir])->replica_hills_file+
"\": will try again in "+
cvm::to_str(replica_update_freq)+" steps.\n");
(replicas[ir])->update_status++;
// 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+"\""+
" failed to read completely the output of replica \""+
(replicas[ir])->replica_id+
"\" after more than "+
cvm::to_str((replicas[ir])->update_status * replica_update_freq)+
" steps. Ensure that it is still running.\n");
}
}
}
// **********************************************************************
// input functions
// **********************************************************************
std::istream & colvarbias_meta::read_restart(std::istream& is)
{
size_t const start_pos = is.tellg();
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+"\""+
".\n");
}
std::string key, brace, conf;
if ( !(is >> key) || !(key == "metadynamics") ||
!(is >> brace) || !(brace == "{") ||
!(is >> colvarparse::read_block("configuration", conf)) ) {
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 useGrids was off when the state file was written, "
"enable rebinGrids now to regenerate the grids.\n");
else {
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 (!rebin_grids) {
if (hills_energy_backup == NULL)
cvm::fatal_error("Error: couldn't read the free energy gradients grid for metadynamics bias \""+
this->name+"\""+
((comm != single_replica) ? ", replica \""+replica_id+"\"" : "")+
"; if useGrids was off when the state file was written, "
"enable rebinGrids now to regenerate the grids.\n");
else {
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.empty();
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.empty())
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.empty())) {
// 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, true);
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.empty())
recount_hills_off_grid(hills.begin(), hills.end(), hills_energy);
}
if (use_grids) {
if (!hills_off_grid.empty()) {
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());
}
{
// 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, true);
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
// **********************************************************************
int colvarbias_meta::setup_output()
{
if (comm == multiple_replicas) {
// 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 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.is_open()) { // 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();
}
// 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
if (!replica_hills_os.is_open()) {
cvm::backup_file(replica_hills_file.c_str());
replica_hills_os.open(replica_hills_file.c_str());
if (!replica_hills_os.is_open())
cvm::error("Error: in opening file \""+
replica_hills_file+"\" for writing.\n", FILE_ERROR);
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.is_open())
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.is_open())
cvm::error("Error: in opening file \""+
replicas_registry_file+"\" for writing.\n", FILE_ERROR);
reg_os << replica_id << " " << replica_list_file << "\n";
reg_os.close();
}
}
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.is_open())
cvm::error("Error: in opening hills output file \"" +
traj_file_name+"\".\n", FILE_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
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(*hills_energy);
pmf->setup();
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 (well_tempered) {
cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
pmf->multiply_constant(well_temper_scale);
}
{
std::string const fes_file_name(fes_file_name_prefix +
((comm != single_replica) ? ".partial" : "") +
(dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file(fes_file_name.c_str());
cvm::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 (well_tempered) {
cvm::real const well_temper_scale = (bias_temperature + cvm::temperature()) / bias_temperature;
pmf->multiply_constant(well_temper_scale);
}
std::string const fes_file_name(fes_file_name_prefix +
(dump_fes_save ?
"."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file(fes_file_name.c_str());
cvm::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());
cvm::ofstream rep_state_os(replica_state_file.c_str());
if (!rep_state_os.is_open())
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();
cvm::backup_file(replica_hills_file.c_str());
replica_hills_os.open(replica_hills_file.c_str());
if (!replica_hills_os.is_open())
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);
size_t i;
os << " ";
for (i = 0; i < centers.size(); i++) {
os << " ";
os << std::setprecision(cvm::cv_prec)
<< std::setw(cvm::cv_width) << centers[i];
}
os << " ";
for (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";
size_t i;
os << " centers ";
for (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 (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;
}
Event Timeline
Log In to Comment