Add Colvars module version 2024-03-12

CMakeLists.txt

@@ -676,6 +676,8 @@ include(gmxManageLmfit)
 
 include(gmxManageMuparser)
 
+include(gmxManageLepton)
+
 include(gmxManageColvars)
 
 ##################################################

cmake/gmxManageLepton.cmake

@@ -0,0 +1,36 @@
+# Add Lepton library, which is developed and distributed as part of OpenMM:
+# https://github.com/openmm/openmm
+
+gmx_option_multichoice(GMX_USE_LEPTON
+    "Build the Lepton library interfaced with GROMACS"
+    INTERNAL
+    INTERNAL NONE)
+mark_as_advanced(GMX_USE_LEPTON)
+
+function(gmx_manage_lepton)
+  if(GMX_USE_LEPTON STREQUAL "INTERNAL")
+
+    set(LEPTON_DIR "${CMAKE_SOURCE_DIR}/src/external/lepton")
+
+    file(GLOB LEPTON_SOURCES ${LEPTON_DIR}/src/*.cpp)
+    add_library(lepton_objlib OBJECT ${LEPTON_SOURCES})
+    target_include_directories(lepton_objlib PRIVATE ${LEPTON_DIR}/include)
+    # TODO support building as shared library
+    target_compile_options(lepton_objlib PRIVATE -DLEPTON_BUILDING_STATIC_LIBRARY)
+    set_target_properties(lepton_objlib PROPERTIES POSITION_INDEPENDENT_CODE ON)
+
+    add_library(lepton INTERFACE)
+    target_sources(lepton INTERFACE $<TARGET_OBJECTS:lepton_objlib>)
+    target_include_directories(lepton SYSTEM INTERFACE $<BUILD_INTERFACE:${LEPTON_DIR}>)
+
+    # Set flags so that Colvars can leverage Lepton functionality
+    # TODO handle the case when Lepton is built without Colvars?
+    target_include_directories(colvars_objlib PRIVATE ${LEPTON_DIR}/include)
+    target_compile_options(colvars_objlib PRIVATE -DLEPTON -DLEPTON_USE_STATIC_LIBRARIES)
+
+  else()
+
+    # Dummy target
+    add_library(lepton INTERFACE)
+  endif()
+endfunction()

src/external/colvars/colvar.cpp

@@ -22,10 +22,8 @@
 #include "colvars_memstream.h"
 
 
-
-std::map<std::string, std::function<colvar::cvc *(const std::string &conf)>>
-    colvar::global_cvc_map =
-        std::map<std::string, std::function<colvar::cvc *(const std::string &conf)>>();
+std::map<std::string, std::function<colvar::cvc *()>> colvar::global_cvc_map =
+    std::map<std::string, std::function<colvar::cvc *()>>();
 
 std::map<std::string, std::string> colvar::global_cvc_desc_map =
     std::map<std::string, std::string>();
@@ -37,6 +35,7 @@ colvar::colvar()
   after_restart = false;
   kinetic_energy = 0.0;
   potential_energy = 0.0;
+  period = 0.0;
 
 #ifdef LEPTON
   dev_null = 0.0;
@@ -139,7 +138,14 @@ int colvar::init(std::string const &conf)
 
     // Sort array of cvcs based on their names
     // Note: default CVC names are in input order for same type of CVC
-    std::sort(cvcs.begin(), cvcs.end(), colvar::compare_cvc);
+    std::sort(cvcs.begin(), cvcs.end(),
+              [](std::shared_ptr<colvar::cvc> const &cvc1,
+                 std::shared_ptr<colvar::cvc> const &cvc2) -> bool {
+                if (cvc1 && cvc2) {
+                  return cvc1->name < cvc2->name;
+                }
+                return false;
+              });
 
     if(cvcs.size() > 1) {
       cvm::log("Sorted list of components for this scripted colvar:\n");
@@ -194,9 +200,9 @@ int colvar::init(std::string const &conf)
 
         if ((cvcs[i])->sup_np < 0) {
           cvm::log("Warning: you chose a negative exponent in the combination; "
-                    "if you apply forces, the simulation may become unstable "
-                    "when the component \""+
-                    (cvcs[i])->function_type+"\" approaches zero.\n");
+                   "if you apply forces, the simulation may become unstable "
+                   "when the component \""+
+                   (cvcs[i])->function_type()+"\" approaches zero.\n");
         }
       }
     }
@@ -303,7 +309,7 @@ int colvar::init(std::string const &conf)
   error_code |= init_grid_parameters(conf);
 
   // Detect if we have a single component that is an alchemical lambda
-  if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type == "alchLambda") {
+  if (is_enabled(f_cv_single_cvc) && cvcs[0]->function_type() == "alchLambda") {
     enable(f_cv_external);
   }
 
@@ -755,8 +761,8 @@ template <typename def_class_name>
 void colvar::add_component_type(char const *def_description, char const *def_config_key)
 {
   if (global_cvc_map.count(def_config_key) == 0) {
-    global_cvc_map[def_config_key] = [](const std::string &cvc_conf) {
-      return new def_class_name(cvc_conf);
+    global_cvc_map[def_config_key] = []() {
+      return new def_class_name();
     };
     global_cvc_desc_map[def_config_key] = std::string(def_description);
   }
@@ -767,48 +773,39 @@ int colvar::init_components_type(const std::string& conf, const char* def_config
   size_t def_count = 0;
   std::string def_conf = "";
   size_t pos = 0;
+  int error_code = COLVARS_OK;
   while ( this->key_lookup(conf,
                            def_config_key,
                            &def_conf,
                            &pos) ) {
-    if (!def_conf.size()) continue;
+
     cvm::log("Initializing "
              "a new \""+std::string(def_config_key)+"\" component"+
              (cvm::debug() ? ", with configuration:\n"+def_conf
               : ".\n"));
-    cvc *cvcp = global_cvc_map[def_config_key](def_conf);
-    cvm::increase_depth();
-    if (cvcp) {
-      int error_code = cvcp->init_code;
-      cvcs.push_back(cvcp);
-      error_code |= cvcp->set_function_type(def_config_key);
-      if (error_code == COLVARS_OK) {
-        error_code |= cvcp->check_keywords(def_conf, def_config_key);
-      }
-      if (error_code != COLVARS_OK) {
-        cvm::decrease_depth();
-        return cvm::error("Error: in setting up component \"" + std::string(def_config_key) +
-                              "\".\n",
-                          COLVARS_INPUT_ERROR);
-      }
-    } else {
-      cvm::decrease_depth();
-      return cvm::error("Error: in allocating component \"" + std::string(def_config_key) + "\".\n",
+    cvc *cvcp = global_cvc_map[def_config_key]();
+    if (!cvcp) {
+      return cvm::error("Error: in creating object of type \"" + std::string(def_config_key) +
+                            "\".\n",
                         COLVARS_MEMORY_ERROR);
     }
+    cvcs.push_back(std::shared_ptr<colvar::cvc>(cvcp));
 
-    if ((cvcp->period != 0.0) || (cvcp->wrap_center != 0.0)) {
-      if (!cvcp->is_enabled(f_cvc_periodic)) {
-        cvm::decrease_depth();
-        return cvm::error("Error: invalid use of period and/or "
-                          "wrapAround in a \"" +
-                              std::string(def_config_key) + "\" component.\n" +
-                              "Period: " + cvm::to_str(cvcp->period) +
-                              " wrapAround: " + cvm::to_str(cvcp->wrap_center),
-                          COLVARS_INPUT_ERROR);
-      }
+    cvm::increase_depth();
+    int error_code_this = cvcp->init(def_conf);
+    if (error_code_this == COLVARS_OK) {
+      // Checking for invalid keywords only if the parsing was successful, otherwise any
+      // early-returns due to errors would raise false positives
+      error_code_this |= cvcp->check_keywords(def_conf, def_config_key);
+    }
+    cvm::decrease_depth();
+    if (error_code_this != COLVARS_OK) {
+      error_code |=
+          cvm::error("Error: in setting up component \"" + std::string(def_config_key) + "\".\n",
+                     COLVARS_INPUT_ERROR);
     }
 
+    // Set default name if it doesn't have one
     if ( ! cvcs.back()->name.size()) {
       std::ostringstream s;
       s << def_config_key << std::setfill('0') << std::setw(4) << ++def_count;
@@ -822,15 +819,13 @@ int colvar::init_components_type(const std::string& conf, const char* def_config
                (cvm::debug() ? ", named \"" + cvcs.back()->name + "\"" : "") + ".\n");
     }
 
-    cvm::decrease_depth();
-
     def_conf = "";
     if (cvm::debug()) {
       cvm::log("Parsed " + cvm::to_str(cvcs.size()) + " components at this time.\n");
     }
   }
 
-  return COLVARS_OK;
+  return error_code;
 }
 
 
@@ -944,7 +939,7 @@ int colvar::init_components(std::string const &conf)
   if (error_code == COLVARS_OK) {
     // Store list of children cvcs for dependency checking purposes
     for (i = 0; i < cvcs.size(); i++) {
-      add_child(cvcs[i]);
+      add_child(cvcs[i].get());
     }
     // By default all CVCs are active at the start
     n_active_cvcs = cvcs.size();
@@ -1221,7 +1216,7 @@ int colvar::init_dependencies() {
 
   // Initialize feature_states for each instance
   feature_states.reserve(f_cv_ntot);
-  for (i = 0; i < f_cv_ntot; i++) {
+  for (i = feature_states.size(); i < f_cv_ntot; i++) {
     feature_states.push_back(feature_state(true, false));
     // Most features are available, so we set them so
     // and list exceptions below
@@ -1284,14 +1279,10 @@ colvar::~colvar()
   // for dependency purposes
   remove_all_children();
 
-  for (std::vector<cvc *>::reverse_iterator ci = cvcs.rbegin();
-      ci != cvcs.rend();
-      ++ci) {
-    // clear all children of this cvc (i.e. its atom groups)
-    // because the cvc base class destructor can't do it early enough
-    // and we don't want to have each cvc derived class do it separately
+  for (auto ci = cvcs.rbegin(); ci != cvcs.rend(); ++ci) {
+    // Clear all children of this cvc (i.e. its atom groups), because the cvc base class destructor
+    // can't do it early enough and we don't want to have each cvc derived class do it separately
     (*ci)->remove_all_children();
-    delete *ci;
   }
   cvcs.clear();
 
@@ -1513,6 +1504,7 @@ int colvar::collect_cvc_values()
               cvm::to_str(x, cvm::cv_width, cvm::cv_prec)+".\n");
 
   if (after_restart) {
+    x_old = x_restart;
     if (cvm::proxy->simulation_running()) {
       cvm::real const jump2 = dist2(x, x_restart) / (width*width);
       if (jump2 > 0.25) {
@@ -1707,12 +1699,13 @@ int colvar::calc_colvar_properties()
     // Do the same if no simulation is running (eg. VMD postprocessing)
     if ((cvm::step_relative() == 0 && !after_restart) || x_ext.type() == colvarvalue::type_notset || !cvm::proxy->simulation_running()) {
       x_ext = x;
+      cvm::log("Initializing extended coordinate to colvar value.\n");
       if (is_enabled(f_cv_reflecting_lower_boundary) && x_ext < lower_boundary) {
-        cvm::log("Warning: initializing extended coordinate to reflective lower boundary, as colvar value is below.");
+        cvm::log("Warning: initializing extended coordinate to reflective lower boundary, as colvar value is below.\n");
         x_ext = lower_boundary;
       }
       if (is_enabled(f_cv_reflecting_upper_boundary) && x_ext > upper_boundary) {
-        cvm::log("Warning: initializing extended coordinate to reflective upper boundary, as colvar value is above.");
+        cvm::log("Warning: initializing extended coordinate to reflective upper boundary, as colvar value is above.\n");
         x_ext = upper_boundary;
       }
 
@@ -1722,8 +1715,18 @@ int colvar::calc_colvar_properties()
     // Special case of a repeated timestep (eg. multiple NAMD "run" statements)
     // revert values of the extended coordinate and velocity prior to latest integration
     if (cvm::proxy->simulation_running() && cvm::step_relative() == prev_timestep) {
-      x_ext = prev_x_ext;
-      v_ext = prev_v_ext;
+      // Detect jumps due to discrete changes in coordinates (eg. in replica exchange schemes)
+      cvm::real const jump2 = dist2(x, x_old) / (width*width);
+      if (jump2 > 0.25) {
+        cvm::log("Detected discrete jump in colvar value from "
+            + cvm::to_str(x_old) + " to " + cvm::to_str(x) + ".\n");
+        cvm::log("Reinitializing extended coordinate to colvar value.\n");
+        x_ext = x;
+      } else {
+        cvm::log("Reinitializing extended coordinate to last value.\n");
+        x_ext = prev_x_ext;
+        v_ext = prev_v_ext;
+      }
     }
     // report the restraint center as "value"
     // These position and velocities come from integration at the _previous timestep_ in update_forces_energy()
@@ -1939,9 +1942,8 @@ int colvar::end_of_step()
   if (cvm::debug())
     cvm::log("End of step for colvar \""+this->name+"\".\n");
 
-  if (is_enabled(f_cv_fdiff_velocity)) {
-    x_old = x;
-  }
+  // Used for fdiff_velocity and for detecting jumps for extended Lagrangian colvars
+  x_old = x;
 
   if (is_enabled(f_cv_subtract_applied_force)) {
     f_old = f;

src/external/colvars/colvar.h

@@ -14,6 +14,7 @@
 #include <list>
 #include <iosfwd>
 #include <map>
+#include <memory>
 
 #include "colvarmodule.h"
 #include "colvarvalue.h"
@@ -617,7 +618,6 @@ class colvar : public colvarparse, public colvardeps {
   class dihedPC;
   class alch_lambda;
   class alch_Flambda;
-  class componentDisabled;
   class CartesianBasedPath;
   class aspath;
   class azpath;
@@ -645,8 +645,7 @@ class colvar : public colvarparse, public colvardeps {
   class map_total;
 
   /// A global mapping of cvc names to the cvc constructors
-  static const std::map<std::string, std::function<colvar::cvc *(const std::string &subcv_conf)>> &
-  get_global_cvc_map()
+  static const std::map<std::string, std::function<colvar::cvc *()>> &get_global_cvc_map()
   {
     return global_cvc_map;
   }
@@ -656,8 +655,8 @@ class colvar : public colvarparse, public colvardeps {
 
 protected:
 
-  /// \brief Array of \link colvar::cvc \endlink objects
-  std::vector<cvc *> cvcs;
+  /// Array of components objects
+  std::vector<std::shared_ptr<colvar::cvc>> cvcs;
 
   /// \brief Flags to enable or disable cvcs at next colvar evaluation
   std::vector<bool> cvc_flags;
@@ -689,8 +688,7 @@ class colvar : public colvarparse, public colvardeps {
 #endif
 
   /// A global mapping of cvc names to the cvc constructors
-  static std::map<std::string, std::function<colvar::cvc *(const std::string &conf)>>
-      global_cvc_map;
+  static std::map<std::string, std::function<colvar::cvc *()>> global_cvc_map;
 
   /// A global mapping of cvc names to the corresponding descriptions
   static std::map<std::string, std::string> global_cvc_desc_map;

src/external/colvars/colvaratoms.cpp

@@ -286,7 +286,7 @@ int cvm::atom_group::init_dependencies() {
   // Initialize feature_states for each instance
   // default as unavailable, not enabled
   feature_states.reserve(f_ag_ntot);
-  for (i = 0; i < colvardeps::f_ag_ntot; i++) {
+  for (i = feature_states.size(); i < colvardeps::f_ag_ntot; i++) {
     feature_states.push_back(feature_state(false, false));
   }
 
@@ -532,8 +532,8 @@ int cvm::atom_group::parse(std::string const &group_conf)
         cvm::error("Error: atomsColValue, if provided, must be non-zero.\n", COLVARS_INPUT_ERROR);
       }
 
-      // NOTE: calls to add_atom() and/or add_atom_id() are in the proxy-implemented function
-      error_code |= cvm::load_atoms(atoms_file_name.c_str(), *this, atoms_col, atoms_col_value);
+      error_code |= cvm::main()->proxy->load_atoms_pdb(atoms_file_name.c_str(), *this, atoms_col,
+                                                       atoms_col_value);
     }
   }