2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2008,2009,2010,2012,2013,2014,2015,2016,2017,2018, 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.
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.
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.
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.
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.
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.
37 #include "mtop_util.h"
45 #include "gromacs/math/vectypes.h"
46 #include "gromacs/topology/atoms.h"
47 #include "gromacs/topology/block.h"
48 #include "gromacs/topology/idef.h"
49 #include "gromacs/topology/ifunc.h"
50 #include "gromacs/topology/topology.h"
51 #include "gromacs/topology/topsort.h"
52 #include "gromacs/utility/arrayref.h"
53 #include "gromacs/utility/fatalerror.h"
54 #include "gromacs/utility/real.h"
55 #include "gromacs/utility/smalloc.h"
57 static int gmx_mtop_maxresnr(const gmx_mtop_t *mtop, int maxres_renum)
61 for (const gmx_moltype_t &moltype : mtop->moltype)
63 const t_atoms &atoms = moltype.atoms;
64 if (atoms.nres > maxres_renum)
66 for (int r = 0; r < atoms.nres; r++)
68 if (atoms.resinfo[r].nr > maxresnr)
70 maxresnr = atoms.resinfo[r].nr;
79 static void buildMolblockIndices(gmx_mtop_t *mtop)
81 mtop->moleculeBlockIndices.resize(mtop->molblock.size());
85 int residueNumberStart = mtop->maxresnr + 1;
86 int moleculeIndexStart = 0;
87 for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
89 const gmx_molblock_t &molb = mtop->molblock[mb];
90 MoleculeBlockIndices &indices = mtop->moleculeBlockIndices[mb];
91 const int numResPerMol = mtop->moltype[molb.type].atoms.nres;
93 indices.numAtomsPerMolecule = mtop->moltype[molb.type].atoms.nr;
94 indices.globalAtomStart = atomIndex;
95 indices.globalResidueStart = residueIndex;
96 atomIndex += molb.nmol*indices.numAtomsPerMolecule;
97 residueIndex += molb.nmol*numResPerMol;
98 indices.globalAtomEnd = atomIndex;
99 indices.residueNumberStart = residueNumberStart;
100 if (numResPerMol <= mtop->maxres_renum)
102 residueNumberStart += molb.nmol*numResPerMol;
104 indices.moleculeIndexStart = moleculeIndexStart;
105 moleculeIndexStart += molb.nmol;
109 void gmx_mtop_finalize(gmx_mtop_t *mtop)
113 if (mtop->molblock.size() == 1 && mtop->molblock[0].nmol == 1)
115 /* We have a single molecule only, no renumbering needed.
116 * This case also covers an mtop converted from pdb/gro/... input,
117 * so we retain the original residue numbering.
119 mtop->maxres_renum = 0;
123 /* We only renumber single residue molecules. Their intra-molecular
124 * residue numbering is anyhow irrelevant.
126 mtop->maxres_renum = 1;
129 env = getenv("GMX_MAXRESRENUM");
132 sscanf(env, "%d", &mtop->maxres_renum);
134 if (mtop->maxres_renum == -1)
136 /* -1 signals renumber residues in all molecules */
137 mtop->maxres_renum = INT_MAX;
140 mtop->maxresnr = gmx_mtop_maxresnr(mtop, mtop->maxres_renum);
142 buildMolblockIndices(mtop);
145 void gmx_mtop_count_atomtypes(const gmx_mtop_t *mtop, int state, int typecount[])
147 for (int i = 0; i < mtop->ffparams.atnr; ++i)
151 for (const gmx_molblock_t &molb : mtop->molblock)
153 const t_atoms &atoms = mtop->moltype[molb.type].atoms;
154 for (int i = 0; i < atoms.nr; ++i)
159 tpi = atoms.atom[i].type;
163 tpi = atoms.atom[i].typeB;
165 typecount[tpi] += molb.nmol;
170 int ncg_mtop(const gmx_mtop_t *mtop)
173 for (const gmx_molblock_t &molb : mtop->molblock)
175 ncg += molb.nmol*mtop->moltype[molb.type].cgs.nr;
181 int gmx_mtop_num_molecules(const gmx_mtop_t &mtop)
183 int numMolecules = 0;
184 for (const gmx_molblock_t &molb : mtop.molblock)
186 numMolecules += molb.nmol;
191 int gmx_mtop_nres(const gmx_mtop_t *mtop)
194 for (const gmx_molblock_t &molb : mtop->molblock)
196 nres += molb.nmol*mtop->moltype[molb.type].atoms.nres;
201 void gmx_mtop_remove_chargegroups(gmx_mtop_t *mtop)
203 for (gmx_moltype_t &molt : mtop->moltype)
205 t_block &cgs = molt.cgs;
206 if (cgs.nr < molt.atoms.nr)
208 cgs.nr = molt.atoms.nr;
209 srenew(cgs.index, cgs.nr + 1);
210 for (int i = 0; i < cgs.nr + 1; i++)
218 typedef struct gmx_mtop_atomloop_all
220 const gmx_mtop_t *mtop;
222 const t_atoms *atoms;
227 } t_gmx_mtop_atomloop_all;
229 gmx_mtop_atomloop_all_t
230 gmx_mtop_atomloop_all_init(const gmx_mtop_t *mtop)
232 struct gmx_mtop_atomloop_all *aloop;
239 &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
241 aloop->maxresnr = mtop->maxresnr;
242 aloop->at_local = -1;
243 aloop->at_global = -1;
248 static void gmx_mtop_atomloop_all_destroy(gmx_mtop_atomloop_all_t aloop)
253 gmx_bool gmx_mtop_atomloop_all_next(gmx_mtop_atomloop_all_t aloop,
254 int *at_global, const t_atom **atom)
256 if (aloop == nullptr)
258 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
264 if (aloop->at_local >= aloop->atoms->nr)
266 if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
268 /* Single residue molecule, increase the count with one */
269 aloop->maxresnr += aloop->atoms->nres;
273 if (aloop->mol >= aloop->mtop->molblock[aloop->mblock].nmol)
276 if (aloop->mblock >= aloop->mtop->molblock.size())
278 gmx_mtop_atomloop_all_destroy(aloop);
281 aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
286 *at_global = aloop->at_global;
287 *atom = &aloop->atoms->atom[aloop->at_local];
292 void gmx_mtop_atomloop_all_names(gmx_mtop_atomloop_all_t aloop,
293 char **atomname, int *resnr, char **resname)
297 *atomname = *(aloop->atoms->atomname[aloop->at_local]);
298 resind_mol = aloop->atoms->atom[aloop->at_local].resind;
299 *resnr = aloop->atoms->resinfo[resind_mol].nr;
300 if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
302 *resnr = aloop->maxresnr + 1 + resind_mol;
304 *resname = *(aloop->atoms->resinfo[resind_mol].name);
307 void gmx_mtop_atomloop_all_moltype(gmx_mtop_atomloop_all_t aloop,
308 const gmx_moltype_t **moltype,
311 *moltype = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type];
312 *at_mol = aloop->at_local;
315 typedef struct gmx_mtop_atomloop_block
317 const gmx_mtop_t *mtop;
319 const t_atoms *atoms;
321 } t_gmx_mtop_atomloop_block;
323 gmx_mtop_atomloop_block_t
324 gmx_mtop_atomloop_block_init(const gmx_mtop_t *mtop)
326 struct gmx_mtop_atomloop_block *aloop;
332 aloop->atoms = &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
333 aloop->at_local = -1;
338 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop)
343 gmx_bool gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop,
344 const t_atom **atom, int *nmol)
346 if (aloop == nullptr)
348 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
353 if (aloop->at_local >= aloop->atoms->nr)
356 if (aloop->mblock >= aloop->mtop->molblock.size())
358 gmx_mtop_atomloop_block_destroy(aloop);
361 aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
365 *atom = &aloop->atoms->atom[aloop->at_local];
366 *nmol = aloop->mtop->molblock[aloop->mblock].nmol;
371 typedef struct gmx_mtop_ilistloop
373 const gmx_mtop_t *mtop;
378 gmx_mtop_ilistloop_init(const gmx_mtop_t *mtop)
380 struct gmx_mtop_ilistloop *iloop;
391 gmx_mtop_ilistloop_init(const gmx_mtop_t &mtop)
393 return gmx_mtop_ilistloop_init(&mtop);
396 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
401 gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
402 const t_ilist **ilist_mol, int *nmol)
404 if (iloop == nullptr)
406 gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
410 if (iloop->mblock >= static_cast<int>(iloop->mtop->molblock.size()))
412 if (iloop->mblock == static_cast<int>(iloop->mtop->molblock.size()) &&
413 iloop->mtop->bIntermolecularInteractions)
415 *ilist_mol = iloop->mtop->intermolecular_ilist;
420 gmx_mtop_ilistloop_destroy(iloop);
425 iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
427 *nmol = iloop->mtop->molblock[iloop->mblock].nmol;
431 typedef struct gmx_mtop_ilistloop_all
433 const gmx_mtop_t *mtop;
437 } t_gmx_mtop_ilist_all;
439 gmx_mtop_ilistloop_all_t
440 gmx_mtop_ilistloop_all_init(const gmx_mtop_t *mtop)
442 struct gmx_mtop_ilistloop_all *iloop;
454 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop)
459 gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
460 const t_ilist **ilist_mol, int *atnr_offset)
463 if (iloop == nullptr)
465 gmx_incons("gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
470 iloop->a_offset += iloop->mtop->moleculeBlockIndices[iloop->mblock].numAtomsPerMolecule;
475 /* Inter-molecular interactions, if present, are indexed with
476 * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
477 * check for this value in this conditional.
479 if (iloop->mblock == iloop->mtop->molblock.size() ||
480 iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
484 if (iloop->mblock >= iloop->mtop->molblock.size())
486 if (iloop->mblock == iloop->mtop->molblock.size() &&
487 iloop->mtop->bIntermolecularInteractions)
489 *ilist_mol = iloop->mtop->intermolecular_ilist;
494 gmx_mtop_ilistloop_all_destroy(iloop);
500 iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
502 *atnr_offset = iloop->a_offset;
507 int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype)
509 gmx_mtop_ilistloop_t iloop;
515 iloop = gmx_mtop_ilistloop_init(mtop);
516 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
518 n += nmol*il[ftype].nr/(1+NRAL(ftype));
521 if (mtop->bIntermolecularInteractions)
523 n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype));
529 int gmx_mtop_ftype_count(const gmx_mtop_t &mtop, int ftype)
531 return gmx_mtop_ftype_count(&mtop, ftype);
534 t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
537 /* In most cases this is too much, but we realloc at the end */
538 snew(cgs_gl.index, mtop->natoms+1);
542 for (const gmx_molblock_t &molb : mtop->molblock)
544 const t_block &cgs_mol = mtop->moltype[molb.type].cgs;
545 for (int mol = 0; mol < molb.nmol; mol++)
547 for (int cg = 0; cg < cgs_mol.nr; cg++)
549 cgs_gl.index[cgs_gl.nr + 1] =
550 cgs_gl.index[cgs_gl.nr] +
551 cgs_mol.index[cg + 1] - cgs_mol.index[cg];
556 cgs_gl.nalloc_index = cgs_gl.nr + 1;
557 srenew(cgs_gl.index, cgs_gl.nalloc_index);
562 static void atomcat(t_atoms *dest, const t_atoms *src, int copies,
563 int maxres_renum, int *maxresnr)
567 int destnr = dest->nr;
571 dest->haveMass = src->haveMass;
572 dest->haveType = src->haveType;
573 dest->haveCharge = src->haveCharge;
574 dest->haveBState = src->haveBState;
575 dest->havePdbInfo = src->havePdbInfo;
579 dest->haveMass = dest->haveMass && src->haveMass;
580 dest->haveType = dest->haveType && src->haveType;
581 dest->haveCharge = dest->haveCharge && src->haveCharge;
582 dest->haveBState = dest->haveBState && src->haveBState;
583 dest->havePdbInfo = dest->havePdbInfo && src->havePdbInfo;
588 size = destnr+copies*srcnr;
589 srenew(dest->atom, size);
590 srenew(dest->atomname, size);
593 srenew(dest->atomtype, size);
594 if (dest->haveBState)
596 srenew(dest->atomtypeB, size);
599 if (dest->havePdbInfo)
601 srenew(dest->pdbinfo, size);
606 size = dest->nres+copies*src->nres;
607 srenew(dest->resinfo, size);
610 /* residue information */
611 for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
613 memcpy((char *) &(dest->resinfo[l]), (char *) &(src->resinfo[0]),
614 (size_t)(src->nres*sizeof(src->resinfo[0])));
617 for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
619 memcpy((char *) &(dest->atom[l]), (char *) &(src->atom[0]),
620 (size_t)(srcnr*sizeof(src->atom[0])));
621 memcpy((char *) &(dest->atomname[l]), (char *) &(src->atomname[0]),
622 (size_t)(srcnr*sizeof(src->atomname[0])));
625 memcpy((char *) &(dest->atomtype[l]), (char *) &(src->atomtype[0]),
626 (size_t)(srcnr*sizeof(src->atomtype[0])));
627 if (dest->haveBState)
629 memcpy((char *) &(dest->atomtypeB[l]), (char *) &(src->atomtypeB[0]),
630 (size_t)(srcnr*sizeof(src->atomtypeB[0])));
633 if (dest->havePdbInfo)
635 memcpy((char *) &(dest->pdbinfo[l]), (char *) &(src->pdbinfo[0]),
636 (size_t)(srcnr*sizeof(src->pdbinfo[0])));
640 /* Increment residue indices */
641 for (l = destnr, j = 0; (j < copies); j++)
643 for (i = 0; (i < srcnr); i++, l++)
645 dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;
649 if (src->nres <= maxres_renum)
651 /* Single residue molecule, continue counting residues */
652 for (j = 0; (j < copies); j++)
654 for (l = 0; l < src->nres; l++)
657 dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;
662 dest->nres += copies*src->nres;
663 dest->nr += copies*src->nr;
666 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
670 init_t_atoms(&atoms, 0, FALSE);
672 int maxresnr = mtop->maxresnr;
673 for (const gmx_molblock_t &molb : mtop->molblock)
675 atomcat(&atoms, &mtop->moltype[molb.type].atoms, molb.nmol,
676 mtop->maxres_renum, &maxresnr);
683 * The cat routines below are old code from src/kernel/topcat.c
686 static void blockcat(t_block *dest, const t_block *src, int copies)
688 int i, j, l, nra, size;
692 size = (dest->nr+copies*src->nr+1);
693 srenew(dest->index, size);
696 nra = dest->index[dest->nr];
697 for (l = dest->nr, j = 0; (j < copies); j++)
699 for (i = 0; (i < src->nr); i++)
701 dest->index[l++] = nra + src->index[i];
705 nra += src->index[src->nr];
708 dest->nr += copies*src->nr;
709 dest->index[dest->nr] = nra;
712 static void blockacat(t_blocka *dest, const t_blocka *src, int copies,
716 int destnr = dest->nr;
717 int destnra = dest->nra;
721 size = (dest->nr+copies*src->nr+1);
722 srenew(dest->index, size);
726 size = (dest->nra+copies*src->nra);
727 srenew(dest->a, size);
730 for (l = destnr, j = 0; (j < copies); j++)
732 for (i = 0; (i < src->nr); i++)
734 dest->index[l++] = dest->nra+src->index[i];
736 dest->nra += src->nra;
738 for (l = destnra, j = 0; (j < copies); j++)
740 for (i = 0; (i < src->nra); i++)
742 dest->a[l++] = dnum+src->a[i];
747 dest->index[dest->nr] = dest->nra;
750 static void ilistcat(int ftype, t_ilist *dest, const t_ilist *src, int copies,
757 dest->nalloc = dest->nr + copies*src->nr;
758 srenew(dest->iatoms, dest->nalloc);
760 for (c = 0; c < copies; c++)
762 for (i = 0; i < src->nr; )
764 dest->iatoms[dest->nr++] = src->iatoms[i++];
765 for (a = 0; a < nral; a++)
767 dest->iatoms[dest->nr++] = dnum + src->iatoms[i++];
774 static void set_posres_params(t_idef *idef, const gmx_molblock_t *molb,
775 int i0, int a_offset)
781 il = &idef->il[F_POSRES];
783 idef->iparams_posres_nalloc = i1;
784 srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
785 for (i = i0; i < i1; i++)
787 ip = &idef->iparams_posres[i];
788 /* Copy the force constants */
789 *ip = idef->iparams[il->iatoms[i*2]];
790 a_molb = il->iatoms[i*2+1] - a_offset;
791 if (molb->posres_xA.empty())
793 gmx_incons("Position restraint coordinates are missing");
795 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
796 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
797 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
798 if (!molb->posres_xB.empty())
800 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
801 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
802 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
806 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
807 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
808 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
810 /* Set the parameter index for idef->iparams_posre */
815 static void set_fbposres_params(t_idef *idef, const gmx_molblock_t *molb,
816 int i0, int a_offset)
822 il = &idef->il[F_FBPOSRES];
824 idef->iparams_fbposres_nalloc = i1;
825 srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
826 for (i = i0; i < i1; i++)
828 ip = &idef->iparams_fbposres[i];
829 /* Copy the force constants */
830 *ip = idef->iparams[il->iatoms[i*2]];
831 a_molb = il->iatoms[i*2+1] - a_offset;
832 if (molb->posres_xA.empty())
834 gmx_incons("Position restraint coordinates are missing");
836 /* Take flat-bottom posres reference from normal position restraints */
837 ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
838 ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
839 ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
840 /* Note: no B-type for flat-bottom posres */
842 /* Set the parameter index for idef->iparams_posre */
847 static void gen_local_top(const gmx_mtop_t *mtop,
848 bool freeEnergyInteractionsAtEnd,
852 int srcnr, destnr, ftype, natoms, nposre_old, nfbposre_old;
853 const gmx_ffparams_t *ffp;
856 gmx_mtop_atomloop_all_t aloop;
859 /* no copying of pointers possible now */
861 top->atomtypes.nr = mtop->atomtypes.nr;
862 if (mtop->atomtypes.atomnumber)
864 snew(top->atomtypes.atomnumber, top->atomtypes.nr);
865 std::copy(mtop->atomtypes.atomnumber, mtop->atomtypes.atomnumber + top->atomtypes.nr, top->atomtypes.atomnumber);
869 top->atomtypes.atomnumber = nullptr;
872 ffp = &mtop->ffparams;
875 idef->ntypes = ffp->ntypes;
876 idef->atnr = ffp->atnr;
877 /* we can no longer copy the pointers to the mtop members,
878 * because they will become invalid as soon as mtop gets free'd.
879 * We also need to make sure to only operate on valid data!
884 snew(idef->functype, ffp->ntypes);
885 std::copy(ffp->functype, ffp->functype + ffp->ntypes, idef->functype);
889 idef->functype = nullptr;
893 snew(idef->iparams, ffp->ntypes);
894 std::copy(ffp->iparams, ffp->iparams + ffp->ntypes, idef->iparams);
898 idef->iparams = nullptr;
900 idef->iparams_posres = nullptr;
901 idef->iparams_posres_nalloc = 0;
902 idef->iparams_fbposres = nullptr;
903 idef->iparams_fbposres_nalloc = 0;
904 idef->fudgeQQ = ffp->fudgeQQ;
905 idef->cmap_grid.ngrid = ffp->cmap_grid.ngrid;
906 idef->cmap_grid.grid_spacing = ffp->cmap_grid.grid_spacing;
907 if (ffp->cmap_grid.cmapdata)
909 snew(idef->cmap_grid.cmapdata, ffp->cmap_grid.ngrid);
910 std::copy(ffp->cmap_grid.cmapdata, ffp->cmap_grid.cmapdata + ffp->cmap_grid.ngrid, idef->cmap_grid.cmapdata);
914 idef->cmap_grid.cmapdata = nullptr;
916 idef->ilsort = ilsortUNKNOWN;
918 init_block(&top->cgs);
919 init_blocka(&top->excls);
920 for (ftype = 0; ftype < F_NRE; ftype++)
922 idef->il[ftype].nr = 0;
923 idef->il[ftype].nalloc = 0;
924 idef->il[ftype].iatoms = nullptr;
928 for (const gmx_molblock_t &molb : mtop->molblock)
930 const gmx_moltype_t &molt = mtop->moltype[molb.type];
932 srcnr = molt.atoms.nr;
935 blockcat(&top->cgs, &molt.cgs, molb.nmol);
937 blockacat(&top->excls, &molt.excls, molb.nmol, destnr, srcnr);
939 nposre_old = idef->il[F_POSRES].nr;
940 nfbposre_old = idef->il[F_FBPOSRES].nr;
941 for (ftype = 0; ftype < F_NRE; ftype++)
944 ftype == F_CONSTR && molt.ilist[F_CONSTRNC].nr > 0)
946 /* Merge all constrains into one ilist.
947 * This simplifies the constraint code.
949 for (int mol = 0; mol < molb.nmol; mol++)
951 ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTR],
952 1, destnr + mol*srcnr, srcnr);
953 ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTRNC],
954 1, destnr + mol*srcnr, srcnr);
957 else if (!(bMergeConstr && ftype == F_CONSTRNC))
959 ilistcat(ftype, &idef->il[ftype], &molt.ilist[ftype],
960 molb.nmol, destnr, srcnr);
963 if (idef->il[F_POSRES].nr > nposre_old)
965 /* Executing this line line stops gmxdump -sys working
966 * correctly. I'm not aware there's an elegant fix. */
967 set_posres_params(idef, &molb, nposre_old/2, natoms);
969 if (idef->il[F_FBPOSRES].nr > nfbposre_old)
971 set_fbposres_params(idef, &molb, nfbposre_old/2, natoms);
974 natoms += molb.nmol*srcnr;
977 if (mtop->bIntermolecularInteractions)
979 for (ftype = 0; ftype < F_NRE; ftype++)
981 ilistcat(ftype, &idef->il[ftype], &mtop->intermolecular_ilist[ftype],
986 if (freeEnergyInteractionsAtEnd && gmx_mtop_bondeds_free_energy(mtop))
988 snew(qA, mtop->natoms);
989 snew(qB, mtop->natoms);
990 aloop = gmx_mtop_atomloop_all_init(mtop);
992 while (gmx_mtop_atomloop_all_next(aloop, &ag, &atom))
997 gmx_sort_ilist_fe(&top->idef, qA, qB);
1003 top->idef.ilsort = ilsortNO_FE;
1008 gmx_mtop_generate_local_top(const gmx_mtop_t *mtop,
1009 bool freeEnergyInteractionsAtEnd)
1011 gmx_localtop_t *top;
1015 gen_local_top(mtop, freeEnergyInteractionsAtEnd, true, top);
1020 /*! \brief Fills an array with molecule begin/end atom indices
1022 * \param[in] mtop The global topology
1023 * \param[out] index Array of size nr. of molecules + 1 to be filled with molecule begin/end indices
1025 static void fillMoleculeIndices(const gmx_mtop_t &mtop,
1026 gmx::ArrayRef<int> index)
1028 int globalAtomIndex = 0;
1029 int globalMolIndex = 0;
1030 index[globalMolIndex] = globalAtomIndex;
1031 for (const gmx_molblock_t &molb : mtop.molblock)
1033 int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
1034 for (int mol = 0; mol < molb.nmol; mol++)
1036 globalAtomIndex += numAtomsPerMolecule;
1037 globalMolIndex += 1;
1038 index[globalMolIndex] = globalAtomIndex;
1043 gmx::RangePartitioning gmx_mtop_molecules(const gmx_mtop_t &mtop)
1045 gmx::RangePartitioning mols;
1047 for (const gmx_molblock_t &molb : mtop.molblock)
1049 int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
1050 for (int mol = 0; mol < molb.nmol; mol++)
1052 mols.appendBlock(numAtomsPerMolecule);
1059 /*! \brief Creates and returns a deprecated t_block struct with molecule indices
1061 * \param[in] mtop The global topology
1063 static t_block gmx_mtop_molecules_t_block(const gmx_mtop_t &mtop)
1067 mols.nr = gmx_mtop_num_molecules(mtop);
1068 mols.nalloc_index = mols.nr + 1;
1069 snew(mols.index, mols.nalloc_index);
1071 fillMoleculeIndices(mtop, gmx::arrayRefFromArray(mols.index, mols.nr + 1));
1076 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop, bool freeMTop)
1078 gmx_localtop_t ltop;
1081 gen_local_top(mtop, false, FALSE, <op);
1082 ltop.idef.ilsort = ilsortUNKNOWN;
1084 top.name = mtop->name;
1085 top.idef = ltop.idef;
1086 top.atomtypes = ltop.atomtypes;
1088 top.excls = ltop.excls;
1089 top.atoms = gmx_mtop_global_atoms(mtop);
1090 top.mols = gmx_mtop_molecules_t_block(*mtop);
1091 top.bIntermolecularInteractions = mtop->bIntermolecularInteractions;
1092 top.symtab = mtop->symtab;
1096 // Clear pointers and counts, such that the pointers copied to top
1097 // keep pointing to valid data after destroying mtop.
1098 mtop->symtab.symbuf = nullptr;
1099 mtop->symtab.nr = 0;
1105 std::vector<size_t> get_atom_index(const gmx_mtop_t *mtop)
1108 std::vector<size_t> atom_index;
1109 gmx_mtop_atomloop_block_t aloopb = gmx_mtop_atomloop_block_init(mtop);
1112 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
1114 if (atom->ptype == eptAtom)
1116 atom_index.push_back(j);
1123 void convertAtomsToMtop(t_symtab *symtab,
1128 mtop->symtab = *symtab;
1132 mtop->moltype.clear();
1133 mtop->moltype.resize(1);
1134 mtop->moltype.back().atoms = *atoms;
1136 mtop->molblock.resize(1);
1137 mtop->molblock[0].type = 0;
1138 mtop->molblock[0].nmol = 1;
1140 mtop->bIntermolecularInteractions = FALSE;
1142 mtop->natoms = atoms->nr;
1144 mtop->haveMoleculeIndices = false;
1146 gmx_mtop_finalize(mtop);