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

update colvars to state of 2014-10-21 with some additional refactoring

Browse Source
This commit is contained in:
Axel Kohlmeyer
2014-10-21 14:00:52 -04:00
parent 35d8036e32
commit 2267cc7839
13 changed files with 504 additions and 179 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/PDF/colvars-refman-lammps.pdf
View File

Binary file not shown.

105
lib/colvars/colvar.cpp
View File

@ -159,39 +159,34 @@ colvar::colvar (std::string const &conf)
enable(task_scripted);
cvm::log("This colvar uses scripted function \"" + scripted_function + "\".");
std::string type_str;
get_keyval (conf, "scriptedFunctionType", type_str, "scalar");
// Only accept scalar scripted colvars
// might accept other types when the infrastructure is in place
// for derivatives of vectors wrt vectors
x.type(colvarvalue::type_scalar);
x_reported.type(x.type());
x.type(colvarvalue::type_notset);
for (i = 0; i < colvarvalue::type_all; i++) {
if (type_str == colvarvalue::type_keyword[i]) {
x.type(colvarvalue::Type(i));
break;
}
}
if (x.type() == colvarvalue::type_notset) {
cvm::error("Could not parse scripted colvar type.");
return;
}
x_reported.type (x.type());
cvm::log(std::string("Expecting colvar value of type ")
+ colvarvalue::type_desc[x.type()]);
// Build ordered list of component values that will be
// passed to the script
// Sort array of cvcs based on values of componentExp
std::vector<cvc *> temp_vec;
for (i = 1; i <= cvcs.size(); i++) {
for (j = 0; j < cvcs.size(); j++) {
if (cvcs[j]->sup_np == int(i)) {
sorted_cvc_values.push_back(cvcs[j]->p_value());
temp_vec.push_back(cvcs[j]);
break;
}
}
}
if (sorted_cvc_values.size() != cvcs.size()) {
cvm::error("Could not find order numbers for all components"
if (temp_vec.size() != cvcs.size()) {
cvm::error("Could not find order numbers for all components "
"in componentExp values.");
return;
}
cvcs = temp_vec;
// Build ordered list of component values that will be
// passed to the script
for (j = 0; j < cvcs.size(); j++) {
sorted_cvc_values.push_back(cvcs[j]->p_value());
}
}
if (!tasks[task_scripted]) {
@ -203,19 +198,23 @@ colvar::colvar (std::string const &conf)
b_Jacobian_force = true;
}
// Test whether this is a single-component variable
// Decide whether the colvar is periodic
// Used to wrap extended DOF if extendedLagrangian is on
if (cvcs.size() == 1 && (cvcs[0])->b_periodic && (cvcs[0])->sup_np == 1
&& (cvcs[0])->sup_coeff == 1.0
if (cvcs.size() == 1 && (cvcs[0])->sup_np == 1
&& (cvcs[0])->sup_coeff == 1.0
&& !tasks[task_scripted]) {
this->b_periodic = true;
this->period = (cvcs[0])->period;
b_single_cvc = true;
b_periodic = (cvcs[0])->b_periodic;
period = (cvcs[0])->period;
// TODO write explicit wrap() function for colvars to allow for
// sup_coeff different from 1
// this->period = (cvcs[0])->period * (cvcs[0])->sup_coeff;
} else {
this->b_periodic = false;
this->period = 0.0;
b_single_cvc = false;
b_periodic = false;
period = 0.0;
}
// check the available features of each cvc
@ -817,13 +816,14 @@ void colvar::calc()
// each atom group will take care of its own ref_pos_group, if defined
}
}
if (tasks[task_output_velocity]) {
for (i = 0; i < cvcs.size(); i++) {
for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
cvcs[i]->atom_groups[ig]->read_velocities();
}
}
}
//// Don't try to get atom velocities, as no back-end currently implements it
// if (tasks[task_output_velocity] && !tasks[task_fdiff_velocity]) {
// for (i = 0; i < cvcs.size(); i++) {
// for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
// cvcs[i]->atom_groups[ig]->read_velocities();
// }
// }
// }
if (tasks[task_system_force]) {
for (i = 0; i < cvcs.size(); i++) {
for (ig = 0; ig < cvcs[i]->atom_groups.size(); ig++) {
@ -1162,8 +1162,11 @@ void colvar::communicate_forces()
for (i = 0; i < cvcs.size(); i++) {
cvm::increase_depth();
// Note: we need a dot product here
(cvcs[i])->apply_force (f * func_grads[i]);
// Force is scalar times colvarvalue (scalar or vector)
// Note: this can only handle scalar colvars (scalar values of f)
// A non-scalar colvar would need the gradient to be expressed
// as an order-2 tensor
(cvcs[i])->apply_force (f.real_value * func_grads[i]);
cvm::decrease_depth();
}
} else if (x.type() == colvarvalue::type_scalar) {
@ -1232,30 +1235,38 @@ bool colvar::periodic_boundaries() const
cvm::real colvar::dist2 (colvarvalue const &x1,
colvarvalue const &x2) const
{
return (cvcs[0])->dist2 (x1, x2);
if (b_single_cvc) {
return (cvcs[0])->dist2(x1, x2);
} else {
return x1.dist2(x2);
}
}
colvarvalue colvar::dist2_lgrad (colvarvalue const &x1,
colvarvalue const &x2) const
{
return (cvcs[0])->dist2_lgrad (x1, x2);
if (b_single_cvc) {
return (cvcs[0])->dist2_lgrad (x1, x2);
} else {
return x1.dist2_grad(x2);
}
}
colvarvalue colvar::dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const
{
return (cvcs[0])->dist2_rgrad (x1, x2);
}
cvm::real colvar::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
return (cvcs[0])->compare (x1, x2);
if (b_single_cvc) {
return (cvcs[0])->dist2_rgrad (x1, x2);
} else {
return x2.dist2_grad(x1);
}
}
void colvar::wrap (colvarvalue &x) const
{
(cvcs[0])->wrap (x);
if (b_single_cvc) {
(cvcs[0])->wrap (x);
}
return;
}

11
lib/colvars/colvar.h
View File

@ -83,6 +83,10 @@ public:
/// combination of \link cvc \endlink elements
bool b_linear;
/// \brief True if this \link colvar \endlink is equal to
/// its only constituent cvc
bool b_single_cvc;
/// \brief True if all \link cvc \endlink objects are capable
/// of calculating inverse gradients
bool b_inverse_gradients;
@ -344,13 +348,6 @@ public:
colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to compare colvar values
///
/// Handles correctly symmetries and periodic boundary conditions
cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to wrap a value into a standard interval
///

18
lib/colvars/colvarbias.cpp
View File

@ -124,6 +124,24 @@ cvm::real colvarbias::energy_difference(std::string const &conf)
}
// So far, these are only implemented in colvarsbias_abf
int colvarbias::bin_num()
{
cvm::error ("Error: bin_num() not implemented.\n");
return -1;
}
int colvarbias::current_bin()
{
cvm::error ("Error: current_bin() not implemented.\n");
return -1;
}
int colvarbias::bin_count(int bin_index)
{
cvm::error ("Error: bin_count() not implemented.\n");
return -1;
}
std::ostream & colvarbias::write_traj_label (std::ostream &os)
{
os << " ";

9
lib/colvars/colvarbias.h
View File

@ -29,6 +29,15 @@ public:
/// Calculate change in energy from using alternate configuration
virtual cvm::real energy_difference(std::string const &conf);
/// Give the total number of bins for a given bias.
virtual int bin_num();
/// Calculate the bin index for a given bias.
virtual int current_bin();
//// Give the count at a given bin index.
virtual int bin_count(int bin_index);
//// Share information between replicas, whatever it may be.
virtual void replica_share() {};
/// Perform analysis tasks
virtual inline void analyse() {}

145
lib/colvars/colvarbias_abf.cpp
View File

@ -44,6 +44,18 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
get_keyval (conf, "historyFreq", history_freq, 0);
b_history_files = (history_freq > 0);
// shared ABF
get_keyval (conf, "shared", shared_on, false);
if (shared_on) {
if (!cvm::replica_enabled || cvm::replica_num() <= 1)
cvm::error ("Error: shared ABF requires more than one replica.");
else
cvm::log ("shared ABF will be applied among "+ cvm::to_str(cvm::replica_num()) + " replicas.\n");
// If shared_freq is not set, we default to output_freq
get_keyval (conf, "sharedFreq", shared_freq, output_freq);
}
// ************* checking the associated colvars *******************
if (colvars.size() == 0) {
@ -109,6 +121,15 @@ colvarbias_abf::colvarbias_abf (std::string const &conf, char const *key)
gradients->samples = samples;
samples->has_parent_data = true;
// 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_gradients = new colvar_grid_gradient (colvars);
last_gradients->samples = last_samples;
last_samples->has_parent_data = true;
shared_last_step = -1;
// If custom grids are provided, read them
if ( input_prefix.size() > 0 ) {
read_gradients_samples ();
@ -130,6 +151,19 @@ colvarbias_abf::~colvarbias_abf()
gradients = NULL;
}
// shared ABF
// We used to only do this if "shared" was defined,
// but now we can call shared externally
if (last_samples) {
delete last_samples;
last_samples = NULL;
}
if (last_gradients) {
delete last_gradients;
last_gradients = NULL;
}
delete [] force;
if (cvm::n_abf_biases > 0)
@ -225,13 +259,103 @@ cvm::real colvarbias_abf::update()
write_gradients_samples (output_prefix);
}
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
// otherwise, backup and replace
write_gradients_samples (output_prefix + ".hist", (cvm::step_absolute() > 0));
}
if (shared_on && shared_last_step >= 0 && cvm::step_absolute() % shared_freq == 0) {
// Share gradients and samples for shared ABF.
replica_share();
}
// Prepare for the first sharing.
if (shared_last_step < 0) {
// Copy the current gradient and count values into last.
last_gradients->copy_grid(*gradients);
last_samples->copy_grid(*samples);
shared_last_step = cvm::step_absolute();
cvm::log ("Prepared sample and gradient buffers at step "+cvm::to_str(cvm::step_absolute())+".");
}
return 0.0;
}
void colvarbias_abf::replica_share () {
int p;
if( !cvm::replica_enabled() ) {
cvm::error ("Error: shared ABF: No replicas.\n");
return;
}
// We must have stored the last_gradients and last_samples.
if (shared_last_step < 0 ) {
cvm::error ("Error: shared ABF: Tried to apply shared ABF before any sampling had occurred.\n");
return;
}
// Share gradients for shared ABF.
cvm::log ("shared ABF: Sharing gradient and samples among replicas at step "+cvm::to_str(cvm::step_absolute()) );
// Count of data items.
size_t data_n = gradients->raw_data_num();
size_t samp_start = data_n*sizeof(cvm::real);
size_t msg_total = data_n*sizeof(size_t) + samp_start;
char* msg_data = new char[msg_total];
if (cvm::replica_index() == 0) {
// Replica 0 collects the delta gradient and count from the others.
for (p = 1; p < cvm::replica_num(); p++) {
// Receive the deltas.
cvm::replica_comm_recv(msg_data, msg_total, p);
// Map the deltas from the others into the grids.
last_gradients->raw_data_in((cvm::real*)(&msg_data[0]));
last_samples->raw_data_in((size_t*)(&msg_data[samp_start]));
// Combine the delta gradient and count of the other replicas
// with Replica 0's current state (including its delta).
gradients->add_grid( *last_gradients );
samples->add_grid( *last_samples );
}
// Now we must send the combined gradient to the other replicas.
gradients->raw_data_out((cvm::real*)(&msg_data[0]));
samples->raw_data_out((size_t*)(&msg_data[samp_start]));
for (p = 1; p < cvm::replica_num(); p++) {
cvm::replica_comm_send(msg_data, msg_total, p);
}
} else {
// All other replicas send their delta gradient and count.
// Calculate the delta gradient and count.
last_gradients->delta_grid (*gradients);
last_samples->delta_grid (*samples);
// Cast the raw char data to the gradient and samples.
last_gradients->raw_data_out((cvm::real*)(&msg_data[0]));
last_samples->raw_data_out((size_t*)(&msg_data[samp_start]));
cvm::replica_comm_send(msg_data, msg_total, 0);
// We now receive the combined gradient from Replica 0.
cvm::replica_comm_recv(msg_data, msg_total, 0);
// We sync to the combined gradient computed by Replica 0.
gradients->raw_data_in((cvm::real*)(&msg_data[0]));
samples->raw_data_in((size_t*)(&msg_data[samp_start]));
}
// Without a barrier it's possible that one replica starts
// share 2 when other replicas haven't finished share 1.
cvm::replica_comm_barrier();
// Done syncing the replicas.
delete[] msg_data;
// Copy the current gradient and count values into last.
last_gradients->copy_grid(*gradients);
last_samples->copy_grid(*samples);
shared_last_step = cvm::step_absolute();
}
void colvarbias_abf::write_gradients_samples (const std::string &prefix, bool append)
{
@ -268,6 +392,27 @@ void colvarbias_abf::write_gradients_samples (const std::string &prefix, bool ap
return;
}
// For Tcl implementation of selection rules.
/// Give the total number of bins for a given bias.
int colvarbias_abf::bin_num() {
return samples->number_of_points(0);
}
/// Calculate the bin index for a given bias.
int colvarbias_abf::current_bin() {
return samples->current_bin_scalar(0);
}
/// Give the count at a given bin index.
int colvarbias_abf::bin_count(int bin_index) {
if (bin_index < 0 || bin_index >= bin_num()) {
cvm::error ("Error: Tried to get bin count from invalid bin index "+cvm::to_str(bin_index));
return -1;
}
std::vector<int> ix(1,(int)bin_index);
return samples->value(ix);
}
void colvarbias_abf::read_gradients_samples ()
{
std::string samples_in_name, gradients_in_name;

23
lib/colvars/colvarbias_abf.h
View File

@ -1,5 +1,3 @@
/// -*- c++ -*-
/************************************************************************
* Headers for the ABF and histogram biases *
************************************************************************/
@ -11,7 +9,6 @@
#include <list>
#include <sstream>
#include <iomanip>
//#include <cmath>
#include "colvarbias.h"
#include "colvargrid.h"
@ -62,8 +59,28 @@ private:
/// n-dim grid of number of samples
colvar_grid_count *samples;
// shared ABF
bool shared_on;
size_t shared_freq;
int shared_last_step;
// Share between replicas -- may be called independently of update
virtual void replica_share();
// Store the last set for shared ABF
colvar_grid_gradient *last_gradients;
colvar_grid_count *last_samples;
// For Tcl implementation of selection rules.
/// Give the total number of bins for a given bias.
virtual int bin_num();
/// Calculate the bin index for a given bias.
virtual int current_bin();
//// Give the count at a given bin index.
virtual int bin_count(int bin_index);
/// Write human-readable FE gradients and sample count
void write_gradients_samples (const std::string &prefix, bool append = false);
void write_last_gradients_samples (const std::string &prefix, bool append = false);
/// Read human-readable FE gradients and sample count (if not using restart)
void read_gradients_samples ();

1
lib/colvars/colvarbias_meta.cpp
View File

@ -723,6 +723,7 @@ void colvarbias_meta::calc_hills_force (size_t const &i,
}
break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
for (h = h_first; h != h_last; h++) {
if (h->value() == 0.0) continue;

116
lib/colvars/colvarcomp.h
View File

@ -219,14 +219,6 @@ public:
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Return a positive number if x2>x1, zero if x2==x1,
/// negative otherwise (can be redefined to transparently implement
/// constraints, symmetries and periodicities) \b Note: \b it \b
/// only \b works \b with \b scalar \b variables, otherwise raises
/// an error
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Wrapp value (for periodic/symmetric cvcs)
virtual void wrap (colvarvalue &x) const;
@ -300,18 +292,6 @@ inline colvarvalue colvar::cvc::dist2_rgrad (colvarvalue const &x1,
return x2.dist2_grad (x1);
}
inline cvm::real colvar::cvc::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
if (this->type() == colvarvalue::type_scalar) {
return cvm::real (x1 - x2);
} else {
cvm::error ("Error: you requested an operation which requires "
"comparison between two non-scalar values.\n");
return 0.0;
}
}
inline void colvar::cvc::wrap (colvarvalue &x) const
{
return;
@ -356,8 +336,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -383,9 +361,6 @@ public:
/// Redefined to handle the box periodicity
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the box periodicity
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -408,8 +383,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -453,8 +426,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Redefined to make use of the user-provided period
virtual void wrap (colvarvalue &x) const;
};
@ -485,8 +456,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -511,8 +480,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -541,8 +508,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -565,8 +530,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -592,8 +555,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -635,8 +596,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -685,8 +644,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -736,9 +693,6 @@ public:
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
virtual void wrap (colvarvalue &x) const;
};
@ -796,8 +750,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
/// \brief Colvar component: self-coordination number within a group
@ -835,8 +787,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
/// \brief Colvar component: hydrogen bond, defined as the product of
@ -872,8 +822,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -915,8 +863,6 @@ public:
// colvarvalue const &x2) const;
// virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
// colvarvalue const &x2) const;
// virtual cvm::real compare (colvarvalue const &x1,
// colvarvalue const &x2) const;
// };
@ -958,8 +904,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
/// \brief Colvar component: dihedPC
@ -988,8 +932,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
/// \brief Colvar component: orientation in space of an atom group,
@ -1030,8 +972,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -1055,8 +995,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -1080,8 +1018,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -1108,8 +1044,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -1141,9 +1075,6 @@ public:
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
/// Redefined to handle the 2*PI periodicity
virtual void wrap (colvarvalue &x) const;
};
@ -1180,8 +1111,6 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad (colvarvalue const &x1,
colvarvalue const &x2) const;
virtual cvm::real compare (colvarvalue const &x1,
colvarvalue const &x2) const;
};
@ -1214,12 +1143,6 @@ public:
return this->dist2_lgrad (x2, x1); \
} \
\
inline cvm::real colvar::TYPE::compare (colvarvalue const &x1, \
colvarvalue const &x2) const \
{ \
return this->dist2_lgrad (x1, x2); \
} \
\
simple_scalar_dist_functions (distance)
// NOTE: distance_z has explicit functions, see below
@ -1268,12 +1191,6 @@ inline colvarvalue colvar::dihedral::dist2_rgrad (colvarvalue const &x1,
return (-2.0) * diff;
}
inline cvm::real colvar::dihedral::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
return dist2_lgrad (x1, x2);
}
inline void colvar::dihedral::wrap (colvarvalue &x) const
{
if ((x.real_value - wrap_center) >= 180.0) {
@ -1313,12 +1230,6 @@ inline colvarvalue colvar::spin_angle::dist2_rgrad (colvarvalue const &x1,
return (-2.0) * diff;
}
inline cvm::real colvar::spin_angle::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
return dist2_lgrad (x1, x2);
}
inline void colvar::spin_angle::wrap (colvarvalue &x) const
{
if ((x.real_value - wrap_center) >= 180.0) {
@ -1370,12 +1281,6 @@ inline colvarvalue colvar::distance_z::dist2_rgrad (colvarvalue const &x1,
return (-2.0) * diff;
}
inline cvm::real colvar::distance_z::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
return dist2_lgrad (x1, x2);
}
inline void colvar::distance_z::wrap (colvarvalue &x) const
{
if (! this->b_periodic) {
@ -1412,13 +1317,6 @@ inline colvarvalue colvar::distance_vec::dist2_rgrad (colvarvalue const &x1,
return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value);
}
inline cvm::real colvar::distance_vec::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
cvm::error ("Error: cannot compare() two distance vectors.\n");
return 0.0;
}
inline cvm::real colvar::distance_dir::dist2 (colvarvalue const &x1,
colvarvalue const &x2) const
{
@ -1437,13 +1335,6 @@ inline colvarvalue colvar::distance_dir::dist2_rgrad (colvarvalue const &x1,
return colvarvalue ((x2.rvector_value - x1.rvector_value), colvarvalue::type_unitvector);
}
inline cvm::real colvar::distance_dir::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
cvm::error ("Error: cannot compare() two distance directions.\n");
return 0.0;
}
// distance between quaternions
inline cvm::real colvar::orientation::dist2 (colvarvalue const &x1,
@ -1464,12 +1355,5 @@ inline colvarvalue colvar::orientation::dist2_rgrad (colvarvalue const &x1,
return x2.quaternion_value.dist2_grad (x1);
}
inline cvm::real colvar::orientation::compare (colvarvalue const &x1,
colvarvalue const &x2) const
{
cvm::error ("Error: cannot compare() two quaternions.\n");
return 0.0;
}
#endif

64
lib/colvars/colvargrid.h
View File

@ -362,6 +362,70 @@ public:
has_data = true;
}
/// \brief Get the change from this to other_grid
/// and store the result in this.
/// this_grid := other_grid - this_grid
/// Grids must have the same dimensions.
void delta_grid (colvar_grid<T> const &other_grid)
{
if (other_grid.multiplicity() != this->multiplicity()) {
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 "
"different size.\n");
return;
}
for (size_t i = 0; i < data.size(); i++) {
data[i] = other_grid.data[i] - data[i];
}
has_data = true;
}
/// \brief Copy data from another grid of the same type, AND
/// identical definition (boundaries, widths)
/// Added for shared ABF.
void copy_grid (colvar_grid<T> const &other_grid)
{
if (other_grid.multiplicity() != this->multiplicity()) {
cvm::error ("Error: trying to copy two grids with "
"different multiplicity.\n");
return;
}
if (other_grid.data.size() != this->data.size()) {
cvm::error ("Error: trying to copy two grids with "
"different size.\n");
return;
}
for (size_t i = 0; i < data.size(); i++) {
data[i] = other_grid.data[i];
}
has_data = true;
}
/// \brief Extract the grid data as they are represented in memory.
/// Put the results in "out_data".
void raw_data_out (T* out_data) const
{
for (size_t i = 0; i < data.size(); i++) out_data[i] = data[i];
}
/// \brief Input the data as they are represented in memory.
void raw_data_in (const T* in_data)
{
for (size_t i = 0; i < data.size(); i++) data[i] = in_data[i];
has_data = true;
}
/// \brief Size of the data as they are represented in memory.
size_t raw_data_num() const { return data.size(); }
/// \brief Get the binned value indexed by ix, or the first of them
/// if the multiplicity is larger than 1

120
lib/colvars/colvarmodule.cpp
View File

@ -435,7 +435,6 @@ std::string colvarmodule::read_colvar(std::string const &name)
return ss.str();
}
cvm::real colvarmodule::energy_difference (std::string const &bias_name,
std::string const &conf)
{
@ -449,6 +448,125 @@ cvm::real colvarmodule::energy_difference (std::string const &bias_name,
return energy_diff;
}
// For shared ABF
cvm::real colvarmodule::read_width(std::string const &name)
{
cvm::increase_depth();
int found = 0;
real width = -1.0;
for (std::vector<colvar *>::iterator cvi = colvars.begin();
cvi != colvars.end();
cvi++) {
if ( (*cvi)->name == name ) {
++found;
width = (*cvi)->width;
}
}
if (found < 1) {
cvm::error ("Error: bias not found.\n");
} else if (found > 1) {
cvm::error ("Error: duplicate bias name.\n");
} else {
cvm::decrease_depth();
}
return width;
}
size_t colvarmodule::bias_current_bin (std::string const &bias_name)
{
cvm::increase_depth();
int found = 0;
size_t ret = 0; // N.B.: size_t is unsigned, so returning -1 would be a problem.
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ( (*bi)->name == bias_name ) {
++found;
ret = (*bi)->current_bin ();
}
}
if (found < 1) {
cvm::error ("Error: bias not found.\n");
} else if (found > 1) {
cvm::error ("Error: duplicate bias name.\n");
} else {
cvm::decrease_depth();
}
return ret;
}
size_t colvarmodule::bias_bin_num (std::string const &bias_name)
{
cvm::increase_depth();
int found = 0;
size_t ret = 0; // N.B.: size_t is unsigned, so returning -1 would be a problem.
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ( (*bi)->name == bias_name ) {
++found;
ret = (*bi)->bin_num ();
}
}
if (found < 1) {
cvm::error ("Error: bias not found.\n");
} else if (found > 1) {
cvm::error ("Error: duplicate bias name.\n");
} else {
cvm::decrease_depth();
}
return ret;
}
size_t colvarmodule::bias_bin_count (std::string const &bias_name, size_t bin_index)
{
cvm::increase_depth();
int found = 0;
size_t ret = 0; // N.B.: size_t is unsigned, so returning -1 would be a problem.
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ( (*bi)->name == bias_name ) {
++found;
ret = (*bi)->bin_count (bin_index);
}
}
if (found < 1) {
cvm::error ("Error: bias not found.\n");
} else if (found > 1) {
cvm::error ("Error: duplicate bias name.\n");
} else {
cvm::decrease_depth();
}
return ret;
}
void colvarmodule::bias_share (std::string const &bias_name)
{
cvm::increase_depth();
int found = 0;
for (std::vector<colvarbias *>::iterator bi = biases.begin();
bi != biases.end();
bi++) {
if ( (*bi)->name == bias_name ) {
++found;
(*bi)->replica_share ();
}
}
if (found < 1) {
cvm::error ("Error: bias not found.\n");
return;
}
if (found > 1) {
cvm::error ("Error: duplicate bias name.\n");
return;
}
cvm::decrease_depth();
}
int colvarmodule::calc() {

41
lib/colvars/colvarmodule.h
View File

@ -4,7 +4,7 @@
#define COLVARMODULE_H
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2014-10-13"
#define COLVARS_VERSION "2014-10-21"
#endif
#ifndef COLVARS_DEBUG
@ -273,6 +273,17 @@ public:
/// currently works for harmonic (force constant and/or centers)
real energy_difference (std::string const &bias_name, std::string const &conf);
/// Give the bin width in the units of the colvar.
real read_width(std::string const &name);
/// Give the total number of bins for a given bias.
size_t bias_bin_num(std::string const &bias_name);
/// Calculate the bin index for a given bias.
size_t bias_current_bin(std::string const &bias_name);
//// Give the count at a given bin index.
size_t bias_bin_count(std::string const &bias_name, size_t bin_index);
//// Share among replicas.
void bias_share(std::string const &bias_name);
/// Calculate collective variables and biases
int calc();
@ -349,6 +360,13 @@ public:
/// Print a message to the main log and exit normally
static void exit (std::string const &message);
// Replica exchange commands.
static bool replica_enabled();
static int replica_index();
static int replica_num();
static void replica_comm_barrier();
static int replica_comm_recv(char* msg_data, int buf_len, int src_rep);
static int replica_comm_send(char* msg_data, int msg_len, int dest_rep);
/// \brief Get the distance between two atomic positions with pbcs handled
/// correctly
@ -536,6 +554,27 @@ inline cvm::real cvm::dt()
return proxy->dt();
}
// Replica exchange commands
inline bool cvm::replica_enabled() {
return proxy->replica_enabled();
}
inline int cvm::replica_index() {
return proxy->replica_index();
}
inline int cvm::replica_num() {
return proxy->replica_num();
}
inline void cvm::replica_comm_barrier() {
return proxy->replica_comm_barrier();
}
inline int cvm::replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
return proxy->replica_comm_recv(msg_data,buf_len,src_rep);
}
inline int cvm::replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
return proxy->replica_comm_send(msg_data,msg_len,dest_rep);
}
inline void cvm::request_system_force()
{
proxy->request_system_force (true);

30
lib/colvars/colvarproxy.h
View File

@ -5,7 +5,7 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2014-09-19"
#define COLVARPROXY_VERSION "2014-10-21"
#endif
@ -68,8 +68,31 @@ public:
virtual int frame (int) { return COLVARS_NOT_IMPLEMENTED; }
// **************** SIMULATION PARAMETERS ****************
// Replica exchange commands:
/// \brief Indicate if multi-replica support is available and active
virtual bool replica_enabled() { return false; }
/// \brief Index of this replica
virtual int replica_index() { return 0; }
/// \brief Total number of replica
virtual int replica_num() { return 1; }
/// \brief Synchronize replica
virtual void replica_comm_barrier() {}
/// \brief Receive data from other replica
virtual int replica_comm_recv(char* msg_data, int buf_len, int src_rep) {
return COLVARS_NOT_IMPLEMENTED;
}
/// \brief Send data to other replica
virtual int replica_comm_send(char* msg_data, int msg_len, int dest_rep) {
return COLVARS_NOT_IMPLEMENTED;
}
// **************** SIMULATION PARAMETERS ****************
/// \brief Prefix to be used for input files (restarts, not
/// configuration)
@ -93,7 +116,7 @@ public:
return output_prefix_str;
}
/// \brief Restarts will be fritten each time this number of steps has passed
/// \brief Restarts will be written each time this number of steps has passed
virtual size_t restart_frequency() = 0;
@ -212,4 +235,3 @@ inline cvm::real colvarproxy::position_dist2 (cvm::atom_pos const &pos1,
}
#endif