Fix legacy topology memory error
[gromacs.git] / src / gromacs / topology / mtop_util.cpp
blobcc439f21ab83ee5ba17dc632ca4ee6e9d08508df
1 /*
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.
35 #include "gmxpre.h"
37 #include "mtop_util.h"
39 #include <limits.h>
40 #include <stddef.h>
41 #include <stdio.h>
42 #include <stdlib.h>
43 #include <string.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)
59 int maxresnr = 0;
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;
76 return maxresnr;
79 static void buildMolblockIndices(gmx_mtop_t *mtop)
81 mtop->moleculeBlockIndices.resize(mtop->molblock.size());
83 int atomIndex = 0;
84 int residueIndex = 0;
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)
111 char *env;
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;
121 else
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");
130 if (env != nullptr)
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)
149 typecount[i] = 0;
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)
156 int tpi;
157 if (state == 0)
159 tpi = atoms.atom[i].type;
161 else
163 tpi = atoms.atom[i].typeB;
165 typecount[tpi] += molb.nmol;
170 int ncg_mtop(const gmx_mtop_t *mtop)
172 int ncg = 0;
173 for (const gmx_molblock_t &molb : mtop->molblock)
175 ncg += molb.nmol*mtop->moltype[molb.type].cgs.nr;
178 return ncg;
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;
188 return numMolecules;
191 int gmx_mtop_nres(const gmx_mtop_t *mtop)
193 int nres = 0;
194 for (const gmx_molblock_t &molb : mtop->molblock)
196 nres += molb.nmol*mtop->moltype[molb.type].atoms.nres;
198 return 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++)
212 cgs.index[i] = i;
218 typedef struct gmx_mtop_atomloop_all
220 const gmx_mtop_t *mtop;
221 size_t mblock;
222 const t_atoms *atoms;
223 int mol;
224 int maxresnr;
225 int at_local;
226 int at_global;
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;
234 snew(aloop, 1);
236 aloop->mtop = mtop;
237 aloop->mblock = 0;
238 aloop->atoms =
239 &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
240 aloop->mol = 0;
241 aloop->maxresnr = mtop->maxresnr;
242 aloop->at_local = -1;
243 aloop->at_global = -1;
245 return aloop;
248 static void gmx_mtop_atomloop_all_destroy(gmx_mtop_atomloop_all_t aloop)
250 sfree(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");
261 aloop->at_local++;
262 aloop->at_global++;
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;
271 aloop->mol++;
272 aloop->at_local = 0;
273 if (aloop->mol >= aloop->mtop->molblock[aloop->mblock].nmol)
275 aloop->mblock++;
276 if (aloop->mblock >= aloop->mtop->molblock.size())
278 gmx_mtop_atomloop_all_destroy(aloop);
279 return FALSE;
281 aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
282 aloop->mol = 0;
286 *at_global = aloop->at_global;
287 *atom = &aloop->atoms->atom[aloop->at_local];
289 return TRUE;
292 void gmx_mtop_atomloop_all_names(gmx_mtop_atomloop_all_t aloop,
293 char **atomname, int *resnr, char **resname)
295 int resind_mol;
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,
309 int *at_mol)
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;
318 size_t mblock;
319 const t_atoms *atoms;
320 int at_local;
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;
328 snew(aloop, 1);
330 aloop->mtop = mtop;
331 aloop->mblock = 0;
332 aloop->atoms = &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
333 aloop->at_local = -1;
335 return aloop;
338 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop)
340 sfree(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");
351 aloop->at_local++;
353 if (aloop->at_local >= aloop->atoms->nr)
355 aloop->mblock++;
356 if (aloop->mblock >= aloop->mtop->molblock.size())
358 gmx_mtop_atomloop_block_destroy(aloop);
359 return FALSE;
361 aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
362 aloop->at_local = 0;
365 *atom = &aloop->atoms->atom[aloop->at_local];
366 *nmol = aloop->mtop->molblock[aloop->mblock].nmol;
368 return TRUE;
371 typedef struct gmx_mtop_ilistloop
373 const gmx_mtop_t *mtop;
374 int mblock;
375 } t_gmx_mtop_ilist;
377 gmx_mtop_ilistloop_t
378 gmx_mtop_ilistloop_init(const gmx_mtop_t *mtop)
380 struct gmx_mtop_ilistloop *iloop;
382 snew(iloop, 1);
384 iloop->mtop = mtop;
385 iloop->mblock = -1;
387 return iloop;
390 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
392 sfree(iloop);
395 gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
396 const t_ilist **ilist_mol, int *nmol)
398 if (iloop == nullptr)
400 gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
403 iloop->mblock++;
404 if (iloop->mblock >= static_cast<int>(iloop->mtop->molblock.size()))
406 if (iloop->mblock == static_cast<int>(iloop->mtop->molblock.size()) &&
407 iloop->mtop->bIntermolecularInteractions)
409 *ilist_mol = iloop->mtop->intermolecular_ilist;
410 *nmol = 1;
411 return TRUE;
414 gmx_mtop_ilistloop_destroy(iloop);
415 return FALSE;
418 *ilist_mol =
419 iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
421 *nmol = iloop->mtop->molblock[iloop->mblock].nmol;
423 return TRUE;
425 typedef struct gmx_mtop_ilistloop_all
427 const gmx_mtop_t *mtop;
428 size_t mblock;
429 int mol;
430 int a_offset;
431 } t_gmx_mtop_ilist_all;
433 gmx_mtop_ilistloop_all_t
434 gmx_mtop_ilistloop_all_init(const gmx_mtop_t *mtop)
436 struct gmx_mtop_ilistloop_all *iloop;
438 snew(iloop, 1);
440 iloop->mtop = mtop;
441 iloop->mblock = 0;
442 iloop->mol = -1;
443 iloop->a_offset = 0;
445 return iloop;
448 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop)
450 sfree(iloop);
453 gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
454 const t_ilist **ilist_mol, int *atnr_offset)
457 if (iloop == nullptr)
459 gmx_incons("gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
462 if (iloop->mol >= 0)
464 iloop->a_offset += iloop->mtop->moleculeBlockIndices[iloop->mblock].numAtomsPerMolecule;
467 iloop->mol++;
469 /* Inter-molecular interactions, if present, are indexed with
470 * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
471 * check for this value in this conditional.
473 if (iloop->mblock == iloop->mtop->molblock.size() ||
474 iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
476 iloop->mblock++;
477 iloop->mol = 0;
478 if (iloop->mblock >= iloop->mtop->molblock.size())
480 if (iloop->mblock == iloop->mtop->molblock.size() &&
481 iloop->mtop->bIntermolecularInteractions)
483 *ilist_mol = iloop->mtop->intermolecular_ilist;
484 *atnr_offset = 0;
485 return TRUE;
488 gmx_mtop_ilistloop_all_destroy(iloop);
489 return FALSE;
493 *ilist_mol =
494 iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
496 *atnr_offset = iloop->a_offset;
498 return TRUE;
501 int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype)
503 gmx_mtop_ilistloop_t iloop;
504 const t_ilist *il;
505 int n, nmol;
507 n = 0;
509 iloop = gmx_mtop_ilistloop_init(mtop);
510 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
512 n += nmol*il[ftype].nr/(1+NRAL(ftype));
515 if (mtop->bIntermolecularInteractions)
517 n += mtop->intermolecular_ilist[ftype].nr/(1+NRAL(ftype));
520 return n;
523 t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
525 t_block cgs_gl;
526 /* In most cases this is too much, but we realloc at the end */
527 snew(cgs_gl.index, mtop->natoms+1);
529 cgs_gl.nr = 0;
530 cgs_gl.index[0] = 0;
531 for (const gmx_molblock_t &molb : mtop->molblock)
533 const t_block &cgs_mol = mtop->moltype[molb.type].cgs;
534 for (int mol = 0; mol < molb.nmol; mol++)
536 for (int cg = 0; cg < cgs_mol.nr; cg++)
538 cgs_gl.index[cgs_gl.nr + 1] =
539 cgs_gl.index[cgs_gl.nr] +
540 cgs_mol.index[cg + 1] - cgs_mol.index[cg];
541 cgs_gl.nr++;
545 cgs_gl.nalloc_index = cgs_gl.nr + 1;
546 srenew(cgs_gl.index, cgs_gl.nalloc_index);
548 return cgs_gl;
551 static void atomcat(t_atoms *dest, const t_atoms *src, int copies,
552 int maxres_renum, int *maxresnr)
554 int i, j, l, size;
555 int srcnr = src->nr;
556 int destnr = dest->nr;
558 if (dest->nr == 0)
560 dest->haveMass = src->haveMass;
561 dest->haveType = src->haveType;
562 dest->haveCharge = src->haveCharge;
563 dest->haveBState = src->haveBState;
564 dest->havePdbInfo = src->havePdbInfo;
566 else
568 dest->haveMass = dest->haveMass && src->haveMass;
569 dest->haveType = dest->haveType && src->haveType;
570 dest->haveCharge = dest->haveCharge && src->haveCharge;
571 dest->haveBState = dest->haveBState && src->haveBState;
572 dest->havePdbInfo = dest->havePdbInfo && src->havePdbInfo;
575 if (srcnr)
577 size = destnr+copies*srcnr;
578 srenew(dest->atom, size);
579 srenew(dest->atomname, size);
580 if (dest->haveType)
582 srenew(dest->atomtype, size);
583 if (dest->haveBState)
585 srenew(dest->atomtypeB, size);
588 if (dest->havePdbInfo)
590 srenew(dest->pdbinfo, size);
593 if (src->nres)
595 size = dest->nres+copies*src->nres;
596 srenew(dest->resinfo, size);
599 /* residue information */
600 for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
602 memcpy((char *) &(dest->resinfo[l]), (char *) &(src->resinfo[0]),
603 (size_t)(src->nres*sizeof(src->resinfo[0])));
606 for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
608 memcpy((char *) &(dest->atom[l]), (char *) &(src->atom[0]),
609 (size_t)(srcnr*sizeof(src->atom[0])));
610 memcpy((char *) &(dest->atomname[l]), (char *) &(src->atomname[0]),
611 (size_t)(srcnr*sizeof(src->atomname[0])));
612 if (dest->haveType)
614 memcpy((char *) &(dest->atomtype[l]), (char *) &(src->atomtype[0]),
615 (size_t)(srcnr*sizeof(src->atomtype[0])));
616 if (dest->haveBState)
618 memcpy((char *) &(dest->atomtypeB[l]), (char *) &(src->atomtypeB[0]),
619 (size_t)(srcnr*sizeof(src->atomtypeB[0])));
622 if (dest->havePdbInfo)
624 memcpy((char *) &(dest->pdbinfo[l]), (char *) &(src->pdbinfo[0]),
625 (size_t)(srcnr*sizeof(src->pdbinfo[0])));
629 /* Increment residue indices */
630 for (l = destnr, j = 0; (j < copies); j++)
632 for (i = 0; (i < srcnr); i++, l++)
634 dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;
638 if (src->nres <= maxres_renum)
640 /* Single residue molecule, continue counting residues */
641 for (j = 0; (j < copies); j++)
643 for (l = 0; l < src->nres; l++)
645 (*maxresnr)++;
646 dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;
651 dest->nres += copies*src->nres;
652 dest->nr += copies*src->nr;
655 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
657 t_atoms atoms;
659 init_t_atoms(&atoms, 0, FALSE);
661 int maxresnr = mtop->maxresnr;
662 for (const gmx_molblock_t &molb : mtop->molblock)
664 atomcat(&atoms, &mtop->moltype[molb.type].atoms, molb.nmol,
665 mtop->maxres_renum, &maxresnr);
668 return atoms;
672 * The cat routines below are old code from src/kernel/topcat.c
675 static void blockcat(t_block *dest, const t_block *src, int copies)
677 int i, j, l, nra, size;
679 if (src->nr)
681 size = (dest->nr+copies*src->nr+1);
682 srenew(dest->index, size);
685 nra = dest->index[dest->nr];
686 for (l = dest->nr, j = 0; (j < copies); j++)
688 for (i = 0; (i < src->nr); i++)
690 dest->index[l++] = nra + src->index[i];
692 if (src->nr > 0)
694 nra += src->index[src->nr];
697 dest->nr += copies*src->nr;
698 dest->index[dest->nr] = nra;
701 static void blockacat(t_blocka *dest, const t_blocka *src, int copies,
702 int dnum, int snum)
704 int i, j, l, size;
705 int destnr = dest->nr;
706 int destnra = dest->nra;
708 if (src->nr)
710 size = (dest->nr+copies*src->nr+1);
711 srenew(dest->index, size);
713 if (src->nra)
715 size = (dest->nra+copies*src->nra);
716 srenew(dest->a, size);
719 for (l = destnr, j = 0; (j < copies); j++)
721 for (i = 0; (i < src->nr); i++)
723 dest->index[l++] = dest->nra+src->index[i];
725 dest->nra += src->nra;
727 for (l = destnra, j = 0; (j < copies); j++)
729 for (i = 0; (i < src->nra); i++)
731 dest->a[l++] = dnum+src->a[i];
733 dnum += snum;
734 dest->nr += src->nr;
736 dest->index[dest->nr] = dest->nra;
739 static void ilistcat(int ftype, t_ilist *dest, const t_ilist *src, int copies,
740 int dnum, int snum)
742 int nral, c, i, a;
744 nral = NRAL(ftype);
746 dest->nalloc = dest->nr + copies*src->nr;
747 srenew(dest->iatoms, dest->nalloc);
749 for (c = 0; c < copies; c++)
751 for (i = 0; i < src->nr; )
753 dest->iatoms[dest->nr++] = src->iatoms[i++];
754 for (a = 0; a < nral; a++)
756 dest->iatoms[dest->nr++] = dnum + src->iatoms[i++];
759 dnum += snum;
763 static void set_posres_params(t_idef *idef, const gmx_molblock_t *molb,
764 int i0, int a_offset)
766 t_ilist *il;
767 int i1, i, a_molb;
768 t_iparams *ip;
770 il = &idef->il[F_POSRES];
771 i1 = il->nr/2;
772 idef->iparams_posres_nalloc = i1;
773 srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
774 for (i = i0; i < i1; i++)
776 ip = &idef->iparams_posres[i];
777 /* Copy the force constants */
778 *ip = idef->iparams[il->iatoms[i*2]];
779 a_molb = il->iatoms[i*2+1] - a_offset;
780 if (molb->posres_xA.empty())
782 gmx_incons("Position restraint coordinates are missing");
784 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
785 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
786 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
787 if (!molb->posres_xB.empty())
789 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
790 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
791 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
793 else
795 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
796 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
797 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
799 /* Set the parameter index for idef->iparams_posre */
800 il->iatoms[i*2] = i;
804 static void set_fbposres_params(t_idef *idef, const gmx_molblock_t *molb,
805 int i0, int a_offset)
807 t_ilist *il;
808 int i1, i, a_molb;
809 t_iparams *ip;
811 il = &idef->il[F_FBPOSRES];
812 i1 = il->nr/2;
813 idef->iparams_fbposres_nalloc = i1;
814 srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
815 for (i = i0; i < i1; i++)
817 ip = &idef->iparams_fbposres[i];
818 /* Copy the force constants */
819 *ip = idef->iparams[il->iatoms[i*2]];
820 a_molb = il->iatoms[i*2+1] - a_offset;
821 if (molb->posres_xA.empty())
823 gmx_incons("Position restraint coordinates are missing");
825 /* Take flat-bottom posres reference from normal position restraints */
826 ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
827 ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
828 ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
829 /* Note: no B-type for flat-bottom posres */
831 /* Set the parameter index for idef->iparams_posre */
832 il->iatoms[i*2] = i;
836 static void gen_local_top(const gmx_mtop_t *mtop,
837 bool freeEnergyInteractionsAtEnd,
838 bool bMergeConstr,
839 gmx_localtop_t *top)
841 int srcnr, destnr, ftype, natoms, nposre_old, nfbposre_old;
842 const gmx_ffparams_t *ffp;
843 t_idef *idef;
844 real *qA, *qB;
845 gmx_mtop_atomloop_all_t aloop;
846 int ag;
848 /* no copying of pointers possible now */
850 top->atomtypes.nr = mtop->atomtypes.nr;
851 if (mtop->atomtypes.atomnumber)
853 snew(top->atomtypes.atomnumber, top->atomtypes.nr);
854 std::copy(mtop->atomtypes.atomnumber, mtop->atomtypes.atomnumber + top->atomtypes.nr, top->atomtypes.atomnumber);
856 else
858 top->atomtypes.atomnumber = nullptr;
861 ffp = &mtop->ffparams;
863 idef = &top->idef;
864 idef->ntypes = ffp->ntypes;
865 idef->atnr = ffp->atnr;
866 /* we can no longer copy the pointers to the mtop members,
867 * because they will become invalid as soon as mtop gets free'd.
868 * We also need to make sure to only operate on valid data!
871 if (ffp->functype)
873 snew(idef->functype, ffp->ntypes);
874 std::copy(ffp->functype, ffp->functype + ffp->ntypes, idef->functype);
876 else
878 idef->functype = nullptr;
880 if (ffp->iparams)
882 snew(idef->iparams, ffp->ntypes);
883 std::copy(ffp->iparams, ffp->iparams + ffp->ntypes, idef->iparams);
885 else
887 idef->iparams = nullptr;
889 idef->iparams_posres = nullptr;
890 idef->iparams_posres_nalloc = 0;
891 idef->iparams_fbposres = nullptr;
892 idef->iparams_fbposres_nalloc = 0;
893 idef->fudgeQQ = ffp->fudgeQQ;
894 idef->cmap_grid.ngrid = ffp->cmap_grid.ngrid;
895 idef->cmap_grid.grid_spacing = ffp->cmap_grid.grid_spacing;
896 if (ffp->cmap_grid.cmapdata)
898 snew(idef->cmap_grid.cmapdata, ffp->cmap_grid.ngrid);
899 std::copy(ffp->cmap_grid.cmapdata, ffp->cmap_grid.cmapdata + ffp->cmap_grid.ngrid, idef->cmap_grid.cmapdata);
901 else
903 idef->cmap_grid.cmapdata = nullptr;
905 idef->ilsort = ilsortUNKNOWN;
907 init_block(&top->cgs);
908 init_blocka(&top->excls);
909 for (ftype = 0; ftype < F_NRE; ftype++)
911 idef->il[ftype].nr = 0;
912 idef->il[ftype].nalloc = 0;
913 idef->il[ftype].iatoms = nullptr;
916 natoms = 0;
917 for (const gmx_molblock_t &molb : mtop->molblock)
919 const gmx_moltype_t &molt = mtop->moltype[molb.type];
921 srcnr = molt.atoms.nr;
922 destnr = natoms;
924 blockcat(&top->cgs, &molt.cgs, molb.nmol);
926 blockacat(&top->excls, &molt.excls, molb.nmol, destnr, srcnr);
928 nposre_old = idef->il[F_POSRES].nr;
929 nfbposre_old = idef->il[F_FBPOSRES].nr;
930 for (ftype = 0; ftype < F_NRE; ftype++)
932 if (bMergeConstr &&
933 ftype == F_CONSTR && molt.ilist[F_CONSTRNC].nr > 0)
935 /* Merge all constrains into one ilist.
936 * This simplifies the constraint code.
938 for (int mol = 0; mol < molb.nmol; mol++)
940 ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTR],
941 1, destnr + mol*srcnr, srcnr);
942 ilistcat(ftype, &idef->il[F_CONSTR], &molt.ilist[F_CONSTRNC],
943 1, destnr + mol*srcnr, srcnr);
946 else if (!(bMergeConstr && ftype == F_CONSTRNC))
948 ilistcat(ftype, &idef->il[ftype], &molt.ilist[ftype],
949 molb.nmol, destnr, srcnr);
952 if (idef->il[F_POSRES].nr > nposre_old)
954 /* Executing this line line stops gmxdump -sys working
955 * correctly. I'm not aware there's an elegant fix. */
956 set_posres_params(idef, &molb, nposre_old/2, natoms);
958 if (idef->il[F_FBPOSRES].nr > nfbposre_old)
960 set_fbposres_params(idef, &molb, nfbposre_old/2, natoms);
963 natoms += molb.nmol*srcnr;
966 if (mtop->bIntermolecularInteractions)
968 for (ftype = 0; ftype < F_NRE; ftype++)
970 ilistcat(ftype, &idef->il[ftype], &mtop->intermolecular_ilist[ftype],
971 1, 0, mtop->natoms);
975 if (freeEnergyInteractionsAtEnd && gmx_mtop_bondeds_free_energy(mtop))
977 snew(qA, mtop->natoms);
978 snew(qB, mtop->natoms);
979 aloop = gmx_mtop_atomloop_all_init(mtop);
980 const t_atom *atom;
981 while (gmx_mtop_atomloop_all_next(aloop, &ag, &atom))
983 qA[ag] = atom->q;
984 qB[ag] = atom->qB;
986 gmx_sort_ilist_fe(&top->idef, qA, qB);
987 sfree(qA);
988 sfree(qB);
990 else
992 top->idef.ilsort = ilsortNO_FE;
996 gmx_localtop_t *
997 gmx_mtop_generate_local_top(const gmx_mtop_t *mtop,
998 bool freeEnergyInteractionsAtEnd)
1000 gmx_localtop_t *top;
1002 snew(top, 1);
1004 gen_local_top(mtop, freeEnergyInteractionsAtEnd, true, top);
1006 return top;
1009 /*! \brief Fills an array with molecule begin/end atom indices
1011 * \param[in] mtop The global topology
1012 * \param[out] index Array of size nr. of molecules + 1 to be filled with molecule begin/end indices
1014 static void fillMoleculeIndices(const gmx_mtop_t &mtop,
1015 gmx::ArrayRef<int> index)
1017 int globalAtomIndex = 0;
1018 int globalMolIndex = 0;
1019 index[globalMolIndex] = globalAtomIndex;
1020 for (const gmx_molblock_t &molb : mtop.molblock)
1022 int numAtomsPerMolecule = mtop.moltype[molb.type].atoms.nr;
1023 for (int mol = 0; mol < molb.nmol; mol++)
1025 globalAtomIndex += numAtomsPerMolecule;
1026 globalMolIndex += 1;
1027 index[globalMolIndex] = globalAtomIndex;
1032 gmx::BlockRanges gmx_mtop_molecules(const gmx_mtop_t &mtop)
1034 gmx::BlockRanges mols;
1036 mols.index.resize(gmx_mtop_num_molecules(mtop) + 1);
1038 fillMoleculeIndices(mtop, mols.index);
1040 return mols;
1043 /*! \brief Creates and returns a deprecated t_block struct with molecule indices
1045 * \param[in] mtop The global topology
1047 static t_block gmx_mtop_molecules_t_block(const gmx_mtop_t &mtop)
1049 t_block mols;
1051 mols.nr = gmx_mtop_num_molecules(mtop);
1052 mols.nalloc_index = mols.nr + 1;
1053 snew(mols.index, mols.nalloc_index);
1055 fillMoleculeIndices(mtop, gmx::arrayRefFromArray(mols.index, mols.nr + 1));
1057 return mols;
1060 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop, bool freeMTop)
1062 gmx_localtop_t ltop;
1063 t_topology top;
1065 gen_local_top(mtop, false, FALSE, &ltop);
1066 ltop.idef.ilsort = ilsortUNKNOWN;
1068 top.name = mtop->name;
1069 top.idef = ltop.idef;
1070 top.atomtypes = ltop.atomtypes;
1071 top.cgs = ltop.cgs;
1072 top.excls = ltop.excls;
1073 top.atoms = gmx_mtop_global_atoms(mtop);
1074 top.mols = gmx_mtop_molecules_t_block(*mtop);
1075 top.bIntermolecularInteractions = mtop->bIntermolecularInteractions;
1076 top.symtab = mtop->symtab;
1078 if (freeMTop)
1080 // Clear pointers and counts, such that the pointers copied to top
1081 // keep pointing to valid data after destroying mtop.
1082 mtop->symtab.symbuf = nullptr;
1083 mtop->symtab.nr = 0;
1086 return top;
1089 std::vector<size_t> get_atom_index(const gmx_mtop_t *mtop)
1092 std::vector<size_t> atom_index;
1093 gmx_mtop_atomloop_block_t aloopb = gmx_mtop_atomloop_block_init(mtop);
1094 const t_atom *atom;
1095 int nmol, j = 0;
1096 while (gmx_mtop_atomloop_block_next(aloopb, &atom, &nmol))
1098 if (atom->ptype == eptAtom)
1100 atom_index.push_back(j);
1102 j++;
1104 return atom_index;
1107 void convertAtomsToMtop(t_symtab *symtab,
1108 char **name,
1109 t_atoms *atoms,
1110 gmx_mtop_t *mtop)
1112 mtop->symtab = *symtab;
1114 mtop->name = name;
1116 mtop->moltype.clear();
1117 mtop->moltype.resize(1);
1118 mtop->moltype.back().atoms = *atoms;
1120 mtop->molblock.resize(1);
1121 mtop->molblock[0].type = 0;
1122 mtop->molblock[0].nmol = 1;
1124 mtop->bIntermolecularInteractions = FALSE;
1126 mtop->natoms = atoms->nr;
1128 mtop->haveMoleculeIndices = false;
1130 gmx_mtop_finalize(mtop);