jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

Merge pull request #31 from lammps/master 

rebase
This commit is contained in:
Jacob Gissinger 2018-11-05 20:25:38 -07:00 committed by GitHub
parent 20b9c7fd79 2635f7d160
commit ed77701e56
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
24 changed files with 800 additions and 376 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

Binary file not shown.

122
lib/colvars/colvar.cpp
View File

@ -426,8 +426,8 @@ int colvar::init_custom_function(std::string const &conf)
if (x.size() != value_evaluators.size()) {
cvm::error("Error: based on custom function type, expected "
+ cvm::to_str(x.size()) + " scalar expressions, but "
+ cvm::to_str(value_evaluators.size() + " were found.\n"));
+ cvm::to_str(x.size()) + " scalar expressions, but "
+ cvm::to_str(value_evaluators.size()) + " were found.\n");
return INPUT_ERROR;
}
@ -455,36 +455,42 @@ int colvar::init_grid_parameters(std::string const &conf)
}
lower_boundary.type(value());
upper_boundary.type(value());
upper_wall.type(value());
if (is_enabled(f_cv_scalar)) {
if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
enable(f_cv_lower_boundary);
}
std::string lw_conf, uw_conf;
if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) {
cvm::log("Warning: lowerWallConstant and lowerWall are deprecated, "
"please define a harmonicWalls bias instead.\n");
lower_wall.type(value());
get_keyval(conf, "lowerWall", lower_wall, lower_boundary);
lw_conf = std::string("\n\
lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
lowerWalls "+cvm::to_str(lower_wall)+"\n");
}
if (get_keyval(conf, "upperBoundary", upper_boundary, upper_boundary)) {
enable(f_cv_upper_boundary);
}
if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) {
cvm::log("Warning: upperWallConstant and upperWall are deprecated, "
"please define a harmonicWalls bias instead.\n");
std::string lw_conf, uw_conf;
if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0,
parse_silent)) {
cvm::log("Reading legacy options lowerWall and lowerWallConstant: "
"consider using a harmonicWalls restraint.\n");
lower_wall.type(value());
if (!get_keyval(conf, "lowerWall", lower_wall, lower_boundary)) {
cvm::log("Warning: lowerWall will need to be "
"defined explicitly in the next release.\n");
}
lw_conf = std::string("\n\
lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\
lowerWalls "+cvm::to_str(lower_wall)+"\n");
}
if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0,
parse_silent)) {
cvm::log("Reading legacy options upperWall and upperWallConstant: "
"consider using a harmonicWalls restraint.\n");
upper_wall.type(value());
get_keyval(conf, "upperWall", upper_wall, upper_boundary);
if (!get_keyval(conf, "upperWall", upper_wall, upper_boundary)) {
cvm::log("Warning: upperWall will need to be "
"defined explicitly in the next release.\n");
}
uw_conf = std::string("\n\
upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\
upperWalls "+cvm::to_str(upper_wall)+"\n");
@ -677,6 +683,7 @@ template<typename def_class_name> int colvar::init_components_type(std::string c
if (cvcp != NULL) {
cvcs.push_back(cvcp);
cvcp->check_keywords(def_conf, def_config_key);
cvcp->config_key = def_config_key;
if (cvm::get_error()) {
cvm::error("Error: in setting up component \""+
std::string(def_config_key)+"\".\n", INPUT_ERROR);
@ -1022,13 +1029,13 @@ int colvar::calc()
int colvar::calc_cvcs(int first_cvc, size_t num_cvcs)
{
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\", components "+
cvm::to_str(first_cvc)+" through "+cvm::to_str(first_cvc+num_cvcs)+".\n");
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
error_code |= check_cvc_range(first_cvc, num_cvcs);
if (error_code != COLVARS_OK) {
return error_code;
@ -1059,9 +1066,10 @@ int colvar::collect_cvc_data()
if (cvm::debug())
cvm::log("Calculating colvar \""+this->name+"\"'s properties.\n");
colvarproxy *proxy = cvm::main()->proxy;
int error_code = COLVARS_OK;
if (cvm::step_relative() > 0) {
if ((cvm::step_relative() > 0) && (!proxy->total_forces_same_step())){
// Total force depends on Jacobian derivative from previous timestep
// collect_cvc_total_forces() uses the previous value of jd
error_code |= collect_cvc_total_forces();
@ -1069,6 +1077,10 @@ int colvar::collect_cvc_data()
error_code |= collect_cvc_values();
error_code |= collect_cvc_gradients();
error_code |= collect_cvc_Jacobians();
if (proxy->total_forces_same_step()){
// Use Jacobian derivative from this timestep
error_code |= collect_cvc_total_forces();
}
error_code |= calc_colvar_properties();
if (cvm::debug())
@ -1394,6 +1406,13 @@ int colvar::calc_colvar_properties()
vr.reset(); // (already 0; added for clarity)
}
// Special case of a repeated timestep (eg. multiple NAMD "run" statements)
// revert values of the extended coordinate and velocity prior to latest integration
if (cvm::step_relative() == prev_timestep) {
xr = prev_xr;
vr = prev_vr;
}
// report the restraint center as "value"
x_reported = xr;
v_reported = vr;
@ -1450,7 +1469,6 @@ cvm::real colvar::update_forces_energy()
// extended variable if there is one
if (is_enabled(f_cv_extended_Lagrangian)) {
if (cvm::debug()) {
cvm::log("Updating extended-Lagrangian degree of freedom.\n");
}
@ -1502,6 +1520,11 @@ cvm::real colvar::update_forces_energy()
ft_reported = f_ext;
}
// backup in case we need to revert this integration timestep
// if the same MD timestep is re-run
prev_xr = xr;
prev_vr = vr;
// leapfrog: starting from x_i, f_i, v_(i-1/2)
vr += (0.5 * dt) * f_ext / ext_mass;
// Because of leapfrog, kinetic energy at time i is approximate
@ -1638,18 +1661,26 @@ int colvar::set_cvc_flags(std::vector<bool> const &flags)
}
void colvar::update_active_cvc_square_norm()
{
active_cvc_square_norm = 0.0;
for (size_t i = 0; i < cvcs.size(); i++) {
if (cvcs[i]->is_enabled()) {
active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
}
}
}
int colvar::update_cvc_flags()
{
// Update the enabled/disabled status of cvcs if necessary
if (cvc_flags.size()) {
n_active_cvcs = 0;
active_cvc_square_norm = 0.;
for (size_t i = 0; i < cvcs.size(); i++) {
cvcs[i]->set_enabled(f_cvc_active, cvc_flags[i]);
if (cvcs[i]->is_enabled()) {
n_active_cvcs++;
active_cvc_square_norm += cvcs[i]->sup_coeff * cvcs[i]->sup_coeff;
}
}
if (!n_active_cvcs) {
@ -1657,12 +1688,49 @@ int colvar::update_cvc_flags()
return COLVARS_ERROR;
}
cvc_flags.clear();
update_active_cvc_square_norm();
}
return COLVARS_OK;
}
int colvar::update_cvc_config(std::vector<std::string> const &confs)
{
cvm::log("Updating configuration for colvar \""+name+"\"");
if (confs.size() != cvcs.size()) {
return cvm::error("Error: Wrong number of CVC config strings. "
"For those CVCs that are not being changed, try passing "
"an empty string.", INPUT_ERROR);
}
int error_code = COLVARS_OK;
int num_changes = 0;
for (size_t i = 0; i < cvcs.size(); i++) {
if (confs[i].size()) {
std::string conf(confs[i]);
cvm::increase_depth();
error_code |= cvcs[i]->colvar::cvc::init(conf);
error_code |= cvcs[i]->check_keywords(conf,
cvcs[i]->config_key.c_str());
cvm::decrease_depth();
num_changes++;
}
}
if (num_changes == 0) {
cvm::log("Warning: no changes were applied through modifycvcs; "
"please check that its argument is a list of strings.\n");
}
update_active_cvc_square_norm();
return error_code;
}
// ******************** METRIC FUNCTIONS ********************
// Use the metrics defined by \link cvc \endlink objects

20
lib/colvars/colvar.h
View File

@ -171,8 +171,12 @@ protected:
// Options for extended_lagrangian
/// Restraint center
colvarvalue xr;
/// Previous value of the restraint center;
colvarvalue prev_xr;
/// Velocity of the restraint center
colvarvalue vr;
/// Previous velocity of the restraint center
colvarvalue prev_vr;
/// Mass of the restraint center
cvm::real ext_mass;
/// Restraint force constant
@ -352,6 +356,9 @@ public:
/// \brief Updates the flags in the CVC objects, and their number
int update_cvc_flags();
/// \brief Modify the configuration of CVCs (currently, only base class data)
int update_cvc_config(std::vector<std::string> const &confs);
protected:
/// \brief Number of CVC objects with an active flag
size_t n_active_cvcs;
@ -359,10 +366,17 @@ protected:
/// Sum of square coefficients for active cvcs
cvm::real active_cvc_square_norm;
/// Update the sum of square coefficients for active cvcs
void update_active_cvc_square_norm();
/// \brief Absolute timestep number when this colvar was last updated
int prev_timestep;
public:
/// \brief Return the number of CVC objects defined
inline size_t num_cvcs() const { return cvcs.size(); }
/// \brief Return the number of CVC objects with an active flag (as set by update_cvc_flags)
inline size_t num_active_cvcs() const { return n_active_cvcs; }
@ -371,21 +385,21 @@ public:
///
/// Handles correctly symmetries and periodic boundary conditions
cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to calculate square distances and gradients
///
/// Handles correctly symmetries and periodic boundary conditions
colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
colvarvalue const &x2) const;
/// \brief Use the internal metrics (as from \link cvc
/// \endlink objects) to wrap a value into a standard interval

4
lib/colvars/colvaratoms.cpp
View File

@ -21,6 +21,8 @@ cvm::atom::atom()
{
index = -1;
id = -1;
mass = 1.0;
charge = 1.0;
reset_data();
}
@ -395,7 +397,7 @@ int cvm::atom_group::parse(std::string const &group_conf)
}
// NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
parse_error |= cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
}
}

2
lib/colvars/colvaratoms.h
View File

@ -90,7 +90,7 @@ public:
/// Destructor
~atom();
/// Set mutable data (everything except id and mass) to zero; update mass
/// Set mutable data (everything except id and mass) to zero
inline void reset_data()
{
pos = cvm::atom_pos(0.0);

58
lib/colvars/colvarbias_abf.cpp
View File

@ -564,6 +564,8 @@ int colvarbias_abf::replica_share() {
return COLVARS_OK;
}
void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool append)
{
std::string samples_out_name = prefix + ".count";
@ -572,10 +574,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *samples_os =
cvm::proxy->output_stream(samples_out_name, mode);
if (!samples_os) {
cvm::error("Error opening ABF samples file " + samples_out_name + " for writing");
return;
}
if (!samples_os) return;
samples->write_multicol(*samples_os);
cvm::proxy->close_output_stream(samples_out_name);
@ -583,10 +582,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
if (num_variables() > 2) {
std::string samples_dx_out_name = prefix + ".count.dx";
std::ostream *samples_dx_os = cvm::proxy->output_stream(samples_dx_out_name, mode);
if (!samples_os) {
cvm::error("Error opening samples file " + samples_dx_out_name + " for writing");
return;
}
if (!samples_os) return;
samples->write_opendx(*samples_dx_os);
*samples_dx_os << std::endl;
cvm::proxy->close_output_stream(samples_dx_out_name);
@ -594,10 +590,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *gradients_os =
cvm::proxy->output_stream(gradients_out_name, mode);
if (!gradients_os) {
cvm::error("Error opening ABF gradient file " + gradients_out_name + " for writing");
return;
}
if (!gradients_os) return;
gradients->write_multicol(*gradients_os);
cvm::proxy->close_output_stream(gradients_out_name);
@ -609,20 +602,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::string pmf_out_name = prefix + ".pmf";
std::ostream *pmf_os = cvm::proxy->output_stream(pmf_out_name, mode);
if (!pmf_os) {
cvm::error("Error opening pmf file " + pmf_out_name + " for writing");
return;
}
if (!pmf_os) return;
pmf->write_multicol(*pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string pmf_dx_out_name = prefix + ".pmf.dx";
std::ostream *pmf_dx_os = cvm::proxy->output_stream(pmf_dx_out_name, mode);
if (!pmf_dx_os) {
cvm::error("Error opening pmf file " + pmf_dx_out_name + " for writing");
return;
}
if (!pmf_dx_os) return;
pmf->write_opendx(*pmf_dx_os);
*pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(pmf_dx_out_name);
@ -639,10 +626,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *z_samples_os =
cvm::proxy->output_stream(z_samples_out_name, mode);
if (!z_samples_os) {
cvm::error("Error opening eABF z-histogram file " + z_samples_out_name + " for writing");
return;
}
if (!z_samples_os) return;
z_samples->write_multicol(*z_samples_os);
cvm::proxy->close_output_stream(z_samples_out_name);
@ -651,10 +635,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *z_gradients_os =
cvm::proxy->output_stream(z_gradients_out_name, mode);
if (!z_gradients_os) {
cvm::error("Error opening eABF z-gradient file " + z_gradients_out_name + " for writing");
return;
}
if (!z_gradients_os) return;
z_gradients->write_multicol(*z_gradients_os);
cvm::proxy->close_output_stream(z_gradients_out_name);
}
@ -672,10 +653,7 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::ostream *czar_gradients_os =
cvm::proxy->output_stream(czar_gradients_out_name, mode);
if (!czar_gradients_os) {
cvm::error("Error opening CZAR gradient file " + czar_gradients_out_name + " for writing");
return;
}
if (!czar_gradients_os) return;
czar_gradients->write_multicol(*czar_gradients_os);
cvm::proxy->close_output_stream(czar_gradients_out_name);
@ -688,20 +666,14 @@ void colvarbias_abf::write_gradients_samples(const std::string &prefix, bool app
std::string czar_pmf_out_name = prefix + ".czar.pmf";
std::ostream *czar_pmf_os = cvm::proxy->output_stream(czar_pmf_out_name, mode);
if (!czar_pmf_os) {
cvm::error("Error opening CZAR pmf file " + czar_pmf_out_name + " for writing");
return;
}
if (!czar_pmf_os) return;
czar_pmf->write_multicol(*czar_pmf_os);
// In dimension higher than 2, dx is easier to handle and visualize
if (num_variables() > 2) {
std::string czar_pmf_dx_out_name = prefix + ".czar.pmf.dx";
std::ostream *czar_pmf_dx_os = cvm::proxy->output_stream(czar_pmf_dx_out_name, mode);
if (!czar_pmf_dx_os) {
cvm::error("Error opening CZAR pmf file " + czar_pmf_dx_out_name + " for writing");
return;
}
if (!czar_pmf_dx_os) return;
czar_pmf->write_opendx(*czar_pmf_dx_os);
*czar_pmf_dx_os << std::endl;
cvm::proxy->close_output_stream(czar_pmf_dx_out_name);
@ -854,3 +826,9 @@ std::istream & colvarbias_abf::read_state_data(std::istream& is)
return is;
}
int colvarbias_abf::write_output_files()
{
write_gradients_samples(output_prefix);
return COLVARS_OK;
}

6
lib/colvars/colvarbias_abf.h
View File

@ -136,13 +136,13 @@ private:
/// Write human-readable FE gradients and sample count, and DX file in dim > 2
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();
std::istream& read_state_data(std::istream&);
std::ostream& write_state_data(std::ostream&);
virtual std::istream& read_state_data(std::istream&);
virtual std::ostream& write_state_data(std::ostream&);
virtual int write_output_files();
};
#endif

22
lib/colvars/colvarbias_restraint.cpp
View File

@ -996,11 +996,16 @@ int colvarbias_restraint_harmonic_walls::init(std::string const &conf)
if ((lower_walls.size() > 0) && (upper_walls.size() > 0)) {
for (i = 0; i < num_variables(); i++) {
if (lower_walls[i] >= upper_walls[i]) {
cvm::error("Error: one upper wall, "+
cvm::to_str(upper_walls[i])+
", is not higher than the lower wall, "+
cvm::to_str(lower_walls[i])+".\n",
INPUT_ERROR);
return cvm::error("Error: one upper wall, "+
cvm::to_str(upper_walls[i])+
", is not higher than the lower wall, "+
cvm::to_str(lower_walls[i])+".\n",
INPUT_ERROR);
}
if (variables(i)->dist2(lower_walls[i], upper_walls[i]) < 1.0e-12) {
return cvm::error("Error: lower wall and upper wall are equal "
"in the domain of the variable \""+
variables(i)->name+"\".\n", INPUT_ERROR);
}
}
if (lower_wall_k * upper_wall_k == 0.0) {
@ -1279,13 +1284,16 @@ cvm::real colvarbias_restraint_linear::energy_difference(std::string const &conf
cvm::real colvarbias_restraint_linear::restraint_potential(size_t i) const
{
return force_k / variables(i)->width * (variables(i)->value() - colvar_centers[i]);
return force_k / variables(i)->width * (variables(i)->value() -
colvar_centers[i]).sum();
}
colvarvalue const colvarbias_restraint_linear::restraint_force(size_t i) const
{
return -1.0 * force_k / variables(i)->width;
colvarvalue dummy(variables(i)->value());
dummy.set_ones();
return -1.0 * force_k / variables(i)->width * dummy;
}

25
lib/colvars/colvarcomp.cpp
View File

@ -43,16 +43,27 @@ colvar::cvc::cvc(std::string const &conf)
int colvar::cvc::init(std::string const &conf)
{
int error_code = COLVARS_OK;
if (cvm::debug())
cvm::log("Initializing cvc base object.\n");
get_keyval(conf, "name", this->name, this->name);
std::string const old_name(name);
if (name.size() > 0) {
// Temporary description until child object is initialized
description = "cvc " + name;
} else {
description = "uninitialized cvc";
cvm::log("Updating configuration for component \""+name+"\"");
}
if (get_keyval(conf, "name", name, name)) {
if (name.size() > 0) {
description = "cvc \"" + name + "\" of type " + function_type;
} else {
description = "unnamed cvc";
}
if ((name != old_name) && (old_name.size() > 0)) {
cvm::error("Error: cannot rename component \""+old_name+
"\" after initialization (new name = \""+name+"\")",
INPUT_ERROR);
name = old_name;
}
}
get_keyval(conf, "componentCoeff", sup_coeff, sup_coeff);
@ -78,7 +89,7 @@ int colvar::cvc::init(std::string const &conf)
if (cvm::debug())
cvm::log("Done initializing cvc base object.\n");
return error_code;
return cvm::get_error();
}

118
lib/colvars/colvarcomp.h
View File

@ -82,6 +82,9 @@ public:
/// this variable definition should be set within the constructor.
std::string function_type;
/// Keyword used in the input to denote this CVC
std::string config_key;
/// \brief Coefficient in the polynomial combination (default: 1.0)
cvm::real sup_coeff;
/// \brief Exponent in the polynomial combination (default: 1)
@ -834,49 +837,70 @@ protected:
cvm::real r0;
/// \brief "Cutoff vector" for anisotropic calculation
cvm::rvector r0_vec;
/// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be
/// used
/// \brief Whether r/r0 or \vec{r}*\vec{1/r0_vec} should be used
bool b_anisotropic;
/// Integer exponent of the function numerator
int en;
/// Integer exponent of the function denominator
int ed;
/// \brief If true, group2 will be treated as a single atom
/// (default: loop over all pairs of atoms in group1 and group2)
bool b_group2_center_only;
/// \brief If true, group2 will be treated as a single atom, stored in this
/// accessory group
cvm::atom_group *group2_center;
/// Tolerance for the pair list
cvm::real tolerance;
/// Frequency of update of the pair list
int pairlist_freq;
/// Pair list
bool *pairlist;
public:
/// Constructor
coordnum(std::string const &conf);
coordnum();
virtual ~coordnum() {}
~coordnum();
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
/// coordination number \param exp_num \i n exponent \param exp_den
/// \i m exponent \param A1 atom \param A2 atom
static cvm::real switching_function(cvm::real const &r0,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
/// vector of different cutoffs in the three directions \param
/// exp_num \i n exponent \param exp_den \i m exponent \param A1
/// atom \param A2 atom
static cvm::real switching_function(cvm::rvector const &r0_vec,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
virtual colvarvalue dist2_lgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
enum {
ef_null = 0,
ef_gradients = 1,
ef_anisotropic = (1<<8),
ef_use_pairlist = (1<<9),
ef_rebuild_pairlist = (1<<10)
};
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), where x = |A1-A2|/r0 \param r0, r0_vec "cutoff" for
/// the coordination number (scalar or vector depending on user choice)
/// \param en Numerator exponent \param ed Denominator exponent \param First
/// atom \param Second atom \param pairlist_elem pointer to pair flag for
/// this pair \param tolerance A pair is defined as having a larger
/// coordination than this number
template<int flags>
static cvm::real switching_function(cvm::real const &r0,
cvm::rvector const &r0_vec,
int en,
int ed,
cvm::atom &A1,
cvm::atom &A2,
bool **pairlist_elem,
cvm::real tolerance);
/// Main workhorse function
template<int flags> int compute_coordnum();
};
@ -887,7 +911,8 @@ class colvar::selfcoordnum
: public colvar::cvc
{
protected:
/// First atom group
/// Selected atoms
cvm::atom_group *group1;
/// \brief "Cutoff" for isotropic calculation (default)
cvm::real r0;
@ -895,22 +920,18 @@ protected:
int en;
/// Integer exponent of the function denominator
int ed;
cvm::real tolerance;
int pairlist_freq;
bool *pairlist;
public:
/// Constructor
selfcoordnum(std::string const &conf);
selfcoordnum();
virtual ~selfcoordnum() {}
~selfcoordnum();
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
/// coordination number \param exp_num \i n exponent \param exp_den
/// \i m exponent \param A1 atom \param A2 atom
static cvm::real switching_function(cvm::real const &r0,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;
@ -918,6 +939,9 @@ public:
colvarvalue const &x2) const;
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// Main workhorse function
template<int flags> int compute_selfcoordnum();
};
@ -947,26 +971,6 @@ public:
virtual void calc_value();
virtual void calc_gradients();
virtual void apply_force(colvarvalue const &force);
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |A1-A2|/r0 \param r0 "cutoff" for the
/// coordination number \param exp_num \i n exponent \param exp_den
/// \i m exponent \param A1 atom \param A2 atom
static cvm::real switching_function(cvm::real const &r0,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
/*
template<bool b_gradients>
/// \brief Calculate a coordination number through the function
/// (1-x**n)/(1-x**m), x = |(A1-A2)*(r0_vec)^-|1 \param r0_vec
/// vector of different cutoffs in the three directions \param
/// exp_num \i n exponent \param exp_den \i m exponent \param A1
/// atom \param A2 atom
static cvm::real switching_function(cvm::rvector const &r0_vec,
int const &exp_num, int const &exp_den,
cvm::atom &A1, cvm::atom &A2);
*/
virtual cvm::real dist2(colvarvalue const &x1,
colvarvalue const &x2) const;

521
lib/colvars/colvarcomp_coordnums.cpp
View File

@ -18,45 +18,36 @@
template<bool calculate_gradients>
template<int flags>
cvm::real colvar::coordnum::switching_function(cvm::real const &r0,
int const &en,
int const &ed,
cvm::rvector const &r0_vec,
int en,
int ed,
cvm::atom &A1,
cvm::atom &A2)
cvm::atom &A2,
bool **pairlist_elem,
cvm::real pairlist_tol)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::real const l2 = diff.norm2()/(r0*r0);
// Assume en and ed are even integers, and avoid sqrt in the following
int const en2 = en/2;
int const ed2 = ed/2;
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
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*r0))*diff;
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
if ((flags & ef_use_pairlist) && !(flags & ef_rebuild_pairlist)) {
bool const within = **pairlist_elem;
(*pairlist_elem)++;
if (!within) {
return 0.0;
}
}
return func;
}
cvm::rvector const r0sq_vec(r0_vec.x*r0_vec.x,
r0_vec.y*r0_vec.y,
r0_vec.z*r0_vec.z);
template<bool calculate_gradients>
cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
int const &en,
int const &ed,
cvm::atom &A1,
cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
cvm::rvector const scal_diff(diff.x/((flags & ef_anisotropic) ?
r0_vec.x : r0),
diff.y/((flags & ef_anisotropic) ?
r0_vec.y : r0),
diff.z/((flags & ef_anisotropic) ?
r0_vec.z : r0));
cvm::real const l2 = scal_diff.norm2();
// Assume en and ed are even integers, and avoid sqrt in the following
@ -65,22 +56,45 @@ cvm::real colvar::coordnum::switching_function(cvm::rvector const &r0_vec,
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
//The subtraction and division stretches the function back to the range of [0,1] from [pairlist_tol,1]
cvm::real const func = (((1.0-xn)/(1.0-xd)) - pairlist_tol) / (1.0-pairlist_tol);
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);
if (flags & ef_rebuild_pairlist) {
//Particles just outside of the cutoff also are considered if they come near.
**pairlist_elem = (func > (-pairlist_tol * 0.5)) ? true : false;
(*pairlist_elem)++;
}
//If the value is too small, we need to exclude it, rather than let it contribute to the sum or the gradients.
if (func < 0)
return 0;
if (flags & ef_gradients) {
//This is the old, completely correct expression for dFdl2:
//cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2) -
// func*ed2*(xd/l2))*(-1.0);
//This can become:
//cvm::real const dFdl2 = (1.0/(1.0-xd))*(en2*(xn/l2)*(1.0-xn)/(1.0-xn) -
// func*ed2*(xd/l2))*(-1.0);
//Recognizing that func = (1.0-xn)/(1.0-xd), we can group together the "func" and get a version of dFdl2 that is 0
//when func=0, which lets us skip this gradient calculation when func=0.
cvm::real const dFdl2 = func * ((ed2*xd/((1.0-xd)*l2)) - (en2*xn/((1.0-xn)*l2)));
cvm::rvector const dl2dx((2.0/((flags & ef_anisotropic) ? r0sq_vec.x :
r0*r0)) * diff.x,
(2.0/((flags & ef_anisotropic) ? r0sq_vec.y :
r0*r0)) * diff.y,
(2.0/((flags & ef_anisotropic) ? r0sq_vec.z :
r0*r0)) * diff.z);
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
}
return func;
}
colvar::coordnum::coordnum(std::string const &conf)
: cvc(conf), b_anisotropic(false), b_group2_center_only(false)
: cvc(conf), b_anisotropic(false), group2_center(NULL), pairlist(NULL)
{
function_type = "coordnum";
x.type(colvarvalue::type_scalar);
@ -90,23 +104,26 @@ colvar::coordnum::coordnum(std::string const &conf)
if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) {
cvm::error("Error: group1 and group2 share a common atom (number: " +
cvm::to_str(atom_number) + ")\n");
cvm::to_str(atom_number) + ")\n", INPUT_ERROR);
return;
}
if (group1->b_dummy) {
cvm::error("Error: only group2 is allowed to be a dummy atom\n");
cvm::error("Error: only group2 is allowed to be a dummy atom\n",
INPUT_ERROR);
return;
}
bool const b_isotropic = get_keyval(conf, "cutoff", r0,
cvm::real(4.0 * cvm::unit_angstrom()));
if (get_keyval(conf, "cutoff3", r0_vec, cvm::rvector(4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom()))) {
if (get_keyval(conf, "cutoff3", r0_vec,
cvm::rvector(4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom(),
4.0 * cvm::unit_angstrom()))) {
if (b_isotropic) {
cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n",
cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" "
"at the same time.\n",
INPUT_ERROR);
return;
}
@ -135,86 +152,178 @@ colvar::coordnum::coordnum(std::string const &conf)
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
bool b_group2_center_only = false;
get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy);
if (b_group2_center_only) {
if (!group2_center) {
group2_center = new cvm::atom_group();
group2_center->add_atom(cvm::atom());
}
}
get_keyval(conf, "tolerance", tolerance, 0.0);
if (tolerance > 0) {
get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
if ( ! (pairlist_freq > 0) ) {
cvm::error("Error: non-positive pairlistfrequency provided.\n",
INPUT_ERROR);
return; // and do not allocate the pairlists below
}
if (b_group2_center_only) {
pairlist = new bool[group1->size()];
}
else {
pairlist = new bool[group1->size() * group2->size()];
}
}
}
colvar::coordnum::coordnum()
: b_anisotropic(false), b_group2_center_only(false)
: b_anisotropic(false), group2_center(NULL), pairlist(NULL)
{
function_type = "coordnum";
x.type(colvarvalue::type_scalar);
}
void colvar::coordnum::calc_value()
colvar::coordnum::~coordnum()
{
x.real_value = 0.0;
if (pairlist != NULL) {
delete [] pairlist;
}
if (group2_center != NULL) {
delete group2_center;
}
}
if (b_group2_center_only) {
// create a fake atom to hold the group2 com coordinates
cvm::atom group2_com_atom;
group2_com_atom.pos = group2->center_of_mass();
template<int compute_flags> int colvar::coordnum::compute_coordnum()
{
if (group2_center) {
(*group2_center)[0].pos = group2->center_of_mass();
group2_center->calc_required_properties();
}
cvm::atom_group *group2p = group2_center ? group2_center : group2;
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);
} 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);
bool const use_pairlist = (pairlist != NULL);
bool const rebuild_pairlist = (pairlist != NULL) &&
(cvm::step_relative() % pairlist_freq == 0);
bool *pairlist_elem = use_pairlist ? pairlist : NULL;
cvm::atom_iter ai1 = group1->begin(), ai2 = group2p->begin();
cvm::atom_iter const ai1_end = group1->end();
cvm::atom_iter const ai2_end = group2p->end();
if (b_anisotropic) {
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist |
ef_rebuild_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
} else {
int const flags = compute_flags | ef_anisotropic | ef_use_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
}
} else { // if (use_pairlist) {
int const flags = compute_flags | ef_anisotropic;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
NULL, 0.0);
}
}
}
} else {
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);
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
} 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);
} else {
int const flags = compute_flags | ef_use_pairlist;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
&pairlist_elem,
tolerance);
}
}
}
} else { // if (use_pairlist) {
int const flags = compute_flags;
for (ai1 = group1->begin(); ai1 != ai1_end; ai1++) {
for (ai2 = group2->begin(); ai2 != ai2_end; ai2++) {
x.real_value += switching_function<flags>(r0, r0_vec, en, ed,
*ai1, *ai2,
NULL, 0.0);
}
}
}
}
if (compute_flags & ef_gradients) {
if (group2_center) {
group2->set_weighted_gradient((*group2_center)[0].grad);
}
}
return COLVARS_OK;
}
void colvar::coordnum::calc_value()
{
x.real_value = 0.0;
if (is_enabled(f_cvc_gradient)) {
compute_coordnum<ef_gradients>();
} else {
compute_coordnum<ef_null>();
}
}
void colvar::coordnum::calc_gradients()
{
if (b_group2_center_only) {
// create a fake atom to hold the group2 com coordinates
cvm::atom group2_com_atom;
group2_com_atom.pos = group2->center_of_mass();
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);
} else {
for (cvm::atom_iter ai1 = group1->begin(); ai1 != group1->end(); ai1++)
switching_function<true>(r0, en, ed, *ai1, group2_com_atom);
}
group2->set_weighted_gradient(group2_com_atom.grad);
} else {
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);
}
} 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);
}
}
}
// Gradients are computed by calc_value() if f_cvc_gradients is enabled
}
@ -235,7 +344,7 @@ simple_scalar_dist_functions(coordnum)
// h_bond member functions
colvar::h_bond::h_bond(std::string const &conf)
: cvc(conf)
: cvc(conf)
{
if (cvm::debug())
cvm::log("Initializing h_bond object.\n");
@ -307,13 +416,24 @@ colvar::h_bond::~h_bond()
void colvar::h_bond::calc_value()
{
x.real_value = colvar::coordnum::switching_function<false>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
int const flags = coordnum::ef_null;
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
x.real_value =
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*atom_groups[0])[0],
(*atom_groups[0])[1],
NULL, 0.0);
}
void colvar::h_bond::calc_gradients()
{
colvar::coordnum::switching_function<true>(r0, en, ed, (*atom_groups[0])[0], (*atom_groups[0])[1]);
int const flags = coordnum::ef_gradients;
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*atom_groups[0])[0],
(*atom_groups[0])[1],
NULL, 0.0);
}
@ -328,7 +448,7 @@ simple_scalar_dist_functions(h_bond)
colvar::selfcoordnum::selfcoordnum(std::string const &conf)
: cvc(conf)
: cvc(conf), pairlist(NULL)
{
function_type = "selfcoordnum";
x.type(colvarvalue::type_scalar);
@ -353,36 +473,115 @@ colvar::selfcoordnum::selfcoordnum(std::string const &conf)
if (!is_enabled(f_cvc_pbc_minimum_image)) {
cvm::log("Warning: only minimum-image distances are used by this variable.\n");
}
get_keyval(conf, "tolerance", tolerance, 0.0);
if (tolerance > 0) {
get_keyval(conf, "pairListFrequency", pairlist_freq, 100);
if ( ! (pairlist_freq > 0) ) {
cvm::error("Error: non-positive pairlistfrequency provided.\n",
INPUT_ERROR);
return;
}
pairlist = new bool[(group1->size()-1) * (group1->size()-1)];
}
}
colvar::selfcoordnum::selfcoordnum()
: pairlist(NULL)
{
function_type = "selfcoordnum";
x.type(colvarvalue::type_scalar);
}
colvar::selfcoordnum::~selfcoordnum()
{
if (pairlist != NULL) {
delete [] pairlist;
}
}
template<int compute_flags> int colvar::selfcoordnum::compute_selfcoordnum()
{
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
bool const use_pairlist = (pairlist != NULL);
bool const rebuild_pairlist = (pairlist != NULL) &&
(cvm::step_relative() % pairlist_freq == 0);
bool *pairlist_elem = use_pairlist ? pairlist : NULL;
size_t i = 0, j = 0;
size_t const n = group1->size();
// Always isotropic (TODO: enable the ellipsoid?)
if (use_pairlist) {
if (rebuild_pairlist) {
int const flags = compute_flags | coordnum::ef_use_pairlist |
coordnum::ef_rebuild_pairlist;
for (i = 0; i < n - 1; i++) {
for (j = i + 1; j < n; j++) {
x.real_value +=
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*group1)[i],
(*group1)[j],
&pairlist_elem,
tolerance);
}
}
} else {
int const flags = compute_flags | coordnum::ef_use_pairlist;
for (i = 0; i < n - 1; i++) {
for (j = i + 1; j < n; j++) {
x.real_value +=
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*group1)[i],
(*group1)[j],
&pairlist_elem,
tolerance);
}
}
}
} else { // if (use_pairlist) {
int const flags = compute_flags | coordnum::ef_null;
for (i = 0; i < n - 1; i++) {
for (j = i + 1; j < n; j++) {
x.real_value +=
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
(*group1)[i],
(*group1)[j],
&pairlist_elem,
tolerance);
}
}
}
return COLVARS_OK;
}
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]);
}
if (is_enabled(f_cvc_gradient)) {
compute_selfcoordnum<coordnum::ef_gradients>();
} else {
compute_selfcoordnum<coordnum::ef_null>();
}
}
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]);
}
}
// Gradients are computed by calc_value() if f_cvc_gradients is enabled
}
void colvar::selfcoordnum::apply_force(colvarvalue const &force)
{
if (!group1->noforce) {
@ -394,6 +593,7 @@ void colvar::selfcoordnum::apply_force(colvarvalue const &force)
simple_scalar_dist_functions(selfcoordnum)
// groupcoordnum member functions
colvar::groupcoordnum::groupcoordnum(std::string const &conf)
: distance(conf), b_anisotropic(false)
@ -415,7 +615,7 @@ colvar::groupcoordnum::groupcoordnum(std::string const &conf)
if (b_scale) {
cvm::error("Error: cannot specify \"scale\" and "
"\"scale3\" at the same time.\n");
"\"scale3\" at the same time.\n");
return;
}
b_anisotropic = true;
@ -453,95 +653,56 @@ colvar::groupcoordnum::groupcoordnum()
}
template<bool calculate_gradients>
cvm::real colvar::groupcoordnum::switching_function(cvm::real const &r0,
int const &en,
int const &ed,
cvm::atom &A1,
cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::real const l2 = diff.norm2()/(r0*r0);
// Assume en and ed are even integers, and avoid sqrt in the following
int const en2 = en/2;
int const ed2 = ed/2;
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
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*r0))*diff;
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
}
return func;
}
#if 0 // AMG: I don't think there's any reason to support anisotropic,
// and I don't have those flags below in calc_value, but
// if I need them, I'll also need to uncomment this method
template<bool calculate_gradients>
cvm::real colvar::groupcoordnum::switching_function(cvm::rvector const &r0_vec,
int const &en,
int const &ed,
cvm::atom &A1,
cvm::atom &A2)
{
cvm::rvector const diff = cvm::position_distance(A1.pos, A2.pos);
cvm::rvector const scal_diff(diff.x/r0_vec.x, diff.y/r0_vec.y, diff.z/r0_vec.z);
cvm::real const l2 = scal_diff.norm2();
// Assume en and ed are even integers, and avoid sqrt in the following
int const en2 = en/2;
int const ed2 = ed/2;
cvm::real const xn = cvm::integer_power(l2, en2);
cvm::real const xd = cvm::integer_power(l2, ed2);
cvm::real const func = (1.0-xn)/(1.0-xd);
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);
A1.grad += (-1.0)*dFdl2*dl2dx;
A2.grad += dFdl2*dl2dx;
}
return func;
}
#endif
void colvar::groupcoordnum::calc_value()
{
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
// create fake atoms to hold the com coordinates
cvm::atom group1_com_atom;
cvm::atom group2_com_atom;
group1_com_atom.pos = group1->center_of_mass();
group2_com_atom.pos = group2->center_of_mass();
x.real_value = coordnum::switching_function<false>(r0, en, ed,
group1_com_atom, group2_com_atom);
if (b_anisotropic) {
int const flags = coordnum::ef_anisotropic;
x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
} else {
int const flags = coordnum::ef_null;
x.real_value = coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
}
}
void colvar::groupcoordnum::calc_gradients()
{
cvm::rvector const r0_vec(0.0); // TODO enable the flag?
cvm::atom group1_com_atom;
cvm::atom group2_com_atom;
group1_com_atom.pos = group1->center_of_mass();
group2_com_atom.pos = group2->center_of_mass();
coordnum::switching_function<true>(r0, en, ed, group1_com_atom, group2_com_atom);
if (b_anisotropic) {
int const flags = coordnum::ef_gradients | coordnum::ef_anisotropic;
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
} else {
int const flags = coordnum::ef_gradients;
coordnum::switching_function<flags>(r0, r0_vec, en, ed,
group1_com_atom,
group2_com_atom,
NULL, 0.0);
}
group1->set_weighted_gradient(group1_com_atom.grad);
group2->set_weighted_gradient(group2_com_atom.grad);
}

