2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2008,2009,2010,
5 * Copyright (c) 2012,2013,2014,2015,2016,2017,2018,2019, by the GROMACS development team, led by
6 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
7 * and including many others, as listed in the AUTHORS file in the
8 * top-level source directory and at http://www.gromacs.org.
10 * GROMACS is free software; you can redistribute it and/or
11 * modify it under the terms of the GNU Lesser General Public License
12 * as published by the Free Software Foundation; either version 2.1
13 * of the License, or (at your option) any later version.
15 * GROMACS is distributed in the hope that it will be useful,
16 * but WITHOUT ANY WARRANTY; without even the implied warranty of
17 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
18 * Lesser General Public License for more details.
20 * You should have received a copy of the GNU Lesser General Public
21 * License along with GROMACS; if not, see
22 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
23 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
25 * If you want to redistribute modifications to GROMACS, please
26 * consider that scientific software is very special. Version
27 * control is crucial - bugs must be traceable. We will be happy to
28 * consider code for inclusion in the official distribution, but
29 * derived work must not be called official GROMACS. Details are found
30 * in the README & COPYING files - if they are missing, get the
31 * official version at http://www.gromacs.org.
33 * To help us fund GROMACS development, we humbly ask that you cite
34 * the research papers on the package. Check out http://www.gromacs.org.
38 #include "mtop_util.h"
46 #include "gromacs/math/vectypes.h"
47 #include "gromacs/topology/atoms.h"
48 #include "gromacs/topology/block.h"
49 #include "gromacs/topology/exclusionblocks.h"
50 #include "gromacs/topology/idef.h"
51 #include "gromacs/topology/ifunc.h"
52 #include "gromacs/topology/topology.h"
53 #include "gromacs/topology/topsort.h"
54 #include "gromacs/utility/arrayref.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/real.h"
57 #include "gromacs/utility/smalloc.h"
59 static int gmx_mtop_maxresnr(const gmx_mtop_t* mtop, int maxres_renum)
63 for (const gmx_moltype_t& moltype : mtop->moltype)
65 const t_atoms& atoms = moltype.atoms;
66 if (atoms.nres > maxres_renum)
68 for (int r = 0; r < atoms.nres; r++)
70 if (atoms.resinfo[r].nr > maxresnr)
72 maxresnr = atoms.resinfo[r].nr;
81 static void buildMolblockIndices(gmx_mtop_t* mtop)
83 mtop->moleculeBlockIndices.resize(mtop->molblock.size());
87 int residueNumberStart = mtop->maxresnr + 1;
88 int moleculeIndexStart = 0;
89 for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
91 const gmx_molblock_t& molb = mtop->molblock[mb];
92 MoleculeBlockIndices& indices = mtop->moleculeBlockIndices[mb];
93 const int numResPerMol = mtop->moltype[molb.type].atoms.nres;
95 indices.numAtomsPerMolecule = mtop->moltype[molb.type].atoms.nr;
96 indices.globalAtomStart = atomIndex;
97 indices.globalResidueStart = residueIndex;
98 atomIndex += molb.nmol * indices.numAtomsPerMolecule;
99 residueIndex += molb.nmol * numResPerMol;
100 indices.globalAtomEnd = atomIndex;
101 indices.residueNumberStart = residueNumberStart;
102 if (numResPerMol <= mtop->maxres_renum)
104 residueNumberStart += molb.nmol * numResPerMol;
106 indices.moleculeIndexStart = moleculeIndexStart;
107 moleculeIndexStart += molb.nmol;
111 void gmx_mtop_finalize(gmx_mtop_t* mtop)
115 if (mtop->molblock.size() == 1 && mtop->molblock[0].nmol == 1)
117 /* We have a single molecule only, no renumbering needed.
118 * This case also covers an mtop converted from pdb/gro/... input,
119 * so we retain the original residue numbering.
121 mtop->maxres_renum = 0;
125 /* We only renumber single residue molecules. Their intra-molecular
126 * residue numbering is anyhow irrelevant.
128 mtop->maxres_renum = 1;
131 env = getenv("GMX_MAXRESRENUM");
134 sscanf(env, "%d", &mtop->maxres_renum);
136 if (mtop->maxres_renum == -1)
138 /* -1 signals renumber residues in all molecules */
139 mtop->maxres_renum = INT_MAX;
142 mtop->maxresnr = gmx_mtop_maxresnr(mtop, mtop->maxres_renum);
144 buildMolblockIndices(mtop);
147 void gmx_mtop_count_atomtypes(const gmx_mtop_t* mtop, int state, int typecount[])
149 for (int i = 0; i < mtop->ffparams.atnr; ++i)
153 for (const gmx_molblock_t& molb : mtop->molblock)
155 const t_atoms& atoms = mtop->moltype[molb.type].atoms;
156 for (int i = 0; i < atoms.nr; ++i)
161 tpi = atoms.atom[i].type;
165 tpi = atoms.atom[i].typeB;
167 typecount[tpi] += molb.nmol;
172 int gmx_mtop_num_molecules(const gmx_mtop_t& mtop)
174 int numMolecules = 0;
175 for (const gmx_molblock_t& molb : mtop.molblock)
177 numMolecules += molb.nmol;
182 int gmx_mtop_nres(const gmx_mtop_t* mtop)
185 for (const gmx_molblock_t& molb : mtop->molblock)
187 nres += molb.nmol * mtop->moltype[molb.type].atoms.nres;
192 AtomIterator::AtomIterator(const gmx_mtop_t& mtop, int globalAtomNumber) :
195 atoms_(&mtop.moltype[mtop.molblock[0].type].atoms),
197 highestResidueNumber_(mtop.maxresnr),
199 globalAtomNumber_(globalAtomNumber)
201 GMX_ASSERT(globalAtomNumber == 0 || globalAtomNumber == mtop.natoms,
202 "Starting at other atoms not implemented yet");
205 AtomIterator& AtomIterator::operator++()
210 if (localAtomNumber_ >= atoms_->nr)
212 if (atoms_->nres <= mtop_->maxresnr)
214 /* Single residue molecule, increase the count with one */
215 highestResidueNumber_ += atoms_->nres;
218 localAtomNumber_ = 0;
219 if (currentMolecule_ >= mtop_->molblock[mblock_].nmol)
222 if (mblock_ >= mtop_->molblock.size())
226 atoms_ = &mtop_->moltype[mtop_->molblock[mblock_].type].atoms;
227 currentMolecule_ = 0;
233 AtomIterator AtomIterator::operator++(int)
235 AtomIterator temp = *this;
240 bool AtomIterator::operator==(const AtomIterator& o) const
242 return mtop_ == o.mtop_ && globalAtomNumber_ == o.globalAtomNumber_;
245 bool AtomIterator::operator!=(const AtomIterator& o) const
247 return !(*this == o);
250 const t_atom& AtomProxy::atom() const
252 return it_->atoms_->atom[it_->localAtomNumber_];
255 int AtomProxy::globalAtomNumber() const
257 return it_->globalAtomNumber_;
260 const char* AtomProxy::atomName() const
262 return *(it_->atoms_->atomname[it_->localAtomNumber_]);
265 const char* AtomProxy::residueName() const
267 int residueIndexInMolecule = it_->atoms_->atom[it_->localAtomNumber_].resind;
268 return *(it_->atoms_->resinfo[residueIndexInMolecule].name);
271 int AtomProxy::residueNumber() const
273 int residueIndexInMolecule = it_->atoms_->atom[it_->localAtomNumber_].resind;
274 if (it_->atoms_->nres <= it_->mtop_->maxres_renum)
276 return it_->highestResidueNumber_ + 1 + residueIndexInMolecule;
280 return it_->atoms_->resinfo[residueIndexInMolecule].nr;
284 const gmx_moltype_t& AtomProxy::moleculeType() const
286 return it_->mtop_->moltype[it_->mtop_->molblock[it_->mblock_].type];
289 int AtomProxy::atomNumberInMol() const
291 return it_->localAtomNumber_;
294 typedef struct gmx_mtop_atomloop_block
296 const gmx_mtop_t* mtop;
298 const t_atoms* atoms;
300 } t_gmx_mtop_atomloop_block;
302 gmx_mtop_atomloop_block_t gmx_mtop_atomloop_block_init(const gmx_mtop_t* mtop)
304 struct gmx_mtop_atomloop_block* aloop;
310 aloop->atoms = &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
311 aloop->at_local = -1;
316 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop)
321 gmx_bool gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop, const t_atom** atom, int* nmol)
323 if (aloop == nullptr)
325 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
330 if (aloop->at_local >= aloop->atoms->nr)
333 if (aloop->mblock >= aloop->mtop->molblock.size())
335 gmx_mtop_atomloop_block_destroy(aloop);
338 aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
342 *atom = &aloop->atoms->atom[aloop->at_local];
343 *nmol = aloop->mtop->molblock[aloop->mblock].nmol;
348 typedef struct gmx_mtop_ilistloop
350 const gmx_mtop_t* mtop;
354 gmx_mtop_ilistloop_t gmx_mtop_ilistloop_init(const gmx_mtop_t* mtop)
356 struct gmx_mtop_ilistloop* iloop;
366 gmx_mtop_ilistloop_t gmx_mtop_ilistloop_init(const gmx_mtop_t& mtop)
368 return gmx_mtop_ilistloop_init(&mtop);
371 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
376 const InteractionLists* gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop, int* nmol)
378 if (iloop == nullptr)
380 gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
384 if (iloop->mblock >= gmx::ssize(iloop->mtop->molblock))
386 if (iloop->mblock == gmx::ssize(iloop->mtop->molblock) && iloop->mtop->bIntermolecularInteractions)
389 return iloop->mtop->intermolecular_ilist.get();
392 gmx_mtop_ilistloop_destroy(iloop);
396 *nmol = iloop->mtop->molblock[iloop->mblock].nmol;
398 return &iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
400 typedef struct gmx_mtop_ilistloop_all
402 const gmx_mtop_t* mtop;
406 } t_gmx_mtop_ilist_all;
408 gmx_mtop_ilistloop_all_t gmx_mtop_ilistloop_all_init(const gmx_mtop_t* mtop)
410 struct gmx_mtop_ilistloop_all* iloop;
422 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop)
427 const InteractionLists* gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop, int* atnr_offset)
430 if (iloop == nullptr)
433 "gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
438 iloop->a_offset += iloop->mtop->moleculeBlockIndices[iloop->mblock].numAtomsPerMolecule;
443 /* Inter-molecular interactions, if present, are indexed with
444 * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
445 * check for this value in this conditional.
447 if (iloop->mblock == iloop->mtop->molblock.size()
448 || iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
452 if (iloop->mblock >= iloop->mtop->molblock.size())
454 if (iloop->mblock == iloop->mtop->molblock.size() && iloop->mtop->bIntermolecularInteractions)
457 return iloop->mtop->intermolecular_ilist.get();
460 gmx_mtop_ilistloop_all_destroy(iloop);
465 *atnr_offset = iloop->a_offset;
467 return &iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
470 int gmx_mtop_ftype_count(const gmx_mtop_t* mtop, int ftype)
472 gmx_mtop_ilistloop_t iloop;
477 iloop = gmx_mtop_ilistloop_init(mtop);
478 while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol))
480 n += nmol * (*il)[ftype].size() / (1 + NRAL(ftype));
483 if (mtop->bIntermolecularInteractions)
485 n += (*mtop->intermolecular_ilist)[ftype].size() / (1 + NRAL(ftype));
491 int gmx_mtop_ftype_count(const gmx_mtop_t& mtop, int ftype)
493 return gmx_mtop_ftype_count(&mtop, ftype);
496 int gmx_mtop_interaction_count(const gmx_mtop_t& mtop, const int unsigned if_flags)
500 gmx_mtop_ilistloop_t iloop = gmx_mtop_ilistloop_init(mtop);
502 while (const InteractionLists* il = gmx_mtop_ilistloop_next(iloop, &nmol))
504 for (int ftype = 0; ftype < F_NRE; ftype++)
506 if ((interaction_function[ftype].flags & if_flags) == if_flags)
508 n += nmol * (*il)[ftype].size() / (1 + NRAL(ftype));
513 if (mtop.bIntermolecularInteractions)
515 for (int ftype = 0; ftype < F_NRE; ftype++)
517 if ((interaction_function[ftype].flags & if_flags) == if_flags)
519 n += (*mtop.intermolecular_ilist)[ftype].size() / (1 + NRAL(ftype));
527 static void atomcat(t_atoms* dest, const t_atoms* src, int copies, int maxres_renum, int* maxresnr)
531 int destnr = dest->nr;
535 dest->haveMass = src->haveMass;
536 dest->haveType = src->haveType;
537 dest->haveCharge = src->haveCharge;
538 dest->haveBState = src->haveBState;
539 dest->havePdbInfo = src->havePdbInfo;
543 dest->haveMass = dest->haveMass && src->haveMass;
544 dest->haveType = dest->haveType && src->haveType;
545 dest->haveCharge = dest->haveCharge && src->haveCharge;
546 dest->haveBState = dest->haveBState && src->haveBState;
547 dest->havePdbInfo = dest->havePdbInfo && src->havePdbInfo;
552 size = destnr + copies * srcnr;
553 srenew(dest->atom, size);
554 srenew(dest->atomname, size);
557 srenew(dest->atomtype, size);
558 if (dest->haveBState)
560 srenew(dest->atomtypeB, size);
563 if (dest->havePdbInfo)
565 srenew(dest->pdbinfo, size);
570 size = dest->nres + copies * src->nres;
571 srenew(dest->resinfo, size);
574 /* residue information */
575 for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
577 memcpy(reinterpret_cast<char*>(&(dest->resinfo[l])), reinterpret_cast<char*>(&(src->resinfo[0])),
578 static_cast<size_t>(src->nres * sizeof(src->resinfo[0])));
581 for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
583 memcpy(reinterpret_cast<char*>(&(dest->atom[l])), reinterpret_cast<char*>(&(src->atom[0])),
584 static_cast<size_t>(srcnr * sizeof(src->atom[0])));
585 memcpy(reinterpret_cast<char*>(&(dest->atomname[l])),
586 reinterpret_cast<char*>(&(src->atomname[0])),
587 static_cast<size_t>(srcnr * sizeof(src->atomname[0])));
590 memcpy(reinterpret_cast<char*>(&(dest->atomtype[l])),
591 reinterpret_cast<char*>(&(src->atomtype[0])),
592 static_cast<size_t>(srcnr * sizeof(src->atomtype[0])));
593 if (dest->haveBState)
595 memcpy(reinterpret_cast<char*>(&(dest->atomtypeB[l])),
596 reinterpret_cast<char*>(&(src->atomtypeB[0])),
597 static_cast<size_t>(srcnr * sizeof(src->atomtypeB[0])));
600 if (dest->havePdbInfo)
602 memcpy(reinterpret_cast<char*>(&(dest->pdbinfo[l])),
603 reinterpret_cast<char*>(&(src->pdbinfo[0])),
604 static_cast<size_t>(srcnr * sizeof(src->pdbinfo[0])));
608 /* Increment residue indices */
609 for (l = destnr, j = 0; (j < copies); j++)
611 for (i = 0; (i < srcnr); i++, l++)
613 dest->atom[l].resind = dest->nres + j * src->nres + src->atom[i].resind;
617 if (src->nres <= maxres_renum)
619 /* Single residue molecule, continue counting residues */
620 for (j = 0; (j < copies); j++)
622 for (l = 0; l < src->nres; l++)
625 dest->resinfo[dest->nres + j * src->nres + l].nr = *maxresnr;
630 dest->nres += copies * src->nres;
631 dest->nr += copies * src->nr;
634 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t* mtop)
638 init_t_atoms(&atoms, 0, FALSE);
640 int maxresnr = mtop->maxresnr;
641 for (const gmx_molblock_t& molb : mtop->molblock)
643 atomcat(&atoms, &mtop->moltype[molb.type].atoms, molb.nmol, mtop->maxres_renum, &maxresnr);
650 * The cat routines below are old code from src/kernel/topcat.c
653 static void blockacat(t_blocka* dest, const t_blocka* src, int copies, int dnum, int snum)
656 int destnr = dest->nr;
657 int destnra = dest->nra;
661 size = (dest->nr + copies * src->nr + 1);
662 srenew(dest->index, size);
666 size = (dest->nra + copies * src->nra);
667 srenew(dest->a, size);
670 for (l = destnr, j = 0; (j < copies); j++)
672 for (i = 0; (i < src->nr); i++)
674 dest->index[l++] = dest->nra + src->index[i];
676 dest->nra += src->nra;
678 for (l = destnra, j = 0; (j < copies); j++)
680 for (i = 0; (i < src->nra); i++)
682 dest->a[l++] = dnum + src->a[i];
687 dest->index[dest->nr] = dest->nra;
688 dest->nalloc_index = dest->nr;
689 dest->nalloc_a = dest->nra;
692 static void ilistcat(int ftype, t_ilist* dest, const InteractionList& src, int copies, int dnum, int snum)
698 dest->nalloc = dest->nr + copies * src.size();
699 srenew(dest->iatoms, dest->nalloc);
701 for (c = 0; c < copies; c++)
703 for (i = 0; i < src.size();)
705 dest->iatoms[dest->nr++] = src.iatoms[i++];
706 for (a = 0; a < nral; a++)
708 dest->iatoms[dest->nr++] = dnum + src.iatoms[i++];
715 static void set_posres_params(t_idef* idef, const gmx_molblock_t* molb, int i0, int a_offset)
721 il = &idef->il[F_POSRES];
723 idef->iparams_posres_nalloc = i1;
724 srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
725 for (i = i0; i < i1; i++)
727 ip = &idef->iparams_posres[i];
728 /* Copy the force constants */
729 *ip = idef->iparams[il->iatoms[i * 2]];
730 a_molb = il->iatoms[i * 2 + 1] - a_offset;
731 if (molb->posres_xA.empty())
733 gmx_incons("Position restraint coordinates are missing");
735 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
736 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
737 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
738 if (!molb->posres_xB.empty())
740 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
741 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
742 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
746 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
747 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
748 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
750 /* Set the parameter index for idef->iparams_posre */
751 il->iatoms[i * 2] = i;
755 static void set_fbposres_params(t_idef* idef, const gmx_molblock_t* molb, int i0, int a_offset)
761 il = &idef->il[F_FBPOSRES];
763 idef->iparams_fbposres_nalloc = i1;
764 srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
765 for (i = i0; i < i1; i++)
767 ip = &idef->iparams_fbposres[i];
768 /* Copy the force constants */
769 *ip = idef->iparams[il->iatoms[i * 2]];
770 a_molb = il->iatoms[i * 2 + 1] - a_offset;
771 if (molb->posres_xA.empty())
773 gmx_incons("Position restraint coordinates are missing");
775 /* Take flat-bottom posres reference from normal position restraints */
776 ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
777 ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
778 ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
779 /* Note: no B-type for flat-bottom posres */
781 /* Set the parameter index for idef->iparams_posre */
782 il->iatoms[i * 2] = i;
786 /*! \brief Copy idef structure from mtop.
788 * Makes a deep copy of an idef data structure from a gmx_mtop_t.
789 * Used to initialize legacy topology types.
791 * \param[in] mtop Reference to input mtop.
792 * \param[in] idef Pointer to idef to populate.
793 * \param[in] mergeConstr Decide if constraints will be merged.
794 * \param[in] freeEnergyInteractionsAtEnd Decide if free energy stuff should
795 * be added at the end.
797 static void copyIdefFromMtop(const gmx_mtop_t& mtop, t_idef* idef, bool freeEnergyInteractionsAtEnd, bool mergeConstr)
799 const gmx_ffparams_t* ffp = &mtop.ffparams;
801 idef->ntypes = ffp->numTypes();
802 idef->atnr = ffp->atnr;
803 /* we can no longer copy the pointers to the mtop members,
804 * because they will become invalid as soon as mtop gets free'd.
805 * We also need to make sure to only operate on valid data!
808 if (!ffp->functype.empty())
810 snew(idef->functype, ffp->functype.size());
811 std::copy(ffp->functype.data(), ffp->functype.data() + ffp->functype.size(), idef->functype);
815 idef->functype = nullptr;
817 if (!ffp->iparams.empty())
819 snew(idef->iparams, ffp->iparams.size());
820 std::copy(ffp->iparams.data(), ffp->iparams.data() + ffp->iparams.size(), idef->iparams);
824 idef->iparams = nullptr;
826 idef->iparams_posres = nullptr;
827 idef->iparams_posres_nalloc = 0;
828 idef->iparams_fbposres = nullptr;
829 idef->iparams_fbposres_nalloc = 0;
830 idef->fudgeQQ = ffp->fudgeQQ;
831 idef->cmap_grid = new gmx_cmap_t;
832 *idef->cmap_grid = ffp->cmap_grid;
833 idef->ilsort = ilsortUNKNOWN;
835 for (int ftype = 0; ftype < F_NRE; ftype++)
837 idef->il[ftype].nr = 0;
838 idef->il[ftype].nalloc = 0;
839 idef->il[ftype].iatoms = nullptr;
843 for (const gmx_molblock_t& molb : mtop.molblock)
845 const gmx_moltype_t& molt = mtop.moltype[molb.type];
847 int srcnr = molt.atoms.nr;
850 int nposre_old = idef->il[F_POSRES].nr;
851 int nfbposre_old = idef->il[F_FBPOSRES].nr;
852 for (int ftype = 0; ftype < F_NRE; ftype++)
854 if (mergeConstr && ftype == F_CONSTR && molt.ilist[F_CONSTRNC].size() > 0)
856 /* Merge all constrains into one ilist.
857 * This simplifies the constraint code.
859 for (int mol = 0; mol < molb.nmol; mol++)
861 ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTR], 1,
862 destnr + mol * srcnr, srcnr);
863 ilistcat(ftype, &idef->il[F_CONSTR], molt.ilist[F_CONSTRNC], 1,
864 destnr + mol * srcnr, srcnr);
867 else if (!(mergeConstr && ftype == F_CONSTRNC))
869 ilistcat(ftype, &idef->il[ftype], molt.ilist[ftype], molb.nmol, destnr, srcnr);
872 if (idef->il[F_POSRES].nr > nposre_old)
874 /* Executing this line line stops gmxdump -sys working
875 * correctly. I'm not aware there's an elegant fix. */
876 set_posres_params(idef, &molb, nposre_old / 2, natoms);
878 if (idef->il[F_FBPOSRES].nr > nfbposre_old)
880 set_fbposres_params(idef, &molb, nfbposre_old / 2, natoms);
883 natoms += molb.nmol * srcnr;
886 if (mtop.bIntermolecularInteractions)
888 for (int ftype = 0; ftype < F_NRE; ftype++)
890 ilistcat(ftype, &idef->il[ftype], (*mtop.intermolecular_ilist)[ftype], 1, 0, mtop.natoms);
894 if (freeEnergyInteractionsAtEnd && gmx_mtop_bondeds_free_energy(&mtop))
896 std::vector<real> qA(mtop.natoms);
897 std::vector<real> qB(mtop.natoms);
898 for (const AtomProxy atomP : AtomRange(mtop))
900 const t_atom& local = atomP.atom();
901 int index = atomP.globalAtomNumber();
903 qB[index] = local.qB;
905 gmx_sort_ilist_fe(idef, qA.data(), qB.data());
909 idef->ilsort = ilsortNO_FE;
913 /*! \brief Copy atomtypes from mtop
915 * Makes a deep copy of t_atomtypes from gmx_mtop_t.
916 * Used to initialize legacy topology types.
918 * \param[in] mtop Reference to input mtop.
919 * \param[in] atomtypes Pointer to atomtypes to populate.
921 static void copyAtomtypesFromMtop(const gmx_mtop_t& mtop, t_atomtypes* atomtypes)
923 atomtypes->nr = mtop.atomtypes.nr;
924 if (mtop.atomtypes.atomnumber)
926 snew(atomtypes->atomnumber, mtop.atomtypes.nr);
927 std::copy(mtop.atomtypes.atomnumber, mtop.atomtypes.atomnumber + mtop.atomtypes.nr,
928 atomtypes->atomnumber);
932 atomtypes->atomnumber = nullptr;
936 /*! \brief Copy excls from mtop.
938 * Makes a deep copy of excls(t_blocka) from gmx_mtop_t.
939 * Used to initialize legacy topology types.
941 * \param[in] mtop Reference to input mtop.
942 * \param[in] excls Pointer to final excls data structure.
944 static void copyExclsFromMtop(const gmx_mtop_t& mtop, t_blocka* excls)
948 for (const gmx_molblock_t& molb : mtop.molblock)
950 const gmx_moltype_t& molt = mtop.moltype[molb.type];
952 int srcnr = molt.atoms.nr;
955 blockacat(excls, &molt.excls, molb.nmol, destnr, srcnr);
957 natoms += molb.nmol * srcnr;
961 /*! \brief Updates inter-molecular exclusion lists
963 * This function updates inter-molecular exclusions to exclude all
964 * non-bonded interactions between a given list of atoms
966 * \param[inout] excls existing exclusions in local topology
967 * \param[in] ids list of global IDs of atoms
969 static void addMimicExclusions(t_blocka* excls, const gmx::ArrayRef<const int> ids)
971 t_blocka inter_excl{};
972 init_blocka(&inter_excl);
973 size_t n_q = ids.size();
975 inter_excl.nr = excls->nr;
976 inter_excl.nra = n_q * n_q;
978 size_t total_nra = n_q * n_q;
980 snew(inter_excl.index, excls->nr + 1);
981 snew(inter_excl.a, total_nra);
983 for (int i = 0; i < excls->nr; ++i)
985 inter_excl.index[i] = 0;
988 /* Here we loop over the list of QM atom ids
989 * and create exclusions between all of them resulting in
990 * n_q * n_q sized exclusion list
993 for (int k = 0; k < inter_excl.nr; ++k)
995 inter_excl.index[k] = prev_index;
996 for (long i = 0; i < ids.ssize(); i++)
1002 size_t index = n_q * i;
1003 inter_excl.index[ids[i]] = index;
1004 prev_index = index + n_q;
1005 for (size_t j = 0; j < n_q; ++j)
1007 inter_excl.a[n_q * i + j] = ids[j];
1011 inter_excl.index[ids[n_q - 1] + 1] = n_q * n_q;
1013 inter_excl.index[inter_excl.nr] = n_q * n_q;
1015 std::vector<gmx::ExclusionBlock> qmexcl2(excls->nr);
1016 gmx::blockaToExclusionBlocks(&inter_excl, qmexcl2);
1018 // Merge the created exclusion list with the existing one
1019 gmx::mergeExclusions(excls, qmexcl2);
1022 static void gen_local_top(const gmx_mtop_t& mtop,
1023 bool freeEnergyInteractionsAtEnd,
1025 gmx_localtop_t* top)
1027 copyAtomtypesFromMtop(mtop, &top->atomtypes);
1028 copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
1029 copyExclsFromMtop(mtop, &top->excls);
1030 if (!mtop.intermolecularExclusionGroup.empty())
1032 addMimicExclusions(&top->excls, mtop.intermolecularExclusionGroup);
1036 void gmx_mtop_generate_local_top(const gmx_mtop_t& mtop, gmx_localtop_t* top, bool freeEnergyInteractionsAtEnd)
1038 gen_local_top(mtop, freeEnergyInteractionsAtEnd, true, top);
1041 /*! \brief Fills an array with molecule begin/end atom indices
1043 * \param[in] mtop The global topology
1044 * \param[out] index Array of size nr. of molecules + 1 to be filled with molecule begin/end indices
1046 static void fillMoleculeIndices(const gmx_mtop_t& mtop, gmx::ArrayRef<int> index)
1048 int globalAtomIndex = 0;
1049 int globalMolIndex = 0;
1050 index[globalMolIndex] = globalAtomIndex;
1051 for (const gmx_molblock_t& molb : mtop.molblock)
1053 int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
1054 for (int mol = 0; mol < molb.nmol; mol++)
1056 globalAtomIndex += numAtomsPerMolecule;
1057 globalMolIndex += 1;
1058 index[globalMolIndex] = globalAtomIndex;
1063 gmx::RangePartitioning gmx_mtop_molecules(const gmx_mtop_t& mtop)
1065 gmx::RangePartitioning mols;
1067 for (const gmx_molblock_t& molb : mtop.molblock)
1069 int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
1070 for (int mol = 0; mol < molb.nmol; mol++)
1072 mols.appendBlock(numAtomsPerMolecule);
1079 /*! \brief Creates and returns a deprecated t_block struct with molecule indices
1081 * \param[in] mtop The global topology
1083 static t_block gmx_mtop_molecules_t_block(const gmx_mtop_t& mtop)
1087 mols.nr = gmx_mtop_num_molecules(mtop);
1088 mols.nalloc_index = mols.nr + 1;
1089 snew(mols.index, mols.nalloc_index);
1091 fillMoleculeIndices(mtop, gmx::arrayRefFromArray(mols.index, mols.nr + 1));
1096 static void gen_t_topology(const gmx_mtop_t& mtop,
1097 bool freeEnergyInteractionsAtEnd,
1101 copyAtomtypesFromMtop(mtop, &top->atomtypes);
1102 copyIdefFromMtop(mtop, &top->idef, freeEnergyInteractionsAtEnd, bMergeConstr);
1103 copyExclsFromMtop(mtop, &top->excls);
1105 top->name = mtop.name;
1106 top->atoms = gmx_mtop_global_atoms(&mtop);
1107 top->mols = gmx_mtop_molecules_t_block(mtop);
1108 top->bIntermolecularInteractions = mtop.bIntermolecularInteractions;
1109 top->symtab = mtop.symtab;
1112 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t* mtop, bool freeMTop)
1116 gen_t_topology(*mtop, false, false, &top);
1120 // Clear pointers and counts, such that the pointers copied to top
1121 // keep pointing to valid data after destroying mtop.
1122 mtop->symtab.symbuf = nullptr;
1123 mtop->symtab.nr = 0;
1128 std::vector<int> get_atom_index(const gmx_mtop_t* mtop)
1131 std::vector<int> atom_index;
1132 for (const AtomProxy atomP : AtomRange(*mtop))
1134 const t_atom& local = atomP.atom();
1135 int index = atomP.globalAtomNumber();
1136 if (local.ptype == eptAtom)
1138 atom_index.push_back(index);
1144 void convertAtomsToMtop(t_symtab* symtab, char** name, t_atoms* atoms, gmx_mtop_t* mtop)
1146 mtop->symtab = *symtab;
1150 mtop->moltype.clear();
1151 mtop->moltype.resize(1);
1152 mtop->moltype.back().atoms = *atoms;
1154 mtop->molblock.resize(1);
1155 mtop->molblock[0].type = 0;
1156 mtop->molblock[0].nmol = 1;
1158 mtop->bIntermolecularInteractions = FALSE;
1160 mtop->natoms = atoms->nr;
1162 mtop->haveMoleculeIndices = false;
1164 gmx_mtop_finalize(mtop);