src/gromacs/applied_forces/qmmm/qmmmtopologypreprocessor.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2021, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35 /*! \internal \file
  36  * \brief
  37  * Implements QMMMTopologyPreprocessor
  38  *
  39  * \author Dmitry Morozov <dmitry.morozov@jyu.fi>
  40  * \ingroup module_applied_forces
  41  */
  42 #include "gmxpre.h"
  43
  44 #include "qmmmtopologypreprocessor.h"
  45
  46 #include "gromacs/selection/indexutil.h"
  47 #include "gromacs/topology/mtop_lookup.h"
  48 #include "gromacs/topology/mtop_util.h"
  49 #include "gromacs/topology/topology.h"
  50 #include "gromacs/utility/exceptions.h"
  51 #include "gromacs/utility/fatalerror.h"
  52
  53 namespace gmx
  54 {
  55
  56 QMMMTopologyPreprocessor::QMMMTopologyPreprocessor(ArrayRef<const index> qmIndices) :
  57     qmIndices_(qmIndices.begin(), qmIndices.end())
  58 {
  59 }
  60
  61 void QMMMTopologyPreprocessor::preprocess(gmx_mtop_t* mtop)
  62 {
  63     // 1) Split QM-containing molecules from other molecules in blocks
  64     splitQMblocks(mtop);
  65
  66     // 2) Nullify charges on all virtual sites consisting of QM only atoms
  67     modifyQMMMVirtualSites(mtop);
  68
  69     // 3) Nullify charges on all QM atoms
  70     removeQMClassicalCharges(mtop);
  71
  72     // 4) Exclude LJ interactions between QM atoms
  73     addQMLJExclusions(mtop);
  74
  75     // 5) Build atomNumbers vector with atomic numbers of all atoms
  76     buildQMMMAtomNumbers(mtop);
  77
  78     // 6) Make F_CONNBOND between atoms within QM region
  79     modifyQMMMTwoCenterInteractions(mtop);
  80
  81     // 7) Remove angles and settles containing 2 or more QM atoms
  82     modifyQMMMThreeCenterInteractions(mtop);
  83
  84     // 8) Remove dihedrals containing 3 or more QM atoms
  85     modifyQMMMFourCenterInteractions(mtop);
  86
  87     // 9) Build vector containing pairs of bonded QM - MM atoms (Link frontier)
  88     buildQMMMLink(mtop);
  89
  90     // finalize topology
  91     mtop->finalize();
  92 }
  93
  94 const QMMMTopologyInfo& QMMMTopologyPreprocessor::topInfo() const
  95 {
  96     return topInfo_;
  97 }
  98
  99 ArrayRef<const int> QMMMTopologyPreprocessor::atomNumbers() const
 100 {
 101     return atomNumbers_;
 102 }
 103
 104 ArrayRef<const real> QMMMTopologyPreprocessor::atomCharges() const
 105 {
 106     return atomCharges_;
 107 }
 108
 109 ArrayRef<const LinkFrontier> QMMMTopologyPreprocessor::linkFrontier() const
 110 {
 111     return linkFrontier_;
 112 }
 113
 114 bool QMMMTopologyPreprocessor::isQMAtom(index globalAtomIndex)
 115 {
 116     return (qmIndices_.find(globalAtomIndex) != qmIndices_.end());
 117 }
 118
 119 void QMMMTopologyPreprocessor::splitQMblocks(gmx_mtop_t* mtop)
 120 {
 121
 122     // Global counter of atoms
 123     index iAt = 0;
 124
 125     /* Counter of molecules point to the specific moltype
 126      * i.e molblock 0 has 2 molecules have moltype 0 and molblock 2 has 1 additional molecule of type 0
 127      * then will be numMoleculesOfType[0] = 2 + 1 = 3
 128      * That counter is needed to decide whether we should create a new moltype for the molecule contatining QM atoms
 129      * or we could modify existing moltype if there is only one molecule of that type
 130      */
 131     std::vector<int> numMoleculesOfType(mtop->moltype.size());
 132     for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++)
 133     {
 134         numMoleculesOfType[mtop->molblock[molBlockIndex].type] += mtop->molblock[molBlockIndex].nmol;
 135     }
 136
 137     // Loop over all blocks in topology
 138     // molBlockIndex - current index of block in mtop
 139     for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++)
 140     {
 141         // Initialize block as non-QM first
 142         bQMBlock_.push_back(false);
 143
 144         // Pointer to current block
 145         gmx_molblock_t* molBlock = &mtop->molblock[molBlockIndex];
 146
 147         // Number of atoms in all molecules of current block
 148         const int numAtomsInMolecule = mtop->moltype[molBlock->type].atoms.nr;
 149
 150         // Loop over all molecules in molBlock
 151         // mol - index of current molecule in molBlock
 152         for (int mol = 0; mol < molBlock->nmol; mol++)
 153         {
 154
 155             // search for QM atoms in current molecule
 156             bool bQMMM = false;
 157             for (int i = 0; i < numAtomsInMolecule; i++)
 158             {
 159                 if (isQMAtom(iAt))
 160                 {
 161                     bQMMM = true;
 162                 }
 163                 iAt++;
 164             }
 165
 166             // Apparently current molecule (molBlockIndex, mol) contains QM atoms
 167             // We should split it from the current block and create new blocks
 168             // For that molecule and all molecules after it new block will be created
 169             if (bQMMM)
 170             {
 171                 // if this block contains only 1 molecule, then splitting not needed
 172                 if (molBlock->nmol > 1)
 173                 {
 174                     // If current molecule is not the first in the block,
 175                     // then we need to first create block for it and all molecule before it
 176                     if (mol > 0)
 177                     {
 178                         // Split the molblock at this molecule
 179                         auto pos = mtop->molblock.begin() + molBlockIndex + 1;
 180                         mtop->molblock.insert(pos, mtop->molblock[molBlockIndex]);
 181                         mtop->molblock[molBlockIndex].nmol = mol;
 182                         mtop->molblock[molBlockIndex + 1].nmol -= mol;
 183                         bQMBlock_[molBlockIndex] = false;
 184                         bQMBlock_.push_back(true);
 185                         molBlockIndex++;
 186                         molBlock = &mtop->molblock[molBlockIndex];
 187                     }
 188
 189                     // If current molecule is not the only one in new block,
 190                     // Then we split new block after that molecule
 191                     if (molBlock->nmol > 1)
 192                     {
 193                         auto pos = mtop->molblock.begin() + molBlockIndex + 1;
 194                         mtop->molblock.insert(pos, mtop->molblock[molBlockIndex]);
 195                         molBlock                           = &mtop->molblock[molBlockIndex];
 196                         mtop->molblock[molBlockIndex].nmol = 1;
 197                         mtop->molblock[molBlockIndex + 1].nmol -= 1;
 198                         bQMBlock_[molBlockIndex] = true;
 199                     }
 200                 }
 201                 else
 202                 {
 203                     bQMBlock_[molBlockIndex] = true;
 204                 }
 205
 206                 // Create a copy of a moltype for a molecule
 207                 // containing QM atoms and append it in the end of the moltype vector
 208                 // if there is only 1 molecule pointing to that type then skip that step
 209                 if (numMoleculesOfType[molBlock->type] > 1)
 210                 {
 211                     // Here comes a huge piece of "not so good" code, because of deleted operator= from gmx_moltype_t
 212                     std::vector<gmx_moltype_t> temp(mtop->moltype.size());
 213                     for (size_t i = 0; i < mtop->moltype.size(); ++i)
 214                     {
 215                         copy_moltype(&mtop->moltype[i], &temp[i]);
 216                     }
 217                     mtop->moltype.resize(temp.size() + 1);
 218                     for (size_t i = 0; i < temp.size(); ++i)
 219                     {
 220                         copy_moltype(&temp[i], &mtop->moltype[i]);
 221                     }
 222                     copy_moltype(&mtop->moltype[molBlock->type], &mtop->moltype.back());
 223
 224                     // Set the molecule type for the QMMM molblock
 225                     molBlock->type = mtop->moltype.size() - 1;
 226                 }
 227             }
 228         }
 229     }
 230
 231     // Call finalize() to rebuild Block Indicies or else atoms lookup will fail
 232     mtop->finalize();
 233 }
 234
 235 void QMMMTopologyPreprocessor::removeQMClassicalCharges(gmx_mtop_t* mtop)
 236 {
 237     // Loop over all atoms and remove charge if they are QM atoms.
 238     // Sum-up total removed charge and remaning charge on MM atoms
 239     // Build atomCharges_ vector
 240     int molBlockIndex = 0;
 241     for (int i = 0; i < mtop->natoms; i++)
 242     {
 243         int indexInMolecule;
 244         mtopGetMolblockIndex(*mtop, i, &molBlockIndex, nullptr, &indexInMolecule);
 245         t_atom* atom = &mtop->moltype[mtop->molblock[molBlockIndex].type].atoms.atom[indexInMolecule];
 246         if (isQMAtom(i))
 247         {
 248             topInfo_.totalClassicalChargeOfQMAtoms += atom->q;
 249             atom->q  = 0.0;
 250             atom->qB = 0.0;
 251         }
 252         else
 253         {
 254             topInfo_.remainingMMCharge += atom->q;
 255         }
 256
 257         atomCharges_.push_back(atom->q);
 258     }
 259 }
 260
 261 void QMMMTopologyPreprocessor::addQMLJExclusions(gmx_mtop_t* mtop)
 262 {
 263     // Add all QM atoms to the mtop->intermolecularExclusionGroup
 264     mtop->intermolecularExclusionGroup.reserve(mtop->intermolecularExclusionGroup.size()
 265                                                + qmIndices_.size());
 266     for (auto i : qmIndices_)
 267     {
 268         mtop->intermolecularExclusionGroup.push_back(i);
 269         topInfo_.numExclusionsMade++;
 270     }
 271 }
 272
 273 void QMMMTopologyPreprocessor::buildQMMMAtomNumbers(gmx_mtop_t* mtop)
 274 {
 275     // Save to atomNumbers_ atom numbers of all atoms
 276     AtomIterator atoms(*mtop);
 277     while (atoms->globalAtomNumber() < mtop->natoms)
 278     {
 279         // Check if we have valid atomnumbers
 280         if (atoms->atom().atomnumber < 0)
 281         {
 282             gmx_fatal(FARGS,
 283                       "Atoms %d does not have atomic number needed for QMMM. Check atomtypes "
 284                       "section in your topology or forcefield.",
 285                       atoms->globalAtomNumber());
 286         }
 287
 288         atomNumbers_.push_back(atoms->atom().atomnumber);
 289         atoms++;
 290     }
 291
 292     // Save in topInfo_ number of QM and MM atoms
 293     topInfo_.numQMAtoms += gmx::ssize(qmIndices_);
 294     topInfo_.numMMAtoms += mtop->natoms - gmx::ssize(qmIndices_);
 295 }
 296
 297 void QMMMTopologyPreprocessor::modifyQMMMTwoCenterInteractions(gmx_mtop_t* mtop)
 298 {
 299     // Loop over all blocks in topology
 300     // molBlockIndex - index of current block in mtop
 301     for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++)
 302     {
 303         // check if current block contains QM atoms
 304         if (bQMBlock_[molBlockIndex])
 305         {
 306             // molType - strucutre with current block type
 307             gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type];
 308
 309             // start - global index if the first atom in the current block
 310             int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart;
 311
 312             // loop over all interaction types
 313             for (int ftype = 0; ftype < F_NRE; ftype++)
 314             {
 315                 // If not bonded interaction or F_CONNBONDS, or some form of Restraints,
 316                 // or not pair interaction, or no interactions of that type: then go the next type
 317                 if (!(interaction_function[ftype].flags & IF_BOND) || ftype == F_CONNBONDS
 318                     || ftype == F_RESTRBONDS || ftype == F_HARMONIC || ftype == F_DISRES || ftype == F_ORIRES
 319                     || ftype == F_ANGRESZ || NRAL(ftype) != 2 || molType->ilist[ftype].empty())
 320                 {
 321                     continue;
 322                 }
 323
 324                 // Number of elements in the iatoms array for the current ftype
 325                 const int numInteractionElements = NRAL(ftype) + 1;
 326
 327                 // Buffer for preserving interactions that are retained
 328                 AtomGroupIndices iatomsBuf;
 329
 330                 // Loop over all interactions of ftype
 331                 for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements)
 332                 {
 333
 334                     // If both atoms are QM and it is IF_CHEMBOND then convert it to F_CONNBONDS
 335                     if (isQMAtom(molType->ilist[ftype].iatoms[j + 1] + start)
 336                         && isQMAtom(molType->ilist[ftype].iatoms[j + 2] + start))
 337                     {
 338
 339                         // Add chemical bond to the F_CONNBONDS (bond type 5)
 340                         if (IS_CHEMBOND(ftype))
 341                         {
 342                             // Bond type is not used in F_CONNBONDS, so for generated bonds we set it to -1
 343                             const int connBondsType = -1;
 344
 345                             // Add new CONNBOND between atoms
 346                             molType->ilist[F_CONNBONDS].iatoms.push_back(connBondsType);
 347                             molType->ilist[F_CONNBONDS].iatoms.push_back(
 348                                     molType->ilist[ftype].iatoms[j + 1]);
 349                             molType->ilist[F_CONNBONDS].iatoms.push_back(
 350                                     molType->ilist[ftype].iatoms[j + 2]);
 351
 352                             topInfo_.numConnBondsAdded++;
 353                         }
 354
 355                         // Since the current bond is not copied into iatomsBuf it will be removed from the topology
 356                         topInfo_.numBondsRemoved++;
 357                     }
 358                     else
 359                     {
 360                         // If one of atoms is not QM then copy interaction into iatomsBuf
 361                         for (int k = 0; k < numInteractionElements; k++)
 362                         {
 363                             iatomsBuf.push_back(molType->ilist[ftype].iatoms[k + j]);
 364                         }
 365                     }
 366                 }
 367
 368                 // Swap iatomsBuf and molType->ilist[ftype].iatoms
 369                 molType->ilist[ftype].iatoms.swap(iatomsBuf);
 370             }
 371         }
 372     }
 373 }
 374
 375 void QMMMTopologyPreprocessor::buildQMMMLink(gmx_mtop_t* mtop)
 376 {
 377     // Loop over all blocks in topology
 378     // molBlockIndex - index of current block in mtop
 379     for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++)
 380     {
 381         // check if current block contains QM atoms
 382         if (bQMBlock_[molBlockIndex])
 383         {
 384             // molType - strucutre with current block type
 385             gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type];
 386
 387             // start - gloabl index of the first atom in the current block
 388             int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart;
 389
 390             // loop over all interaction types
 391             for (int ftype = 0; ftype < F_NRE; ftype++)
 392             {
 393                 // If not chemical bond interaction or not pair interaction
 394                 // or no interactions of that type: then skip current ftype
 395                 if (!(interaction_function[ftype].flags & IF_CHEMBOND) || NRAL(ftype) != 2
 396                     || molType->ilist[ftype].empty())
 397                 {
 398                     continue;
 399                 }
 400
 401                 // Number of elements in the iatoms array for the current ftype
 402                 const int numInteractionElements = NRAL(ftype) + 1;
 403
 404                 // Loop over all interactions of ftype
 405                 for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements)
 406                 {
 407                     // Global indexes of atoms involved into the interaction
 408                     int a1 = molType->ilist[ftype].iatoms[j + 1] + start;
 409                     int a2 = molType->ilist[ftype].iatoms[j + 2] + start;
 410
 411                     // Update Link Frontier List if one of the atoms QM and one MM
 412                     if (isQMAtom(a1) && !isQMAtom(a2))
 413                     {
 414                         linkFrontier_.push_back({ a1, a2 });
 415                         topInfo_.numLinkBonds++;
 416                     }
 417                     if (isQMAtom(a2) && !isQMAtom(a1))
 418                     {
 419                         linkFrontier_.push_back({ a2, a1 });
 420                         topInfo_.numLinkBonds++;
 421                     }
 422
 423                     // Check if it is constrained bond within QM subsystem
 424                     if (isQMAtom(a2) && isQMAtom(a1) && (interaction_function[ftype].flags & IF_CONSTRAINT))
 425                     {
 426                         topInfo_.numConstrainedBondsInQMSubsystem++;
 427                     }
 428                 }
 429             }
 430         }
 431     }
 432 }
 433
 434 void QMMMTopologyPreprocessor::modifyQMMMThreeCenterInteractions(gmx_mtop_t* mtop)
 435 {
 436     // Loop over all blocks in topology
 437     // molBlockIndex - index of current block in mtop
 438     for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++)
 439     {
 440         // check if current block contains QM atoms
 441         if (bQMBlock_[molBlockIndex])
 442         {
 443             // molType - strucutre with current block type
 444             gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type];
 445
 446             // start - global index of the first atom in the current block
 447             int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart;
 448
 449             // loop over all interaction types
 450             for (int ftype = 0; ftype < F_NRE; ftype++)
 451             {
 452                 // If not bonded interaction or Restraints
 453                 // or not three-particle interaction or no interactions of that type
 454                 // and not Settle: then go the next type
 455                 if ((!(interaction_function[ftype].flags & IF_BOND) || ftype == F_RESTRANGLES
 456                      || NRAL(ftype) != 3 || molType->ilist[ftype].empty())
 457                     && (ftype != F_SETTLE))
 458                 {
 459                     continue;
 460                 }
 461
 462                 // Number of elements in the iatoms array for the current ftype
 463                 const int numInteractionElements = NRAL(ftype) + 1;
 464
 465                 // Buffer for preserving interactions that are retained
 466                 AtomGroupIndices iatomsBuf;
 467
 468                 // Loop over all interactions of ftype
 469                 for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements)
 470                 {
 471                     // Calculate number of qm atoms in the interaction
 472                     int numQm = 0;
 473                     for (int k = 1; k <= NRAL(ftype); k++)
 474                     {
 475                         if (isQMAtom(molType->ilist[ftype].iatoms[j + k] + start))
 476                         {
 477                             numQm++;
 478                         }
 479                     }
 480
 481                     // If at least 2 atoms are QM then remove interaction
 482                     if (numQm >= 2)
 483                     {
 484                         // If this is SETTLE then replace it with two F_CONNBONDS
 485                         if (ftype == F_SETTLE)
 486                         {
 487                             // Bond type is not used in F_CONNBONDS, so for generated bonds we set it to -1
 488                             const int connBondsType = -1;
 489
 490                             // Add CONNBOND between atoms 1 and 2 first
 491                             molType->ilist[F_CONNBONDS].iatoms.push_back(connBondsType);
 492                             molType->ilist[F_CONNBONDS].iatoms.push_back(
 493                                     molType->ilist[ftype].iatoms[j + 1]);
 494                             molType->ilist[F_CONNBONDS].iatoms.push_back(
 495                                     molType->ilist[ftype].iatoms[j + 2]);
 496                             // Then between atoms 1 and 3
 497                             molType->ilist[F_CONNBONDS].iatoms.push_back(connBondsType);
 498                             molType->ilist[F_CONNBONDS].iatoms.push_back(
 499                                     molType->ilist[ftype].iatoms[j + 1]);
 500                             molType->ilist[F_CONNBONDS].iatoms.push_back(
 501                                     molType->ilist[ftype].iatoms[j + 3]);
 502
 503                             topInfo_.numConnBondsAdded += 2;
 504                             topInfo_.numSettleRemoved++;
 505                         }
 506                         else
 507                         {
 508                             // If it is normal angle then it should be just removed
 509                             topInfo_.numAnglesRemoved++;
 510                         }
 511                     }
 512                     else
 513                     {
 514                         // If several (>1) atoms is not QM then preserve interaction in the iatomsBuf
 515                         for (int k = 0; k < numInteractionElements; k++)
 516                         {
 517                             iatomsBuf.push_back(molType->ilist[ftype].iatoms[k + j]);
 518                         }
 519                     }
 520                 }
 521
 522                 // Swap iatomsBuf and molType->ilist[ftype].iatoms
 523                 molType->ilist[ftype].iatoms.swap(iatomsBuf);
 524             }
 525         }
 526     }
 527 }
 528
 529
 530 void QMMMTopologyPreprocessor::modifyQMMMFourCenterInteractions(gmx_mtop_t* mtop)
 531 {
 532     // Loop over all blocks in topology
 533     // molBlockIndex - index of current block in mtop
 534     for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++)
 535     {
 536         // check if current block contains QM atoms
 537         if (bQMBlock_[molBlockIndex])
 538         {
 539             // molType - strucutre with current block type
 540             gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type];
 541
 542             // start - global index of the first atom in the current block
 543             int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart;
 544
 545             // loop over all interaction types
 546             for (int ftype = 0; ftype < F_NRE; ftype++)
 547             {
 548                 // If not bonded interaction or Restraints
 549                 // or not four-particle interaction or no interactions of that type: then go the next type
 550                 if (!(interaction_function[ftype].flags & IF_BOND) || ftype == F_RESTRDIHS
 551                     || NRAL(ftype) != 4 || molType->ilist[ftype].empty())
 552                 {
 553                     continue;
 554                 }
 555
 556                 // Number of elements in the iatoms array for the current ftype
 557                 const int numInteractionElements = NRAL(ftype) + 1;
 558
 559                 // Buffer for preserving interactions that are retained
 560                 AtomGroupIndices iatomsBuf;
 561
 562                 // Loop over all interactions of ftype
 563                 for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements)
 564                 {
 565                     // Calculate number of qm atoms in the interaction
 566                     int numQm = 0;
 567                     for (int k = 1; k <= NRAL(ftype); k++)
 568                     {
 569                         if (isQMAtom(molType->ilist[ftype].iatoms[j + k] + start))
 570                         {
 571                             numQm++;
 572                         }
 573                     }
 574
 575                     // If at least 3 atoms are QM then remove interaction
 576                     if (numQm >= 3)
 577                     {
 578                         topInfo_.numDihedralsRemoved++;
 579                     }
 580                     else
 581                     {
 582                         // If several (>1) atoms is not QM then preserve interaction in the iatomsBuf
 583                         for (int k = 0; k < numInteractionElements; k++)
 584                         {
 585                             iatomsBuf.push_back(molType->ilist[ftype].iatoms[k + j]);
 586                         }
 587                     }
 588                 }
 589
 590                 // Swap iatomsBuf and molType->ilist[ftype].iatoms
 591                 molType->ilist[ftype].iatoms.swap(iatomsBuf);
 592             }
 593         }
 594     }
 595 }
 596
 597 void QMMMTopologyPreprocessor::modifyQMMMVirtualSites(gmx_mtop_t* mtop)
 598 {
 599     // Loop over all blocks in topology
 600     // molBlockIndex - index of current block in mtop
 601     for (size_t molBlockIndex = 0; molBlockIndex < mtop->molblock.size(); molBlockIndex++)
 602     {
 603         // check if current block contains QM atoms
 604         if (bQMBlock_[molBlockIndex])
 605         {
 606             // molType - strucutre with current block type
 607             gmx_moltype_t* molType = &mtop->moltype[mtop->molblock[molBlockIndex].type];
 608
 609             // start - global index of the first atom in the current block
 610             int start = mtop->moleculeBlockIndices[molBlockIndex].globalAtomStart;
 611
 612             // loop over all interaction types
 613             for (int ftype = 0; ftype < F_NRE; ftype++)
 614             {
 615                 // If this is VSite interaction and ilist is not empty
 616                 if (IS_VSITE(ftype) && !molType->ilist[ftype].empty())
 617                 {
 618                     // Number of elements in the iatoms array for the current ftype
 619                     const int numInteractionElements = NRAL(ftype) + 1;
 620
 621                     // Loop over all interactions of ftype
 622                     for (int j = 0; j < molType->ilist[ftype].size(); j += numInteractionElements)
 623                     {
 624                         // Calculate number of qm atoms in the interaction
 625                         int numQm = 0;
 626                         // Here k starts from 2 because first atom in the interaction is an actuall vsite index
 627                         for (int k = 2; k <= NRAL(ftype); k++)
 628                         {
 629                             if (isQMAtom(molType->ilist[ftype].iatoms[j + k] + start))
 630                             {
 631                                 numQm++;
 632                             }
 633                         }
 634
 635                         // If all atoms froming that virtual site are QM atoms
 636                         // then remove classical charge from that virtual site
 637                         if (numQm == (NRAL(ftype) - 1))
 638                         {
 639                             topInfo_.numVirtualSitesModified++;
 640                             topInfo_.totalClassicalChargeOfQMAtoms +=
 641                                     molType->atoms.atom[molType->ilist[ftype].iatoms[j + 1]].q;
 642                             molType->atoms.atom[molType->ilist[ftype].iatoms[j + 1]].q  = 0.0;
 643                             molType->atoms.atom[molType->ilist[ftype].iatoms[j + 1]].qB = 0.0;
 644                         }
 645                     }
 646                 }
 647             }
 648         }
 649     }
 650 }
 651
 652 } // namespace gmx