diff --git a/.github/workflows/backend-template.yml b/.github/workflows/backend-template.yml index 84e31e4f5..3be97e6b1 100644 --- a/.github/workflows/backend-template.yml +++ b/.github/workflows/backend-template.yml @@ -121,9 +121,6 @@ jobs: ref: 'master' path: 'devel-tools/packages' - - name: Test build recipes - run: bash devel-tools/check_build_recipes - - name: Build and test standalone library run: cmake -P devel-tools/build_test_library.cmake diff --git a/.github/workflows/test-library.yml b/.github/workflows/test-library.yml index d034195f2..ba09aac96 100644 --- a/.github/workflows/test-library.yml +++ b/.github/workflows/test-library.yml @@ -39,9 +39,6 @@ jobs: # run: | # sudo apt -y install ccache - - name: Test build recipes - run: bash devel-tools/check_build_recipes - - name: Convert BibTeX references to code shell: bash working-directory: doc diff --git a/devel-tools/check_build_recipes b/devel-tools/check_build_recipes deleted file mode 100755 index ec1197039..000000000 --- a/devel-tools/check_build_recipes +++ /dev/null @@ -1,26 +0,0 @@ -#!/bin/bash - -# TODO check later for headers in the dependencies; conditional compilation -# may give false results - -echo "Checking Colvars build recipes for VMD" - -pushd $(git rev-parse --show-toplevel) 1> /dev/null -retcode=0 -for f in src/colvar*.cpp src/colvar*.h ; do - f=$(basename $f) - extension=${f#*.} - if [ ${extension} = 'cpp' ] ; then - prefix=${f%.cpp} - vmd_source=${prefix}.C - else - prefix=${f%.h} - vmd_source=${f} - fi - if ! grep -q ${vmd_source} vmd/src/colvars_files.pl ; then - echo "Error: ${vmd_source} missing from vmd/src/colvars_files.pl" >&2 - retcode=1 - fi -done -popd 1> /dev/null -exit ${retcode} diff --git a/devel-tools/compile-gromacs.sh b/devel-tools/compile-gromacs.sh index e95ead5aa..396fd7d3e 100755 --- a/devel-tools/compile-gromacs.sh +++ b/devel-tools/compile-gromacs.sh @@ -45,8 +45,8 @@ compile_gromacs_target() { while [ $# -ge 1 ]; do if [ "${1,,}" = "debug" ]; then GMX_BUILD_TYPE=Debug - GMX_BUILD_OPTS+=(-DCMAKE_BUILD_TYPE=Debug -DCMAKE_VERBOSE_MAKEFILE=yes) - GMX_BUILD_OPTS+=(-DCOLVARS_DEBUG=ON) + GMX_BUILD_OPTS+=(-DCMAKE_VERBOSE_MAKEFILE=yes) + GMX_BUILD_OPTS+=(-DCMAKE_CXX_FLAGS="-DCOLVARS_DEBUG=1") else GMX_INSTALL_DIR=${1} fi diff --git a/devel-tools/generate-namd-makefile.sh b/devel-tools/generate-namd-makefile.sh index d504a0d71..b9676f25b 100755 --- a/devel-tools/generate-namd-makefile.sh +++ b/devel-tools/generate-namd-makefile.sh @@ -1,9 +1,11 @@ #!/bin/bash +set -e + # Generate a file suitable for inclusion by the NAMD Makefile if [ -f "${2}" ] ; then - echo "Usage: ${0} [... source files ...]" >& 1 + echo "Usage: ${0} [... source files ...]" >& 2 exit 1 fi diff --git a/devel-tools/generate-vmd-makefile.sh b/devel-tools/generate-vmd-makefile.sh new file mode 100644 index 000000000..8b19547ec --- /dev/null +++ b/devel-tools/generate-vmd-makefile.sh @@ -0,0 +1,87 @@ +#!/bin/bash + +set -e + +# Generate the file colvars_files.pl used by VMD (i.e. a poor person's version of a dependency file) + +if [ ! -f "${1}" ] ; then + echo "Usage: ${0} [... source files ...] [... header files ...]" >& 2 + exit 1 +fi + +all_files=("${@}") + + +canonize_filename() { + local file=$1 + if [ "${file%.cpp}" != "${file}" ] ; then + basename "${file%.cpp}.C" + else + basename "${file}" + fi +} + +source_files=() +header_files=() +cuda_files=() +for file in "${all_files[@]}" ; do + if [ "${file%.cpp}" != "${file}" ] ; then + source_files+=($(canonize_filename "${file}")) + elif [ "${file%.C}" != "${file}" ] ; then + source_files+=($(canonize_filename "${file}")) + elif [ "${file%.h}" != "${file}" ] ; then + header_files+=($(canonize_filename "${file}")) + elif [ "${file%.cu}" != "${file}" ] ; then + echo "Error: building of Colvars CUDA files with VMD is not supported yet" >&2 + # cuda_files+=($(canonize_filename "${file}")) + else + echo "Error: file ${file} has an unsupported extension" >&2 + exit 1 + fi +done + +cat < 0$ ameliorates this to some degree, although every pair is still checked to regenerate the pair list. +This component returns a dimensionless number, which ranges from 0 (all interatomic distances are much larger than the cutoff) to $N_{\mathtt{group1}} \times N_{\mathtt{group2}}$ (all possible pairwise distances are much smaller than the cutoff). +The computational cost also grows with the number of pairwise distances calculated, unless either \texttt{group1CenterOnly} or \texttt{group2CenterOnly} are used. +Setting $\mathtt{tolerance} > 0$ could mitigate this by looping over pairs less frequently than every step. @@ -8201,6 +8212,9 @@ The legacy keyword \texttt{disableForces} for atom groups is now deprecated and will be discontinued in a future release. Atom groups now have an automated way to save computation if forces are not used, and enabling this option otherwise would lead to incorrect behavior. +\item \textbf{Colvars version 2026-XX-XX or later.}\\ + Periodic boundary conditions are handled internally, rather than using code by \MDENGINE{} to compute the minimum-image convention (see \ref{sec:cvc_pbcs}). + \end{itemize} diff --git a/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch b/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch index 00a64b7b8..b694bb766 100644 --- a/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch +++ b/gromacs/src/applied_forces/colvars/colvarproxygromacs.cpp.patch @@ -1,14 +1,24 @@ diff --git a/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp b/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp -index 76b350c611..37e5b53647 100644 +index 76b350c611..3876022dd3 100644 --- a/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp +++ b/src/gromacs/applied_forces/colvars/colvarproxygromacs.cpp @@ -79,9 +79,6 @@ ColvarProxyGromacs::ColvarProxyGromacs(const std::string& colvarsConfigString, // Retrieve masses and charges from input file updated_masses_ = updated_charges_ = true; - + - // User-scripted forces are not available in GROMACS - have_scripts = false; - angstrom_value_ = 0.1; - + // From Gnu units +@@ -281,7 +278,8 @@ ColvarProxyGromacs::~ColvarProxyGromacs() + } + + +-cvm::rvector ColvarProxyGromacs::position_distance(cvm::atom_pos const& pos1, cvm::atom_pos const& pos2) const ++cvm::rvector ColvarProxyGromacs::position_distance_engine(cvm::atom_pos const& pos1, ++ cvm::atom_pos const& pos2) const + { + rvec r1, r2, dr; + r1[0] = pos1.x; diff --git a/gromacs/src/applied_forces/colvars/colvarproxygromacs.h.diff b/gromacs/src/applied_forces/colvars/colvarproxygromacs.h.diff deleted file mode 100644 index 3a37a89be..000000000 --- a/gromacs/src/applied_forces/colvars/colvarproxygromacs.h.diff +++ /dev/null @@ -1,13 +0,0 @@ -diff --git a/src/gromacs/applied_forces/colvars/colvarproxygromacs.h b/src/gromacs/applied_forces/colvars/colvarproxygromacs.h -index 31ca7daf6e..32b1e8d0d7 100644 ---- a/src/gromacs/applied_forces/colvars/colvarproxygromacs.h -+++ b/src/gromacs/applied_forces/colvars/colvarproxygromacs.h -@@ -97,8 +97,6 @@ protected: - - - public: -- friend class cvm::atom; -- - /*! \brief Construct ColvarProxyGromacs from its parameters - * - * \param[in] colvarsConfigString Content of the colvars input file. diff --git a/gromacs/src/applied_forces/colvars/colvarproxygromacs.h.patch b/gromacs/src/applied_forces/colvars/colvarproxygromacs.h.patch new file mode 100644 index 000000000..7567c7c1b --- /dev/null +++ b/gromacs/src/applied_forces/colvars/colvarproxygromacs.h.patch @@ -0,0 +1,15 @@ +diff --git a/src/gromacs/applied_forces/colvars/colvarproxygromacs.h b/src/gromacs/applied_forces/colvars/colvarproxygromacs.h +index 32b1e8d0d7..3c0a0dee9a 100644 +--- a/src/gromacs/applied_forces/colvars/colvarproxygromacs.h ++++ b/src/gromacs/applied_forces/colvars/colvarproxygromacs.h +@@ -146,8 +146,8 @@ public: + */ + int check_atom_id(int atomNumber) override; + +- //! Compute the minimum distance with respect to the PBC between 2 atoms. +- cvm::rvector position_distance(cvm::atom_pos const& pos1, cvm::atom_pos const& pos2) const override; ++ //! Compute the minimum distance with respect to the PBCs between 2 positions using GROMACS code ++ cvm::rvector position_distance_engine(cvm::atom_pos const& pos1, cvm::atom_pos const& pos2) const override; + }; + + } // namespace gmx diff --git a/gromacs/src/applied_forces/colvars/colvarsMDModule.cpp.diff b/gromacs/src/applied_forces/colvars/colvarsMDModule.cpp.diff deleted file mode 100644 index 72ed7a57b..000000000 --- a/gromacs/src/applied_forces/colvars/colvarsMDModule.cpp.diff +++ /dev/null @@ -1,12 +0,0 @@ -diff --git a/src/gromacs/applied_forces/colvars/colvarsMDModule.cpp b/src/gromacs/applied_forces/colvars/colvarsMDModule.cpp -index 46df82f723..91e624c03c 100644 ---- a/src/gromacs/applied_forces/colvars/colvarsMDModule.cpp -+++ b/src/gromacs/applied_forces/colvars/colvarsMDModule.cpp -@@ -56,6 +56,7 @@ - #include "gromacs/utility/mpicomm.h" - #include "gromacs/utility/stringutil.h" - -+#include "colvars_version.h" - #include "colvarsforceprovider.h" - #include "colvarsoptions.h" - #include "colvarssimulationsparameters.h" diff --git a/gromacs/src/applied_forces/colvars/colvarsforceprovider.cpp.patch b/gromacs/src/applied_forces/colvars/colvarsforceprovider.cpp.patch new file mode 100644 index 000000000..6844d47c9 --- /dev/null +++ b/gromacs/src/applied_forces/colvars/colvarsforceprovider.cpp.patch @@ -0,0 +1,55 @@ +diff --git a/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp b/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp +index 5860eed790..bc34e6ccaa 100644 +--- a/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp ++++ b/src/gromacs/applied_forces/colvars/colvarsforceprovider.cpp +@@ -337,15 +337,47 @@ void ColvarsForceProvider::calculateForces(const ForceProviderInput& forceProvid + ForceProviderOutput* forceProviderOutput) + { + +- // Construct t_pbc struct +- set_pbc(&gmxPbc_, pbcType_, forceProviderInput.box_); ++ auto const& box = forceProviderInput.box_; ++ ++ // Set PBC struct for use in ColvarProxyGromacs::position_distance() ++ set_pbc(&gmxPbc_, pbcType_, box); ++ ++ // Share PBCs with the library as well ++ switch (pbcType_) ++ { ++ case PbcType::Unset: ++ GMX_RELEASE_ASSERT(false, "PBC type not set when being used by ColvarsForceProvider."); ++ break; ++ case PbcType::XY: ++ boundaries_.set_boundaries(true, ++ true, ++ false, ++ cvm::rvector{ box[0][0], box[0][1], box[0][2] }, ++ cvm::rvector{ box[1][0], box[1][1], box[1][2] }, ++ cvm::rvector{ 0.0, 0.0, 0.0 }); ++ break; ++ case PbcType::Xyz: ++ boundaries_.set_boundaries(true, ++ true, ++ true, ++ cvm::rvector{ box[0][0], box[0][1], box[0][2] }, ++ cvm::rvector{ box[1][0], box[1][1], box[1][2] }, ++ cvm::rvector{ box[2][0], box[2][1], box[2][2] }); ++ break; ++ case PbcType::Screw: ++ boundaries_.reset(); ++ boundaries_.set_type(cvm::system_boundary_conditions::types::unsupported); ++ break; ++ case PbcType::No: boundaries_.reset(); break; ++ default: ++ GMX_RELEASE_ASSERT(false, "Invalid pbcType in ColvarsForceProvider::calculateForces"); ++ } + + const MpiComm& mpiComm = forceProviderInput.mpiComm_; + // Local atom coords + const gmx::ArrayRef x = forceProviderInput.x_; + // Local atom coords (coerced into into old gmx type) + const rvec* xPointer = &(x.data()->as_vec()); +- const auto& box = forceProviderInput.box_; + + colvars->it = forceProviderInput.step_; + diff --git a/gromacs/update_patches_from_git.sh b/gromacs/update_patches_from_git.sh new file mode 100644 index 000000000..a60f7ba3d --- /dev/null +++ b/gromacs/update_patches_from_git.sh @@ -0,0 +1,21 @@ +#!/bin/bash + +GROMACS_SOURCE=${1} +if [ -n "${GROMACS_SOURCE}" ] ; then + GROMACS_SOURCE=$(git -C ${GROMACS_SOURCE} rev-parse --show-toplevel) +fi + +if [ $? != 0 ] || [ ! -s ${GROMACS_SOURCE}/src/gromacs/version.h.cmakein ] ; then + echo >&2 + echo "Usage: $0 " >& 2 + exit 1 +fi + +COLVARS_SOURCE=$(dirname $(realpath $0)) +COLVARS_SOURCE=$(git -C ${COLVARS_SOURCE} rev-parse --show-toplevel) + +modified_files=($(git -C ${GROMACS_SOURCE} ls-files -m src/gromacs/applied_forces/colvars)) + +for file in ${modified_files[@]} ; do + git -C ${GROMACS_SOURCE} diff ${file} > ${COLVARS_SOURCE}/gromacs/src/${file#src\/gromacs\/}.patch +done diff --git a/lammps/src/COLVARS/colvarproxy_lammps.cpp b/lammps/src/COLVARS/colvarproxy_lammps.cpp index be216b0bf..cc2b3b2c9 100644 --- a/lammps/src/COLVARS/colvarproxy_lammps.cpp +++ b/lammps/src/COLVARS/colvarproxy_lammps.cpp @@ -1,4 +1,3 @@ -// clang-format off // -*- c++ -*- // This file is part of the Collective Variables module (Colvars). @@ -133,24 +132,12 @@ double colvarproxy_lammps::compute() } previous_step = _lmp->update->ntimestep; - unit_cell_x.set(_lmp->domain->xprd, 0.0, 0.0); - unit_cell_y.set(0.0, _lmp->domain->yprd, 0.0); - unit_cell_z.set(0.0, 0.0, _lmp->domain->zprd); - - if (_lmp->domain->xperiodic == 0 && _lmp->domain->yperiodic == 0 && - _lmp->domain->zperiodic == 0) { - boundaries_type = boundaries_non_periodic; - reset_pbc_lattice(); - } else if ((_lmp->domain->nonperiodic == 0) && - (_lmp->domain->dimension == 3) && - (_lmp->domain->triclinic == 0)) { - // Orthogonal unit cell - boundaries_type = boundaries_pbc_ortho; - colvarproxy_system::update_pbc_lattice(); - // It is safer to let LAMMPS deal with high-tilt triclinic boxes - } else { - boundaries_type = boundaries_unsupported; - } + boundaries_.set_boundaries(_lmp->domain->xperiodic, + _lmp->domain->yperiodic, + _lmp->domain->zperiodic, + cvm::rvector{_lmp->domain->xprd, 0.0, 0.0}, + cvm::rvector{_lmp->domain->xy, _lmp->domain->yprd, 0.0}, + cvm::rvector{_lmp->domain->xz, _lmp->domain->yz, _lmp->domain->zprd}); if (cvm::debug()) { cvm::log(std::string(cvm::line_marker) + @@ -187,14 +174,13 @@ double colvarproxy_lammps::compute() /* ---------------------------------------------------------------------- */ -cvm::rvector colvarproxy_lammps::position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) - const +cvm::rvector colvarproxy_lammps::position_distance_engine(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const { double xtmp = pos2.x - pos1.x; double ytmp = pos2.y - pos1.y; double ztmp = pos2.z - pos1.z; - _lmp->domain->minimum_image_big(FLERR, xtmp,ytmp,ztmp); + _lmp->domain->minimum_image_big(FLERR, xtmp, ytmp, ztmp); return {xtmp, ytmp, ztmp}; } diff --git a/lammps/src/COLVARS/colvarproxy_lammps.h b/lammps/src/COLVARS/colvarproxy_lammps.h index 9acdfc343..a11f81a35 100644 --- a/lammps/src/COLVARS/colvarproxy_lammps.h +++ b/lammps/src/COLVARS/colvarproxy_lammps.h @@ -79,8 +79,8 @@ class colvarproxy_lammps : public colvarproxy { void log(std::string const &message) override; void error(std::string const &message) override; - [[nodiscard]] cvm::rvector position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) const override; + [[nodiscard]] cvm::rvector position_distance_engine(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const override; cvm::real rand_gaussian() override; diff --git a/misc_interfaces/stubs/colvarproxy_stub.cpp b/misc_interfaces/stubs/colvarproxy_stub.cpp index df2daabd7..d6b02bfb6 100644 --- a/misc_interfaces/stubs/colvarproxy_stub.cpp +++ b/misc_interfaces/stubs/colvarproxy_stub.cpp @@ -24,6 +24,9 @@ colvarproxy_stub::colvarproxy_stub() { version_int = get_version_from_string(COLVARPROXY_VERSION); + + use_internal_pbc_ = true; + b_simulation_running = false; // both fields are taken from data structures already available @@ -49,8 +52,6 @@ colvarproxy_stub::~colvarproxy_stub() int colvarproxy_stub::setup() { - boundaries_type = boundaries_non_periodic; - reset_pbc_lattice(); colvars->it = colvars->it_restart = 0; if (colvars) { diff --git a/namd/src/colvarproxy_namd.C b/namd/src/colvarproxy_namd.C index 7f6d0f570..5877d5a50 100644 --- a/namd/src/colvarproxy_namd.C +++ b/namd/src/colvarproxy_namd.C @@ -388,29 +388,10 @@ void colvarproxy_namd::calculate() if (accelMDOn) update_accelMD_info(); auto *lattice = globalmaster->get_lattice(); - - { - Vector const a = lattice->a(); - Vector const b = lattice->b(); - Vector const c = lattice->c(); - unit_cell_x.set(a.x, a.y, a.z); - unit_cell_y.set(b.x, b.y, b.z); - unit_cell_z.set(c.x, c.y, c.z); - } - - if (!lattice->a_p() && !lattice->b_p() && !lattice->c_p()) { - boundaries_type = boundaries_non_periodic; - reset_pbc_lattice(); - } else if (lattice->a_p() && lattice->b_p() && lattice->c_p()) { - if (lattice->orthogonal()) { - boundaries_type = boundaries_pbc_ortho; - } else { - boundaries_type = boundaries_pbc_triclinic; - } - colvarproxy_system::update_pbc_lattice(); - } else { - boundaries_type = boundaries_unsupported; - } + boundaries_.set_boundaries(lattice->a_p(), lattice->b_p(), lattice->c_p(), + cvm::rvector{lattice->a().x, lattice->a().y, lattice->a().z}, + cvm::rvector{lattice->b().x, lattice->b().y, lattice->b().z}, + cvm::rvector{lattice->c().x, lattice->c().y, lattice->c().z}); if (cvm::debug()) { cvm::log(std::string(cvm::line_marker)+ @@ -889,19 +870,17 @@ void colvarproxy_namd::update_atom_properties(int index) } -cvm::rvector colvarproxy_namd::position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) - const +cvm::rvector colvarproxy_namd::position_distance_engine(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const { - Position const p1(pos1.x, pos1.y, pos1.z); - Position const p2(pos2.x, pos2.y, pos2.z); + Position const p1{pos1.x, pos1.y, pos1.z}; + Position const p2{pos2.x, pos2.y, pos2.z}; // return p2 - p1 Vector const d = globalmaster->get_lattice()->delta(p2, p1); - return cvm::rvector(d.x, d.y, d.z); + return {d.x, d.y, d.z}; } - enum e_pdb_field { e_pdb_none, e_pdb_occ, diff --git a/namd/src/colvarproxy_namd.h b/namd/src/colvarproxy_namd.h index 57f114b58..b58589b8e 100644 --- a/namd/src/colvarproxy_namd.h +++ b/namd/src/colvarproxy_namd.h @@ -203,8 +203,8 @@ class colvarproxy_namd : public colvarproxy { void update_atom_properties(int index); - cvm::rvector position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) const; + cvm::rvector position_distance_engine(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const override; int load_atoms_pdb(char const *filename, cvm::atom_group &atoms, diff --git a/src/.clang-format b/src/.clang-format index 66e78c8b5..e3f508c76 100644 --- a/src/.clang-format +++ b/src/.clang-format @@ -1,14 +1,16 @@ --- -Language: Cpp -ColumnLimit: 100 -UseTab: Never -TabWidth: 2 -IndentWidth: 2 -ObjCBlockIndentWidth: 2 -PenaltyBreakAssignment: 2 AccessModifierOffset: -2 +BinPackParameters: false BreakBeforeBraces: WebKit +ColumnLimit: 100 EmptyLineAfterAccessModifier: Leave -NamespaceIndentation: None +IndentWidth: 2 +Language: Cpp MaxEmptyLinesToKeep: 2 +NamespaceIndentation: None +ObjCBlockIndentWidth: 2 +PenaltyBreakAssignment: 2 +PenaltyReturnTypeOnItsOwnLine: 10000 +TabWidth: 2 +UseTab: Never ... diff --git a/src/colvar.cpp b/src/colvar.cpp index 90f288112..48979e161 100644 --- a/src/colvar.cpp +++ b/src/colvar.cpp @@ -21,6 +21,7 @@ #include "colvarbias.h" #include "colvars_memstream.h" #include "colvarcomp_torchann.h" +#include "colvarcomp_coordnums.h" std::map> colvar::global_cvc_map = std::map>(); diff --git a/src/colvarcomp.cpp b/src/colvarcomp.cpp index 53854799b..666a39c9a 100644 --- a/src/colvarcomp.cpp +++ b/src/colvarcomp.cpp @@ -446,6 +446,14 @@ int colvar::cvc::set_param(std::string const ¶m_name, void colvar::cvc::read_data() { + if (!is_enabled(f_cvc_pbc_minimum_image)) { + // Copy boundary conditions from the proxy + boundary_conditions = cvm::main()->proxy->get_system_boundaries(); + } else { + // Set as non-periodic boundary conditions (default) for this CVC + boundary_conditions.reset(); + } + if (is_enabled(f_cvc_explicit_atom_groups)) { for (auto agi = atom_groups.begin(); agi != atom_groups.end(); agi++) { auto &atoms = *(*agi); diff --git a/src/colvarcomp.h b/src/colvarcomp.h index 2f378effa..528e790df 100644 --- a/src/colvarcomp.h +++ b/src/colvarcomp.h @@ -347,6 +347,9 @@ class colvar::cvc /// \brief CVC-specific default colvar width (default: not provided) cvm::real width = 0.0; + + /// Boundary conditions for the system, copied from the engine when needed + cvm::system_boundary_conditions boundary_conditions; }; @@ -795,180 +798,6 @@ class colvar::dihedral }; - -/// \brief Colvar component: coordination number between two groups -/// (colvarvalue::type_scalar type, range [0:N1*N2]) -class colvar::coordnum - : public colvar::cvc -{ -protected: - /// First atom group - cvm::atom_group *group1 = nullptr; - /// Second atom group - cvm::atom_group *group2 = nullptr; - /// \brief "Cutoff" for isotropic calculation (default) - cvm::real r0; - /// \brief "Cutoff vector" for anisotropic calculation - cvm::rvector r0_vec; - /// \brief Whether r/r0 or \vec{r}*\vec{1/r0_vec} should be used - bool b_anisotropic = false; - /// Integer exponent of the function numerator - int en = 6; - /// Integer exponent of the function denominator - int ed = 12; - - /// If true, group2 will be treated as a single atom - bool b_group2_center_only = false; - - /// Tolerance for the pair list - cvm::real tolerance = 0.0; - - /// Frequency of update of the pair list - int pairlist_freq = 100; - - /// Pair list - bool *pairlist = nullptr; - -public: - - coordnum(); - virtual ~coordnum(); - virtual int init(std::string const &conf); - virtual void calc_value(); - virtual void calc_gradients(); - - enum { - ef_null = 0, - ef_gradients = 1, - ef_anisotropic = (1<<8), - ef_use_pairlist = (1<<9), - ef_rebuild_pairlist = (1<<10) - }; - - /// \brief Calculate a coordination number through the function - /// (1-x**n)/(1-x**m), where x = |A1-A2|/r0 \param r0, r0_vec "cutoff" for - /// the coordination number (scalar or vector depending on user choice) - /// \param en Numerator exponent \param ed Denominator exponent \param First - /// atom \param Second atom \param pairlist_elem pointer to pair flag for - /// this pair \param tolerance A pair is defined as having a larger - /// coordination than this number - template - static cvm::real switching_function(cvm::real const &r0, - cvm::rvector const &inv_r0_vec, - cvm::rvector const &inv_r0sq_vec, - int en, - int ed, - const cvm::real a1x, - const cvm::real a1y, - const cvm::real a1z, - const cvm::real a2x, - const cvm::real a2y, - const cvm::real a2z, - cvm::real& g1x, - cvm::real& g1y, - cvm::real& g1z, - cvm::real& g2x, - cvm::real& g2y, - cvm::real& g2z, - bool **pairlist_elem, - cvm::real tolerance); - - /// Workhorse function - template int compute_coordnum(); - - /// Workhorse function - template void main_loop(bool **pairlist_elem); - -}; - - - -/// \brief Colvar component: self-coordination number within a group -/// (colvarvalue::type_scalar type, range [0:N*(N-1)/2]) -class colvar::selfcoordnum - : public colvar::cvc -{ -protected: - /// Selected atoms - cvm::atom_group *group1 = nullptr; - /// \brief "Cutoff" for isotropic calculation (default) - cvm::real r0; - /// Integer exponent of the function numerator - int en = 6; - /// Integer exponent of the function denominator - int ed = 12; - cvm::real tolerance = 0.0; - int pairlist_freq = 100; - - bool *pairlist = nullptr; - -public: - - selfcoordnum(); - ~selfcoordnum(); - virtual int init(std::string const &conf); - virtual void calc_value(); - virtual void calc_gradients(); - - /// Main workhorse function - template int compute_selfcoordnum(); -}; - - - -/// \brief Colvar component: coordination number between two groups -/// (colvarvalue::type_scalar type, range [0:N1*N2]) -class colvar::groupcoordnum - : public colvar::distance -{ -protected: - /// \brief "Cutoff" for isotropic calculation (default) - cvm::real r0; - /// \brief "Cutoff vector" for anisotropic calculation - cvm::rvector r0_vec; - /// \brief Wheter dist/r0 or \vec{dist}*\vec{1/r0_vec} should ne be - /// used - bool b_anisotropic = false; - /// Integer exponent of the function numerator - int en = 6; - /// Integer exponent of the function denominator - int ed = 12; -public: - groupcoordnum(); - virtual ~groupcoordnum() {} - virtual int init(std::string const &conf); - virtual void calc_value(); - virtual void calc_gradients(); -}; - - - -/// \brief Colvar component: hydrogen bond, defined as the product of -/// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol)) -/// (colvarvalue::type_scalar type, range [0:1]) -class colvar::h_bond - : public colvar::cvc -{ -protected: - /// \brief "Cutoff" distance between acceptor and donor - cvm::real r0; - /// Integer exponent of the function numerator - int en = 6; - /// Integer exponent of the function denominator - int ed = 8; -public: - /// Constructor for atoms already allocated - h_bond(cvm::atom_group::simple_atom const &acceptor, - cvm::atom_group::simple_atom const &donor, - cvm::real r0, int en, int ed); - h_bond(); - virtual ~h_bond() {} - virtual int init(std::string const &conf); - virtual void calc_value(); - virtual void calc_gradients(); -}; - - /// \brief Colvar component: alpha helix content of a contiguous /// segment of 5 or more residues, implemented as a sum of Ca-Ca-Ca /// angles and hydrogen bonds (colvarvalue::type_scalar type, range diff --git a/src/colvarcomp_angles.cpp b/src/colvarcomp_angles.cpp index ddd960721..2f32fc239 100644 --- a/src/colvarcomp_angles.cpp +++ b/src/colvarcomp_angles.cpp @@ -63,13 +63,9 @@ void colvar::angle::calc_value() cvm::atom_pos const g2_pos = group2->center_of_mass(); cvm::atom_pos const g3_pos = group3->center_of_mass(); - r21 = is_enabled(f_cvc_pbc_minimum_image) ? - cvm::position_distance(g2_pos, g1_pos) : - g1_pos - g2_pos; + r21 = boundary_conditions.position_distance(g2_pos, g1_pos); r21l = r21.norm(); - r23 = is_enabled(f_cvc_pbc_minimum_image) ? - cvm::position_distance(g2_pos, g3_pos) : - g3_pos - g2_pos; + r23 = boundary_conditions.position_distance(g2_pos, g3_pos); r23l = r23.norm(); cvm::real const cos_theta = (r21*r23)/(r21l*r23l); @@ -158,9 +154,7 @@ void colvar::dipole_angle::calc_value() r21 = group1->dipole(); r21l = r21.norm(); - r23 = is_enabled(f_cvc_pbc_minimum_image) ? - cvm::position_distance(g2_pos, g3_pos) : - g3_pos - g2_pos; + r23 = boundary_conditions.position_distance(g2_pos, g3_pos); r23l = r23.norm(); cvm::real denom = (r21l*r23l); @@ -281,15 +275,9 @@ void colvar::dihedral::calc_value() cvm::atom_pos const g4_pos = group4->center_of_mass(); // Usual sign convention: r12 = r2 - r1 - r12 = is_enabled(f_cvc_pbc_minimum_image) ? - cvm::position_distance(g1_pos, g2_pos) : - g2_pos - g1_pos; - r23 = is_enabled(f_cvc_pbc_minimum_image) ? - cvm::position_distance(g2_pos, g3_pos) : - g3_pos - g2_pos; - r34 = is_enabled(f_cvc_pbc_minimum_image) ? - cvm::position_distance(g3_pos, g4_pos) : - g4_pos - g3_pos; + r12 = boundary_conditions.position_distance(g1_pos, g2_pos); + r23 = boundary_conditions.position_distance(g2_pos, g3_pos); + r34 = boundary_conditions.position_distance(g3_pos, g4_pos); cvm::rvector const n1 = cvm::rvector::outer(r12, r23); cvm::rvector const n2 = cvm::rvector::outer(r23, r34); diff --git a/src/colvarcomp_coordnums.cpp b/src/colvarcomp_coordnums.cpp index 99d3f3e69..410cc38f9 100644 --- a/src/colvarcomp_coordnums.cpp +++ b/src/colvarcomp_coordnums.cpp @@ -12,124 +12,34 @@ #include "colvarvalue.h" #include "colvar.h" #include "colvarcomp.h" - -template -inline -cvm::real colvar::coordnum::switching_function(cvm::real const &r0, - cvm::rvector const &inv_r0_vec, - cvm::rvector const &inv_r0sq_vec, - int en, - int ed, - const cvm::real a1x, - const cvm::real a1y, - const cvm::real a1z, - const cvm::real a2x, - const cvm::real a2y, - const cvm::real a2z, - cvm::real& g1x, - cvm::real& g1y, - cvm::real& g1z, - cvm::real& g2x, - cvm::real& g2y, - cvm::real& g2z, - bool **pairlist_elem, - cvm::real pairlist_tol) -{ - if ((flags & ef_use_pairlist) && !(flags & ef_rebuild_pairlist)) { - bool const within = **pairlist_elem; - (*pairlist_elem)++; - if (!within) { - return 0.0; - } - } - - const cvm::atom_pos pos1{a1x, a1y, a1z}; - const cvm::atom_pos pos2{a2x, a2y, a2z}; - cvm::rvector const diff = cvm::position_distance(pos1, pos2); - cvm::rvector const scal_diff(diff.x * inv_r0_vec.x, - diff.y * inv_r0_vec.y, - diff.z * inv_r0_vec.z); - cvm::real const l2 = scal_diff.norm2(); - - // Assume en and ed are even integers, and avoid sqrt in the following - int const en2 = en/2; - int const ed2 = ed/2; - - cvm::real const xn = cvm::integer_power(l2, en2); - cvm::real const xd = cvm::integer_power(l2, ed2); - cvm::real const eps_l2 = 1.0e-7; - cvm::real const h = l2 - 1.0; - cvm::real const en2_r = (cvm::real) en2; - cvm::real const ed2_r = (cvm::real) ed2; - cvm::real func_no_pairlist; - - if (std::abs(h) < eps_l2) { - // Order-2 Taylor expansion: c0 + c1*h + c2*h^2 - cvm::real const c0 = en2_r / ed2_r; - cvm::real const c1 = (en2_r * (en2_r - ed2_r)) / (2.0 * ed2_r); - cvm::real const c2 = (en2_r * (en2_r - ed2_r) * (2.0 * en2_r - ed2_r - 3.0)) / (12.0 * ed2_r); - func_no_pairlist = c0 + h * (c1 + h * c2); - } else { - func_no_pairlist = (1.0 - xn) / (1.0 - xd); - } - - cvm::real func, inv_one_pairlist_tol; - if (flags & ef_use_pairlist) { - inv_one_pairlist_tol = 1 / (1.0-pairlist_tol); - func = (func_no_pairlist - pairlist_tol) * inv_one_pairlist_tol; - } else { - func = func_no_pairlist; - } - - if (flags & ef_rebuild_pairlist) { - //Particles just outside of the cutoff also are considered if they come near. - **pairlist_elem = (func > (-pairlist_tol * 0.5)) ? true : false; - (*pairlist_elem)++; - } - //If the value is too small, we need to exclude it, rather than let it contribute to the sum or the gradients. - if (func < 0) - return 0; - - if (flags & ef_gradients) { - // Logarithmic derivative: 1st-order Taylor expansion around l2 = 1 - cvm::real log_deriv; - if (std::abs(h) < eps_l2) { - cvm::real const g0 = 0.5 * (en2_r - ed2_r); - cvm::real const g1 = ((en2_r - ed2_r) * (en2_r + ed2_r - 6.0)) / 12.0; - log_deriv = g0 + h * g1; - } else { - log_deriv = (ed2_r * xd / ((1.0 - xd) * l2)) - (en2_r * xn / ((1.0 - xn) * l2)); - } - cvm::real const dFdl2 = (flags & ef_use_pairlist) ? - func_no_pairlist * inv_one_pairlist_tol * log_deriv : - func * log_deriv; - - cvm::rvector const dl2dx((2.0 * inv_r0sq_vec.x) * diff.x, - (2.0 * inv_r0sq_vec.y) * diff.y, - (2.0 * inv_r0sq_vec.z) * diff.z); - - const cvm::rvector G = dFdl2*dl2dx; - g1x += -1.0*G.x; - g1y += -1.0*G.y; - g1z += -1.0*G.z; - g2x += G.x; - g2y += G.y; - g2z += G.z; - } - - return func; -} +#include "colvarcomp_coordnums.h" colvar::coordnum::coordnum() { set_function_type("coordNum"); x.type(colvarvalue::type_scalar); - colvarproxy *proxy = cvm::main()->proxy; - r0 = proxy->angstrom_to_internal(4.0); - r0_vec = cvm::rvector(proxy->angstrom_to_internal(4.0), - proxy->angstrom_to_internal(4.0), - proxy->angstrom_to_internal(4.0)); + cvm::real const r0 = cvm::main()->proxy->angstrom_to_internal(4.0); + update_cutoffs({r0, r0, r0}); + // Boundaries will be set later, when the number of pairs is known +} + + +void colvar::coordnum::update_cutoffs(cvm::rvector const &r0_vec_i) +{ + r0_vec = r0_vec_i; + + inv_r0_vec = { + 1.0 / r0_vec.x, + 1.0 / r0_vec.y, + 1.0 / r0_vec.z + }; + + inv_r0sq_vec = { + inv_r0_vec.x * inv_r0_vec.x, + inv_r0_vec.y * inv_r0_vec.y, + inv_r0_vec.z * inv_r0_vec.z + }; } @@ -138,37 +48,74 @@ int colvar::coordnum::init(std::string const &conf) int error_code = cvc::init(conf); group1 = parse_group(conf, "group1"); - group2 = parse_group(conf, "group2"); - if (!group1 || !group2) { + if (!group1) { return error_code | COLVARS_INPUT_ERROR; } - if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) { - error_code |= cvm::error( - "Error: group1 and group2 share a common atom (number: " + cvm::to_str(atom_number) + ")\n", - COLVARS_INPUT_ERROR); + if (group1->b_dummy) { + error_code |= cvm::error("Error: group1 may not be a dummy atom\n", COLVARS_INPUT_ERROR); } - if (group1->b_dummy) { - error_code |= - cvm::error("Error: only group2 is allowed to be a dummy atom\n", COLVARS_INPUT_ERROR); + if (function_type() != "selfCoordNum") { + + group2 = parse_group(conf, "group2"); + if (!group2) { + return error_code | COLVARS_INPUT_ERROR; + } + + if (int atom_number = cvm::atom_group::overlap(*group1, *group2)) { + error_code |= cvm::error("Error: group1 and group2 share a common atom (number: " + + cvm::to_str(atom_number) + ")\n", + COLVARS_INPUT_ERROR); + } + + if (function_type() == "coordNum") { + get_keyval(conf, "group1CenterOnly", b_group1_center_only, group2->b_dummy); + get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy); + } + + if (function_type() == "groupCoord") { + // In groupCoord, these flags are hard-coded + b_group1_center_only = true; + b_group2_center_only = true; + } + + size_t const group1_num_coords = b_group1_center_only ? 1 : group1->size(); + size_t const group2_num_coords = b_group2_center_only ? 1 : group2->size(); + + num_pairs = group1_num_coords * group2_num_coords; + + } else { + + // selfCoordNum case + num_pairs = (group1->size() * (group1->size() - 1)) / 2; } - bool const b_isotropic = get_keyval(conf, "cutoff", r0, r0); + init_scalar_boundaries(0.0, num_pairs); + + // Get the default value from r0_vec to report it + cvm::real r0 = r0_vec[0]; + bool const b_redefined_cutoff = get_keyval(conf, "cutoff", r0, r0); if (get_keyval(conf, "cutoff3", r0_vec, r0_vec)) { - if (b_isotropic) { - error_code |= cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" " - "at the same time.\n", - COLVARS_INPUT_ERROR); + if (b_redefined_cutoff) { + error_code |= + cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n", + COLVARS_INPUT_ERROR); } - b_anisotropic = true; // remove meaningless negative signs if (r0_vec.x < 0.0) r0_vec.x *= -1.0; if (r0_vec.y < 0.0) r0_vec.y *= -1.0; if (r0_vec.z < 0.0) r0_vec.z *= -1.0; + + update_cutoffs(r0_vec); + + } else { + if (b_redefined_cutoff) { + update_cutoffs({r0, r0, r0}); + } } get_keyval(conf, "expNumer", en, en); @@ -184,150 +131,138 @@ int colvar::coordnum::init(std::string const &conf) COLVARS_INPUT_ERROR); } - if (!is_enabled(f_cvc_pbc_minimum_image)) { - cvm::log("Warning: only minimum-image distances are used by this variable.\n"); - } - - get_keyval(conf, "group2CenterOnly", b_group2_center_only, group2->b_dummy); - - get_keyval(conf, "tolerance", tolerance, tolerance); - if (tolerance > 0) { - cvm::main()->cite_feature("coordNum pairlist"); - get_keyval(conf, "pairListFrequency", pairlist_freq, pairlist_freq); - if ( ! (pairlist_freq > 0) ) { - return cvm::error("Error: non-positive pairlistfrequency provided.\n", - COLVARS_INPUT_ERROR); - // return and do not allocate the pairlists below - } - size_t n; - if (b_group2_center_only) { - n = group1->size(); - } else { - n = group1->size() * group2->size(); + if (function_type() != "groupCoord") { + // All coordNum variables may benefit from a pairlist, except groupCoord + get_keyval(conf, "tolerance", tolerance, tolerance); + if (tolerance > 0) { + cvm::main()->cite_feature("coordNum pairlist"); + compute_tolerance_l2_max(); + get_keyval(conf, "pairListFrequency", pairlist_freq, pairlist_freq); + if ( ! (pairlist_freq > 0) ) { + return cvm::error("Error: non-positive pairlistfrequency provided.\n", + COLVARS_INPUT_ERROR); + // return and do not allocate the pairlists below + } + pairlist.reset(new bool[num_pairs]); + auto *pairlist_elem = pairlist.get(); + for (size_t ip = 0; ip < num_pairs; ip++, pairlist_elem++) { + *pairlist_elem = true; + } } - pairlist = new bool[n]; - for (size_t i = 0; i < n; i++) pairlist[i] = true; } - init_scalar_boundaries(0.0, b_group2_center_only ? - static_cast(group1->size()) : - static_cast(group1->size() * - group2->size())); - return error_code; } -colvar::coordnum::~coordnum() -{ - if (pairlist) { - delete [] pairlist; - } -} +colvar::coordnum::~coordnum() {} -template void colvar::coordnum::main_loop(bool **pairlist_elem) +void colvar::coordnum::compute_tolerance_l2_max() { - const cvm::rvector inv_r0_vec( - 1.0 / ((flags & ef_anisotropic) ? r0_vec.x : r0), - 1.0 / ((flags & ef_anisotropic) ? r0_vec.y : r0), - 1.0 / ((flags & ef_anisotropic) ? r0_vec.z : r0)); - cvm::rvector const inv_r0sq_vec( - inv_r0_vec.x*inv_r0_vec.x, - inv_r0_vec.y*inv_r0_vec.y, - inv_r0_vec.z*inv_r0_vec.z); - if (b_group2_center_only) { - const cvm::atom_pos group2_com = group2->center_of_mass(); - cvm::rvector group2_com_grad(0, 0, 0); - for (size_t i = 0; i < group1->size(); ++i) { - x.real_value += switching_function(r0, inv_r0_vec, - inv_r0sq_vec, en, ed, - group1->pos_x(i), - group1->pos_y(i), - group1->pos_z(i), - group2_com.x, - group2_com.y, - group2_com.z, - group1->grad_x(i), - group1->grad_y(i), - group1->grad_z(i), - group2_com_grad.x, - group2_com_grad.y, - group2_com_grad.z, - pairlist_elem, - tolerance); - } - if (b_group2_center_only) { - group2->set_weighted_gradient(group2_com_grad); - } - } else { - for (size_t i = 0; i < group1->size(); ++i) { - for (size_t j = 0; j < group2->size(); ++j) { - x.real_value += switching_function(r0, inv_r0_vec, - inv_r0sq_vec, en, ed, - group1->pos_x(i), - group1->pos_y(i), - group1->pos_z(i), - group2->pos_x(j), - group2->pos_y(j), - group2->pos_z(j), - group1->grad_x(i), - group1->grad_y(i), - group1->grad_z(i), - group2->grad_x(j), - group2->grad_y(j), - group2->grad_z(j), - pairlist_elem, - tolerance); - } + cvm::real l2 = 1.001; + cvm::real F = 0.0; + cvm::real dFdl2 = 0.0; + constexpr size_t num_iters_max = 1000000; + constexpr cvm::real result_tol = 1.0e-6; + constexpr cvm::real dF_tol = 1.0e-9; + size_t i; + // Find the value of l2 such that F(l2) = 0 using the Newton method + for (i = 0; i < num_iters_max; i++) { + F = switching_function(l2, dFdl2, en, ed, tolerance); + if ((std::fabs(F) < result_tol) || (std::fabs(dFdl2) < dF_tol)) { + break; } + l2 -= F / dFdl2; + } + tolerance_l2_max = l2; + if (cvm::debug()) { + cvm::log("Found max valid l2 in " + cvm::to_str(i+1) + " iterations, result = " + cvm::to_str(l2) + " f(result) = " + cvm::to_str(F)); } } -template int colvar::coordnum::compute_coordnum() +template +void inline colvar::coordnum::main_loop() { - bool const use_pairlist = (pairlist != NULL); - bool const rebuild_pairlist = (pairlist != NULL) && - (cvm::step_relative() % pairlist_freq == 0); + size_t const group1_num_coords = use_group1_com ? 1 : group1->size(); + size_t const group2_num_coords = use_group2_com ? 1 : group2->size(); - bool *pairlist_elem = use_pairlist ? pairlist : NULL; + cvm::atom_pos const group1_com = group1->center_of_mass(); + cvm::atom_pos const group2_com = group2->center_of_mass(); + cvm::rvector group1_com_grad, group2_com_grad; - if (b_anisotropic) { + bool *pairlist_elem = pairlist.get(); - if (use_pairlist) { - if (rebuild_pairlist) { - int const flags = compute_flags | ef_anisotropic | ef_use_pairlist | - ef_rebuild_pairlist; - main_loop(&pairlist_elem); - } else { - int const flags = compute_flags | ef_anisotropic | ef_use_pairlist; - main_loop(&pairlist_elem); + for (size_t i = 0; i < group1_num_coords; ++i) { + + cvm::real const x1 = use_group1_com ? group1_com.x : group1->pos_x(i); + cvm::real const y1 = use_group1_com ? group1_com.y : group1->pos_y(i); + cvm::real const z1 = use_group1_com ? group1_com.z : group1->pos_z(i); + + cvm::real &gx1 = use_group1_com ? group1_com_grad.x : group1->grad_x(i); + cvm::real &gy1 = use_group1_com ? group1_com_grad.y : group1->grad_y(i); + cvm::real &gz1 = use_group1_com ? group1_com_grad.z : group1->grad_z(i); + + for (size_t j = 0; j < group2_num_coords; ++j) { + + cvm::real const x2 = use_group2_com ? group2_com.x : group2->pos_x(j); + cvm::real const y2 = use_group2_com ? group2_com.y : group2->pos_y(j); + cvm::real const z2 = use_group2_com ? group2_com.z : group2->pos_z(j); + + cvm::real &gx2 = use_group2_com ? group2_com_grad.x : group2->grad_x(j); + cvm::real &gy2 = use_group2_com ? group2_com_grad.y : group2->grad_y(j); + cvm::real &gz2 = use_group2_com ? group2_com_grad.z : group2->grad_z(j); + + bool const within = + ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) || + !(flags & ef_use_pairlist); + + cvm::real const partial = + within ? compute_pair_coordnum(inv_r0_vec, inv_r0sq_vec, en, ed, + x1, y1, z1, x2, y2, z2, + gx1, gy1, gz1, gx2, gy2, gz2, + tolerance, tolerance_l2_max, boundary_conditions) + : 0.0; + + if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) { + *pairlist_elem = partial > 0.0 ? true : false; } - } else { + x.real_value += partial; - int const flags = compute_flags | ef_anisotropic; - main_loop(NULL); + if (flags & ef_use_pairlist) { + pairlist_elem++; + } } + } - } else { + if (use_group1_com) { + group1->set_weighted_gradient(group1_com_grad); + } + if (use_group2_com) { + group2->set_weighted_gradient(group2_com_grad); + } +} - if (use_pairlist) { - if (rebuild_pairlist) { - int const flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist; - main_loop(&pairlist_elem); - } else { - int const flags = compute_flags | ef_use_pairlist; - main_loop(&pairlist_elem); - } +template +int colvar::coordnum::compute_coordnum() +{ + bool const use_pairlist = pairlist.get(); + bool const rebuild_pairlist = use_pairlist && (cvm::step_relative() % pairlist_freq == 0); + if (use_pairlist) { + if (rebuild_pairlist) { + constexpr int flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist; + main_loop(); } else { - - int const flags = compute_flags; - main_loop(NULL); + constexpr int flags = compute_flags | ef_use_pairlist; + main_loop(); } + } else { + constexpr int flags = compute_flags; + main_loop(); } return COLVARS_OK; @@ -338,9 +273,40 @@ void colvar::coordnum::calc_value() { x.real_value = 0.0; if (is_enabled(f_cvc_gradient)) { - compute_coordnum(); + + constexpr int flags = ef_gradients; + + if (b_group1_center_only) { + if (b_group2_center_only) { + compute_coordnum(); + } else { + compute_coordnum(); + } + } else { + if (b_group2_center_only) { + compute_coordnum(); + } else { + compute_coordnum(); + } + } + } else { - compute_coordnum(); + + constexpr int flags = ef_null; + + if (b_group1_center_only) { + if (b_group2_center_only) { + compute_coordnum(); + } else { + compute_coordnum(); + } + } else { + if (b_group2_center_only) { + compute_coordnum(); + } else { + compute_coordnum(); + } + } } } @@ -356,8 +322,8 @@ void colvar::coordnum::calc_gradients() colvar::h_bond::h_bond() { - colvarproxy *proxy = cvm::main()->proxy; - r0 = proxy->angstrom_to_internal(3.3); + cvm::real const r0 = cvm::main()->proxy->angstrom_to_internal(3.3); + r0_vec = {r0, r0, r0}; set_function_type("hBond"); x.type(colvarvalue::type_scalar); init_scalar_boundaries(0.0, 1.0); @@ -391,7 +357,11 @@ int colvar::h_bond::init(std::string const &conf) modify_atom.add_atom(cvm::atom_group::init_atom_from_proxy(p, d_num)); } - get_keyval(conf, "cutoff", r0, r0); + cvm::real r0 = r0_vec[0]; + bool const b_redefined_cutoff = get_keyval(conf, "cutoff", r0, r0); + if (b_redefined_cutoff) { + r0_vec = {r0, r0, r0}; + } get_keyval(conf, "expNumer", en, en); get_keyval(conf, "expDenom", ed, ed); @@ -415,7 +385,7 @@ colvar::h_bond::h_bond(cvm::atom_group::simple_atom const &acceptor, cvm::real r0_i, int en_i, int ed_i) : h_bond() { - r0 = r0_i; + r0_vec = {r0_i, r0_i, r0_i}; en = en_i; ed = ed_i; register_atom_group(new cvm::atom_group); @@ -427,8 +397,7 @@ colvar::h_bond::h_bond(cvm::atom_group::simple_atom const &acceptor, void colvar::h_bond::calc_value() { - int const flags = coordnum::ef_null; - cvm::rvector const r0_vec(0.0); // TODO enable the flag? + constexpr int flags = coordnum::ef_null; cvm::rvector G1, G2; const cvm::atom_pos A1{atom_groups[0]->pos_x(0), atom_groups[0]->pos_y(0), @@ -436,184 +405,139 @@ void colvar::h_bond::calc_value() const cvm::atom_pos A2{atom_groups[0]->pos_x(1), atom_groups[0]->pos_y(1), atom_groups[0]->pos_z(1)}; - const cvm::rvector inv_r0_vec( - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); - cvm::rvector const inv_r0sq_vec( - inv_r0_vec.x*inv_r0_vec.x, - inv_r0_vec.y*inv_r0_vec.y, - inv_r0_vec.z*inv_r0_vec.z); - x.real_value = - coordnum::switching_function(r0, inv_r0_vec, inv_r0sq_vec, en, ed, - atom_groups[0]->pos_x(0), - atom_groups[0]->pos_y(0), - atom_groups[0]->pos_z(0), - atom_groups[0]->pos_x(1), - atom_groups[0]->pos_y(1), - atom_groups[0]->pos_z(1), - atom_groups[0]->grad_x(0), - atom_groups[0]->grad_y(0), - atom_groups[0]->grad_z(0), - atom_groups[0]->grad_x(1), - atom_groups[0]->grad_y(1), - atom_groups[0]->grad_z(1), - NULL, 0.0); + + const cvm::rvector inv_r0_vec{ + 1.0 / r0_vec.x, + 1.0 / r0_vec.y, + 1.0 / r0_vec.z + }; + cvm::rvector const inv_r0sq_vec{ + inv_r0_vec.x * inv_r0_vec.x, + inv_r0_vec.y * inv_r0_vec.y, + inv_r0_vec.z * inv_r0_vec.z + }; + + x.real_value = coordnum::compute_pair_coordnum(inv_r0_vec, inv_r0sq_vec, en, ed, + atom_groups[0]->pos_x(0), + atom_groups[0]->pos_y(0), + atom_groups[0]->pos_z(0), + atom_groups[0]->pos_x(1), + atom_groups[0]->pos_y(1), + atom_groups[0]->pos_z(1), + atom_groups[0]->grad_x(0), + atom_groups[0]->grad_y(0), + atom_groups[0]->grad_z(0), + atom_groups[0]->grad_x(1), + atom_groups[0]->grad_y(1), + atom_groups[0]->grad_z(1), + 0.0, 1.0e20, + boundary_conditions); // Skip the gradient } void colvar::h_bond::calc_gradients() { - int const flags = coordnum::ef_gradients; - cvm::rvector const r0_vec(0.0); // TODO enable the flag? - const cvm::rvector inv_r0_vec( - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); - cvm::rvector const inv_r0sq_vec( + int constexpr flags = coordnum::ef_gradients; + const cvm::rvector inv_r0_vec{ + 1.0 / r0_vec.x, + 1.0 / r0_vec.y, + 1.0 / r0_vec.z + }; + cvm::rvector const inv_r0sq_vec{ inv_r0_vec.x*inv_r0_vec.x, inv_r0_vec.y*inv_r0_vec.y, - inv_r0_vec.z*inv_r0_vec.z); - coordnum::switching_function(r0, inv_r0_vec, inv_r0sq_vec, en, ed, - atom_groups[0]->pos_x(0), - atom_groups[0]->pos_y(0), - atom_groups[0]->pos_z(0), - atom_groups[0]->pos_x(1), - atom_groups[0]->pos_y(1), - atom_groups[0]->pos_z(1), - atom_groups[0]->grad_x(0), - atom_groups[0]->grad_y(0), - atom_groups[0]->grad_z(0), - atom_groups[0]->grad_x(1), - atom_groups[0]->grad_y(1), - atom_groups[0]->grad_z(1), - NULL, 0.0); + inv_r0_vec.z*inv_r0_vec.z + }; + coordnum::compute_pair_coordnum(inv_r0_vec, inv_r0sq_vec, en, ed, + atom_groups[0]->pos_x(0), + atom_groups[0]->pos_y(0), + atom_groups[0]->pos_z(0), + atom_groups[0]->pos_x(1), + atom_groups[0]->pos_y(1), + atom_groups[0]->pos_z(1), + atom_groups[0]->grad_x(0), + atom_groups[0]->grad_y(0), + atom_groups[0]->grad_z(0), + atom_groups[0]->grad_x(1), + atom_groups[0]->grad_y(1), + atom_groups[0]->grad_z(1), + 0.0, 1.0e20, + boundary_conditions); } - colvar::selfcoordnum::selfcoordnum() { set_function_type("selfCoordNum"); - x.type(colvarvalue::type_scalar); - r0 = cvm::main()->proxy->angstrom_to_internal(4.0); } -int colvar::selfcoordnum::init(std::string const &conf) +template inline void colvar::selfcoordnum::selfcoordnum_sequential_loop() { - int error_code = cvc::init(conf); - - group1 = parse_group(conf, "group1"); - - if (!group1 || group1->size() == 0) { - return error_code | COLVARS_INPUT_ERROR; - } - - get_keyval(conf, "cutoff", r0, r0); - get_keyval(conf, "expNumer", en, en); - get_keyval(conf, "expDenom", ed, ed); - - - if ((en % 2) || (ed % 2)) { - error_code |= cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", - COLVARS_INPUT_ERROR); - } - - if ((en <= 0) || (ed <= 0)) { - error_code |= cvm::error("Error: negative exponent(s) provided.\n", COLVARS_INPUT_ERROR); - } + size_t const n = group1->size(); + bool *pairlist_elem = pairlist.get(); + + for (size_t i = 0; i < n - 1; i++) { + + cvm::real const x1 = group1->pos_x(i); + cvm::real const y1 = group1->pos_y(i); + cvm::real const z1 = group1->pos_z(i); + cvm::real &gx1 = group1->grad_x(i); + cvm::real &gy1 = group1->grad_y(i); + cvm::real &gz1 = group1->grad_z(i); + + for (size_t j = i + 1; j < n; j++) { + + cvm::real const x2 = group1->pos_x(j); + cvm::real const y2 = group1->pos_y(j); + cvm::real const z2 = group1->pos_z(j); + cvm::real &gx2 = group1->grad_x(j); + cvm::real &gy2 = group1->grad_y(j); + cvm::real &gz2 = group1->grad_z(j); + + bool const within = + ((flags & ef_use_pairlist) && (*pairlist_elem || (flags & ef_rebuild_pairlist))) || + !(flags & ef_use_pairlist); + + cvm::real const partial = + within ? compute_pair_coordnum(inv_r0_vec, inv_r0sq_vec, en, ed, + x1, y1, z1, x2, y2, z2, + gx1, gy1, gz1, gx2, gy2, gz2, + tolerance, tolerance_l2_max, boundary_conditions) + : 0.0; + + if ((flags & ef_use_pairlist) && (flags & ef_rebuild_pairlist)) { + *pairlist_elem = partial > 0.0 ? true : false; + } - if (!is_enabled(f_cvc_pbc_minimum_image)) { - cvm::log("Warning: only minimum-image distances are used by this variable.\n"); - } + x.real_value += partial; - get_keyval(conf, "tolerance", tolerance, tolerance); - if (tolerance > 0) { - get_keyval(conf, "pairListFrequency", pairlist_freq, pairlist_freq); - if ( ! (pairlist_freq > 0) ) { - error_code |= cvm::error("Error: non-positive pairlistfrequency provided.\n", - COLVARS_INPUT_ERROR); + if (flags & ef_use_pairlist) { + pairlist_elem++; + } } - size_t n = (group1->size()-1) * (group1->size()-1); - pairlist = new bool[n]; - for (size_t i = 0; i < n; i++) pairlist[i] = true; - } - - init_scalar_boundaries(0.0, static_cast((group1->size()-1) * - (group1->size()-1))); - - return error_code; -} - - -colvar::selfcoordnum::~selfcoordnum() -{ - if (pairlist) { - delete [] pairlist; } } template int colvar::selfcoordnum::compute_selfcoordnum() { - cvm::rvector const r0_vec(0.0); // TODO enable the flag? - - bool const use_pairlist = (pairlist != NULL); - bool const rebuild_pairlist = (pairlist != NULL) && - (cvm::step_relative() % pairlist_freq == 0); - - bool *pairlist_elem = use_pairlist ? pairlist : NULL; - size_t i = 0, j = 0; - size_t const n = group1->size(); - - // Always isotropic (TODO: enable the ellipsoid?) -#define CALL_KERNEL(flags) do { \ - const cvm::rvector inv_r0_vec( \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); \ - cvm::rvector const inv_r0sq_vec( \ - inv_r0_vec.x*inv_r0_vec.x, \ - inv_r0_vec.y*inv_r0_vec.y, \ - inv_r0_vec.z*inv_r0_vec.z); \ - for (i = 0; i < n - 1; i++) { \ - for (j = i + 1; j < n; j++) { \ - x.real_value += \ - coordnum::switching_function( \ - r0, inv_r0_vec, inv_r0sq_vec, en, ed, \ - group1->pos_x(i), \ - group1->pos_y(i), \ - group1->pos_z(i), \ - group1->pos_x(j), \ - group1->pos_y(j), \ - group1->pos_z(j), \ - group1->grad_x(i), \ - group1->grad_y(i), \ - group1->grad_z(i), \ - group1->grad_x(j), \ - group1->grad_y(j), \ - group1->grad_z(j), \ - &pairlist_elem, tolerance); \ - } \ - } \ -} while (0); + bool const use_pairlist = pairlist.get(); + bool const rebuild_pairlist = use_pairlist && (cvm::step_relative() % pairlist_freq == 0); if (use_pairlist) { if (rebuild_pairlist) { - int const flags = compute_flags | coordnum::ef_use_pairlist | - coordnum::ef_rebuild_pairlist; - CALL_KERNEL(flags); + int constexpr flags = compute_flags | ef_use_pairlist | ef_rebuild_pairlist; + selfcoordnum_sequential_loop(); } else { - int const flags = compute_flags | coordnum::ef_use_pairlist; - CALL_KERNEL(flags); + int constexpr flags = compute_flags | ef_use_pairlist; + selfcoordnum_sequential_loop(); } - } else { // if (use_pairlist) { - int const flags = compute_flags | coordnum::ef_null; - CALL_KERNEL(flags); + } else { + int constexpr flags = compute_flags | ef_null; + selfcoordnum_sequential_loop(); } -#undef CALL_KERNEL return COLVARS_OK; } @@ -635,125 +559,23 @@ void colvar::selfcoordnum::calc_gradients() } - -colvar::groupcoordnum::groupcoordnum() -{ - set_function_type("groupCoord"); - x.type(colvarvalue::type_scalar); - init_scalar_boundaries(0.0, 1.0); - colvarproxy *proxy = cvm::main()->proxy; - r0 = proxy->angstrom_to_internal(4.0); - r0_vec = cvm::rvector(proxy->angstrom_to_internal(4.0), - proxy->angstrom_to_internal(4.0), - proxy->angstrom_to_internal(4.0)); -} - - -int colvar::groupcoordnum::init(std::string const &conf) -{ - int error_code = distance::init(conf); - - // group1 and group2 are already initialized by distance() - if (group1->b_dummy || group2->b_dummy) { - return cvm::error("Error: neither group can be a dummy atom\n", COLVARS_INPUT_ERROR); - } - - bool const b_scale = get_keyval(conf, "cutoff", r0, r0); - - if (get_keyval(conf, "cutoff3", r0_vec, r0_vec)) { - if (b_scale) { - error_code |= - cvm::error("Error: cannot specify \"cutoff\" and \"cutoff3\" at the same time.\n", - COLVARS_INPUT_ERROR); - } - b_anisotropic = true; - // remove meaningless negative signs - if (r0_vec.x < 0.0) r0_vec.x *= -1.0; - if (r0_vec.y < 0.0) r0_vec.y *= -1.0; - if (r0_vec.z < 0.0) r0_vec.z *= -1.0; - } - - get_keyval(conf, "expNumer", en, en); - get_keyval(conf, "expDenom", ed, ed); - - if ((en % 2) || (ed % 2)) { - error_code |= cvm::error("Error: odd exponent(s) provided, can only use even ones.\n", - COLVARS_INPUT_ERROR); - } - - if ((en <= 0) || (ed <= 0)) { - error_code |= cvm::error("Error: negative exponent(s) provided.\n", COLVARS_INPUT_ERROR); - } - - if (!is_enabled(f_cvc_pbc_minimum_image)) { - cvm::log("Warning: only minimum-image distances are used by this variable.\n"); - } - - return error_code; -} +colvar::groupcoordnum::groupcoordnum() { set_function_type("groupCoord"); } void colvar::groupcoordnum::calc_value() { - const cvm::atom_pos A1 = group1->center_of_mass(); - const cvm::atom_pos A2 = group2->center_of_mass(); -#define CALL_KERNEL(flags) do { \ - const cvm::rvector inv_r0_vec( \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); \ - cvm::rvector const inv_r0sq_vec( \ - inv_r0_vec.x*inv_r0_vec.x, \ - inv_r0_vec.y*inv_r0_vec.y, \ - inv_r0_vec.z*inv_r0_vec.z); \ - cvm::rvector G1, G2; \ - const cvm::rvector r0sq_vec(r0_vec.x*r0_vec.x, \ - r0_vec.y*r0_vec.y, \ - r0_vec.z*r0_vec.z); \ - x.real_value = coordnum::switching_function(r0, inv_r0_vec, inv_r0sq_vec, en, ed, \ - A1.x, A1.y, A1.z, \ - A2.x, A2.y, A2.z, \ - G1.x, G1.y, G1.z, \ - G2.x, G2.y, G2.z, NULL, 0.0); \ -} while (0); - if (b_anisotropic) { - int const flags = coordnum::ef_anisotropic; - CALL_KERNEL(flags); + x.real_value = 0.0; + if (is_enabled(f_cvc_gradient)) { + constexpr int flags = ef_gradients; + compute_coordnum(); } else { - int const flags = coordnum::ef_null; - CALL_KERNEL(flags); + constexpr int flags = ef_null; + compute_coordnum(); } -#undef CALL_KERNEL } void colvar::groupcoordnum::calc_gradients() { - const cvm::atom_pos A1 = group1->center_of_mass(); - const cvm::atom_pos A2 = group2->center_of_mass(); - cvm::rvector G1(0, 0, 0), G2(0, 0, 0); -#define CALL_KERNEL(flags) do { \ - const cvm::rvector inv_r0_vec( \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.x : r0), \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.y : r0), \ - 1.0 / ((flags & coordnum::ef_anisotropic) ? r0_vec.z : r0)); \ - cvm::rvector const inv_r0sq_vec( \ - inv_r0_vec.x*inv_r0_vec.x, \ - inv_r0_vec.y*inv_r0_vec.y, \ - inv_r0_vec.z*inv_r0_vec.z); \ - coordnum::switching_function(r0, inv_r0_vec, inv_r0sq_vec, en, ed, \ - A1.x, A1.y, A1.z, \ - A2.x, A2.y, A2.z, \ - G1.x, G1.y, G1.z, \ - G2.x, G2.y, G2.z, NULL, 0.0); \ -} while (0); - if (b_anisotropic) { - int const flags = coordnum::ef_gradients | coordnum::ef_anisotropic; - CALL_KERNEL(flags); - } else { - int const flags = coordnum::ef_gradients; - CALL_KERNEL(flags); - } - group1->set_weighted_gradient(G1); - group2->set_weighted_gradient(G2); + // Gradients are computed by calc_value() if f_cvc_gradients is enabled } diff --git a/src/colvarcomp_coordnums.h b/src/colvarcomp_coordnums.h new file mode 100644 index 000000000..60a3cac47 --- /dev/null +++ b/src/colvarcomp_coordnums.h @@ -0,0 +1,283 @@ +// -*- c++ -*- + +// This file is part of the Collective Variables module (Colvars). +// The original version of Colvars and its updates are located at: +// https://github.com/Colvars/colvars +// Please update all Colvars source files before making any changes. +// If you wish to distribute your changes, please submit them to the +// Colvars repository at GitHub. +// + +#ifndef COLVARCOMP_COORDNUM_H +#define COLVARCOMP_COORDNUM_H + +#include + +#include "colvar.h" +#include "colvarcomp.h" +#include "colvarmodule.h" + + +/// \brief Colvar component: coordination number between two groups +/// (colvarvalue::type_scalar type, range [0:N1*N2]) +class colvar::coordnum : public colvar::cvc { +public: + + coordnum(); + virtual ~coordnum(); + virtual int init(std::string const &conf); + virtual void calc_value(); + virtual void calc_gradients(); + + enum { + ef_null = 0, + ef_gradients = 1, + ef_use_pairlist = (1 << 9), + ef_rebuild_pairlist = (1 << 10) + }; + + /// Compute the switching function (1-l2**(en/2))/(1-l2**(ed/2)) for a given squared scaled distance + /// @param[in] l2 Square norm of (Dx/r0x, Dy/r0y, Dz/r0z) + /// @param[out] dFdl2 Derivative of the result with respect to l2 + /// @param en Numerator exponent + /// @param ed Denominator exponent + /// @param pairlist_tol Pairlist tolerance + template + static cvm::real switching_function(cvm::real const &l2, cvm::real &dFdl2, int en, int ed, + cvm::real pairlist_tol); + + /// Main kernel for the coordination number + template + static cvm::real compute_pair_coordnum(cvm::rvector const &inv_r0_vec, + cvm::rvector const &inv_r0sq_vec, int en, int ed, + const cvm::real a1x, const cvm::real a1y, const cvm::real a1z, + const cvm::real a2x, const cvm::real a2y, const cvm::real a2z, + cvm::real &g1x, cvm::real &g1y, cvm::real &g1z, + cvm::real &g2x, cvm::real &g2y, cvm::real &g2z, + cvm::real pairlist_tol, cvm::real pairlist_tol_l2_max, + cvm::system_boundary_conditions const &bc); + + /// Workhorse function + template int compute_coordnum(); + + /// Workhorse function + template void main_loop(); + +protected: + /// First atom group + cvm::atom_group *group1 = nullptr; + /// Second atom group + cvm::atom_group *group2 = nullptr; + + /// Cutoff distances along each dimension + cvm::rvector r0_vec; + + /// Inverse of r0_vec + cvm::rvector inv_r0_vec; + + /// Square of inv_r0_vec + cvm::rvector inv_r0sq_vec; + + /// Set r0_vec and related fields + void update_cutoffs(cvm::rvector const &r0_vec_i); + + /// Integer exponent of the function numerator + int en = 6; + /// Integer exponent of the function denominator + int ed = 12; + + /// The number of pairwise distances being calculated + size_t num_pairs = 0; + + /// If true, group1 will be treated as a single atom + bool b_group1_center_only = false; + + /// If true, group2 will be treated as a single atom + bool b_group2_center_only = false; + + /// Tolerance for the pair list + cvm::real tolerance = 0.0; + + /// Value of the squared scaled distance (l^2) that matches the given tolerance + cvm::real tolerance_l2_max = 1.0e20; + + /// Recompute the value of tolerance_l2_max + void compute_tolerance_l2_max(); + + /// Frequency of update of the pair list + int pairlist_freq = 100; + + /// Pair list + std::unique_ptr pairlist; + +}; + + +/// \brief Colvar component: self-coordination number within a group +/// (colvarvalue::type_scalar type, range [0:N*(N-1)/2]) +class colvar::selfcoordnum : public colvar::coordnum { +public: + + selfcoordnum(); + virtual void calc_value(); + virtual void calc_gradients(); + + /// Workhorse function + template void selfcoordnum_sequential_loop(); + + /// Main workhorse function + template int compute_selfcoordnum(); +}; + + +/// \brief Colvar component: coordination number between two groups +/// (colvarvalue::type_scalar type, range [0:N1*N2]) +class colvar::groupcoordnum : public colvar::coordnum { +public: + groupcoordnum(); + virtual ~groupcoordnum() {} + virtual void calc_value(); + virtual void calc_gradients(); +}; + + +/// \brief Colvar component: hydrogen bond, defined as the product of +/// a colvar::coordnum and 1/2*(1-cos((180-ang)/ang_tol)) +/// (colvarvalue::type_scalar type, range [0:1]) +class colvar::h_bond : public colvar::cvc { +public: + /// Constructor for atoms already allocated + h_bond(cvm::atom_group::simple_atom const &acceptor, cvm::atom_group::simple_atom const &donor, + cvm::real r0, int en, int ed); + h_bond(); + virtual ~h_bond() {} + virtual int init(std::string const &conf); + virtual void calc_value(); + virtual void calc_gradients(); + +protected: + /// Cutoff distances along each dimension + cvm::rvector r0_vec; + /// Integer exponent of the function numerator + int en = 6; + /// Integer exponent of the function denominator + int ed = 8; +}; + + +template +inline cvm::real colvar::coordnum::switching_function(cvm::real const &l2, cvm::real &dFdl2, + int en, int ed, + cvm::real pairlist_tol) +{ + // Assume en and ed are even integers, and avoid sqrt in the following + int const en2 = en/2; + int const ed2 = ed/2; + + cvm::real const xn = cvm::integer_power(l2, en2); + cvm::real const xd = cvm::integer_power(l2, ed2); + cvm::real const eps_l2 = 1.0e-7; + cvm::real const h = l2 - 1.0; + cvm::real const en2_r = (cvm::real) en2; + cvm::real const ed2_r = (cvm::real) ed2; + cvm::real func_no_pairlist; + + if (std::abs(h) < eps_l2) { + // Order-2 Taylor expansion: c0 + c1*h + c2*h^2 + cvm::real const c0 = en2_r / ed2_r; + cvm::real const c1 = (en2_r * (en2_r - ed2_r)) / (2.0 * ed2_r); + cvm::real const c2 = (en2_r * (en2_r - ed2_r) * (2.0 * en2_r - ed2_r - 3.0)) / (12.0 * ed2_r); + func_no_pairlist = c0 + h * (c1 + h * c2); + } else { + func_no_pairlist = (1.0 - xn) / (1.0 - xd); + } + + cvm::real func, inv_one_pairlist_tol; + if (flags & ef_use_pairlist) { + inv_one_pairlist_tol = 1 / (1.0-pairlist_tol); + func = (func_no_pairlist - pairlist_tol) * inv_one_pairlist_tol; + } else { + func = func_no_pairlist; + } + + // If the value is too small and we are correcting for the tolerance, the result is negative + // and we need to exclude it rather than let it contribute to the sum or the gradients. + if (func < 0) + return 0; + + if (flags & ef_gradients) { + // Logarithmic derivative: 1st-order Taylor expansion around l2 = 1 + cvm::real log_deriv; + if (std::abs(h) < eps_l2) { + cvm::real const g0 = 0.5 * (en2_r - ed2_r); + cvm::real const g1 = ((en2_r - ed2_r) * (en2_r + ed2_r - 6.0)) / 12.0; + log_deriv = g0 + h * g1; + } else { + log_deriv = (ed2_r * xd / ((1.0 - xd) * l2)) - (en2_r * xn / ((1.0 - xn) * l2)); + } + dFdl2 = (flags & ef_use_pairlist) ? + func_no_pairlist * inv_one_pairlist_tol * log_deriv : + func * log_deriv; + } + + return func; +} + + +template +inline cvm::real colvar::coordnum::compute_pair_coordnum(cvm::rvector const &inv_r0_vec, + cvm::rvector const &inv_r0sq_vec, + int en, + int ed, + const cvm::real a1x, + const cvm::real a1y, + const cvm::real a1z, + const cvm::real a2x, + const cvm::real a2y, + const cvm::real a2z, + cvm::real& g1x, + cvm::real& g1y, + cvm::real& g1z, + cvm::real& g2x, + cvm::real& g2y, + cvm::real& g2z, + cvm::real pairlist_tol, + cvm::real pairlist_tol_l2_max, + cvm::system_boundary_conditions const &bc) +{ + const cvm::atom_pos pos1{a1x, a1y, a1z}; + const cvm::atom_pos pos2{a2x, a2y, a2z}; + cvm::rvector const diff = bc.position_distance(pos1, pos2); + + cvm::rvector const scal_diff(diff.x * inv_r0_vec.x, + diff.y * inv_r0_vec.y, + diff.z * inv_r0_vec.z); + cvm::real const l2 = scal_diff.norm2(); + if (flags & ef_use_pairlist) { + if (l2 > pairlist_tol_l2_max) { + // Exit if the distance is such that F(l2) < pairlist_tol + return 0.0; + } + } + + cvm::real dFdl2 = 0.0; + cvm::real F = switching_function(l2, dFdl2, en, ed, pairlist_tol); + + if ((flags & ef_gradients) && (F > 0.0)) { + cvm::rvector const dl2dx((2.0 * inv_r0sq_vec.x) * diff.x, + (2.0 * inv_r0sq_vec.y) * diff.y, + (2.0 * inv_r0sq_vec.z) * diff.z); + + const cvm::rvector G = dFdl2*dl2dx; + g1x += -1.0*G.x; + g1y += -1.0*G.y; + g1z += -1.0*G.z; + g2x += G.x; + g2y += G.y; + g2z += G.z; + } + + return F; +} + +#endif // COLVARCOMP_COORDNUM_H diff --git a/src/colvarcomp_distances.cpp b/src/colvarcomp_distances.cpp index f0af812b1..88c8444a9 100644 --- a/src/colvarcomp_distances.cpp +++ b/src/colvarcomp_distances.cpp @@ -51,12 +51,7 @@ int colvar::distance::init(std::string const &conf) void colvar::distance::calc_value() { - if (!is_enabled(f_cvc_pbc_minimum_image)) { - dist_v = group2->center_of_mass() - group1->center_of_mass(); - } else { - dist_v = cvm::position_distance(group1->center_of_mass(), - group2->center_of_mass()); - } + dist_v = boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass()); x.real_value = dist_v.norm(); } @@ -98,12 +93,8 @@ colvar::distance_vec::distance_vec() void colvar::distance_vec::calc_value() { - if (!is_enabled(f_cvc_pbc_minimum_image)) { - x.rvector_value = group2->center_of_mass() - group1->center_of_mass(); - } else { - x.rvector_value = cvm::position_distance(group1->center_of_mass(), - group2->center_of_mass()); - } + x.rvector_value = + boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass()); } @@ -126,19 +117,13 @@ void colvar::distance_vec::apply_force(colvarvalue const &force) cvm::real colvar::distance_vec::dist2(colvarvalue const &x1, colvarvalue const &x2) const { - if (is_enabled(f_cvc_pbc_minimum_image)) { - return (cvm::position_distance(x1.rvector_value, x2.rvector_value)).norm2(); - } - return (x2.rvector_value - x1.rvector_value).norm2(); + return (boundary_conditions.position_distance(x1.rvector_value, x2.rvector_value)).norm2(); } colvarvalue colvar::distance_vec::dist2_lgrad(colvarvalue const &x1, colvarvalue const &x2) const { - if (is_enabled(f_cvc_pbc_minimum_image)) { - return 2.0 * cvm::position_distance(x2.rvector_value, x1.rvector_value); - } - return 2.0 * (x2.rvector_value - x1.rvector_value); + return 2.0 * boundary_conditions.position_distance(x2.rvector_value, x1.rvector_value); } @@ -200,21 +185,12 @@ void colvar::distance_z::calc_value() cvm::rvector const M = main->center_of_mass(); cvm::rvector const R1 = ref1->center_of_mass(); if (fixed_axis) { - if (!is_enabled(f_cvc_pbc_minimum_image)) { - dist_v = M - R1; - } else { - dist_v = cvm::position_distance(R1, M); - } + dist_v = boundary_conditions.position_distance(R1, M); } else { cvm::rvector const R2 = ref2->center_of_mass(); cvm::rvector const C = 0.5 * (R1 + R2); - if (!is_enabled(f_cvc_pbc_minimum_image)) { - dist_v = M - C; - axis = R2 - R1; - } else { - dist_v = cvm::position_distance(C, M); - axis = cvm::position_distance(R1, R2); - } + dist_v = boundary_conditions.position_distance(C, M); + axis = boundary_conditions.position_distance(R1, R2); axis_norm = axis.norm(); axis = axis.unit(); } @@ -272,19 +248,9 @@ colvar::distance_xy::distance_xy() void colvar::distance_xy::calc_value() { - if (!is_enabled(f_cvc_pbc_minimum_image)) { - dist_v = main->center_of_mass() - ref1->center_of_mass(); - } else { - dist_v = cvm::position_distance(ref1->center_of_mass(), - main->center_of_mass()); - } + dist_v = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass()); if (!fixed_axis) { - if (!is_enabled(f_cvc_pbc_minimum_image)) { - v12 = ref2->center_of_mass() - ref1->center_of_mass(); - } else { - v12 = cvm::position_distance(ref1->center_of_mass(), - ref2->center_of_mass()); - } + v12 = boundary_conditions.position_distance(ref1->center_of_mass(), ref2->center_of_mass()); axis_norm = v12.norm(); axis = v12.unit(); } @@ -308,12 +274,7 @@ void colvar::distance_xy::calc_gradients() ref1->set_weighted_gradient(-1.0 * x_inv * dist_v_ortho); main->set_weighted_gradient( x_inv * dist_v_ortho); } else { - if (!is_enabled(f_cvc_pbc_minimum_image)) { - v13 = main->center_of_mass() - ref1->center_of_mass(); - } else { - v13 = cvm::position_distance(ref1->center_of_mass(), - main->center_of_mass()); - } + v13 = boundary_conditions.position_distance(ref1->center_of_mass(), main->center_of_mass()); A = (dist_v * axis) / axis_norm; ref1->set_weighted_gradient( (A - 1.0) * x_inv * dist_v_ortho); @@ -354,12 +315,8 @@ colvar::distance_dir::distance_dir() void colvar::distance_dir::calc_value() { - if (!is_enabled(f_cvc_pbc_minimum_image)) { - dist_v = group2->center_of_mass() - group1->center_of_mass(); - } else { - dist_v = cvm::position_distance(group1->center_of_mass(), - group2->center_of_mass()); - } + dist_v = + boundary_conditions.position_distance(group1->center_of_mass(), group2->center_of_mass()); x.rvector_value = dist_v.unit(); } @@ -452,7 +409,7 @@ int colvar::distance_inv::init(std::string const &conf) void colvar::distance_inv::calc_value() { -#define CALL_KERNEL(USE_PBC_MINIMUM_IMAGE) do { \ +#define CALL_KERNEL() do { \ const int factor = -1*(exponent/2); \ for (size_t i = 0; i < group1->size(); ++i) { \ const cvm::atom_pos pos1(group1->pos_x(i), \ @@ -464,11 +421,7 @@ void colvar::distance_inv::calc_value() group2->pos_y(j), \ group2->pos_z(j)); \ cvm::rvector dv; \ - if (USE_PBC_MINIMUM_IMAGE) { \ - dv = cvm::position_distance(pos1, pos2); \ - } else { \ - dv = pos2 - pos1; \ - } \ + dv = boundary_conditions.position_distance(pos1, pos2); \ cvm::real const d2 = dv.norm2(); \ cvm::real const dinv = cvm::integer_power(d2, factor); \ x.real_value += dinv; \ @@ -484,11 +437,7 @@ void colvar::distance_inv::calc_value() } \ } while (0); x.real_value = 0.0; - if (!is_enabled(f_cvc_pbc_minimum_image)) { - CALL_KERNEL(false); - } else { - CALL_KERNEL(true); - } + CALL_KERNEL(); x.real_value *= 1.0 / cvm::real(group1->size() * group2->size()); x.real_value = cvm::pow(x.real_value, -1.0/cvm::real(exponent)); @@ -537,7 +486,7 @@ int colvar::distance_pairs::init(std::string const &conf) void colvar::distance_pairs::calc_value() { x.vector1d_value.resize(group1->size() * group2->size()); -#define CALL_KERNEL(USE_PBC_MINIMUM_IMAGE) do { \ +#define CALL_KERNEL() do { \ for (size_t i1 = 0; i1 < group1->size(); ++i1) { \ const cvm::atom_pos pos1(group1->pos_x(i1), \ group1->pos_y(i1), \ @@ -547,9 +496,7 @@ void colvar::distance_pairs::calc_value() const cvm::atom_pos pos2(group2->pos_x(i2), \ group2->pos_y(i2), \ group2->pos_z(i2)); \ - const cvm::rvector dv = USE_PBC_MINIMUM_IMAGE ? \ - cvm::position_distance(pos1, pos2) : \ - pos2 - pos1; \ + const cvm::rvector dv = boundary_conditions.position_distance(pos1, pos2); \ cvm::real const d = dv.norm(); \ x.vector1d_value[i1*group2->size() + i2] = d; \ const cvm::rvector g2 = dv.unit(); \ @@ -563,11 +510,7 @@ void colvar::distance_pairs::calc_value() group1->grad_z(i1) += g1.z;*/ \ } \ } while (0); - if (!is_enabled(f_cvc_pbc_minimum_image)) { - CALL_KERNEL(false); - } else { - CALL_KERNEL(true); - } + CALL_KERNEL(); #undef CALL_KERNEL } @@ -580,7 +523,7 @@ void colvar::distance_pairs::calc_gradients() void colvar::distance_pairs::apply_force(colvarvalue const &force) { -#define CALL_KERNEL(USE_PBC_MINIMUM_IMAGE) do { \ +#define CALL_KERNEL() do { \ auto group1_force_obj = group1->get_group_force_object(); \ auto group2_force_obj = group2->get_group_force_object(); \ for (size_t i1 = 0; i1 < group1->size(); i1++) { \ @@ -592,9 +535,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force) const cvm::atom_pos pos2(group2->pos_x(i2), \ group2->pos_y(i2), \ group2->pos_z(i2)); \ - const cvm::rvector dv = USE_PBC_MINIMUM_IMAGE ? \ - cvm::position_distance(pos1, pos2) : \ - pos2 - pos1; \ + const cvm::rvector dv = boundary_conditions.position_distance(pos1, pos2); \ cvm::real const d = dv.norm(); \ x.vector1d_value[i1*group2->size() + i2] = d; \ const cvm::rvector f2 = force[i1*group2->size() + i2] * dv.unit();\ @@ -604,11 +545,7 @@ void colvar::distance_pairs::apply_force(colvarvalue const &force) group1_force_obj.add_atom_force(i1, f1); \ } \ } while (0); - if (!is_enabled(f_cvc_pbc_minimum_image)) { - CALL_KERNEL(false); - } else { - CALL_KERNEL(true); - } + CALL_KERNEL(); #undef CALL_KERNEL } diff --git a/src/colvarcomp_protein.cpp b/src/colvarcomp_protein.cpp index d5b37ff9d..d58edf2d7 100644 --- a/src/colvarcomp_protein.cpp +++ b/src/colvarcomp_protein.cpp @@ -13,6 +13,7 @@ #include "colvarvalue.h" #include "colvar.h" #include "colvarcomp.h" +#include "colvarcomp_coordnums.h" colvar::alpha_angles::alpha_angles() diff --git a/src/colvarmodule.cpp b/src/colvarmodule.cpp index a7e3111eb..fc50b53ab 100644 --- a/src/colvarmodule.cpp +++ b/src/colvarmodule.cpp @@ -2484,13 +2484,6 @@ void cvm::request_total_force() } -cvm::rvector cvm::position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) -{ - return proxy->position_distance(pos1, pos2); -} - - cvm::real cvm::rand_gaussian(void) { return proxy->rand_gaussian(); diff --git a/src/colvarmodule.h b/src/colvarmodule.h index 0470dd942..8fb7c4b92 100644 --- a/src/colvarmodule.h +++ b/src/colvarmodule.h @@ -227,6 +227,7 @@ static inline real acos(real const &x) template class matrix2d; class quaternion; class rotation; + class system_boundary_conditions; class usage; class memory_stream; @@ -373,9 +374,13 @@ static inline real acos(real const &x) std::vector *biases_active(); /// \brief Whether debug output should be enabled (compile-time option) - static inline bool debug() + static inline bool constexpr debug() { +#if (defined(__HIP_DEVICE_COMPILE__)) || (defined(__CUDA_ARCH__)) + return false; +#else return COLVARS_DEBUG; +#endif } /// How many objects (variables and biases) are configured yet? @@ -793,11 +798,6 @@ static inline real acos(real const &x) return 5; } - /// \brief Get the distance between two atomic positions with pbcs handled - /// correctly - static rvector position_distance(atom_pos const &pos1, - atom_pos const &pos2); - /// \brief Names of .ndx files that have been loaded std::vector index_file_names; diff --git a/src/colvarproxy_system.cpp b/src/colvarproxy_system.cpp index 6b082b8ed..2ffef4149 100644 --- a/src/colvarproxy_system.cpp +++ b/src/colvarproxy_system.cpp @@ -22,12 +22,10 @@ colvarproxy_system::colvarproxy_system() timestep_ = 1.0; target_temperature_ = 0.0; boltzmann_ = 0.001987191; // Default: kcal/mol/K - boundaries_type = boundaries_unsupported; total_force_requested = false; indirect_lambda_biasing_force = 0.0; cached_alch_lambda_changed = false; cached_alch_lambda = -1.0; - reset_pbc_lattice(); } @@ -90,73 +88,10 @@ bool colvarproxy_system::total_forces_same_step() const } -inline int round_to_integer(cvm::real x) +cvm::rvector colvarproxy_system::position_distance_engine(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const { - return int(cvm::floor(x+0.5)); -} - - -void colvarproxy_system::update_pbc_lattice() -{ - // Periodicity is assumed in all directions - - if (boundaries_type == boundaries_unsupported || - boundaries_type == boundaries_non_periodic) { - cvm::error("Error: setting PBC lattice with unsupported boundaries.\n", - COLVARS_BUG_ERROR); - return; - } - - { - cvm::rvector const v = cvm::rvector::outer(unit_cell_y, unit_cell_z); - reciprocal_cell_x = v/(v*unit_cell_x); - } - { - cvm::rvector const v = cvm::rvector::outer(unit_cell_z, unit_cell_x); - reciprocal_cell_y = v/(v*unit_cell_y); - } - { - cvm::rvector const v = cvm::rvector::outer(unit_cell_x, unit_cell_y); - reciprocal_cell_z = v/(v*unit_cell_z); - } -} - - -void colvarproxy_system::reset_pbc_lattice() -{ - unit_cell_x.reset(); - unit_cell_y.reset(); - unit_cell_z.reset(); - reciprocal_cell_x.reset(); - reciprocal_cell_y.reset(); - reciprocal_cell_z.reset(); -} - - -cvm::rvector colvarproxy_system::position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) - const -{ - if (boundaries_type == boundaries_unsupported) { - cvm::error("Error: unsupported boundary conditions.\n", COLVARS_INPUT_ERROR); - } - - cvm::rvector diff = (pos2 - pos1); - - if (boundaries_type == boundaries_non_periodic) return diff; - - cvm::real const x_shift = round_to_integer(reciprocal_cell_x*diff); - cvm::real const y_shift = round_to_integer(reciprocal_cell_y*diff); - cvm::real const z_shift = round_to_integer(reciprocal_cell_z*diff); - - diff.x -= x_shift*unit_cell_x.x + y_shift*unit_cell_y.x + - z_shift*unit_cell_z.x; - diff.y -= x_shift*unit_cell_x.y + y_shift*unit_cell_y.y + - z_shift*unit_cell_z.y; - diff.z -= x_shift*unit_cell_x.z + y_shift*unit_cell_y.z + - z_shift*unit_cell_z.z; - - return diff; + return boundaries_.position_distance(pos1, pos2); } diff --git a/src/colvarproxy_system.h b/src/colvarproxy_system.h index 38f3b5a64..fe93fa3c8 100644 --- a/src/colvarproxy_system.h +++ b/src/colvarproxy_system.h @@ -10,6 +10,7 @@ #ifndef COLVARPROXY_SYSTEM_H #define COLVARPROXY_SYSTEM_H +#include "colvars_system.h" /// Methods for accessing the simulation system (PBCs, integrator, etc) class colvarproxy_system { @@ -85,15 +86,17 @@ class colvarproxy_system { /// Pass restraint energy value for current timestep to MD engine virtual void add_energy(cvm::real energy); - /// \brief Get the PBC-aware distance vector between two positions - virtual cvm::rvector position_distance(cvm::atom_pos const &pos1, - cvm::atom_pos const &pos2) const; + /// Account for system boundaries within the Colvars library (as opposed to using the MD engine) + inline bool & use_internal_pbc() { return use_internal_pbc_; } - /// Recompute PBC reciprocal lattice (assumes XYZ periodicity) - void update_pbc_lattice(); + /// Get the PBC-aware distance vector between two positions using the MD engine (for VMD, this is internal) + virtual cvm::rvector position_distance_engine(cvm::atom_pos const &pos1, + cvm::atom_pos const &pos2) const; - /// Set the lattice vectors to zero - void reset_pbc_lattice(); + /// Get the current system boundary conditions + inline cvm::system_boundary_conditions const &get_system_boundaries() const { + return boundaries_; + } /// \brief Tell the proxy whether total forces are needed (they may not /// always be available) @@ -179,26 +182,12 @@ class colvarproxy_system { /// Whether the total forces have been requested bool total_force_requested; - /// \brief Type of boundary conditions - /// - /// Orthogonal and triclinic cells are made available to objects. - /// For any other conditions (mixed periodicity, triclinic cells in LAMMPS) - /// minimum-image distances are computed by the host engine regardless. - enum Boundaries_type { - boundaries_non_periodic, - boundaries_pbc_ortho, - boundaries_pbc_triclinic, - boundaries_unsupported - }; - - /// Type of boundary conditions - Boundaries_type boundaries_type; - - /// Bravais lattice vectors - cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z; - - /// Reciprocal lattice vectors - cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z; + /// Use the PBC functions from the Colvars library (as opposed to MD engine) + bool use_internal_pbc_ = false; + + /// Current system boundary conditions + cvm::system_boundary_conditions boundaries_; }; + #endif diff --git a/src/colvars_system.h b/src/colvars_system.h new file mode 100644 index 000000000..986a7cb55 --- /dev/null +++ b/src/colvars_system.h @@ -0,0 +1,230 @@ +// -*- c++ -*- + +// This file is part of the Collective Variables module (Colvars). +// The original version of Colvars and its updates are located at: +// https://github.com/Colvars/colvars +// Please update all Colvars source files before making any changes. +// If you wish to distribute your changes, please submit them to the +// Colvars repository at GitHub. + +#ifndef COLVARS_SYSTEM_H +#define COLVARS_SYSTEM_H + +#include "colvartypes.h" + +/// Class to store the system's boundary conditions +class colvarmodule::system_boundary_conditions { +public: + + /// Type of boundary conditions defined for the current computation + enum class types { + non_periodic, /// All three dimensions are non-periodic + mixed, /// Some dimensions are periodic, but others are not + pbc_orthogonal, /// All three dimensions are periodic, lattice vectors are orthogonal + pbc_triclinic, /// All three dimensions are periodic, lattice vectors are *not* orthogonal + unsupported /// Unsupported boundary conditions + }; + + /// Default constructor + inline COLVARS_HOST_DEVICE system_boundary_conditions() { reset(); } + + /// Copy constructor + system_boundary_conditions(system_boundary_conditions const &) = default; + + /// Type of boundary conditions in the current computation + inline COLVARS_HOST_DEVICE types type() const { return type_; } + + /// Set the type of boundary explicitly + inline COLVARS_HOST_DEVICE void set_type(types t) { type_ = t; } + + /// Compute the distance between two positions + cvm::rvector position_distance(cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) const; + + /// Compute a shift vector that accounts for tilt factors up to 0.5 + cvm::rvector get_triclinic_shift(cvm::rvector const &diff) const; + + /// Reset to defaults (non-periodic) + inline COLVARS_HOST_DEVICE void reset() { + periodic_x = periodic_y = periodic_z = false; + type_ = types::non_periodic; + unit_cell_x.reset(); + unit_cell_y.reset(); + unit_cell_z.reset(); + reciprocal_cell_x.reset(); + reciprocal_cell_y.reset(); + reciprocal_cell_z.reset(); + } + + /// Set from explicit boundary configuration + void set_boundaries(bool periodic_x_in, + bool periodic_y_in, + bool periodic_z_in, + cvm::rvector const &A, + cvm::rvector const &B, + cvm::rvector const &C); + +protected: + + /// Type of boundary conditions in the current computation + types type_ = types::non_periodic; + + /// Bravais lattice vectors + cvm::rvector unit_cell_x, unit_cell_y, unit_cell_z; + + /// Reciprocal lattice vectors + cvm::rvector reciprocal_cell_x, reciprocal_cell_y, reciprocal_cell_z; + + /// Periodic flags in each dimension + bool periodic_x = false, periodic_y = false, periodic_z = false; +}; + + +/// Set from explicit boundary configuration +inline COLVARS_HOST_DEVICE void cvm::system_boundary_conditions::set_boundaries( + bool periodic_x_in, + bool periodic_y_in, + bool periodic_z_in, + cvm::rvector const &A, + cvm::rvector const &B, + cvm::rvector const &C) +{ + constexpr double diagonal_tol2 = 1.0e-10; + + periodic_x = periodic_x_in; + periodic_y = periodic_y_in; + periodic_z = periodic_z_in; + + if ((periodic_x == periodic_y) && (periodic_x == periodic_z)) { + if (periodic_x) { + // Temporarily set as fully-periodic & orthogonal; will check below for triclinic + set_type(types::pbc_orthogonal); + } else { + set_type(types::non_periodic); + return; + } + } else { + set_type(types::mixed); + } + + bool off_diagonal = false; + + if (periodic_x) { + unit_cell_x = A; + if ((A.y * A.y) > diagonal_tol2 || (A.z * A.z) > diagonal_tol2) { + off_diagonal = true; + } else { + reciprocal_cell_x = unit_cell_x/unit_cell_x.norm2(); + } + } else { + unit_cell_x.reset(); + reciprocal_cell_x.reset(); + } + + if (periodic_y) { + unit_cell_y = B; + if ((B.x * B.x) > diagonal_tol2 || (B.z * B.z) > diagonal_tol2) { + off_diagonal = true; + } else { + reciprocal_cell_y = unit_cell_y/unit_cell_y.norm2(); + } + } else { + unit_cell_y.reset(); + reciprocal_cell_y.reset(); + } + + if (periodic_z) { + unit_cell_z = C; + if ((C.x * C.x) > diagonal_tol2 || (C.y * C.y) > diagonal_tol2) { + off_diagonal = true; + } else { + reciprocal_cell_z = unit_cell_z/unit_cell_z.norm2(); + } + } else { + unit_cell_z.reset(); + reciprocal_cell_z.reset(); + } + + if (type() == types::pbc_orthogonal && off_diagonal) { + set_type(types::pbc_triclinic); + } + + if (type() == types::pbc_triclinic) { + cvm::rvector const v_yz = cvm::rvector::outer(unit_cell_y, unit_cell_z); + reciprocal_cell_x = v_yz / (v_yz * unit_cell_x); + cvm::rvector const v_zx = cvm::rvector::outer(unit_cell_z, unit_cell_x); + reciprocal_cell_y = v_zx / (v_zx * unit_cell_y); + cvm::rvector const v_xy = cvm::rvector::outer(unit_cell_x, unit_cell_y); + reciprocal_cell_z = v_xy / (v_xy * unit_cell_z); + } + + if (cvm::debug()) { + cvm::log("periodic_x = " + cvm::to_str(periodic_x)); + cvm::log("periodic_y = " + cvm::to_str(periodic_y)); + cvm::log("periodic_z = " + cvm::to_str(periodic_z)); + + cvm::log("unit_cell_x = " + cvm::to_str(unit_cell_x)); + cvm::log("unit_cell_y = " + cvm::to_str(unit_cell_y)); + cvm::log("unit_cell_z = " + cvm::to_str(unit_cell_z)); + + cvm::log("reciprocal_cell_x = " + cvm::to_str(reciprocal_cell_x)); + cvm::log("reciprocal_cell_y = " + cvm::to_str(reciprocal_cell_y)); + cvm::log("reciprocal_cell_z = " + cvm::to_str(reciprocal_cell_z)); + } +} + + +inline COLVARS_HOST_DEVICE cvm::rvector cvm::system_boundary_conditions::position_distance( + cvm::atom_pos const &pos1, cvm::atom_pos const &pos2) const +{ + cvm::rvector diff = (pos2 - pos1); + + if (type() == types::non_periodic || type() == types::unsupported) { + return diff; + } + + cvm::real const x_shift = ::floor(reciprocal_cell_x * diff + 0.5); + cvm::real const y_shift = ::floor(reciprocal_cell_y * diff + 0.5); + cvm::real const z_shift = ::floor(reciprocal_cell_z * diff + 0.5); + + diff.x -= x_shift * unit_cell_x.x + y_shift * unit_cell_y.x + z_shift * unit_cell_z.x; + diff.y -= x_shift * unit_cell_x.y + y_shift * unit_cell_y.y + z_shift * unit_cell_z.y; + diff.z -= x_shift * unit_cell_x.z + y_shift * unit_cell_y.z + z_shift * unit_cell_z.z; + + if (type() != types::pbc_orthogonal) { + // Matches both "mixed" and "pbc_triclinic", because reciprocal cell vectors are not used + diff += get_triclinic_shift(diff); + } + + return diff; +} + + +inline COLVARS_HOST_DEVICE cvm::rvector cvm::system_boundary_conditions::get_triclinic_shift( + cvm::rvector const &diff) const +{ + cvm::real min_dist2 = diff.norm2(); + cvm::rvector result{0.0, 0.0, 0.0}; + + int const nx = periodic_x ? 1 : 0; + int const ny = periodic_y ? 1 : 0; + int const nz = periodic_z ? 1 : 0; + + // Loop over neighboring cells to find a shorter distance + for (int ix = -nx; ix <= nx; ix++) { + for (int iy = -ny; iy <= ny; iy++) { + for (int iz = -nz; iz <= nz; iz++) { + cvm::rvector const shift = ix * unit_cell_x + iy * unit_cell_y + iz * unit_cell_z; + cvm::real const this_dist2 = (diff + shift).norm2(); + if (this_dist2 < min_dist2) { + result = shift; + min_dist2 = this_dist2; + } + } + } + } + + return result; +} + + +#endif diff --git a/tests/input_files/compare_test.sh b/tests/input_files/compare_test.sh index 6ee2073fe..93a45a07d 100755 --- a/tests/input_files/compare_test.sh +++ b/tests/input_files/compare_test.sh @@ -22,13 +22,17 @@ fi echo "Top level directory is ${TOPDIR}" -SPIFF=$(${TOPDIR}/devel-tools/get_spiff) -if [ $? != 0 ] ; then +if hash spiff >& /dev/null ; then + SPIFF=spiff +else + SPIFF=$(${TOPDIR}/devel-tools/get_spiff) + if [ $? != 0 ] ; then echo "Error: spiff is not available and could not be downloaded/built." >& 2 exit 1 -else + else echo "Using spiff executable from $SPIFF" hash -p ${SPIFF} spiff + fi fi dir=$1 diff --git a/update-colvars-code.sh b/update-colvars-code.sh index 990cf8532..b544ad241 100755 --- a/update-colvars-code.sh +++ b/update-colvars-code.sh @@ -494,7 +494,12 @@ then condcopy "${src}" "${target}/src/${tgt}" done - condcopy "${source}/vmd/src/colvars_files.pl" "${target}/src/colvars_files.pl" + bash "${source}/devel-tools/generate-vmd-makefile.sh" \ + "${source}/src/"*.cpp \ + "${source}/src/"*.h \ + "${source}/vmd/src/"*.C \ + "${source}/vmd/src/"*.h \ + > "${target}/src/colvars_files.pl" echo ' done.' diff --git a/vmd/src/colvarproxy_vmd.C b/vmd/src/colvarproxy_vmd.C index f5f310807..78347943e 100644 --- a/vmd/src/colvarproxy_vmd.C +++ b/vmd/src/colvarproxy_vmd.C @@ -80,6 +80,9 @@ colvarproxy_vmd::colvarproxy_vmd(Tcl_Interp *interp, VMDApp *v, int molid) { engine_name_ = "VMD"; version_int = get_version_from_string(COLVARPROXY_VERSION); + + use_internal_pbc_ = true; + b_simulation_running = false; // both fields are taken from data structures already available @@ -221,36 +224,18 @@ int colvarproxy_vmd::update_input() angstrom_to_internal(vmdpos[atoms_ids[i]*3+2])); } - Timestep const *ts = vmdmol->get_frame(vmdmol_frame); - { - // Get lattice vectors - float A[3]; - float B[3]; - float C[3]; - ts->get_transform_vectors(A, B, C); - unit_cell_x.set(angstrom_to_internal(A[0]), angstrom_to_internal(A[1]), angstrom_to_internal(A[2])); - unit_cell_y.set(angstrom_to_internal(B[0]), angstrom_to_internal(B[1]), angstrom_to_internal(B[2])); - unit_cell_z.set(angstrom_to_internal(C[0]), angstrom_to_internal(C[1]), angstrom_to_internal(C[2])); - } - - if (ts->a_length + ts->b_length + ts->c_length < 1.0e-6) { - boundaries_type = boundaries_non_periodic; - reset_pbc_lattice(); - } else if ((ts->a_length > 1.0e-6) && - (ts->b_length > 1.0e-6) && - (ts->c_length > 1.0e-6)) { - if (((ts->alpha-90.0)*(ts->alpha-90.0)) + - ((ts->beta-90.0)*(ts->beta-90.0)) + - ((ts->gamma-90.0)*(ts->gamma-90.0)) < 1.0e-6) { - boundaries_type = boundaries_pbc_ortho; - } else { - boundaries_type = boundaries_pbc_triclinic; - } - colvarproxy_system::update_pbc_lattice(); - } else { - boundaries_type = boundaries_unsupported; - } + + float A[3]; + float B[3]; + float C[3]; + ts->get_transform_vectors(A, B, C); + boundaries_.set_boundaries((ts->a_length > 0.0), + (ts->b_length > 0.0), + (ts->c_length > 0.0), + cvm::rvector{A[0], A[1], A[2]}, + cvm::rvector{B[0], B[1], B[2]}, + cvm::rvector{C[0], C[1], C[2]}); return error_code; } @@ -809,7 +794,7 @@ void colvarproxy_vmd::compute_voldata(VolumetricData const *voldata, cvm::atom_pos const origin(0.0, 0.0, 0.0); for (; i < atoms->size(); ++i) { // Wrap around the origin - cvm::rvector const wrapped_pos = position_distance( + cvm::rvector const wrapped_pos = position_distance_engine( origin, cvm::atom_pos(atoms->pos_x(i), atoms->pos_y(i), atoms->pos_z(i))); coord[0] = internal_to_angstrom(wrapped_pos.x); coord[1] = internal_to_angstrom(wrapped_pos.y); diff --git a/vmd/src/colvars_files.pl b/vmd/src/colvars_files.pl deleted file mode 100644 index 819cd0f27..000000000 --- a/vmd/src/colvars_files.pl +++ /dev/null @@ -1,111 +0,0 @@ -# List of files for the Colvars module - -our ( $colvars_defines ); - -our ( @colvars_cc ); -our ( @colvars_cu ); -our ( @colvars_ccpp ); -our ( @colvars_h ); - -$colvars_defines = " -DVMDCOLVARS"; - -@colvars_cc = (); -@colvars_cu = (); -@colvars_ccpp = ('colvaratoms.C', - 'colvaratoms_gpu.C', - 'colvarbias.C', - 'colvarbias_abf.C', - 'colvarbias_abmd.C', - 'colvarbias_alb.C', - 'colvarbias_histogram.C', - 'colvarbias_meta.C', - 'colvarbias_restraint.C', - 'colvarbias_histogram_reweight_amd.C', - 'colvarbias_opes.C', - 'colvar.C', - 'colvarcomp_neuralnetwork.C', - 'colvar_neuralnetworkcompute.C', - 'colvarcomp_torchann.C', - 'colvarcomp.C', - 'colvarcomp_alchlambda.C', - 'colvarcomp_angles.C', - 'colvarcomp_apath.C', - 'colvarcomp_combination.C', - 'colvarcomp_coordnums.C', - 'colvarcomp_distances.C', - 'colvarcomp_gpath.C', - 'colvarcomp_protein.C', - 'colvarcomp_rotations.C', - 'colvarcomp_volmaps.C', - 'colvardeps.C', - 'colvargrid.C', - 'colvargrid_integrate.C', - 'colvarmodule.C', - 'colvarparams.C', - 'colvarparse.C', - 'colvarproxy.C', - 'colvarproxy_io.C', - 'colvarproxy_replicas.C', - 'colvarproxy_system.C', - 'colvarproxy_tcl.C', - 'colvarproxy_vmd.C', - 'colvarproxy_gpu.C', - 'colvarproxy_volmaps.C', - 'colvarscript.C', - 'colvarscript_commands.C', - 'colvarscript_commands_bias.C', - 'colvarscript_commands_colvar.C', - 'colvars_memstream.C', - 'colvartypes.C', - 'colvarvalue.C', - 'colvar_gpu_support.C', - 'colvar_gpu_calc.C', - 'colvar_rotation_derivative.C', - 'nr_jacobi.C'); -@colvars_h = ('colvar_UIestimator.h', - 'colvar_arithmeticpath.h', - 'colvar_geometricpath.h', - 'colvar_neuralnetworkcompute.h', - 'colvar_rotation_derivative.h', - 'colvaratoms.h', - 'colvaratoms_gpu.h', - 'colvarbias.h', - 'colvarbias_abf.h', - 'colvarbias_abmd.h', - 'colvarbias_alb.h', - 'colvarbias_histogram.h', - 'colvarbias_histogram_reweight_amd.h', - 'colvarbias_meta.h', - 'colvarbias_restraint.h', - 'colvarbias_opes.h', - 'colvarcomp.h', - 'colvarcomp_torchann.h', - 'colvardeps.h', - 'colvargrid.h', - 'colvargrid_def.h', - 'colvargrid_integrate.h', - 'colvar.h', - 'colvarmodule.h', - 'colvarmodule_refs.h', - 'colvarmodule_utils.h', - 'colvarparams.h', - 'colvarparse.h', - 'colvarproxy.h', - 'colvarproxy_io.h', - 'colvarproxy_replicas.h', - 'colvarproxy_system.h', - 'colvarproxy_tcl.h', - 'colvarproxy_vmd.h', - 'colvarproxy_gpu.h', - 'colvarproxy_volmaps.h', - 'colvarscript.h', - 'colvarscript_commands.h', - 'colvarscript_commands_bias.h', - 'colvarscript_commands_colvar.h', - 'colvars_memstream.h', - 'colvars_version.h', - 'colvartypes.h', - 'colvarvalue.h', - 'colvar_gpu_support.h', - 'colvar_gpu_calc.h', - 'nr_jacobi.h');