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
)
64 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
66 atoms
= &mtop
->moltype
[mt
].atoms
;
67 if (atoms
->nres
> maxres_renum
)
69 for (r
= 0; r
< atoms
->nres
; r
++)
71 if (atoms
->resinfo
[r
].nr
> maxresnr
)
73 maxresnr
= atoms
->resinfo
[r
].nr
;
82 static void finalizeMolblocks(gmx_mtop_t
*mtop
)
86 int residueNumberStart
= mtop
->maxresnr
+ 1;
87 int moleculeIndexStart
= 0;
88 for (int mb
= 0; mb
< mtop
->nmolblock
; mb
++)
90 gmx_molblock_t
*molb
= &mtop
->molblock
[mb
];
91 int numResPerMol
= mtop
->moltype
[molb
->type
].atoms
.nres
;
92 molb
->globalAtomStart
= atomIndex
;
93 molb
->globalResidueStart
= residueIndex
;
94 atomIndex
+= molb
->nmol
*molb
->natoms_mol
;
95 residueIndex
+= molb
->nmol
*numResPerMol
;
96 molb
->globalAtomEnd
= atomIndex
;
97 molb
->residueNumberStart
= residueNumberStart
;
98 if (numResPerMol
<= mtop
->maxres_renum
)
100 residueNumberStart
+= molb
->nmol
*numResPerMol
;
102 molb
->moleculeIndexStart
= moleculeIndexStart
;
103 moleculeIndexStart
+= molb
->nmol
;
107 void gmx_mtop_finalize(gmx_mtop_t
*mtop
)
111 if (mtop
->nmolblock
== 1 && mtop
->molblock
[0].nmol
== 1)
113 /* We have a single molecule only, no renumbering needed.
114 * This case also covers an mtop converted from pdb/gro/... input,
115 * so we retain the original residue numbering.
117 mtop
->maxres_renum
= 0;
121 /* We only renumber single residue molecules. Their intra-molecular
122 * residue numbering is anyhow irrelevant.
124 mtop
->maxres_renum
= 1;
127 env
= getenv("GMX_MAXRESRENUM");
130 sscanf(env
, "%d", &mtop
->maxres_renum
);
132 if (mtop
->maxres_renum
== -1)
134 /* -1 signals renumber residues in all molecules */
135 mtop
->maxres_renum
= INT_MAX
;
138 mtop
->maxresnr
= gmx_mtop_maxresnr(mtop
, mtop
->maxres_renum
);
140 finalizeMolblocks(mtop
);
143 void gmx_mtop_count_atomtypes(const gmx_mtop_t
*mtop
, int state
, int typecount
[])
145 int i
, mb
, nmol
, tpi
;
148 for (i
= 0; i
< mtop
->ffparams
.atnr
; ++i
)
152 for (mb
= 0; mb
< mtop
->nmolblock
; ++mb
)
154 nmol
= mtop
->molblock
[mb
].nmol
;
155 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
156 for (i
= 0; i
< atoms
->nr
; ++i
)
160 tpi
= atoms
->atom
[i
].type
;
164 tpi
= atoms
->atom
[i
].typeB
;
166 typecount
[tpi
] += nmol
;
171 int ncg_mtop(const gmx_mtop_t
*mtop
)
177 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
180 mtop
->molblock
[mb
].nmol
*
181 mtop
->moltype
[mtop
->molblock
[mb
].type
].cgs
.nr
;
187 int gmx_mtop_num_molecules(const gmx_mtop_t
&mtop
)
189 int numMolecules
= 0;
190 for (int mb
= 0; mb
< mtop
.nmolblock
; mb
++)
192 numMolecules
+= mtop
.molblock
[mb
].nmol
;
197 int gmx_mtop_nres(const gmx_mtop_t
*mtop
)
200 for (int mb
= 0; mb
< mtop
->nmolblock
; ++mb
)
203 mtop
->molblock
[mb
].nmol
*
204 mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
.nres
;
209 void gmx_mtop_remove_chargegroups(gmx_mtop_t
*mtop
)
215 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
217 cgs
= &mtop
->moltype
[mt
].cgs
;
218 if (cgs
->nr
< mtop
->moltype
[mt
].atoms
.nr
)
220 cgs
->nr
= mtop
->moltype
[mt
].atoms
.nr
;
221 srenew(cgs
->index
, cgs
->nr
+1);
222 for (i
= 0; i
< cgs
->nr
+1; i
++)
230 typedef struct gmx_mtop_atomloop_all
232 const gmx_mtop_t
*mtop
;
239 } t_gmx_mtop_atomloop_all
;
241 gmx_mtop_atomloop_all_t
242 gmx_mtop_atomloop_all_init(const gmx_mtop_t
*mtop
)
244 struct gmx_mtop_atomloop_all
*aloop
;
251 &mtop
->moltype
[mtop
->molblock
[aloop
->mblock
].type
].atoms
;
253 aloop
->maxresnr
= mtop
->maxresnr
;
254 aloop
->at_local
= -1;
255 aloop
->at_global
= -1;
260 static void gmx_mtop_atomloop_all_destroy(gmx_mtop_atomloop_all_t aloop
)
265 gmx_bool
gmx_mtop_atomloop_all_next(gmx_mtop_atomloop_all_t aloop
,
266 int *at_global
, const t_atom
**atom
)
268 if (aloop
== nullptr)
270 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
276 if (aloop
->at_local
>= aloop
->atoms
->nr
)
278 if (aloop
->atoms
->nres
<= aloop
->mtop
->maxres_renum
)
280 /* Single residue molecule, increase the count with one */
281 aloop
->maxresnr
+= aloop
->atoms
->nres
;
285 if (aloop
->mol
>= aloop
->mtop
->molblock
[aloop
->mblock
].nmol
)
288 if (aloop
->mblock
>= aloop
->mtop
->nmolblock
)
290 gmx_mtop_atomloop_all_destroy(aloop
);
293 aloop
->atoms
= &aloop
->mtop
->moltype
[aloop
->mtop
->molblock
[aloop
->mblock
].type
].atoms
;
298 *at_global
= aloop
->at_global
;
299 *atom
= &aloop
->atoms
->atom
[aloop
->at_local
];
304 void gmx_mtop_atomloop_all_names(gmx_mtop_atomloop_all_t aloop
,
305 char **atomname
, int *resnr
, char **resname
)
309 *atomname
= *(aloop
->atoms
->atomname
[aloop
->at_local
]);
310 resind_mol
= aloop
->atoms
->atom
[aloop
->at_local
].resind
;
311 *resnr
= aloop
->atoms
->resinfo
[resind_mol
].nr
;
312 if (aloop
->atoms
->nres
<= aloop
->mtop
->maxres_renum
)
314 *resnr
= aloop
->maxresnr
+ 1 + resind_mol
;
316 *resname
= *(aloop
->atoms
->resinfo
[resind_mol
].name
);
319 void gmx_mtop_atomloop_all_moltype(gmx_mtop_atomloop_all_t aloop
,
320 gmx_moltype_t
**moltype
, int *at_mol
)
322 *moltype
= &aloop
->mtop
->moltype
[aloop
->mtop
->molblock
[aloop
->mblock
].type
];
323 *at_mol
= aloop
->at_local
;
326 typedef struct gmx_mtop_atomloop_block
328 const gmx_mtop_t
*mtop
;
332 } t_gmx_mtop_atomloop_block
;
334 gmx_mtop_atomloop_block_t
335 gmx_mtop_atomloop_block_init(const gmx_mtop_t
*mtop
)
337 struct gmx_mtop_atomloop_block
*aloop
;
343 aloop
->atoms
= &mtop
->moltype
[mtop
->molblock
[aloop
->mblock
].type
].atoms
;
344 aloop
->at_local
= -1;
349 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop
)
354 gmx_bool
gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop
,
355 const t_atom
**atom
, int *nmol
)
357 if (aloop
== nullptr)
359 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
364 if (aloop
->at_local
>= aloop
->atoms
->nr
)
367 if (aloop
->mblock
>= aloop
->mtop
->nmolblock
)
369 gmx_mtop_atomloop_block_destroy(aloop
);
372 aloop
->atoms
= &aloop
->mtop
->moltype
[aloop
->mtop
->molblock
[aloop
->mblock
].type
].atoms
;
376 *atom
= &aloop
->atoms
->atom
[aloop
->at_local
];
377 *nmol
= aloop
->mtop
->molblock
[aloop
->mblock
].nmol
;
382 typedef struct gmx_mtop_ilistloop
384 const gmx_mtop_t
*mtop
;
389 gmx_mtop_ilistloop_init(const gmx_mtop_t
*mtop
)
391 struct gmx_mtop_ilistloop
*iloop
;
401 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop
)
406 gmx_bool
gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop
,
407 t_ilist
**ilist_mol
, int *nmol
)
409 if (iloop
== nullptr)
411 gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
415 if (iloop
->mblock
>= iloop
->mtop
->nmolblock
)
417 if (iloop
->mblock
== iloop
->mtop
->nmolblock
&&
418 iloop
->mtop
->bIntermolecularInteractions
)
420 *ilist_mol
= iloop
->mtop
->intermolecular_ilist
;
425 gmx_mtop_ilistloop_destroy(iloop
);
430 iloop
->mtop
->moltype
[iloop
->mtop
->molblock
[iloop
->mblock
].type
].ilist
;
432 *nmol
= iloop
->mtop
->molblock
[iloop
->mblock
].nmol
;
436 typedef struct gmx_mtop_ilistloop_all
438 const gmx_mtop_t
*mtop
;
442 } t_gmx_mtop_ilist_all
;
444 gmx_mtop_ilistloop_all_t
445 gmx_mtop_ilistloop_all_init(const gmx_mtop_t
*mtop
)
447 struct gmx_mtop_ilistloop_all
*iloop
;
459 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop
)
464 gmx_bool
gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop
,
465 t_ilist
**ilist_mol
, int *atnr_offset
)
468 if (iloop
== nullptr)
470 gmx_incons("gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
475 iloop
->a_offset
+= iloop
->mtop
->molblock
[iloop
->mblock
].natoms_mol
;
480 /* Inter-molecular interactions, if present, are indexed with
481 * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
482 * check for this value in this conditional.
484 if (iloop
->mblock
== iloop
->mtop
->nmolblock
||
485 iloop
->mol
>= iloop
->mtop
->molblock
[iloop
->mblock
].nmol
)
489 if (iloop
->mblock
>= iloop
->mtop
->nmolblock
)
491 if (iloop
->mblock
== iloop
->mtop
->nmolblock
&&
492 iloop
->mtop
->bIntermolecularInteractions
)
494 *ilist_mol
= iloop
->mtop
->intermolecular_ilist
;
499 gmx_mtop_ilistloop_all_destroy(iloop
);
505 iloop
->mtop
->moltype
[iloop
->mtop
->molblock
[iloop
->mblock
].type
].ilist
;
507 *atnr_offset
= iloop
->a_offset
;
512 int gmx_mtop_ftype_count(const gmx_mtop_t
*mtop
, int ftype
)
514 gmx_mtop_ilistloop_t iloop
;
520 iloop
= gmx_mtop_ilistloop_init(mtop
);
521 while (gmx_mtop_ilistloop_next(iloop
, &il
, &nmol
))
523 n
+= nmol
*il
[ftype
].nr
/(1+NRAL(ftype
));
526 if (mtop
->bIntermolecularInteractions
)
528 n
+= mtop
->intermolecular_ilist
[ftype
].nr
/(1+NRAL(ftype
));
534 t_block
gmx_mtop_global_cgs(const gmx_mtop_t
*mtop
)
536 t_block cgs_gl
, *cgs_mol
;
538 gmx_molblock_t
*molb
;
540 /* In most cases this is too much, but we realloc at the end */
541 snew(cgs_gl
.index
, mtop
->natoms
+1);
545 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
547 molb
= &mtop
->molblock
[mb
];
548 cgs_mol
= &mtop
->moltype
[molb
->type
].cgs
;
549 for (mol
= 0; mol
< molb
->nmol
; mol
++)
551 for (cg
= 0; cg
< cgs_mol
->nr
; cg
++)
553 cgs_gl
.index
[cgs_gl
.nr
+1] =
554 cgs_gl
.index
[cgs_gl
.nr
] +
555 cgs_mol
->index
[cg
+1] - cgs_mol
->index
[cg
];
560 cgs_gl
.nalloc_index
= cgs_gl
.nr
+ 1;
561 srenew(cgs_gl
.index
, cgs_gl
.nalloc_index
);
566 static void atomcat(t_atoms
*dest
, t_atoms
*src
, int copies
,
567 int maxres_renum
, int *maxresnr
)
571 int destnr
= dest
->nr
;
575 dest
->haveMass
= src
->haveMass
;
576 dest
->haveType
= src
->haveType
;
577 dest
->haveCharge
= src
->haveCharge
;
578 dest
->haveBState
= src
->haveBState
;
579 dest
->havePdbInfo
= src
->havePdbInfo
;
583 dest
->haveMass
= dest
->haveMass
&& src
->haveMass
;
584 dest
->haveType
= dest
->haveType
&& src
->haveType
;
585 dest
->haveCharge
= dest
->haveCharge
&& src
->haveCharge
;
586 dest
->haveBState
= dest
->haveBState
&& src
->haveBState
;
587 dest
->havePdbInfo
= dest
->havePdbInfo
&& src
->havePdbInfo
;
592 size
= destnr
+copies
*srcnr
;
593 srenew(dest
->atom
, size
);
594 srenew(dest
->atomname
, size
);
597 srenew(dest
->atomtype
, size
);
598 if (dest
->haveBState
)
600 srenew(dest
->atomtypeB
, size
);
603 if (dest
->havePdbInfo
)
605 srenew(dest
->pdbinfo
, size
);
610 size
= dest
->nres
+copies
*src
->nres
;
611 srenew(dest
->resinfo
, size
);
614 /* residue information */
615 for (l
= dest
->nres
, j
= 0; (j
< copies
); j
++, l
+= src
->nres
)
617 memcpy((char *) &(dest
->resinfo
[l
]), (char *) &(src
->resinfo
[0]),
618 (size_t)(src
->nres
*sizeof(src
->resinfo
[0])));
621 for (l
= destnr
, j
= 0; (j
< copies
); j
++, l
+= srcnr
)
623 memcpy((char *) &(dest
->atom
[l
]), (char *) &(src
->atom
[0]),
624 (size_t)(srcnr
*sizeof(src
->atom
[0])));
625 memcpy((char *) &(dest
->atomname
[l
]), (char *) &(src
->atomname
[0]),
626 (size_t)(srcnr
*sizeof(src
->atomname
[0])));
629 memcpy((char *) &(dest
->atomtype
[l
]), (char *) &(src
->atomtype
[0]),
630 (size_t)(srcnr
*sizeof(src
->atomtype
[0])));
631 if (dest
->haveBState
)
633 memcpy((char *) &(dest
->atomtypeB
[l
]), (char *) &(src
->atomtypeB
[0]),
634 (size_t)(srcnr
*sizeof(src
->atomtypeB
[0])));
637 if (dest
->havePdbInfo
)
639 memcpy((char *) &(dest
->pdbinfo
[l
]), (char *) &(src
->pdbinfo
[0]),
640 (size_t)(srcnr
*sizeof(src
->pdbinfo
[0])));
644 /* Increment residue indices */
645 for (l
= destnr
, j
= 0; (j
< copies
); j
++)
647 for (i
= 0; (i
< srcnr
); i
++, l
++)
649 dest
->atom
[l
].resind
= dest
->nres
+j
*src
->nres
+src
->atom
[i
].resind
;
653 if (src
->nres
<= maxres_renum
)
655 /* Single residue molecule, continue counting residues */
656 for (j
= 0; (j
< copies
); j
++)
658 for (l
= 0; l
< src
->nres
; l
++)
661 dest
->resinfo
[dest
->nres
+j
*src
->nres
+l
].nr
= *maxresnr
;
666 dest
->nres
+= copies
*src
->nres
;
667 dest
->nr
+= copies
*src
->nr
;
670 t_atoms
gmx_mtop_global_atoms(const gmx_mtop_t
*mtop
)
674 gmx_molblock_t
*molb
;
676 init_t_atoms(&atoms
, 0, FALSE
);
678 maxresnr
= mtop
->maxresnr
;
679 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
681 molb
= &mtop
->molblock
[mb
];
682 atomcat(&atoms
, &mtop
->moltype
[molb
->type
].atoms
, molb
->nmol
,
683 mtop
->maxres_renum
, &maxresnr
);
690 * The cat routines below are old code from src/kernel/topcat.c
693 static void blockcat(t_block
*dest
, t_block
*src
, int copies
)
695 int i
, j
, l
, nra
, size
;
699 size
= (dest
->nr
+copies
*src
->nr
+1);
700 srenew(dest
->index
, size
);
703 nra
= dest
->index
[dest
->nr
];
704 for (l
= dest
->nr
, j
= 0; (j
< copies
); j
++)
706 for (i
= 0; (i
< src
->nr
); i
++)
708 dest
->index
[l
++] = nra
+ src
->index
[i
];
710 nra
+= src
->index
[src
->nr
];
712 dest
->nr
+= copies
*src
->nr
;
713 dest
->index
[dest
->nr
] = nra
;
716 static void blockacat(t_blocka
*dest
, t_blocka
*src
, int copies
,
720 int destnr
= dest
->nr
;
721 int destnra
= dest
->nra
;
725 size
= (dest
->nr
+copies
*src
->nr
+1);
726 srenew(dest
->index
, size
);
730 size
= (dest
->nra
+copies
*src
->nra
);
731 srenew(dest
->a
, size
);
734 for (l
= destnr
, j
= 0; (j
< copies
); j
++)
736 for (i
= 0; (i
< src
->nr
); i
++)
738 dest
->index
[l
++] = dest
->nra
+src
->index
[i
];
740 dest
->nra
+= src
->nra
;
742 for (l
= destnra
, j
= 0; (j
< copies
); j
++)
744 for (i
= 0; (i
< src
->nra
); i
++)
746 dest
->a
[l
++] = dnum
+src
->a
[i
];
751 dest
->index
[dest
->nr
] = dest
->nra
;
754 static void ilistcat(int ftype
, t_ilist
*dest
, t_ilist
*src
, int copies
,
761 dest
->nalloc
= dest
->nr
+ copies
*src
->nr
;
762 srenew(dest
->iatoms
, dest
->nalloc
);
764 for (c
= 0; c
< copies
; c
++)
766 for (i
= 0; i
< src
->nr
; )
768 dest
->iatoms
[dest
->nr
++] = src
->iatoms
[i
++];
769 for (a
= 0; a
< nral
; a
++)
771 dest
->iatoms
[dest
->nr
++] = dnum
+ src
->iatoms
[i
++];
778 static void set_posres_params(t_idef
*idef
, gmx_molblock_t
*molb
,
779 int i0
, int a_offset
)
785 il
= &idef
->il
[F_POSRES
];
787 idef
->iparams_posres_nalloc
= i1
;
788 srenew(idef
->iparams_posres
, idef
->iparams_posres_nalloc
);
789 for (i
= i0
; i
< i1
; i
++)
791 ip
= &idef
->iparams_posres
[i
];
792 /* Copy the force constants */
793 *ip
= idef
->iparams
[il
->iatoms
[i
*2]];
794 a_molb
= il
->iatoms
[i
*2+1] - a_offset
;
795 if (molb
->nposres_xA
== 0)
797 gmx_incons("Position restraint coordinates are missing");
799 ip
->posres
.pos0A
[XX
] = molb
->posres_xA
[a_molb
][XX
];
800 ip
->posres
.pos0A
[YY
] = molb
->posres_xA
[a_molb
][YY
];
801 ip
->posres
.pos0A
[ZZ
] = molb
->posres_xA
[a_molb
][ZZ
];
802 if (molb
->nposres_xB
> 0)
804 ip
->posres
.pos0B
[XX
] = molb
->posres_xB
[a_molb
][XX
];
805 ip
->posres
.pos0B
[YY
] = molb
->posres_xB
[a_molb
][YY
];
806 ip
->posres
.pos0B
[ZZ
] = molb
->posres_xB
[a_molb
][ZZ
];
810 ip
->posres
.pos0B
[XX
] = ip
->posres
.pos0A
[XX
];
811 ip
->posres
.pos0B
[YY
] = ip
->posres
.pos0A
[YY
];
812 ip
->posres
.pos0B
[ZZ
] = ip
->posres
.pos0A
[ZZ
];
814 /* Set the parameter index for idef->iparams_posre */
819 static void set_fbposres_params(t_idef
*idef
, gmx_molblock_t
*molb
,
820 int i0
, int a_offset
)
826 il
= &idef
->il
[F_FBPOSRES
];
828 idef
->iparams_fbposres_nalloc
= i1
;
829 srenew(idef
->iparams_fbposres
, idef
->iparams_fbposres_nalloc
);
830 for (i
= i0
; i
< i1
; i
++)
832 ip
= &idef
->iparams_fbposres
[i
];
833 /* Copy the force constants */
834 *ip
= idef
->iparams
[il
->iatoms
[i
*2]];
835 a_molb
= il
->iatoms
[i
*2+1] - a_offset
;
836 if (molb
->nposres_xA
== 0)
838 gmx_incons("Position restraint coordinates are missing");
840 /* Take flat-bottom posres reference from normal position restraints */
841 ip
->fbposres
.pos0
[XX
] = molb
->posres_xA
[a_molb
][XX
];
842 ip
->fbposres
.pos0
[YY
] = molb
->posres_xA
[a_molb
][YY
];
843 ip
->fbposres
.pos0
[ZZ
] = molb
->posres_xA
[a_molb
][ZZ
];
844 /* Note: no B-type for flat-bottom posres */
846 /* Set the parameter index for idef->iparams_posre */
851 static void gen_local_top(const gmx_mtop_t
*mtop
,
852 bool freeEnergyInteractionsAtEnd
,
856 int mb
, srcnr
, destnr
, ftype
, natoms
, mol
, nposre_old
, nfbposre_old
;
857 gmx_molblock_t
*molb
;
859 const gmx_ffparams_t
*ffp
;
862 gmx_mtop_atomloop_all_t aloop
;
865 top
->atomtypes
= mtop
->atomtypes
;
867 ffp
= &mtop
->ffparams
;
870 idef
->ntypes
= ffp
->ntypes
;
871 idef
->atnr
= ffp
->atnr
;
872 idef
->functype
= ffp
->functype
;
873 idef
->iparams
= ffp
->iparams
;
874 idef
->iparams_posres
= nullptr;
875 idef
->iparams_posres_nalloc
= 0;
876 idef
->iparams_fbposres
= nullptr;
877 idef
->iparams_fbposres_nalloc
= 0;
878 idef
->fudgeQQ
= ffp
->fudgeQQ
;
879 idef
->cmap_grid
= ffp
->cmap_grid
;
880 idef
->ilsort
= ilsortUNKNOWN
;
882 init_block(&top
->cgs
);
883 init_blocka(&top
->excls
);
884 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
886 idef
->il
[ftype
].nr
= 0;
887 idef
->il
[ftype
].nalloc
= 0;
888 idef
->il
[ftype
].iatoms
= nullptr;
892 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
894 molb
= &mtop
->molblock
[mb
];
895 molt
= &mtop
->moltype
[molb
->type
];
897 srcnr
= molt
->atoms
.nr
;
900 blockcat(&top
->cgs
, &molt
->cgs
, molb
->nmol
);
902 blockacat(&top
->excls
, &molt
->excls
, molb
->nmol
, destnr
, srcnr
);
904 nposre_old
= idef
->il
[F_POSRES
].nr
;
905 nfbposre_old
= idef
->il
[F_FBPOSRES
].nr
;
906 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
909 ftype
== F_CONSTR
&& molt
->ilist
[F_CONSTRNC
].nr
> 0)
911 /* Merge all constrains into one ilist.
912 * This simplifies the constraint code.
914 for (mol
= 0; mol
< molb
->nmol
; mol
++)
916 ilistcat(ftype
, &idef
->il
[F_CONSTR
], &molt
->ilist
[F_CONSTR
],
917 1, destnr
+mol
*srcnr
, srcnr
);
918 ilistcat(ftype
, &idef
->il
[F_CONSTR
], &molt
->ilist
[F_CONSTRNC
],
919 1, destnr
+mol
*srcnr
, srcnr
);
922 else if (!(bMergeConstr
&& ftype
== F_CONSTRNC
))
924 ilistcat(ftype
, &idef
->il
[ftype
], &molt
->ilist
[ftype
],
925 molb
->nmol
, destnr
, srcnr
);
928 if (idef
->il
[F_POSRES
].nr
> nposre_old
)
930 /* Executing this line line stops gmxdump -sys working
931 * correctly. I'm not aware there's an elegant fix. */
932 set_posres_params(idef
, molb
, nposre_old
/2, natoms
);
934 if (idef
->il
[F_FBPOSRES
].nr
> nfbposre_old
)
936 set_fbposres_params(idef
, molb
, nfbposre_old
/2, natoms
);
939 natoms
+= molb
->nmol
*srcnr
;
942 if (mtop
->bIntermolecularInteractions
)
944 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
946 ilistcat(ftype
, &idef
->il
[ftype
], &mtop
->intermolecular_ilist
[ftype
],
951 if (freeEnergyInteractionsAtEnd
&& gmx_mtop_bondeds_free_energy(mtop
))
953 snew(qA
, mtop
->natoms
);
954 snew(qB
, mtop
->natoms
);
955 aloop
= gmx_mtop_atomloop_all_init(mtop
);
957 while (gmx_mtop_atomloop_all_next(aloop
, &ag
, &atom
))
962 gmx_sort_ilist_fe(&top
->idef
, qA
, qB
);
968 top
->idef
.ilsort
= ilsortNO_FE
;
973 gmx_mtop_generate_local_top(const gmx_mtop_t
*mtop
,
974 bool freeEnergyInteractionsAtEnd
)
980 gen_local_top(mtop
, freeEnergyInteractionsAtEnd
, true, top
);
985 /*! \brief Fills an array with molecule begin/end atom indices
987 * \param[in] mtop The global topology
988 * \param[out] index Array of size nr. of molecules + 1 to be filled with molecule begin/end indices
990 static void fillMoleculeIndices(const gmx_mtop_t
&mtop
,
991 gmx::ArrayRef
<int> index
)
993 int globalAtomIndex
= 0;
994 int globalMolIndex
= 0;
995 index
[globalMolIndex
] = globalAtomIndex
;
996 for (int mb
= 0; mb
< mtop
.nmolblock
; mb
++)
998 const gmx_molblock_t
&molb
= mtop
.molblock
[mb
];
999 for (int mol
= 0; mol
< molb
.nmol
; mol
++)
1001 globalAtomIndex
+= molb
.natoms_mol
;
1002 globalMolIndex
+= 1;
1003 index
[globalMolIndex
] = globalAtomIndex
;
1008 gmx::BlockRanges
gmx_mtop_molecules(const gmx_mtop_t
&mtop
)
1010 gmx::BlockRanges mols
;
1012 mols
.index
.resize(gmx_mtop_num_molecules(mtop
) + 1);
1014 fillMoleculeIndices(mtop
, mols
.index
);
1019 /*! \brief Creates and returns a deprecated t_block struct with molecule indices
1021 * \param[in] mtop The global topology
1023 static t_block
gmx_mtop_molecules_t_block(const gmx_mtop_t
&mtop
)
1027 mols
.nr
= gmx_mtop_num_molecules(mtop
);
1028 mols
.nalloc_index
= mols
.nr
+ 1;
1029 snew(mols
.index
, mols
.nalloc_index
);
1031 fillMoleculeIndices(mtop
, gmx::arrayRefFromArray(mols
.index
, mols
.nr
+ 1));
1036 t_topology
gmx_mtop_t_to_t_topology(gmx_mtop_t
*mtop
, bool freeMTop
)
1039 gmx_localtop_t ltop
;
1042 gen_local_top(mtop
, false, FALSE
, <op
);
1043 ltop
.idef
.ilsort
= ilsortUNKNOWN
;
1045 top
.name
= mtop
->name
;
1046 top
.idef
= ltop
.idef
;
1047 top
.atomtypes
= ltop
.atomtypes
;
1049 top
.excls
= ltop
.excls
;
1050 top
.atoms
= gmx_mtop_global_atoms(mtop
);
1051 top
.mols
= gmx_mtop_molecules_t_block(*mtop
);
1052 top
.bIntermolecularInteractions
= mtop
->bIntermolecularInteractions
;
1053 top
.symtab
= mtop
->symtab
;
1057 // Free pointers that have not been copied to top.
1058 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
1060 done_moltype(&mtop
->moltype
[mt
]);
1062 sfree(mtop
->moltype
);
1064 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
1066 done_molblock(&mtop
->molblock
[mb
]);
1068 sfree(mtop
->molblock
);
1070 done_gmx_groups_t(&mtop
->groups
);
1076 std::vector
<size_t> get_atom_index(const gmx_mtop_t
*mtop
)
1079 std::vector
<size_t> atom_index
;
1080 gmx_mtop_atomloop_block_t aloopb
= gmx_mtop_atomloop_block_init(mtop
);
1083 while (gmx_mtop_atomloop_block_next(aloopb
, &atom
, &nmol
))
1085 if (atom
->ptype
== eptAtom
)
1087 atom_index
.push_back(j
);
1094 void convertAtomsToMtop(t_symtab
*symtab
,
1099 mtop
->symtab
= *symtab
;
1104 // This snew clears all entries, we should replace it by an initializer
1105 snew(mtop
->moltype
, mtop
->nmoltype
);
1106 mtop
->moltype
[0].atoms
= *atoms
;
1107 init_block(&mtop
->moltype
[0].cgs
);
1108 init_blocka(&mtop
->moltype
[0].excls
);
1110 mtop
->nmolblock
= 1;
1111 // This snew clears all entries, we should replace it by an initializer
1112 snew(mtop
->molblock
, mtop
->nmolblock
);
1113 mtop
->molblock
[0].type
= 0;
1114 mtop
->molblock
[0].nmol
= 1;
1115 mtop
->molblock
[0].natoms_mol
= atoms
->nr
;
1117 mtop
->bIntermolecularInteractions
= FALSE
;
1119 mtop
->natoms
= atoms
->nr
;
1121 mtop
->haveMoleculeIndices
= false;
1123 gmx_mtop_finalize(mtop
);