From 250d67b01d26e602169d8415daa121abd729b56b Mon Sep 17 00:00:00 2001 From: Dmitry Morozov Date: Tue, 14 Sep 2021 07:57:27 +0000 Subject: [PATCH] Added class for making QMMM-related topology modifications --- .../applied_forces/qmmm/CMakeLists.txt | 1 + .../qmmm/qmmmtopologypreprocessor.cpp | 652 ++++++++++++++++++ .../qmmm/qmmmtopologypreprocessor.h | 206 ++++++ .../applied_forces/qmmm/tests/CMakeLists.txt | 1 + .../qmmm/tests/qmmmtopologypreprocessor.cpp | 244 +++++++ ...AlanineDipeptideWithLinksNoConstraints.xml | 18 + ...anineDipeptideWithLinksWithConstraints.xml | 18 + ...eprocessorTest_FourWatersFirstQMNoLink.xml | 18 + ...rocessorTest_FourWatersFirstQMWithLink.xml | 18 + ...orTest_FourWatersSeondAndForthQMNoLink.xml | 18 + ...ologyPreprocessorTest_RemovingQMVsites.xml | 18 + 11 files changed, 1212 insertions(+) create mode 100644 src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.cpp create mode 100644 src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.h create mode 100644 src/gromacs/applied_forces/qmmm/tests/qmmmtopologypreprocessor.cpp create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksNoConstraints.xml create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksWithConstraints.xml create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMNoLink.xml create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMWithLink.xml create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersSeondAndForthQMNoLink.xml create mode 100644 src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_RemovingQMVsites.xml diff --git a/src/gromacs/applied_forces/qmmm/CMakeLists.txt b/src/gromacs/applied_forces/qmmm/CMakeLists.txt index c64cadff3f..ba20be3824 100644 --- a/src/gromacs/applied_forces/qmmm/CMakeLists.txt +++ b/src/gromacs/applied_forces/qmmm/CMakeLists.txt @@ -34,6 +34,7 @@ gmx_add_libgromacs_sources( qmmminputgenerator.cpp + qmmmtopologypreprocessor.cpp ) if (BUILD_TESTING) diff --git a/src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.cpp b/src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.cpp new file mode 100644 index 0000000000..6f4b7be73a --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.cpp @@ -0,0 +1,652 @@ +/* + * 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 QMMMTopologyPreprocessor + * + * \author Dmitry Morozov + * \ingroup module_applied_forces + */ +#include "gmxpre.h" + +#include "qmmmtopologypreprocessor.h" + +#include "gromacs/selection/indexutil.h" +#include "gromacs/topology/mtop_lookup.h" +#include "gromacs/topology/mtop_util.h" +#include "gromacs/topology/topology.h" +#include "gromacs/utility/exceptions.h" +#include "gromacs/utility/fatalerror.h" + +namespace gmx +{ + +QMMMTopologyPreprocessor::QMMMTopologyPreprocessor(ArrayRef qmIndices) : + qmIndices_(qmIndices.begin(), qmIndices.end()) +{ +} + +void QMMMTopologyPreprocessor::preprocess(gmx_mtop_t* mtop) +{ + // 1) Split QM-containing molecules from other molecules in blocks + splitQMblocks(mtop); + + // 2) Nullify charges on all virtual sites consisting of QM only atoms + modifyQMMMVirtualSites(mtop); + + // 3) Nullify charges on all QM atoms + removeQMClassicalCharges(mtop); + + // 4) Exclude LJ interactions between QM atoms + addQMLJExclusions(mtop); + + // 5) Build atomNumbers vector with atomic numbers of all atoms + buildQMMMAtomNumbers(mtop); + + // 6) Make F_CONNBOND between atoms within QM region + modifyQMMMTwoCenterInteractions(mtop); + + // 7) Remove angles and settles containing 2 or more QM atoms + modifyQMMMThreeCenterInteractions(mtop); + + // 8) Remove dihedrals containing 3 or more QM atoms + modifyQMMMFourCenterInteractions(mtop); + + // 9) Build vector containing pairs of bonded QM - MM atoms (Link frontier) + buildQMMMLink(mtop); + + // finalize topology + mtop->finalize(); +} + +const QMMMTopologyInfo& QMMMTopologyPreprocessor::topInfo() const +{ + return topInfo_; +} + +ArrayRef QMMMTopologyPreprocessor::atomNumbers() const +{ + return atomNumbers_; +} + +ArrayRef QMMMTopologyPreprocessor::atomCharges() const +{ + return atomCharges_; +} + +ArrayRef QMMMTopologyPreprocessor::linkFrontier() const +{ + return linkFrontier_; +} + +bool QMMMTopologyPreprocessor::isQMAtom(index globalAtomIndex) +{ + return (qmIndices_.find(globalAtomIndex) != qmIndices_.end()); +} + +void QMMMTopologyPreprocessor::splitQMblocks(gmx_mtop_t* mtop) +{ + + // Global counter of atoms + index iAt = 0; + + /* Counter of molecules point to the specific moltype + * i.e molblock 0 has 2 molecules have moltype 0 and molblock 2 has 1 additional molecule of type 0 + * then will be numMoleculesOfType[0] = 2 + 1 = 3 + * That counter is needed to decide whether we should create a new moltype for the molecule contatining QM atoms + * or we could modify existing moltype if there is only one molecule of that type + */ + std::vector numMoleculesOfType(mtop->moltype.size()); + for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++) + { + numMoleculesOfType[mtop->molblock[molBlockIndex].type] += mtop->molblock[molBlockIndex].nmol; + } + + // Loop over all blocks in topology + // molBlockIndex - current index of block in mtop + for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++) + { + // Initialize block as non-QM first + bQMBlock_.push_back(false); + + // Pointer to current block + gmx_molblock_t* molBlock = &mtop->molblock[molBlockIndex]; + + // Number of atoms in all molecules of current block + const int numAtomsInMolecule = mtop->moltype[molBlock->type].atoms.nr; + + // Loop over all molecules in molBlock + // mol - index of current molecule in molBlock + for (int mol = 0; mol < molBlock->nmol; mol++) + { + + // search for QM atoms in current molecule + bool bQMMM = false; + for (int i = 0; i < numAtomsInMolecule; i++) + { + if (isQMAtom(iAt)) + { + bQMMM = true; + } + iAt++; + } + + // Apparently current molecule (molBlockIndex, mol) contains QM atoms + // We should split it from the current block and create new blocks + // For that molecule and all molecules after it new block will be created + if (bQMMM) + { + // if this block contains only 1 molecule, then splitting not needed + if (molBlock->nmol > 1) + { + // If current molecule is not the first in the block, + // then we need to first create block for it and all molecule before it + if (mol > 0) + { + // Split the molblock at this molecule + auto pos = mtop->molblock.begin() + molBlockIndex + 1; + mtop->molblock.insert(pos, mtop->molblock[molBlockIndex]); + mtop->molblock[molBlockIndex].nmol = mol; + mtop->molblock[molBlockIndex + 1].nmol -= mol; + bQMBlock_[molBlockIndex] = false; + bQMBlock_.push_back(true); + molBlockIndex++; + molBlock = &mtop->molblock[molBlockIndex]; + } + + // If current molecule is not the only one in new block, + // Then we split new block after that molecule + if (molBlock->nmol > 1) + { + auto pos = mtop->molblock.begin() + molBlockIndex + 1; + mtop->molblock.insert(pos, mtop->molblock[molBlockIndex]); + molBlock = &mtop->molblock[molBlockIndex]; + mtop->molblock[molBlockIndex].nmol = 1; + mtop->molblock[molBlockIndex + 1].nmol -= 1; + bQMBlock_[molBlockIndex] = true; + } + } + else + { + bQMBlock_[molBlockIndex] = true; + } + + // Create a copy of a moltype for a molecule + // containing QM atoms and append it in the end of the moltype vector + // if there is only 1 molecule pointing to that type then skip that step + if (numMoleculesOfType[molBlock->type] > 1) + { + // Here comes a huge piece of "not so good" code, because of deleted operator= from gmx_moltype_t + std::vector temp(mtop->moltype.size()); + for (size_t i = 0; i < mtop->moltype.size(); ++i) + { + copy_moltype(&mtop->moltype[i], &temp[i]); + } + mtop->moltype.resize(temp.size() + 1); + for (size_t i = 0; i < temp.size(); ++i) + { + copy_moltype(&temp[i], &mtop->moltype[i]); + } + copy_moltype(&mtop->moltype[molBlock->type], &mtop->moltype.back()); + + // Set the molecule type for the QMMM molblock + molBlock->type = mtop->moltype.size() - 1; + } + } + } + } + + // Call finalize() to rebuild Block Indicies or else atoms lookup will fail + mtop->finalize(); +} + +void QMMMTopologyPreprocessor::removeQMClassicalCharges(gmx_mtop_t* mtop) +{ + // Loop over all atoms and remove charge if they are QM atoms. + // Sum-up total removed charge and remaning charge on MM atoms + // Build atomCharges_ vector + int molBlockIndex = 0; + for (int i = 0; i < mtop->natoms; i++) + { + int indexInMolecule; + mtopGetMolblockIndex(*mtop, i, &molBlockIndex, nullptr, &indexInMolecule); + t_atom* atom = &mtop->moltype[mtop->molblock[molBlockIndex].type].atoms.atom[indexInMolecule]; + if (isQMAtom(i)) + { + topInfo_.totalClassicalChargeOfQMAtoms += atom->q; + atom->q = 0.0; + atom->qB = 0.0; + } + else + { + topInfo_.remainingMMCharge += atom->q; + } + + atomCharges_.push_back(atom->q); + } +} + +void QMMMTopologyPreprocessor::addQMLJExclusions(gmx_mtop_t* mtop) +{ + // Add all QM atoms to the mtop->intermolecularExclusionGroup + mtop->intermolecularExclusionGroup.reserve(mtop->intermolecularExclusionGroup.size() + + qmIndices_.size()); + for (auto i : qmIndices_) + { + mtop->intermolecularExclusionGroup.push_back(i); + topInfo_.numExclusionsMade++; + } +} + +void QMMMTopologyPreprocessor::buildQMMMAtomNumbers(gmx_mtop_t* mtop) +{ + // Save to atomNumbers_ atom numbers of all atoms + AtomIterator atoms(*mtop); + while (atoms->globalAtomNumber() < mtop->natoms) + { + // Check if we have valid atomnumbers + if (atoms->atom().atomnumber < 0) + { + gmx_fatal(FARGS, + "Atoms %d does not have atomic number needed for QMMM. Check atomtypes " + "section in your topology or forcefield.", + atoms->globalAtomNumber()); + } + + atomNumbers_.push_back(atoms->atom().atomnumber); + atoms++; + } + + // Save in topInfo_ number of QM and MM atoms + topInfo_.numQMAtoms += gmx::ssize(qmIndices_); + topInfo_.numMMAtoms += mtop->natoms - gmx::ssize(qmIndices_); +} + +void QMMMTopologyPreprocessor::modifyQMMMTwoCenterInteractions(gmx_mtop_t* mtop) +{ + // Loop over all blocks in topology + // molBlockIndex - index of current block in mtop + for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++) + { + // check if current block contains QM atoms + if (bQMBlock_[molBlockIndex]) + { + // molType - strucutre with current block type + gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type]; + + // start - global index if the first atom in the current block + int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart; + + // loop over all interaction types + for (int ftype = 0; ftype < F_NRE; ftype++) + { + // If not bonded interaction or F_CONNBONDS, or some form of Restraints, + // or not pair interaction, or no interactions of that type: then go the next type + if (!(interaction_function[ftype].flags & IF_BOND) || ftype == F_CONNBONDS + || ftype == F_RESTRBONDS || ftype == F_HARMONIC || ftype == F_DISRES || ftype == F_ORIRES + || ftype == F_ANGRESZ || NRAL(ftype) != 2 || molType->ilist[ftype].empty()) + { + continue; + } + + // Number of elements in the iatoms array for the current ftype + const int numInteractionElements = NRAL(ftype) + 1; + + // Buffer for preserving interactions that are retained + AtomGroupIndices iatomsBuf; + + // Loop over all interactions of ftype + for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements) + { + + // If both atoms are QM and it is IF_CHEMBOND then convert it to F_CONNBONDS + if (isQMAtom(molType->ilist[ftype].iatoms[j + 1] + start) + && isQMAtom(molType->ilist[ftype].iatoms[j + 2] + start)) + { + + // Add chemical bond to the F_CONNBONDS (bond type 5) + if (IS_CHEMBOND(ftype)) + { + // Bond type is not used in F_CONNBONDS, so for generated bonds we set it to -1 + const int connBondsType = -1; + + // Add new CONNBOND between atoms + molType->ilist[F_CONNBONDS].iatoms.push_back(connBondsType); + molType->ilist[F_CONNBONDS].iatoms.push_back( + molType->ilist[ftype].iatoms[j + 1]); + molType->ilist[F_CONNBONDS].iatoms.push_back( + molType->ilist[ftype].iatoms[j + 2]); + + topInfo_.numConnBondsAdded++; + } + + // Since the current bond is not copied into iatomsBuf it will be removed from the topology + topInfo_.numBondsRemoved++; + } + else + { + // If one of atoms is not QM then copy interaction into iatomsBuf + for (int k = 0; k < numInteractionElements; k++) + { + iatomsBuf.push_back(molType->ilist[ftype].iatoms[k + j]); + } + } + } + + // Swap iatomsBuf and molType->ilist[ftype].iatoms + molType->ilist[ftype].iatoms.swap(iatomsBuf); + } + } + } +} + +void QMMMTopologyPreprocessor::buildQMMMLink(gmx_mtop_t* mtop) +{ + // Loop over all blocks in topology + // molBlockIndex - index of current block in mtop + for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++) + { + // check if current block contains QM atoms + if (bQMBlock_[molBlockIndex]) + { + // molType - strucutre with current block type + gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type]; + + // start - gloabl index of the first atom in the current block + int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart; + + // loop over all interaction types + for (int ftype = 0; ftype < F_NRE; ftype++) + { + // If not chemical bond interaction or not pair interaction + // or no interactions of that type: then skip current ftype + if (!(interaction_function[ftype].flags & IF_CHEMBOND) || NRAL(ftype) != 2 + || molType->ilist[ftype].empty()) + { + continue; + } + + // Number of elements in the iatoms array for the current ftype + const int numInteractionElements = NRAL(ftype) + 1; + + // Loop over all interactions of ftype + for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements) + { + // Global indexes of atoms involved into the interaction + int a1 = molType->ilist[ftype].iatoms[j + 1] + start; + int a2 = molType->ilist[ftype].iatoms[j + 2] + start; + + // Update Link Frontier List if one of the atoms QM and one MM + if (isQMAtom(a1) && !isQMAtom(a2)) + { + linkFrontier_.push_back({ a1, a2 }); + topInfo_.numLinkBonds++; + } + if (isQMAtom(a2) && !isQMAtom(a1)) + { + linkFrontier_.push_back({ a2, a1 }); + topInfo_.numLinkBonds++; + } + + // Check if it is constrained bond within QM subsystem + if (isQMAtom(a2) && isQMAtom(a1) && (interaction_function[ftype].flags & IF_CONSTRAINT)) + { + topInfo_.numConstrainedBondsInQMSubsystem++; + } + } + } + } + } +} + +void QMMMTopologyPreprocessor::modifyQMMMThreeCenterInteractions(gmx_mtop_t* mtop) +{ + // Loop over all blocks in topology + // molBlockIndex - index of current block in mtop + for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++) + { + // check if current block contains QM atoms + if (bQMBlock_[molBlockIndex]) + { + // molType - strucutre with current block type + gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type]; + + // start - global index of the first atom in the current block + int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart; + + // loop over all interaction types + for (int ftype = 0; ftype < F_NRE; ftype++) + { + // If not bonded interaction or Restraints + // or not three-particle interaction or no interactions of that type + // and not Settle: then go the next type + if ((!(interaction_function[ftype].flags & IF_BOND) || ftype == F_RESTRANGLES + || NRAL(ftype) != 3 || molType->ilist[ftype].empty()) + && (ftype != F_SETTLE)) + { + continue; + } + + // Number of elements in the iatoms array for the current ftype + const int numInteractionElements = NRAL(ftype) + 1; + + // Buffer for preserving interactions that are retained + AtomGroupIndices iatomsBuf; + + // Loop over all interactions of ftype + for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements) + { + // Calculate number of qm atoms in the interaction + int numQm = 0; + for (int k = 1; k <= NRAL(ftype); k++) + { + if (isQMAtom(molType->ilist[ftype].iatoms[j + k] + start)) + { + numQm++; + } + } + + // If at least 2 atoms are QM then remove interaction + if (numQm >= 2) + { + // If this is SETTLE then replace it with two F_CONNBONDS + if (ftype == F_SETTLE) + { + // Bond type is not used in F_CONNBONDS, so for generated bonds we set it to -1 + const int connBondsType = -1; + + // Add CONNBOND between atoms 1 and 2 first + molType->ilist[F_CONNBONDS].iatoms.push_back(connBondsType); + molType->ilist[F_CONNBONDS].iatoms.push_back( + molType->ilist[ftype].iatoms[j + 1]); + molType->ilist[F_CONNBONDS].iatoms.push_back( + molType->ilist[ftype].iatoms[j + 2]); + // Then between atoms 1 and 3 + molType->ilist[F_CONNBONDS].iatoms.push_back(connBondsType); + molType->ilist[F_CONNBONDS].iatoms.push_back( + molType->ilist[ftype].iatoms[j + 1]); + molType->ilist[F_CONNBONDS].iatoms.push_back( + molType->ilist[ftype].iatoms[j + 3]); + + topInfo_.numConnBondsAdded += 2; + topInfo_.numSettleRemoved++; + } + else + { + // If it is normal angle then it should be just removed + topInfo_.numAnglesRemoved++; + } + } + else + { + // If several (>1) atoms is not QM then preserve interaction in the iatomsBuf + for (int k = 0; k < numInteractionElements; k++) + { + iatomsBuf.push_back(molType->ilist[ftype].iatoms[k + j]); + } + } + } + + // Swap iatomsBuf and molType->ilist[ftype].iatoms + molType->ilist[ftype].iatoms.swap(iatomsBuf); + } + } + } +} + + +void QMMMTopologyPreprocessor::modifyQMMMFourCenterInteractions(gmx_mtop_t* mtop) +{ + // Loop over all blocks in topology + // molBlockIndex - index of current block in mtop + for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++) + { + // check if current block contains QM atoms + if (bQMBlock_[molBlockIndex]) + { + // molType - strucutre with current block type + gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type]; + + // start - global index of the first atom in the current block + int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart; + + // loop over all interaction types + for (int ftype = 0; ftype < F_NRE; ftype++) + { + // If not bonded interaction or Restraints + // or not four-particle interaction or no interactions of that type: then go the next type + if (!(interaction_function[ftype].flags & IF_BOND) || ftype == F_RESTRDIHS + || NRAL(ftype) != 4 || molType->ilist[ftype].empty()) + { + continue; + } + + // Number of elements in the iatoms array for the current ftype + const int numInteractionElements = NRAL(ftype) + 1; + + // Buffer for preserving interactions that are retained + AtomGroupIndices iatomsBuf; + + // Loop over all interactions of ftype + for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements) + { + // Calculate number of qm atoms in the interaction + int numQm = 0; + for (int k = 1; k <= NRAL(ftype); k++) + { + if (isQMAtom(molType->ilist[ftype].iatoms[j + k] + start)) + { + numQm++; + } + } + + // If at least 3 atoms are QM then remove interaction + if (numQm >= 3) + { + topInfo_.numDihedralsRemoved++; + } + else + { + // If several (>1) atoms is not QM then preserve interaction in the iatomsBuf + for (int k = 0; k < numInteractionElements; k++) + { + iatomsBuf.push_back(molType->ilist[ftype].iatoms[k + j]); + } + } + } + + // Swap iatomsBuf and molType->ilist[ftype].iatoms + molType->ilist[ftype].iatoms.swap(iatomsBuf); + } + } + } +} + +void QMMMTopologyPreprocessor::modifyQMMMVirtualSites(gmx_mtop_t* mtop) +{ + // Loop over all blocks in topology + // molBlockIndex - index of current block in mtop + for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++) + { + // check if current block contains QM atoms + if (bQMBlock_[molBlockIndex]) + { + // molType - strucutre with current block type + gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type]; + + // start - global index of the first atom in the current block + int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart; + + // loop over all interaction types + for (int ftype = 0; ftype < F_NRE; ftype++) + { + // If this is VSite interaction and ilist is not empty + if (IS_VSITE(ftype) && !molType->ilist[ftype].empty()) + { + // Number of elements in the iatoms array for the current ftype + const int numInteractionElements = NRAL(ftype) + 1; + + // Loop over all interactions of ftype + for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements) + { + // Calculate number of qm atoms in the interaction + int numQm = 0; + // Here k starts from 2 because first atom in the interaction is an actuall vsite index + for (int k = 2; k <= NRAL(ftype); k++) + { + if (isQMAtom(molType->ilist[ftype].iatoms[j + k] + start)) + { + numQm++; + } + } + + // If all atoms froming that virtual site are QM atoms + // then remove classical charge from that virtual site + if (numQm == (NRAL(ftype) - 1)) + { + topInfo_.numVirtualSitesModified++; + topInfo_.totalClassicalChargeOfQMAtoms += + molType->atoms.atom[molType->ilist[ftype].iatoms[j + 1]].q; + molType->atoms.atom[molType->ilist[ftype].iatoms[j + 1]].q = 0.0; + molType->atoms.atom[molType->ilist[ftype].iatoms[j + 1]].qB = 0.0; + } + } + } + } + } + } +} + +} // namespace gmx diff --git a/src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.h b/src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.h new file mode 100644 index 0000000000..64a93ffc24 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.h @@ -0,0 +1,206 @@ +/* + * 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 + * QMMMTopologyPrepocessor class responsible for + * all modificatios of the topology during input pre-processing + * + * \author Dmitry Morozov + * \ingroup module_applied_forces + */ +#ifndef GMX_APPLIED_FORCES_QMMMTOPOLOGYPREPROCESSOR_H +#define GMX_APPLIED_FORCES_QMMMTOPOLOGYPREPROCESSOR_H + +#include +#include + +#include "gromacs/utility/basedefinitions.h" +#include "gromacs/utility/arrayref.h" + +#include "qmmmtypes.h" + +struct gmx_mtop_t; + +namespace gmx +{ + +/*! \internal + * \brief Contains various information about topology modifications + * Used for statistics during topology pre-processing within QMMMTopologyPreprocessor class + */ +struct QMMMTopologyInfo +{ + //! Total number of MM atoms + int numMMAtoms = 0; + //! Total number of QM atoms + int numQMAtoms = 0; + //! Total remaining charge of MM part + real remainingMMCharge = 0.0; + //! Total classical charge removed from QM atoms + real totalClassicalChargeOfQMAtoms = 0.0; + //! Total number of Non-bonded (LJ) exclusions made for QM-QM interactions + int numExclusionsMade = 0; + //! Total number of removed classical Bonds between QM-QM atoms + int numBondsRemoved = 0; + //! Total number of removed classical Angles between QM-QM atoms + int numAnglesRemoved = 0; + //! Total number of removed classical Dihedrals between QM-QM atoms + int numDihedralsRemoved = 0; + //! Total number of removed F_SETTLE between QM-QM atoms + int numSettleRemoved = 0; + //! Total number of empty chemical bonds (F_CONNBONDS) added between QM-QM atoms + int numConnBondsAdded = 0; + //! Total number of virtual sites, that consisting of QM atoms only, which charge has been removed + int numVirtualSitesModified = 0; + //! Total number of constrained bonds within QM subsystem + int numConstrainedBondsInQMSubsystem = 0; + //! Total number of broken bonds between QM and MM atoms (Link Frontier) + int numLinkBonds = 0; +}; + +/*! \internal + * \brief Class implementing gmx_mtop_t QMMM modifications during preprocessing + * 1) Split QM-containing molecules from other molecules in blocks + * 2) Nullify charges on all virtual sites consisting of QM only atoms + * 3) Nullifies charges on all QM atoms + * 4) Excludes LJ interactions between QM atoms + * 5) Builds vector with atomic numbers of all atoms + * 6) Makes F_CONNBOND between atoms within QM region + * 7) Removes angles and settles containing 2 or more QM atoms + * 8) Removes dihedrals containing 3 or more QM atoms + * 9) Builds vector containing pairs of bonded QM - MM atoms (Link frontier) + */ +class QMMMTopologyPreprocessor +{ +public: + /*! \brief Constructor for QMMMTopologyPreprocessor from its parameters + * + * \param[in] qmIndices Array with global indicies of QM atoms + */ + QMMMTopologyPreprocessor(ArrayRef qmIndices); + + /*! \brief Pocesses mtop topology and prepares atomNumbers_ and linkFrontier_ vectors + * Builds topInfo_ containing information about topology modifications + * + * \param[in,out] mtop Topology that needs to be modified + */ + void preprocess(gmx_mtop_t* mtop); + + //! \brief Returns data about modifications made via QMMMTopologyInfo + const QMMMTopologyInfo& topInfo() const; + + //! \brief Returns view of atomic numbers for all atoms in the processed topology + ArrayRef atomNumbers() const; + + //! \brief Returns view of point charges for all atoms in the processed topology + ArrayRef atomCharges() const; + + //! \brief Returns view of the whole Link Frontier for the processed topology + ArrayRef linkFrontier() const; + +private: + //! Retruns true if globalAtomIndex belongs to QM region + bool isQMAtom(index globalAtomIndex); + + /*! \brief Splits QM containing molecules out of MM blocks in topology + * Modifies blocks in topology + * Updates bQMBlock vector containing QM flags of all blocks in modified mtop + */ + void splitQMblocks(gmx_mtop_t* mtop); + + /*! \brief Removes classical charges from QM atoms + * Provides data about removed charge via topInfo_ + */ + void removeQMClassicalCharges(gmx_mtop_t* mtop); + + //! \brief Build exlusion list for LJ interactions between QM atoms + void addQMLJExclusions(gmx_mtop_t* mtop); + + /*! \brief Builds atomNumbers_ vector + * Provides data about total number of QM and MM atoms via topInfo_ + */ + void buildQMMMAtomNumbers(gmx_mtop_t* mtop); + + /*! \brief Modifies pairwise bonded interactions + * Removes any other pairwise bonded interactions between QM-QM atoms + * Creates F_CONNBOND between QM atoms + * Any restraints and constraints will be kept + * Provides data about modifications via topInfo_ + */ + void modifyQMMMTwoCenterInteractions(gmx_mtop_t* mtop); + + /*! \brief Builds link_ vector with pairs of atoms indicting broken QM - MM chemical bonds. + * Also performs search of constrained bonds within QM subsystem. + */ + void buildQMMMLink(gmx_mtop_t* mtop); + + /*! \brief Modifies three-centers interactions (i.e. Angles, Settles) + * Removes any other three-centers bonded interactions including 2 or more QM atoms + * Any restraints and constraints will be kept + * Any F_SETTLE containing QM atoms will be converted to the pair of F_CONNBONDS + * Provides data about modifications via topInfo_ + */ + void modifyQMMMThreeCenterInteractions(gmx_mtop_t* mtop); + + /*! \brief Modifies four-centers interactions + * Removes any other four-centers bonded interactions including 3 or more QM atoms + * Any restraints and constraints will be kept + * Provides data about modifications via topInfo_ + */ + void modifyQMMMFourCenterInteractions(gmx_mtop_t* mtop); + + //! \brief Removes charge from all virtual sites which are consists of only QM atoms + void modifyQMMMVirtualSites(gmx_mtop_t* mtop); + + //! Vector indicating which molblocks have QM atoms + std::vector bQMBlock_; + /*! \brief Global indices of QM atoms; + * The dominant operation is search and we also expect the set of qm atoms to be very small + * relative to the rest, so set should outperform unordered set, i.e. unsorted std::vector. + */ + std::set qmIndices_; + //! Vector with atom numbers for the whole system + std::vector atomNumbers_; + //! Vector with atom point charges for the whole system + std::vector atomCharges_; + //! Vector with pairs of indices defining broken bonds in QMMM + std::vector linkFrontier_; + //! Structure with information about modifications made + QMMMTopologyInfo topInfo_; +}; + +} // namespace gmx + +#endif diff --git a/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt b/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt index 93db354c8b..7c879e151a 100644 --- a/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt +++ b/src/gromacs/applied_forces/qmmm/tests/CMakeLists.txt @@ -35,4 +35,5 @@ gmx_add_unit_test(QMMMAppliedForcesUnitTest qmmm_applied_forces-test CPP_SOURCE_FILES qmmminputgenerator.cpp + qmmmtopologypreprocessor.cpp ) diff --git a/src/gromacs/applied_forces/qmmm/tests/qmmmtopologypreprocessor.cpp b/src/gromacs/applied_forces/qmmm/tests/qmmmtopologypreprocessor.cpp new file mode 100644 index 0000000000..ed5814c040 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/qmmmtopologypreprocessor.cpp @@ -0,0 +1,244 @@ +/* + * 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/qmmmtopologypreprocessor.h" + +#include + +#include + +#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 QMMMTopologyPreprocessorTest : public ::testing::Test +{ +public: + /*! \brief Generates tpr file from *.top and *.gro existing in the simulation database directory + * and loads gmx_mtop_t from it + */ + void makeMtopFromFile(const std::string& fileName, const std::string& mdpContent) + { + const std::string simData = gmx::test::TestFileManager::getTestSimulationDatabaseDirectory(); + + // Generate empty mdp file + const std::string mdpInputFileName = fileManager_.getTemporaryFilePath(fileName + ".mdp"); + gmx::TextWriter::writeFileFromString(mdpInputFileName, mdpContent); + + // Generate tpr file + const std::string tprName = fileManager_.getTemporaryFilePath(fileName + ".tpr"); + { + gmx::test::CommandLine caller; + caller.append("grompp"); + caller.addOption("-f", mdpInputFileName); + caller.addOption("-p", gmx::Path::join(simData, fileName + ".top")); + caller.addOption("-c", gmx::Path::join(simData, fileName + ".gro")); + caller.addOption("-o", tprName); + ASSERT_EQ(0, gmx_grompp(caller.argc(), caller.argv())); + } + + // Load topology + bool fullTopology; + PbcType pbcType; + matrix box; + readConfAndTopology(tprName.c_str(), &fullTopology, &mtop_, &pbcType, nullptr, nullptr, box); + } + +protected: + gmx::test::TestFileManager fileManager_; + std::vector qmIndices_; + gmx_mtop_t mtop_; +}; + +class QMMMTopologyPreprocessorChecker : public gmx::test::TestReferenceChecker +{ +public: + QMMMTopologyPreprocessorChecker(const gmx::test::TestReferenceChecker& rootChecker) : + gmx::test::TestReferenceChecker(rootChecker) + { + // Tolerance of all charges and vectors should be 1E-3 + setDefaultTolerance(gmx::test::absoluteTolerance(0.001)); + } + + void checkAll(const QMMMTopologyInfo& info, const gmx_mtop_t& mtop) + { + checkInteger(info.numMMAtoms, "Number of MM atoms"); + checkInteger(info.numQMAtoms, "Number of QM atoms"); + checkReal(info.remainingMMCharge, "MM charge"); + checkReal(info.totalClassicalChargeOfQMAtoms, "QM charge"); + checkInteger(info.numExclusionsMade, "Exclusions Made"); + checkInteger(info.numBondsRemoved, "Removed bonds"); + checkInteger(info.numAnglesRemoved, "Removed angles"); + checkInteger(info.numDihedralsRemoved, "Removed dihedrals"); + checkInteger(info.numSettleRemoved, "Removed settles"); + checkInteger(info.numConnBondsAdded, "Generated CONNBONDS"); + checkInteger(info.numVirtualSitesModified, "Removed vsites"); + checkInteger(info.numConstrainedBondsInQMSubsystem, "QM Constraints Found"); + checkInteger(info.numLinkBonds, "Link-atom sites"); + checkInteger(mtop.molblock.size(), "Molblocks in topology"); + } +}; + +TEST_F(QMMMTopologyPreprocessorTest, CanConstruct) +{ + qmIndices_ = { 0, 1, 2 }; + EXPECT_NO_THROW(QMMMTopologyPreprocessor topPrep(qmIndices_)); +} + +TEST_F(QMMMTopologyPreprocessorTest, FourWatersFirstQMNoLink) +{ + // Reference input 4x SPCE waters from database 4waters.top (first one QM) and no Link atoms + makeMtopFromFile("4water", ""); + qmIndices_ = { 0, 1, 2 }; + + QMMMTopologyPreprocessor topPrep(qmIndices_); + topPrep.preprocess(&mtop_); + + // Get data about changes and check it + QMMMTopologyInfo info = topPrep.topInfo(); + + gmx::test::TestReferenceData data; + QMMMTopologyPreprocessorChecker checker(data.rootChecker()); + checker.checkAll(info, mtop_); +} + +TEST_F(QMMMTopologyPreprocessorTest, FourWatersSeondAndForthQMNoLink) +{ + // Reference input 4x SPCE waters from database 4waters.top (second and forth are QM) and no Link atoms + makeMtopFromFile("4water", ""); + qmIndices_ = { 3, 4, 5, 9, 10, 11 }; + + QMMMTopologyPreprocessor topPrep(qmIndices_); + topPrep.preprocess(&mtop_); + + // Get data about changes and check it + QMMMTopologyInfo info = topPrep.topInfo(); + + gmx::test::TestReferenceData data; + QMMMTopologyPreprocessorChecker checker(data.rootChecker()); + checker.checkAll(info, mtop_); +} + +TEST_F(QMMMTopologyPreprocessorTest, FourWatersFirstQMWithLink) +{ + // Reference input 4x SPCE waters from database 4waters.top (first one QM) with Link atom + makeMtopFromFile("4water", ""); + qmIndices_ = { 0, 1 }; + + QMMMTopologyPreprocessor topPrep(qmIndices_); + topPrep.preprocess(&mtop_); + + // Get data about changes and check it + QMMMTopologyInfo info = topPrep.topInfo(); + + gmx::test::TestReferenceData data; + QMMMTopologyPreprocessorChecker checker(data.rootChecker()); + checker.checkAll(info, mtop_); +} + +TEST_F(QMMMTopologyPreprocessorTest, AlanineDipeptideWithLinksNoConstraints) +{ + // Reference input alanine_vacuo.top + makeMtopFromFile("alanine_vacuo", ""); + qmIndices_ = { 8, 9, 10, 11, 12, 13 }; + + QMMMTopologyPreprocessor topPrep(qmIndices_); + topPrep.preprocess(&mtop_); + + // Get data about changes and check it + QMMMTopologyInfo info = topPrep.topInfo(); + + gmx::test::TestReferenceData data; + QMMMTopologyPreprocessorChecker checker(data.rootChecker()); + checker.checkAll(info, mtop_); +} + +TEST_F(QMMMTopologyPreprocessorTest, AlanineDipeptideWithLinksWithConstraints) +{ + // Reference input alanine_vacuo.top with constraints=all-bonds + makeMtopFromFile("alanine_vacuo", "constraints = all-bonds"); + qmIndices_ = { 8, 9, 10, 11, 12, 13 }; + + QMMMTopologyPreprocessor topPrep(qmIndices_); + topPrep.preprocess(&mtop_); + + // Get data about changes and check it + QMMMTopologyInfo info = topPrep.topInfo(); + + gmx::test::TestReferenceData data; + QMMMTopologyPreprocessorChecker checker(data.rootChecker()); + checker.checkAll(info, mtop_); +} + +TEST_F(QMMMTopologyPreprocessorTest, RemovingQMVsites) +{ + // Reference input vistes_test.top + makeMtopFromFile("vsite_test", ""); + qmIndices_ = { 0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15 }; + + QMMMTopologyPreprocessor topPrep(qmIndices_); + topPrep.preprocess(&mtop_); + + // Get data about changes and check it + QMMMTopologyInfo info = topPrep.topInfo(); + + gmx::test::TestReferenceData data; + QMMMTopologyPreprocessorChecker checker(data.rootChecker()); + checker.checkAll(info, mtop_); +} + +} // namespace gmx diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksNoConstraints.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksNoConstraints.xml new file mode 100644 index 0000000000..4816ff4069 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksNoConstraints.xml @@ -0,0 +1,18 @@ + + + + 16 + 6 + -0.11439999999999996 + 0.1144 + 6 + 8 + 11 + 9 + 0 + 5 + 0 + 0 + 2 + 1 + diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksWithConstraints.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksWithConstraints.xml new file mode 100644 index 0000000000..76fa2b49b3 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_AlanineDipeptideWithLinksWithConstraints.xml @@ -0,0 +1,18 @@ + + + + 16 + 6 + -0.11439999999999996 + 0.1144 + 6 + 3 + 11 + 9 + 0 + 0 + 0 + 5 + 2 + 1 + diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMNoLink.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMNoLink.xml new file mode 100644 index 0000000000..82a39dcf50 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMNoLink.xml @@ -0,0 +1,18 @@ + + + + 9 + 3 + 0 + 0 + 3 + 0 + 0 + 0 + 1 + 2 + 0 + 0 + 0 + 2 + diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMWithLink.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMWithLink.xml new file mode 100644 index 0000000000..7c15142d28 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersFirstQMWithLink.xml @@ -0,0 +1,18 @@ + + + + 10 + 2 + 0.40999999999999998 + -0.40999999999999998 + 2 + 0 + 0 + 0 + 1 + 2 + 0 + 0 + 1 + 2 + diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersSeondAndForthQMNoLink.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersSeondAndForthQMNoLink.xml new file mode 100644 index 0000000000..37806c3141 --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_FourWatersSeondAndForthQMNoLink.xml @@ -0,0 +1,18 @@ + + + + 6 + 6 + 0 + 0 + 6 + 0 + 0 + 0 + 2 + 4 + 0 + 0 + 0 + 4 + diff --git a/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_RemovingQMVsites.xml b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_RemovingQMVsites.xml new file mode 100644 index 0000000000..a532f1a5ed --- /dev/null +++ b/src/gromacs/applied_forces/qmmm/tests/refdata/QMMMTopologyPreprocessorTest_RemovingQMVsites.xml @@ -0,0 +1,18 @@ + + + + 8 + 16 + 0 + 0 + 16 + 28 + 14 + 13 + 0 + 15 + 8 + 0 + 0 + 1 + -- 2.22.0