diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index c662982080..128762cb06 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -98,8 +98,8 @@ colvar::colvar (std::string const &conf) "on an axis", "distanceZ", distance_z); initialize_components ("distance projection " "on a plane", "distanceXY", distance_xy); - initialize_components ("average distance weighted by inverse sixth power", - "distance6", distance6); + initialize_components ("average distance weighted by inverse power", + "distanceInv", distance_inv); initialize_components ("coordination " "number", "coordNum", coordnum); diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index a9542a606c..396c63dc5a 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -471,7 +471,7 @@ public: class distance; class distance_z; class distance_xy; - class distance6; + class distance_inv; class angle; class dihedral; class coordnum; diff --git a/lib/colvars/colvaratoms.h b/lib/colvars/colvaratoms.h index 481a93a9bd..782ae6b9c8 100644 --- a/lib/colvars/colvaratoms.h +++ b/lib/colvars/colvaratoms.h @@ -1,3 +1,5 @@ +// -*- c++ -*- + #ifndef COLVARATOMS_H #define COLVARATOMS_H @@ -304,15 +306,8 @@ public: /// the colvar has not a scalar value) or the biases require to /// micromanage the forces. void apply_forces (std::vector const &forces); + }; - - #endif - - -// Emacs -// Local Variables: -// mode: C++ -// End: diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index 00d333dc89..7954fb5c1d 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -470,16 +470,16 @@ public: /// \brief Colvar component: average distance between two groups of atoms, weighted as the sixth power, /// as in NMR refinements (colvarvalue::type_scalar type, range (0:*)) -class colvar::distance6 +class colvar::distance_inv : public colvar::distance { protected: /// Components of the distance vector orthogonal to the axis - cvm::real smoothing; + int exponent; public: - distance6 (std::string const &conf); - distance6(); - virtual inline ~distance6() {} + distance_inv (std::string const &conf); + distance_inv(); + virtual inline ~distance_inv() {} virtual void calc_value(); virtual void calc_gradients(); virtual void apply_force (colvarvalue const &force); @@ -1296,7 +1296,7 @@ inline void colvar::spin_angle::wrap (colvarvalue &x) const simple_scalar_dist_functions (distance) // NOTE: distance_z has explicit functions, see below simple_scalar_dist_functions (distance_xy) - simple_scalar_dist_functions (distance6) + simple_scalar_dist_functions (distance_inv) simple_scalar_dist_functions (angle) simple_scalar_dist_functions (coordnum) simple_scalar_dist_functions (selfcoordnum) diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index 2b46538564..35d29efe4b 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -454,25 +454,34 @@ void colvar::distance_dir::apply_force (colvarvalue const &force) -colvar::distance6::distance6 (std::string const &conf) +colvar::distance_inv::distance_inv (std::string const &conf) : distance (conf) { - function_type = "distance6"; + function_type = "distance_inv"; + get_keyval (conf, "exponent", exponent, 6); + if (exponent%2) { + cvm::fatal_error ("Error: odd exponent provided, can only use even ones.\n"); + } + if (exponent <= 0) { + cvm::fatal_error ("Error: negative or zero exponent provided.\n"); + } + b_inverse_gradients = false; b_Jacobian_derivative = false; x.type (colvarvalue::type_scalar); } -colvar::distance6::distance6() +colvar::distance_inv::distance_inv() { - function_type = "distance6"; + function_type = "distance_inv"; + exponent = 6; b_inverse_gradients = false; b_Jacobian_derivative = false; b_1site_force = false; x.type (colvarvalue::type_scalar); } -void colvar::distance6::calc_value() +void colvar::distance_inv::calc_value() { group1.reset_atoms_data(); group2.reset_atoms_data(); @@ -486,8 +495,11 @@ void colvar::distance6::calc_value() for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { cvm::rvector const dv = ai2->pos - ai1->pos; cvm::real const d2 = dv.norm2(); - x.real_value += 1.0/(d2*d2*d2); - cvm::rvector const dsumddv = -6.0/(d2*d2*d2*d2) * dv; + cvm::real dinv = 1.0; + for (int ne = 0; ne < exponent/2; ne++) + dinv *= 1.0/d2; + x.real_value += dinv; + cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv; ai1->grad += -1.0 * dsumddv; ai2->grad += dsumddv; } @@ -497,8 +509,11 @@ void colvar::distance6::calc_value() for (cvm::atom_iter ai2 = group2.begin(); ai2 != group2.end(); ai2++) { cvm::rvector const dv = cvm::position_distance (ai1->pos, ai2->pos); cvm::real const d2 = dv.norm2(); - x.real_value += 1.0/(d2*d2*d2); - cvm::rvector const dsumddv = -6.0/(d2*d2*d2*d2) * dv; + cvm::real dinv = 1.0; + for (int ne = 0; ne < exponent/2; ne++) + dinv *= 1.0/d2; + x.real_value += dinv; + cvm::rvector const dsumddv = -(cvm::real (exponent)) * dinv/d2 * dv; ai1->grad += -1.0 * dsumddv; ai2->grad += dsumddv; } @@ -506,12 +521,12 @@ void colvar::distance6::calc_value() } x.real_value *= 1.0 / cvm::real (group1.size() * group2.size()); - x.real_value = std::pow (x.real_value, -1.0/6.0); + x.real_value = std::pow (x.real_value, -1.0/(cvm::real (exponent))); } -void colvar::distance6::calc_gradients() +void colvar::distance_inv::calc_gradients() { - cvm::real const dxdsum = (-1.0/6.0) * std::pow (x.real_value, 7.0) / cvm::real (group1.size() * group2.size()); + cvm::real const dxdsum = (-1.0/(cvm::real (exponent))) * std::pow (x.real_value, exponent+1) / cvm::real (group1.size() * group2.size()); for (cvm::atom_iter ai1 = group1.begin(); ai1 != group1.end(); ai1++) { ai1->grad *= dxdsum; } @@ -520,7 +535,7 @@ void colvar::distance6::calc_gradients() } } -void colvar::distance6::apply_force (colvarvalue const &force) +void colvar::distance_inv::apply_force (colvarvalue const &force) { if (!group1.noforce) group1.apply_colvar_force (force.real_value); diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index f587f6a321..2dfc1bd2ac 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -2,7 +2,7 @@ #define COLVARMODULE_H #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2012-07-17" +#define COLVARS_VERSION "2012-08-08" #endif #ifndef COLVARS_DEBUG