diff --git a/cmake/CMakeLists.txt b/cmake/CMakeLists.txt index 13e1d1539f..44c864941c 100644 --- a/cmake/CMakeLists.txt +++ b/cmake/CMakeLists.txt @@ -2,7 +2,7 @@ # CMake build system # This file is part of LAMMPS # Created by Christoph Junghans and Richard Berger -cmake_minimum_required(VERSION 2.8.12) +cmake_minimum_required(VERSION 3.10) project(lammps CXX) set(SOVERSION 0) @@ -55,12 +55,9 @@ if(${CMAKE_CXX_COMPILER_ID} STREQUAL "Intel") set (CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} -restrict -std=c++11") endif() -option(DISABLE_CXX11_REQUIREMENT "Disable check that requires C++11 for compiling LAMMPS" OFF) -if(DISABLE_CXX11_REQUIREMENT) - add_definitions(-DLAMMPS_CXX98) -# else() -# set(CMAKE_CXX_STANDARD 11) -endif() +# we require C++11 +set(CMAKE_CXX_STANDARD 11) +set(CMAKE_CXX_STANDARD_REQUIRED ON) # GNU compiler features if (${CMAKE_CXX_COMPILER_ID} STREQUAL "GNU") @@ -355,9 +352,6 @@ endforeach(HEADER) set(MATH_LIBRARIES "m" CACHE STRING "math library") mark_as_advanced( MATH_LIBRARIES ) include(CheckLibraryExists) -if (CMAKE_VERSION VERSION_LESS "3.4") - enable_language(C) # check_library_exists isn't supported without a C compiler before v3.4 -endif() # RB: disabled this check because it breaks with KOKKOS CUDA enabled #foreach(FUNC sin cos) # check_library_exists(${MATH_LIBRARIES} ${FUNC} "" FOUND_${FUNC}_${MATH_LIBRARIES}) diff --git a/cmake/Modules/Packages/GPU.cmake b/cmake/Modules/Packages/GPU.cmake index dab9d51a3f..200d8fb51e 100644 --- a/cmake/Modules/Packages/GPU.cmake +++ b/cmake/Modules/Packages/GPU.cmake @@ -1,7 +1,4 @@ if(PKG_GPU) - if (CMAKE_VERSION VERSION_LESS "3.1") - message(FATAL_ERROR "For the GPU package you need at least cmake-3.1") - endif() set(GPU_SOURCES_DIR ${LAMMPS_SOURCE_DIR}/GPU) set(GPU_SOURCES ${GPU_SOURCES_DIR}/gpu_extra.h ${GPU_SOURCES_DIR}/fix_gpu.h diff --git a/cmake/Modules/Packages/LATTE.cmake b/cmake/Modules/Packages/LATTE.cmake index de7116780b..8bcda84cdc 100644 --- a/cmake/Modules/Packages/LATTE.cmake +++ b/cmake/Modules/Packages/LATTE.cmake @@ -8,9 +8,6 @@ if(PKG_LATTE) endif() option(DOWNLOAD_LATTE "Download the LATTE library instead of using an already installed one" ${DOWNLOAD_LATTE_DEFAULT}) if(DOWNLOAD_LATTE) - if (CMAKE_VERSION VERSION_LESS "3.7") # due to SOURCE_SUBDIR - message(FATAL_ERROR "For downlading LATTE you need at least cmake-3.7") - endif() if(CMAKE_GENERATOR STREQUAL "Ninja") message(FATAL_ERROR "Cannot build downloaded LATTE library with Ninja build tool") endif() diff --git a/cmake/Modules/Packages/MSCG.cmake b/cmake/Modules/Packages/MSCG.cmake index b442580583..99d98659ee 100644 --- a/cmake/Modules/Packages/MSCG.cmake +++ b/cmake/Modules/Packages/MSCG.cmake @@ -8,9 +8,6 @@ if(PKG_MSCG) endif() option(DOWNLOAD_MSCG "Download MSCG library instead of using an already installed one)" ${DOWNLOAD_MSCG_DEFAULT}) if(DOWNLOAD_MSCG) - if (CMAKE_VERSION VERSION_LESS "3.7") # due to SOURCE_SUBDIR - message(FATAL_ERROR "For downlading MSCG you need at least cmake-3.7") - endif() if(CMAKE_GENERATOR STREQUAL "Ninja") message(FATAL_ERROR "Cannot build downloaded MSCG library with Ninja build tool") endif() diff --git a/cmake/Modules/Packages/USER-MOLFILE.cmake b/cmake/Modules/Packages/USER-MOLFILE.cmake index b8c4234d26..16ffc34994 100644 --- a/cmake/Modules/Packages/USER-MOLFILE.cmake +++ b/cmake/Modules/Packages/USER-MOLFILE.cmake @@ -1,8 +1,4 @@ if(PKG_USER-MOLFILE) - if (CMAKE_VERSION VERSION_LESS "3.10") # due to INTERFACE without a library - message(FATAL_ERROR "For configuring USER-MOLFILE you need CMake 3.10 or later") - endif() - set(MOLFILE_INCLUDE_DIRS "${LAMMPS_LIB_SOURCE_DIR}/molfile" CACHE STRING "Path to VMD molfile plugin headers") add_library(molfile INTERFACE) target_include_directories(molfile INTERFACE ${MOLFILE_INCLUDE_DIRS}) diff --git a/doc/src/Build_cmake.rst b/doc/src/Build_cmake.rst index 8314afaa0e..939431c6b3 100644 --- a/doc/src/Build_cmake.rst +++ b/doc/src/Build_cmake.rst @@ -95,10 +95,8 @@ this directory or sub-directories within it that CMake creates. directory to un-install all packages. The purge removes all the \*.h files auto-generated by make. -You must have CMake version 2.8 or later on your system to build -LAMMPS. A handful of LAMMPS packages (KOKKOS, LATTE, MSCG) require a -later version. CMake will print a message telling you if a later -version is required. Installation instructions for CMake are below. +You must have CMake version 3.10 or later on your system to build +LAMMPS. Installation instructions for CMake are below. After the initial build, if you edit LAMMPS source files, or add your own new files to the source directory, you can just re-type make from diff --git a/doc/src/Build_settings.rst b/doc/src/Build_settings.rst index 1e668ee95d..5f835f76a0 100644 --- a/doc/src/Build_settings.rst +++ b/doc/src/Build_settings.rst @@ -23,11 +23,11 @@ explain how to do this for building both with CMake and make. C++11 standard compliance ------------------------------------------ -The LAMMPS developers plan to transition to make the C++11 standard the -minimum requirement for compiling LAMMPS. Currently this only applies to -some packages like KOKKOS while the rest aims to be compatible with the C++98 -standard. Most currently used compilers are compatible with C++11; some need -to set extra flags to enable C++11 compliance. Example for GNU c++: +A C++11 standard compatible compiler is a requirement for compiling LAMMPS. +LAMMPS version 3 March 2020 is the last version compatible with the previous +C++98 standard for the core code and most packages. Most currently used +C++ compilers are compatible with C++11, but some older ones may need extra +flags to enable C++11 compliance. Example for GNU c++ 4.8.x: .. code-block:: make diff --git a/doc/src/PDF/colvars-refman-lammps.pdf b/doc/src/PDF/colvars-refman-lammps.pdf index dd028ae65e..c9a67d6a1d 100644 Binary files a/doc/src/PDF/colvars-refman-lammps.pdf and b/doc/src/PDF/colvars-refman-lammps.pdf differ diff --git a/lib/colvars/Makefile.deps b/lib/colvars/Makefile.deps index 33a9dd10ed..2c0dabea6a 100644 --- a/lib/colvars/Makefile.deps +++ b/lib/colvars/Makefile.deps @@ -27,46 +27,58 @@ $(COLVARS_OBJ_DIR)colvarbias_restraint.o: colvarbias_restraint.cpp \ $(COLVARS_OBJ_DIR)colvarcomp_angles.o: colvarcomp_angles.cpp \ colvarmodule.h colvars_version.h colvar.h colvarvalue.h colvartypes.h \ colvarparse.h colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvar_arithmeticpath.h -$(COLVARS_OBJ_DIR)colvarcomp_apath.o: colvarcomp_apath.cpp + colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h +$(COLVARS_OBJ_DIR)colvarcomp_apath.o: colvarcomp_apath.cpp colvarmodule.h \ + colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ + colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_coordnums.o: colvarcomp_coordnums.cpp \ colvarmodule.h colvars_version.h colvarparse.h colvarvalue.h \ colvartypes.h colvarparams.h colvaratoms.h colvarproxy.h colvardeps.h \ - colvar.h colvarcomp.h colvar_arithmeticpath.h + colvar.h colvarcomp.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp.o: colvarcomp.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvar.h colvarparse.h \ colvarparams.h colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \ - colvar_arithmeticpath.h + colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_distances.o: colvarcomp_distances.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ - colvaratoms.h colvarproxy.h colvar_arithmeticpath.h -$(COLVARS_OBJ_DIR)colvarcomp_gpath.o: colvarcomp_gpath.cpp + colvaratoms.h colvarproxy.h colvar_arithmeticpath.h \ + colvar_geometricpath.h +$(COLVARS_OBJ_DIR)colvarcomp_gpath.o: colvarcomp_gpath.cpp colvarmodule.h \ + colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ + colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ + colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_protein.o: colvarcomp_protein.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ - colvaratoms.h colvarproxy.h colvar_arithmeticpath.h + colvaratoms.h colvarproxy.h colvar_arithmeticpath.h \ + colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarcomp_rotations.o: colvarcomp_rotations.cpp \ colvarmodule.h colvars_version.h colvarvalue.h colvartypes.h \ colvarparse.h colvarparams.h colvar.h colvardeps.h colvarcomp.h \ - colvaratoms.h colvarproxy.h colvar_arithmeticpath.h + colvaratoms.h colvarproxy.h colvar_arithmeticpath.h \ + colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvar.o: colvar.cpp colvarmodule.h colvars_version.h \ colvarvalue.h colvartypes.h colvarparse.h colvarparams.h colvar.h \ colvardeps.h colvarcomp.h colvaratoms.h colvarproxy.h \ - colvar_arithmeticpath.h colvarscript.h colvarbias.h + colvar_arithmeticpath.h colvar_geometricpath.h colvarscript.h \ + colvarbias.h $(COLVARS_OBJ_DIR)colvardeps.o: colvardeps.cpp colvarmodule.h \ colvars_version.h colvarproxy.h colvartypes.h colvarvalue.h colvardeps.h \ colvarparse.h colvarparams.h $(COLVARS_OBJ_DIR)colvargrid.o: colvargrid.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvarparse.h \ colvarparams.h colvar.h colvardeps.h colvarcomp.h colvaratoms.h \ - colvarproxy.h colvar_arithmeticpath.h colvargrid.h + colvarproxy.h colvar_arithmeticpath.h colvar_geometricpath.h \ + colvargrid.h $(COLVARS_OBJ_DIR)colvarmodule.o: colvarmodule.cpp colvarmodule.h \ colvars_version.h colvarparse.h colvarvalue.h colvartypes.h \ colvarparams.h colvarproxy.h colvar.h colvardeps.h colvarbias.h \ colvarbias_abf.h colvargrid.h colvar_UIestimator.h colvarbias_alb.h \ colvarbias_histogram.h colvarbias_meta.h colvarbias_restraint.h \ - colvarscript.h colvaratoms.h colvarcomp.h colvar_arithmeticpath.h + colvarscript.h colvaratoms.h colvarcomp.h colvar_arithmeticpath.h \ + colvar_geometricpath.h $(COLVARS_OBJ_DIR)colvarparams.o: colvarparams.cpp colvarmodule.h \ colvars_version.h colvarvalue.h colvartypes.h colvarparams.h $(COLVARS_OBJ_DIR)colvarparse.o: colvarparse.cpp colvarmodule.h \ diff --git a/lib/colvars/colvar.cpp b/lib/colvars/colvar.cpp index b853aaa1b8..fb95bda5ec 100644 --- a/lib/colvars/colvar.cpp +++ b/lib/colvars/colvar.cpp @@ -33,6 +33,8 @@ colvar::colvar() dev_null = 0.0; #endif + expand_boundaries = false; + description = "uninitialized colvar"; init_dependencies(); } @@ -474,6 +476,8 @@ int colvar::init_custom_function(std::string const &conf) int colvar::init_grid_parameters(std::string const &conf) { + int error_code = COLVARS_OK; + colvarmodule *cv = cvm::main(); cvm::real default_width = width; @@ -527,15 +531,18 @@ int colvar::init_grid_parameters(std::string const &conf) disable(f_cv_hard_upper_boundary); } + // Parse legacy wall options and set up a harmonicWalls bias if needed + cvm::real lower_wall_k = 0.0, upper_wall_k = 0.0; + cvm::real lower_wall = 0.0, upper_wall = 0.0; std::string lw_conf, uw_conf; + if (get_keyval(conf, "lowerWallConstant", lower_wall_k, 0.0, parse_silent)) { cvm::log("Reading legacy options lowerWall and lowerWallConstant: " - "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n"); - lower_wall.type(value()); - if (!get_keyval(conf, "lowerWall", lower_wall, lower_boundary)) { - cvm::log("Warning: lowerWall will need to be " - "defined explicitly in the next release.\n"); + "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n"); + if (!get_keyval(conf, "lowerWall", lower_wall)) { + error_code != cvm::error("Error: the value of lowerWall must be set " + "explicitly.\n", INPUT_ERROR); } lw_conf = std::string("\n\ lowerWallConstant "+cvm::to_str(lower_wall_k*width*width)+"\n\ @@ -545,11 +552,10 @@ int colvar::init_grid_parameters(std::string const &conf) if (get_keyval(conf, "upperWallConstant", upper_wall_k, 0.0, parse_silent)) { cvm::log("Reading legacy options upperWall and upperWallConstant: " - "consider using a harmonicWalls restraint\n(caution: force constant would then be scaled by width^2).\n"); - upper_wall.type(value()); - if (!get_keyval(conf, "upperWall", upper_wall, upper_boundary)) { - cvm::log("Warning: upperWall will need to be " - "defined explicitly in the next release.\n"); + "consider using a harmonicWalls restraint (caution: force constant would then be scaled by width^2).\n"); + if (!get_keyval(conf, "upperWall", upper_wall)) { + error_code != cvm::error("Error: the value of upperWall must be set " + "explicitly.\n", INPUT_ERROR); } uw_conf = std::string("\n\ upperWallConstant "+cvm::to_str(upper_wall_k*width*width)+"\n\ @@ -558,12 +564,11 @@ int colvar::init_grid_parameters(std::string const &conf) if (lw_conf.size() && uw_conf.size()) { if (lower_wall >= upper_wall) { - cvm::error("Error: the upper wall, "+ - cvm::to_str(upper_wall)+ - ", is not higher than the lower wall, "+ - cvm::to_str(lower_wall)+".\n", - INPUT_ERROR); - return INPUT_ERROR; + error_code |= cvm::error("Error: the upper wall, "+ + cvm::to_str(upper_wall)+ + ", is not higher than the lower wall, "+ + cvm::to_str(lower_wall)+".\n", + INPUT_ERROR); } } @@ -575,7 +580,7 @@ harmonicWalls {\n\ colvars "+this->name+"\n"+lw_conf+uw_conf+"\ timeStepFactor "+cvm::to_str(time_step_factor)+"\n"+ "}\n"); - cv->append_new_config(walls_conf); + error_code |= cv->append_new_config(walls_conf); } } @@ -588,29 +593,27 @@ harmonicWalls {\n\ // consistency checks for boundaries and walls if (is_enabled(f_cv_lower_boundary) && is_enabled(f_cv_upper_boundary)) { if (lower_boundary >= upper_boundary) { - cvm::error("Error: the upper boundary, "+ - cvm::to_str(upper_boundary)+ - ", is not higher than the lower boundary, "+ - cvm::to_str(lower_boundary)+".\n", - INPUT_ERROR); - return INPUT_ERROR; + error_code |= cvm::error("Error: the upper boundary, "+ + cvm::to_str(upper_boundary)+ + ", is not higher than the lower boundary, "+ + cvm::to_str(lower_boundary)+".\n", + INPUT_ERROR); } } - get_keyval(conf, "expandBoundaries", expand_boundaries, false); + get_keyval(conf, "expandBoundaries", expand_boundaries, expand_boundaries); if (expand_boundaries && periodic_boundaries()) { - cvm::error("Error: trying to expand boundaries that already " - "cover a whole period of a periodic colvar.\n", - INPUT_ERROR); - return INPUT_ERROR; + error_code |= cvm::error("Error: trying to expand boundaries that already " + "cover a whole period of a periodic colvar.\n", + INPUT_ERROR); } + if (expand_boundaries && is_enabled(f_cv_hard_lower_boundary) && is_enabled(f_cv_hard_upper_boundary)) { - cvm::error("Error: inconsistent configuration " - "(trying to expand boundaries with both " - "hardLowerBoundary and hardUpperBoundary enabled).\n", - INPUT_ERROR); - return INPUT_ERROR; + error_code |= cvm::error("Error: inconsistent configuration " + "(trying to expand boundaries, but both " + "hardLowerBoundary and hardUpperBoundary " + "are enabled).\n", INPUT_ERROR); } return COLVARS_OK; @@ -2077,12 +2080,12 @@ void colvar::wrap(colvarvalue &x_unwrapped) const // ******************** INPUT FUNCTIONS ******************** -std::istream & colvar::read_restart(std::istream &is) +std::istream & colvar::read_state(std::istream &is) { size_t const start_pos = is.tellg(); std::string conf; - if ( !(is >> colvarparse::read_block("colvar", conf)) ) { + if ( !(is >> colvarparse::read_block("colvar", &conf)) ) { // this is not a colvar block is.clear(); is.seekg(start_pos, std::ios::beg); @@ -2092,15 +2095,24 @@ std::istream & colvar::read_restart(std::istream &is) { std::string check_name = ""; - if ( (get_keyval(conf, "name", check_name, - std::string(""), colvarparse::parse_silent)) && - (check_name != name) ) { - cvm::error("Error: the state file does not match the " - "configuration file, at colvar \""+name+"\".\n"); - } + get_keyval(conf, "name", check_name, + std::string(""), colvarparse::parse_silent); if (check_name.size() == 0) { cvm::error("Error: Collective variable in the " - "restart file without any identifier.\n"); + "restart file without any identifier.\n", INPUT_ERROR); + is.clear(); + is.seekg(start_pos, std::ios::beg); + is.setstate(std::ios::failbit); + return is; + } + + if (check_name != name) { + if (cvm::debug()) { + cvm::log("Ignoring state of colvar \""+check_name+ + "\": this colvar is named \""+name+"\".\n"); + } + is.seekg(start_pos, std::ios::beg); + return is; } } @@ -2199,7 +2211,7 @@ std::istream & colvar::read_traj(std::istream &is) // ******************** OUTPUT FUNCTIONS ******************** -std::ostream & colvar::write_restart(std::ostream &os) { +std::ostream & colvar::write_state(std::ostream &os) { os << "colvar {\n" << " name " << name << "\n" diff --git a/lib/colvars/colvar.h b/lib/colvars/colvar.h index f7d5a55627..e6fafbdec3 100644 --- a/lib/colvars/colvar.h +++ b/lib/colvars/colvar.h @@ -229,17 +229,9 @@ public: /// \brief Location of the lower boundary colvarvalue lower_boundary; - /// \brief Location of the lower wall - colvarvalue lower_wall; - /// \brief Force constant for the lower boundary potential (|x-xb|^2) - cvm::real lower_wall_k; /// \brief Location of the upper boundary colvarvalue upper_boundary; - /// \brief Location of the upper wall - colvarvalue upper_wall; - /// \brief Force constant for the upper boundary potential (|x-xb|^2) - cvm::real upper_wall_k; /// \brief Is the interval defined by the two boundaries periodic? bool periodic_boundaries() const; @@ -453,9 +445,9 @@ public: std::ostream & write_traj_label(std::ostream &os); /// Read the collective variable from a restart file - std::istream & read_restart(std::istream &is); + std::istream & read_state(std::istream &is); /// Write the collective variable to a restart file - std::ostream & write_restart(std::ostream &os); + std::ostream & write_state(std::ostream &os); /// Write output files (if defined, e.g. in analysis mode) int write_output_files(); diff --git a/lib/colvars/colvarbias.cpp b/lib/colvars/colvarbias.cpp index ff7bbad19b..09bcca01b5 100644 --- a/lib/colvars/colvarbias.cpp +++ b/lib/colvars/colvarbias.cpp @@ -15,8 +15,10 @@ colvarbias::colvarbias(char const *key) - : bias_type(to_lower_cppstr(key)) { + bias_type = to_lower_cppstr(key); + state_keyword = bias_type; + description = "uninitialized " + cvm::to_str(key) + " bias"; init_dependencies(); rank = 1; @@ -25,6 +27,7 @@ colvarbias::colvarbias(char const *key) b_output_energy = false; reset(); state_file_step = 0L; + matching_state = false; } @@ -80,10 +83,21 @@ int colvarbias::init(std::string const &conf) cvm::log("Reinitializing bias \""+name+"\".\n"); } + colvar_values.resize(num_variables()); + + for (i = 0; i < num_variables(); i++) { + colvar_values[i].type(colvars[i]->value().type()); + colvar_forces[i].type(colvar_values[i].type()); + previous_colvar_forces[i].type(colvar_values[i].type()); + } + output_prefix = cvm::output_prefix(); get_keyval(conf, "outputEnergy", b_output_energy, b_output_energy); + // Disabled by default in base class; default value can be overridden by derived class constructor + get_keyval_feature(this, conf, "bypassExtendedLagrangian", f_cvb_bypass_ext_lagrangian, is_enabled(f_cvb_bypass_ext_lagrangian), parse_silent); + get_keyval(conf, "timeStepFactor", time_step_factor, 1); if (time_step_factor < 1) { cvm::error("Error: timeStepFactor must be 1 or greater.\n"); @@ -114,6 +128,11 @@ int colvarbias::init_dependencies() { init_feature(f_cvb_apply_force, "apply force", f_type_user); require_feature_children(f_cvb_apply_force, f_cv_gradient); + init_feature(f_cvb_bypass_ext_lagrangian, "bypass extended-Lagrangian coordinates", f_type_user); + // The exclusion below prevents the inconsistency where biasing forces are applied onto + // the actual colvar, while total forces are measured on the extended coordinate + exclude_feature_self(f_cvb_bypass_ext_lagrangian, f_cvb_get_total_force); + init_feature(f_cvb_get_total_force, "obtain total force", f_type_dynamic); require_feature_children(f_cvb_get_total_force, f_cv_total_force); @@ -158,6 +177,12 @@ int colvarbias::init_dependencies() { // only compute TI samples when deriving from colvarbias_ti feature_states[f_cvb_calc_ti_samples].available = false; + // The feature f_cvb_bypass_ext_lagrangian is only implemented by some derived classes + // (initially, harmonicWalls) + feature_states[f_cvb_bypass_ext_lagrangian].available = false; + // disabled by default; can be changed by derived classes that implement it + feature_states[f_cvb_bypass_ext_lagrangian].enabled = false; + return COLVARS_OK; } @@ -265,6 +290,11 @@ int colvarbias::update() has_data = true; + // Update the cached colvar values + for (size_t i = 0; i < num_variables(); i++) { + colvar_values[i] = colvars[i]->value(); + } + error_code |= calc_energy(NULL); error_code |= calc_forces(NULL); @@ -304,9 +334,11 @@ void colvarbias::communicate_forces() // may send forces to the same colvar // which is why rescaling has to happen now: the colvar is not // aware of this bias' time_step_factor - variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i]); - } - for (i = 0; i < num_variables(); i++) { + if (is_enabled(f_cvb_bypass_ext_lagrangian)) { + variables(i)->add_bias_force_actual_value(cvm::real(time_step_factor) * colvar_forces[i]); + } else { + variables(i)->add_bias_force(cvm::real(time_step_factor) * colvar_forces[i]); + } previous_colvar_forces[i] = colvar_forces[i]; } } @@ -368,20 +400,27 @@ std::string const colvarbias::get_state_params() const int colvarbias::set_state_params(std::string const &conf) { - std::string new_name = ""; - if (colvarparse::get_keyval(conf, "name", new_name, - std::string(""), colvarparse::parse_silent) && - (new_name != this->name)) { - cvm::error("Error: in the state file, the " - "\""+bias_type+"\" block has a different name, \""+new_name+ - "\": different system?\n", INPUT_ERROR); - } + matching_state = false; - if (name.size() == 0) { + std::string check_name = ""; + colvarparse::get_keyval(conf, "name", check_name, + std::string(""), colvarparse::parse_silent); + + if (check_name.size() == 0) { cvm::error("Error: \""+bias_type+"\" block within the restart file " "has no identifiers.\n", INPUT_ERROR); } + if (check_name != this->name) { + if (cvm::debug()) { + cvm::log("Ignoring state of bias \""+check_name+ + "\": this bias is named \""+name+"\".\n"); + } + return COLVARS_OK; + } + + matching_state = true; + colvarparse::get_keyval(conf, "step", state_file_step, cvm::step_absolute(), colvarparse::parse_silent); @@ -396,7 +435,7 @@ std::ostream & colvarbias::write_state(std::ostream &os) } os.setf(std::ios::scientific, std::ios::floatfield); os.precision(cvm::cv_prec); - os << bias_type << " {\n" + os << state_keyword << " {\n" << " configuration {\n"; std::istringstream is(get_state_params()); std::string line; @@ -415,13 +454,12 @@ std::istream & colvarbias::read_state(std::istream &is) size_t const start_pos = is.tellg(); std::string key, brace, conf; - if ( !(is >> key) || !(key == bias_type) || + if ( !(is >> key) || !(key == state_keyword || key == bias_type) || !(is >> brace) || !(brace == "{") || - !(is >> colvarparse::read_block("configuration", conf)) || + !(is >> colvarparse::read_block("configuration", &conf)) || (set_state_params(conf) != COLVARS_OK) ) { - if (key != bias_type) - cvm::log("Found key \"" + key + "\" instead of \"" + bias_type + "\"\n"); - cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+ + cvm::error("Error: in reading state configuration for \""+bias_type+ + "\" bias \""+ this->name+"\" at position "+ cvm::to_str(static_cast(is.tellg()))+ " in stream.\n", INPUT_ERROR); @@ -431,6 +469,15 @@ std::istream & colvarbias::read_state(std::istream &is) return is; } + if (matching_state == false) { + // This state is not for this bias + is.seekg(start_pos, std::ios::beg); + return is; + } + + cvm::log("Restarting "+bias_type+" bias \""+name+"\" from step number "+ + cvm::to_str(state_file_step)+".\n"); + if (!read_state_data(is)) { cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+ this->name+"\" at position "+ diff --git a/lib/colvars/colvarbias.h b/lib/colvars/colvarbias.h index ed72c6063d..5179c42853 100644 --- a/lib/colvars/colvarbias.h +++ b/lib/colvars/colvarbias.h @@ -23,9 +23,12 @@ public: /// Name of this bias std::string name; - /// Type of this bias + /// Keyword indicating the type of this bias std::string bias_type; + /// Keyword used in state files (== bias_type most of the time) + std::string state_keyword; + /// If there is more than one bias of this type, record its rank int rank; @@ -65,7 +68,7 @@ public: virtual int calc_forces(std::vector const *values); /// Send forces to the collective variables - virtual void communicate_forces(); + void communicate_forces(); /// Carry out operations needed before next step is run virtual int end_of_step(); @@ -144,10 +147,10 @@ public: std::istream & read_state_data_key(std::istream &is, char const *key); /// Write the bias configuration to a restart file or other stream - virtual std::ostream & write_state(std::ostream &os); + std::ostream & write_state(std::ostream &os); /// Read the bias configuration from a restart file or other stream - virtual std::istream & read_state(std::istream &is); + std::istream & read_state(std::istream &is); /// Write a label to the trajectory file (comment line) virtual std::ostream & write_traj_label(std::ostream &os); @@ -207,6 +210,9 @@ protected: /// through each colvar object std::vector colvars; + /// \brief Up to date value of each colvar + std::vector colvar_values; + /// \brief Current forces from this bias to the variables std::vector colvar_forces; @@ -226,6 +232,9 @@ protected: /// \brief Step number read from the last state file cvm::step_number state_file_step; + /// Flag used to tell if the state string being read is for this bias + bool matching_state; + }; diff --git a/lib/colvars/colvarbias_abf.cpp b/lib/colvars/colvarbias_abf.cpp index bcd532e680..5382a828f7 100644 --- a/lib/colvars/colvarbias_abf.cpp +++ b/lib/colvars/colvarbias_abf.cpp @@ -115,7 +115,7 @@ int colvarbias_abf::init(std::string const &conf) if (colvars[i]->value().type() != colvarvalue::type_scalar) { cvm::error("Error: ABF bias can only use scalar-type variables.\n"); } - colvars[i]->enable(f_cv_grid); + colvars[i]->enable(f_cv_grid); // Could be a child dependency of a f_cvb_use_grids feature if (hide_Jacobian) { colvars[i]->enable(f_cv_hide_Jacobian); } diff --git a/lib/colvars/colvarbias_alb.cpp b/lib/colvars/colvarbias_alb.cpp index 1d20245fb2..3508270bf6 100644 --- a/lib/colvars/colvarbias_alb.cpp +++ b/lib/colvars/colvarbias_alb.cpp @@ -246,6 +246,12 @@ int colvarbias_alb::update() int colvarbias_alb::set_state_params(std::string const &conf) { + int error_code = colvarbias::set_state_params(conf); + + if (error_code != COLVARS_OK) { + return error_code; + } + if (!get_keyval(conf, "setCoupling", set_coupling)) cvm::fatal_error("Error: current setCoupling is missing from the restart.\n"); diff --git a/lib/colvars/colvarbias_histogram.cpp b/lib/colvars/colvarbias_histogram.cpp index b59a63f6e3..92985b0f61 100644 --- a/lib/colvars/colvarbias_histogram.cpp +++ b/lib/colvars/colvarbias_histogram.cpp @@ -17,6 +17,7 @@ colvarbias_histogram::colvarbias_histogram(char const *key) : colvarbias(key), grid(NULL), out_name("") { + provide(f_cvb_bypass_ext_lagrangian); // Allow histograms of actual cv for extended-Lagrangian } @@ -84,6 +85,10 @@ int colvarbias_histogram::init(std::string const &conf) grid = new colvar_grid_scalar(); grid->init_from_colvars(colvars); + if (is_enabled(f_cvb_bypass_ext_lagrangian)) { + grid->request_actual_value(); + } + { std::string grid_conf; if (key_lookup(conf, "histogramGrid", &grid_conf)) { diff --git a/lib/colvars/colvarbias_meta.cpp b/lib/colvars/colvarbias_meta.cpp index 267d556be6..05f293d872 100644 --- a/lib/colvars/colvarbias_meta.cpp +++ b/lib/colvars/colvarbias_meta.cpp @@ -38,7 +38,12 @@ colvarbias_meta::colvarbias_meta(char const *key) new_hills_begin = hills.end(); hills_traj_os = NULL; + hill_width = 0.0; + + new_hill_freq = 1000; + use_grids = true; + grids_freq = 0; rebin_grids = false; hills_energy = NULL; hills_energy_gradients = NULL; @@ -48,6 +53,8 @@ colvarbias_meta::colvarbias_meta(char const *key) dump_fes_save = false; dump_replica_fes = false; + b_hills_traj = false; + ebmeta_equil_steps = 0L; replica_update_freq = 0; @@ -58,6 +65,7 @@ colvarbias_meta::colvarbias_meta(char const *key) int colvarbias_meta::init(std::string const &conf) { int error_code = COLVARS_OK; + size_t i = 0; error_code |= colvarbias::init(conf); error_code |= colvarbias_ti::init(conf); @@ -71,16 +79,38 @@ int colvarbias_meta::init(std::string const &conf) cvm::error("Error: hillWeight must be provided, and a positive number.\n", INPUT_ERROR); } - get_keyval(conf, "newHillFrequency", new_hill_freq, 1000); + get_keyval(conf, "newHillFrequency", new_hill_freq, new_hill_freq); if (new_hill_freq > 0) { enable(f_cvb_history_dependent); + if (grids_freq == 0) { + grids_freq = new_hill_freq; + } } - get_keyval(conf, "hillWidth", hill_width, cvm::sqrt(2.0 * PI) / 2.0); - cvm::log("Half-widths of the Gaussian hills (sigma's):\n"); - for (size_t i = 0; i < num_variables(); i++) { - cvm::log(variables(i)->name+std::string(": ")+ - cvm::to_str(0.5 * variables(i)->width * hill_width)); + get_keyval(conf, "gaussianSigmas", colvar_sigmas, colvar_sigmas); + + get_keyval(conf, "hillWidth", hill_width, hill_width); + + if ((colvar_sigmas.size() > 0) && (hill_width > 0.0)) { + error_code |= cvm::error("Error: hillWidth and gaussianSigmas are " + "mutually exclusive.", INPUT_ERROR); + } + + if (hill_width > 0.0) { + colvar_sigmas.resize(num_variables()); + // Print the calculated sigma parameters + cvm::log("Half-widths of the Gaussian hills (sigma's):\n"); + for (i = 0; i < num_variables(); i++) { + colvar_sigmas[i] = variables(i)->width * hill_width / 2.0; + cvm::log(variables(i)->name+std::string(": ")+ + cvm::to_str(colvar_sigmas[i])); + } + } + + if (colvar_sigmas.size() == 0) { + error_code |= cvm::error("Error: positive values are required for " + "either hillWidth or gaussianSigmas.", + INPUT_ERROR); } { @@ -95,11 +125,18 @@ int colvarbias_meta::init(std::string const &conf) get_keyval(conf, "useGrids", use_grids, use_grids); if (use_grids) { - get_keyval(conf, "gridsUpdateFrequency", grids_freq, new_hill_freq); + + for (i = 0; i < num_variables(); i++) { + if (2.0*colvar_sigmas[i] < variables(i)->width) { + cvm::log("Warning: gaussianSigmas is too narrow for the grid " + "spacing along "+variables(i)->name+"."); + } + } + + get_keyval(conf, "gridsUpdateFrequency", grids_freq, grids_freq); get_keyval(conf, "rebinGrids", rebin_grids, rebin_grids); expand_grids = false; - size_t i; for (i = 0; i < num_variables(); i++) { variables(i)->enable(f_cv_grid); // Could be a child dependency of a f_cvb_use_grids feature if (variables(i)->expand_boundaries) { @@ -116,15 +153,17 @@ int colvarbias_meta::init(std::string const &conf) get_keyval(conf, "keepHills", keep_hills, keep_hills); get_keyval(conf, "keepFreeEnergyFiles", dump_fes_save, dump_fes_save); - hills_energy = new colvar_grid_scalar(colvars); - hills_energy_gradients = new colvar_grid_gradient(colvars); + if (hills_energy == NULL) { + hills_energy = new colvar_grid_scalar(colvars); + hills_energy_gradients = new colvar_grid_gradient(colvars); + } } else { dump_fes = false; } - get_keyval(conf, "writeHillsTrajectory", b_hills_traj, false); + get_keyval(conf, "writeHillsTrajectory", b_hills_traj, b_hills_traj); error_code |= init_replicas_params(conf); error_code |= init_well_tempered_params(conf); @@ -323,7 +362,7 @@ int colvarbias_meta::clear_state_data() // ********************************************************************** std::list::const_iterator -colvarbias_meta::create_hill(colvarbias_meta::hill const &h) +colvarbias_meta::add_hill(colvarbias_meta::hill const &h) { hill_iter const hills_end = hills.end(); hills.push_back(h); @@ -563,12 +602,14 @@ int colvarbias_meta::update_bias() case single_replica: - create_hill(hill(hill_weight*hills_scale, colvars, hill_width)); + add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale, + colvar_values, colvar_sigmas)); break; case multiple_replicas: - create_hill(hill(hill_weight*hills_scale, colvars, hill_width, replica_id)); + add_hill(hill(cvm::step_absolute(), hill_weight*hills_scale, + colvar_values, colvar_sigmas, replica_id)); std::ostream *replica_hills_os = cvm::proxy->get_output_stream(replica_hills_file); if (replica_hills_os) { @@ -730,35 +771,21 @@ void colvarbias_meta::calc_hills(colvarbias_meta::hill_iter h_first, std::vector const *values) { size_t i = 0; - std::vector curr_values(num_variables()); - for (i = 0; i < num_variables(); i++) { - curr_values[i].type(variables(i)->value()); - } - - if (values) { - for (i = 0; i < num_variables(); i++) { - curr_values[i] = (*values)[i]; - } - } else { - for (i = 0; i < num_variables(); i++) { - curr_values[i] = variables(i)->value(); - } - } for (hill_iter h = h_first; h != h_last; h++) { // compute the gaussian exponent cvm::real cv_sqdev = 0.0; for (i = 0; i < num_variables(); i++) { - colvarvalue const &x = curr_values[i]; + colvarvalue const &x = values ? (*values)[i] : colvar_values[i]; colvarvalue const ¢er = h->centers[i]; - cvm::real const half_width = 0.5 * h->widths[i]; - cv_sqdev += (variables(i)->dist2(x, center)) / (half_width*half_width); + cvm::real const sigma = h->sigmas[i]; + cv_sqdev += (variables(i)->dist2(x, center)) / (sigma*sigma); } // compute the gaussian if (cv_sqdev > 23.0) { - // set it to zero if the exponent is more negative than log(1.0E-05) + // set it to zero if the exponent is more negative than log(1.0E-06) h->value(0.0); } else { h->value(cvm::exp(-0.5*cv_sqdev)); @@ -775,22 +802,22 @@ void colvarbias_meta::calc_hills_force(size_t const &i, std::vector const *values) { // Retrieve the value of the colvar - colvarvalue const x(values ? (*values)[i] : variables(i)->value()); + colvarvalue const x(values ? (*values)[i] : colvar_values[i]); // 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 (variables(i)->value().type()) { + switch (x.type()) { case colvarvalue::type_scalar: for (h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; - cvm::real const half_width = 0.5 * h->widths[i]; + cvm::real const sigma = h->sigmas[i]; forces[i].real_value += - ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * + ( h->weight() * h->value() * (0.5 / (sigma*sigma)) * (variables(i)->dist2_lgrad(x, center)).real_value ); } break; @@ -801,9 +828,9 @@ void colvarbias_meta::calc_hills_force(size_t const &i, for (h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; - cvm::real const half_width = 0.5 * h->widths[i]; + cvm::real const sigma = h->sigmas[i]; forces[i].rvector_value += - ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * + ( h->weight() * h->value() * (0.5 / (sigma*sigma)) * (variables(i)->dist2_lgrad(x, center)).rvector_value ); } break; @@ -813,9 +840,9 @@ void colvarbias_meta::calc_hills_force(size_t const &i, for (h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; - cvm::real const half_width = 0.5 * h->widths[i]; + cvm::real const sigma = h->sigmas[i]; forces[i].quaternion_value += - ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * + ( h->weight() * h->value() * (0.5 / (sigma*sigma)) * (variables(i)->dist2_lgrad(x, center)).quaternion_value ); } break; @@ -824,9 +851,9 @@ void colvarbias_meta::calc_hills_force(size_t const &i, for (h = h_first; h != h_last; h++) { if (h->value() == 0.0) continue; colvarvalue const ¢er = h->centers[i]; - cvm::real const half_width = 0.5 * h->widths[i]; + cvm::real const sigma = h->sigmas[i]; forces[i].vector1d_value += - ( h->weight() * h->value() * (0.5 / (half_width*half_width)) * + ( h->weight() * h->value() * (0.5 / (sigma*sigma)) * (variables(i)->dist2_lgrad(x, center)).vector1d_value ); } break; @@ -1246,14 +1273,20 @@ void colvarbias_meta::read_replica_files() int colvarbias_meta::set_state_params(std::string const &state_conf) { - std::string new_replica = ""; - if (colvarparse::get_keyval(state_conf, "replicaID", new_replica, + int error_code = colvarbias::set_state_params(state_conf); + + if (error_code != COLVARS_OK) { + return error_code; + } + + std::string check_replica = ""; + if (colvarparse::get_keyval(state_conf, "replicaID", check_replica, std::string(""), colvarparse::parse_silent) && - (new_replica != this->replica_id)) { - cvm::error("Error: in the state file, the " - "\"metadynamics\" block has a different replicaID: different system?\n", - INPUT_ERROR); - return INPUT_ERROR; + (check_replica != this->replica_id)) { + return cvm::error("Error: in the state file , the " + "\"metadynamics\" block has a different replicaID ("+ + check_replica+" instead of "+replica_id+").\n", + INPUT_ERROR); } return COLVARS_OK; @@ -1479,17 +1512,18 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) if (!is) return is; // do nothing if failbit is set size_t const start_pos = is.tellg(); + size_t i = 0; std::string data; - if ( !(is >> read_block("hill", data)) ) { + if ( !(is >> read_block("hill", &data)) ) { is.clear(); is.seekg(start_pos, std::ios::beg); is.setstate(std::ios::failbit); return is; } - cvm::step_number h_it; - get_keyval(data, "step", h_it, 0L, parse_silent); + cvm::step_number h_it = 0L; + get_keyval(data, "step", h_it, h_it, parse_restart); if (h_it <= state_file_step) { if (cvm::debug()) cvm::log("Skipping a hill older than the state file for metadynamics bias \""+ @@ -1499,31 +1533,24 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) } cvm::real h_weight; - get_keyval(data, "weight", h_weight, hill_weight, parse_silent); + get_keyval(data, "weight", h_weight, hill_weight, parse_restart); std::vector h_centers(num_variables()); - for (size_t i = 0; i < num_variables(); i++) { + for (i = 0; i < num_variables(); i++) { h_centers[i].type(variables(i)->value()); } - { - // it is safer to read colvarvalue objects one at a time; - // TODO: change this it later - std::string centers_input; - key_lookup(data, "centers", ¢ers_input); - std::istringstream centers_is(centers_input); - for (size_t i = 0; i < num_variables(); i++) { - centers_is >> h_centers[i]; - } - } + get_keyval(data, "centers", h_centers, h_centers, parse_restart); - std::vector h_widths(num_variables()); - get_keyval(data, "widths", h_widths, - std::vector(num_variables(), (cvm::sqrt(2.0 * PI) / 2.0)), - parse_silent); + std::vector h_sigmas(num_variables()); + get_keyval(data, "widths", h_sigmas, h_sigmas, parse_restart); + for (i = 0; i < num_variables(); i++) { + // For backward compatibility, read the widths instead of the sigmas + h_sigmas[i] /= 2.0; + } std::string h_replica = ""; if (comm != single_replica) { - get_keyval(data, "replicaID", h_replica, replica_id, parse_silent); + get_keyval(data, "replicaID", h_replica, replica_id, parse_restart); if (h_replica != replica_id) cvm::fatal_error("Error: trying to read a hill created by replica \""+h_replica+ "\" for replica \""+replica_id+ @@ -1531,7 +1558,7 @@ std::istream & colvarbias_meta::read_hill(std::istream &is) } hill_iter const hills_end = hills.end(); - hills.push_back(hill(h_it, h_weight, h_centers, h_widths, h_replica)); + hills.push_back(hill(h_it, h_weight, h_centers, h_sigmas, h_replica)); if (new_hills_begin == hills_end) { // if new_hills_begin is unset, set it for the first time new_hills_begin = hills.end(); @@ -1753,8 +1780,8 @@ void colvarbias_meta::write_pmf() if (ebmeta) { int nt_points=pmf->number_of_points(); for (int i = 0; i < nt_points; i++) { - cvm:: real pmf_val=0.0; - cvm:: real target_val=target_dist->value(i); + cvm::real pmf_val=0.0; + cvm::real target_val=target_dist->value(i); if (target_val>0) { pmf_val=pmf->value(i); pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val); @@ -1794,8 +1821,8 @@ void colvarbias_meta::write_pmf() if (ebmeta) { int nt_points=pmf->number_of_points(); for (int i = 0; i < nt_points; i++) { - cvm:: real pmf_val=0.0; - cvm:: real target_val=target_dist->value(i); + cvm::real pmf_val=0.0; + cvm::real target_val=target_dist->value(i); if (target_val>0) { pmf_val=pmf->value(i); pmf_val=pmf_val+cvm::temperature() * cvm::boltzmann() * cvm::logn(target_val); @@ -1889,10 +1916,10 @@ std::string colvarbias_meta::hill::output_traj() } os << " "; - for (i = 0; i < widths.size(); i++) { + for (i = 0; i < sigmas.size(); i++) { os << " "; os << std::setprecision(cvm::cv_prec) - << std::setw(cvm::cv_width) << widths[i]; + << std::setw(cvm::cv_width) << sigmas[i]; } os << " "; @@ -1903,6 +1930,47 @@ std::string colvarbias_meta::hill::output_traj() } +colvarbias_meta::hill::hill(cvm::step_number it_in, + cvm::real W_in, + std::vector const &cv_values, + std::vector const &cv_sigmas, + std::string const &replica_in) + : it(it_in), + sW(1.0), + W(W_in), + centers(cv_values.size()), + sigmas(cv_values.size()), + replica(replica_in) +{ + for (size_t i = 0; i < cv_values.size(); i++) { + centers[i].type(cv_values[i]); + centers[i] = cv_values[i]; + sigmas[i] = cv_sigmas[i]; + } + if (cvm::debug()) { + cvm::log("New hill, applied to "+cvm::to_str(cv_values.size())+ + " collective variables, with centers "+ + cvm::to_str(centers)+", sigmas "+ + cvm::to_str(sigmas)+" and weight "+ + cvm::to_str(W)+".\n"); + } +} + + +colvarbias_meta::hill::hill(colvarbias_meta::hill const &h) + : sW(1.0), + W(h.W), + centers(h.centers), + sigmas(h.sigmas), + it(h.it), + replica(h.replica) +{} + + +colvarbias_meta::hill::~hill() +{} + + std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h) { os.setf(std::ios::scientific, std::ios::floatfield); @@ -1927,12 +1995,13 @@ std::ostream & operator << (std::ostream &os, colvarbias_meta::hill const &h) } os << "\n"; + // For backward compatibility, write the widths instead of the sigmas os << " widths "; - for (i = 0; i < (h.widths).size(); i++) { + for (i = 0; i < (h.sigmas).size(); i++) { os << " " << std::setprecision(cvm::cv_prec) << std::setw(cvm::cv_width) - << h.widths[i]; + << 2.0 * h.sigmas[i]; } os << "\n"; diff --git a/lib/colvars/colvarbias_meta.h b/lib/colvars/colvarbias_meta.h index abf1417c32..9f9fc63fc7 100644 --- a/lib/colvars/colvarbias_meta.h +++ b/lib/colvars/colvarbias_meta.h @@ -71,11 +71,14 @@ public: protected: - /// \brief width of a hill + /// Width of a hill in number of grid points /// /// The local width of each collective variable, multiplied by this /// number, provides the hill width along that direction - cvm::real hill_width; + cvm::real hill_width; + + /// The sigma parameters of the Gaussian hills + std::vector colvar_sigmas; /// \brief Number of simulation steps between two hills size_t new_hill_freq; @@ -113,11 +116,11 @@ protected: /// \brief Add a new hill; if a .hills trajectory is written, /// write it there; if there is more than one replica, communicate /// it to the others - virtual std::list::const_iterator create_hill(hill const &h); + std::list::const_iterator add_hill(hill const &h); /// \brief Remove a previously saved hill (returns an iterator for /// the next hill in the list) - virtual std::list::const_iterator delete_hill(hill_iter &h); + std::list::const_iterator delete_hill(hill_iter &h); /// \brief Calculate the values of the hills, incrementing /// bias_energy @@ -266,6 +269,9 @@ class colvarbias_meta::hill { protected: + /// Time step at which this hill was added + cvm::step_number it; + /// Value of the hill function (ranges between 0 and 1) cvm::real hill_value; @@ -275,83 +281,35 @@ protected: /// Maximum height in energy of the hill cvm::real W; - /// Center of the hill in the collective variable space - std::vector centers; + /// Centers of the hill in the collective variable space + std::vector centers; - /// Widths of the hill in the collective variable space - std::vector widths; + /// Half-widths of the hill in the collective variable space + std::vector sigmas; + + /// Identity of the replica who added this hill + std::string replica; public: friend class colvarbias_meta; - /// Time step at which this hill was added - cvm::step_number it; - - /// Identity of the replica who added this hill (only in multiple replica simulations) - std::string replica; - - /// \brief Runtime constructor: data are read directly from - /// collective variables \param weight Weight of the hill \param - /// cv Pointer to the array of collective variables involved \param - /// replica (optional) Identity of the replica which creates the - /// hill - inline hill(cvm::real const &W_in, - std::vector &cv, - cvm::real const &hill_width, - std::string const &replica_in = "") - : sW(1.0), - W(W_in), - centers(cv.size()), - widths(cv.size()), - it(cvm::step_absolute()), - replica(replica_in) - { - for (size_t i = 0; i < cv.size(); i++) { - centers[i].type(cv[i]->value()); - centers[i] = cv[i]->value(); - widths[i] = cv[i]->width * hill_width; - } - if (cvm::debug()) - cvm::log("New hill, applied to "+cvm::to_str(cv.size())+ - " collective variables, with centers "+ - cvm::to_str(centers)+", widths "+ - cvm::to_str(widths)+" and weight "+ - cvm::to_str(W)+".\n"); - } - - /// \brief General constructor: all data are explicitly passed as - /// arguments (used for instance when reading hills saved on a - /// file) \param it Time step of creation of the hill \param - /// weight Weight of the hill \param centers Center of the hill - /// \param widths Width of the hill around centers \param replica - /// (optional) Identity of the replica which creates the hill - inline hill(cvm::step_number const &it_in, - cvm::real const &W_in, - std::vector const ¢ers_in, - std::vector const &widths_in, - std::string const &replica_in = "") - : sW(1.0), - W(W_in), - centers(centers_in), - widths(widths_in), - it(it_in), - replica(replica_in) - {} + /// Constructor of a hill object + /// \param it Step number at which the hill was added + /// \param W Weight of the hill (energy units) + /// \param cv_values Array of collective variable values + /// \param cv_sigmas Array of collective variable values + /// \param replica ID of the replica that creates the hill (optional) + hill(cvm::step_number it, cvm::real W, + std::vector const &cv_values, + std::vector const &cv_sigmas, + std::string const &replica = ""); /// Copy constructor - inline hill(colvarbias_meta::hill const &h) - : sW(1.0), - W(h.W), - centers(h.centers), - widths(h.widths), - it(h.it), - replica(h.replica) - {} + hill(colvarbias_meta::hill const &h); /// Destructor - inline ~hill() - {} + ~hill(); /// Get the energy inline cvm::real energy() @@ -440,8 +398,7 @@ public: std::string output_traj(); /// Write the hill to an output stream - inline friend std::ostream & operator << (std::ostream &os, - hill const &h); + friend std::ostream & operator << (std::ostream &os, hill const &h); }; diff --git a/lib/colvars/colvarbias_restraint.cpp b/lib/colvars/colvarbias_restraint.cpp index 1fd9c84e98..486d77688b 100644 --- a/lib/colvars/colvarbias_restraint.cpp +++ b/lib/colvars/colvarbias_restraint.cpp @@ -17,6 +17,7 @@ colvarbias_restraint::colvarbias_restraint(char const *key) : colvarbias(key), colvarbias_ti(key) { + state_keyword = "restraint"; } @@ -705,67 +706,6 @@ std::ostream & colvarbias_restraint_k_moving::write_traj(std::ostream &os) -// redefined due to legacy state file keyword "harmonic" -std::istream & colvarbias_restraint::read_state(std::istream &is) -{ - size_t const start_pos = is.tellg(); - - std::string key, brace, conf; - if ( !(is >> key) || !(key == "restraint" || key == "harmonic") || - !(is >> brace) || !(brace == "{") || - !(is >> colvarparse::read_block("configuration", conf)) || - (set_state_params(conf) != COLVARS_OK) ) { - cvm::error("Error: in reading state configuration for \""+bias_type+"\" bias \""+ - this->name+"\" at position "+ - cvm::to_str(static_cast(is.tellg()))+ - " in stream.\n", INPUT_ERROR); - is.clear(); - is.seekg(start_pos, std::ios::beg); - is.setstate(std::ios::failbit); - return is; - } - - if (!read_state_data(is)) { - cvm::error("Error: in reading state data for \""+bias_type+"\" bias \""+ - this->name+"\" at position "+ - cvm::to_str(static_cast(is.tellg()))+ - " in stream.\n", INPUT_ERROR); - is.clear(); - is.seekg(start_pos, std::ios::beg); - is.setstate(std::ios::failbit); - } - - is >> brace; - if (brace != "}") { - cvm::log("brace = "+brace+"\n"); - cvm::error("Error: corrupt restart information for \""+bias_type+"\" bias \""+ - this->name+"\": no matching brace at position "+ - cvm::to_str(static_cast(is.tellg()))+" in stream.\n"); - is.setstate(std::ios::failbit); - } - - return is; -} - - -std::ostream & colvarbias_restraint::write_state(std::ostream &os) -{ - os.setf(std::ios::scientific, std::ios::floatfield); - os << "restraint {\n" - << " configuration {\n"; - std::istringstream is(get_state_params()); - std::string line; - while (std::getline(is, line)) { - os << " " << line << "\n"; - } - os << " }\n"; - write_state_data(os); - os << "}\n\n"; - return os; -} - - - colvarbias_restraint_harmonic::colvarbias_restraint_harmonic(char const *key) : colvarbias(key), colvarbias_ti(key), @@ -928,6 +868,9 @@ colvarbias_restraint_harmonic_walls::colvarbias_restraint_harmonic_walls(char co { lower_wall_k = -1.0; upper_wall_k = -1.0; + // This bias implements the bias_actual_colvars feature (most others do not) + provide(f_cvb_bypass_ext_lagrangian); + set_enabled(f_cvb_bypass_ext_lagrangian); // Defaults to enabled } @@ -1075,23 +1018,13 @@ int colvarbias_restraint_harmonic_walls::update() } -void colvarbias_restraint_harmonic_walls::communicate_forces() -{ - for (size_t i = 0; i < num_variables(); i++) { - if (cvm::debug()) { - cvm::log("Communicating a force to colvar \""+ - variables(i)->name+"\".\n"); - } - // Impulse-style multiple timestep - variables(i)->add_bias_force_actual_value(cvm::real(time_step_factor) * colvar_forces[i]); - } -} - - cvm::real colvarbias_restraint_harmonic_walls::colvar_distance(size_t i) const { colvar *cv = variables(i); - colvarvalue const &cvv = variables(i)->actual_value(); + + colvarvalue const &cvv = is_enabled(f_cvb_bypass_ext_lagrangian) ? + variables(i)->actual_value() : + variables(i)->value(); // For a periodic colvar, both walls may be applicable at the same time // in which case we pick the closer one diff --git a/lib/colvars/colvarbias_restraint.h b/lib/colvars/colvarbias_restraint.h index 644835f129..ac25e35646 100644 --- a/lib/colvars/colvarbias_restraint.h +++ b/lib/colvars/colvarbias_restraint.h @@ -35,8 +35,6 @@ public: virtual int set_state_params(std::string const &conf); // virtual std::ostream & write_state_data(std::ostream &os); // virtual std::istream & read_state_data(std::istream &os); - virtual std::ostream & write_state(std::ostream &os); - virtual std::istream & read_state(std::istream &is); virtual std::ostream & write_traj_label(std::ostream &os); virtual std::ostream & write_traj(std::ostream &os); @@ -258,7 +256,6 @@ public: colvarbias_restraint_harmonic_walls(char const *key); virtual int init(std::string const &conf); virtual int update(); - virtual void communicate_forces(); virtual std::string const get_state_params() const; virtual int set_state_params(std::string const &conf); virtual std::ostream & write_state_data(std::ostream &os); diff --git a/lib/colvars/colvardeps.cpp b/lib/colvars/colvardeps.cpp index 8425982be4..f83458b51c 100644 --- a/lib/colvars/colvardeps.cpp +++ b/lib/colvars/colvardeps.cpp @@ -428,11 +428,11 @@ void colvardeps::require_feature_alt(int f, int g, int h, int i, int j) { void colvardeps::print_state() { size_t i; - cvm::log("Enabled features of \"" + description + "\" (with reference count)"); + cvm::log("Features of \"" + description + "\" ON/OFF (refcount)"); for (i = 0; i < feature_states.size(); i++) { - if (is_enabled(i)) - cvm::log("- " + features()[i]->description + " (" - + cvm::to_str(feature_states[i].ref_count) + ")"); + std::string onoff = is_enabled(i) ? "ON" : "OFF"; + cvm::log("- " + features()[i]->description + " " + onoff + " (" + + cvm::to_str(feature_states[i].ref_count) + ")"); } cvm::increase_depth(); for (i=0; i> key) && (key == std::string("grid_parameters"))) { is.seekg(start_pos, std::ios::beg); - is >> colvarparse::read_block("grid_parameters", conf); + is >> colvarparse::read_block("grid_parameters", &conf); parse_params(conf, colvarparse::parse_silent); } else { cvm::log("Grid parameters are missing in the restart file, using those from the configuration.\n"); diff --git a/lib/colvars/colvarmodule.cpp b/lib/colvars/colvarmodule.cpp index ecf0ae800d..83632d9458 100644 --- a/lib/colvars/colvarmodule.cpp +++ b/lib/colvars/colvarmodule.cpp @@ -1284,7 +1284,7 @@ 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 (is >> colvarparse::read_block("configuration", &restart_conf)) { parse->get_keyval(restart_conf, "step", it_restart, static_cast(0), @@ -1329,30 +1329,84 @@ std::istream & colvarmodule::read_restart(std::istream &is) parse->clear_keyword_registry(); } - // colvars restart - cvm::increase_depth(); - for (std::vector::iterator cvi = colvars.begin(); - cvi != colvars.end(); - cvi++) { - if ( !((*cvi)->read_restart(is)) ) { - cvm::error("Error: in reading restart configuration for collective variable \""+ - (*cvi)->name+"\".\n", - INPUT_ERROR); + print_total_forces_errning(warn_total_forces); + + read_objects_state(is); + + return is; +} + + + +std::istream & colvarmodule::read_objects_state(std::istream &is) +{ + size_t pos = 0; + std::string word; + + while (is.good()) { + pos = is.tellg(); + word.clear(); + is >> word; + + if (word.size()) { + + is.seekg(pos, std::ios::beg); + + if (word == "colvar") { + + cvm::increase_depth(); + for (std::vector::iterator cvi = colvars.begin(); + cvi != colvars.end(); + cvi++) { + if ( !((*cvi)->read_state(is)) ) { + // Here an error signals that the variable is a match, but the + // state is corrupt; otherwise, the variable rewinds is silently + cvm::error("Error: in reading restart configuration for " + "collective variable \""+(*cvi)->name+"\".\n", + INPUT_ERROR); + } + if (static_cast(is.tellg()) > pos) break; // found it + } + cvm::decrease_depth(); + + } else { + + cvm::increase_depth(); + for (std::vector::iterator bi = biases.begin(); + bi != biases.end(); + bi++) { + if (((*bi)->state_keyword != word) && (*bi)->bias_type != word) { + // Skip biases with different type; state_keyword is used to + // support different versions of the state file format + continue; + } + if (!((*bi)->read_state(is))) { + // Same as above, an error means a match but the state is incorrect + cvm::error("Error: in reading restart configuration for bias \""+ + (*bi)->name+"\".\n", + INPUT_ERROR); + } + if (static_cast(is.tellg()) > pos) break; // found it + } + cvm::decrease_depth(); + } } + + if (static_cast(is.tellg()) == pos) { + // This block has not been read by any object: discard it and move on + // to the next one + is >> colvarparse::read_block(word, NULL); + } + + if (!is) break; } - // biases restart - for (std::vector::iterator bi = biases.begin(); - bi != biases.end(); - bi++) { - if (!((*bi)->read_state(is))) { - cvm::error("Error: in reading restart configuration for bias \""+ - (*bi)->name+"\".\n", - INPUT_ERROR); - } - } - cvm::decrease_depth(); + return is; +} + +int colvarmodule::print_total_forces_errning(bool warn_total_forces) +{ if (warn_total_forces) { cvm::log(cvm::line_marker); cvm::log("WARNING: The definition of system forces has changed. Please see:\n"); @@ -1367,10 +1421,11 @@ to:\n\ and load it to continue this simulation.\n"); output_prefix() = output_prefix()+".tmp"; write_restart_file(output_prefix()+".colvars.state"); - cvm::error("Exiting with error until issue is addressed.\n", INPUT_ERROR); + return cvm::error("Exiting with error until issue is addressed.\n", + INPUT_ERROR); } - return is; + return COLVARS_OK; } @@ -1499,7 +1554,7 @@ std::ostream & colvarmodule::write_restart(std::ostream &os) for (std::vector::iterator cvi = colvars.begin(); cvi != colvars.end(); cvi++) { - (*cvi)->write_restart(os); + (*cvi)->write_state(os); } for (std::vector::iterator bi = biases.begin(); diff --git a/lib/colvars/colvarmodule.h b/lib/colvars/colvarmodule.h index d3bc5e1b80..9346e8e1d6 100644 --- a/lib/colvars/colvarmodule.h +++ b/lib/colvars/colvarmodule.h @@ -429,8 +429,15 @@ public: /// (Re)initialize the output trajectory and state file (does not write it yet) int setup_output(); - /// Read the input restart file + /// Read a restart file std::istream & read_restart(std::istream &is); + + /// Read the states of individual objects; allows for changes + std::istream & read_objects_state(std::istream &is); + + /// If needed (old restart file), print the warning that cannot be ignored + int print_total_forces_errning(bool warn_total_forces); + /// Write the output restart file std::ostream & write_restart(std::ostream &os); diff --git a/lib/colvars/colvarparse.cpp b/lib/colvars/colvarparse.cpp index 77f03793a3..a09d6be9ef 100644 --- a/lib/colvars/colvarparse.cpp +++ b/lib/colvars/colvarparse.cpp @@ -840,6 +840,17 @@ bool colvarparse::key_lookup(std::string const &conf, } +colvarparse::read_block::read_block(std::string const &key_in, + std::string *data_in) + : key(key_in), data(data_in) +{ +} + + +colvarparse::read_block::~read_block() +{} + + std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb) { size_t start_pos = is.tellg(); @@ -856,7 +867,9 @@ std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb) } if (next != "{") { - (*rb.data) = next; + if (rb.data) { + *(rb.data) = next; + } return is; } @@ -870,9 +883,15 @@ std::istream & operator>> (std::istream &is, colvarparse::read_block const &rb) br_old = br; br++; } - if (brace_count) (*rb.data).append(line + "\n"); + if (brace_count) { + if (rb.data) { + (rb.data)->append(line + "\n"); + } + } else { - (*rb.data).append(line, 0, br_old); + if (rb.data) { + (rb.data)->append(line, 0, br_old); + } break; } } diff --git a/lib/colvars/colvarparse.h b/lib/colvars/colvarparse.h index f43f6527c8..b7d42fdffa 100644 --- a/lib/colvars/colvarparse.h +++ b/lib/colvars/colvarparse.h @@ -272,14 +272,18 @@ public: /// skipping other blocks class read_block { - std::string const key; + /// The keyword that identifies the block + std::string const key; + + /// Where to keep the data (may be NULL) std::string * const data; public: - inline read_block(std::string const &key_in, std::string &data_in) - : key(key_in), data(&data_in) - {} - inline ~read_block() {} + + read_block(std::string const &key_in, std::string *data_in = NULL); + + ~read_block(); + friend std::istream & operator >> (std::istream &is, read_block const &rb); }; diff --git a/lib/colvars/colvars_version.h b/lib/colvars/colvars_version.h index aa985bedd8..cc8744d62b 100644 --- a/lib/colvars/colvars_version.h +++ b/lib/colvars/colvars_version.h @@ -1,3 +1,3 @@ #ifndef COLVARS_VERSION -#define COLVARS_VERSION "2020-01-27" +#define COLVARS_VERSION "2020-02-25" #endif diff --git a/src/USER-COLVARS/colvarproxy_lammps.cpp b/src/USER-COLVARS/colvarproxy_lammps.cpp index 2129461c9a..dabb0bf923 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.cpp +++ b/src/USER-COLVARS/colvarproxy_lammps.cpp @@ -160,14 +160,20 @@ void colvarproxy_lammps::init(const char *conf_file) } } -void colvarproxy_lammps::add_config_file(const char *conf_file) +int colvarproxy_lammps::add_config_file(const char *conf_file) { - colvars->read_config_file(conf_file); + return colvars->read_config_file(conf_file); } -void colvarproxy_lammps::add_config_string(const std::string &conf) +int colvarproxy_lammps::add_config_string(const std::string &conf) { - colvars->read_config_string(conf); + return colvars->read_config_string(conf); +} + +int colvarproxy_lammps::read_state_file(char const *state_filename) +{ + input_prefix() = std::string(state_filename); + return colvars->setup_input(); } colvarproxy_lammps::~colvarproxy_lammps() diff --git a/src/USER-COLVARS/colvarproxy_lammps.h b/src/USER-COLVARS/colvarproxy_lammps.h index 29ba5023ed..95e7d4fb32 100644 --- a/src/USER-COLVARS/colvarproxy_lammps.h +++ b/src/USER-COLVARS/colvarproxy_lammps.h @@ -86,10 +86,13 @@ class colvarproxy_lammps : public colvarproxy { void write_output_files(); // read additional config from file - void add_config_file(char const *config_filename); + int add_config_file(char const *config_filename); // read additional config from string - void add_config_string(const std::string &config); + int add_config_string(const std::string &config); + + // load a state file + int read_state_file(char const *state_filename); // Request to set the units used internally by Colvars int set_unit_system(std::string const &units_in, bool check_only); diff --git a/src/USER-COLVARS/fix_colvars.cpp b/src/USER-COLVARS/fix_colvars.cpp index 9863e8c20f..7ce2a93680 100644 --- a/src/USER-COLVARS/fix_colvars.cpp +++ b/src/USER-COLVARS/fix_colvars.cpp @@ -502,7 +502,7 @@ int FixColvars::modify_param(int narg, char **arg) if (me == 0) { if (! proxy) error->one(FLERR,"Cannot use fix_modify before initialization"); - proxy->add_config_file(arg[1]); + return proxy->add_config_file(arg[1]) == COLVARS_OK ? 0 : 2; } return 2; } else if (strcmp(arg[0],"config") == 0) { @@ -510,8 +510,16 @@ int FixColvars::modify_param(int narg, char **arg) if (me == 0) { if (! proxy) error->one(FLERR,"Cannot use fix_modify before initialization"); - std::string conf(arg[1]); - proxy->add_config_string(conf); + std::string const conf(arg[1]); + return proxy->add_config_string(conf) == COLVARS_OK ? 0 : 2; + } + return 2; + } else if (strcmp(arg[0],"load") == 0) { + if (narg < 2) error->all(FLERR,"Illegal fix_modify command"); + if (me == 0) { + if (! proxy) + error->one(FLERR,"Cannot use fix_modify before initialization"); + return proxy->read_state_file(arg[1]) == COLVARS_OK ? 0 : 2; } return 2; } diff --git a/src/lmptype.h b/src/lmptype.h index 65e46535fc..5c19ab7483 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -28,6 +28,12 @@ #ifndef LMP_LMPTYPE_H #define LMP_LMPTYPE_H +// C++11 check + +#if __cplusplus < 201103L +#error LAMMPS requires a C++11 (or later) compliant compiler. Enable C++11 compatibility or upgrade the compiler. +#endif + #ifndef __STDC_LIMIT_MACROS #define __STDC_LIMIT_MACROS #endif @@ -38,8 +44,8 @@ #include #include -#include // requires C++-11 -#include // requires C++-11 +#include +#include // grrr - IBM Power6 does not provide this def in their system header files