Added QMMMForceProvider class for QMMM MdModule
authorDmitry Morozov <aracsmd@gmail.com>
Thu, 23 Sep 2021 08:03:08 +0000 (08:03 +0000)
committerJoe Jordan <ejjordan12@gmail.com>
Thu, 23 Sep 2021 08:03:08 +0000 (08:03 +0000)
cmake/gmxManageCP2K.cmake
src/gromacs/applied_forces/qmmm/CMakeLists.txt
src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp [new file with mode: 0644]
src/gromacs/applied_forces/qmmm/qmmmforceprovider.h [new file with mode: 0644]
src/gromacs/applied_forces/qmmm/qmmmforceprovider_stub.cpp [new file with mode: 0644]
src/gromacs/applied_forces/qmmm/qmmmtypes.h
src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt
src/gromacs/applied_forces/qmmm/tests/qmmmforceprovider.cpp [new file with mode: 0644]

index 99f5b70c040d571f75972a96cef642f33f178e18..d23ca3392e1c8a3a1f53c9287d15644bbdfe23a4 100644 (file)
@@ -41,6 +41,9 @@ if(GMX_CP2K)
         message(FATAL_ERROR "To build GROMACS with CP2K Interface both CP2K_DIR and CP2K_LINKER_FLAGS should be defined")
     endif()
 
+    # Add directory with libcp2k.h into system include directories
+    include_directories(SYSTEM "${CP2K_DIR}/../../../src/start")
+
     # Add libcp2k and DBCSR for linking 
     set(CMAKE_CXX_STANDARD_LIBRARIES "${CMAKE_CXX_STANDARD_LIBRARIES} -Wl,--allow-multiple-definition -L${CP2K_DIR} -lcp2k -L${CP2K_DIR}/exts/dbcsr -ldbcsr")
 
index ba20be38243be8a7c794c28f613adeda94e27c12..3707aa1188f206a9b67e9af5e538d663338875bc 100644 (file)
@@ -37,6 +37,18 @@ gmx_add_libgromacs_sources(
     qmmmtopologypreprocessor.cpp
     )
 
+# If we have libcp2k linked then compile qmmmforceprovider.cpp
+# In other case compile stub implementation qmmmforceprovider_stub.cpp
+if (GMX_CP2K)
+gmx_add_libgromacs_sources(
+    qmmmforceprovider.cpp
+    )
+else()
+gmx_add_libgromacs_sources(
+    qmmmforceprovider_stub.cpp
+    )
+endif()
+
 if (BUILD_TESTING)
     add_subdirectory(tests)
 endif()
