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

Merge pull request #883 from giacomofiorin/colvars-update 

Fixes to Colvars library (version 2018-04-29)
This commit is contained in:
Steve Plimpton 2018-05-03 11:25:32 -06:00 committed by GitHub
parent 83f5c331cd 0c005f5cb0
commit 35294dafbc
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
21 changed files with 605 additions and 452 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.

53
lib/colvars/colvar.cpp
View File

@ -7,6 +7,10 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <list>
#include <vector>
#include <algorithm>
#include "colvarmodule.h"
#include "colvarvalue.h"
#include "colvarparse.h"
@ -14,15 +18,6 @@
#include "colvarcomp.h"
#include "colvarscript.h"
// used in build_atom_list()
#include <algorithm>
/// Compare two cvcs using their names
/// Used to sort CVC array in scripted coordinates
bool compare(colvar::cvc *i, colvar::cvc *j) {
return i->name < j->name;
}
colvar::colvar()
@ -34,6 +29,15 @@ colvar::colvar()
}
namespace {
/// Compare two cvcs using their names
/// Used to sort CVC array in scripted coordinates
bool compare(colvar::cvc *i, colvar::cvc *j)
{
return i->name < j->name;
}
}
int colvar::init(std::string const &conf)
{
cvm::log("Initializing a new collective variable.\n");
@ -143,6 +147,9 @@ int colvar::init(std::string const &conf)
x.type(cvc_value);
x_reported.type(cvc_value);
}
set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
// If using scripted biases, any colvar may receive bias forces
// and will need its gradient
if (cvm::scripted_forces()) {
@ -198,6 +205,7 @@ int colvar::init(std::string const &conf)
if (is_enabled(f_cv_homogeneous) && cvcs[0]->b_periodic) { // TODO make this a CVC feature
bool b_periodic = true;
period = cvcs[0]->period;
wrap_center = cvcs[0]->wrap_center;
for (i = 1; i < cvcs.size(); i++) {
if (!cvcs[i]->b_periodic || cvcs[i]->period != period) {
b_periodic = false;
@ -211,6 +219,14 @@ int colvar::init(std::string const &conf)
set_enabled(f_cv_periodic, b_periodic);
}
// Allow scripted/custom functions to be defined as periodic
if ( (is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function)) && is_enabled(f_cv_scalar) ) {
if (get_keyval(conf, "period", period, 0.)) {
set_enabled(f_cv_periodic, true);
get_keyval(conf, "wrapAround", wrap_center, 0.);
}
}
// check that cvcs are compatible
for (i = 0; i < cvcs.size(); i++) {
@ -443,8 +459,6 @@ int colvar::init_grid_parameters(std::string const &conf)
upper_boundary.type(value());
upper_wall.type(value());
set_enabled(f_cv_scalar, (value().type() == colvarvalue::type_scalar));
if (is_enabled(f_cv_scalar)) {
if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) {
@ -1503,7 +1517,7 @@ cvm::real colvar::update_forces_energy()
vr += (0.5 * dt) * f_ext / ext_mass;
xr += dt * vr;
xr.apply_constraints();
if (this->is_enabled(f_cv_periodic)) this->wrap(xr);
this->wrap(xr);
}
// Now adding the force on the actual colvar (for those biases that
@ -1714,9 +1728,18 @@ colvarvalue colvar::dist2_rgrad(colvarvalue const &x1,
void colvar::wrap(colvarvalue &x) const
{
if (is_enabled(f_cv_homogeneous)) {
(cvcs[0])->wrap(x);
if ( !is_enabled(f_cv_periodic) ) {
return;
}
if ( is_enabled(f_cv_scripted) || is_enabled(f_cv_custom_function) ) {
// Scripted functions do their own wrapping, as cvcs might not be periodic
cvm::real shift = std::floor((x.real_value - wrap_center) / period + 0.5);
x.real_value -= shift * period;
} else {
cvcs[0]->wrap(x);
}
return;
}
@ -2244,7 +2267,7 @@ void colvar::calc_runave()
runave.type(value().type());
runave.reset();
// first-step operations
// first-step operationsf
if (cvm::debug())
cvm::log("Colvar \""+this->name+

1
lib/colvars/colvar.h
View File

@ -216,6 +216,7 @@ public:
/// Period, if this variable is periodic
cvm::real period;
cvm::real wrap_center;
/// \brief Expand the boundaries of multiples of width, to keep the

59
lib/colvars/colvaratoms.cpp
View File

@ -7,6 +7,10 @@
// If you wish to distribute your changes, please submit them to the
// Colvars repository at GitHub.
#include <list>
#include <vector>
#include <algorithm>
#include "colvarmodule.h"
#include "colvarproxy.h"
#include "colvarparse.h"
@ -436,7 +440,7 @@ int cvm::atom_group::parse(std::string const &group_conf)
}
bool b_print_atom_ids = false;
get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false, colvarparse::parse_silent);
get_keyval(group_conf, "printAtomIDs", b_print_atom_ids, false);
// Calculate all required properties (such as total mass)
setup();
@ -715,13 +719,10 @@ int cvm::atom_group::parse_fitting_options(std::string const &group_conf)
"if provided, must be non-zero.\n", INPUT_ERROR);
return COLVARS_ERROR;
}
} else {
// if not, rely on existing atom indices for the group
group_for_fit->create_sorted_ids();
ref_pos.resize(group_for_fit->size());
}
cvm::load_coords(ref_pos_file.c_str(), ref_pos, group_for_fit->sorted_ids,
ref_pos.resize(group_for_fit->size());
cvm::load_coords(ref_pos_file.c_str(), &ref_pos, group_for_fit,
ref_pos_col, ref_pos_col_value);
}
@ -789,33 +790,39 @@ void cvm::atom_group::do_feature_side_effects(int id)
}
int cvm::atom_group::create_sorted_ids(void)
int cvm::atom_group::create_sorted_ids()
{
// Only do the work if the vector is not yet populated
if (sorted_ids.size())
if (sorted_atoms_ids.size())
return COLVARS_OK;
std::list<int> temp_id_list;
for (cvm::atom_iter ai = this->begin(); ai != this->end(); ai++) {
temp_id_list.push_back(ai->id);
// Sort the internal IDs
std::list<int> sorted_atoms_ids_list;
for (size_t i = 0; i < this->size(); i++) {
sorted_atoms_ids_list.push_back(atoms_ids[i]);
}
temp_id_list.sort();
temp_id_list.unique();
if (temp_id_list.size() != this->size()) {
cvm::error("Error: duplicate atom IDs in atom group? (found " +
cvm::to_str(temp_id_list.size()) +
" unique atom IDs instead of" +
cvm::to_str(this->size()) + ").\n");
return COLVARS_ERROR;
sorted_atoms_ids_list.sort();
sorted_atoms_ids_list.unique();
if (sorted_atoms_ids_list.size() != this->size()) {
return cvm::error("Error: duplicate atom IDs in atom group? (found " +
cvm::to_str(sorted_atoms_ids_list.size()) +
" unique atom IDs instead of " +
cvm::to_str(this->size()) + ").\n", BUG_ERROR);
}
sorted_ids = std::vector<int> (temp_id_list.size());
unsigned int id_i = 0;
std::list<int>::iterator li;
for (li = temp_id_list.begin(); li != temp_id_list.end(); ++li) {
sorted_ids[id_i] = *li;
id_i++;
// Compute map between sorted and unsorted elements
sorted_atoms_ids.resize(this->size());
sorted_atoms_ids_map.resize(this->size());
std::list<int>::iterator lsii = sorted_atoms_ids_list.begin();
size_t ii = 0;
for ( ; ii < this->size(); lsii++, ii++) {
sorted_atoms_ids[ii] = *lsii;
size_t const pos = std::find(atoms_ids.begin(), atoms_ids.end(), *lsii) -
atoms_ids.begin();
sorted_atoms_ids_map[ii] = pos;
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
return COLVARS_OK;
}

38
lib/colvars/colvaratoms.h
View File

@ -227,9 +227,16 @@ protected:
/// \brief Array of atom objects
std::vector<cvm::atom> atoms;
/// \brief Array of atom identifiers for the MD program (0-based)
/// \brief Internal atom IDs for host code
std::vector<int> atoms_ids;
/// Sorted list of internal atom IDs (populated on-demand by
/// create_sorted_ids); used to read coordinate files
std::vector<int> sorted_atoms_ids;
/// Map entries of sorted_atoms_ids onto the original positions in the group
std::vector<int> sorted_atoms_ids_map;
/// \brief Dummy atom position
cvm::atom_pos dummy_atom_pos;
@ -273,19 +280,34 @@ public:
return atoms.size();
}
std::string const print_atom_ids() const;
/// \brief If this option is on, this group merely acts as a wrapper
/// for a fixed position; any calls to atoms within or to
/// functions that return disaggregated data will fail
bool b_dummy;
/// Sorted list of zero-based (internal) atom ids
/// (populated on-demand by create_sorted_ids)
std::vector<int> sorted_ids;
/// Internal atom IDs (populated during initialization)
inline std::vector<int> const &ids() const
{
return atoms_ids;
}
/// Allocates and populates the sorted list of atom ids
int create_sorted_ids(void);
std::string const print_atom_ids() const;
/// Allocates and populates sorted_ids and sorted_ids_map
int create_sorted_ids();
/// Sorted internal atom IDs (populated on-demand by create_sorted_ids);
/// used to read coordinate files
inline std::vector<int> const &sorted_ids() const
{
return sorted_atoms_ids;
}
/// Map entries of sorted_atoms_ids onto the original positions in the group
inline std::vector<int> const &sorted_ids_map() const
{
return sorted_atoms_ids_map;
}
/// Detect whether two groups share atoms
/// If yes, returns 1-based number of a common atom; else, returns 0

17
lib/colvars/colvarbias_abf.cpp
View File

@ -17,17 +17,17 @@ colvarbias_abf::colvarbias_abf(char const *key)
: colvarbias(key),
b_UI_estimator(false),
b_CZAR_estimator(false),
pabf_freq(0),
system_force(NULL),
gradients(NULL),
pmf(NULL),
samples(NULL),
pmf(NULL),
z_gradients(NULL),
z_samples(NULL),
czar_gradients(NULL),
czar_pmf(NULL),
last_gradients(NULL),
last_samples(NULL),
pabf_freq(0)
last_samples(NULL)
{
}
@ -315,8 +315,6 @@ colvarbias_abf::~colvarbias_abf()
int colvarbias_abf::update()
{
int iter;
if (cvm::debug()) cvm::log("Updating ABF bias " + this->name);
size_t i;
@ -368,7 +366,12 @@ int colvarbias_abf::update()
if ( b_integrate ) {
if ( pabf_freq && cvm::step_relative() % pabf_freq == 0 ) {
cvm::real err;
iter = pmf->integrate(integrate_steps, integrate_tol, err);
int iter = pmf->integrate(integrate_steps, integrate_tol, err);
if ( iter == integrate_steps ) {
cvm::log("Warning: PMF integration did not converge to " + cvm::to_str(integrate_tol)
+ " in " + cvm::to_str(integrate_steps)
+ " steps. Residual error: " + cvm::to_str(err));
}
pmf->set_zero_minimum(); // TODO: do this only when necessary
}
}
@ -485,7 +488,6 @@ int colvarbias_abf::update()
int colvarbias_abf::replica_share() {
int p;
if ( !cvm::replica_enabled() ) {
cvm::error("Error: shared ABF: No replicas.\n");
@ -507,6 +509,7 @@ int colvarbias_abf::replica_share() {
char* msg_data = new char[msg_total];
if (cvm::replica_index() == 0) {
int p;
// Replica 0 collects the delta gradient and count from the others.
for (p = 1; p < cvm::replica_num(); p++) {
// Receive the deltas.

46
lib/colvars/colvarcomp.cpp
View File

@ -21,6 +21,9 @@ colvar::cvc::cvc()
b_try_scalable(true)
{
init_cvc_requires();
sup_coeff = 1.0;
period = 0.0;
wrap_center = 0.0;
}
@ -30,40 +33,47 @@ colvar::cvc::cvc(std::string const &conf)
b_periodic(false),
b_try_scalable(true)
{
init_cvc_requires();
sup_coeff = 1.0;
period = 0.0;
wrap_center = 0.0;
init(conf);
}
int colvar::cvc::init(std::string const &conf)
{
int error_code = COLVARS_OK;
if (cvm::debug())
cvm::log("Initializing cvc base object.\n");
init_cvc_requires();
if (get_keyval(conf, "name", this->name, std::string(""), parse_silent)) {
get_keyval(conf, "name", this->name, this->name);
if (name.size() > 0) {
// Temporary description until child object is initialized
description = "cvc " + name;
} else {
description = "uninitialized cvc";
}
get_keyval(conf, "componentCoeff", sup_coeff, 1.0);
get_keyval(conf, "componentExp", sup_np, 1);
get_keyval(conf, "componentCoeff", sup_coeff, sup_coeff);
get_keyval(conf, "componentExp", sup_np, sup_np);
get_keyval(conf, "period", period, 0.0);
get_keyval(conf, "wrapAround", wrap_center, 0.0);
get_keyval(conf, "period", period, period);
get_keyval(conf, "wrapAround", wrap_center, wrap_center);
get_keyval_feature((colvarparse *)this, conf, "debugGradients",
get_keyval_feature(dynamic_cast<colvarparse *>(this), conf, "debugGradients",
f_cvc_debug_gradient, false, parse_silent);
{
bool b_no_PBC = false;
get_keyval(conf, "forceNoPBC", b_no_PBC, false);
if (b_no_PBC) {
disable(f_cvc_pbc_minimum_image);
} else {
enable(f_cvc_pbc_minimum_image);
}
// this does not use get_keyval_feature() only for backward compatibility
bool b_no_PBC = !is_enabled(f_cvc_pbc_minimum_image); // Enabled by default
get_keyval(conf, "forceNoPBC", b_no_PBC, b_no_PBC);
if (b_no_PBC) {
disable(f_cvc_pbc_minimum_image);
} else {
enable(f_cvc_pbc_minimum_image);
}
// Attempt scalable calculations when in parallel? (By default yes, if available)
get_keyval(conf, "scalable", b_try_scalable, true);
get_keyval(conf, "scalable", b_try_scalable, b_try_scalable);
if (cvm::debug())
cvm::log("Done initializing cvc base object.\n");

12
lib/colvars/colvarcomp.h
View File

@ -98,12 +98,14 @@ public:
/// \brief Constructor
///
/// At least one constructor which reads a string should be defined
/// for every class inheriting from cvc \param conf Contents
/// of the configuration file pertaining to this \link cvc
/// \endlink
/// Calls the init() function of the class
cvc(std::string const &conf);
/// An init function should be defined for every class inheriting from cvc
/// \param conf Contents of the configuration file pertaining to this \link
/// cvc \endlink
virtual int init(std::string const &conf);
/// \brief Within the constructor, make a group parse its own
/// options from the provided configuration string
/// Returns reference to new group
@ -231,7 +233,7 @@ public:
virtual colvarvalue dist2_rgrad(colvarvalue const &x1,
colvarvalue const &x2) const;
/// \brief Wrapp value (for periodic/symmetric cvcs)
/// \brief Wrap value (for periodic/symmetric cvcs)
virtual void wrap(colvarvalue &x) const;
/// \brief Pointers to all atom groups, to let colvars collect info

30
lib/colvars/colvarcomp_distances.cpp
View File

@ -148,7 +148,7 @@ void colvar::distance_vec::apply_force(colvarvalue const &force)
cvm::real colvar::distance_vec::dist2(colvarvalue const &x1,
colvarvalue const &x2) const
{
return cvm::position_dist2(x1.rvector_value, x2.rvector_value);
return (cvm::position_distance(x1.rvector_value, x2.rvector_value)).norm2();
}
@ -192,7 +192,7 @@ colvar::distance_z::distance_z(std::string const &conf)
// this group is optional
ref2 = parse_group(conf, "ref2", true);
if (ref2 && ref2->size()) {
if ( ref2 ) {
cvm::log("Using axis joining the centers of mass of groups \"ref\" and \"ref2\"");
fixed_axis = false;
if (key_lookup(conf, "axis"))
@ -306,7 +306,7 @@ void colvar::distance_z::apply_force(colvarvalue const &force)
if (!ref1->noforce)
ref1->apply_colvar_force(force.real_value);
if (ref2 && ref2->size() && !ref2->noforce)
if (ref2 && !ref2->noforce)
ref2->apply_colvar_force(force.real_value);
if (!main->noforce)
@ -464,7 +464,7 @@ void colvar::distance_xy::apply_force(colvarvalue const &force)
if (!ref1->noforce)
ref1->apply_colvar_force(force.real_value);
if (ref2 && ref2->size() && !ref2->noforce)
if (ref2 && !ref2->noforce)
ref2->apply_colvar_force(force.real_value);
if (!main->noforce)
@ -979,14 +979,12 @@ colvar::rmsd::rmsd(std::string const &conf)
"if provided, must be non-zero.\n");
return;
}
} else {
// if not, rely on existing atom indices for the group
atoms->create_sorted_ids();
ref_pos.resize(atoms->size());
}
cvm::load_coords(ref_pos_file.c_str(), ref_pos, atoms->sorted_ids,
ref_pos_col, ref_pos_col_value);
ref_pos.resize(atoms->size());
cvm::load_coords(ref_pos_file.c_str(), &ref_pos, atoms,
ref_pos_col, ref_pos_col_value);
}
}
@ -1172,13 +1170,11 @@ colvar::eigenvector::eigenvector(std::string const &conf)
"if provided, must be non-zero.\n");
return;
}
} else {
// if not, use atom indices
atoms->create_sorted_ids();
}
ref_pos.resize(atoms->size());
cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value);
cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
file_col, file_col_value);
}
}
@ -1249,13 +1245,11 @@ colvar::eigenvector::eigenvector(std::string const &conf)
cvm::error("Error: vectorColValue, if provided, must be non-zero.\n");
return;
}
} else {
// if not, use atom indices
atoms->create_sorted_ids();
}
eigenvec.resize(atoms->size());
cvm::load_coords(file_name.c_str(), eigenvec, atoms->sorted_ids, file_col, file_col_value);
cvm::load_coords(file_name.c_str(), &eigenvec, atoms,
file_col, file_col_value);
}
}