2
lib/colvars/colvarcomp_distances.cpp
View File

@ -939,7 +939,7 @@ colvar::rmsd::rmsd(std::string const &conf)
bool b_Jacobian_derivative = true;
if (atoms->fitting_group != NULL && b_Jacobian_derivative) {
cvm::log("The option \"refPositionsGroup\" (alternative group for fitting) was enabled: "
cvm::log("The option \"fittingGroup\" (alternative group for fitting) was enabled: "
"Jacobian derivatives of the RMSD will not be calculated.\n");
b_Jacobian_derivative = false;
}

31
lib/colvars/colvarmodule.cpp
View File

@ -139,7 +139,7 @@ int colvarmodule::read_config_file(char const *config_filename)
// read the config file into a string
std::string conf = "";
std::string line;
while (colvarparse::getline_nocomments(config_s, line)) {
while (parse->read_config_line(config_s, line)) {
// Delete lines that contain only white space after removing comments
if (line.find_first_not_of(colvarparse::white_space) != std::string::npos)
conf.append(line+"\n");
@ -159,11 +159,12 @@ int colvarmodule::read_config_string(std::string const &config_str)
// strip the comments away
std::string conf = "";
std::string line;
while (colvarparse::getline_nocomments(config_s, line)) {
while (parse->read_config_line(config_s, line)) {
// Delete lines that contain only white space after removing comments
if (line.find_first_not_of(colvarparse::white_space) != std::string::npos)
conf.append(line+"\n");
}
return parse_config(conf);
}
@ -191,6 +192,12 @@ int colvarmodule::parse_config(std::string &conf)
{
extra_conf.clear();
// Check that the input has matching braces
if (colvarparse::check_braces(conf, 0) != COLVARS_OK) {
return cvm::error("Error: unmatched curly braces in configuration.\n",
INPUT_ERROR);
}
// Parse global options
if (catch_input_errors(parse_global_params(conf))) {
return get_error();
@ -235,6 +242,12 @@ int colvarmodule::parse_config(std::string &conf)
}
std::string const & colvarmodule::get_config() const
{
return parse->get_config();
}
int colvarmodule::append_new_config(std::string const &new_conf)
{
extra_conf += new_conf;
@ -246,9 +259,13 @@ int colvarmodule::parse_global_params(std::string const &conf)
{
colvarmodule *cvm = cvm::main();
std::string index_file_name;
if (parse->get_keyval(conf, "indexFile", index_file_name)) {
cvm->read_index_file(index_file_name.c_str());
{
std::string index_file_name;
size_t pos = 0;
while (parse->key_lookup(conf, "indexFile", &index_file_name, &pos)) {
cvm->read_index_file(index_file_name.c_str());
index_file_name.clear();
}
}
if (parse->get_keyval(conf, "smp", proxy->b_smp_active, proxy->b_smp_active)) {
@ -1073,10 +1090,10 @@ colvarmodule::~colvarmodule()
int colvarmodule::reset()
{
parse->init();
cvm::log("Resetting the Collective Variables module.\n");
parse->init();
// Iterate backwards because we are deleting the elements as we go
for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
bi != biases.rend();

8
lib/colvars/colvarmodule.h
View File

@ -131,7 +131,7 @@ public:
/// Module-wide error state
/// see constants at the top of this file
protected:
private:
static int errorCode;
@ -274,6 +274,9 @@ public:
/// \brief Parse a "clean" config string (no comments)
int parse_config(std::string &conf);
/// Get the configuration string read so far (includes comments)
std::string const & get_config() const;
// Parse functions (setup internal data based on a string)
/// Allow reading from Windows text files using using std::getline
@ -296,6 +299,9 @@ public:
private:
/// Configuration string read so far by the module (includes comments)
std::string config_string;
/// Auto-generated configuration during parsing (e.g. to implement
/// back-compatibility)
std::string extra_conf;

55
lib/colvars/colvarparse.cpp
View File

@ -43,9 +43,10 @@ template<typename TYPE> bool colvarparse::_get_keyval_scalar_(std::string const
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
std::istringstream is(data);
@ -98,9 +99,10 @@ bool colvarparse::_get_keyval_scalar_string_(std::string const &conf,
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
std::istringstream is(data);
@ -162,9 +164,10 @@ template<typename TYPE> bool colvarparse::_get_keyval_vector_(std::string const
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
std::istringstream is(data);
@ -319,9 +322,10 @@ bool colvarparse::get_keyval(std::string const &conf,
}
} while (b_found);
if (found_count > 1)
cvm::log("Warning: found more than one instance of \""+
std::string(key)+"\".\n");
if (found_count > 1) {
cvm::error("Error: found more than one instance of \""+
std::string(key)+"\".\n", INPUT_ERROR);
}
if (data.size()) {
if ( (data == std::string("on")) ||
@ -535,6 +539,19 @@ int colvarparse::check_keywords(std::string &conf, char const *key)
}
std::istream & colvarparse::read_config_line(std::istream &is,
std::string &line)
{
cvm::getline(is, line);
config_string += line+'\n';
size_t const comment = line.find('#');
if (comment != std::string::npos) {
line.erase(comment);
}
return is;
}
std::istream & colvarparse::getline_nocomments(std::istream &is,
std::string &line)
{
@ -607,7 +624,7 @@ bool colvarparse::key_lookup(std::string const &conf,
}
// check that there are matching braces between here and the end of conf
bool const b_not_within_block = brace_check(conf, pos);
bool const b_not_within_block = (check_braces(conf, pos) == COLVARS_OK);
bool const b_isolated = (b_isolated_left && b_isolated_right &&
b_not_within_block);
@ -781,19 +798,15 @@ std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb)
}
bool colvarparse::brace_check(std::string const &conf,
int colvarparse::check_braces(std::string const &conf,
size_t const start_pos)
{
size_t brace_count = 0;
int brace_count = 0;
size_t brace = start_pos;
while ( (brace = conf.find_first_of("{}", brace)) != std::string::npos) {
while ((brace = conf.find_first_of("{}", brace)) != std::string::npos) {
if (conf[brace] == '{') brace_count++;
if (conf[brace] == '}') brace_count--;
brace++;
}
if (brace_count != 0)
return false;
else
return true;
return (brace_count != 0) ? INPUT_ERROR : COLVARS_OK;
}

24
lib/colvars/colvarparse.h
View File

@ -24,7 +24,7 @@
/// need to parse input inherit from this
class colvarparse {
private:
protected:
/// \brief List of legal keywords for this object: this is updated
/// by each call to colvarparse::get_keyval() or
@ -47,7 +47,7 @@ private:
/// \brief Remove all the values from the config string
void strip_values(std::string &conf);
/// \brief Configuration string of the object
/// \brief Configuration string of the object (includes comments)
std::string config_string;
public:
@ -72,7 +72,7 @@ public:
}
/// Set a new config string for this object
inline void init(const std::string& conf)
inline void init(std::string const &conf)
{
if (! config_string.size()) {
init();
@ -80,7 +80,8 @@ public:
}
}
inline const std::string& get_config()
/// Get the configuration string (includes comments)
inline std::string const & get_config() const
{
return config_string;
}
@ -284,14 +285,19 @@ public:
std::string *data = NULL,
size_t *save_pos = NULL);
/// \brief Reads a configuration line, adds it to config_string, and returns
/// the stream \param is Input stream \param s String that will hold the
/// configuration line, with comments stripped
std::istream & read_config_line(std::istream &is, std::string &line);
/// \brief Works as std::getline() but also removes everything
/// between a comment character and the following newline
static std::istream & getline_nocomments(std::istream &is,
std::string &s);
static std::istream & getline_nocomments(std::istream &is, std::string &s);
/// Check if the content of the file has matching braces
bool brace_check(std::string const &conf,
size_t const start_pos = 0);
/// \brief Check if the content of a config string has matching braces
/// \param conf The configuration string \param start_pos Start the count
/// from this position
static int check_braces(std::string const &conf, size_t const start_pos);
};

20
lib/colvars/colvarproxy.cpp
View File

@ -289,13 +289,23 @@ colvarproxy_smp::colvarproxy_smp()
omp_lock_state = NULL;
#if defined(_OPENMP)
if (smp_thread_id() == 0) {
omp_lock_state = reinterpret_cast<void *>(new omp_lock_t);
omp_init_lock(reinterpret_cast<omp_lock_t *>(omp_lock_state));
}
#endif
}
colvarproxy_smp::~colvarproxy_smp() {}
colvarproxy_smp::~colvarproxy_smp()
{
#if defined(_OPENMP)
if (smp_thread_id() == 0) {
if (omp_lock_state) {
delete reinterpret_cast<omp_lock_t *>(omp_lock_state);
}
}
#endif
}
int colvarproxy_smp::smp_enabled()
@ -499,6 +509,14 @@ char const *colvarproxy_script::script_obj_to_str(unsigned char *obj)
}
std::vector<std::string> colvarproxy_script::script_obj_to_str_vector(unsigned char *obj)
{
cvm::error("Error: trying to print a script object without a scripting "
"language interface.\n", BUG_ERROR);
return std::vector<std::string>();
}
int colvarproxy_script::run_force_callback()
{
return COLVARS_NOT_IMPLEMENTED;

3
lib/colvars/colvarproxy.h
View File

@ -461,6 +461,9 @@ public:
/// Convert a script object (Tcl or Python call argument) to a C string
virtual char const *script_obj_to_str(unsigned char *obj);
/// Convert a script object (Tcl or Python call argument) to a vector of strings
virtual std::vector<std::string> script_obj_to_str_vector(unsigned char *obj);
/// Pointer to the scripting interface object
/// (does not need to be allocated in a new interface)
colvarscript *script;

2
lib/colvars/colvars_version.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2018-04-29"
#define COLVARS_VERSION "2018-10-16"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars

35
lib/colvars/colvarscript.cpp
View File

@ -82,6 +82,14 @@ int colvarscript::run(int objc, unsigned char *const objv[])
int error_code = COLVARS_OK;
// If command is found in map, execute it
std::string const cmd_key("cv_"+cmd);
if (comm_str_map.count(cmd_key) > 0) {
error_code |= (*(comm_fns[comm_str_map[cmd_key]]))(
reinterpret_cast<void *>(this), objc, objv);
return error_code;
}
if (cmd == "colvar") {
if (objc < 3) {
result = "Missing parameters\n" + help_string();
@ -295,11 +303,12 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
}
if (subcmd == "delete") {
size_t i;
for (i = 0; i < cv->biases.size(); i++) {
while (cv->biases.size() > 0) {
size_t i = cv->biases.size()-1;
cvm::log("Warning: before deleting colvar " + cv->name
+ ", deleting parent bias " + cv->biases[i]->name);
delete cv->biases[i];
}
cv->biases.clear();
// colvar destructor is tasked with the cleanup
delete cv;
// TODO this could be done by the destructors
@ -373,6 +382,23 @@ int colvarscript::proc_colvar(colvar *cv, int objc, unsigned char *const objv[])
return COLVARS_OK;
}
if (subcmd == "modifycvcs") {
if (objc < 4) {
result = "cvcflags: missing parameter: vector of strings";
return COLVARSCRIPT_ERROR;
}
std::vector<std::string> const confs(proxy->script_obj_to_str_vector(objv[3]));
cvm::increase_depth();
int res = cv->update_cvc_config(confs);
cvm::decrease_depth();
if (res != COLVARS_OK) {
result = "Error setting CVC flags";
return COLVARSCRIPT_ERROR;
}
result = "0";
return COLVARS_OK;
}
if ((subcmd == "get") || (subcmd == "set") || (subcmd == "state")) {
return proc_features(cv, objc, objv);
}
@ -547,6 +573,8 @@ std::string colvarscript::help_string() const
Managing the Colvars module:\n\
configfile <file name> -- read configuration from a file\n\
config <string> -- read configuration from the given string\n\
getconfig -- get the module's configuration string\n\
resetindexgroups -- clear the index groups loaded so far\n\
reset -- delete all internal configuration\n\
delete -- delete this Colvars module instance\n\
version -- return version of Colvars code\n\
@ -579,6 +607,7 @@ Accessing collective variables:\n\
colvar <name> gettotalforce -- return total force of colvar <name>\n\
colvar <name> getconfig -- return config string of colvar <name>\n\
colvar <name> cvcflags <fl> -- enable or disable cvcs according to 0/1 flags\n\
colvar <name> modifycvcs <str> -- pass new config strings to each CVC\n\
colvar <name> get <f> -- get the value of the colvar feature <f>\n\
colvar <name> set <f> <val> -- set the value of the colvar feature <f>\n\
\n\

19
lib/colvars/colvarscript.h
View File

@ -71,8 +71,10 @@ public:
cv_help,
cv_version,
cv_config,
cv_getconfig,
cv_configfile,
cv_reset,
cv_resetindexgroups,
cv_delete,
cv_list,
cv_list_biases,
@ -254,6 +256,23 @@ extern "C" {
return COLVARSCRIPT_ERROR;
)
CVSCRIPT(cv_getconfig,
"Get the module's configuration string read so far",
0, 0,
{ },
script->set_str_result(cvm::main()->get_config());
return COLVARS_OK;
)
CVSCRIPT(cv_resetindexgroups,
"Clear the index groups loaded so far, allowing to replace them",
0, 0,
{ },
cvm::main()->index_group_names.clear();
cvm::main()->index_groups.clear();
return COLVARS_OK;
)
CVSCRIPT(cv_addenergy,
"Add an energy to the MD engine",
1, 1,

34
lib/colvars/colvarvalue.cpp
View File

@ -379,6 +379,40 @@ void colvarvalue::set_random()
}
void colvarvalue::set_ones(cvm::real assigned_value)
{
size_t ic;
switch (this->type()) {
case colvarvalue::type_scalar:
this->real_value = assigned_value;
break;
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
this->rvector_value.x = assigned_value;
this->rvector_value.y = assigned_value;
this->rvector_value.z = assigned_value;
break;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
this->quaternion_value.q0 = assigned_value;
this->quaternion_value.q1 = assigned_value;
this->quaternion_value.q2 = assigned_value;
this->quaternion_value.q3 = assigned_value;
break;
case colvarvalue::type_vector:
for (ic = 0; ic < this->vector1d_value.size(); ic++) {
this->vector1d_value[ic] = assigned_value;
}
break;
case colvarvalue::type_notset:
default:
undef_op();
break;
}
}
void colvarvalue::undef_op() const
{
cvm::error("Error: Undefined operation on a colvar of type \""+

43
lib/colvars/colvarvalue.h
View File

@ -187,6 +187,12 @@ public:
return std::sqrt(this->norm2());
}
/// Sum of the components of this colvarvalue (if more than one dimension)
cvm::real sum() const;
/// Return a colvarvalue object of the same type and all components set to 1
colvarvalue ones() const;
/// Square distance between this \link colvarvalue \endlink and another
cvm::real dist2(colvarvalue const &x2) const;
@ -272,17 +278,21 @@ public:
/// Get a single colvarvalue out of elements of the vector
colvarvalue const get_elem(int const i_begin, int const i_end, Type const vt) const;
/// Get a single colvarvalue out of elements of the vector
colvarvalue const get_elem(int const icv) const;
/// Set elements of the vector from a single colvarvalue (uses the rank of x
/// to compute the length)
void set_elem(int const icv, colvarvalue const &x);
/// Set elements of the vector from a single colvarvalue
void set_elem(int const i_begin, int const i_end, colvarvalue const &x);
/// Make each element a random number in N(0,1)
void set_random();
/// Get a single colvarvalue out of elements of the vector
colvarvalue const get_elem(int const icv) const;
/// Set elements of the vector from a single colvarvalue
void set_elem(int const icv, colvarvalue const &x);
/// Make each element equal to the given argument
void set_ones(cvm::real assigned_value = 1.0);
/// Get a scalar number out of an element of the vector
cvm::real operator [] (int const i) const;
@ -683,6 +693,29 @@ inline cvm::real colvarvalue::norm2() const
}
inline cvm::real colvarvalue::sum() const
{
switch (value_type) {
case colvarvalue::type_scalar:
return (this->real_value);
case colvarvalue::type_3vector:
case colvarvalue::type_unit3vector:
case colvarvalue::type_unit3vectorderiv:
return (this->rvector_value).x + (this->rvector_value).y +
(this->rvector_value).z;
case colvarvalue::type_quaternion:
case colvarvalue::type_quaternionderiv:
return (this->quaternion_value).q0 + (this->quaternion_value).q1 +
(this->quaternion_value).q2 + (this->quaternion_value).q3;
case colvarvalue::type_vector:
return (this->vector1d_value).sum();
case colvarvalue::type_notset:
default:
return 0.0;
}
}
inline cvm::real colvarvalue::dist2(colvarvalue const &x2) const
{
colvarvalue::check_types(*this, x2);

2
src/USER-COLVARS/colvarproxy_lammps_version.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2018-04-29"
#define COLVARPROXY_VERSION "2018-08-29"
// This file is part of the Collective Variables module (Colvars).
// The original version of Colvars and its updates are located at:
// https://github.com/colvars/colvars