From c8170794de75246116850afcd3a752509327de8b Mon Sep 17 00:00:00 2001 From: Dmitry Morozov Date: Thu, 12 Aug 2021 09:40:38 +0000 Subject: [PATCH] Added class for generating standard CP2K input --- src/gromacs/applied_forces/CMakeLists.txt | 3 +- .../applied_forces/qmmm/CMakeLists.txt | 41 ++ .../qmmm/qmmminputgenerator.cpp | 507 ++++++++++++++++++ .../applied_forces/qmmm/qmmminputgenerator.h | 160 ++++++ src/gromacs/applied_forces/qmmm/qmmmtypes.h | 140 +++++ .../applied_forces/qmmm/tests/.clang-tidy | 9 + .../applied_forces/qmmm/tests/CMakeLists.txt | 38 ++ .../qmmm/tests/qmmminputgenerator.cpp | 147 +++++ ...MInputGeneratorTest_TwoWatersPBENoLink.xml | 149 +++++ ...nputGeneratorTest_TwoWatersPBEWithLink.xml | 153 ++++++ 10 files changed, 1346 insertions(+), 1 deletion(-) create mode 100644 src/gromacs/applied_forces/qmmm/CMakeLists.txt create mode 100644 src/gromacs/applied_forces/qmmm/qmmminputgenerator.cpp create mode 100644 src/gromacs/applied_forces/qmmm/qmmminputgenerator.h create mode 100644 src/gromacs/applied_forces/qmmm/qmmmtypes.h create mode 100644 src/gromacs/applied_forces/qmmm/tests/.clang-tidy create mode 100644 src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt create mode 100644 src/gromacs/applied_forces/qmmm/tests/qmmminputgenerator.cpp create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBENoLink.xml create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBEWithLink.xml diff --git a/src/gromacs/applied_forces/CMakeLists.txt b/src/gromacs/applied_forces/CMakeLists.txt index 92e40d3895..5397c64584 100644 --- a/src/gromacs/applied_forces/CMakeLists.txt +++ b/src/gromacs/applied_forces/CMakeLists.txt @@ -1,7 +1,7 @@ # # This file is part of the GROMACS molecular simulation package. # -# Copyright (c) 2015,2016,2019,2020, by the GROMACS development team, led by +# Copyright (c) 2015,2016,2019,2020,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. @@ -73,6 +73,7 @@ gmx_add_libgromacs_sources( add_subdirectory(awh) add_subdirectory(densityfitting) +add_subdirectory(qmmm) if (BUILD_TESTING) add_subdirectory(tests) diff --git a/src/gromacs/applied_forces/qmmm/CMakeLists.txt b/src/gromacs/applied_forces/qmmm/CMakeLists.txt new file mode 100644 index 0000000000..c64cadff3f --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/CMakeLists.txt @@ -0,0 +1,41 @@ +# +# 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. + +gmx_add_libgromacs_sources( + qmmminputgenerator.cpp + ) + +if (BUILD_TESTING) + add_subdirectory(tests) +endif() diff --git a/src/gromacs/applied_forces/qmmm/qmmminputgenerator.cpp b/src/gromacs/applied_forces/qmmm/qmmminputgenerator.cpp new file mode 100644 index 0000000000..3bcca1ceb2 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmminputgenerator.cpp @@ -0,0 +1,507 @@ +/* + * 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 input generator class for CP2K QMMM + * + * \author Dmitry Morozov + * \ingroup module_applied_forces + */ + +#include "gmxpre.h" + +#include "qmmminputgenerator.h" + +#include "gromacs/math/units.h" +#include "gromacs/math/vec.h" +#include "gromacs/utility/stringutil.h" + +namespace gmx +{ + +QMMMInputGenerator::QMMMInputGenerator(const QMMMParameters& parameters, + PbcType pbcType, + const matrix box, + ArrayRef q, + ArrayRef x) : + parameters_(parameters), + pbc_(pbcType), + qmBox_{ { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } }, + qmCenter_{ 0.0, 0.0, 0.0 }, + qmTrans_{ 0.0, 0.0, 0.0 }, + q_(q), + x_(x) +{ + // Copy box + copy_mat(box, box_); + + // Fill qmAtoms_ set + for (const auto& val : parameters_.qmIndices_) + { + qmAtoms_.emplace(val); + } + + // Compute QM box + computeQMBox(sc_qmBoxScale, sc_qmBoxMinLength); +} + +bool QMMMInputGenerator::isQMAtom(index globalAtomIndex) const +{ + return (qmAtoms_.find(globalAtomIndex) != qmAtoms_.end()); +} + +void QMMMInputGenerator::computeQMBox(real scale, real minNorm) +{ + // Init atom numbers + size_t nQm = parameters_.qmIndices_.size(); + + // If there is only one QM atom, then just copy the box_ + if (nQm < 2) + { + copy_mat(box_, qmBox_); + qmCenter_ = x_[parameters_.qmIndices_[0]]; + qmTrans_ = RVec(qmBox_[0]) / 2 + RVec(qmBox_[1]) / 2 + RVec(qmBox_[2]) / 2 - qmCenter_; + return; + } + + // Initialize pbc + t_pbc pbc; + set_pbc(&pbc, pbc_, box_); + + /* To compute qmBox_: + * 1) Compute maximum dimension - maxDist betweeen QM atoms within PBC + * 2) Make projection of the each box_ vector onto the cross prod of other two vectors + * 3) Calculate scales so that norm will be maxDist for each box_ vector + * 4) Apply scales to the box_ to get qmBox_ + */ + RVec dx(0.0, 0.0, 0.0); + real maxDist = 0.0; + + // Search for the maxDist - maximum distance within QM system + for (size_t i = 0; i < nQm - 1; i++) + { + for (size_t j = i + 1; j < nQm; j++) + { + pbc_dx(&pbc, x_[parameters_.qmIndices_[i]], x_[parameters_.qmIndices_[j]], dx); + maxDist = dx.norm() > maxDist ? dx.norm() : maxDist; + } + } + + // Apply scale factor: qmBox_ should be *scale times bigger than maxDist + maxDist *= scale; + + // Compute QM Box vectors + RVec vec0 = computeQMBoxVec(box_[0], box_[1], box_[2], maxDist, minNorm, norm(box_[0])); + RVec vec1 = computeQMBoxVec(box_[1], box_[0], box_[2], maxDist, minNorm, norm(box_[1])); + RVec vec2 = computeQMBoxVec(box_[2], box_[0], box_[1], maxDist, minNorm, norm(box_[2])); + + // Form qmBox_ + copy_rvec(vec0, qmBox_[0]); + copy_rvec(vec1, qmBox_[1]); + copy_rvec(vec2, qmBox_[2]); + + /* Now we need to also compute translation vector. + * In order to center QM atoms in the computed qmBox_ + * + * First compute center of QM system by averaging PBC-aware distance vectors + * with respect to the first QM atom. + */ + for (size_t i = 1; i < nQm; i++) + { + // compute pbc-aware distance vector between QM atom 0 and QM atom i + pbc_dx(&pbc, x_[parameters_.qmIndices_[i]], x_[parameters_.qmIndices_[0]], dx); + + // add that to the qmCenter_ + qmCenter_ += dx; + } + + // Average over all nQm atoms and add first atom coordiantes + qmCenter_ = qmCenter_ / nQm + x_[parameters_.qmIndices_[0]]; + + // Translation vector will be center of qmBox_ - qmCenter_ + qmTrans_ = vec0 / 2 + vec1 / 2 + vec2 / 2 - qmCenter_; +} + +std::string QMMMInputGenerator::generateGlobalSection() +{ + std::string res; + + res += "&GLOBAL\n"; + res += " PRINT_LEVEL LOW\n"; + res += " PROJECT GROMACS\n"; + res += " RUN_TYPE ENERGY_FORCE\n"; + res += "&END GLOBAL\n"; + + return res; +} + +std::string QMMMInputGenerator::generateDFTSection() const +{ + std::string res; + + res += " &DFT\n"; + + // write charge and multiplicity + res += formatString(" CHARGE %d\n", parameters_.qmCharge_); + res += formatString(" MULTIPLICITY %d\n", parameters_.qmMult_); + + // If multiplicity is not 1 then we should use unrestricted Kohn-Sham + if (parameters_.qmMult_ > 1) + { + res += " UKS\n"; + } + + // Basis files, Grid setup and SCF parameters + res += " BASIS_SET_FILE_NAME BASIS_MOLOPT\n"; + res += " POTENTIAL_FILE_NAME POTENTIAL\n"; + res += " &MGRID\n"; + res += " NGRIDS 5\n"; + res += " CUTOFF 450\n"; + res += " REL_CUTOFF 50\n"; + res += " COMMENSURATE\n"; + res += " &END MGRID\n"; + res += " &SCF\n"; + res += " SCF_GUESS RESTART\n"; + res += " EPS_SCF 5.0E-8\n"; + res += " MAX_SCF 20\n"; + res += " &OT T\n"; + res += " MINIMIZER DIIS\n"; + res += " STEPSIZE 0.15\n"; + res += " PRECONDITIONER FULL_ALL\n"; + res += " &END OT\n"; + res += " &OUTER_SCF T\n"; + res += " MAX_SCF 20\n"; + res += " EPS_SCF 5.0E-8\n"; + res += " &END OUTER_SCF\n"; + res += " &END SCF\n"; + + // DFT functional parameters + res += " &XC\n"; + res += " DENSITY_CUTOFF 1.0E-12\n"; + res += " GRADIENT_CUTOFF 1.0E-12\n"; + res += " TAU_CUTOFF 1.0E-12\n"; + res += formatString(" &XC_FUNCTIONAL %s\n", c_qmmmQMMethodNames[parameters_.qmMethod_]); + res += " &END XC_FUNCTIONAL\n"; + res += " &END XC\n"; + res += " &QS\n"; + res += " METHOD GPW\n"; + res += " EPS_DEFAULT 1.0E-10\n"; + res += " EXTRAPOLATION ASPC\n"; + res += " EXTRAPOLATION_ORDER 4\n"; + res += " &END QS\n"; + res += " &END DFT\n"; + + return res; +} + +std::string QMMMInputGenerator::generateQMMMSection() const +{ + std::string res; + + // Init some numbers + const size_t nQm = parameters_.qmIndices_.size(); + + // Count the numbers of individual QM atoms per type + std::vector num_atoms(periodic_system.size(), 0); + for (size_t i = 0; i < nQm; i++) + { + num_atoms[parameters_.atomNumbers_[parameters_.qmIndices_[i]]]++; + } + + res += " &QMMM\n"; + + // QM cell + res += " &CELL\n"; + res += formatString( + " A %.3lf %.3lf %.3lf\n", qmBox_[0][0] * 10, qmBox_[0][1] * 10, qmBox_[0][2] * 10); + res += formatString( + " B %.3lf %.3lf %.3lf\n", qmBox_[1][0] * 10, qmBox_[1][1] * 10, qmBox_[1][2] * 10); + res += formatString( + " C %.3lf %.3lf %.3lf\n", qmBox_[2][0] * 10, qmBox_[2][1] * 10, qmBox_[2][2] * 10); + res += " PERIODIC XYZ\n"; + res += " &END CELL\n"; + + res += " CENTER EVERY_STEP\n"; + res += " CENTER_GRID TRUE\n"; + res += " &WALLS\n"; + res += " TYPE REFLECTIVE\n"; + res += " &END WALLS\n"; + + res += " ECOUPL GAUSS\n"; + res += " USE_GEEP_LIB 12\n"; + res += " &PERIODIC\n"; + res += " GMAX 1.0E+00\n"; + res += " &MULTIPOLE ON\n"; + res += " RCUT 1.0E+01\n"; + res += " EWALD_PRECISION 1.0E-06\n"; + res += " &END\n"; + res += " &END PERIODIC\n"; + + // Print indicies of QM atoms + // Loop over counter of QM atom types + for (size_t i = 0; i < num_atoms.size(); i++) + { + if (num_atoms[i] > 0) + { + res += formatString(" &QM_KIND %3s\n", periodic_system[i].c_str()); + res += " MM_INDEX"; + // Loop over all QM atoms indexes + for (size_t j = 0; j < nQm; j++) + { + if (parameters_.atomNumbers_[parameters_.qmIndices_[j]] == static_cast(i)) + { + res += formatString(" %d", static_cast(parameters_.qmIndices_[j] + 1)); + } + } + + res += "\n"; + res += " &END QM_KIND\n"; + } + } + + // Print &LINK groups + // Loop over parameters_.link_ + for (size_t i = 0; i < parameters_.link_.size(); i++) + { + res += " &LINK\n"; + res += formatString(" QM_INDEX %d\n", static_cast(parameters_.link_[i].qm) + 1); + res += formatString(" MM_INDEX %d\n", static_cast(parameters_.link_[i].mm) + 1); + res += " &END LINK\n"; + } + + res += " &END QMMM\n"; + + return res; +} + +std::string QMMMInputGenerator::generateMMSection() +{ + std::string res; + + res += " &MM\n"; + res += " &FORCEFIELD\n"; + res += " DO_NONBONDED FALSE\n"; + res += " &END FORCEFIELD\n"; + res += " &POISSON\n"; + res += " &EWALD\n"; + res += " EWALD_TYPE NONE\n"; + res += " &END EWALD\n"; + res += " &END POISSON\n"; + res += " &END MM\n"; + + return res; +} + +std::string QMMMInputGenerator::generateSubsysSection() const +{ + std::string res; + + // Init some numbers + size_t nQm = parameters_.qmIndices_.size(); + size_t nMm = parameters_.mmIndices_.size(); + size_t nAt = nQm + nMm; + + // Count the numbers of individual QM atoms per type + std::vector num_atoms(periodic_system.size(), 0); + for (size_t i = 0; i < nQm; i++) + { + num_atoms[parameters_.atomNumbers_[parameters_.qmIndices_[i]]]++; + } + + res += " &SUBSYS\n"; + + // Print cell parameters + res += " &CELL\n"; + res += formatString(" A %.3lf %.3lf %.3lf\n", box_[0][0] * 10, box_[0][1] * 10, box_[0][2] * 10); + res += formatString(" B %.3lf %.3lf %.3lf\n", box_[1][0] * 10, box_[1][1] * 10, box_[1][2] * 10); + res += formatString(" C %.3lf %.3lf %.3lf\n", box_[2][0] * 10, box_[2][1] * 10, box_[2][2] * 10); + res += " PERIODIC XYZ\n"; + res += " &END CELL\n"; + + // Print topology section + res += " &TOPOLOGY\n"; + + // Placeholder for PDB file name that will be replaced with actuall name during mdrun + res += " COORD_FILE_NAME %s\n"; + + res += " COORD_FILE_FORMAT PDB\n"; + res += " CHARGE_EXTENDED TRUE\n"; + res += " CONNECTIVITY OFF\n"; + res += " &GENERATE\n"; + res += " &ISOLATED_ATOMS\n"; + res += formatString(" LIST %d..%d\n", 1, static_cast(nAt)); + res += " &END\n"; + res += " &END GENERATE\n"; + res += " &END TOPOLOGY\n"; + + // Now we will print basises for all types of QM atoms + // Loop over counter of QM atom types + for (size_t i = 0; i < num_atoms.size(); i++) + { + if (num_atoms[i] > 0) + { + res += " &KIND " + periodic_system[i] + "\n"; + res += " ELEMENT " + periodic_system[i] + "\n"; + res += " BASIS_SET DZVP-MOLOPT-GTH\n"; + res += " POTENTIAL GTH-" + std::string(c_qmmmQMMethodNames[parameters_.qmMethod_]); + res += "\n"; + res += " &END KIND\n"; + } + } + + // Add element kind X - they are represents virtual sites + res += " &KIND X\n"; + res += " ELEMENT H\n"; + res += " &END KIND\n"; + res += " &END SUBSYS\n"; + + return res; +} + + +std::string QMMMInputGenerator::generateCP2KInput() const +{ + + std::string inp; + + // Begin CP2K input generation + inp += generateGlobalSection(); + + inp += "&FORCE_EVAL\n"; + inp += " METHOD QMMM\n"; + + // Add DFT parameters + inp += generateDFTSection(); + + // Add QMMM parameters + inp += generateQMMMSection(); + + // Add MM parameters + inp += generateMMSection(); + + // Add SUBSYS parameters + inp += generateSubsysSection(); + + inp += "&END FORCE_EVAL\n"; + + return inp; +} + +std::string QMMMInputGenerator::generateCP2KPdb() const +{ + + std::string pdb; + + /* Generate *.pdb formatted lines + * and append to std::string pdb + */ + for (size_t i = 0; i < x_.size(); i++) + { + pdb += "ATOM "; + + // Here we need to print i % 100000 because atom counter in *.pdb has only 5 digits + pdb += formatString("%5d ", static_cast(i % 100000)); + + // Atom name + pdb += formatString(" %3s ", periodic_system[parameters_.atomNumbers_[i]].c_str()); + + // Lable atom as QM or MM residue + if (isQMAtom(i)) + { + pdb += " QM 1 "; + } + else + { + pdb += " MM 2 "; + } + + // Coordinates + pdb += formatString("%7.3lf %7.3lf %7.3lf 1.00 0.00 ", + (x_[i][XX] + qmTrans_[XX]) * 10, + (x_[i][YY] + qmTrans_[YY]) * 10, + (x_[i][ZZ] + qmTrans_[ZZ]) * 10); + + // Atom symbol + pdb += formatString(" %3s ", periodic_system[parameters_.atomNumbers_[i]].c_str()); + + // Point charge for MM atoms or 0 for QM atoms + pdb += formatString("%lf\n", q_[i]); + } + + return pdb; +} + +const RVec& QMMMInputGenerator::qmTrans() const +{ + return qmTrans_; +} + +const matrix& QMMMInputGenerator::qmBox() const +{ + return qmBox_; +} + + +RVec computeQMBoxVec(const RVec& a, const RVec& b, const RVec& c, real h, real minNorm, real maxNorm) +{ + RVec vec0(a); + RVec vec1(b); + RVec vec2(c); + RVec dx(0.0, 0.0, 0.0); + RVec res(0.0, 0.0, 0.0); + + // Compute scales sc that will transform a + dx = vec1.cross(vec2); + dx /= dx.norm(); + + // Transform a + res = h / static_cast(fabs(vec0.dot(dx))) * a; + + // If vector is smaller than minL then scale it up + if (res.norm() < minNorm) + { + res *= minNorm / res.norm(); + } + + // If vector is longer than maxL then scale it down + if (res.norm() > maxNorm) + { + res *= maxNorm / res.norm(); + } + + return res; +} + +} // namespace gmx diff --git a/src/gromacs/applied_forces/qmmm/qmmminputgenerator.h b/src/gromacs/applied_forces/qmmm/qmmminputgenerator.h new file mode 100644 index 0000000000..a0c9c91d56 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmminputgenerator.h @@ -0,0 +1,160 @@ +/* + * 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 input generator class for CP2K QMMM + * + * \author Dmitry Morozov + * \ingroup module_applied_forces + */ +#ifndef GMX_APPLIED_FORCES_QMMMINPUTGENERATOR_H +#define GMX_APPLIED_FORCES_QMMMINPUTGENERATOR_H + +#include + +#include "gromacs/pbcutil/pbc.h" +#include "gromacs/utility/arrayref.h" + +#include "qmmmtypes.h" + +namespace gmx +{ + +/*! \internal \brief + * Class that takes QMMMParameters, Coordinates, Point charges, Box dimensions, pbcType. + * Generates QM/MM sample input parameters and pdb-style coordinates for CP2K. + * Input are generated as std::string objects which can be stored in tpr KVT + * and/or flushed into the files. + */ +class QMMMInputGenerator +{ +public: + /*! \brief Construct QMMMInputGenerator from its parameters + * + * \param[in] parameters Structure with QMMM parameters + * \param[in] pbcType Periodic boundary conditions + * \param[in] box Matrix with full box of the system + * \param[in] q Point charges of each atom in the system + * \param[in] x Coordinates of each atom in the system + */ + QMMMInputGenerator(const QMMMParameters& parameters, + PbcType pbcType, + const matrix box, + ArrayRef q, + ArrayRef x); + + /*! \brief Generates sample CP2K input file + * + */ + std::string generateCP2KInput() const; + + /*! \brief Generates PDB file suitable for usage with CP2K. + * In that PDB file Point Charges of MM atoms are provided with Extended Beta field + */ + std::string generateCP2KPdb() const; + + //! \brief Returns computed QM box dimensions + const matrix& qmBox() const; + + //! \brief Returns computed translation vector in order to center QM atoms inside QM box + const RVec& qmTrans() const; + +private: + //! \brief Check if atom belongs to the global index of qmAtoms_ + bool isQMAtom(index globalAtomIndex) const; + + /*!\brief Calculates dimensions and center of the QM box. + * Also evaluates translation for the system in order to center QM atoms inside QM box + * + * \param[in] scale Factor of how much QM box would be bigger than the radius of QM system + * \param[in] minNorm Minimum norm of the QM box vector + */ + void computeQMBox(real scale, real minNorm); + + //! \brief Generates &GLOBAL section of CP2K Input + static std::string generateGlobalSection(); + + //! \brief Generates &DFT section of CP2K Input + std::string generateDFTSection() const; + + //! \brief Generates &QMMM section of CP2K Input + std::string generateQMMMSection() const; + + //! \brief Generates &MM section of CP2K Input + static std::string generateMMSection(); + + //! \brief Generates &SUBSYS section of CP2K Input + std::string generateSubsysSection() const; + + //! QMMM Parameters structure + const QMMMParameters& parameters_; + //! Simulation PbcType + PbcType pbc_; + //! Simulation Box + matrix box_; + //! QM box + matrix qmBox_; + //! PBC-aware center of QM subsystem + RVec qmCenter_; + //! Translation that shifts qmCenter_ to the center of qmBox_ + RVec qmTrans_; + //! Set containing indexes of all QM atoms + std::set qmAtoms_; + //! Atoms point charges + ArrayRef q_; + //! Atoms coordinates + ArrayRef x_; + + //! Scale of the generated QM box with respect to the QM-subsystem size + static constexpr real sc_qmBoxScale = 1.5; + //! Minimum length of the generated QM box vectors in nm + static constexpr real sc_qmBoxMinLength = 1.0; +}; + +/*! \brief Transforms vector a such as distance from it to the plane defined by vectors b and c + * will be h minimum length will be milL and maximum length maxL + * + * \param[in] a Vector which should be scaled + * \param[in] b First vector that forms the plane + * \param[in] c Second vector that forms the plane + * \param[in] h Distance from the end of a to the plane of (b,c) + * \param[in] minNorm Minimum norm of vector + * \param[in] maxNorm Maximum norm of vector + */ +RVec computeQMBoxVec(const RVec& a, const RVec& b, const RVec& c, real h, real minNorm, real maxNorm); + +} // namespace gmx + +#endif // GMX_APPLIED_FORCES_QMMMINPUTGENERATOR_H diff --git a/src/gromacs/applied_forces/qmmm/qmmmtypes.h b/src/gromacs/applied_forces/qmmm/qmmmtypes.h new file mode 100644 index 0000000000..beaf8b1441 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmmtypes.h @@ -0,0 +1,140 @@ +/* + * 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 structers and types needed to evaluate forces and energies for QM/MM + * + * \author Dmitry Morozov + * \author Christian Blau + * \ingroup module_applied_forces + */ +#ifndef GMX_APPLIED_FORCES_QMMMTYPES_H +#define GMX_APPLIED_FORCES_QMMMTYPES_H + +#include +#include +#include + +#include "gromacs/math/vectypes.h" +#include "gromacs/utility/classhelpers.h" +#include "gromacs/utility/enumerationhelpers.h" +#include "gromacs/utility/real.h" + +namespace gmx +{ + +/*! \internal + * \brief Helper structure with indexes of broken bonds between QM and MM + * Used to determine and store pair of QM and MM atoms between which chemical bond is broken + */ +struct LinkFrontier +{ + //! Global index of QM atom at Frontier + index qm; + //! Global index of MM atom at Frontier + index mm; +}; + +/*! \brief Enumerator for supported QM methods + * Also could be INPUT which means external input file provided + * with the name determined by QMMMParameters::qminputfilename_ + */ +enum class QMMMQMMethod +{ + PBE, //!< DFT with PBE functional + BLYP, //!< DFT with BLYP functional + INPUT, //!< User provides suitable input file for QM package + Count +}; + +//! The names of the supported QM methods +static const EnumerationArray c_qmmmQMMethodNames = { + { "PBE", "BLYP", "INPUT" } +}; + +//! 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 " +}; + +/*! \internal + * \brief Holding all parameters needed for QM/MM simulation. + * Also used for setting all default parameter values. + */ +struct QMMMParameters +{ + //! Indicate if QM/MM is active (default false) + bool active_ = false; + //! Indices of the atoms that are part of the QM region (default whole System) + std::vector qmIndices_; + //! Indices of the atoms that are part of the MM region (default no MM atoms) + std::vector mmIndices_; + //! Vector with pairs of indicies defining broken bonds in QMMM (default determined from topology) + std::vector link_; + //! Vector with atomic numbers of all atoms in the system (default determined from topology) + std::vector atomNumbers_; + //! Total charge of QM system (default 0) + int qmCharge_ = 0; + //! Total multiplicity of QM system (default 1) + 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"; + //! 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 + RVec qmTrans_; + + //! Constructor with default initializers for arrays + QMMMParameters() : + qmBox_{ { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } }, qmTrans_{ 0.0, 0.0, 0.0 } + { + } + + GMX_DISALLOW_COPY_AND_ASSIGN(QMMMParameters); +}; + +} // namespace gmx + +#endif // GMX_APPLIED_FORCES_QMMMTYPES_H diff --git a/src/gromacs/applied_forces/qmmm/tests/.clang-tidy b/src/gromacs/applied_forces/qmmm/tests/.clang-tidy new file mode 100644 index 0000000000..458c9c0eed --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/.clang-tidy @@ -0,0 +1,9 @@ +# List of rationales for check suppressions (where known). +# This have to precede the list because inline comments are not +# supported by clang-tidy. +# +# -cppcoreguidelines-avoid-non-const-global-variables +# There are quite a lot of static variables in the test code that +# can not be replaced. +Checks: -cppcoreguidelines-avoid-non-const-global-variables, +InheritParentConfig: true diff --git a/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt b/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt new file mode 100644 index 0000000000..93db354c8b --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt @@ -0,0 +1,38 @@ +# +# 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. + +gmx_add_unit_test(QMMMAppliedForcesUnitTest qmmm_applied_forces-test + CPP_SOURCE_FILES + qmmminputgenerator.cpp + ) diff --git a/src/gromacs/applied_forces/qmmm/tests/qmmminputgenerator.cpp b/src/gromacs/applied_forces/qmmm/tests/qmmminputgenerator.cpp new file mode 100644 index 0000000000..ae2f8c310b --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/qmmminputgenerator.cpp @@ -0,0 +1,147 @@ +/* + * 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 QMMMInputGenerator class for QMMM MDModule + * + * \author Dmitry Morozov + * \ingroup module_applied_forces + */ +#include "gmxpre.h" + +#include "gromacs/applied_forces/qmmm/qmmminputgenerator.h" + +#include + +#include + +#include "gromacs/math/vec.h" +#include "gromacs/utility/arrayref.h" +#include "testutils/refdata.h" +#include "testutils/testasserts.h" +#include "testutils/testfilemanager.h" + +namespace gmx +{ + +class QMMMInputGeneratorTest : public ::testing::Test +{ +public: + //! Reference input 2x TIP3P waters (first one QM) and no Link atoms + void initParameters2WatersNoLink() + { + parameters_.qmIndices_ = { 0, 1, 2 }; + parameters_.mmIndices_ = { 3, 4, 5 }; + parameters_.atomNumbers_ = { 8, 1, 1, 8, 1, 1 }; + q_ = { 0.0, 0.0, 0.0, -0.834, 0.417, 0.417 }; + pbc_ = PbcType::Xyz; + matrix box1_ = { { 2.0, 0.0, 0.0 }, { 0.0, 2.0, 0.0 }, { 0.0, 0.0, 2.0 } }; + copy_mat(box1_, box_); + coords_ = { RVec(0.005, 0.600, 0.244), RVec(-0.017, 0.690, 0.270), + RVec(0.051, 0.610, 0.161), RVec(0.975, 0.631, 1.569), + RVec(0.951, 0.615, 1.661), RVec(0.916, 0.701, 1.541) }; + } + + //! Reference input 2x TIP3P waters (first two atoms are QM) and one link should be generated + void initParameters2WatersWithLink() + { + parameters_.qmIndices_ = { 0, 1 }; + parameters_.mmIndices_ = { 2, 3, 4, 5 }; + parameters_.atomNumbers_ = { 8, 1, 1, 8, 1, 1 }; + parameters_.link_ = { { 0, 2 } }; + q_ = { 0.0, 0.0, 0.417, -0.834, 0.417, 0.417 }; + pbc_ = PbcType::Xyz; + matrix box1_ = { { 2.0, 0.0, 0.0 }, { 0.0, 2.0, 0.0 }, { 0.0, 0.0, 2.0 } }; + copy_mat(box1_, box_); + coords_ = { RVec(0.005, 0.600, 0.244), RVec(-0.017, 0.690, 0.270), + RVec(0.051, 0.610, 0.161), RVec(0.975, 0.631, 1.569), + RVec(0.951, 0.615, 1.661), RVec(0.916, 0.701, 1.541) }; + } + +protected: + QMMMParameters parameters_; + PbcType pbc_; + matrix box_ = { { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 }, { 0.0, 0.0, 0.0 } }; + + std::vector q_; + std::vector coords_; +}; + + +TEST_F(QMMMInputGeneratorTest, CanConstruct) +{ + initParameters2WatersNoLink(); + EXPECT_NO_THROW(QMMMInputGenerator inpGen(parameters_, pbc_, box_, q_, coords_)); +} + +TEST_F(QMMMInputGeneratorTest, TwoWatersPBENoLink) +{ + // Create Input Generator + initParameters2WatersNoLink(); + QMMMInputGenerator inpGen(parameters_, pbc_, box_, q_, coords_); + + gmx::test::TestReferenceData data; + gmx::test::TestReferenceChecker checker(data.rootChecker()); + + // Tolerance of all coordinates and vectors should be 1E-3 (as in gro or pdb files) + checker.setDefaultTolerance(gmx::test::absoluteTolerance(0.001)); + checker.checkVector(inpGen.qmTrans(), "QM Translation"); + checker.checkVector(inpGen.qmBox()[0], "QM Box Vector 1"); + checker.checkVector(inpGen.qmBox()[1], "QM Box Vector 2"); + checker.checkVector(inpGen.qmBox()[2], "QM Box Vector 3"); + checker.checkString(inpGen.generateCP2KInput(), "Input"); + checker.checkString(inpGen.generateCP2KPdb(), "PDB"); +} + +TEST_F(QMMMInputGeneratorTest, TwoWatersPBEWithLink) +{ + // Create Input Generator + initParameters2WatersWithLink(); + QMMMInputGenerator inpGen(parameters_, pbc_, box_, q_, coords_); + + gmx::test::TestReferenceData data; + gmx::test::TestReferenceChecker checker(data.rootChecker()); + + // Tolerance of all coordinates and vectors should be 1E-3 (as in gro or pdb files) + checker.setDefaultTolerance(gmx::test::absoluteTolerance(0.001)); + checker.checkVector(inpGen.qmTrans(), "QM Translation"); + checker.checkVector(inpGen.qmBox()[0], "QM Box Vector 1"); + checker.checkVector(inpGen.qmBox()[1], "QM Box Vector 2"); + checker.checkVector(inpGen.qmBox()[2], "QM Box Vector 3"); + checker.checkString(inpGen.generateCP2KInput(), "Input"); + checker.checkString(inpGen.generateCP2KPdb(), "PDB"); +} + +} // namespace gmx diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBENoLink.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBENoLink.xml new file mode 100644 index 0000000000..d922f8a3c4 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBENoLink.xml @@ -0,0 +1,149 @@ + + + + + 0.48699999 + -0.13333333 + 0.27499998 + + + 1 + 0 + 0 + + + 0 + 1 + 0 + + + 0 + 0 + 1 + + &GLOBAL + PRINT_LEVEL LOW + PROJECT GROMACS + RUN_TYPE ENERGY_FORCE +&END GLOBAL +&FORCE_EVAL + METHOD QMMM + &DFT + CHARGE 0 + MULTIPLICITY 1 + BASIS_SET_FILE_NAME BASIS_MOLOPT + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + NGRIDS 5 + CUTOFF 450 + REL_CUTOFF 50 + COMMENSURATE + &END MGRID + &SCF + SCF_GUESS RESTART + EPS_SCF 5.0E-8 + MAX_SCF 20 + &OT T + MINIMIZER DIIS + STEPSIZE 0.15 + PRECONDITIONER FULL_ALL + &END OT + &OUTER_SCF T + MAX_SCF 20 + EPS_SCF 5.0E-8 + &END OUTER_SCF + &END SCF + &XC + DENSITY_CUTOFF 1.0E-12 + GRADIENT_CUTOFF 1.0E-12 + TAU_CUTOFF 1.0E-12 + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &QS + METHOD GPW + EPS_DEFAULT 1.0E-10 + EXTRAPOLATION ASPC + EXTRAPOLATION_ORDER 4 + &END QS + &END DFT + &QMMM + &CELL + A 10.000 0.000 0.000 + B 0.000 10.000 0.000 + C 0.000 0.000 10.000 + PERIODIC XYZ + &END CELL + CENTER EVERY_STEP + CENTER_GRID TRUE + &WALLS + TYPE REFLECTIVE + &END WALLS + ECOUPL GAUSS + USE_GEEP_LIB 12 + &PERIODIC + GMAX 1.0E+00 + &MULTIPOLE ON + RCUT 1.0E+01 + EWALD_PRECISION 1.0E-06 + &END + &END PERIODIC + &QM_KIND H + MM_INDEX 2 3 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &END QMMM + &MM + &FORCEFIELD + DO_NONBONDED FALSE + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE NONE + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + A 20.000 0.000 0.000 + B 0.000 20.000 0.000 + C 0.000 0.000 20.000 + PERIODIC XYZ + &END CELL + &TOPOLOGY + COORD_FILE_NAME %s + COORD_FILE_FORMAT PDB + CHARGE_EXTENDED TRUE + CONNECTIVITY OFF + &GENERATE + &ISOLATED_ATOMS + LIST 1..6 + &END + &END GENERATE + &END TOPOLOGY + &KIND H + ELEMENT H + BASIS_SET DZVP-MOLOPT-GTH + POTENTIAL GTH-PBE + &END KIND + &KIND O + ELEMENT O + BASIS_SET DZVP-MOLOPT-GTH + POTENTIAL GTH-PBE + &END KIND + &KIND X + ELEMENT H + &END KIND + &END SUBSYS +&END FORCE_EVAL + + ATOM 0 O QM 1 4.920 4.667 5.190 1.00 0.00 O 0.000000 +ATOM 1 H QM 1 4.700 5.567 5.450 1.00 0.00 H 0.000000 +ATOM 2 H QM 1 5.380 4.767 4.360 1.00 0.00 H 0.000000 +ATOM 3 O MM 2 14.620 4.977 18.440 1.00 0.00 O -0.834000 +ATOM 4 H MM 2 14.380 4.817 19.360 1.00 0.00 H 0.417000 +ATOM 5 H MM 2 14.030 5.677 18.160 1.00 0.00 H 0.417000 + + diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBEWithLink.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBEWithLink.xml new file mode 100644 index 0000000000..282019cb39 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMInputGeneratorTest_TwoWatersPBEWithLink.xml @@ -0,0 +1,153 @@ + + + + + 0.50599998 + -0.14499998 + 0.243 + + + 1 + 0 + 0 + + + 0 + 1 + 0 + + + 0 + 0 + 1 + + &GLOBAL + PRINT_LEVEL LOW + PROJECT GROMACS + RUN_TYPE ENERGY_FORCE +&END GLOBAL +&FORCE_EVAL + METHOD QMMM + &DFT + CHARGE 0 + MULTIPLICITY 1 + BASIS_SET_FILE_NAME BASIS_MOLOPT + POTENTIAL_FILE_NAME POTENTIAL + &MGRID + NGRIDS 5 + CUTOFF 450 + REL_CUTOFF 50 + COMMENSURATE + &END MGRID + &SCF + SCF_GUESS RESTART + EPS_SCF 5.0E-8 + MAX_SCF 20 + &OT T + MINIMIZER DIIS + STEPSIZE 0.15 + PRECONDITIONER FULL_ALL + &END OT + &OUTER_SCF T + MAX_SCF 20 + EPS_SCF 5.0E-8 + &END OUTER_SCF + &END SCF + &XC + DENSITY_CUTOFF 1.0E-12 + GRADIENT_CUTOFF 1.0E-12 + TAU_CUTOFF 1.0E-12 + &XC_FUNCTIONAL PBE + &END XC_FUNCTIONAL + &END XC + &QS + METHOD GPW + EPS_DEFAULT 1.0E-10 + EXTRAPOLATION ASPC + EXTRAPOLATION_ORDER 4 + &END QS + &END DFT + &QMMM + &CELL + A 10.000 0.000 0.000 + B 0.000 10.000 0.000 + C 0.000 0.000 10.000 + PERIODIC XYZ + &END CELL + CENTER EVERY_STEP + CENTER_GRID TRUE + &WALLS + TYPE REFLECTIVE + &END WALLS + ECOUPL GAUSS + USE_GEEP_LIB 12 + &PERIODIC + GMAX 1.0E+00 + &MULTIPOLE ON + RCUT 1.0E+01 + EWALD_PRECISION 1.0E-06 + &END + &END PERIODIC + &QM_KIND H + MM_INDEX 2 + &END QM_KIND + &QM_KIND O + MM_INDEX 1 + &END QM_KIND + &LINK + QM_INDEX 1 + MM_INDEX 3 + &END LINK + &END QMMM + &MM + &FORCEFIELD + DO_NONBONDED FALSE + &END FORCEFIELD + &POISSON + &EWALD + EWALD_TYPE NONE + &END EWALD + &END POISSON + &END MM + &SUBSYS + &CELL + A 20.000 0.000 0.000 + B 0.000 20.000 0.000 + C 0.000 0.000 20.000 + PERIODIC XYZ + &END CELL + &TOPOLOGY + COORD_FILE_NAME %s + COORD_FILE_FORMAT PDB + CHARGE_EXTENDED TRUE + CONNECTIVITY OFF + &GENERATE + &ISOLATED_ATOMS + LIST 1..6 + &END + &END GENERATE + &END TOPOLOGY + &KIND H + ELEMENT H + BASIS_SET DZVP-MOLOPT-GTH + POTENTIAL GTH-PBE + &END KIND + &KIND O + ELEMENT O + BASIS_SET DZVP-MOLOPT-GTH + POTENTIAL GTH-PBE + &END KIND + &KIND X + ELEMENT H + &END KIND + &END SUBSYS +&END FORCE_EVAL + + ATOM 0 O QM 1 5.110 4.550 4.870 1.00 0.00 O 0.000000 +ATOM 1 H QM 1 4.890 5.450 5.130 1.00 0.00 H 0.000000 +ATOM 2 H MM 2 5.570 4.650 4.040 1.00 0.00 H 0.417000 +ATOM 3 O MM 2 14.810 4.860 18.120 1.00 0.00 O -0.834000 +ATOM 4 H MM 2 14.570 4.700 19.040 1.00 0.00 H 0.417000 +ATOM 5 H MM 2 14.220 5.560 17.840 1.00 0.00 H 0.417000 + + -- 2.22.0