7
lib/colvars/colvarcomp_rotations.cpp
View File

@ -50,12 +50,11 @@ colvar::orientation::orientation(std::string const &conf)
"if provided, must be non-zero.\n");
return;
}
} else {
// if not, use atom indices
atoms->create_sorted_ids();
}
ref_pos.resize(atoms->size());
cvm::load_coords(file_name.c_str(), ref_pos, atoms->sorted_ids, file_col, file_col_value);
cvm::load_coords(file_name.c_str(), &ref_pos, atoms,
file_col, file_col_value);
}
}

2
lib/colvars/colvardeps.cpp
View File

@ -549,7 +549,7 @@ void colvardeps::init_cv_requires() {
f_req_exclude(f_cv_custom_function, f_cv_scripted);
init_feature(f_cv_periodic, "periodic", f_type_static);
f_req_self(f_cv_periodic, f_cv_homogeneous);
f_req_self(f_cv_periodic, f_cv_scalar);
init_feature(f_cv_scalar, "scalar", f_type_static);
init_feature(f_cv_linear, "linear", f_type_static);
init_feature(f_cv_homogeneous, "homogeneous", f_type_static);

15
lib/colvars/colvargrid.cpp
View File

@ -329,7 +329,6 @@ void integrate_potential::update_div_local(const std::vector<int> &ix0)
const int linear_index = address(ix0);
int i, j, k;
std::vector<int> ix = ix0;
const cvm::real * g;
if (nd == 2) {
// gradients at grid points surrounding the current scalar grid point
@ -783,7 +782,7 @@ void integrate_potential::atimes(const std::vector<cvm::real> &A, std::vector<cv
/// Inversion of preconditioner matrix (e.g. diagonal of the Laplacian)
void integrate_potential::asolve(const std::vector<cvm::real> &b, std::vector<cvm::real> &x)
{
for (size_t i=0; i<nt; i++) {
for (size_t i=0; i<int(nt); i++) {
x[i] = b[i] * inv_lap_diag[i]; // Jacobi preconditioner - little benefit in tests so far
}
return;
@ -803,7 +802,7 @@ void integrate_potential::nr_linbcg_sym(const std::vector<cvm::real> &b, std::ve
iter=0;
atimes(x,r);
for (j=0;j<nt;j++) {
for (j=0;j<int(nt);j++) {
r[j]=b[j]-r[j];
}
bnrm=l2norm(b);
@ -814,26 +813,26 @@ void integrate_potential::nr_linbcg_sym(const std::vector<cvm::real> &b, std::ve
bkden = 1.0;
while (iter < itmax) {
++iter;
for (bknum=0.0,j=0;j<nt;j++) {
for (bknum=0.0,j=0;j<int(nt);j++) {
bknum += r[j]*r[j]; // precon: z[j]*r[j]
}
if (iter == 1) {
for (j=0;j<nt;j++) {
for (j=0;j<int(nt);j++) {
p[j] = r[j]; // precon: p[j] = z[j]
}
} else {
bk=bknum/bkden;
for (j=0;j<nt;j++) {
for (j=0;j<int(nt);j++) {
p[j] = bk*p[j] + r[j]; // precon: bk*p[j] + z[j]
}
}
bkden = bknum;
atimes(p,z);
for (akden=0.0,j=0;j<nt;j++) {
for (akden=0.0,j=0;j<int(nt);j++) {
akden += z[j]*p[j];
}
ak = bknum/akden;
for (j=0;j<nt;j++) {
for (j=0;j<int(nt);j++) {
x[j] += ak*p[j];
r[j] -= ak*z[j];
}

157
lib/colvars/colvarmodule.cpp
View File

@ -30,10 +30,10 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
depth_s = 0;
cv_traj_os = NULL;
// pointer to the proxy object
if (proxy == NULL) {
proxy = proxy_in;
parse = new colvarparse();
proxy = proxy_in; // Pointer to the proxy object
parse = new colvarparse(); // Parsing object for global options
version_int = proxy->get_version_from_string(COLVARS_VERSION);
} else {
// TODO relax this error to handle multiple molecules in VMD
// once the module is not static anymore
@ -58,22 +58,24 @@ colvarmodule::colvarmodule(colvarproxy *proxy_in)
// "it_restart" will be set by the input state file, if any;
// "it" should be updated by the proxy
colvarmodule::it = colvarmodule::it_restart = 0;
colvarmodule::it_restart_from_state_file = true;
colvarmodule::use_scripted_forces = false;
use_scripted_forces = false;
scripting_after_biases = false;
colvarmodule::b_analysis = false;
b_analysis = false;
colvarmodule::debug_gradients_step_size = 1.0e-07;
colvarmodule::rotation::monitor_crossings = false;
colvarmodule::rotation::crossing_threshold = 1.0e-02;
colvarmodule::cv_traj_freq = 100;
colvarmodule::restart_out_freq = proxy->restart_frequency();
cv_traj_freq = 100;
restart_out_freq = proxy->restart_frequency();
// by default overwrite the existing trajectory file
colvarmodule::cv_traj_append = false;
cv_traj_append = false;
cv_traj_write_labels = true;
}
@ -189,26 +191,27 @@ int colvarmodule::parse_config(std::string &conf)
{
extra_conf.clear();
// parse global options
// Parse global options
if (catch_input_errors(parse_global_params(conf))) {
return get_error();
}
// parse the options for collective variables
// Parse the options for collective variables
if (catch_input_errors(parse_colvars(conf))) {
return get_error();
}
// parse the options for biases
// Parse the options for biases
if (catch_input_errors(parse_biases(conf))) {
return get_error();
}
// done parsing known keywords, check that all keywords found were valid ones
// Done parsing known keywords, check that all keywords found were valid ones
if (catch_input_errors(parse->check_keywords(conf, "colvarmodule"))) {
return get_error();
}
// Parse auto-generated configuration (e.g. for back-compatibility)
if (extra_conf.size()) {
catch_input_errors(parse_global_params(extra_conf));
catch_input_errors(parse_colvars(extra_conf));
@ -222,13 +225,11 @@ int colvarmodule::parse_config(std::string &conf)
cvm::log("Collective variables module (re)initialized.\n");
cvm::log(cvm::line_marker);
// update any necessary proxy data
// Update any necessary proxy data
proxy->setup();
if (cv_traj_os != NULL) {
// configuration might have changed, better redo the labels
write_traj_label(*cv_traj_os);
}
// configuration might have changed, better redo the labels
cv_traj_write_labels = true;
return get_error();
}
@ -279,15 +280,18 @@ int colvarmodule::parse_global_params(std::string const &conf)
parse->get_keyval(conf, "colvarsTrajAppend",
cv_traj_append, cv_traj_append, colvarparse::parse_silent);
parse->get_keyval(conf, "scriptedColvarForces", use_scripted_forces, false);
parse->get_keyval(conf, "scriptedColvarForces",
use_scripted_forces, use_scripted_forces);
parse->get_keyval(conf, "scriptingAfterBiases", scripting_after_biases, true);
parse->get_keyval(conf, "scriptingAfterBiases",
scripting_after_biases, scripting_after_biases);
if (use_scripted_forces && !proxy->force_script_defined) {
cvm::error("User script for scripted colvar forces not found.", INPUT_ERROR);
return cvm::error("User script for scripted colvar forces not found.",
INPUT_ERROR);
}
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
return cvm::get_error();
}
@ -432,13 +436,13 @@ int colvarmodule::parse_biases(std::string const &conf)
}
int colvarmodule::num_variables() const
size_t colvarmodule::num_variables() const
{
return colvars.size();
}
int colvarmodule::num_variables_feature(int feature_id) const
size_t colvarmodule::num_variables_feature(int feature_id) const
{
size_t n = 0;
for (std::vector<colvar *>::const_iterator cvi = colvars.begin();
@ -452,13 +456,13 @@ int colvarmodule::num_variables_feature(int feature_id) const
}
int colvarmodule::num_biases() const
size_t colvarmodule::num_biases() const
{
return biases.size();
}
int colvarmodule::num_biases_feature(int feature_id) const
size_t colvarmodule::num_biases_feature(int feature_id) const
{
size_t n = 0;
for (std::vector<colvarbias *>::const_iterator bi = biases.begin();
@ -472,7 +476,7 @@ int colvarmodule::num_biases_feature(int feature_id) const
}
int colvarmodule::num_biases_type(std::string const &type) const
size_t colvarmodule::num_biases_type(std::string const &type) const
{
size_t n = 0;
for (std::vector<colvarbias *>::const_iterator bi = biases.begin();
@ -971,13 +975,20 @@ int colvarmodule::write_restart_file(std::string const &out_name)
int colvarmodule::write_traj_files()
{
if (cv_traj_os == NULL) {
open_traj_file(cv_traj_name);
if (open_traj_file(cv_traj_name) != COLVARS_OK) {
return cvm::get_error();
} else {
cv_traj_write_labels = true;
}
}
// write labels in the traj file every 1000 lines and at first timestep
if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 || cvm::step_relative() == 0) {
if ((cvm::step_absolute() % (cv_traj_freq * 1000)) == 0 ||
cvm::step_relative() == 0 ||
cv_traj_write_labels) {
write_traj_label(*cv_traj_os);
}
cv_traj_write_labels = false;
if ((cvm::step_absolute() % cv_traj_freq) == 0) {
write_traj(*cv_traj_os);
@ -1064,7 +1075,7 @@ int colvarmodule::reset()
{
parse->init();
cvm::log("Resetting the Collective Variables Module.\n");
cvm::log("Resetting the Collective Variables module.\n");
// Iterate backwards because we are deleting the elements as we go
for (std::vector<colvarbias *>::reverse_iterator bi = biases.rbegin();
@ -1073,6 +1084,7 @@ int colvarmodule::reset()
delete *bi; // the bias destructor updates the biases array
}
biases.clear();
biases_active_.clear();
// Iterate backwards because we are deleting the elements as we go
for (std::vector<colvar *>::reverse_iterator cvi = colvars.rbegin();
@ -1088,7 +1100,7 @@ int colvarmodule::reset()
proxy->reset();
if (cv_traj_os != NULL) {
// Do not close file here, as we might not be done with it yet.
// Do not close traj file here, as we might not be done with it yet.
proxy->flush_output_stream(cv_traj_os);
}
@ -1188,12 +1200,10 @@ std::istream & colvarmodule::read_restart(std::istream &is)
// read global restart information
std::string restart_conf;
if (is >> colvarparse::read_block("configuration", restart_conf)) {
if (it_restart_from_state_file) {
parse->get_keyval(restart_conf, "step",
it_restart, (size_t) 0,
colvarparse::parse_silent);
parse->get_keyval(restart_conf, "step",
it_restart, (size_t) 0,
colvarparse::parse_silent);
it = it_restart;
}
std::string restart_version;
parse->get_keyval(restart_conf, "version",
restart_version, std::string(""),
@ -1688,40 +1698,59 @@ int cvm::read_index_file(char const *filename)
return (cvm::get_error() ? COLVARS_ERROR : COLVARS_OK);
}
int cvm::load_atoms(char const *file_name,
cvm::atom_group &atoms,
std::string const &pdb_field,
double const pdb_field_value)
double pdb_field_value)
{
return proxy->load_atoms(file_name, atoms, pdb_field, pdb_field_value);
}
int cvm::load_coords(char const *file_name,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::string const &pdb_field,
double const pdb_field_value)
{
// Differentiate between PDB and XYZ files
// for XYZ files, use CVM internal parser
// otherwise call proxy function for PDB
std::string const ext(strlen(file_name) > 4 ? (file_name + (strlen(file_name) - 4)) : file_name);
int cvm::load_coords(char const *file_name,
std::vector<cvm::rvector> *pos,
cvm::atom_group *atoms,
std::string const &pdb_field,
double pdb_field_value)
{
int error_code = COLVARS_OK;
std::string const ext(strlen(file_name) > 4 ?
(file_name + (strlen(file_name) - 4)) :
file_name);
atoms->create_sorted_ids();
std::vector<cvm::rvector> sorted_pos(atoms->size(), cvm::rvector(0.0));
// Differentiate between PDB and XYZ files
if (colvarparse::to_lower_cppstr(ext) == std::string(".xyz")) {
if ( pdb_field.size() > 0 ) {
cvm::error("Error: PDB column may not be specified for XYZ coordinate file.\n", INPUT_ERROR);
return COLVARS_ERROR;
if (pdb_field.size() > 0) {
return cvm::error("Error: PDB column may not be specified "
"for XYZ coordinate files.\n", INPUT_ERROR);
}
return cvm::load_coords_xyz(file_name, pos, indices);
// For XYZ files, use internal parser
error_code |= cvm::load_coords_xyz(file_name, &sorted_pos, atoms);
} else {
return proxy->load_coords(file_name, pos, indices, pdb_field, pdb_field_value);
// Otherwise, call proxy function for PDB
error_code |= proxy->load_coords(file_name,
sorted_pos, atoms->sorted_ids(),
pdb_field, pdb_field_value);
}
std::vector<int> const &map = atoms->sorted_ids_map();
for (size_t i = 0; i < atoms->size(); i++) {
(*pos)[map[i]] = sorted_pos[i];
}
return error_code;
}
int cvm::load_coords_xyz(char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices)
std::vector<rvector> *pos,
cvm::atom_group *atoms)
{
std::ifstream xyz_is(filename);
unsigned int natoms;
@ -1736,12 +1765,12 @@ int cvm::load_coords_xyz(char const *filename,
cvm::getline(xyz_is, line);
cvm::getline(xyz_is, line);
xyz_is.width(255);
std::vector<atom_pos>::iterator pos_i = pos.begin();
std::vector<atom_pos>::iterator pos_i = pos->begin();
if (pos.size() != natoms) { // Use specified indices
if (pos->size() != natoms) { // Use specified indices
int next = 0; // indices are zero-based
std::vector<int>::const_iterator index = indices.begin();
for ( ; pos_i != pos.end() ; pos_i++, index++) {
std::vector<int>::const_iterator index = atoms->sorted_ids().begin();
for ( ; pos_i != pos->end() ; pos_i++, index++) {
while ( next < *index ) {
cvm::getline(xyz_is, line);
@ -1751,7 +1780,7 @@ int cvm::load_coords_xyz(char const *filename,
xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
}
} else { // Use all positions
for ( ; pos_i != pos.end() ; pos_i++) {
for ( ; pos_i != pos->end() ; pos_i++) {
xyz_is >> symbol;
xyz_is >> (*pos_i)[0] >> (*pos_i)[1] >> (*pos_i)[2];
}
@ -1792,17 +1821,13 @@ void cvm::request_total_force()
proxy->request_total_force(true);
}
cvm::rvector cvm::position_distance(atom_pos const &pos1,
atom_pos const &pos2)
cvm::rvector cvm::position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
return proxy->position_distance(pos1, pos2);
}
cvm::real cvm::position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
return proxy->position_dist2(pos1, pos2);
}
cvm::real cvm::rand_gaussian(void)
{

71
lib/colvars/colvarmodule.h
View File

@ -73,8 +73,17 @@ private:
/// Impossible to initialize the main object without arguments
colvarmodule();
/// Integer representing the version string (allows comparisons)
int version_int;
public:
/// Get the version number (higher = more recent)
int version_number() const
{
return version_int;
}
friend class colvarproxy;
// TODO colvarscript should be unaware of colvarmodule's internals
friend class colvarscript;
@ -158,10 +167,6 @@ public:
return it;
}
/// If true, get it_restart from the state file; if set to false,
/// the MD program is providing it
bool it_restart_from_state_file;
/// \brief Finite difference step size (if there is no dynamics, or
/// if gradients need to be tested independently from the size of
/// dt)
@ -306,19 +311,19 @@ private:
public:
/// Return how many variables are defined
int num_variables() const;
size_t num_variables() const;
/// Return how many variables have this feature enabled
int num_variables_feature(int feature_id) const;
size_t num_variables_feature(int feature_id) const;
/// Return how many biases are defined
int num_biases() const;
size_t num_biases() const;
/// Return how many biases have this feature enabled
int num_biases_feature(int feature_id) const;
size_t num_biases_feature(int feature_id) const;
/// Return how many biases of this type are defined
int num_biases_type(std::string const &type) const;
size_t num_biases_type(std::string const &type) const;
/// Return the names of time-dependent biases with forces enabled (ABF,
/// metadynamics, etc)
@ -479,9 +484,6 @@ public:
/// Print a message to the main log and set global error code
static int error(std::string const &message, int code = COLVARS_ERROR);
/// 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();
@ -495,15 +497,6 @@ public:
static rvector position_distance(atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Get the square distance between two positions (with
/// periodic boundary conditions handled transparently)
///
/// Note: in the case of periodic boundary conditions, this provides
/// an analytical square distance (while taking the square of
/// position_distance() would produce leads to a cusp)
static real position_dist2(atom_pos const &pos1,
atom_pos const &pos2);
/// \brief Names of groups from a Gromacs .ndx file to be read at startup
std::list<std::string> index_group_names;
@ -513,29 +506,36 @@ public:
/// \brief Read a Gromacs .ndx file
int read_index_file(char const *filename);
/// \brief Create atoms from a file \param filename name of the file
/// (usually a PDB) \param atoms array of the atoms to be allocated
/// \param pdb_field (optiona) if "filename" is a PDB file, use this
/// field to determine which are the atoms to be set
/// \brief Select atom IDs from a file (usually PDB) \param filename name of
/// the file \param atoms array into which atoms read from "filename" will be
/// appended \param pdb_field (optional) if the file is a PDB and this
/// string is non-empty, select atoms for which this field is non-zero
/// \param pdb_field_value (optional) if non-zero, select only atoms whose
/// pdb_field equals this
static int load_atoms(char const *filename,
atom_group &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
double pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from a file
/// (PDB or XYZ)
/// \brief Load coordinates for a group of atoms from a file (PDB or XYZ);
/// if "pos" is already allocated, the number of its elements must match the
/// number of entries in "filename" \param filename name of the file \param
/// pos array of coordinates \param atoms group containing the atoms (used
/// to obtain internal IDs) \param pdb_field (optional) if the file is a PDB
/// and this string is non-empty, select atoms for which this field is
/// non-zero \param pdb_field_value (optional) if non-zero, select only
/// atoms whose pdb_field equals this
static int load_coords(char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices,
std::vector<rvector> *pos,
atom_group *atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
double pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from an
/// XYZ file
static int load_coords_xyz(char const *filename,
std::vector<atom_pos> &pos,
const std::vector<int> &indices);
std::vector<rvector> *pos,
atom_group *atoms);
/// Frequency for collective variables trajectory output
static size_t cv_traj_freq;
@ -568,6 +568,9 @@ protected:
/// Appending to the existing trajectory file?
bool cv_traj_append;
/// Write labels at the next iteration
bool cv_traj_write_labels;
private:
/// Counter for the current depth in the object hierarchy (useg e.g. in output)

91
lib/colvars/colvarproxy.cpp
View File

@ -26,7 +26,10 @@
colvarproxy_system::colvarproxy_system() {}
colvarproxy_system::colvarproxy_system()
{
reset_pbc_lattice();
}
colvarproxy_system::~colvarproxy_system() {}
@ -55,10 +58,73 @@ bool colvarproxy_system::total_forces_same_step() const
}
cvm::real colvarproxy_system::position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
inline int round_to_integer(cvm::real x)
{
return (position_distance(pos1, pos2)).norm2();
return std::floor(x+0.5);
}
void colvarproxy_system::update_pbc_lattice()
{
// Periodicity is assumed in all directions
if (boundaries_type == boundaries_unsupported ||
boundaries_type == boundaries_non_periodic) {
cvm::error("Error: setting PBC lattice with unsupported boundaries.\n",
BUG_ERROR);
return;
}
{
cvm::rvector const v = cvm::rvector::outer(unit_cell_y, unit_cell_z);
reciprocal_cell_x = v/(v*unit_cell_x);
}
{
cvm::rvector const v = cvm::rvector::outer(unit_cell_z, unit_cell_x);
reciprocal_cell_y = v/(v*unit_cell_y);
}
{
cvm::rvector const v = cvm::rvector::outer(unit_cell_x, unit_cell_y);
reciprocal_cell_z = v/(v*unit_cell_z);
}
}
void colvarproxy_system::reset_pbc_lattice()
{
unit_cell_x.reset();
unit_cell_y.reset();
unit_cell_z.reset();
reciprocal_cell_x.reset();
reciprocal_cell_y.reset();
reciprocal_cell_z.reset();
}
cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
const
{
if (boundaries_type == boundaries_unsupported) {
cvm::error("Error: unsupported boundary conditions.\n", INPUT_ERROR);
}
cvm::rvector diff = (pos2 - pos1);
if (boundaries_type == boundaries_non_periodic) return diff;
cvm::real const x_shift = round_to_integer(reciprocal_cell_x*diff);
cvm::real const y_shift = round_to_integer(reciprocal_cell_y*diff);
cvm::real const z_shift = round_to_integer(reciprocal_cell_z*diff);
diff.x -= x_shift*unit_cell_x.x + y_shift*unit_cell_y.x +
z_shift*unit_cell_z.x;
diff.y -= x_shift*unit_cell_x.y + y_shift*unit_cell_y.y +
z_shift*unit_cell_z.y;
diff.z -= x_shift*unit_cell_x.z + y_shift*unit_cell_y.z +
z_shift*unit_cell_z.z;
return diff;
}
@ -132,7 +198,7 @@ void colvarproxy_atoms::clear_atom(int index)
int colvarproxy_atoms::load_atoms(char const *filename,
cvm::atom_group &atoms,
std::string const &pdb_field,
double const)
double)
{
return cvm::error("Error: loading atom identifiers from a file "
"is currently not implemented.\n",
@ -142,9 +208,9 @@ int colvarproxy_atoms::load_atoms(char const *filename,
int colvarproxy_atoms::load_coords(char const *filename,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::vector<int> const &sorted_ids,
std::string const &pdb_field,
double const)
double)
{
return cvm::error("Error: loading atomic coordinates from a file "
"is currently not implemented.\n",
@ -661,6 +727,7 @@ int colvarproxy_io::close_output_stream(std::string const &output_name)
for ( ; osi != output_files.end(); osi++, osni++) {
if (*osni == output_name) {
((std::ofstream *) (*osi))->close();
delete *osi;
output_files.erase(osi);
output_stream_names.erase(osni);
return COLVARS_OK;
@ -729,3 +796,13 @@ size_t colvarproxy::restart_frequency()
return 0;
}
int colvarproxy::get_version_from_string(char const *version_string)
{
std::string const v(version_string);
std::istringstream is(v.substr(0, 4) + v.substr(5, 2) + v.substr(8, 2));
int newint;
is >> newint;
return newint;
}

80
lib/colvars/colvarproxy.h
View File

@ -14,6 +14,7 @@
#include <list>
#include "colvarmodule.h"
#include "colvartypes.h"
#include "colvarvalue.h"
@ -29,7 +30,7 @@
/// To interface to a new MD engine, the simplest solution is to derive a new
/// class from \link colvarproxy \endlink. Currently implemented are: \link
/// colvarproxy_lammps, \endlink, \link colvarproxy_namd, \endlink, \link
/// colvarproxy_vmd, \endlink.
/// colvarproxy_vmd \endlink.
// forward declarations
@ -68,14 +69,16 @@ public:
/// \brief Get the PBC-aware distance vector between two positions
virtual cvm::rvector position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2) = 0;
cvm::atom_pos const &pos2) const;
/// \brief Get the PBC-aware square distance between two positions;
/// may need to be reimplemented independently from position_distance() for optimization purposes
virtual cvm::real position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2);
/// Recompute PBC reciprocal lattice (assumes XYZ periodicity)
void update_pbc_lattice();
/// Tell the proxy whether total forces are needed (may not always be available)
/// Set the lattice vectors to zero
void reset_pbc_lattice();
/// \brief Tell the proxy whether total forces are needed (they may not
/// always be available)
virtual void request_total_force(bool yesno);
/// Are total forces being used?
@ -83,6 +86,29 @@ public:
/// Are total forces from the current step available?
virtual bool total_forces_same_step() const;
protected:
/// \brief Type of boundary conditions
///
/// Orthogonal and triclinic cells are made available to objects.
/// For any other conditions (mixed periodicity, triclinic cells in LAMMPS)
/// minimum-image distances are computed by the host engine regardless.
enum Boundaries_type {
boundaries_non_periodic,
boundaries_pbc_ortho,
boundaries_pbc_triclinic,
boundaries_unsupported
};
/// Type of boundary conditions
Boundaries_type boundaries_type;
/// Bravais lattice vectors
cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z;
/// Reciprocal lattice vectors
cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z;
};
@ -121,24 +147,30 @@ public:
/// (costly) set the corresponding atoms_ncopies to zero
virtual void clear_atom(int index);
/// \brief Read atom identifiers from a file \param filename name of
/// the file (usually a PDB) \param atoms array to which atoms read
/// from "filename" will be appended \param pdb_field (optiona) if
/// "filename" is a PDB file, use this field to determine which are
/// the atoms to be set
/// \brief Select atom IDs from a file (usually PDB) \param filename name of
/// the file \param atoms array to which atoms read from "filename" will be
/// appended \param pdb_field (optional) if the file is a PDB and this
/// string is non-empty, select atoms for which this field is non-zero
/// \param pdb_field_value (optional) if non-zero, select only atoms whose
/// pdb_field equals this
virtual int load_atoms(char const *filename,
cvm::atom_group &atoms,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
double pdb_field_value = 0.0);
/// \brief Load the coordinates for a group of atoms from a file
/// (usually a PDB); if "pos" is already allocated, the number of its
/// elements must match the number of atoms in "filename"
/// \brief Load a set of coordinates from a file (usually PDB); if "pos" is
/// already allocated, the number of its elements must match the number of
/// entries in "filename" \param filename name of the file \param pos array
/// of coordinates \param sorted_ids array of sorted internal IDs, used to
/// loop through the file only once \param pdb_field (optional) if the file
/// is a PDB and this string is non-empty, select atoms for which this field
/// is non-zero \param pdb_field_value (optional) if non-zero, select only
/// atoms whose pdb_field equals this
virtual int load_coords(char const *filename,
std::vector<cvm::atom_pos> &pos,
const std::vector<int> &indices,
std::vector<int> const &sorted_ids,
std::string const &pdb_field,
double const pdb_field_value = 0.0);
double pdb_field_value = 0.0);
/// Clear atomic data
int reset();
@ -636,6 +668,15 @@ public:
return b_simulation_running;
}
/// Convert a version string "YYYY-MM-DD" into an integer
int get_version_from_string(char const *version_string);
/// Get the version number (higher = more recent)
int version_number() const
{
return version_int;
}
protected:
/// Whether a simulation is running (warn against irrecovarable errors)
@ -644,6 +685,9 @@ protected:
/// Whether the entire module should be deallocated by the host engine
bool b_delete_requested;
/// Integer representing the version string (allows comparisons)
int version_int;
};

2
lib/colvars/colvars_version.h
View File

@ -1,5 +1,5 @@
#ifndef COLVARS_VERSION
#define COLVARS_VERSION "2018-01-17"
#define COLVARS_VERSION "2018-04-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

308
lib/colvars/colvartypes.h
View File

@ -130,7 +130,7 @@ public:
}
}
inline void operator *= (cvm::real const &a)
inline void operator *= (cvm::real a)
{
size_t i;
for (i = 0; i < this->size(); i++) {
@ -138,7 +138,7 @@ public:
}
}
inline void operator /= (cvm::real const &a)
inline void operator /= (cvm::real a)
{
size_t i;
for (i = 0; i < this->size(); i++) {
@ -146,7 +146,8 @@ public:
}
}
inline friend vector1d<T> operator + (vector1d<T> const &v1, vector1d<T> const &v2)
inline friend vector1d<T> operator + (vector1d<T> const &v1,
vector1d<T> const &v2)
{
check_sizes(v1.size(), v2.size());
vector1d<T> result(v1.size());
@ -157,7 +158,8 @@ public:
return result;
}
inline friend vector1d<T> operator - (vector1d<T> const &v1, vector1d<T> const &v2)
inline friend vector1d<T> operator - (vector1d<T> const &v1,
vector1d<T> const &v2)
{
check_sizes(v1.size(), v2.size());
vector1d<T> result(v1.size());
@ -168,7 +170,7 @@ public:
return result;
}
inline friend vector1d<T> operator * (vector1d<T> const &v, cvm::real const &a)
inline friend vector1d<T> operator * (vector1d<T> const &v, cvm::real a)
{
vector1d<T> result(v.size());
size_t i;
@ -178,12 +180,12 @@ public:
return result;
}
inline friend vector1d<T> operator * (cvm::real const &a, vector1d<T> const &v)
inline friend vector1d<T> operator * (cvm::real a, vector1d<T> const &v)
{
return v * a;
}
inline friend vector1d<T> operator / (vector1d<T> const &v, cvm::real const &a)
inline friend vector1d<T> operator / (vector1d<T> const &v, cvm::real a)
{
vector1d<T> result(v.size());
size_t i;
@ -246,7 +248,8 @@ public:
}
/// Assign a vector to a slice of this vector
inline void sliceassign(size_t const i1, size_t const i2, vector1d<T> const &v)
inline void sliceassign(size_t const i1, size_t const i2,
vector1d<T> const &v)
{
if ((i2 < i1) || (i2 >= this->size())) {
cvm::error("Error: trying to slice a vector using incorrect boundaries.\n");
@ -259,12 +262,13 @@ public:
/// Formatted output
inline size_t output_width(size_t const &real_width) const
inline size_t output_width(size_t real_width) const
{
return real_width*(this->size()) + 3*(this->size()-1) + 4;
}
inline friend std::istream & operator >> (std::istream &is, cvm::vector1d<T> &v)
inline friend std::istream & operator >> (std::istream &is,
cvm::vector1d<T> &v)
{
if (v.size() == 0) return is;
size_t const start_pos = is.tellg();
@ -288,7 +292,8 @@ public:
return is;
}
inline friend std::ostream & operator << (std::ostream &os, cvm::vector1d<T> const &v)
inline friend std::ostream & operator << (std::ostream &os,
cvm::vector1d<T> const &v)
{
std::streamsize const w = os.width();
std::streamsize const p = os.precision();
@ -377,6 +382,15 @@ protected:
{
return vector1d<T>(length, data);
}
inline int set(cvm::vector1d<T> const &v) const
{
if (v.size() != length) {
return cvm::error("Error: setting a matrix row from a vector of "
"incompatible size.\n", BUG_ERROR);
}
for (size_t i = 0; i < length; i++) data[i] = v[i];
return COLVARS_OK;
}
};
std::vector<T> data;
@ -515,9 +529,12 @@ public:
{
if ((m1.outer_length != m2.outer_length) ||
(m1.inner_length != m2.inner_length)) {
cvm::error("Error: trying to perform an operation between matrices of different sizes, "+
cvm::to_str(m1.outer_length)+"x"+cvm::to_str(m1.inner_length)+" and "+
cvm::to_str(m2.outer_length)+"x"+cvm::to_str(m2.inner_length)+".\n");
cvm::error("Error: trying to perform an operation between "
"matrices of different sizes, "+
cvm::to_str(m1.outer_length)+"x"+
cvm::to_str(m1.inner_length)+" and "+
cvm::to_str(m2.outer_length)+"x"+
cvm::to_str(m2.inner_length)+".\n");
}
}
@ -539,7 +556,7 @@ public:
}
}
inline void operator *= (cvm::real const &a)
inline void operator *= (cvm::real a)
{
size_t i;
for (i = 0; i < data.size(); i++) {
@ -547,7 +564,7 @@ public:
}
}
inline void operator /= (cvm::real const &a)
inline void operator /= (cvm::real a)
{
size_t i;
for (i = 0; i < data.size(); i++) {
@ -555,7 +572,8 @@ public:
}
}
inline friend matrix2d<T> operator + (matrix2d<T> const &m1, matrix2d<T> const &m2)
inline friend matrix2d<T> operator + (matrix2d<T> const &m1,
matrix2d<T> const &m2)
{
check_sizes(m1, m2);
matrix2d<T> result(m1.outer_length, m1.inner_length);
@ -566,7 +584,8 @@ public:
return result;
}
inline friend matrix2d<T> operator - (matrix2d<T> const &m1, matrix2d<T> const &m2)
inline friend matrix2d<T> operator - (matrix2d<T> const &m1,
matrix2d<T> const &m2)
{
check_sizes(m1, m2);
matrix2d<T> result(m1.outer_length, m1.inner_length);
@ -577,7 +596,7 @@ public:
return result;
}
inline friend matrix2d<T> operator * (matrix2d<T> const &m, cvm::real const &a)
inline friend matrix2d<T> operator * (matrix2d<T> const &m, cvm::real a)
{
matrix2d<T> result(m.outer_length, m.inner_length);
size_t i;
@ -587,12 +606,12 @@ public:
return result;
}
inline friend matrix2d<T> operator * (cvm::real const &a, matrix2d<T> const &m)
inline friend matrix2d<T> operator * (cvm::real a, matrix2d<T> const &m)
{
return m * a;
}
inline friend matrix2d<T> operator / (matrix2d<T> const &m, cvm::real const &a)
inline friend matrix2d<T> operator / (matrix2d<T> const &m, cvm::real a)
{
matrix2d<T> result(m.outer_length, m.inner_length);
size_t i;
@ -602,34 +621,17 @@ public:
return result;
}
/// Matrix multiplication
// inline friend matrix2d<T> const & operator * (matrix2d<T> const &m1, matrix2d<T> const &m2)
// {
// matrix2d<T> result(m1.outer_length, m2.inner_length);
// if (m1.inner_length != m2.outer_length) {
// cvm::error("Error: trying to multiply two matrices of incompatible sizes, "+
// cvm::to_str(m1.outer_length)+"x"+cvm::to_str(m1.inner_length)+" and "+
// cvm::to_str(m2.outer_length)+"x"+cvm::to_str(m2.inner_length)+".\n");
// } else {
// size_t i, j, k;
// for (i = 0; i < m1.outer_length; i++) {
// for (j = 0; j < m2.inner_length; j++) {
// for (k = 0; k < m1.inner_length; k++) {
// result[i][j] += m1[i][k] * m2[k][j];
// }
// }
// }
// }
// return result;
// }
/// vector-matrix multiplication
inline friend vector1d<T> operator * (vector1d<T> const &v, matrix2d<T> const &m)
inline friend vector1d<T> operator * (vector1d<T> const &v,
matrix2d<T> const &m)
{
vector1d<T> result(m.inner_length);
if (m.outer_length != v.size()) {
cvm::error("Error: trying to multiply a vector and a matrix of incompatible sizes, "+
cvm::to_str(v.size()) + " and " + cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length) + ".\n");
cvm::error("Error: trying to multiply a vector and a matrix "
"of incompatible sizes, "+
cvm::to_str(v.size()) + " and " +
cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length) +
".\n");
} else {
size_t i, k;
for (i = 0; i < m.inner_length; i++) {
@ -641,25 +643,6 @@ public:
return result;
}
// /// matrix-vector multiplication (unused for now)
// inline friend vector1d<T> const & operator * (matrix2d<T> const &m, vector1d<T> const &v)
// {
// vector1d<T> result(m.outer_length);
// if (m.inner_length != v.size()) {
// cvm::error("Error: trying to multiply a matrix and a vector of incompatible sizes, "+
// cvm::to_str(m.outer_length)+"x"+cvm::to_str(m.inner_length)
// + " and " + cvm::to_str(v.length) + ".\n");
// } else {
// size_t i, k;
// for (i = 0; i < m.outer_length; i++) {
// for (k = 0; k < m.inner_length; k++) {
// result[i] += m[i][k] * v[k];
// }
// }
// }
// return result;
// }
/// Formatted output
friend std::ostream & operator << (std::ostream &os,
matrix2d<T> const &m)
@ -725,49 +708,52 @@ public:
cvm::real x, y, z;
inline rvector()
: x(0.0), y(0.0), z(0.0)
{}
{
reset();
}
inline rvector(cvm::real const &x_i,
cvm::real const &y_i,
cvm::real const &z_i)
: x(x_i), y(y_i), z(z_i)
{}
/// \brief Set all components to zero
inline void reset()
{
set(0.0);
}
inline rvector(cvm::real x_i, cvm::real y_i, cvm::real z_i)
{
set(x_i, y_i, z_i);
}
inline rvector(cvm::vector1d<cvm::real> const &v)
: x(v[0]), y(v[1]), z(v[2])
{}
{
set(v[0], v[1], v[2]);
}
inline rvector(cvm::real t)
: x(t), y(t), z(t)
{}
{
set(t);
}
/// \brief Set all components to a scalar value
inline void set(cvm::real const &value) {
/// \brief Set all components to a scalar
inline void set(cvm::real value)
{
x = y = z = value;
}
/// \brief Assign all components
inline void set(cvm::real const &x_i,
cvm::real const &y_i,
cvm::real const &z_i) {
inline void set(cvm::real x_i, cvm::real y_i, cvm::real z_i)
{
x = x_i;
y = y_i;
z = z_i;
}
/// \brief Set all components to zero
inline void reset() {
x = y = z = 0.0;
}
/// \brief Access cartesian components by index
inline cvm::real & operator [] (int const &i) {
inline cvm::real & operator [] (int i) {
return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : x;
}
/// \brief Access cartesian components by index
inline cvm::real const & operator [] (int const &i) const {
inline cvm::real operator [] (int i) const {
return (i == 0) ? x : (i == 1) ? y : (i == 2) ? z : x;
}
@ -780,14 +766,6 @@ public:
return result;
}
inline cvm::rvector & operator = (cvm::real const &v)
{
x = v;
y = v;
z = v;
return *this;
}
inline void operator += (cvm::rvector const &v)
{
x += v.x;
@ -802,7 +780,7 @@ public:
z -= v.z;
}
inline void operator *= (cvm::real const &v)
inline void operator *= (cvm::real v)
{
x *= v;
y *= v;
@ -832,13 +810,14 @@ public:
return (n > 0. ? cvm::rvector(x, y, z)/n : cvm::rvector(1., 0., 0.));
}
static inline size_t output_width(size_t const &real_width)
static inline size_t output_width(size_t real_width)
{
return 3*real_width + 10;
}
static inline cvm::rvector outer(cvm::rvector const &v1, cvm::rvector const &v2)
static inline cvm::rvector outer(cvm::rvector const &v1,
cvm::rvector const &v2)
{
return cvm::rvector( v1.y*v2.z - v2.y*v1.z,
-v1.x*v2.z + v2.x*v1.z,
@ -850,41 +829,35 @@ public:
return cvm::rvector(-v.x, -v.y, -v.z);
}
friend inline int operator == (cvm::rvector const &v1, cvm::rvector const &v2)
{
return (v1.x == v2.x) && (v1.y == v2.y) && (v1.z == v2.z);
}
friend inline int operator != (cvm::rvector const &v1, cvm::rvector const &v2)
{
return (v1.x != v2.x) || (v1.y != v2.y) || (v1.z != v2.z);
}
friend inline cvm::rvector operator + (cvm::rvector const &v1, cvm::rvector const &v2)
friend inline cvm::rvector operator + (cvm::rvector const &v1,
cvm::rvector const &v2)
{
return cvm::rvector(v1.x + v2.x, v1.y + v2.y, v1.z + v2.z);
}
friend inline cvm::rvector operator - (cvm::rvector const &v1, cvm::rvector const &v2)
friend inline cvm::rvector operator - (cvm::rvector const &v1,
cvm::rvector const &v2)
{
return cvm::rvector(v1.x - v2.x, v1.y - v2.y, v1.z - v2.z);
}
friend inline cvm::real operator * (cvm::rvector const &v1, cvm::rvector const &v2)
/// Inner (dot) product
friend inline cvm::real operator * (cvm::rvector const &v1,
cvm::rvector const &v2)
{
return v1.x * v2.x + v1.y * v2.y + v1.z * v2.z;
}
friend inline cvm::rvector operator * (cvm::real const &a, cvm::rvector const &v)
friend inline cvm::rvector operator * (cvm::real a, cvm::rvector const &v)
{
return cvm::rvector(a*v.x, a*v.y, a*v.z);
}
friend inline cvm::rvector operator * (cvm::rvector const &v, cvm::real const &a)
friend inline cvm::rvector operator * (cvm::rvector const &v, cvm::real a)
{
return cvm::rvector(a*v.x, a*v.y, a*v.z);
}
friend inline cvm::rvector operator / (cvm::rvector const &v, cvm::real const &a)
friend inline cvm::rvector operator / (cvm::rvector const &v, cvm::real a)
{
return cvm::rvector(v.x/a, v.y/a, v.z/a);
}
@ -946,15 +919,9 @@ public:
{}
/// Constructor component by component
inline rmatrix(cvm::real const &xxi,
cvm::real const &xyi,
cvm::real const &xzi,
cvm::real const &yxi,
cvm::real const &yyi,
cvm::real const &yzi,
cvm::real const &zxi,
cvm::real const &zyi,
cvm::real const &zzi)
inline rmatrix(cvm::real xxi, cvm::real xyi, cvm::real xzi,
cvm::real yxi, cvm::real yyi, cvm::real yzi,
cvm::real zxi, cvm::real zyi, cvm::real zzi)
: cvm::matrix2d<cvm::real>(3, 3)
{
this->xx() = xxi;
@ -983,31 +950,13 @@ public:
inline cvm::rmatrix transpose() const
{
return cvm::rmatrix(this->xx(),
this->yx(),
this->zx(),
this->xy(),
this->yy(),
this->zy(),
this->xz(),
this->yz(),
this->zz());
return cvm::rmatrix(this->xx(), this->yx(), this->zx(),
this->xy(), this->yy(), this->zy(),
this->xz(), this->yz(), this->zz());
}
friend cvm::rvector operator * (cvm::rmatrix const &m, cvm::rvector const &r);
// friend inline cvm::rmatrix const operator * (cvm::rmatrix const &m1, cvm::rmatrix const &m2) {
// return cvm::rmatrix (m1.xx()*m2.xx() + m1.xy()*m2.yx() + m1.xz()*m2.yz(),
// m1.xx()*m2.xy() + m1.xy()*m2.yy() + m1.xz()*m2.zy(),
// m1.xx()*m2.xz() + m1.xy()*m2.yz() + m1.xz()*m2.zz(),
// m1.yx()*m2.xx() + m1.yy()*m2.yx() + m1.yz()*m2.yz(),
// m1.yx()*m2.xy() + m1.yy()*m2.yy() + m1.yz()*m2.yy(),
// m1.yx()*m2.xz() + m1.yy()*m2.yz() + m1.yz()*m2.yz(),
// m1.zx()*m2.xx() + m1.zy()*m2.yx() + m1.zz()*m2.yz(),
// m1.zx()*m2.xy() + m1.zy()*m2.yy() + m1.zz()*m2.yy(),
// m1.zx()*m2.xz() + m1.zy()*m2.yz() + m1.zz()*m2.yz());
// }
};
@ -1031,7 +980,7 @@ public:
cvm::real q0, q1, q2, q3;
/// Constructor from a 3-d vector
inline quaternion(cvm::real const &x, cvm::real const &y, cvm::real const &z)
inline quaternion(cvm::real x, cvm::real y, cvm::real z)
: q0(0.0), q1(x), q2(y), q3(z)
{}
@ -1041,10 +990,10 @@ public:
{}
/// Constructor component by component
inline quaternion(cvm::real const &q0i,
cvm::real const &q1i,
cvm::real const &q2i,
cvm::real const &q3i)
inline quaternion(cvm::real q0i,
cvm::real q1i,
cvm::real q2i,
cvm::real q3i)
: q0(q0i), q1(q1i), q2(q2i), q3(q3i)
{}
@ -1055,9 +1004,9 @@ public:
/// "Constructor" after Euler angles (in radians)
///
/// http://en.wikipedia.org/wiki/Conversion_between_quaternions_and_Euler_angles
inline void set_from_euler_angles(cvm::real const &phi_in,
cvm::real const &theta_in,
cvm::real const &psi_in)
inline void set_from_euler_angles(cvm::real phi_in,
cvm::real theta_in,
cvm::real psi_in)
{
q0 = ( (std::cos(phi_in/2.0)) * (std::cos(theta_in/2.0)) * (std::cos(psi_in/2.0)) +
(std::sin(phi_in/2.0)) * (std::sin(theta_in/2.0)) * (std::sin(psi_in/2.0)) );
@ -1079,7 +1028,7 @@ public:
}
/// \brief Set all components to a scalar
inline void set(cvm::real const &value = 0.0)
inline void set(cvm::real value)
{
q0 = q1 = q2 = q3 = value;
}
@ -1087,7 +1036,7 @@ public:
/// \brief Set all components to zero (null quaternion)
inline void reset()
{
q0 = q1 = q2 = q3 = 0.0;
set(0.0);
}
/// \brief Set the q0 component to 1 and the others to 0 (quaternion
@ -1099,7 +1048,7 @@ public:
}
/// Tell the number of characters required to print a quaternion, given that of a real number
static inline size_t output_width(size_t const &real_width)
static inline size_t output_width(size_t real_width)
{
return 4*real_width + 13;
}
@ -1113,7 +1062,7 @@ public:
friend std::istream & operator >> (std::istream &is, cvm::quaternion &q);
/// Access the quaternion as a 4-d array (return a reference)
inline cvm::real & operator [] (int const &i) {
inline cvm::real & operator [] (int i) {
switch (i) {
case 0:
return this->q0;
@ -1130,7 +1079,7 @@ public:
}
/// Access the quaternion as a 4-d array (return a value)
inline cvm::real operator [] (int const &i) const {
inline cvm::real operator [] (int i) const {
switch (i) {
case 0:
return this->q0;
@ -1175,12 +1124,12 @@ public:
return cvm::quaternion(q0, -q1, -q2, -q3);
}
inline void operator *= (cvm::real const &a)
inline void operator *= (cvm::real a)
{
q0 *= a; q1 *= a; q2 *= a; q3 *= a;
}
inline void operator /= (cvm::real const &a)
inline void operator /= (cvm::real a)
{
q0 /= a; q1 /= a; q2 /= a; q3 /= a;
}
@ -1215,19 +1164,22 @@ public:
}
friend inline cvm::quaternion operator + (cvm::quaternion const &h, cvm::quaternion const &q)
friend inline cvm::quaternion operator + (cvm::quaternion const &h,
cvm::quaternion const &q)
{
return cvm::quaternion(h.q0+q.q0, h.q1+q.q1, h.q2+q.q2, h.q3+q.q3);
}
friend inline cvm::quaternion operator - (cvm::quaternion const &h, cvm::quaternion const &q)
friend inline cvm::quaternion operator - (cvm::quaternion const &h,
cvm::quaternion const &q)
{
return cvm::quaternion(h.q0-q.q0, h.q1-q.q1, h.q2-q.q2, h.q3-q.q3);
}
/// \brief Provides the quaternion product. \b NOTE: for the inner
/// product use: \code h.inner (q); \endcode
friend inline cvm::quaternion operator * (cvm::quaternion const &h, cvm::quaternion const &q)
friend inline cvm::quaternion operator * (cvm::quaternion const &h,
cvm::quaternion const &q)
{
return cvm::quaternion(h.q0*q.q0 - h.q1*q.q1 - h.q2*q.q2 - h.q3*q.q3,
h.q0*q.q1 + h.q1*q.q0 + h.q2*q.q3 - h.q3*q.q2,
@ -1235,18 +1187,18 @@ public:
h.q0*q.q3 + h.q3*q.q0 + h.q1*q.q2 - h.q2*q.q1);
}
friend inline cvm::quaternion operator * (cvm::real const &c,
friend inline cvm::quaternion operator * (cvm::real c,
cvm::quaternion const &q)
{
return cvm::quaternion(c*q.q0, c*q.q1, c*q.q2, c*q.q3);
}
friend inline cvm::quaternion operator * (cvm::quaternion const &q,
cvm::real const &c)
cvm::real c)
{
return cvm::quaternion(q.q0*c, q.q1*c, q.q2*c, q.q3*c);
}
friend inline cvm::quaternion operator / (cvm::quaternion const &q,
cvm::real const &c)
cvm::real c)
{
return cvm::quaternion(q.q0/c, q.q1/c, q.q2/c, q.q3/c);
}
@ -1407,7 +1359,7 @@ public:
std::vector< cvm::vector1d<cvm::rvector> > dQ0_1, dQ0_2;
/// Allocate space for the derivatives of the rotation
inline void request_group1_gradients(size_t const &n)
inline void request_group1_gradients(size_t n)
{
dS_1.resize(n, cvm::matrix2d<cvm::rvector>(4, 4));
dL0_1.resize(n, cvm::rvector(0.0, 0.0, 0.0));
@ -1415,7 +1367,7 @@ public:
}
/// Allocate space for the derivatives of the rotation
inline void request_group2_gradients(size_t const &n)
inline void request_group2_gradients(size_t n)
{
dS_2.resize(n, cvm::matrix2d<cvm::rvector>(4, 4));
dL0_2.resize(n, cvm::rvector(0.0, 0.0, 0.0));
@ -1448,7 +1400,7 @@ public:
}
/// Constructor after an axis of rotation and an angle (in radians)
inline rotation(cvm::real const &angle, cvm::rvector const &axis)
inline rotation(cvm::real angle, cvm::rvector const &axis)
: b_debug_gradients(false)
{
cvm::rvector const axis_n = axis.unit();
@ -1500,20 +1452,18 @@ public:
if (q.q0 != 0.0) {
// cvm::real const x = iprod/q.q0;
cvm::real const dspindx =
(180.0/PI) * 2.0 * (1.0 / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
cvm::real const dspindx = (180.0/PI) * 2.0 * (1.0 / (1.0 + (iprod*iprod)/(q.q0*q.q0)));
return
cvm::quaternion( dspindx * (iprod * (-1.0) / (q.q0*q.q0)),
dspindx * ((1.0/q.q0) * axis.x),
dspindx * ((1.0/q.q0) * axis.y),
dspindx * ((1.0/q.q0) * axis.z));
return cvm::quaternion( dspindx * (iprod * (-1.0) / (q.q0*q.q0)),
dspindx * ((1.0/q.q0) * axis.x),
dspindx * ((1.0/q.q0) * axis.y),
dspindx * ((1.0/q.q0) * axis.z));
} else {
// (1/(1+x^2)) ~ (1/x)^2
return
cvm::quaternion((180.0/PI) * 2.0 * ((-1.0)/iprod), 0.0, 0.0, 0.0);
// XX TODO: What if iprod == 0? XX
// The documentation of spinAngle discourages its use when q_vec and
// axis are not close
return cvm::quaternion((180.0/PI) * 2.0 * ((-1.0)/iprod), 0.0, 0.0, 0.0);
}
}

59
src/USER-COLVARS/colvarproxy_lammps.cpp
View File

@ -133,6 +133,8 @@ colvarproxy_lammps::colvarproxy_lammps(LAMMPS_NS::LAMMPS *lmp,
void colvarproxy_lammps::init(const char *conf_file)
{
version_int = get_version_from_string(COLVARPROXY_VERSION);
// create the colvarmodule instance
colvars = new colvarmodule(this);
@ -203,6 +205,25 @@ double colvarproxy_lammps::compute()
}
previous_step = _lmp->update->ntimestep;
unit_cell_x.set(_lmp->domain->xprd, 0.0, 0.0);
unit_cell_y.set(0.0, _lmp->domain->yprd, 0.0);
unit_cell_z.set(0.0, 0.0, _lmp->domain->zprd);
if (_lmp->domain->xperiodic == 0 && _lmp->domain->yperiodic == 0 &&
_lmp->domain->zperiodic == 0) {
boundaries_type = boundaries_non_periodic;
reset_pbc_lattice();
} else if ((_lmp->domain->nonperiodic == 0) &&
(_lmp->domain->dimension == 3) &&
(_lmp->domain->triclinic == 0)) {
// Orthogonal unit cell
boundaries_type = boundaries_pbc_ortho;
colvarproxy_system::update_pbc_lattice();
// It is safer to let LAMMPS deal with high-tilt triclinic boxes
} else {
boundaries_type = boundaries_unsupported;
}
if (cvm::debug()) {
cvm::log(std::string(cvm::line_marker)+
"colvarproxy_lammps, step no. "+cvm::to_str(colvars->it)+"\n"+
@ -263,8 +284,10 @@ bool colvarproxy_lammps::deserialize_status(std::string &rst)
}
}
cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
const
{
double xtmp = pos2.x - pos1.x;
double ytmp = pos2.y - pos1.y;
@ -273,28 +296,6 @@ cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1,
return cvm::rvector(xtmp, ytmp, ztmp);
}
cvm::real colvarproxy_lammps::position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2)
{
double xtmp = pos2.x - pos1.x;
double ytmp = pos2.y - pos1.y;
double ztmp = pos2.z - pos1.z;
_lmp->domain->minimum_image(xtmp,ytmp,ztmp);
return cvm::real(xtmp*xtmp + ytmp*ytmp + ztmp*ztmp);
}
void colvarproxy_lammps::select_closest_image(cvm::atom_pos &pos,
cvm::atom_pos const &ref)
{
double xtmp = pos.x - ref.x;
double ytmp = pos.y - ref.y;
double ztmp = pos.z - ref.z;
_lmp->domain->minimum_image(xtmp,ytmp,ztmp);
pos.x = ref.x + xtmp;
pos.y = ref.y + ytmp;
pos.z = ref.z + ztmp;
}
void colvarproxy_lammps::log(std::string const &message)
{
@ -308,12 +309,14 @@ void colvarproxy_lammps::log(std::string const &message)
}
}
void colvarproxy_lammps::error(std::string const &message)
{
// In LAMMPS, all errors are fatal
fatal_error(message);
}
void colvarproxy_lammps::fatal_error(std::string const &message)
{
log(message);
@ -321,13 +324,6 @@ void colvarproxy_lammps::fatal_error(std::string const &message)
"Fatal error in the collective variables module.\n");
}
void colvarproxy_lammps::exit(std::string const &message)
{
log(message);
log("Request to exit the simulation made.\n");
do_exit=true;
}
int colvarproxy_lammps::backup_file(char const *filename)
{
@ -342,10 +338,12 @@ int colvarproxy_lammps::backup_file(char const *filename)
// multi-replica support
void colvarproxy_lammps::replica_comm_barrier() {
void colvarproxy_lammps::replica_comm_barrier()
{
MPI_Barrier(inter_comm);
}
int colvarproxy_lammps::replica_comm_recv(char* msg_data,
int buf_len, int src_rep)
{
@ -359,6 +357,7 @@ int colvarproxy_lammps::replica_comm_recv(char* msg_data,
return retval;
}
int colvarproxy_lammps::replica_comm_send(char* msg_data,
int msg_len, int dest_rep)
{

7
src/USER-COLVARS/colvarproxy_lammps.h
View File

@ -120,14 +120,9 @@ class colvarproxy_lammps : public colvarproxy {
void log(std::string const &message);
void error(std::string const &message);
void fatal_error(std::string const &message);
void exit(std::string const &message);
cvm::rvector position_distance(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2);
cvm::real position_dist2(cvm::atom_pos const &pos1,
cvm::atom_pos const &pos2);
void select_closest_image(cvm::atom_pos &pos,
cvm::atom_pos const &ref_pos);
cvm::atom_pos const &pos2) const;
int backup_file(char const *filename);

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

@ -1,5 +1,5 @@
#ifndef COLVARPROXY_VERSION
#define COLVARPROXY_VERSION "2017-12-01"
#define COLVARPROXY_VERSION "2018-04-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