gmx_add_libgromacs_sources(
qmmminputgenerator.cpp
+ qmmmtopologypreprocessor.cpp
)
if (BUILD_TESTING)
--- /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 QMMMTopologyPreprocessor
+ *
+ * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
+ * \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<const index> 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<const int> QMMMTopologyPreprocessor::atomNumbers() const
+{
+ return atomNumbers_;
+}
+
+ArrayRef<const real> QMMMTopologyPreprocessor::atomCharges() const
+{
+ return atomCharges_;
+}
+
+ArrayRef<const LinkFrontier> 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<int> 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<gmx_moltype_t> 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
--- /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
+ * QMMMTopologyPrepocessor class responsible for
+ * all modificatios of the topology during input pre-processing
+ *
+ * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
+ * \ingroup module_applied_forces
+ */
+#ifndef GMX_APPLIED_FORCES_QMMMTOPOLOGYPREPROCESSOR_H
+#define GMX_APPLIED_FORCES_QMMMTOPOLOGYPREPROCESSOR_H
+
+#include <set>
+#include <string>
+
+#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<const index> 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<const int> atomNumbers() const;
+
+ //! \brief Returns view of point charges for all atoms in the processed topology
+ ArrayRef<const real> atomCharges() const;
+
+ //! \brief Returns view of the whole Link Frontier for the processed topology
+ ArrayRef<const LinkFrontier> 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<bool> 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<int> qmIndices_;
+ //! Vector with atom numbers for the whole system
+ std::vector<int> atomNumbers_;
+ //! Vector with atom point charges for the whole system
+ std::vector<real> atomCharges_;
+ //! Vector with pairs of indices defining broken bonds in QMMM
+ std::vector<LinkFrontier> linkFrontier_;
+ //! Structure with information about modifications made
+ QMMMTopologyInfo topInfo_;
+};
+
+} // namespace gmx
+
+#endif
gmx_add_unit_test(QMMMAppliedForcesUnitTest qmmm_applied_forces-test
CPP_SOURCE_FILES
qmmminputgenerator.cpp
+ qmmmtopologypreprocessor.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/qmmmtopologypreprocessor.h"
+
+#include <vector>
+
+#include <gtest/gtest.h>
+
+#include "gromacs/fileio/confio.h"
+#include "gromacs/gmxpreprocess/grompp.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/topology/mtop_lookup.h"
+#include "gromacs/topology/mtop_util.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/arrayref.h"
+#include "gromacs/utility/path.h"
+#include "gromacs/utility/textwriter.h"
+#include "testutils/cmdlinetest.h"
+#include "testutils/refdata.h"
+#include "testutils/testasserts.h"
+#include "testutils/testfilemanager.h"
+
+namespace gmx
+{
+
+class 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<index> 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
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Int Name="Number of MM atoms">16</Int>
+ <Int Name="Number of QM atoms">6</Int>
+ <Real Name="MM charge">-0.11439999999999996</Real>
+ <Real Name="QM charge">0.1144</Real>
+ <Int Name="Exclusions Made">6</Int>
+ <Int Name="Removed bonds">8</Int>
+ <Int Name="Removed angles">11</Int>
+ <Int Name="Removed dihedrals">9</Int>
+ <Int Name="Removed settles">0</Int>
+ <Int Name="Generated CONNBONDS">5</Int>
+ <Int Name="Removed vsites">0</Int>
+ <Int Name="QM Constraints Found">0</Int>
+ <Int Name="Link-atom sites">2</Int>
+ <Int Name="Molblocks in topology">1</Int>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Int Name="Number of MM atoms">16</Int>
+ <Int Name="Number of QM atoms">6</Int>
+ <Real Name="MM charge">-0.11439999999999996</Real>
+ <Real Name="QM charge">0.1144</Real>
+ <Int Name="Exclusions Made">6</Int>
+ <Int Name="Removed bonds">3</Int>
+ <Int Name="Removed angles">11</Int>
+ <Int Name="Removed dihedrals">9</Int>
+ <Int Name="Removed settles">0</Int>
+ <Int Name="Generated CONNBONDS">0</Int>
+ <Int Name="Removed vsites">0</Int>
+ <Int Name="QM Constraints Found">5</Int>
+ <Int Name="Link-atom sites">2</Int>
+ <Int Name="Molblocks in topology">1</Int>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Int Name="Number of MM atoms">9</Int>
+ <Int Name="Number of QM atoms">3</Int>
+ <Real Name="MM charge">0</Real>
+ <Real Name="QM charge">0</Real>
+ <Int Name="Exclusions Made">3</Int>
+ <Int Name="Removed bonds">0</Int>
+ <Int Name="Removed angles">0</Int>
+ <Int Name="Removed dihedrals">0</Int>
+ <Int Name="Removed settles">1</Int>
+ <Int Name="Generated CONNBONDS">2</Int>
+ <Int Name="Removed vsites">0</Int>
+ <Int Name="QM Constraints Found">0</Int>
+ <Int Name="Link-atom sites">0</Int>
+ <Int Name="Molblocks in topology">2</Int>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Int Name="Number of MM atoms">10</Int>
+ <Int Name="Number of QM atoms">2</Int>
+ <Real Name="MM charge">0.40999999999999998</Real>
+ <Real Name="QM charge">-0.40999999999999998</Real>
+ <Int Name="Exclusions Made">2</Int>
+ <Int Name="Removed bonds">0</Int>
+ <Int Name="Removed angles">0</Int>
+ <Int Name="Removed dihedrals">0</Int>
+ <Int Name="Removed settles">1</Int>
+ <Int Name="Generated CONNBONDS">2</Int>
+ <Int Name="Removed vsites">0</Int>
+ <Int Name="QM Constraints Found">0</Int>
+ <Int Name="Link-atom sites">1</Int>
+ <Int Name="Molblocks in topology">2</Int>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Int Name="Number of MM atoms">6</Int>
+ <Int Name="Number of QM atoms">6</Int>
+ <Real Name="MM charge">0</Real>
+ <Real Name="QM charge">0</Real>
+ <Int Name="Exclusions Made">6</Int>
+ <Int Name="Removed bonds">0</Int>
+ <Int Name="Removed angles">0</Int>
+ <Int Name="Removed dihedrals">0</Int>
+ <Int Name="Removed settles">2</Int>
+ <Int Name="Generated CONNBONDS">4</Int>
+ <Int Name="Removed vsites">0</Int>
+ <Int Name="QM Constraints Found">0</Int>
+ <Int Name="Link-atom sites">0</Int>
+ <Int Name="Molblocks in topology">4</Int>
+</ReferenceData>
--- /dev/null
+<?xml version="1.0"?>
+<?xml-stylesheet type="text/xsl" href="referencedata.xsl"?>
+<ReferenceData>
+ <Int Name="Number of MM atoms">8</Int>
+ <Int Name="Number of QM atoms">16</Int>
+ <Real Name="MM charge">0</Real>
+ <Real Name="QM charge">0</Real>
+ <Int Name="Exclusions Made">16</Int>
+ <Int Name="Removed bonds">28</Int>
+ <Int Name="Removed angles">14</Int>
+ <Int Name="Removed dihedrals">13</Int>
+ <Int Name="Removed settles">0</Int>
+ <Int Name="Generated CONNBONDS">15</Int>
+ <Int Name="Removed vsites">8</Int>
+ <Int Name="QM Constraints Found">0</Int>
+ <Int Name="Link-atom sites">0</Int>
+ <Int Name="Molblocks in topology">1</Int>
+</ReferenceData>
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff