#
# 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.
add_subdirectory(awh)
add_subdirectory(densityfitting)
+add_subdirectory(qmmm)
if (BUILD_TESTING)
add_subdirectory(tests)
--- /dev/null
+#
+# 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()
--- /dev/null
+/*
+ * 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 <dmitry.morozov@jyu.fi>
+ * \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<const real> q,
+ ArrayRef<const RVec> 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<int> 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<index>(i))
+ {
+ res += formatString(" %d", static_cast<int>(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<int>(parameters_.link_[i].qm) + 1);
+ res += formatString(" MM_INDEX %d\n", static_cast<int>(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<int> 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<int>(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<int>(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<real>(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
--- /dev/null
+/*
+ * 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 <dmitry.morozov@jyu.fi>
+ * \ingroup module_applied_forces
+ */
+#ifndef GMX_APPLIED_FORCES_QMMMINPUTGENERATOR_H
+#define GMX_APPLIED_FORCES_QMMMINPUTGENERATOR_H
+
+#include <set>
+
+#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<const real> q,
+ ArrayRef<const RVec> 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<index> qmAtoms_;
+ //! Atoms point charges
+ ArrayRef<const real> q_;
+ //! Atoms coordinates
+ ArrayRef<const RVec> 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
--- /dev/null
+/*
+ * 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 <dmitry.morozov@jyu.fi>
+ * \author Christian Blau <blau@kth.se>
+ * \ingroup module_applied_forces
+ */
+#ifndef GMX_APPLIED_FORCES_QMMMTYPES_H
+#define GMX_APPLIED_FORCES_QMMMTYPES_H
+
+#include <memory>
+#include <string>
+#include <vector>
+
+#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<QMMMQMMethod, const char*> c_qmmmQMMethodNames = {
+ { "PBE", "BLYP", "INPUT" }
+};
+
+//! 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 "
+};
+
+/*! \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<index> qmIndices_;
+ //! Indices of the atoms that are part of the MM region (default no MM atoms)
+ std::vector<index> mmIndices_;
+ //! Vector with pairs of indicies defining broken bonds in QMMM (default determined from topology)
+ std::vector<LinkFrontier> link_;
+ //! Vector with atomic numbers of all atoms in the system (default determined from topology)
+ std::vector<int> 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
--- /dev/null
+# 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
--- /dev/null
+#
+# 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
+ )
--- /dev/null
+/*
+ * 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 <dmitry.morozov@jyu.fi>
+ * \ingroup module_applied_forces
+ */
+#include "gmxpre.h"
+
+#include "gromacs/applied_forces/qmmm/qmmminputgenerator.h"
+
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#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<real> q_;
+ std::vector<RVec> 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
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Vector Name="QM Translation">
+ <Real Name="X">0.48699999</Real>
+ <Real Name="Y">-0.13333333</Real>
+ <Real Name="Z">0.27499998</Real>
+ </Vector>
+ <Vector Name="QM Box Vector 1">
+ <Real Name="X">1</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector Name="QM Box Vector 2">
+ <Real Name="X">0</Real>
+ <Real Name="Y">1</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector Name="QM Box Vector 3">
+ <Real Name="X">0</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">1</Real>
+ </Vector>
+ <String Name="Input">&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
+</String>
+ <String Name="PDB">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
+</String>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Vector Name="QM Translation">
+ <Real Name="X">0.50599998</Real>
+ <Real Name="Y">-0.14499998</Real>
+ <Real Name="Z">0.243</Real>
+ </Vector>
+ <Vector Name="QM Box Vector 1">
+ <Real Name="X">1</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector Name="QM Box Vector 2">
+ <Real Name="X">0</Real>
+ <Real Name="Y">1</Real>
+ <Real Name="Z">0</Real>
+ </Vector>
+ <Vector Name="QM Box Vector 3">
+ <Real Name="X">0</Real>
+ <Real Name="Y">0</Real>
+ <Real Name="Z">1</Real>
+ </Vector>
+ <String Name="Input">&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
+</String>
+ <String Name="PDB">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
+</String>
+</ReferenceData>
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff