diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index fbbbf0ba0d..25b0f20b36 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -107,6 +107,8 @@ colvar::colvar(std::string const &conf) "on a plane", "distanceXY", distance_xy); initialize_components("average distance weighted by inverse power", "distanceInv", distance_inv); + initialize_components("N1xN2-long vector of pairwise distances", + "distancePairs", distance_pairs); initialize_components("coordination " "number", "coordNum", coordnum); @@ -163,7 +165,7 @@ colvar::colvar(std::string const &conf) // might accept other types when the infrastructure is in place // for derivatives of vectors wrt vectors x.type(colvarvalue::type_scalar); - x_reported.type(x.type()); + x_reported.type(x); // Sort array of cvcs based on values of componentExp std::vector temp_vec; @@ -206,14 +208,15 @@ colvar::colvar(std::string const &conf) } if (!tasks[task_scripted]) { - colvarvalue::Type const value_type = (cvcs[0])->type(); + colvarvalue const &cvc_value = (cvcs[0])->value(); if (cvm::debug()) cvm::log ("This collective variable is a "+ - colvarvalue::type_desc(value_type)+", corresponding to "+ - cvm::to_str (colvarvalue::num_df(value_type))+ - " internal degrees of freedom.\n"); - x.type(value_type); - x_reported.type(value_type); + colvarvalue::type_desc(cvc_value.type())+ + ((cvc_value.size() > 1) ? " with "+ + cvm::to_str(cvc_value.size())+" individual components.\n" : + ".\n")); + x.type(cvc_value); + x_reported.type(cvc_value); } // If using scripted biases, any colvar may receive bias forces @@ -292,13 +295,10 @@ colvar::colvar(std::string const &conf) b_Jacobian_force = false; // components may have different types only for scripted functions - if (!tasks[task_scripted] && (cvcs[i])->type() != (cvcs[0])->type() ) { + if (!tasks[task_scripted] && (colvarvalue::check_types(cvcs[i]->value(), + cvcs[0]->value())) ) { cvm::error("ERROR: you are definining this collective variable " - "by using components of different types, \""+ - colvarvalue::type_desc((cvcs[0])->type())+ - "\" and \""+ - colvarvalue::type_desc((cvcs[i])->type())+ - "\". " + "by using components of different types. " "You must use the same type in order to " " sum them together.\n", INPUT_ERROR); return; @@ -310,13 +310,13 @@ colvar::colvar(std::string const &conf) cvm::error("Error: \"width\" must be positive.\n", INPUT_ERROR); } - lower_boundary.type(this->type()); - lower_wall.type(this->type()); + lower_boundary.type(value()); + lower_wall.type(value()); - upper_boundary.type(this->type()); - upper_wall.type(this->type()); + upper_boundary.type(value()); + upper_wall.type(value()); - if (this->type() == colvarvalue::type_scalar) { + if (value().type() == colvarvalue::type_scalar) { if (get_keyval(conf, "lowerBoundary", lower_boundary, lower_boundary)) { enable(task_lower_boundary); @@ -399,9 +399,9 @@ colvar::colvar(std::string const &conf) enable(task_extended_lagrangian); - xr.type(this->type()); - vr.type(this->type()); - fr.type(this->type()); + xr.type(value()); + vr.type(value()); + fr.type(value()); const bool found = get_keyval(conf, "extendedTemp", temp, cvm::temperature()); if (temp <= 0.0) { @@ -645,7 +645,7 @@ int colvar::enable(colvar::task const &t) cvm::log("Enabling calculation of the Jacobian force " "on this colvar.\n"); } - fj.type(this->type()); + fj.type(value()); break; case task_system_force: @@ -658,8 +658,8 @@ int colvar::enable(colvar::task const &t) } cvm::request_system_force(); } - ft.type(this->type()); - ft_reported.type(this->type()); + ft.type(value()); + ft_reported.type(value()); break; case task_output_applied_force: @@ -670,9 +670,9 @@ int colvar::enable(colvar::task const &t) break; case task_fdiff_velocity: - x_old.type(this->type()); - v_fdiff.type(this->type()); - v_reported.type(this->type()); + x_old.type(value()); + v_fdiff.type(value()); + v_reported.type(value()); break; case task_output_velocity: @@ -680,7 +680,7 @@ int colvar::enable(colvar::task const &t) break; case task_grid: - if (this->type() != colvarvalue::type_scalar) { + if (value().type() != colvarvalue::type_scalar) { cvm::error("Cannot calculate a grid for collective variable, \""+ this->name+"\", because its value is not a scalar number.\n", INPUT_ERROR); @@ -689,12 +689,12 @@ int colvar::enable(colvar::task const &t) case task_extended_lagrangian: enable(task_gradients); - v_reported.type(this->type()); + v_reported.type(value()); break; case task_lower_boundary: case task_upper_boundary: - if (this->type() != colvarvalue::type_scalar) { + if (value().type() != colvarvalue::type_scalar) { cvm::error("Error: this colvar is not a scalar value " "and cannot produce a grid.\n", INPUT_ERROR); @@ -711,12 +711,12 @@ int colvar::enable(colvar::task const &t) break; case task_gradients: - f.type(this->type()); - fb.type(this->type()); + f.type(value()); + fb.type(value()); break; case task_collect_gradients: - if (this->type() != colvarvalue::type_scalar) { + if (value().type() != colvarvalue::type_scalar) { cvm::error("Collecting atomic gradients for non-scalar collective variable \""+ this->name+"\" is not yet implemented.\n", INPUT_ERROR); @@ -1056,8 +1056,8 @@ cvm::real colvar::update() if (cvm::debug()) cvm::log("Updating colvar \""+this->name+"\".\n"); - // set to zero the applied force + f.type(value()); f.reset(); // add the biases' force, which at this point should already have @@ -1068,7 +1068,11 @@ cvm::real colvar::update() if (tasks[task_lower_wall] || tasks[task_upper_wall]) { // wall force - colvarvalue fw(this->type()); + colvarvalue fw(value()); + fw.reset(); + + if (cvm::debug()) + cvm::log("Calculating wall forces for colvar \""+this->name+"\".\n"); // if the two walls are applied concurrently, decide which is the // closer one (on a periodic colvar, both walls may be applicable @@ -1313,36 +1317,36 @@ std::istream & colvar::read_restart(std::istream &is) { std::string check_name = ""; if ( (get_keyval(conf, "name", check_name, - std::string(""), colvarparse::parse_silent)) && + std::string(""), colvarparse::parse_silent)) && (check_name != name) ) { cvm::error("Error: the state file does not match the " - "configuration file, at colvar \""+name+"\".\n"); + "configuration file, at colvar \""+name+"\".\n"); } if (check_name.size() == 0) { cvm::error("Error: Collective variable in the " - "restart file without any identifier.\n"); + "restart file without any identifier.\n"); } } if ( !(get_keyval(conf, "x", x, - colvarvalue(x.type()), colvarparse::parse_silent)) ) { + colvarvalue(x.type()), colvarparse::parse_silent)) ) { cvm::log("Error: restart file does not contain " - "the value of the colvar \""+ - name+"\" .\n"); + "the value of the colvar \""+ + name+"\" .\n"); } else { cvm::log("Restarting collective variable \""+name+"\" from value: "+ - cvm::to_str(x)+"\n"); + cvm::to_str(x)+"\n"); } if (tasks[colvar::task_extended_lagrangian]) { if ( !(get_keyval(conf, "extended_x", xr, - colvarvalue(x.type()), colvarparse::parse_silent)) && + colvarvalue(x.type()), colvarparse::parse_silent)) && !(get_keyval(conf, "extended_v", vr, - colvarvalue(x.type()), colvarparse::parse_silent)) ) { + colvarvalue(x.type()), colvarparse::parse_silent)) ) { cvm::log("Error: restart file does not contain " - "\"extended_x\" or \"extended_v\" for the colvar \""+ - name+"\", but you requested \"extendedLagrangian\".\n"); + "\"extended_x\" or \"extended_v\" for the colvar \""+ + name+"\", but you requested \"extendedLagrangian\".\n"); } } @@ -1355,10 +1359,10 @@ std::istream & colvar::read_restart(std::istream &is) if (tasks[task_output_velocity]) { if ( !(get_keyval(conf, "v", v_fdiff, - colvarvalue(x.type()), colvarparse::parse_silent)) ) { + colvarvalue(x.type()), colvarparse::parse_silent)) ) { cvm::log("Error: restart file does not contain " - "the velocity for the colvar \""+ - name+"\", but you requested \"outputVelocity\".\n"); + "the velocity for the colvar \""+ + name+"\", but you requested \"outputVelocity\".\n"); } if (tasks[task_extended_lagrangian]) { @@ -1636,11 +1640,11 @@ int colvar::calc_acf() colvar *cfcv = (acf_colvar_name.size() ? cvm::colvar_by_name(acf_colvar_name) : this); - if (cfcv->type() != this->type()) { + if (colvarvalue::check_types(cfcv->value(), value())) { cvm::error("Error: correlation function between \""+cfcv->name+ - "\" and \""+this->name+"\" cannot be calculated, " - "because their value types are different.\n", - INPUT_ERROR); + "\" and \""+this->name+"\" cannot be calculated, " + "because their value types are different.\n", + INPUT_ERROR); } acf_nframes = 0; @@ -1823,7 +1827,7 @@ void colvar::calc_runave() { if (x_history.empty()) { - runave.type(x.type()); + runave.type(value().type()); runave.reset(); // first-step operations diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index e144fea500..80c5a530a4 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -31,9 +31,8 @@ /// /// Please note that most of its members are \link colvarvalue /// \endlink objects, i.e. they can handle different data types -/// together, and must all be set to the same type of colvar::x by -/// using the colvarvalue::type() member function before using them -/// together in assignments or other operations; this is usually done +/// together, and must all be set to the same type of colvar::value() +/// before using them together in assignments or other operations; this is usually done /// automatically in the constructor. If you add a new member of /// \link colvarvalue \endlink type, you should also add its /// initialization line in the \link colvar \endlink constructor. @@ -45,16 +44,13 @@ public: /// Name std::string name; - /// Type of value - colvarvalue::Type type() const; - - /// \brief Current value (previously obtained from calc() or read_traj()) + /// \brief Current value (previously set by calc() or by read_traj()) colvarvalue const & value() const; /// \brief Current actual value (not extended DOF) colvarvalue const & actual_value() const; - /// \brief Current velocity (previously obtained from calc() or read_traj()) + /// \brief Current velocity (previously set by calc() or by read_traj()) colvarvalue const & velocity() const; /// \brief Current system force (previously obtained from calc() or @@ -486,6 +482,7 @@ public: class distance_z; class distance_xy; class distance_inv; + class distance_pairs; class angle; class dihedral; class coordnum; @@ -540,11 +537,6 @@ public: } }; -inline colvarvalue::Type colvar::type() const -{ - return x.type(); -} - inline colvarvalue const & colvar::value() const { @@ -577,6 +569,7 @@ inline void colvar::add_bias_force(colvarvalue const &force) inline void colvar::reset_bias_force() { + fb.type(value()); fb.reset(); } diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index 08bc82bb08..8f1d85b242 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -83,17 +83,19 @@ colvarbias::~colvarbias() void colvarbias::add_colvar(std::string const &cv_name) { - if (colvar *cvp = cvm::colvar_by_name(cv_name)) { - cvp->enable(colvar::task_gradients); + if (colvar *cv = cvm::colvar_by_name(cv_name)) { + cv->enable(colvar::task_gradients); if (cvm::debug()) cvm::log("Applying this bias to collective variable \""+ - cvp->name+"\".\n"); - colvars.push_back(cvp); - colvar_forces.push_back(colvarvalue(cvp->type())); - cvp->biases.push_back(this); // add back-reference to this bias to colvar + cv->name+"\".\n"); + colvars.push_back(cv); + colvar_forces.push_back(colvarvalue()); + colvar_forces.back().type(cv->value()); // make sure each forces is initialized to zero + colvar_forces.back().reset(); + cv->biases.push_back(this); // add back-reference to this bias to colvar } else { cvm::error("Error: cannot find a colvar named \""+ - cv_name+"\".\n"); + cv_name+"\".\n"); } } @@ -103,8 +105,7 @@ void colvarbias::communicate_forces() for (size_t i = 0; i < colvars.size(); i++) { if (cvm::debug()) { cvm::log("Communicating a force to colvar \""+ - colvars[i]->name+"\", of type \""+ - colvarvalue::type_desc(colvars[i]->type())+"\".\n"); + colvars[i]->name+"\".\n"); } colvars[i]->add_bias_force(colvar_forces[i]); } diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index 4c53de04c7..543be6a95f 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -64,7 +64,7 @@ colvarbias_abf::colvarbias_abf(std::string const &conf, char const *key) for (size_t i = 0; i < colvars.size(); i++) { - if (colvars[i]->type() != colvarvalue::type_scalar) { + if (colvars[i]->value().type() != colvarvalue::type_scalar) { cvm::error("Error: ABF bias can only use scalar-type variables.\n"); } diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp index 605ac39b62..309a089258 100644 --- a/lib/colvars/colvarbias_alb.cpp +++ b/lib/colvars/colvarbias_alb.cpp @@ -31,7 +31,7 @@ colvarbias_alb::colvarbias_alb(std::string const &conf, char const *key) : for (i = 0; i < colvars.size(); i++) { - colvar_centers[i].type(colvars[i]->type()); + colvar_centers[i].type(colvars[i]->value()); //zero moments means[i] = ssd[i] = 0; diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 21c995b311..e8e104366b 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -646,7 +646,7 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, { std::vector curr_values(colvars.size()); for (size_t i = 0; i < colvars.size(); i++) { - curr_values[i].type(colvars[i]->type()); + curr_values[i].type(colvars[i]->value()); } if (colvar_values.size()) { @@ -689,15 +689,14 @@ void colvarbias_meta::calc_hills_force(size_t const &i, std::vector const &values) { // Retrieve the value of the colvar - colvarvalue x(values.size() ? values[i].type() : colvars[i]->type()); - x = (values.size() ? values[i] : colvars[i]->value()); + colvarvalue const x(values.size() ? values[i] : colvars[i]->value()); // do the type check only once (all colvarvalues in the hills series // were already saved with their types matching those in the // colvars) hill_iter h; - switch (colvars[i]->type()) { + switch (colvars[i]->value().type()) { case colvarvalue::type_scalar: for (h = h_first; h != h_last; h++) { @@ -1433,7 +1432,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) std::vector h_centers(colvars.size()); for (size_t i = 0; i < colvars.size(); i++) { - h_centers[i].type((colvars[i]->value()).type()); + h_centers[i].type(colvars[i]->value()); } { // it is safer to read colvarvalue objects one at a time; diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index 81ee6266ef..de5b6d2622 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -280,7 +280,7 @@ public: replica(replica_in) { for (size_t i = 0; i < cv.size(); i++) { - centers[i].type(cv[i]->type()); + centers[i].type(cv[i]->value()); centers[i] = cv[i]->value(); widths[i] = cv[i]->width * hill_width; } diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index 5b8049d6a7..76586291e2 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -6,42 +6,61 @@ colvarbias_restraint::colvarbias_restraint(std::string const &conf, - char const *key) + char const *key) : colvarbias(conf, key), target_nstages(0), target_nsteps(0) { get_keyval(conf, "forceConstant", force_k, 1.0); - // get the initial restraint centers - colvar_centers.resize(colvars.size()); - colvar_centers_raw.resize(colvars.size()); - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers[i].type(colvars[i]->type()); - colvar_centers_raw[i].type(colvars[i]->type()); - } - if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { - for (size_t i = 0; i < colvars.size(); i++) { - colvar_centers[i].apply_constraints(); - colvar_centers_raw[i] = colvar_centers[i]; + { + // get the initial restraint centers + colvar_centers.resize(colvars.size()); + colvar_centers_raw.resize(colvars.size()); + size_t i; + for (i = 0; i < colvars.size(); i++) { + colvar_centers[i].type(colvars[i]->value()); + colvar_centers_raw[i].type(colvars[i]->value()); + if (cvm::debug()) { + cvm::log("colvarbias_restraint: center size = "+ + cvm::to_str(colvar_centers[i].vector1d_value.size())+"\n"); + } + } + if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { + for (i = 0; i < colvars.size(); i++) { + if (cvm::debug()) { + cvm::log("colvarbias_restraint: parsing initial centers, i = "+cvm::to_str(i)+".\n"); + } + + colvar_centers[i].apply_constraints(); + colvar_centers_raw[i] = colvar_centers[i]; + } + } else { + colvar_centers.clear(); + cvm::error("Error: must define the initial centers of the restraints.\n"); + } + + if (colvar_centers.size() != colvars.size()) { + cvm::error("Error: number of centers does not match " + "that of collective variables.\n"); } - } else { - colvar_centers.clear(); - cvm::error("Error: must define the initial centers of the restraints.\n"); } - if (colvar_centers.size() != colvars.size()) - cvm::error("Error: number of centers does not match " - "that of collective variables.\n"); - - if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) { - b_chg_centers = true; - for (size_t i = 0; i < target_centers.size(); i++) { - target_centers[i].apply_constraints(); + { + if (cvm::debug()) { + cvm::log("colvarbias_restraint: parsing target centers.\n"); + } + + size_t i; + if (get_keyval(conf, "targetCenters", target_centers, colvar_centers)) { + b_chg_centers = true; + for (i = 0; i < target_centers.size(); i++) { + target_centers[i].apply_constraints(); + } + } else { + b_chg_centers = false; + target_centers.clear(); } - } else { - b_chg_centers = false; - target_centers.clear(); } if (get_keyval(conf, "targetForceConstant", target_force_k, 0.0)) { @@ -110,7 +129,9 @@ void colvarbias_restraint::change_configuration(std::string const &conf) get_keyval(conf, "forceConstant", force_k, force_k); if (get_keyval(conf, "centers", colvar_centers, colvar_centers)) { for (size_t i = 0; i < colvars.size(); i++) { + colvar_centers[i].type(colvars[i]->value()); colvar_centers[i].apply_constraints(); + colvar_centers_raw[i].type(colvars[i]->value()); colvar_centers_raw[i] = colvar_centers[i]; } } @@ -128,10 +149,11 @@ cvm::real colvarbias_restraint::energy_difference(std::string const &conf) alt_colvar_centers.resize(colvars.size()); size_t i; for (i = 0; i < colvars.size(); i++) { - alt_colvar_centers[i].type(colvars[i]->type()); + alt_colvar_centers[i].type(colvars[i]->value()); } if (get_keyval(conf, "centers", alt_colvar_centers, colvar_centers)) { for (i = 0; i < colvars.size(); i++) { + colvar_centers[i].type(colvars[i]->value()); colvar_centers[i].apply_constraints(); } } @@ -177,7 +199,7 @@ cvm::real colvarbias_restraint::update() // if we are restarting a staged calculation) centers_incr.resize(colvars.size()); for (size_t i = 0; i < colvars.size(); i++) { - centers_incr[i].type(colvars[i]->type()); + centers_incr[i].type(colvars[i]->value()); centers_incr[i] = (target_centers[i] - colvar_centers_raw[i]) / cvm::real( target_nstages ? (target_nstages - stage) : (target_nsteps - cvm::step_absolute())); @@ -285,6 +307,7 @@ cvm::real colvarbias_restraint::update() // Force and energy calculation for (size_t i = 0; i < colvars.size(); i++) { + colvar_forces[i].type(colvars[i]->value()); colvar_forces[i] = -1.0 * restraint_force(restraint_convert_k(force_k, colvars[i]->width), colvars[i], colvar_centers[i]); diff --git a/lib/colvars/colvarcomp.cpp b/lib/colvars/colvarcomp.cpp index 7d0ccff933..fc31c44ab5 100644 --- a/lib/colvars/colvarcomp.cpp +++ b/lib/colvars/colvarcomp.cpp @@ -90,7 +90,8 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group) cvm::rotation const rot_0 = group.rot; cvm::rotation const rot_inv = group.rot.inverse(); - cvm::real const x_0 = x.real_value; + cvm::real x_0 = x.real_value; + if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_0 = x[0]; // cvm::log("gradients = "+cvm::to_str (gradients)+"\n"); @@ -130,7 +131,8 @@ void colvar::cvc::debug_gradients(cvm::atom_group &group) group.calc_apply_roto_translation(); } calc_value(); - cvm::real const x_1 = x.real_value; + cvm::real x_1 = x.real_value; + if ((x.type() == colvarvalue::type_vector) && (x.size() == 1)) x_1 = x[0]; cvm::log("Atom "+cvm::to_str(ia)+", component "+cvm::to_str(id)+":\n"); cvm::log("dx(actual) = "+cvm::to_str(x_1 - x_0, 21, 14)+"\n"); diff --git a/lib/colvars/colvarcomp.h b/lib/colvars/colvarcomp.h index a61a27f888..6137172954 100644 --- a/lib/colvars/colvarcomp.h +++ b/lib/colvars/colvarcomp.h @@ -79,10 +79,6 @@ public: /// this variable definition should be set within the constructor. std::string function_type; - /// \brief Type of \link colvarvalue \endlink that this cvc - /// provides - colvarvalue::Type type() const; - /// \brief Coefficient in the polynomial combination (default: 1.0) cvm::real sup_coeff; /// \brief Exponent in the polynomial combination (default: 1) @@ -248,11 +244,6 @@ protected: -inline colvarvalue::Type colvar::cvc::type() const -{ - return x.type(); -} - inline colvarvalue const & colvar::cvc::value() const { return x; @@ -483,6 +474,34 @@ public: }; +/// \brief Colvar component: N1xN2 vector of pairwise distances +/// (colvarvalue::type_vector type, range (0:*) for each component) +class colvar::distance_pairs + : public colvar::cvc +{ +protected: + /// First atom group + cvm::atom_group group1; + /// Second atom group + cvm::atom_group group2; + /// Use absolute positions, ignoring PBCs when present + bool b_no_PBC; +public: + distance_pairs(std::string const &conf); + distance_pairs(); + virtual inline ~distance_pairs() {} + virtual void calc_value(); + virtual void calc_gradients(); + virtual void apply_force(colvarvalue const &force); + 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; +}; + + /// \brief Colvar component: Radius of gyration of an atom group /// (colvarvalue::type_scalar type, range [0:*)) @@ -1335,6 +1354,25 @@ inline colvarvalue colvar::distance_dir::dist2_rgrad(colvarvalue const &x1, return colvarvalue((x2.rvector_value - x1.rvector_value), colvarvalue::type_unit3vector); } +inline cvm::real colvar::distance_pairs::dist2(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return (x1.vector1d_value - x2.vector1d_value).norm2(); +} + +inline colvarvalue colvar::distance_pairs::dist2_lgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return colvarvalue((x1.vector1d_value - x2.vector1d_value), colvarvalue::type_vector); +} + +inline colvarvalue colvar::distance_pairs::dist2_rgrad(colvarvalue const &x1, + colvarvalue const &x2) const +{ + return colvarvalue((x2.vector1d_value - x1.vector1d_value), colvarvalue::type_vector); +} + + // distance between quaternions inline cvm::real colvar::orientation::dist2(colvarvalue const &x1, diff --git a/lib/colvars/colvarcomp_distances.cpp b/lib/colvars/colvarcomp_distances.cpp index c697cb53a7..c0f032ddd4 100644 --- a/lib/colvars/colvarcomp_distances.cpp +++ b/lib/colvars/colvarcomp_distances.cpp @@ -527,6 +527,99 @@ void colvar::distance_inv::apply_force(colvarvalue const &force) +colvar::distance_pairs::distance_pairs(std::string const &conf) + : cvc(conf) +{ + function_type = "distance_pairs"; + b_inverse_gradients = false; + b_Jacobian_derivative = false; + + if (get_keyval(conf, "forceNoPBC", b_no_PBC, false)) { + cvm::log("Computing distance using absolute positions (not minimal-image)"); + } + + parse_group(conf, "group1", group1); + atom_groups.push_back(&group1); + + parse_group(conf, "group2", group2); + atom_groups.push_back(&group2); + + x.type(colvarvalue::type_vector); + x.vector1d_value.resize(group1.size() * group2.size()); +} + + +colvar::distance_pairs::distance_pairs() +{ + function_type = "distance_pairs"; + b_inverse_gradients = false; + b_Jacobian_derivative = false; + x.type(colvarvalue::type_vector); +} + + +void colvar::distance_pairs::calc_value() +{ + x.vector1d_value.resize(group1.size() * group2.size()); + + if (b_no_PBC) { + size_t i1, i2; + for (i1 = 0; i1 < group1.size(); i1++) { + for (i2 = 0; i2 < group2.size(); i2++) { + cvm::rvector const dv = group2[i2].pos - group1[i1].pos; + cvm::real const d = dv.norm(); + x.vector1d_value[i1*group2.size() + i2] = d; + group1[i1].grad = -1.0 * dv.unit(); + group2[i2].grad = dv.unit(); + } + } + } else { + size_t i1, i2; + for (i1 = 0; i1 < group1.size(); i1++) { + for (i2 = 0; i2 < group2.size(); i2++) { + cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos); + cvm::real const d = dv.norm(); + x.vector1d_value[i1*group2.size() + i2] = d; + group1[i1].grad = -1.0 * dv.unit(); + group2[i2].grad = dv.unit(); + } + } + } +} + +void colvar::distance_pairs::calc_gradients() +{ + // will be calculated on the fly in apply_force() + if (b_debug_gradients) { + cvm::log("Debugging gradients:\n"); + debug_gradients(group1); + } +} + +void colvar::distance_pairs::apply_force(colvarvalue const &force) +{ + if (b_no_PBC) { + size_t i1, i2; + for (i1 = 0; i1 < group1.size(); i1++) { + for (i2 = 0; i2 < group2.size(); i2++) { + cvm::rvector const dv = group2[i2].pos - group1[i1].pos; + group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit()); + group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit()); + } + } + } else { + size_t i1, i2; + for (i1 = 0; i1 < group1.size(); i1++) { + for (i2 = 0; i2 < group2.size(); i2++) { + cvm::rvector const dv = cvm::position_distance(group1[i1].pos, group2[i2].pos); + group1[i1].apply_force(force[i1*group2.size() + i2] * (-1.0) * dv.unit()); + group2[i2].apply_force(force[i1*group2.size() + i2] * dv.unit()); + } + } + } +} + + colvar::gyration::gyration(std::string const &conf) : cvc(conf) { diff --git a/lib/colvars/colvargrid.h b/lib/colvars/colvargrid.h index e058b2f2da..17ff6a4fb2 100644 --- a/lib/colvars/colvargrid.h +++ b/lib/colvars/colvargrid.h @@ -223,7 +223,7 @@ public: for (size_t i = 0; i < cv.size(); i++) { - if (cv[i]->type() != colvarvalue::type_scalar) { + if (cv[i]->value().type() != colvarvalue::type_scalar) { cvm::error("Colvar grids can only be automatically " "constructed for scalar variables. " "ABF and histogram can not be used; metadynamics " diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index 463b23e4de..291c0108e4 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -19,11 +19,11 @@ size_t colvarparse::dummy_pos = 0; #define _get_keyval_scalar_string_(TYPE) \ \ - bool colvarparse::get_keyval(std::string const &conf, \ - char const *key, \ - TYPE &value, \ - TYPE const &def_value, \ - Parse_Mode const parse_mode) \ + bool colvarparse::get_keyval(std::string const &conf, \ + char const *key, \ + TYPE &value, \ + TYPE const &def_value, \ + Parse_Mode const parse_mode) \ { \ std::string data; \ bool b_found = false, b_found_any = false; \ @@ -31,7 +31,7 @@ size_t colvarparse::dummy_pos = 0; \ do { \ std::string data_this = ""; \ - b_found = key_lookup(conf, key, data_this, save_pos); \ + b_found = key_lookup(conf, key, data_this, save_pos); \ if (b_found) { \ if (!b_found_any) \ b_found_any = true; \ @@ -41,30 +41,30 @@ size_t colvarparse::dummy_pos = 0; } while (b_found); \ \ if (found_count > 1) \ - cvm::log("Warning: found more than one instance of \""+ \ - std::string(key)+"\".\n"); \ + cvm::log("Warning: found more than one instance of \""+ \ + std::string(key)+"\".\n"); \ \ if (data.size()) { \ - std::istringstream is(data); \ - TYPE x(def_value); \ + std::istringstream is(data); \ + TYPE x(def_value); \ if (is >> x) \ value = x; \ else \ - cvm::error("Error: in parsing \""+ \ - std::string(key)+"\".\n", INPUT_ERROR); \ + cvm::error("Error: in parsing \""+ \ + std::string(key)+"\".\n", INPUT_ERROR); \ if (parse_mode != parse_silent) { \ - cvm::log("# "+std::string(key)+" = "+ \ - cvm::to_str(value)+"\n"); \ + cvm::log("# "+std::string(key)+" = "+ \ + cvm::to_str(value)+"\n"); \ } \ } else { \ \ if (b_found_any) \ - cvm::error("Error: improper or missing value " \ - "for \""+std::string(key)+"\".\n", INPUT_ERROR); \ + cvm::error("Error: improper or missing value " \ + "for \""+std::string(key)+"\".\n", INPUT_ERROR); \ value = def_value; \ if (parse_mode != parse_silent) { \ - cvm::log("# "+std::string(key)+" = \""+ \ - cvm::to_str(def_value)+"\" [default]\n"); \ + cvm::log("# "+std::string(key)+" = \""+ \ + cvm::to_str(def_value)+"\" [default]\n"); \ } \ } \ \ @@ -74,11 +74,11 @@ size_t colvarparse::dummy_pos = 0; #define _get_keyval_scalar_(TYPE) \ \ - bool colvarparse::get_keyval(std::string const &conf, \ - char const *key, \ - TYPE &value, \ - TYPE const &def_value, \ - Parse_Mode const parse_mode) \ + bool colvarparse::get_keyval(std::string const &conf, \ + char const *key, \ + TYPE &value, \ + TYPE const &def_value, \ + Parse_Mode const parse_mode) \ { \ std::string data; \ bool b_found = false, b_found_any = false; \ @@ -86,7 +86,7 @@ size_t colvarparse::dummy_pos = 0; \ do { \ std::string data_this = ""; \ - b_found = key_lookup(conf, key, data_this, save_pos); \ + b_found = key_lookup(conf, key, data_this, save_pos); \ if (b_found) { \ if (!b_found_any) \ b_found_any = true; \ @@ -96,37 +96,37 @@ size_t colvarparse::dummy_pos = 0; } while (b_found); \ \ if (found_count > 1) \ - cvm::log("Warning: found more than one instance of \""+ \ - std::string(key)+"\".\n"); \ + cvm::log("Warning: found more than one instance of \""+ \ + std::string(key)+"\".\n"); \ \ if (data.size()) { \ - std::istringstream is(data); \ + std::istringstream is(data); \ size_t data_count = 0; \ - TYPE x(def_value); \ + TYPE x(def_value); \ while (is >> x) { \ value = x; \ data_count++; \ } \ if (data_count == 0) \ - cvm::fatal_error("Error: in parsing \""+ \ - std::string(key)+"\".\n"); \ + cvm::fatal_error("Error: in parsing \""+ \ + std::string(key)+"\".\n"); \ if (data_count > 1) \ - cvm::error("Error: multiple values " \ - "are not allowed for keyword \""+ \ - std::string(key)+"\".\n", INPUT_ERROR); \ + cvm::error("Error: multiple values " \ + "are not allowed for keyword \""+ \ + std::string(key)+"\".\n", INPUT_ERROR); \ if (parse_mode != parse_silent) { \ - cvm::log("# "+std::string(key)+" = "+ \ - cvm::to_str(value)+"\n"); \ + cvm::log("# "+std::string(key)+" = "+ \ + cvm::to_str(value)+"\n"); \ } \ } else { \ \ if (b_found_any) \ cvm::error("Error: improper or missing value " \ - "for \""+std::string(key)+"\".\n", INPUT_ERROR); \ + "for \""+std::string(key)+"\".\n", INPUT_ERROR); \ value = def_value; \ if (parse_mode != parse_silent) { \ - cvm::log("# "+std::string(key)+" = "+ \ - cvm::to_str(def_value)+" [default]\n"); \ + cvm::log("# "+std::string(key)+" = "+ \ + cvm::to_str(def_value)+" [default]\n"); \ } \ } \ \ @@ -138,11 +138,11 @@ size_t colvarparse::dummy_pos = 0; #define _get_keyval_vector_(TYPE) \ \ - bool colvarparse::get_keyval(std::string const &conf, \ - char const *key, \ - std::vector &values, \ - std::vector const &def_values, \ - Parse_Mode const parse_mode) \ + bool colvarparse::get_keyval(std::string const &conf, \ + char const *key, \ + std::vector &values, \ + std::vector const &def_values, \ + Parse_Mode const parse_mode) \ { \ std::string data; \ bool b_found = false, b_found_any = false; \ @@ -150,7 +150,7 @@ size_t colvarparse::dummy_pos = 0; \ do { \ std::string data_this = ""; \ - b_found = key_lookup(conf, key, data_this, save_pos); \ + b_found = key_lookup(conf, key, data_this, save_pos); \ if (b_found) { \ if (!b_found_any) \ b_found_any = true; \ @@ -160,11 +160,11 @@ size_t colvarparse::dummy_pos = 0; } while (b_found); \ \ if (found_count > 1) \ - cvm::log("Warning: found more than one instance of \""+ \ - std::string(key)+"\".\n"); \ + cvm::log("Warning: found more than one instance of \""+ \ + std::string(key)+"\".\n"); \ \ if (data.size()) { \ - std::istringstream is(data); \ + std::istringstream is(data); \ \ if (values.size() == 0) { \ \ @@ -172,44 +172,44 @@ size_t colvarparse::dummy_pos = 0; if (def_values.size()) \ x = def_values; \ else \ - x.assign(1, TYPE()); \ + x.assign(1, TYPE()); \ \ for (size_t i = 0; \ ( is >> x[ ((i> x) \ values[i] = x; \ else \ - cvm::error("Error: in parsing \""+ \ - std::string(key)+"\".\n", INPUT_ERROR); \ + cvm::error("Error: in parsing \""+ \ + std::string(key)+"\".\n", INPUT_ERROR); \ } \ } \ \ if (parse_mode != parse_silent) { \ - cvm::log("# "+std::string(key)+" = "+ \ - cvm::to_str(values)+"\n"); \ + cvm::log("# "+std::string(key)+" = "+ \ + cvm::to_str(values)+"\n"); \ } \ \ } else { \ \ if (b_found_any) \ - cvm::error("Error: improper or missing values for \""+ \ - std::string(key)+"\".\n", INPUT_ERROR); \ + cvm::error("Error: improper or missing values for \""+ \ + std::string(key)+"\".\n", INPUT_ERROR); \ \ for (size_t i = 0; i < values.size(); i++) \ values[i] = def_values[ (i > def_values.size()) ? 0 : i ]; \ \ if (parse_mode != parse_silent) { \ - cvm::log("# "+std::string(key)+" = "+ \ - cvm::to_str(def_values)+" [default]\n"); \ + cvm::log("# "+std::string(key)+" = "+ \ + cvm::to_str(def_values)+" [default]\n"); \ } \ } \ \ @@ -241,10 +241,10 @@ _get_keyval_vector_(colvarvalue); bool colvarparse::get_keyval(std::string const &conf, - char const *key, - bool &value, - bool const &def_value, - Parse_Mode const parse_mode) + char const *key, + bool &value, + bool const &def_value, + Parse_Mode const parse_mode) { std::string data; bool b_found = false, b_found_any = false; @@ -263,7 +263,7 @@ bool colvarparse::get_keyval(std::string const &conf, if (found_count > 1) cvm::log("Warning: found more than one instance of \""+ - std::string(key)+"\".\n"); + std::string(key)+"\".\n"); if (data.size()) { if ( (data == std::string("on")) || @@ -271,15 +271,15 @@ bool colvarparse::get_keyval(std::string const &conf, (data == std::string("true")) ) { value = true; } else if ( (data == std::string("off")) || - (data == std::string("no")) || - (data == std::string("false")) ) { + (data == std::string("no")) || + (data == std::string("false")) ) { value = false; } else cvm::fatal_error("Error: boolean values only are allowed " - "for \""+std::string(key)+"\".\n"); + "for \""+std::string(key)+"\".\n"); if (parse_mode != parse_silent) { cvm::log("# "+std::string(key)+" = "+ - (value ? "on" : "off")+"\n"); + (value ? "on" : "off")+"\n"); } } else { @@ -292,7 +292,7 @@ bool colvarparse::get_keyval(std::string const &conf, value = def_value; if (parse_mode != parse_silent) { cvm::log("# "+std::string(key)+" = "+ - (def_value ? "on" : "off")+" [default]\n"); + (def_value ? "on" : "off")+" [default]\n"); } } } @@ -351,7 +351,7 @@ int colvarparse::check_keywords(std::string &conf, char const *key) { if (cvm::debug()) cvm::log("Configuration string for \""+std::string(key)+ - "\": \"\n"+conf+"\".\n"); + "\": \"\n"+conf+"\".\n"); strip_values(conf); // after stripping, the config string has either empty lines, or @@ -383,7 +383,7 @@ int colvarparse::check_keywords(std::string &conf, char const *key) } if (!found_keyword) { cvm::log("Error: keyword \""+uk+"\" is not supported, " - "or not recognized in this context.\n"); + "or not recognized in this context.\n"); cvm::set_error_bits(INPUT_ERROR); return COLVARS_ERROR; } @@ -396,8 +396,8 @@ int colvarparse::check_keywords(std::string &conf, char const *key) std::istream & colvarparse::getline_nocomments(std::istream &is, - std::string &line, - char const delim) + std::string &line, + char const delim) { std::getline(is, line, delim); size_t const comment = line.find('#'); @@ -409,9 +409,9 @@ std::istream & colvarparse::getline_nocomments(std::istream &is, bool colvarparse::key_lookup(std::string const &conf, - char const *key_in, - std::string &data, - size_t &save_pos) + char const *key_in, + std::string &data, + size_t &save_pos) { // add this keyword to the register (in its camelCase version) add_keyword(key_in); @@ -445,7 +445,7 @@ bool colvarparse::key_lookup(std::string const &conf, if (pos > 0) { if ( std::string("\n"+white_space+ - "}").find(conf[pos-1]) == + "}").find(conf[pos-1]) == std::string::npos ) { // none of the valid delimiting characters is on the left of key b_isolated_left = false; @@ -454,7 +454,7 @@ bool colvarparse::key_lookup(std::string const &conf, if (pos < conf.size()-key.size()-1) { if ( std::string("\n"+white_space+ - "{").find(conf[pos+key.size()]) == + "{").find(conf[pos+key.size()]) == std::string::npos ) { // none of the valid delimiting characters is on the right of key b_isolated_right = false; @@ -522,9 +522,9 @@ bool colvarparse::key_lookup(std::string const &conf, // add a new line if (line_end >= conf.size()) { cvm::fatal_error("Parse error: reached the end while " - "looking for closing brace; until now " - "the following was parsed: \"\n"+ - line+"\".\n"); + "looking for closing brace; until now " + "the following was parsed: \"\n"+ + line+"\".\n"); return false; } size_t const old_end = line.size(); @@ -548,11 +548,11 @@ bool colvarparse::key_lookup(std::string const &conf, // set data_begin after the opening brace data_begin = line.find_first_of('{', data_begin) + 1; data_begin = line.find_first_not_of(white_space, - data_begin); + data_begin); // set data_end before the closing brace data_end = brace; data_end = line.find_last_not_of(white_space+"}", - data_end) + 1; + data_end) + 1; // data_end_absolute = line_end; if (data_end > line.size()) @@ -625,7 +625,7 @@ std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb) bool colvarparse::brace_check(std::string const &conf, - size_t const start_pos) + size_t const start_pos) { size_t brace_count = 0; size_t brace = start_pos; diff --git a/lib/colvars/colvarscript.cpp b/lib/colvars/colvarscript.cpp index 8bce761dd4..a550a09915 100644 --- a/lib/colvars/colvarscript.cpp +++ b/lib/colvars/colvarscript.cpp @@ -246,7 +246,7 @@ int colvarscript::proc_colvar(int argc, char const *argv[]) { std::istringstream is(f_str); is.width(cvm::cv_width); is.precision(cvm::cv_prec); - colvarvalue force(cv->type()); + colvarvalue force(cv->value()); force.is_derivative(); if (force.from_simple_string(is.str()) != COLVARS_OK) { result = "addforce : error parsing force value"; diff --git a/lib/colvars/colvartypes.h b/lib/colvars/colvartypes.h index d5c2e39341..9bed239962 100644 --- a/lib/colvars/colvartypes.h +++ b/lib/colvars/colvartypes.h @@ -237,7 +237,7 @@ public: inline size_t output_width(size_t const &real_width) const { - return real_width*this->size() + real_width*(this->size()-1) + 4; + return real_width*(this->size()) + 3*(this->size()-1) + 4; } inline friend std::istream & operator >> (std::istream &is, cvm::vector1d &v) @@ -252,8 +252,9 @@ public: return is; } size_t count = 0; - while ( (is >> v[count]) && (is >> sep) && (sep == ',') && (count < v.size()) ) { - count++; + while ( (is >> v[count]) && + (count < (v.size()-1) ? ((is >> sep) && (sep == ',')) : true) ) { + if (++count == v.size()) break; } if (count < v.size()) { is.clear(); @@ -268,6 +269,7 @@ public: std::streamsize const w = os.width(); std::streamsize const p = os.precision(); + os.width(2); os << "( "; size_t i; for (i = 0; i < v.size()-1; i++) { @@ -582,7 +584,8 @@ public: std::streamsize const w = os.width(); std::streamsize const p = os.precision(); - os << "("; + os.width(2); + os << "( "; size_t i; for (i = 0; i < m.outer_length; i++) { os << " ( "; diff --git a/lib/colvars/colvarvalue.h b/lib/colvars/colvarvalue.h index 1513e42af0..16fc5737ca 100644 --- a/lib/colvars/colvarvalue.h +++ b/lib/colvars/colvarvalue.h @@ -170,23 +170,25 @@ public: { // reset the value based on the previous type reset(); - if (value_type == type_vector) { - // delete allocated data if necessary + if ((value_type == type_vector) && (vti != type_vector)) { vector1d_value.resize(0); } value_type = vti; - // reset the value based on the new type - reset(); } /// Set the type after another \link colvarvalue \endlink inline void type(colvarvalue const &x) { - // reset the value based on the previous type + // reset the value held from based on the previous type reset(); + if (x.type() == type_vector) { + vector1d_value.resize(x.vector1d_value.size()); + } else { + if (value_type == type_vector) { + vector1d_value.resize(0); + } + } value_type = x.type(); - // reset the value based on the new type - reset(); } /// Make the type a derivative of the original type @@ -319,10 +321,10 @@ public: /// Ensure that the two types are the same within a binary operator - void static check_types(colvarvalue const &x1, colvarvalue const &x2); + int static check_types(colvarvalue const &x1, colvarvalue const &x2); /// Ensure that the two types are the same within an assignment, or that the left side is type_notset - void static check_types_assign(Type const &vt1, Type const &vt2); + int static check_types_assign(Type const &vt1, Type const &vt2); /// Undefined operation void undef_op() const; @@ -530,10 +532,10 @@ inline colvarvalue::colvarvalue(cvm::vector1d const &v, Type const &v { // reset the value based on the previous type reset(); - if (v.size() != num_dimensions(value_type)) { + if ((vti != type_vector) && (v.size() != num_dimensions(value_type))) { cvm::error("Error: trying to initialize a variable of type \""+type_desc(vti)+ - "\" using a vector of size \""+cvm::to_str(v.size())+ - "\".\n"); + "\" using a vector of size "+cvm::to_str(v.size())+ + ".\n"); value_type = type_notset; } else { value_type = vti; @@ -561,10 +563,12 @@ inline colvarvalue::colvarvalue(cvm::vector1d const &v, Type const &v } -inline void colvarvalue::check_types(colvarvalue const &x1, - colvarvalue const &x2) +inline int colvarvalue::check_types(colvarvalue const &x1, + colvarvalue const &x2) { - if (!colvarvalue::type_checking()) return; + if (!colvarvalue::type_checking()) { + return COLVARS_OK; + } if (x1.type() != x2.type()) { if (((x1.type() == type_unit3vector) && @@ -575,39 +579,49 @@ inline void colvarvalue::check_types(colvarvalue const &x1, (x2.type() == type_quaternionderiv)) || ((x2.type() == type_quaternion) && (x1.type() == type_quaternionderiv))) { - return; + return COLVARS_OK; + } else { + cvm::error("Trying to perform an operation between two colvar " + "values with different types, \""+ + colvarvalue::type_desc(x1.type())+ + "\" and \""+ + colvarvalue::type_desc(x2.type())+ + "\".\n"); + return COLVARS_ERROR; } - cvm::error("Performing an operation between two colvar " - "values with different types, \""+ - colvarvalue::type_desc(x1.type())+ - "\" and \""+ - colvarvalue::type_desc(x2.type())+ - "\".\n"); } if (x1.type() == type_vector) { if (x1.vector1d_value.size() != x2.vector1d_value.size()) { - cvm::error("Performing an operation between two vector colvar " + cvm::error("Trying to perform an operation between two vector colvar " "values with different sizes, "+ cvm::to_str(x1.vector1d_value.size())+ " and "+ cvm::to_str(x2.vector1d_value.size())+ ".\n"); + return COLVARS_ERROR; } } + return COLVARS_OK; } -inline void colvarvalue::check_types_assign(colvarvalue::Type const &vt1, - colvarvalue::Type const &vt2) +inline int colvarvalue::check_types_assign(colvarvalue::Type const &vt1, + colvarvalue::Type const &vt2) { + if (!colvarvalue::type_checking()) { + return COLVARS_OK; + } + if (vt1 != type_notset) { if (vt1 != vt2) { cvm::error("Trying to assign a colvar value with type \""+ type_desc(vt2)+"\" to one with type \""+ type_desc(vt1)+"\".\n"); + return COLVARS_ERROR; } } + return COLVARS_OK; }