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
;
390 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t 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");
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
;
414 gmx_mtop_ilistloop_destroy(iloop
);
419 iloop
->mtop
->moltype
[iloop
->mtop
->molblock
[iloop
->mblock
].type
].ilist
;
421 *nmol
= iloop
->mtop
->molblock
[iloop
->mblock
].nmol
;
425 typedef struct gmx_mtop_ilistloop_all
427 const gmx_mtop_t
*mtop
;
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
;
448 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t 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");
464 iloop
->a_offset
+= iloop
->mtop
->moleculeBlockIndices
[iloop
->mblock
].numAtomsPerMolecule
;
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
)
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
;
488 gmx_mtop_ilistloop_all_destroy(iloop
);
494 iloop
->mtop
->moltype
[iloop
->mtop
->molblock
[iloop
->mblock
].type
].ilist
;
496 *atnr_offset
= iloop
->a_offset
;
501 int gmx_mtop_ftype_count(const gmx_mtop_t
*mtop
, int ftype
)
503 gmx_mtop_ilistloop_t iloop
;
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
));
523 t_block
gmx_mtop_global_cgs(const gmx_mtop_t
*mtop
)
526 /* In most cases this is too much, but we realloc at the end */
527 snew(cgs_gl
.index
, mtop
->natoms
+1);
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
];
545 cgs_gl
.nalloc_index
= cgs_gl
.nr
+ 1;
546 srenew(cgs_gl
.index
, cgs_gl
.nalloc_index
);
551 static void atomcat(t_atoms
*dest
, const t_atoms
*src
, int copies
,
552 int maxres_renum
, int *maxresnr
)
556 int destnr
= dest
->nr
;
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
;
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
;
577 size
= destnr
+copies
*srcnr
;
578 srenew(dest
->atom
, size
);
579 srenew(dest
->atomname
, size
);
582 srenew(dest
->atomtype
, size
);
583 if (dest
->haveBState
)
585 srenew(dest
->atomtypeB
, size
);
588 if (dest
->havePdbInfo
)
590 srenew(dest
->pdbinfo
, size
);
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])));
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
++)
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
)
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
);
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
;
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
];
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
,
705 int destnr
= dest
->nr
;
706 int destnra
= dest
->nra
;
710 size
= (dest
->nr
+copies
*src
->nr
+1);
711 srenew(dest
->index
, size
);
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
];
736 dest
->index
[dest
->nr
] = dest
->nra
;
739 static void ilistcat(int ftype
, t_ilist
*dest
, const t_ilist
*src
, int copies
,
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
++];
763 static void set_posres_params(t_idef
*idef
, const gmx_molblock_t
*molb
,
764 int i0
, int a_offset
)
770 il
= &idef
->il
[F_POSRES
];
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
];
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 */
804 static void set_fbposres_params(t_idef
*idef
, const gmx_molblock_t
*molb
,
805 int i0
, int a_offset
)
811 il
= &idef
->il
[F_FBPOSRES
];
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 */
836 static void gen_local_top(const gmx_mtop_t
*mtop
,
837 bool freeEnergyInteractionsAtEnd
,
841 int srcnr
, destnr
, ftype
, natoms
, nposre_old
, nfbposre_old
;
842 const gmx_ffparams_t
*ffp
;
845 gmx_mtop_atomloop_all_t aloop
;
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
);
858 top
->atomtypes
.atomnumber
= nullptr;
861 ffp
= &mtop
->ffparams
;
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!
873 snew(idef
->functype
, ffp
->ntypes
);
874 std::copy(ffp
->functype
, ffp
->functype
+ ffp
->ntypes
, idef
->functype
);
878 idef
->functype
= nullptr;
882 snew(idef
->iparams
, ffp
->ntypes
);
883 std::copy(ffp
->iparams
, ffp
->iparams
+ ffp
->ntypes
, idef
->iparams
);
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
);
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;
917 for (const gmx_molblock_t
&molb
: mtop
->molblock
)
919 const gmx_moltype_t
&molt
= mtop
->moltype
[molb
.type
];
921 srcnr
= molt
.atoms
.nr
;
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
++)
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
],
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
);
981 while (gmx_mtop_atomloop_all_next(aloop
, &ag
, &atom
))
986 gmx_sort_ilist_fe(&top
->idef
, qA
, qB
);
992 top
->idef
.ilsort
= ilsortNO_FE
;
997 gmx_mtop_generate_local_top(const gmx_mtop_t
*mtop
,
998 bool freeEnergyInteractionsAtEnd
)
1000 gmx_localtop_t
*top
;
1004 gen_local_top(mtop
, freeEnergyInteractionsAtEnd
, true, 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
);
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
)
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));
1060 t_topology
gmx_mtop_t_to_t_topology(gmx_mtop_t
*mtop
, bool freeMTop
)
1062 gmx_localtop_t ltop
;
1065 gen_local_top(mtop
, false, FALSE
, <op
);
1066 ltop
.idef
.ilsort
= ilsortUNKNOWN
;
1068 top
.name
= mtop
->name
;
1069 top
.idef
= ltop
.idef
;
1070 top
.atomtypes
= ltop
.atomtypes
;
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
;
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;
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
);
1096 while (gmx_mtop_atomloop_block_next(aloopb
, &atom
, &nmol
))
1098 if (atom
->ptype
== eptAtom
)
1100 atom_index
.push_back(j
);
1107 void convertAtomsToMtop(t_symtab
*symtab
,
1112 mtop
->symtab
= *symtab
;
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
);