diff --git a/src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp b/src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp
new file mode 100644 (file)
index 0000000..f685efd
--- /dev/null
@@ -0,0 +1,302 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Implements force provider for QMMM
+ *
+ * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
+ * \author Christian Blau <blau@kth.se>
+ * \ingroup module_applied_forces
+ */
+
+#include "gmxpre.h"
+
+#include "qmmmforceprovider.h"
+
+#include "gromacs/domdec/domdec_struct.h"
+#include "gromacs/gmxlib/network.h"
+#include "gromacs/math/units.h"
+#include "gromacs/mdtypes/commrec.h"
+#include "gromacs/mdtypes/enerdata.h"
+#include "gromacs/utility/exceptions.h"
+#include "gromacs/utility/stringutil.h"
+#include "gromacs/utility/textwriter.h"
+
+#include <libcp2k.h>
+
+namespace gmx
+{
+
+QMMMForceProvider::QMMMForceProvider(const QMMMParameters& parameters,
+                                     const LocalAtomSet&   localQMAtomSet,
+                                     const LocalAtomSet&   localMMAtomSet,
+                                     PbcType               pbcType,
+                                     const MDLogger&       logger) :
+    parameters_(parameters),
+    qmAtoms_(localQMAtomSet),
+    mmAtoms_(localMMAtomSet),
+    pbcType_(pbcType),
+    logger_(logger),
+    box_{ { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } }
+{
+}
+
+QMMMForceProvider::~QMMMForceProvider()
+{
+    if (force_env_ != -1)
+    {
+        cp2k_destroy_force_env(force_env_);
+        if (GMX_LIB_MPI)
+        {
+            cp2k_finalize_without_mpi();
+        }
+        else
+        {
+            cp2k_finalize();
+        }
+    }
+}
+
+bool QMMMForceProvider::isQMAtom(index globalAtomIndex)
+{
+    return std::find(qmAtoms_.globalIndex().begin(), qmAtoms_.globalIndex().end(), globalAtomIndex)
+           != qmAtoms_.globalIndex().end();
+}
+
+void QMMMForceProvider::appendLog(const std::string& msg)
+{
+    GMX_LOG(logger_.info).asParagraph().appendText(msg);
+}
+
+void QMMMForceProvider::initCP2KForceEnvironment(const t_commrec& cr)
+{
+    // Check that we have filename either defined in KVT or deduced from *.tpr name
+    GMX_RELEASE_ASSERT(parameters_.qmFileNameBase_.empty(),
+                       "Names of CP2K input/output files is not defined in QMMMParameters");
+
+
+    // Names of CP2K Input and Output files
+    const std::string cp2kInputName  = parameters_.qmFileNameBase_ + ".inp";
+    const std::string cp2kPdbName    = parameters_.qmFileNameBase_ + ".pdb";
+    const std::string cp2kOutputName = parameters_.qmFileNameBase_ + ".out";
+
+    // Write CP2K input if we are Master
+    if (MASTER(&cr))
+    {
+        // In the CP2K Input we need to substitute placeholder with the actuall *.pdb file name
+        TextWriter::writeFileFromString(
+                cp2kInputName, formatString(parameters_.qmInput_.c_str(), cp2kPdbName.c_str()));
+
+        // Write *.pdb with point charges for CP2K
+        TextWriter::writeFileFromString(cp2kPdbName, parameters_.qmPdb_);
+    }
+
+    // Check if we have external MPI library
+    if (GMX_LIB_MPI)
+    {
+        /* Attempt to init CP2K and create force environment in case we have an external MPI library.
+         * _without_mpi - means that libcp2k will not call MPI_Init() itself,
+         * but rather expects it to be called already
+         */
+        cp2k_init_without_mpi();
+
+        // Here #if (GMX_LIB_MPI) needed because of MPI_Comm_c2f() usage
+        // TODO: Probably there should be more elegant solution
+#if GMX_LIB_MPI
+        cp2k_create_force_env_comm(
+                &force_env_, cp2kInputName.c_str(), cp2kOutputName.c_str(), MPI_Comm_c2f(cr.mpi_comm_mysim));
+#endif
+    }
+    else
+    {
+
+        // If we have thread-MPI or no-MPI then we should initialize CP2P differently
+        if (cr.nnodes > 1)
+        {
+            // CP2K could not use thread-MPI parallelization. In case -ntmpi > 1 throw an error.
+            std::string msg =
+                    "Using GROMACS with CP2K requires the GROMACS build to be configured "
+                    "with an MPI library or to use a single thread-MPI rank (-ntmpi 1). "
+                    "In the latter case, manual use of -ntomp is also advisable when "
+                    "the node has many cores to fill them with threads.";
+            GMX_THROW(NotImplementedError(msg.c_str()));
+        }
+
+        // Attempt to init CP2K and create force environment
+        cp2k_init();
+        cp2k_create_force_env(&force_env_, cp2kInputName.c_str(), cp2kOutputName.c_str());
+    }
+
+    // Set flag of successful initialization
+    isCp2kLibraryInitialized_ = true;
+
+} // namespace gmx
+
+void QMMMForceProvider::calculateForces(const ForceProviderInput& fInput, ForceProviderOutput* fOutput)
+{
+    // Total number of atoms in the system
+    size_t numAtoms = qmAtoms_.numAtomsGlobal() + mmAtoms_.numAtomsGlobal();
+
+    // Save box
+    copy_mat(fInput.box_, box_);
+
+    // Initialize PBC
+    t_pbc pbc;
+    set_pbc(&pbc, pbcType_, box_);
+
+    /*
+     * 1) If calculateForce called first time, then we need to init CP2K,
+     *    as it was not possible during QMMMForceProvider constructor
+     *    due to absence of full communication record at that point.
+     */
+    if (!isCp2kLibraryInitialized_)
+    {
+        try
+        {
+            initCP2KForceEnvironment(fInput.cr_);
+        }
+        GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
+    }
+
+    /*
+     * 2) We need to gather fInput.x_ in case of MPI / DD setup
+     */
+
+    // x - coordinates (gathered across nodes in case of DD)
+    std::vector<RVec> x(numAtoms, RVec({ 0.0, 0.0, 0.0 }));
+
+    // Fill cordinates of local QM atoms and add translation
+    for (size_t i = 0; i < qmAtoms_.numAtomsLocal(); i++)
+    {
+        x[qmAtoms_.globalIndex()[qmAtoms_.collectiveIndex()[i]]] =
+                fInput.x_[qmAtoms_.localIndex()[i]] + parameters_.qmTrans_;
+    }
+
+    // Fill cordinates of local MM atoms and add translation
+    for (size_t i = 0; i < mmAtoms_.numAtomsLocal(); i++)
+    {
+        x[mmAtoms_.globalIndex()[mmAtoms_.collectiveIndex()[i]]] =
+                fInput.x_[mmAtoms_.localIndex()[i]] + parameters_.qmTrans_;
+    }
+
+    // If we are in MPI / DD conditions then gather coordinates over nodes
+    if (havePPDomainDecomposition(&fInput.cr_))
+    {
+        gmx_sum(3 * numAtoms, x.data()->as_vec(), &fInput.cr_);
+    }
+
+    // Put all atoms into the central box (they might be shifted out of it because of the translation)
+    put_atoms_in_box(pbcType_, fInput.box_, ArrayRef<RVec>(x));
+
+    /*
+     * 3) Cast data to double format of libcp2k
+     *    update coordinates and box in CP2K and perform QM calculation
+     */
+    // x_d - coordinates casted to linear dobule vector for CP2K with parameters_.qmTrans_ added
+    std::vector<double> x_d(3 * numAtoms, 0.0);
+    for (size_t i = 0; i < numAtoms; i++)
+    {
+        x_d[3 * i]     = static_cast<double>((x[i][XX]) / c_bohr2Nm);
+        x_d[3 * i + 1] = static_cast<double>((x[i][YY]) / c_bohr2Nm);
+        x_d[3 * i + 2] = static_cast<double>((x[i][ZZ]) / c_bohr2Nm);
+    }
+
+    // box_d - box_ casted to linear dobule vector for CP2K
+    std::vector<double> box_d(9);
+    for (size_t i = 0; i < DIM; i++)
+    {
+        box_d[3 * i]     = static_cast<double>(box_[0][i] / c_bohr2Nm);
+        box_d[3 * i + 1] = static_cast<double>(box_[1][i] / c_bohr2Nm);
+        box_d[3 * i + 2] = static_cast<double>(box_[2][i] / c_bohr2Nm);
+    }
+
+    // Update coordinates and box in CP2K
+    cp2k_set_positions(force_env_, x_d.data(), 3 * numAtoms);
+    cp2k_set_cell(force_env_, box_d.data());
+
+    // Run CP2K calculation
+    cp2k_calc_energy_force(force_env_);
+
+    /*
+     * 4) Get output data
+     * We need to fill only local part into fOutput
+     */
+
+    // Only master process should add QM + QMMM energy
+    if (MASTER(&fInput.cr_))
+    {
+        double qmEner = 0.0;
+        cp2k_get_potential_energy(force_env_, &qmEner);
+        fOutput->enerd_.term[F_EQM] += qmEner * c_hartree2Kj * c_avogadro;
+    }
+
+    // Get Forces they are in Hartree/Bohr and will be converted to kJ/mol/nm
+    std::vector<double> cp2kForce(3 * numAtoms, 0.0);
+    cp2k_get_forces(force_env_, cp2kForce.data(), 3 * numAtoms);
+
+    // Fill forces on QM atoms first
+    for (size_t i = 0; i < qmAtoms_.numAtomsLocal(); i++)
+    {
+        fOutput->forceWithVirial_.force_[qmAtoms_.localIndex()[i]][XX] +=
+                static_cast<real>(cp2kForce[3 * qmAtoms_.globalIndex()[qmAtoms_.collectiveIndex()[i]]])
+                * c_hartreeBohr2Md;
+
+        fOutput->forceWithVirial_.force_[qmAtoms_.localIndex()[i]][YY] +=
+                static_cast<real>(cp2kForce[3 * qmAtoms_.globalIndex()[qmAtoms_.collectiveIndex()[i]] + 1])
+                * c_hartreeBohr2Md;
+
+        fOutput->forceWithVirial_.force_[qmAtoms_.localIndex()[i]][ZZ] +=
+                static_cast<real>(cp2kForce[3 * qmAtoms_.globalIndex()[qmAtoms_.collectiveIndex()[i]] + 2])
+                * c_hartreeBohr2Md;
+    }
+
+    // Filll forces on MM atoms then
+    for (size_t i = 0; i < mmAtoms_.numAtomsLocal(); i++)
+    {
+        fOutput->forceWithVirial_.force_[mmAtoms_.localIndex()[i]][XX] +=
+                static_cast<real>(cp2kForce[3 * mmAtoms_.globalIndex()[mmAtoms_.collectiveIndex()[i]]])
+                * c_hartreeBohr2Md;
+
+        fOutput->forceWithVirial_.force_[mmAtoms_.localIndex()[i]][YY] +=
+                static_cast<real>(cp2kForce[3 * mmAtoms_.globalIndex()[mmAtoms_.collectiveIndex()[i]] + 1])
+                * c_hartreeBohr2Md;
+
+        fOutput->forceWithVirial_.force_[mmAtoms_.localIndex()[i]][ZZ] +=
+                static_cast<real>(cp2kForce[3 * mmAtoms_.globalIndex()[mmAtoms_.collectiveIndex()[i]] + 2])
+                * c_hartreeBohr2Md;
+    }
+};
+
+} // namespace gmx
diff --git a/src/gromacs/applied_forces/qmmm/qmmmforceprovider.h b/src/gromacs/applied_forces/qmmm/qmmmforceprovider.h
new file mode 100644 (file)
index 0000000..e90565c
--- /dev/null
@@ -0,0 +1,123 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Declares force provider for QMMM
+ *
+ * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
+ * \author Christian Blau <blau@kth.se>
+ * \ingroup module_applied_forces
+ */
+#ifndef GMX_APPLIED_FORCES_QMMMFORCEPROVIDER_H
+#define GMX_APPLIED_FORCES_QMMMFORCEPROVIDER_H
+
+#include "gromacs/domdec/localatomset.h"
+#include "gromacs/mdtypes/forceoutput.h"
+#include "gromacs/mdtypes/iforceprovider.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/utility/classhelpers.h"
+#include "gromacs/utility/logger.h"
+
+#include "qmmmtypes.h"
+
+namespace gmx
+{
+
+#ifdef __clang__
+#    pragma clang diagnostic push
+#    pragma clang diagnostic ignored "-Wunused-private-field"
+#endif
+
+//! Type for CP2K force environment handle
+typedef int force_env_t;
+
+/*! \internal \brief
+ * Implements IForceProvider for QM/MM.
+ */
+class QMMMForceProvider final : public IForceProvider
+{
+public:
+    QMMMForceProvider(const QMMMParameters& parameters,
+                      const LocalAtomSet&   localQMAtomSet,
+                      const LocalAtomSet&   localMMAtomSet,
+                      PbcType               pbcType,
+                      const MDLogger&       logger);
+
+    //! Destruct force provider for QMMM and finalize libcp2k
+    ~QMMMForceProvider();
+
+    /*!\brief Calculate forces of QMMM.
+     * \param[in] fInput input for force provider
+     * \param[out] fOutput output for force provider
+     */
+    void calculateForces(const ForceProviderInput& fInput, ForceProviderOutput* fOutput) override;
+
+private:
+    //! Write message to the log
+    void appendLog(const std::string& msg);
+
+    /*!\brief Check if atom belongs to the global index of qmAtoms_
+     * \param[in] globalAtomIndex global index of the atom to check
+     */
+    bool isQMAtom(index globalAtomIndex);
+
+    /*!\brief Initialization of QM program.
+     * \param[in] cr connection record structure
+     */
+    void initCP2KForceEnvironment(const t_commrec& cr);
+
+    const QMMMParameters& parameters_;
+    const LocalAtomSet&   qmAtoms_;
+    const LocalAtomSet&   mmAtoms_;
+    const PbcType         pbcType_;
+    const MDLogger&       logger_;
+
+    //! Internal copy of PBC box
+    matrix box_;
+
+    //! Flag wether initCP2KForceEnvironment() has been called already
+    bool isCp2kLibraryInitialized_ = false;
+
+    //! CP2K force environment handle
+    force_env_t force_env_ = -1;
+};
+
+#ifdef __clang__
+#    pragma clang diagnostic pop
+#endif
+
+} // namespace gmx
+
+#endif // GMX_APPLIED_FORCES_QMMMFORCEPROVIDER_H
diff --git a/src/gromacs/applied_forces/qmmm/qmmmforceprovider_stub.cpp b/src/gromacs/applied_forces/qmmm/qmmmforceprovider_stub.cpp
new file mode 100644 (file)
index 0000000..c45c0f0
--- /dev/null
@@ -0,0 +1,114 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Stub implementation of QMMMForceProvider
+ * Compiled in case if CP2K is not linked
+ *
+ * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
+ * \author Christian Blau <blau@kth.se>
+ * \ingroup module_applied_forces
+ */
+
+#include "gmxpre.h"
+
+#include "qmmmforceprovider.h"
+
+#include "gromacs/utility/exceptions.h"
+
+namespace gmx
+{
+
+
+#ifdef __clang__
+#    pragma clang diagnostic push
+#    pragma clang diagnostic ignored "-Wmissing-noreturn"
+#endif
+
+QMMMForceProvider::QMMMForceProvider(const QMMMParameters& parameters,
+                                     const LocalAtomSet&   localQMAtomSet,
+                                     const LocalAtomSet&   localMMAtomSet,
+                                     PbcType               pbcType,
+                                     const MDLogger&       logger) :
+    parameters_(parameters),
+    qmAtoms_(localQMAtomSet),
+    mmAtoms_(localMMAtomSet),
+    pbcType_(pbcType),
+    logger_(logger),
+    box_{ { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } }
+{
+    GMX_THROW(
+            InternalError("CP2K has not been linked into GROMACS, QMMM simulation is not "
+                          "possible.\nPlease, reconfigure GROMACS with -DGMX_CP2K=ON\n"));
+}
+
+QMMMForceProvider::~QMMMForceProvider() {}
+
+// NOLINTNEXTLINE(readability-convert-member-functions-to-static)
+bool QMMMForceProvider::isQMAtom(index /*globalAtomIndex*/)
+{
+    GMX_THROW(
+            InternalError("CP2K has not been linked into GROMACS, QMMM simulation is not "
+                          "possible.\nPlease, reconfigure GROMACS with -DGMX_CP2K=ON\n"));
+}
+
+// NOLINTNEXTLINE(readability-convert-member-functions-to-static)
+void QMMMForceProvider::appendLog(const std::string& /*msg*/)
+{
+    GMX_THROW(
+            InternalError("CP2K has not been linked into GROMACS, QMMM simulation is not "
+                          "possible.\nPlease, reconfigure GROMACS with -DGMX_CP2K=ON\n"));
+}
+
+// NOLINTNEXTLINE(readability-convert-member-functions-to-static)
+void QMMMForceProvider::initCP2KForceEnvironment(const t_commrec& /*cr*/)
+{
+    GMX_THROW(
+            InternalError("CP2K has not been linked into GROMACS, QMMM simulation is not "
+                          "possible.\nPlease, reconfigure GROMACS with -DGMX_CP2K=ON\n"));
+}
+
+void QMMMForceProvider::calculateForces(const ForceProviderInput& /*fInput*/, ForceProviderOutput* /*fOutput*/)
+{
+    GMX_THROW(
+            InternalError("CP2K has not been linked into GROMACS, QMMM simulation is not "
+                          "possible.\nPlease, reconfigure GROMACS with -DGMX_CP2K=ON\n"));
+};
+
+#ifdef __clang__
+#    pragma clang diagnostic pop
+#endif
+
+} // namespace gmx
index beaf8b144168cb95d991bb0063d2cc18351a2c1c..37eebc77f887a57f117da1bcdc69de07680b666d 100644 (file)
@@ -86,9 +86,16 @@ static const EnumerationArray<QMMMQMMethod, const char*> c_qmmmQMMethodNames = {
 
 //! symbols of the elements in periodic table
 const std::vector<std::string> periodic_system = {
-    "X  ", "H  ", "He ", "Li ", "Be ", "B  ", "C  ", "N  ", "O  ", "F  ", "Ne ", "Na ", "Mg ",
-    "Al ", "Si ", "P  ", "S  ", "Cl ", "Ar ", "K  ", "Ca ", "Sc ", "Ti ", "V  ", "Cr ", "Mn ",
-    "Fe ", "Co ", "Ni ", "Cu ", "Zn ", "Ga ", "Ge ", "As ", "Se ", "Br ", "Kr "
+    "X  ", "H  ", "He ", "Li ", "Be ", "B  ", "C  ", "N  ", "O  ", "F  ", "Ne ", "Na ",
+    "Mg ", "Al ", "Si ", "P  ", "S  ", "Cl ", "Ar ", "K  ", "Ca ", "Sc ", "Ti ", "V  ",
+    "Cr ", "Mn ", "Fe ", "Co ", "Ni ", "Cu ", "Zn ", "Ga ", "Ge ", "As ", "Se ", "Br ",
+    "Kr ", "Rb ", "Sr ", "Y  ", "Zr ", "Nb ", "Mo ", "Tc ", "Ru ", "Rh ", "Pd ", "Ag ",
+    "Cd ", "In ", "Sn ", "Sb ", "Te ", "I  ", "Xe ", "Cs ", "Ba ", "La ", "Ce ", "Pr ",
+    "Nd ", "Pm ", "Sm ", "Eu ", "Gd ", "Tb ", "Dy ", "Ho ", "Er ", "Tm ", "Yb ", "Lu ",
+    "Hf ", "Ta ", "W  ", "Re ", "Os ", "Ir ", "Pt ", "Au ", "Hg ", "Tl ", "Pb ", "Bi ",
+    "Po ", "At ", "Rn ", "Fr ", "Ra ", "Ac ", "Th ", "Pa ", "U  ", "Np ", "Pu ", "Am ",
+    "Cm ", "Bk ", "Cf ", "Es ", "Fm ", "Md ", "No ", "Lr ", "Rf ", "Db ", "Sg ", "Bh ",
+    "Hs ", "Mt ", "Ds ", "Rg ", "Cn ", "Nh ", "Fl ", "Mc ", "Lv ", "Ts ", "Og "
 };
 
 /*! \internal
@@ -113,17 +120,17 @@ struct QMMMParameters
     int qmMult_ = 1;
     //! Method used for QM calculation (default DFT with PBE functional)
     QMMMQMMethod qmMethod_ = QMMMQMMethod::PBE;
-    //! String containing name of the input file for CP2K (default "topol-qmmm.inp")
-    std::string qmInputFileName_ = "topol-qmmm.inp";
-    //! String containing name of the PDB file for CP2K input (default "topol-qmmm.pdb")
-    std::string qmPdbFileName_ = "topol-qmmm.pdb";
+    /*! \brief String containing name of the CP2K files (*.inp, *.out, *.pdb)
+     * default value empty, means will be deduced from *.tpr name during mdrun
+     */
+    std::string qmFileNameBase_;
     //! String containing whole CP2K input which can be stored inside *.tpr
     std::string qmInput_;
     //! String containing PDB file for CP2K input which can be stored inside *.tpr
     std::string qmPdb_;
     //! Matrix that contains vectors defining QM box
     matrix qmBox_;
-    //! Matrix that contains vectors defining QM box
+    //! Translation vector to center QM subsystem inside the QM Box
     RVec qmTrans_;
 
     //! Constructor with default initializers for arrays
index 7c879e151a68606a69a03282ce44a70e3605c733..230b8385c6389773e8e05e02eca42a61ac69aced 100644 (file)
@@ -36,4 +36,12 @@ gmx_add_unit_test(QMMMAppliedForcesUnitTest qmmm_applied_forces-test
     CPP_SOURCE_FILES
         qmmminputgenerator.cpp
         qmmmtopologypreprocessor.cpp
+        qmmmforceprovider.cpp
         )
+
+# For QMMMForceProvider test we need to define wether libcp2k linked
+if (GMX_CP2K)
+    set_source_files_properties(qmmmforceprovider.cpp PROPERTIES COMPILE_DEFINITIONS GMX_CP2K=1)
+else()
+    set_source_files_properties(qmmmforceprovider.cpp PROPERTIES COMPILE_DEFINITIONS GMX_CP2K=0)
+endif()
diff --git a/src/gromacs/applied_forces/qmmm/tests/qmmmforceprovider.cpp b/src/gromacs/applied_forces/qmmm/tests/qmmmforceprovider.cpp
new file mode 100644 (file)
index 0000000..6d648c2
--- /dev/null
@@ -0,0 +1,111 @@
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 2021, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+/*! \internal \file
+ * \brief
+ * Tests for QMMMForceProvider class for QMMM MDModule
+ *
+ * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
+ * \ingroup module_applied_forces
+ */
+#include "gmxpre.h"
+
+#include "gromacs/applied_forces/qmmm/qmmmforceprovider.h"
+
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/domdec/localatomsetmanager.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/gmxpreprocess/grompp.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/mtop_lookup.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/path.h"
+#include "gromacs/utility/textwriter.h"
+#include "testutils/cmdlinetest.h"
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+#include "testutils/testfilemanager.h"
+
+namespace gmx
+{
+
+class QMMMForceProviderTest : public ::testing::Test
+{
+public:
+    void setDefaultParameters()
+    {
+        parameters_.active_                = true;
+        std::vector<gmx::index> qmIndicies = { 0, 1, 2 };
+        std::vector<gmx::index> mmIndicies = { 3, 4, 5 };
+        LocalAtomSet            set1 = atomSetManager_.add(ArrayRef<const gmx::index>(qmIndicies));
+        LocalAtomSet            set2 = atomSetManager_.add(ArrayRef<const gmx::index>(mmIndicies));
+        qmAtomSet_                   = std::make_unique<LocalAtomSet>(set1);
+        mmAtomSet_                   = std::make_unique<LocalAtomSet>(set2);
+    }
+
+protected:
+    QMMMParameters                parameters_;
+    LocalAtomSetManager           atomSetManager_;
+    std::unique_ptr<LocalAtomSet> qmAtomSet_;
+    std::unique_ptr<LocalAtomSet> mmAtomSet_;
+    PbcType                       pbcType_;
+    MDLogger                      logger_;
+};
+
+TEST_F(QMMMForceProviderTest, CanConstructOrNot)
+{
+    setDefaultParameters();
+
+    // GMX_CP2K is defined in CMakeList.txt trough set_source_files_properties()
+    if (GMX_CP2K)
+    {
+        // if libcp2k linked then we do not expect throws
+        EXPECT_NO_THROW(QMMMForceProvider forceProvider(
+                parameters_, *qmAtomSet_, *mmAtomSet_, pbcType_, logger_));
+    }
+    else
+    {
+        // if libcp2k not linked then we expect throw from constructor
+        EXPECT_ANY_THROW(QMMMForceProvider forceProvider(
+                parameters_, *qmAtomSet_, *mmAtomSet_, pbcType_, logger_));
+    }
+}
+
+} // namespace gmx