ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@13126 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2015-02-19 18:20:52 +00:00
parent f300f2367e
commit 067833cf07
21 changed files with 597 additions and 360 deletions
Download Patch File Download Diff File Expand all files Collapse all files

8
lib/colvars/colvar.cpp
View File

@ -177,7 +177,7 @@ colvar::colvar(std::string const &conf)
cvm::error("Could not parse scripted colvar type.");
return;
}
x_reported.type (x.type());
x_reported.type(x.type());
cvm::log(std::string("Expecting colvar value of type ")
+ colvarvalue::type_desc(x.type()));
@ -987,7 +987,9 @@ void colvar::calc()
}
}
if (tasks[task_system_force]) {
if (tasks[task_system_force] && !tasks[task_extended_lagrangian]) {
// If extended Lagrangian is enabled, system force calculation is trivial
// and done together with integration of the extended coordinate.
if (tasks[task_scripted]) {
// TODO see if this could reasonably be done in a generic way
@ -1594,7 +1596,7 @@ int colvar::write_output_files()
if (acf.size()) {
cvm::log("Writing acf to file \""+acf_outfile+"\".\n");
std::ofstream acf_os(acf_outfile.c_str());
cvm::ofstream acf_os(acf_outfile.c_str());
if (! acf_os.good()) {
cvm::error("Cannot open file \""+acf_outfile+"\".\n", FILE_ERROR);
}

2
lib/colvars/colvar.h
View File

@ -457,7 +457,7 @@ protected:
/// Timesteps to skip between two values in the running average series
size_t runave_stride;
/// Name of the file to write the running average
std::ofstream runave_os;
cvm::ofstream runave_os;
/// Current value of the running average
colvarvalue runave;
/// Current value of the square deviation from the running average

9
lib/colvars/colvaratoms.cpp
View File

@ -315,21 +315,24 @@ int cvm::atom_group::parse(std::string const &conf,
#endif
if (!b_dummy) {
// calculate total mass (TODO: this is the step that most needs deferred re-initialization)
this->total_mass = 0.0;
for (cvm::atom_iter ai = this->begin();
ai != this->end(); ai++) {
this->total_mass += ai->mass;
}
}
if (!b_dummy) {
// whether these atoms will ever receive forces or not
bool enable_forces = true;
// disableForces is deprecated
if (get_keyval(group_conf, "enableForces", enable_forces, true, mode)) {
noforce = !enable_forces;
} else {
get_keyval(group_conf, "disableForces", noforce, false, mode);
get_keyval(group_conf, "disableForces", noforce, false, colvarparse::parse_silent);
}
get_keyval(group_conf, "weights", weights, weights, colvarparse::parse_silent);
}
// FITTING OPTIONS

3
lib/colvars/colvaratoms.h
View File

@ -139,6 +139,9 @@ public:
/// Allocates and populates the sorted list of atom ids
int create_sorted_ids(void);
/// List of user-defined weights to be used by certain CVCs
std::vector<cvm::real> weights;
/// \brief When updating atomic coordinates, translate them to align with the
/// center of mass of the reference coordinates

115
lib/colvars/colvarbias_abf.cpp
View File

@ -76,6 +76,8 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
if (!colvars[i]->tasks[colvar::task_extended_lagrangian]) {
// request computation of Jacobian force
// ultimately, will be done regardless of extended Lagrangian
// and colvar should then just report zero Jacobian force
colvars[i]->enable(colvar::task_Jacobian_force);
// request Jacobian force as part as system force
@ -116,7 +118,11 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
force = new cvm::real [colvars.size()];
// Construct empty grids based on the colvars
samples = new colvar_grid_count (colvars);
if (cvm::debug()) {
cvm::log("Allocating count and free energy gradient grids.\n");
}
samples = new colvar_grid_count(colvars);
gradients = new colvar_grid_gradient(colvars);
gradients->samples = samples;
samples->has_parent_data = true;
@ -124,7 +130,7 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key)
// For shared ABF, we store a second set of grids.
// This used to be only if "shared" was defined,
// but now we allow calling share externally (e.g. from Tcl).
last_samples = new colvar_grid_count (colvars);
last_samples = new colvar_grid_count(colvars);
last_gradients = new colvar_grid_gradient(colvars);
last_gradients->samples = last_samples;
last_samples->has_parent_data = true;
@ -260,9 +266,9 @@ cvm::real colvarbias_abf::update()
}
if (b_history_files && (cvm::step_absolute() % history_freq) == 0) {
cvm::log("ABFHISTORYFILE "+cvm::to_str(cvm::step_absolute()));
// append to existing file only if cvm::step_absolute() > 0
// file already exists iff cvm::step_relative() > 0
// otherwise, backup and replace
write_gradients_samples(output_prefix + ".hist", (cvm::step_absolute() > 0));
write_gradients_samples(output_prefix + ".hist", (cvm::step_relative() > 0));
}
if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) {
@ -365,28 +371,32 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::string gradients_out_name = prefix + ".grad";
std::ios::openmode mode = (append ? std::ios::app : std::ios::out);
std::ofstream samples_os;
std::ofstream gradients_os;
cvm::ofstream samples_os;
cvm::ofstream gradients_os;
if (!append) cvm::backup_file(samples_out_name.c_str());
samples_os.open(samples_out_name.c_str(), mode);
if (!samples_os.good()) cvm::error("Error opening ABF samples file " + samples_out_name + " for writing");
if (!samples_os.is_open()) {
cvm::error("Error opening ABF samples file " + samples_out_name + " for writing");
}
samples->write_multicol(samples_os);
samples_os.close();
if (!append) cvm::backup_file(gradients_out_name.c_str());
gradients_os.open(gradients_out_name.c_str(), mode);
if (!gradients_os.good()) cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing");
if (!gradients_os.is_open()) {
cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing");
}
gradients->write_multicol(gradients_os);
gradients_os.close();
if (colvars.size() == 1) {
std::string pmf_out_name = prefix + ".pmf";
if (!append) cvm::backup_file(pmf_out_name.c_str());
std::ofstream pmf_os;
cvm::ofstream pmf_os;
// Do numerical integration and output a PMF
pmf_os.open(pmf_out_name.c_str(), mode);
if (!pmf_os.good()) cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
if (!pmf_os.is_open()) cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
gradients->write_1D_integral(pmf_os);
pmf_os << std::endl;
pmf_os.close();
@ -428,13 +438,13 @@ void colvarbias_abf::read_gradients_samples()
cvm::log("Reading sample count from " + samples_in_name + " and gradients from " + gradients_in_name);
is.open(samples_in_name.c_str());
if (!is.good()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");
if (!is.is_open()) cvm::error("Error opening ABF samples file " + samples_in_name + " for reading");
samples->read_multicol(is, true);
is.close();
is.clear();
is.open(gradients_in_name.c_str());
if (!is.good()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading");
if (!is.is_open()) cvm::error("Error opening ABF gradient file " + gradients_in_name + " for reading");
gradients->read_multicol(is, true);
is.close();
}
@ -550,13 +560,41 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
: colvarbias(conf, key),
grid(NULL), out_name("")
{
get_keyval(conf, "outputfreq", output_freq, cvm::restart_out_freq);
get_keyval(conf, "outputFreq", output_freq, cvm::restart_out_freq);
/// with VMD, this may not be an error
// if ( output_freq == 0 ) {
// cvm::error("User required histogram with zero output frequency");
// }
if ( output_freq == 0 ) {
cvm::error("User required histogram with zero output frequency");
{
colvar_array_size = 1;
bool colvar_array = false;
if (get_keyval(conf, "sumVectorColvars", colvar_array, colvar_array)) {
size_t i;
for (i = 0; i < colvars.size(); i++) {
if (colvars[i]->value().type() == colvarvalue::type_vector) {
if (colvar_array_size == 1) {
colvar_array_size = colvars[i]->value().size();
} else {
if (colvar_array_size != colvars[i]->value().size()) {
cvm::error("Error: trying to combine vector colvars of different lengths.\n", INPUT_ERROR);
}
}
}
}
}
}
grid = new colvar_grid_count();
{
std::string grid_conf;
if (key_lookup(conf, "grid", grid_conf)) {
grid->parse_params(grid_conf);
} else {
grid->init_from_colvars(colvars);
}
}
grid = new colvar_grid_count (colvars);
bin.assign(colvars.size(), 0);
cvm::log("Finished histogram setup.\n");
@ -565,7 +603,8 @@ colvarbias_histogram::colvarbias_histogram(std::string const &conf, char const *
/// Destructor
colvarbias_histogram::~colvarbias_histogram()
{
if (grid_os.good()) grid_os.close();
if (grid_os.is_open())
grid_os.close();
if (grid) {
delete grid;
@ -579,7 +618,10 @@ colvarbias_histogram::~colvarbias_histogram()
/// Update the grid
cvm::real colvarbias_histogram::update()
{
if (cvm::debug()) cvm::log("Updating Grid bias " + this->name);
if (cvm::debug()) {
cvm::log("Updating histogram bias " + this->name);
}
// At the first timestep, we need to assign out_name since
// output_prefix is unset during the constructor
@ -589,22 +631,49 @@ cvm::real colvarbias_histogram::update()
cvm::log("Histogram " + this->name + " will be written to file \"" + out_name + "\"");
}
for (size_t i=0; i<colvars.size(); i++) {
bin[i] = grid->current_bin_scalar(i);
bin.assign(colvars.size(), 0);
{
// update indices for all scalar values
size_t i;
for (i = 0; i < colvars.size(); i++) {
if (colvars[i]->value().type() == colvarvalue::type_scalar) {
bin[i] = grid->value_to_bin_scalar(colvars[i]->value(), i);
}
}
}
if ( grid->index_ok(bin) ) { // Only within bounds of the grid...
grid->incr_count(bin);
if (colvar_array_size > 1) {
// update indices for all vector/array values
size_t iv, i;
for (iv = 0; iv < colvar_array_size; iv++) {
for (i = 0; i < colvars.size(); i++) {
if (colvars[i]->value().type() == colvarvalue::type_vector) {
bin[i] = grid->value_to_bin_scalar(colvars[i]->value().vector1d_value[iv], i);
}
}
if (grid->index_ok(bin)) {
// Only within bounds of the grid...
grid->incr_count(bin);
}
}
} else {
if (grid->index_ok(bin)) {
// Only within bounds of the grid...
grid->incr_count(bin);
}
}
if (output_freq && (cvm::step_absolute() % output_freq) == 0) {
if (cvm::debug()) cvm::log("Histogram bias trying to write grid to disk");
grid_os.open(out_name.c_str());
if (!grid_os.good()) cvm::error("Error opening histogram file " + out_name + " for writing");
if (!grid_os.is_open()) cvm::error("Error opening histogram file " + out_name + " for writing");
grid->write_multicol(grid_os);
grid_os.close();
}
return 0.0; // no bias energy for histogram
}

9
lib/colvars/colvarbias_abf.h
View File

@ -1,3 +1,4 @@
// -*- c++ -*-
/************************************************************************
* Headers for the ABF and histogram biases *
************************************************************************/
@ -100,7 +101,7 @@ public:
cvm::real update();
private:
protected:
/// n-dim histogram
colvar_grid_count *grid;
@ -108,8 +109,12 @@ private:
std::string out_name;
int output_freq;
/// If one or more of the variables are \link type_vector \endlink, treat them as arrays of this length
size_t colvar_array_size;
void write_grid();
std::ofstream grid_os; /// Stream for writing grid to disk
cvm::ofstream grid_os; /// Stream for writing grid to disk
std::istream& read_restart(std::istream&);
std::ostream& write_restart(std::ostream&);

18
lib/colvars/colvarbias_meta.cpp
View File

@ -86,7 +86,7 @@ colvarbias_meta::colvarbias_meta(std::string const &conf, char const *key)
colvars[i]->enable(colvar::task_grid);
}
hills_energy = new colvar_grid_scalar (colvars);
hills_energy = new colvar_grid_scalar(colvars);
hills_energy_gradients = new colvar_grid_gradient(colvars);
} else {
rebin_grids = false;
@ -430,11 +430,11 @@ cvm::real colvarbias_meta::update()
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->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->create(new_sizes, 0.0, colvars.size());
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);
@ -746,6 +746,8 @@ void colvarbias_meta::calc_hills_force(size_t const &i,
break;
case colvarvalue::type_notset:
case colvarvalue::type_all:
default:
break;
}
}
@ -959,7 +961,7 @@ void colvarbias_meta::update_replicas_registry()
(replicas.back())->comm = multiple_replicas;
if (use_grids) {
(replicas.back())->hills_energy = new colvar_grid_scalar (colvars);
(replicas.back())->hills_energy = new colvar_grid_scalar(colvars);
(replicas.back())->hills_energy_gradients = new colvar_grid_gradient(colvars);
}
}
@ -1552,7 +1554,7 @@ void colvarbias_meta::write_pmf()
{
// allocate a new grid to store the pmf
colvar_grid_scalar *pmf = new colvar_grid_scalar(*hills_energy);
pmf->create();
pmf->setup();
std::string fes_file_name_prefix(cvm::output_prefix);
@ -1582,7 +1584,7 @@ void colvarbias_meta::write_pmf()
"."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file(fes_file_name.c_str());
std::ofstream fes_os(fes_file_name.c_str());
cvm::ofstream fes_os(fes_file_name.c_str());
pmf->write_multicol(fes_os);
fes_os.close();
}
@ -1606,7 +1608,7 @@ void colvarbias_meta::write_pmf()
"."+cvm::to_str(cvm::step_absolute()) : "") +
".pmf");
cvm::backup_file(fes_file_name.c_str());
std::ofstream fes_os(fes_file_name.c_str());
cvm::ofstream fes_os(fes_file_name.c_str());
pmf->write_multicol(fes_os);
fes_os.close();
}
@ -1621,7 +1623,7 @@ 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());
cvm::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");

4
lib/colvars/colvarbias_meta.h
View File

@ -61,7 +61,7 @@ protected:
/// Write the hill logfile
bool b_hills_traj;
/// Logfile of hill management (creation and deletion)
std::ofstream hills_traj_os;
cvm::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
@ -223,7 +223,7 @@ protected:
std::string replica_hills_file;
/// \brief Output stream corresponding to replica_hills_file
std::ofstream replica_hills_os;
cvm::ofstream replica_hills_os;
/// Position within replica_hills_file (when reading it)
int replica_hills_file_pos;

4
lib/colvars/colvarcomp.h
View File

@ -815,8 +815,6 @@ class colvar::h_bond
: public colvar::cvc
{
protected:
/// Atoms involved in the component
cvm::atom acceptor, donor;
/// \brief "Cutoff" distance between acceptor and donor
cvm::real r0;
/// Integer exponent of the function numerator
@ -1141,6 +1139,8 @@ class colvar::cartesian
protected:
/// Atom group
cvm::atom_group atoms;
/// Which Cartesian coordinates to include
std::vector<size_t> axes;
public:
cartesian(std::string const &conf);
cartesian();

9
lib/colvars/colvarcomp_angles.cpp
View File

@ -54,10 +54,6 @@ colvar::angle::angle()
void colvar::angle::calc_value()
{
group1.read_positions();
group2.read_positions();
group3.read_positions();
cvm::atom_pos const g1_pos = group1.center_of_mass();
cvm::atom_pos const g2_pos = group2.center_of_mass();
cvm::atom_pos const g3_pos = group3.center_of_mass();
@ -215,11 +211,6 @@ colvar::dihedral::dihedral()
void colvar::dihedral::calc_value()
{
group1.read_positions();
group2.read_positions();
group3.read_positions();
group4.read_positions();
cvm::atom_pos const g1_pos = group1.center_of_mass();
cvm::atom_pos const g2_pos = group2.center_of_mass();
cvm::atom_pos const g3_pos = group3.center_of_mass();

79
lib/colvars/colvarcomp_coordnums.cpp
View File

@ -62,8 +62,8 @@ cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
if (calculate_gradients) {
cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) - func*ed2*(xd/l2))*(-1.0);
cvm::rvector const dl2dx((2.0/(r0_vec.x*r0_vec.x))*diff.x,
(2.0/(r0_vec.y*r0_vec.y))*diff.y,
(2.0/(r0_vec.z*r0_vec.z))*diff.z);
(2.0/(r0_vec.y*r0_vec.y))*diff.y,
(2.0/(r0_vec.z*r0_vec.z))*diff.z);
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
}
@ -134,10 +134,10 @@ void colvar::coordnum::calc_value()
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
x.real_value += switching_function<false> (r0_vec, en, ed, *ai1, group2_com_atom);
x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, group2_com_atom);
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
x.real_value += switching_function<false> (r0, en, ed, *ai1, group2_com_atom);
x.real_value += switching_function<false>(r0, en, ed, *ai1, group2_com_atom);
}
} else {
@ -145,12 +145,12 @@ void colvar::coordnum::calc_value()
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
x.real_value += switching_function<false> (r0_vec, en, ed, *ai1, *ai2);
x.real_value += switching_function<false>(r0_vec, en, ed, *ai1, *ai2);
}
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
x.real_value += switching_function<false> (r0, en, ed, *ai1, *ai2);
x.real_value += switching_function<false>(r0, en, ed, *ai1, *ai2);
}
}
}
@ -168,10 +168,10 @@ void colvar::coordnum::calc_gradients()
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
switching_function<true> (r0_vec, en, ed, *ai1, group2_com_atom);
switching_function<true>(r0_vec, en, ed, *ai1, group2_com_atom);
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
switching_function<true> (r0, en, ed, *ai1, group2_com_atom);
switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
}
group2.set_weighted_gradient(group2_com_atom.grad);
@ -181,27 +181,15 @@ void colvar::coordnum::calc_gradients()
if (b_anisotropic) {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
switching_function<true> (r0_vec, en, ed, *ai1, *ai2);
switching_function<true>(r0_vec, en, ed, *ai1, *ai2);
}
} else {
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++)
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
switching_function<true> (r0, en, ed, *ai1, *ai2);
switching_function<true>(r0, en, ed, *ai1, *ai2);
}
}
}
// if (cvm::debug()) {
// for (size_t i = 0; i < group1.size(); i++) {
// cvm::log("atom["+cvm::to_str (group1[i].id+1)+"] gradient: "+
// cvm::to_str (group1[i].grad)+"\n");
// }
// for (size_t i = 0; i < group2.size(); i++) {
// cvm::log("atom["+cvm::to_str (group2[i].id+1)+"] gradient: "+
// cvm::to_str (group2[i].grad)+"\n");
// }
// }
}
void colvar::coordnum::apply_force(colvarvalue const &force)
@ -234,8 +222,8 @@ colvar::h_bond::h_bond(std::string const &conf)
cvm::fatal_error("Error: either acceptor or donor undefined.\n");
}
acceptor = cvm::atom(a_num);
donor = cvm::atom(d_num);
cvm::atom acceptor = cvm::atom(a_num);
cvm::atom donor = cvm::atom(d_num);
atom_groups.push_back(new cvm::atom_group);
atom_groups[0]->add_atom(acceptor);
atom_groups[0]->add_atom(donor);
@ -253,12 +241,10 @@ colvar::h_bond::h_bond(std::string const &conf)
}
colvar::h_bond::h_bond(cvm::atom const &acceptor_i,
cvm::atom const &donor_i,
cvm::real r0_i, int en_i, int ed_i)
: acceptor(acceptor_i),
donor(donor_i),
r0(r0_i), en(en_i), ed(ed_i)
colvar::h_bond::h_bond(cvm::atom const &acceptor,
cvm::atom const &donor,
cvm::real r0_i, int en_i, int ed_i)
: r0(r0_i), en(en_i), ed(ed_i)
{
function_type = "h_bond";
x.type(colvarvalue::type_scalar);
@ -278,30 +264,25 @@ colvar::h_bond::h_bond()
colvar::h_bond::~h_bond()
{
for (unsigned int i=0; i<atom_groups.size(); i++) {
delete atom_groups[i];
}
delete atom_groups[0];
}
void colvar::h_bond::calc_value()
{
x.real_value = colvar::coordnum::switching_function<false> (r0, en, ed, acceptor, donor);
x.real_value = colvar::coordnum::switching_function<false>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
}
void colvar::h_bond::calc_gradients()
{
colvar::coordnum::switching_function<true> (r0, en, ed, acceptor, donor);
(*atom_groups[0])[0].grad = acceptor.grad;
(*atom_groups[0])[1].grad = donor.grad;
colvar::coordnum::switching_function<true>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
}
void colvar::h_bond::apply_force(colvarvalue const &force)
{
acceptor.apply_force(force.real_value * acceptor.grad);
donor.apply_force (force.real_value * donor.grad);
(atom_groups[0])->apply_colvar_force(force);
}
@ -339,23 +320,27 @@ colvar::selfcoordnum::selfcoordnum()
void colvar::selfcoordnum::calc_value()
{
x.real_value = 0.0;
for (size_t i = 0; i < group1.size() - 1; i++)
for (size_t j = i + 1; j < group1.size(); j++)
x.real_value += colvar::coordnum::switching_function<false> (r0, en, ed, group1[i], group1[j]);
for (size_t i = 0; i < group1.size() - 1; i++) {
for (size_t j = i + 1; j < group1.size(); j++) {
x.real_value += colvar::coordnum::switching_function<false>(r0, en, ed, group1[i], group1[j]);
}
}
}
void colvar::selfcoordnum::calc_gradients()
{
for (size_t i = 0; i < group1.size() - 1; i++)
for (size_t j = i + 1; j < group1.size(); j++)
colvar::coordnum::switching_function<true> (r0, en, ed, group1[i], group1[j]);
for (size_t i = 0; i < group1.size() - 1; i++) {
for (size_t j = i + 1; j < group1.size(); j++) {
colvar::coordnum::switching_function<true>(r0, en, ed, group1[i], group1[j]);
}
}
}
void colvar::selfcoordnum::apply_force(colvarvalue const &force)
{
if (!group1.noforce)
if (!group1.noforce) {
group1.apply_colvar_force(force.real_value);
}
}

63
lib/colvars/colvarcomp_distances.cpp
View File

@ -447,9 +447,10 @@ colvar::distance_inv::distance_inv(std::string const &conf)
for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) {
for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) {
if (ai1->id == ai2->id)
cvm::error("Error: group1 and group1 have some atoms in common: this is not allowed for distanceInv.\n");
return;
if (ai1->id == ai2->id) {
cvm::error("Error: group1 and group2 have some atoms in common: this is not allowed for distanceInv.\n");
return;
}
}
}
@ -853,10 +854,11 @@ colvar::rmsd::rmsd(std::string const &conf)
if (get_keyval(conf, "refPositionsCol", ref_pos_col, std::string(""))) {
// if provided, use PDB column to select coordinates
bool found = get_keyval(conf, "refPositionsColValue", ref_pos_col_value, 0.0);
if (found && !ref_pos_col_value)
if (found && !ref_pos_col_value) {
cvm::error("Error: refPositionsColValue, "
"if provided, must be non-zero.\n");
"if provided, must be non-zero.\n");
return;
}
} else {
// if not, rely on existing atom indices for the group
atoms.create_sorted_ids();
@ -1294,17 +1296,40 @@ colvar::cartesian::cartesian(std::string const &conf)
parse_group(conf, "atoms", atoms);
atom_groups.push_back(&atoms);
bool use_x, use_y, use_z;
get_keyval(conf, "useX", use_x, true);
get_keyval(conf, "useY", use_y, true);
get_keyval(conf, "useZ", use_z, true);
axes.clear();
if (use_x) axes.push_back(0);
if (use_y) axes.push_back(1);
if (use_z) axes.push_back(2);
if (axes.size() == 0) {
cvm::error("Error: a \"cartesian\" component was defined with all three axes disabled.\n");
}
x.type(colvarvalue::type_vector);
x.vector1d_value.resize(atoms.size() * 3);
x.vector1d_value.resize(atoms.size() * axes.size());
}
void colvar::cartesian::calc_value()
{
int ia, j;
size_t const dim = axes.size();
size_t ia, j;
for (ia = 0; ia < atoms.size(); ia++) {
for (j = 0; j < 3; j++) {
x.vector1d_value[3*ia + j] = atoms[ia].pos[j];
for (j = 0; j < dim; j++) {
x.vector1d_value[dim*ia + j] = atoms[ia].pos[axes[j]];
}
}
if (atoms.weights.size()) {
for (ia = 0; ia < atoms.size(); ia++) {
for (j = 0; j < dim; j++) {
x.vector1d_value[dim*ia + j] *= atoms.weights[ia];
}
}
}
}
@ -1320,14 +1345,24 @@ void colvar::cartesian::calc_gradients()
void colvar::cartesian::apply_force(colvarvalue const &force)
{
int ia, j;
size_t const dim = axes.size();
size_t ia, j;
if (!atoms.noforce) {
cvm::rvector f;
for (ia = 0; ia < atoms.size(); ia++) {
for (j = 0; j < 3; j++) {
f[j] = force.vector1d_value[3*ia + j];
if (atoms.weights.size()) {
for (ia = 0; ia < atoms.size(); ia++) {
for (j = 0; j < dim; j++) {
f[axes[j]] = force.vector1d_value[dim*ia + j] / atoms.weights[ia];
}
atoms[ia].apply_force(f);
}
} else {
for (ia = 0; ia < atoms.size(); ia++) {
for (j = 0; j < dim; j++) {
f[axes[j]] = force.vector1d_value[dim*ia + j];
}
atoms[ia].apply_force(f);
}
atoms[ia].apply_force(f);
}
}
}

12
lib/colvars/colvarcomp_protein.cpp
View File

@ -66,8 +66,11 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
for (size_t i = 0; i < residues.size()-2; i++) {
theta.push_back(new colvar::angle(cvm::atom(r[i ], "CA", sid),
cvm::atom(r[i+1], "CA", sid),
cvm::atom(r[i+2], "CA", sid)));
cvm::atom(r[i+1], "CA", sid),
cvm::atom(r[i+2], "CA", sid)));
atom_groups.push_back(theta.back()->atom_groups[0]);
atom_groups.push_back(theta.back()->atom_groups[1]);
atom_groups.push_back(theta.back()->atom_groups[2]);
}
} else {
@ -85,8 +88,9 @@ colvar::alpha_angles::alpha_angles(std::string const &conf)
for (size_t i = 0; i < residues.size()-4; i++) {
hb.push_back(new colvar::h_bond(cvm::atom(r[i ], "O", sid),
cvm::atom(r[i+4], "N", sid),
r0, en, ed));
cvm::atom(r[i+4], "N", sid),
r0, en, ed));
atom_groups.push_back(hb.back()->atom_groups[0]);
}
} else {

8
lib/colvars/colvargrid.cpp
View File

@ -10,16 +10,18 @@
colvar_grid_count::colvar_grid_count()
: colvar_grid<size_t>()
{}
{
mult = 1;
}
colvar_grid_count::colvar_grid_count(std::vector<int> const &nx_i,
size_t const &def_count)
: colvar_grid<size_t>(nx_i, def_count)
: colvar_grid<size_t>(nx_i, def_count, 1)
{}
colvar_grid_count::colvar_grid_count(std::vector<colvar *> &colvars,
size_t const &def_count)
: colvar_grid<size_t>(colvars, def_count)
: colvar_grid<size_t>(colvars, def_count, 1)
{}
std::istream & colvar_grid_count::read_restart(std::istream &is)

261
lib/colvars/colvargrid.h
View File

@ -15,7 +15,7 @@
/// \brief Grid of values of a function of several collective
/// variables \param T The data type
///
/// Only scalar colvars supported so far
/// Only scalar colvars supported so far: vector colvars are treated as arrays
template <class T> class colvar_grid : public colvarparse {
protected:
@ -56,7 +56,7 @@ protected:
addr += ix[i]*nxc[i];
if (cvm::debug()) {
if (ix[i] >= nx[i]) {
cvm::error("Error: exceeding bounds in colvar_grid.\n");
cvm::error("Error: exceeding bounds in colvar_grid.\n", BUG_ERROR);
return 0;
}
}
@ -125,35 +125,51 @@ public:
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)
/// \brief Allocate data
int setup(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;
if (cvm::debug()) {
cvm::log("Allocating grid: multiplicity = "+cvm::to_str(mult_i)+
", dimensions = "+cvm::to_str(nx_i)+".\n");
}
mult = mult_i;
data.clear();
nx = nx_i;
nd = nx.size();
nxc.resize(nd);
// setup dimensions
nt = mult;
for (int i = nd-1; i >= 0; i--) {
if (nx_i[i] <= 0) {
cvm::error("Error: providing an invalid number of points, "+
cvm::to_str(nx_i[i])+".\n");
return;
if (nx[i] <= 0) {
cvm::error("Error: providing an invalid number of grid points, "+
cvm::to_str(nx[i])+".\n", BUG_ERROR);
return COLVARS_ERROR;
}
nxc[i] = nt;
nt *= nx[i];
}
if (cvm::debug()) {
cvm::log("Total number of grid elements = "+cvm::to_str(nt)+".\n");
}
data.reserve(nt);
data.assign(nt, t);
return COLVARS_OK;
}
/// \brief Allocate data (allow initialization also after construction)
void create()
int setup()
{
create(this->nx, T(), this->mult);
return setup(this->nx, T(), this->mult);
}
/// \brief Reset data (in case the grid is being reused)
@ -168,6 +184,8 @@ public:
{
save_delimiters = false;
nd = nt = 0;
mult = 1;
this->setup();
}
/// Destructor
@ -176,20 +194,20 @@ public:
/// \brief "Almost copy-constructor": only copies configuration
/// parameters from another grid, but doesn't reallocate stuff;
/// create() must be called after that;
/// setup() must be called after that;
colvar_grid(colvar_grid<T> const &g) : nd(g.nd),
nx(g.nx),
mult(g.mult),
data(),
cv(g.cv),
actual_value(g.actual_value),
lower_boundaries(g.lower_boundaries),
upper_boundaries(g.upper_boundaries),
periodic(g.periodic),
hard_lower_boundaries(g.hard_lower_boundaries),
hard_upper_boundaries(g.hard_upper_boundaries),
widths(g.widths),
has_data(false)
nx(g.nx),
mult(g.mult),
data(),
cv(g.cv),
actual_value(g.actual_value),
lower_boundaries(g.lower_boundaries),
upper_boundaries(g.upper_boundaries),
periodic(g.periodic),
hard_lower_boundaries(g.hard_lower_boundaries),
hard_upper_boundaries(g.hard_upper_boundaries),
widths(g.widths),
has_data(false)
{
save_delimiters = false;
}
@ -199,48 +217,68 @@ public:
/// 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)
T const &t = T(),
size_t mult_i = 1)
: has_data(false)
{
save_delimiters = false;
this->create(nx_i, t, mult_i);
this->setup(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)
T const &t = T(),
size_t mult_i = 1,
bool margin = false)
: has_data(false)
{
save_delimiters = false;
this->init_from_colvars(colvars, t, mult_i, margin);
}
std::vector<int> nx_i;
int init_from_colvars(std::vector<colvar *> const &colvars,
T const &t = T(),
size_t mult_i = 1,
bool margin = false)
{
if (cvm::debug()) {
cvm::log("Reading grid configuration from collective variables.\n");
}
if (cvm::debug())
cv = colvars;
nd = colvars.size();
mult = mult_i;
size_t i;
for (i = 0; i < cv.size(); i++) {
if (!cv[i]->tasks[colvar::task_lower_boundary] ||
!cv[i]->tasks[colvar::task_upper_boundary]) {
cvm::error("Tried to initialize a grid on a "
"variable with undefined boundaries.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
}
if (cvm::debug()) {
cvm::log("Allocating a grid for "+cvm::to_str(colvars.size())+
" collective variables.\n");
" collective variables, multiplicity = "+cvm::to_str(mult_i)+".\n");
}
for (size_t i = 0; i < cv.size(); i++) {
for (i = 0; i < cv.size(); i++) {
if (cv[i]->value().type() != colvarvalue::type_scalar) {
cvm::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");
return;
"constructed for scalar variables. "
"ABF and histogram can not be used; metadynamics "
"can be used with useGrids disabled.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
if (cv[i]->width <= 0.0) {
cvm::error("Tried to initialize a grid on a "
"variable with negative or zero width.\n");
return;
}
if (!cv[i]->tasks[colvar::task_lower_boundary] || !cv[i]->tasks[colvar::task_upper_boundary]) {
cvm::error("Tried to initialize a grid on a "
"variable with undefined boundaries.\n");
return;
"variable with negative or zero width.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
widths.push_back(cv[i]->width);
@ -271,35 +309,49 @@ public:
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);
this->init_from_boundaries();
return this->setup();
}
int init_from_boundaries(T const &t = T(),
size_t const &mult_i = 1)
{
if (cvm::debug()) {
cvm::log("Configuring grid dimensions from colvars boundaries.\n");
}
// these will have to be updated
nx.clear();
nxc.clear();
nt = 0;
for (size_t i = 0; i < lower_boundaries.size(); i++) {
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.push_back(nbins_round);
}
return COLVARS_OK;
}
/// Wrap an index vector around periodic boundary conditions
/// also checks validity of non-periodic indices
@ -309,14 +361,16 @@ public:
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::error("Trying to wrap illegal index vector(non-PBC): "
+ cvm::to_str(ix));
if (ix[i] < 0 || ix[i] >= nx[i]) {
cvm::error("Trying to wrap illegal index vector (non-PBC) for a grid point: "
+ cvm::to_str(ix), BUG_ERROR);
return;
}
}
}
}
/// \brief Report the bin corresponding to the current value of variable i
inline int current_bin_scalar(int const i) const
{
@ -373,7 +427,7 @@ public:
cvm::error("Error: trying to subtract two grids with "
"different multiplicity.\n");
return;
}
}
if (other_grid.data.size() != this->data.size()) {
cvm::error("Error: trying to subtract two grids with "
@ -395,8 +449,8 @@ public:
if (other_grid.multiplicity() != this->multiplicity()) {
cvm::error("Error: trying to copy two grids with "
"different multiplicity.\n");
return;
}
return;
}
if (other_grid.data.size() != this->data.size()) {
cvm::error("Error: trying to copy two grids with "
@ -673,9 +727,11 @@ public:
return os;
}
bool parse_params(std::string const &conf)
/// Read a grid definition from a config string
int parse_params(std::string const &conf)
{
if (cvm::debug()) cvm::log("Reading grid configuration from string.\n");
std::vector<int> old_nx = nx;
std::vector<colvarvalue> old_lb = lower_boundaries;
@ -688,36 +744,52 @@ public:
cvm::to_str(nd_in)+") than the grid defined "
"in the configuration file("+cvm::to_str(nd)+
").\n");
return false;
return COLVARS_ERROR;
}
}
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);
// support also camel case
colvarparse::get_keyval(conf, "lowerBoundaries",
lower_boundaries, lower_boundaries, colvarparse::parse_silent);
colvarparse::get_keyval(conf, "upperBoundaries",
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);
if (nd < lower_boundaries.size()) nd = lower_boundaries.size();
if (! actual_value.size()) actual_value.assign(nd, false);
if (! periodic.size()) periodic.assign(nd, false);
if (! widths.size()) widths.assign(nd, 1.0);
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],
if (old_nx.size()) {
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;
new_params = true;
}
}
} else {
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);
init_from_boundaries();
// data.resize(0); // no longer needed: setup calls clear()
return this->setup(nx, T(), mult);
}
return true;
return COLVARS_OK;
}
/// \brief Check that the grid information inside (boundaries,
@ -834,6 +906,7 @@ public:
<< periodic[i] << "\n";
}
for (std::vector<int> ix = new_index(); index_ok(ix); incr(ix) ) {
if (ix.back() == 0) {

152
lib/colvars/colvarmodule.cpp
View File

@ -24,12 +24,12 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
// TODO relax this error to handle multiple molecules in VMD
// once the module is not static anymore
cvm::error("Error: trying to allocate the collective "
"variable module twice.\n");
"variable module twice.\n");
return;
}
cvm::log(cvm::line_marker);
cvm::log("Initializing the collective variables module, version "+
cvm::to_str(COLVARS_VERSION)+".\n");
cvm::to_str(COLVARS_VERSION)+".\n");
// set initial default values
@ -56,14 +56,14 @@ int colvarmodule::config_file(char const *config_filename)
{
cvm::log(cvm::line_marker);
cvm::log("Reading new configuration from file \""+
std::string(config_filename)+"\":\n");
std::string(config_filename)+"\":\n");
// open the configfile
config_s.open(config_filename);
if (!config_s) {
if (!config_s.is_open()) {
cvm::error("Error: in opening configuration file \""+
std::string(config_filename)+"\".\n",
FILE_ERROR);
std::string(config_filename)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
}
@ -134,8 +134,10 @@ int colvarmodule::config(std::string &conf)
cvm::log("Collective variables module (re)initialized.\n");
cvm::log(cvm::line_marker);
// configuration might have changed, better redo the labels
write_traj_label(cv_traj_os);
if (cv_traj_os.is_open()) {
// configuration might have changed, better redo the labels
write_traj_label(cv_traj_os);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
@ -151,24 +153,24 @@ int colvarmodule::parse_global_params(std::string const &conf)
parse->get_keyval(conf, "analysis", b_analysis, b_analysis);
parse->get_keyval(conf, "debugGradientsStepSize", debug_gradients_step_size,
debug_gradients_step_size,
colvarparse::parse_silent);
debug_gradients_step_size,
colvarparse::parse_silent);
parse->get_keyval(conf, "eigenvalueCrossingThreshold",
colvarmodule::rotation::crossing_threshold,
colvarmodule::rotation::crossing_threshold,
colvarparse::parse_silent);
colvarmodule::rotation::crossing_threshold,
colvarmodule::rotation::crossing_threshold,
colvarparse::parse_silent);
parse->get_keyval(conf, "colvarsTrajFrequency", cv_traj_freq, cv_traj_freq);
parse->get_keyval(conf, "colvarsRestartFrequency",
restart_out_freq, restart_out_freq);
restart_out_freq, restart_out_freq);
// if this is true when initializing, it means
// we are continuing after a reset(): default to true
parse->get_keyval(conf, "colvarsTrajAppend", cv_traj_append, cv_traj_append);
parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false,
colvarparse::parse_silent);
colvarparse::parse_silent);
if (use_scripted_forces && !proxy->force_script_defined) {
cvm::fatal_error("User script for scripted colvar forces not found.");
@ -193,8 +195,8 @@ int colvarmodule::parse_colvars(std::string const &conf)
colvars.push_back(new colvar(colvar_conf));
if (cvm::get_error() ||
((colvars.back())->check_keywords(colvar_conf, "colvar") != COLVARS_OK)) {
cvm::log("Error while constructing colvar number " +
cvm::to_str(colvars.size()) + " : deleting.");
cvm::log("Error while constructing colvar number " +
cvm::to_str(colvars.size()) + " : deleting.");
delete colvars.back(); // the colvar destructor updates the colvars array
return COLVARS_ERROR;
}
@ -214,8 +216,8 @@ int colvarmodule::parse_colvars(std::string const &conf)
if (colvars.size())
cvm::log(cvm::line_marker);
cvm::log("Collective variables initialized, "+
cvm::to_str(colvars.size())+
" in total.\n");
cvm::to_str(colvars.size())+
" in total.\n");
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
@ -223,9 +225,9 @@ int colvarmodule::parse_colvars(std::string const &conf)
bool colvarmodule::check_new_bias(std::string &conf, char const *key)
{
if (cvm::get_error() ||
(biases.back()->check_keywords(conf, key) != COLVARS_OK)) {
(biases.back()->check_keywords(conf, key) != COLVARS_OK)) {
cvm::log("Error while constructing bias number " +
cvm::to_str(biases.size()) + " : deleting.\n");
cvm::to_str(biases.size()) + " : deleting.\n");
delete biases.back(); // the bias destructor updates the biases array
return true;
}
@ -330,7 +332,7 @@ colvar *colvarmodule::colvar_by_name(std::string const &name) {
int colvarmodule::change_configuration(std::string const &bias_name,
std::string const &conf)
std::string const &conf)
{
// This is deprecated; supported strategy is to delete the bias
// and parse the new config
@ -357,7 +359,7 @@ std::string colvarmodule::read_colvar(std::string const &name)
}
cvm::real colvarmodule::energy_difference(std::string const &bias_name,
std::string const &conf)
std::string const &conf)
{
cvm::increase_depth();
colvarbias *b;
@ -449,7 +451,7 @@ int colvarmodule::calc() {
if (cvm::debug()) {
cvm::log(cvm::line_marker);
cvm::log("Collective variables module, step no. "+
cvm::to_str(cvm::step_absolute())+"\n");
cvm::to_str(cvm::step_absolute())+"\n");
}
// calculate collective variables and their gradients
@ -533,7 +535,7 @@ int colvarmodule::calc() {
// equation of motion (extended system)
if (cvm::debug())
cvm::log("Updating the internal degrees of freedom "
"of colvars (if they have any).\n");
"of colvars (if they have any).\n");
cvm::increase_depth();
for (cvi = colvars.begin(); cvi != colvars.end(); cvi++) {
total_colvar_energy += (*cvi)->update();
@ -564,10 +566,10 @@ int colvarmodule::calc() {
if ( (cvm::step_relative() > 0) &&
((cvm::step_absolute() % restart_out_freq) == 0) ) {
cvm::log("Writing the state file \""+
restart_out_name+"\".\n");
restart_out_name+"\".\n");
proxy->backup_file(restart_out_name.c_str());
restart_out_os.open(restart_out_name.c_str());
if (!write_restart(restart_out_os))
if (!restart_out_os.is_open() || !write_restart(restart_out_os))
cvm::error("Error: in writing restart file.\n");
restart_out_os.close();
}
@ -594,7 +596,7 @@ int colvarmodule::calc() {
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_name+"\".\n");
cv_traj_os.flush();
}
}
@ -694,8 +696,8 @@ int colvarmodule::setup_input()
std::ifstream input_is(restart_in_name.c_str());
if (!input_is.good()) {
cvm::error("Error: in opening restart file \""+
std::string(restart_in_name)+"\".\n",
FILE_ERROR);
std::string(restart_in_name)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
} else {
cvm::log("Restarting from file \""+restart_in_name+"\".\n");
@ -725,9 +727,9 @@ int colvarmodule::setup_output()
output_prefix = proxy->output_prefix();
if (output_prefix.size()) {
cvm::log("The final output state file will be \""+
(output_prefix.size() ?
std::string(output_prefix+".colvars.state") :
std::string("colvars.state"))+"\".\n");
(output_prefix.size() ?
std::string(output_prefix+".colvars.state") :
std::string("colvars.state"))+"\".\n");
// cvm::log (cvm::line_marker);
}
@ -752,8 +754,8 @@ std::istream & colvarmodule::read_restart(std::istream &is)
if (is >> colvarparse::read_block("configuration", restart_conf)) {
if (it_restart_from_state_file) {
parse->get_keyval(restart_conf, "step",
it_restart, (size_t) 0,
colvarparse::parse_silent);
it_restart, (size_t) 0,
colvarparse::parse_silent);
it = it_restart;
}
}
@ -767,8 +769,8 @@ std::istream & colvarmodule::read_restart(std::istream &is)
cvi++) {
if ( !((*cvi)->read_restart(is)) ) {
cvm::error("Error: in reading restart configuration for collective variable \""+
(*cvi)->name+"\".\n",
INPUT_ERROR);
(*cvi)->name+"\".\n",
INPUT_ERROR);
}
}
@ -776,10 +778,11 @@ std::istream & colvarmodule::read_restart(std::istream &is)
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if (!((*bi)->read_restart(is)))
if (!((*bi)->read_restart(is))) {
cvm::error("Error: in reading restart configuration for bias \""+
(*bi)->name+"\".\n",
INPUT_ERROR);
(*bi)->name+"\".\n",
INPUT_ERROR);
}
}
cvm::decrease_depth();
@ -802,11 +805,11 @@ int colvarmodule::write_output_files()
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();
std::ostream * os = proxy->output_stream(out_name);
os->setf(std::ios::scientific, std::ios::floatfield);
this->write_restart(*os);
proxy->close_output_stream(out_name);
cvm::increase_depth();
for (std::vector<colvar *>::iterator cvi = colvars.begin();
@ -824,11 +827,11 @@ int colvarmodule::write_output_files()
int colvarmodule::read_traj(char const *traj_filename,
size_t traj_read_begin,
size_t traj_read_end)
size_t traj_read_begin,
size_t traj_read_end)
{
cvm::log("Opening trajectory file \""+
std::string(traj_filename)+"\".\n");
std::string(traj_filename)+"\".\n");
std::ifstream traj_is(traj_filename);
while (true) {
@ -839,7 +842,7 @@ int colvarmodule::read_traj(char const *traj_filename,
do {
if (!colvarparse::getline_nocomments(traj_is, line)) {
cvm::log("End of file \""+std::string(traj_filename)+
"\" reached, or corrupted file.\n");
"\" reached, or corrupted file.\n");
traj_is.close();
return false;
}
@ -867,8 +870,8 @@ int colvarmodule::read_traj(char const *traj_filename,
(it > traj_read_end) ) {
std::cerr << "\n";
cvm::error("Reached the end of the trajectory, "
"read_end = "+cvm::to_str(traj_read_end)+"\n",
FILE_ERROR);
"read_end = "+cvm::to_str(traj_read_end)+"\n",
FILE_ERROR);
return COLVARS_ERROR;
}
@ -877,9 +880,9 @@ int colvarmodule::read_traj(char const *traj_filename,
cvi++) {
if (!(*cvi)->read_traj(is)) {
cvm::error("Error: in reading colvar \""+(*cvi)->name+
"\" from trajectory file \""+
std::string(traj_filename)+"\".\n",
FILE_ERROR);
"\" from trajectory file \""+
std::string(traj_filename)+"\".\n",
FILE_ERROR);
return COLVARS_ERROR;
}
}
@ -927,18 +930,18 @@ int colvarmodule::open_traj_file(std::string const &file_name)
// (re)open trajectory file
if (cv_traj_append) {
cvm::log("Appending to colvar trajectory file \""+file_name+
"\".\n");
"\".\n");
cv_traj_os.open(file_name.c_str(), std::ios::app);
} else {
cvm::log("Writing to colvar trajectory file \""+file_name+
"\".\n");
"\".\n");
proxy->backup_file(file_name.c_str());
cv_traj_os.open(file_name.c_str(), std::ios::out);
cv_traj_os.open(file_name.c_str());
}
if (!cv_traj_os.is_open()) {
cvm::error("Error: cannot write to file \""+file_name+"\".\n",
FILE_ERROR);
FILE_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
@ -1046,10 +1049,11 @@ void cvm::exit(std::string const &message)
int cvm::read_index_file(char const *filename)
{
std::ifstream is(filename, std::ios::binary);
if (!is.good())
if (!is.good()) {
cvm::error("Error: in opening index file \""+
std::string(filename)+"\".\n",
FILE_ERROR);
std::string(filename)+"\".\n",
FILE_ERROR);
}
while (is.good()) {
char open, close;
@ -1062,8 +1066,8 @@ int cvm::read_index_file(char const *filename)
names_i++) {
if (*names_i == group_name) {
cvm::error("Error: the group name \""+group_name+
"\" appears in multiple index files.\n",
FILE_ERROR);
"\" appears in multiple index files.\n",
FILE_ERROR);
}
}
cvm::index_group_names.push_back(group_name);
@ -1093,7 +1097,7 @@ int cvm::read_index_file(char const *filename)
}
cvm::log("The following index groups were read from the index file \""+
std::string(filename)+"\":\n");
std::string(filename)+"\":\n");
std::list<std::string>::iterator names_i = index_group_names.begin();
std::list<std::vector<int> >::iterator lists_i = index_groups.begin();
for ( ; names_i != index_group_names.end() ; names_i++, lists_i++) {
@ -1103,18 +1107,18 @@ int cvm::read_index_file(char const *filename)
}
int cvm::load_atoms(char const *file_name,
std::vector<cvm::atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value)
std::vector<cvm::atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value)
{
return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value);
}
int 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)
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value)
{
// Differentiate between PDB and XYZ files
// for XYZ files, use CVM internal parser
@ -1134,8 +1138,8 @@ int cvm::load_coords(char const *file_name,
int cvm::load_coords_xyz(char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices)
std::vector<atom_pos> &pos,
const std::vector<int> &indices)
{
std::ifstream xyz_is(filename);
unsigned int natoms;
@ -1144,7 +1148,7 @@ int cvm::load_coords_xyz(char const *filename,
if ( ! (xyz_is >> natoms) ) {
cvm::error("Error: cannot parse XYZ file "
+ std::string(filename) + ".\n", INPUT_ERROR);
+ std::string(filename) + ".\n", INPUT_ERROR);
}
// skip comment line
std::getline(xyz_is, line);

23
lib/colvars/colvarmodule.h
View File

@ -4,7 +4,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2014-11-21"
#define COLVARS_VERSION "2015-02-04"
#endif
#ifndef COLVARS_DEBUG
@ -45,6 +45,10 @@
#include <vector>
#include <list>
#ifdef NAMD_VERSION
// use Lustre-friendly wrapper to POSIX write()
#include "fstream_namd.h"
#endif
class colvarparse;
class colvar;
@ -239,6 +243,12 @@ public:
/// (Re)initialize the output trajectory and state file (does not write it yet)
int setup_output();
#ifdef NAMD_VERSION
typedef ofstream_namd ofstream;
#else
typedef std::ofstream ofstream;
#endif
/// Read the input restart file
std::istream & read_restart(std::istream &is);
/// Write the output restart file
@ -443,6 +453,7 @@ public:
/// Pseudo-random number with Gaussian distribution
static real rand_gaussian(void);
protected:
/// Configuration file
@ -452,16 +463,16 @@ protected:
colvarparse *parse;
/// Name of the trajectory file
std::string cv_traj_name;
std::string cv_traj_name;
/// Collective variables output trajectory file
std::ofstream cv_traj_os;
colvarmodule::ofstream cv_traj_os;
/// Appending to the existing trajectory file?
bool cv_traj_append;
bool cv_traj_append;
/// Output restart file
std::ofstream restart_out_os;
colvarmodule::ofstream restart_out_os;
/// \brief Counter for the current depth in the object hierarchy (useg e.g. in output
static size_t depth;
@ -470,12 +481,12 @@ protected:
static bool use_scripted_forces;
public:
/// \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;
/// Increase the depth (number of indentations in the output)
static void increase_depth();

36
lib/colvars/colvarparse.cpp
View File

@ -47,20 +47,22 @@ size_t colvarparse::dummy_pos = 0;
if (data.size()) { \
std::istringstream is(data); \
TYPE x(def_value); \
if (is >> x) \
if (is >> x) { \
value = x; \
else \
} else { \
cvm::error("Error: in parsing \""+ \
std::string(key)+"\".\n", INPUT_ERROR); \
} \
if (parse_mode != parse_silent) { \
cvm::log("# "+std::string(key)+" = "+ \
cvm::to_str(value)+"\n"); \
} \
} else { \
\
if (b_found_any) \
if (b_found_any) { \
cvm::error("Error: improper or missing value " \
"for \""+std::string(key)+"\".\n", INPUT_ERROR); \
} \
value = def_value; \
if (parse_mode != parse_silent) { \
cvm::log("# "+std::string(key)+" = \""+ \
@ -110,19 +112,21 @@ size_t colvarparse::dummy_pos = 0;
if (data_count == 0) \
cvm::fatal_error("Error: in parsing \""+ \
std::string(key)+"\".\n"); \
if (data_count > 1) \
if (data_count > 1) { \
cvm::error("Error: multiple values " \
"are not allowed for keyword \""+ \
std::string(key)+"\".\n", INPUT_ERROR); \
} \
if (parse_mode != parse_silent) { \
cvm::log("# "+std::string(key)+" = "+ \
cvm::to_str(value)+"\n"); \
} \
} else { \
\
if (b_found_any) \
if (b_found_any) { \
cvm::error("Error: improper or missing value " \
"for \""+std::string(key)+"\".\n", INPUT_ERROR); \
} \
value = def_value; \
if (parse_mode != parse_silent) { \
cvm::log("# "+std::string(key)+" = "+ \
@ -185,11 +189,12 @@ size_t colvarparse::dummy_pos = 0;
size_t i = 0; \
for ( ; i < values.size(); i++) { \
TYPE x(values[i]); \
if (is >> x) \
if (is >> x) { \
values[i] = x; \
else \
} else { \
cvm::error("Error: in parsing \""+ \
std::string(key)+"\".\n", INPUT_ERROR); \
} \
} \
} \
\
@ -200,9 +205,10 @@ size_t colvarparse::dummy_pos = 0;
\
} else { \
\
if (b_found_any) \
if (b_found_any) { \
cvm::error("Error: improper or missing values for \""+ \
std::string(key)+"\".\n", INPUT_ERROR); \
} \
\
for (size_t i = 0; i < values.size(); i++) \
values[i] = def_values[ (i > def_values.size()) ? 0 : i ]; \
@ -476,16 +482,6 @@ bool colvarparse::key_lookup(std::string const &conf,
}
}
// check it is not between quotes
// if ( (conf.find_last_of ("\"",
// conf.find_last_of (white_space, pos)) !=
// std::string::npos) &&
// (conf.find_first_of ("\"",
// conf.find_first_of (white_space, pos)) !=
// std::string::npos) )
// return false;
// save the pointer for a future call (when iterating over multiple
// valid instances of the same keyword)
save_pos = pos + key.size();
@ -493,7 +489,7 @@ bool colvarparse::key_lookup(std::string const &conf,
// get the remainder of the line
size_t pl = conf.rfind("\n", pos);
size_t line_begin = (pl == std::string::npos) ? 0 : pos;
size_t nl = conf.find ("\n", pos);
size_t nl = conf.find("\n", pos);
size_t line_end = (nl == std::string::npos) ? conf.size() : nl;
std::string line(conf, line_begin, (line_end-line_begin));
@ -563,7 +559,7 @@ bool colvarparse::key_lookup(std::string const &conf,
if (data.size() && save_delimiters) {
data_begin_pos.push_back(conf.find(data, pos+key.size()));
data_end_pos.push_back (data_begin_pos.back()+data.size());
data_end_pos.push_back(data_begin_pos.back()+data.size());
// std::cerr << "key = " << key << ", data = \""
// << data << "\", data_begin, data_end = "
// << data_begin_pos.back() << ", " << data_end_pos.back()

52
lib/colvars/colvarparse.h
View File

@ -92,36 +92,36 @@ public:
/// wrapper class (colvarvalue.h).
#define _get_keyval_scalar_proto_(_type_,_def_value_) \
bool get_keyval(std::string const &conf, \
char const *key, \
_type_ &value, \
_type_ const &def_value = _def_value_, \
Parse_Mode const parse_mode = parse_normal)
bool get_keyval(std::string const &conf, \
char const *key, \
_type_ &value, \
_type_ const &def_value = _def_value_, \
Parse_Mode const parse_mode = parse_normal)
_get_keyval_scalar_proto_(int, (int)0);
_get_keyval_scalar_proto_(size_t, (size_t)0);
_get_keyval_scalar_proto_(std::string, std::string(""));
_get_keyval_scalar_proto_(cvm::real, (cvm::real)0.0);
_get_keyval_scalar_proto_(cvm::rvector, cvm::rvector());
_get_keyval_scalar_proto_(cvm::quaternion, cvm::quaternion());
_get_keyval_scalar_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset));
_get_keyval_scalar_proto_(bool, false);
_get_keyval_scalar_proto_(int, (int)0);
_get_keyval_scalar_proto_(size_t, (size_t)0);
_get_keyval_scalar_proto_(std::string, std::string(""));
_get_keyval_scalar_proto_(cvm::real, (cvm::real)0.0);
_get_keyval_scalar_proto_(cvm::rvector, cvm::rvector());
_get_keyval_scalar_proto_(cvm::quaternion, cvm::quaternion());
_get_keyval_scalar_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset));
_get_keyval_scalar_proto_(bool, false);
#define _get_keyval_vector_proto_(_type_,_def_value_) \
bool get_keyval(std::string const &conf, \
char const *key, \
std::vector<_type_> &values, \
std::vector<_type_> const &def_values = \
std::vector<_type_> (0, static_cast<_type_>(_def_value_)), \
Parse_Mode const parse_mode = parse_normal)
bool get_keyval(std::string const &conf, \
char const *key, \
std::vector<_type_> &values, \
std::vector<_type_> const &def_values = \
std::vector<_type_> (0, static_cast<_type_>(_def_value_)), \
Parse_Mode const parse_mode = parse_normal)
_get_keyval_vector_proto_(int, 0);
_get_keyval_vector_proto_(size_t, 0);
_get_keyval_vector_proto_(std::string, std::string(""));
_get_keyval_vector_proto_(cvm::real, 0.0);
_get_keyval_vector_proto_(cvm::rvector, cvm::rvector());
_get_keyval_vector_proto_(cvm::quaternion, cvm::quaternion());
_get_keyval_vector_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset));
_get_keyval_vector_proto_(int, 0);
_get_keyval_vector_proto_(size_t, 0);
_get_keyval_vector_proto_(std::string, std::string(""));
_get_keyval_vector_proto_(cvm::real, 0.0);
_get_keyval_vector_proto_(cvm::rvector, cvm::rvector());
_get_keyval_vector_proto_(cvm::quaternion, cvm::quaternion());
_get_keyval_vector_proto_(colvarvalue, colvarvalue(colvarvalue::type_notset));
/// \brief Check that all the keywords within "conf" are in the list

89
lib/colvars/colvarproxy.h
View File

@ -3,6 +3,8 @@
#ifndef COLVARPROXY_H
#define COLVARPROXY_H
#include <fstream>
#include <list>
#include "colvarmodule.h"
#include "colvarvalue.h"
@ -128,25 +130,25 @@ public:
/// \brief Get the PBC-aware distance vector between two positions
virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2) = 0;
cvm::atom_pos const &pos2) = 0;
/// \brief Get the PBC-aware square distance between two positions;
/// may be implemented independently from position_distance() for optimization purposes
virtual cvm::real position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2);
cvm::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 The reference position
virtual void select_closest_image(cvm::atom_pos &pos,
cvm::atom_pos const &ref_pos) = 0;
cvm::atom_pos const &ref_pos) = 0;
/// \brief Perform select_closest_image() on a set of atomic positions
///
/// After that, distance vectors can then be calculated directly,
/// without using position_distance()
void select_closest_images(std::vector<cvm::atom_pos> &pos,
cvm::atom_pos const &ref_pos);
cvm::atom_pos const &ref_pos);
// **************** SCRIPTING INTERFACE ****************
@ -164,13 +166,13 @@ public:
virtual int run_force_callback() { return COLVARS_NOT_IMPLEMENTED; }
virtual int run_colvar_callback(std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
colvarvalue &value)
std::vector<const colvarvalue *> const &cvcs,
colvarvalue &value)
{ return COLVARS_NOT_IMPLEMENTED; }
virtual int run_colvar_gradient_callback(std::string const &name,
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient)
std::vector<const colvarvalue *> const &cvcs,
std::vector<cvm::matrix2d<cvm::real> > &gradient)
{ return COLVARS_NOT_IMPLEMENTED; }
// **************** INPUT/OUTPUT ****************
@ -193,27 +195,78 @@ public:
/// "filename" is a PDB file, use this field to determine which are
/// the atoms to be set
virtual int load_atoms(char const *filename,
std::vector<cvm::atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0) = 0;
std::vector<cvm::atom> &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0) = 0;
/// \brief Load the coordinates for a group of atoms from a file
/// (usually a PDB); if "pos" is already allocated, the number of its
/// elements must match the number of atoms in "filename"
virtual int load_coords(char const *filename,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value = 0.0) = 0;
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value = 0.0) = 0;
protected:
/// \brief Open output files: by default, these are ofstream objects.
/// Allows redefinition to implement different output mechanisms
std::list<std::ostream *> output_files;
/// \brief Identifiers for output_stream objects: by default, these are the names of the files
std::list<std::string> output_stream_names;
public:
// TODO the following definitions may be moved to a .cpp file
/// \brief Returns a reference to the given output channel;
/// if this is not open already, then open it
virtual std::ostream * output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
return *osi;
}
}
output_stream_names.push_back(output_name);
std::ofstream * os = new std::ofstream(output_name.c_str());
if (!os->is_open()) {
cvm::error("Error: cannot write to file \""+output_name+"\".\n",
FILE_ERROR);
}
output_files.push_back(os);
return os;
}
/// \brief Closes the given output channel
virtual int close_output_stream(std::string const &output_name)
{
std::list<std::ostream *>::iterator osi = output_files.begin();
std::list<std::string>::iterator osni = output_stream_names.begin();
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
((std::ofstream *) (*osi))->close();
output_files.erase(osi);
output_stream_names.erase(osni);
return COLVARS_OK;
}
}
return COLVARS_ERROR;
}
/// \brief Rename the given file, before overwriting it
virtual int backup_file(char const *filename)
{ return COLVARS_NOT_IMPLEMENTED; }
{
return COLVARS_NOT_IMPLEMENTED;
}
};
inline void colvarproxy::select_closest_images(std::vector<cvm::atom_pos> &pos,
cvm::atom_pos const &ref_pos)
cvm::atom_pos const &ref_pos)
{
for (std::vector<cvm::atom_pos>::iterator pi = pos.begin();
pi != pos.end(); ++pi) {
@ -222,7 +275,7 @@ inline void colvarproxy::select_closest_images(std::vector<cvm::atom_pos> &pos,
}
inline cvm::real colvarproxy::position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
cvm::atom_pos const &pos2)
{
return (position_distance(pos1, pos2)).norm2();
}

1
lib/colvars/colvartypes.h
View File

@ -324,7 +324,6 @@ protected:
public:
T * data;
size_t length;
friend class matrix2d;
inline row(T * const row_data, size_t const inner_length)
: data(row_data), length(inner_length)
{}