From: Dmitry Morozov Date: Thu, 23 Sep 2021 08:03:08 +0000 (+0000) Subject: Added QMMMForceProvider class for QMMM MdModule X-Git-Url: http://biod.pnpi.spb.ru/gitweb/?a=commitdiff_plain;h=21fa29c865259b3e480b9f73370d847615177451;p=alexxy%2Fgromacs.git Added QMMMForceProvider class for QMMM MdModule --- diff --git a/cmake/gmxManageCP2K.cmake b/cmake/gmxManageCP2K.cmake index 99f5b70c04..d23ca3392e 100644 --- a/cmake/gmxManageCP2K.cmake +++ b/cmake/gmxManageCP2K.cmake @@ -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") diff --git a/src/gromacs/applied_forces/qmmm/CMakeLists.txt b/src/gromacs/applied_forces/qmmm/CMakeLists.txt index ba20be3824..3707aa1188 100644 --- a/src/gromacs/applied_forces/qmmm/CMakeLists.txt +++ b/src/gromacs/applied_forces/qmmm/CMakeLists.txt @@ -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 index 0000000000..f685efdcf7 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmmforceprovider.cpp @@ -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 + * \author Christian Blau + * \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 + +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 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(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 x_d(3 * numAtoms, 0.0); + for (size_t i = 0; i < numAtoms; i++) + { + x_d[3 * i] = static_cast((x[i][XX]) / c_bohr2Nm); + x_d[3 * i + 1] = static_cast((x[i][YY]) / c_bohr2Nm); + x_d[3 * i + 2] = static_cast((x[i][ZZ]) / c_bohr2Nm); + } + + // box_d - box_ casted to linear dobule vector for CP2K + std::vector box_d(9); + for (size_t i = 0; i < DIM; i++) + { + box_d[3 * i] = static_cast(box_[0][i] / c_bohr2Nm); + box_d[3 * i + 1] = static_cast(box_[1][i] / c_bohr2Nm); + box_d[3 * i + 2] = static_cast(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 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(cp2kForce[3 * qmAtoms_.globalIndex()[qmAtoms_.collectiveIndex()[i]]]) + * c_hartreeBohr2Md; + + fOutput->forceWithVirial_.force_[qmAtoms_.localIndex()[i]][YY] += + static_cast(cp2kForce[3 * qmAtoms_.globalIndex()[qmAtoms_.collectiveIndex()[i]] + 1]) + * c_hartreeBohr2Md; + + fOutput->forceWithVirial_.force_[qmAtoms_.localIndex()[i]][ZZ] += + static_cast(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(cp2kForce[3 * mmAtoms_.globalIndex()[mmAtoms_.collectiveIndex()[i]]]) + * c_hartreeBohr2Md; + + fOutput->forceWithVirial_.force_[mmAtoms_.localIndex()[i]][YY] += + static_cast(cp2kForce[3 * mmAtoms_.globalIndex()[mmAtoms_.collectiveIndex()[i]] + 1]) + * c_hartreeBohr2Md; + + fOutput->forceWithVirial_.force_[mmAtoms_.localIndex()[i]][ZZ] += + static_cast(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 index 0000000000..e90565c446 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmmforceprovider.h @@ -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 + * \author Christian Blau + * \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 index 0000000000..c45c0f043b --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmmforceprovider_stub.cpp @@ -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 + * \author Christian Blau + * \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 diff --git a/src/gromacs/applied_forces/qmmm/qmmmtypes.h b/src/gromacs/applied_forces/qmmm/qmmmtypes.h index beaf8b1441..37eebc77f8 100644 --- a/src/gromacs/applied_forces/qmmm/qmmmtypes.h +++ b/src/gromacs/applied_forces/qmmm/qmmmtypes.h @@ -86,9 +86,16 @@ static const EnumerationArray c_qmmmQMMethodNames = { //! symbols of the elements in periodic table const std::vector 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 diff --git a/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt b/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt index 7c879e151a..230b8385c6 100644 --- a/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt +++ b/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt @@ -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 index 0000000000..6d648c2513 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/qmmmforceprovider.cpp @@ -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 + * \ingroup module_applied_forces + */ +#include "gmxpre.h" + +#include "gromacs/applied_forces/qmmm/qmmmforceprovider.h" + +#include + +#include + +#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 qmIndicies = { 0, 1, 2 }; + std::vector mmIndicies = { 3, 4, 5 }; + LocalAtomSet set1 = atomSetManager_.add(ArrayRef(qmIndicies)); + LocalAtomSet set2 = atomSetManager_.add(ArrayRef(mmIndicies)); + qmAtomSet_ = std::make_unique(set1); + mmAtomSet_ = std::make_unique(set2); + } + +protected: + QMMMParameters parameters_; + LocalAtomSetManager atomSetManager_; + std::unique_ptr qmAtomSet_; + std::unique_ptr 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