2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 2008,2009,2010,2012,2013,2014,2015,2016, 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/fatalerror.h"
53 #include "gromacs/utility/real.h"
54 #include "gromacs/utility/smalloc.h"
56 static int gmx_mtop_maxresnr(const gmx_mtop_t
*mtop
, int maxres_renum
)
63 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
65 atoms
= &mtop
->moltype
[mt
].atoms
;
66 if (atoms
->nres
> maxres_renum
)
68 for (r
= 0; r
< atoms
->nres
; r
++)
70 if (atoms
->resinfo
[r
].nr
> maxresnr
)
72 maxresnr
= atoms
->resinfo
[r
].nr
;
81 void gmx_mtop_finalize(gmx_mtop_t
*mtop
)
85 if (mtop
->nmolblock
== 1 && mtop
->molblock
[0].nmol
== 1)
87 /* We have a single molecule only, no renumbering needed.
88 * This case also covers an mtop converted from pdb/gro/... input,
89 * so we retain the original residue numbering.
91 mtop
->maxres_renum
= 0;
95 /* We only renumber single residue molecules. Their intra-molecular
96 * residue numbering is anyhow irrelevant.
98 mtop
->maxres_renum
= 1;
101 env
= getenv("GMX_MAXRESRENUM");
104 sscanf(env
, "%d", &mtop
->maxres_renum
);
106 if (mtop
->maxres_renum
== -1)
108 /* -1 signals renumber residues in all molecules */
109 mtop
->maxres_renum
= INT_MAX
;
112 mtop
->maxresnr
= gmx_mtop_maxresnr(mtop
, mtop
->maxres_renum
);
115 void gmx_mtop_count_atomtypes(const gmx_mtop_t
*mtop
, int state
, int typecount
[])
117 int i
, mb
, nmol
, tpi
;
120 for (i
= 0; i
< mtop
->ffparams
.atnr
; ++i
)
124 for (mb
= 0; mb
< mtop
->nmolblock
; ++mb
)
126 nmol
= mtop
->molblock
[mb
].nmol
;
127 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
128 for (i
= 0; i
< atoms
->nr
; ++i
)
132 tpi
= atoms
->atom
[i
].type
;
136 tpi
= atoms
->atom
[i
].typeB
;
138 typecount
[tpi
] += nmol
;
143 int ncg_mtop(const gmx_mtop_t
*mtop
)
149 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
152 mtop
->molblock
[mb
].nmol
*
153 mtop
->moltype
[mtop
->molblock
[mb
].type
].cgs
.nr
;
159 void gmx_mtop_remove_chargegroups(gmx_mtop_t
*mtop
)
165 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
167 cgs
= &mtop
->moltype
[mt
].cgs
;
168 if (cgs
->nr
< mtop
->moltype
[mt
].atoms
.nr
)
170 cgs
->nr
= mtop
->moltype
[mt
].atoms
.nr
;
171 srenew(cgs
->index
, cgs
->nr
+1);
172 for (i
= 0; i
< cgs
->nr
+1; i
++)
188 typedef struct gmx_mtop_atomlookup
190 const gmx_mtop_t
*mtop
;
194 } t_gmx_mtop_atomlookup
;
197 gmx_mtop_atomlookup_t
198 gmx_mtop_atomlookup_init(const gmx_mtop_t
*mtop
)
200 t_gmx_mtop_atomlookup
*alook
;
202 int a_start
, a_end
, na
, na_start
= -1;
207 alook
->nmb
= mtop
->nmolblock
;
209 snew(alook
->mba
, alook
->nmb
);
212 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
214 na
= mtop
->molblock
[mb
].nmol
*mtop
->molblock
[mb
].natoms_mol
;
215 a_end
= a_start
+ na
;
217 alook
->mba
[mb
].a_start
= a_start
;
218 alook
->mba
[mb
].a_end
= a_end
;
219 alook
->mba
[mb
].na_mol
= mtop
->molblock
[mb
].natoms_mol
;
221 /* We start the binary search with the largest block */
222 if (mb
== 0 || na
> na_start
)
224 alook
->mb_start
= mb
;
234 gmx_mtop_atomlookup_t
235 gmx_mtop_atomlookup_settle_init(const gmx_mtop_t
*mtop
)
237 t_gmx_mtop_atomlookup
*alook
;
239 int na
, na_start
= -1;
241 alook
= gmx_mtop_atomlookup_init(mtop
);
243 /* Check if the starting molblock has settle */
244 if (mtop
->moltype
[mtop
->molblock
[alook
->mb_start
].type
].ilist
[F_SETTLE
].nr
== 0)
246 /* Search the largest molblock with settle */
247 alook
->mb_start
= -1;
248 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
250 if (mtop
->moltype
[mtop
->molblock
[mb
].type
].ilist
[F_SETTLE
].nr
> 0)
252 na
= alook
->mba
[mb
].a_end
- alook
->mba
[mb
].a_start
;
253 if (alook
->mb_start
== -1 || na
> na_start
)
255 alook
->mb_start
= mb
;
261 if (alook
->mb_start
== -1)
263 gmx_incons("gmx_mtop_atomlookup_settle_init called without settles");
271 gmx_mtop_atomlookup_destroy(gmx_mtop_atomlookup_t alook
)
277 void gmx_mtop_atomnr_to_atom(const gmx_mtop_atomlookup_t alook
,
282 int a_start
, atnr_mol
;
285 if (atnr_global
< 0 || atnr_global
>= mtop
->natoms
)
287 gmx_fatal(FARGS
, "gmx_mtop_atomnr_to_moltype was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
288 atnr_global
, 0, mtop
->natoms
-1);
294 mb
= alook
->mb_start
;
298 a_start
= alook
->mba
[mb
].a_start
;
299 if (atnr_global
< a_start
)
303 else if (atnr_global
>= alook
->mba
[mb
].a_end
)
311 mb
= ((mb0
+ mb1
+ 1)>>1);
314 atnr_mol
= (atnr_global
- a_start
) % alook
->mba
[mb
].na_mol
;
316 *atom
= &alook
->mtop
->moltype
[alook
->mtop
->molblock
[mb
].type
].atoms
.atom
[atnr_mol
];
319 void gmx_mtop_atomnr_to_ilist(const gmx_mtop_atomlookup_t alook
,
321 t_ilist
**ilist_mol
, int *atnr_offset
)
324 int a_start
, atnr_local
;
327 if (atnr_global
< 0 || atnr_global
>= mtop
->natoms
)
329 gmx_fatal(FARGS
, "gmx_mtop_atomnr_to_moltype was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
330 atnr_global
, 0, mtop
->natoms
-1);
336 mb
= alook
->mb_start
;
340 a_start
= alook
->mba
[mb
].a_start
;
341 if (atnr_global
< a_start
)
345 else if (atnr_global
>= alook
->mba
[mb
].a_end
)
353 mb
= ((mb0
+ mb1
+ 1)>>1);
356 *ilist_mol
= alook
->mtop
->moltype
[alook
->mtop
->molblock
[mb
].type
].ilist
;
358 atnr_local
= (atnr_global
- a_start
) % alook
->mba
[mb
].na_mol
;
360 *atnr_offset
= atnr_global
- atnr_local
;
363 void gmx_mtop_atomnr_to_molblock_ind(const gmx_mtop_atomlookup_t alook
,
365 int *molb
, int *molnr
, int *atnr_mol
)
371 if (atnr_global
< 0 || atnr_global
>= mtop
->natoms
)
373 gmx_fatal(FARGS
, "gmx_mtop_atomnr_to_moltype was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
374 atnr_global
, 0, mtop
->natoms
-1);
380 mb
= alook
->mb_start
;
384 a_start
= alook
->mba
[mb
].a_start
;
385 if (atnr_global
< a_start
)
389 else if (atnr_global
>= alook
->mba
[mb
].a_end
)
397 mb
= ((mb0
+ mb1
+ 1)>>1);
401 *molnr
= (atnr_global
- a_start
) / alook
->mba
[mb
].na_mol
;
402 *atnr_mol
= atnr_global
- a_start
- (*molnr
)*alook
->mba
[mb
].na_mol
;
405 void gmx_mtop_atominfo_global(const gmx_mtop_t
*mtop
, int atnr_global
,
406 char **atomname
, int *resnr
, char **resname
)
408 int mb
, a_start
, a_end
, maxresnr
, at_loc
;
409 t_atoms
*atoms
= NULL
;
411 if (atnr_global
< 0 || atnr_global
>= mtop
->natoms
)
413 gmx_fatal(FARGS
, "gmx_mtop_atominfo_global was called with atnr_global=%d which is not in the atom range of this system (%d-%d)",
414 atnr_global
, 0, mtop
->natoms
-1);
419 maxresnr
= mtop
->maxresnr
;
424 /* cppcheck-suppress nullPointer #6330 will be fixed in cppcheck 1.73 */
425 if (atoms
->nres
<= mtop
->maxres_renum
)
427 /* Single residue molecule, keep counting */
428 /* cppcheck-suppress nullPointer #6330 will be fixed in cppcheck 1.73 */
429 maxresnr
+= mtop
->molblock
[mb
].nmol
*atoms
->nres
;
433 atoms
= &mtop
->moltype
[mtop
->molblock
[mb
].type
].atoms
;
435 a_end
= a_start
+ mtop
->molblock
[mb
].nmol
*atoms
->nr
;
437 while (atnr_global
>= a_end
);
439 at_loc
= (atnr_global
- a_start
) % atoms
->nr
;
440 *atomname
= *(atoms
->atomname
[at_loc
]);
441 if (atoms
->nres
> mtop
->maxres_renum
)
443 *resnr
= atoms
->resinfo
[atoms
->atom
[at_loc
].resind
].nr
;
447 /* Single residue molecule, keep counting */
448 *resnr
= maxresnr
+ 1 + (atnr_global
- a_start
)/atoms
->nr
*atoms
->nres
+ atoms
->atom
[at_loc
].resind
;
450 *resname
= *(atoms
->resinfo
[atoms
->atom
[at_loc
].resind
].name
);
453 typedef struct gmx_mtop_atomloop_all
455 const gmx_mtop_t
*mtop
;
462 } t_gmx_mtop_atomloop_all
;
464 gmx_mtop_atomloop_all_t
465 gmx_mtop_atomloop_all_init(const gmx_mtop_t
*mtop
)
467 struct gmx_mtop_atomloop_all
*aloop
;
474 &mtop
->moltype
[mtop
->molblock
[aloop
->mblock
].type
].atoms
;
476 aloop
->maxresnr
= mtop
->maxresnr
;
477 aloop
->at_local
= -1;
478 aloop
->at_global
= -1;
483 static void gmx_mtop_atomloop_all_destroy(gmx_mtop_atomloop_all_t aloop
)
488 gmx_bool
gmx_mtop_atomloop_all_next(gmx_mtop_atomloop_all_t aloop
,
489 int *at_global
, t_atom
**atom
)
493 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
499 if (aloop
->at_local
>= aloop
->atoms
->nr
)
501 if (aloop
->atoms
->nres
<= aloop
->mtop
->maxres_renum
)
503 /* Single residue molecule, increase the count with one */
504 aloop
->maxresnr
+= aloop
->atoms
->nres
;
508 if (aloop
->mol
>= aloop
->mtop
->molblock
[aloop
->mblock
].nmol
)
511 if (aloop
->mblock
>= aloop
->mtop
->nmolblock
)
513 gmx_mtop_atomloop_all_destroy(aloop
);
516 aloop
->atoms
= &aloop
->mtop
->moltype
[aloop
->mtop
->molblock
[aloop
->mblock
].type
].atoms
;
521 *at_global
= aloop
->at_global
;
522 *atom
= &aloop
->atoms
->atom
[aloop
->at_local
];
527 void gmx_mtop_atomloop_all_names(gmx_mtop_atomloop_all_t aloop
,
528 char **atomname
, int *resnr
, char **resname
)
532 *atomname
= *(aloop
->atoms
->atomname
[aloop
->at_local
]);
533 resind_mol
= aloop
->atoms
->atom
[aloop
->at_local
].resind
;
534 *resnr
= aloop
->atoms
->resinfo
[resind_mol
].nr
;
535 if (aloop
->atoms
->nres
<= aloop
->mtop
->maxres_renum
)
537 *resnr
= aloop
->maxresnr
+ 1 + resind_mol
;
539 *resname
= *(aloop
->atoms
->resinfo
[resind_mol
].name
);
542 void gmx_mtop_atomloop_all_moltype(gmx_mtop_atomloop_all_t aloop
,
543 gmx_moltype_t
**moltype
, int *at_mol
)
545 *moltype
= &aloop
->mtop
->moltype
[aloop
->mtop
->molblock
[aloop
->mblock
].type
];
546 *at_mol
= aloop
->at_local
;
549 typedef struct gmx_mtop_atomloop_block
551 const gmx_mtop_t
*mtop
;
555 } t_gmx_mtop_atomloop_block
;
557 gmx_mtop_atomloop_block_t
558 gmx_mtop_atomloop_block_init(const gmx_mtop_t
*mtop
)
560 struct gmx_mtop_atomloop_block
*aloop
;
566 aloop
->atoms
= &mtop
->moltype
[mtop
->molblock
[aloop
->mblock
].type
].atoms
;
567 aloop
->at_local
= -1;
572 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop
)
577 gmx_bool
gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop
,
578 t_atom
**atom
, int *nmol
)
582 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
587 if (aloop
->at_local
>= aloop
->atoms
->nr
)
590 if (aloop
->mblock
>= aloop
->mtop
->nmolblock
)
592 gmx_mtop_atomloop_block_destroy(aloop
);
595 aloop
->atoms
= &aloop
->mtop
->moltype
[aloop
->mtop
->molblock
[aloop
->mblock
].type
].atoms
;
599 *atom
= &aloop
->atoms
->atom
[aloop
->at_local
];
600 *nmol
= aloop
->mtop
->molblock
[aloop
->mblock
].nmol
;
605 typedef struct gmx_mtop_ilistloop
607 const gmx_mtop_t
*mtop
;
612 gmx_mtop_ilistloop_init(const gmx_mtop_t
*mtop
)
614 struct gmx_mtop_ilistloop
*iloop
;
624 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop
)
629 gmx_bool
gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop
,
630 t_ilist
**ilist_mol
, int *nmol
)
634 gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
638 if (iloop
->mblock
>= iloop
->mtop
->nmolblock
)
640 if (iloop
->mblock
== iloop
->mtop
->nmolblock
&&
641 iloop
->mtop
->bIntermolecularInteractions
)
643 *ilist_mol
= iloop
->mtop
->intermolecular_ilist
;
648 gmx_mtop_ilistloop_destroy(iloop
);
653 iloop
->mtop
->moltype
[iloop
->mtop
->molblock
[iloop
->mblock
].type
].ilist
;
655 *nmol
= iloop
->mtop
->molblock
[iloop
->mblock
].nmol
;
659 typedef struct gmx_mtop_ilistloop_all
661 const gmx_mtop_t
*mtop
;
665 } t_gmx_mtop_ilist_all
;
667 gmx_mtop_ilistloop_all_t
668 gmx_mtop_ilistloop_all_init(const gmx_mtop_t
*mtop
)
670 struct gmx_mtop_ilistloop_all
*iloop
;
682 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop
)
687 gmx_bool
gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop
,
688 t_ilist
**ilist_mol
, int *atnr_offset
)
693 gmx_incons("gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
698 iloop
->a_offset
+= iloop
->mtop
->molblock
[iloop
->mblock
].natoms_mol
;
703 /* Inter-molecular interactions, if present, are indexed with
704 * iloop->mblock == iloop->mtop->nmolblock, thus we should separately
705 * check for this value in this conditional.
707 if (iloop
->mblock
== iloop
->mtop
->nmolblock
||
708 iloop
->mol
>= iloop
->mtop
->molblock
[iloop
->mblock
].nmol
)
712 if (iloop
->mblock
>= iloop
->mtop
->nmolblock
)
714 if (iloop
->mblock
== iloop
->mtop
->nmolblock
&&
715 iloop
->mtop
->bIntermolecularInteractions
)
717 *ilist_mol
= iloop
->mtop
->intermolecular_ilist
;
722 gmx_mtop_ilistloop_all_destroy(iloop
);
728 iloop
->mtop
->moltype
[iloop
->mtop
->molblock
[iloop
->mblock
].type
].ilist
;
730 *atnr_offset
= iloop
->a_offset
;
735 int gmx_mtop_ftype_count(const gmx_mtop_t
*mtop
, int ftype
)
737 gmx_mtop_ilistloop_t iloop
;
743 iloop
= gmx_mtop_ilistloop_init(mtop
);
744 while (gmx_mtop_ilistloop_next(iloop
, &il
, &nmol
))
746 n
+= nmol
*il
[ftype
].nr
/(1+NRAL(ftype
));
749 if (mtop
->bIntermolecularInteractions
)
751 n
+= mtop
->intermolecular_ilist
[ftype
].nr
/(1+NRAL(ftype
));
757 t_block
gmx_mtop_global_cgs(const gmx_mtop_t
*mtop
)
759 t_block cgs_gl
, *cgs_mol
;
761 gmx_molblock_t
*molb
;
763 /* In most cases this is too much, but we realloc at the end */
764 snew(cgs_gl
.index
, mtop
->natoms
+1);
768 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
770 molb
= &mtop
->molblock
[mb
];
771 cgs_mol
= &mtop
->moltype
[molb
->type
].cgs
;
772 for (mol
= 0; mol
< molb
->nmol
; mol
++)
774 for (cg
= 0; cg
< cgs_mol
->nr
; cg
++)
776 cgs_gl
.index
[cgs_gl
.nr
+1] =
777 cgs_gl
.index
[cgs_gl
.nr
] +
778 cgs_mol
->index
[cg
+1] - cgs_mol
->index
[cg
];
783 cgs_gl
.nalloc_index
= cgs_gl
.nr
+ 1;
784 srenew(cgs_gl
.index
, cgs_gl
.nalloc_index
);
789 static void atomcat(t_atoms
*dest
, t_atoms
*src
, int copies
,
790 int maxres_renum
, int *maxresnr
)
794 int destnr
= dest
->nr
;
798 dest
->haveMass
= src
->haveMass
;
799 dest
->haveType
= src
->haveType
;
800 dest
->haveCharge
= src
->haveCharge
;
801 dest
->haveBState
= src
->haveBState
;
802 dest
->havePdbInfo
= src
->havePdbInfo
;
806 dest
->haveMass
= dest
->haveMass
&& src
->haveMass
;
807 dest
->haveType
= dest
->haveType
&& src
->haveType
;
808 dest
->haveCharge
= dest
->haveCharge
&& src
->haveCharge
;
809 dest
->haveBState
= dest
->haveBState
&& src
->haveBState
;
810 dest
->havePdbInfo
= dest
->havePdbInfo
&& src
->havePdbInfo
;
815 size
= destnr
+copies
*srcnr
;
816 srenew(dest
->atom
, size
);
817 srenew(dest
->atomname
, size
);
820 srenew(dest
->atomtype
, size
);
821 if (dest
->haveBState
)
823 srenew(dest
->atomtypeB
, size
);
826 if (dest
->havePdbInfo
)
828 srenew(dest
->pdbinfo
, size
);
833 size
= dest
->nres
+copies
*src
->nres
;
834 srenew(dest
->resinfo
, size
);
837 /* residue information */
838 for (l
= dest
->nres
, j
= 0; (j
< copies
); j
++, l
+= src
->nres
)
840 memcpy((char *) &(dest
->resinfo
[l
]), (char *) &(src
->resinfo
[0]),
841 (size_t)(src
->nres
*sizeof(src
->resinfo
[0])));
844 for (l
= destnr
, j
= 0; (j
< copies
); j
++, l
+= srcnr
)
846 memcpy((char *) &(dest
->atom
[l
]), (char *) &(src
->atom
[0]),
847 (size_t)(srcnr
*sizeof(src
->atom
[0])));
848 memcpy((char *) &(dest
->atomname
[l
]), (char *) &(src
->atomname
[0]),
849 (size_t)(srcnr
*sizeof(src
->atomname
[0])));
852 memcpy((char *) &(dest
->atomtype
[l
]), (char *) &(src
->atomtype
[0]),
853 (size_t)(srcnr
*sizeof(src
->atomtype
[0])));
854 if (dest
->haveBState
)
856 memcpy((char *) &(dest
->atomtypeB
[l
]), (char *) &(src
->atomtypeB
[0]),
857 (size_t)(srcnr
*sizeof(src
->atomtypeB
[0])));
860 if (dest
->havePdbInfo
)
862 memcpy((char *) &(dest
->pdbinfo
[l
]), (char *) &(src
->pdbinfo
[0]),
863 (size_t)(srcnr
*sizeof(src
->pdbinfo
[0])));
867 /* Increment residue indices */
868 for (l
= destnr
, j
= 0; (j
< copies
); j
++)
870 for (i
= 0; (i
< srcnr
); i
++, l
++)
872 dest
->atom
[l
].resind
= dest
->nres
+j
*src
->nres
+src
->atom
[i
].resind
;
876 if (src
->nres
<= maxres_renum
)
878 /* Single residue molecule, continue counting residues */
879 for (j
= 0; (j
< copies
); j
++)
881 for (l
= 0; l
< src
->nres
; l
++)
884 dest
->resinfo
[dest
->nres
+j
*src
->nres
+l
].nr
= *maxresnr
;
889 dest
->nres
+= copies
*src
->nres
;
890 dest
->nr
+= copies
*src
->nr
;
893 t_atoms
gmx_mtop_global_atoms(const gmx_mtop_t
*mtop
)
897 gmx_molblock_t
*molb
;
899 init_t_atoms(&atoms
, 0, FALSE
);
901 maxresnr
= mtop
->maxresnr
;
902 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
904 molb
= &mtop
->molblock
[mb
];
905 atomcat(&atoms
, &mtop
->moltype
[molb
->type
].atoms
, molb
->nmol
,
906 mtop
->maxres_renum
, &maxresnr
);
913 * The cat routines below are old code from src/kernel/topcat.c
916 static void blockcat(t_block
*dest
, t_block
*src
, int copies
)
918 int i
, j
, l
, nra
, size
;
922 size
= (dest
->nr
+copies
*src
->nr
+1);
923 srenew(dest
->index
, size
);
926 nra
= dest
->index
[dest
->nr
];
927 for (l
= dest
->nr
, j
= 0; (j
< copies
); j
++)
929 for (i
= 0; (i
< src
->nr
); i
++)
931 dest
->index
[l
++] = nra
+ src
->index
[i
];
933 nra
+= src
->index
[src
->nr
];
935 dest
->nr
+= copies
*src
->nr
;
936 dest
->index
[dest
->nr
] = nra
;
939 static void blockacat(t_blocka
*dest
, t_blocka
*src
, int copies
,
943 int destnr
= dest
->nr
;
944 int destnra
= dest
->nra
;
948 size
= (dest
->nr
+copies
*src
->nr
+1);
949 srenew(dest
->index
, size
);
953 size
= (dest
->nra
+copies
*src
->nra
);
954 srenew(dest
->a
, size
);
957 for (l
= destnr
, j
= 0; (j
< copies
); j
++)
959 for (i
= 0; (i
< src
->nr
); i
++)
961 dest
->index
[l
++] = dest
->nra
+src
->index
[i
];
963 dest
->nra
+= src
->nra
;
965 for (l
= destnra
, j
= 0; (j
< copies
); j
++)
967 for (i
= 0; (i
< src
->nra
); i
++)
969 dest
->a
[l
++] = dnum
+src
->a
[i
];
974 dest
->index
[dest
->nr
] = dest
->nra
;
977 static void ilistcat(int ftype
, t_ilist
*dest
, t_ilist
*src
, int copies
,
984 dest
->nalloc
= dest
->nr
+ copies
*src
->nr
;
985 srenew(dest
->iatoms
, dest
->nalloc
);
987 for (c
= 0; c
< copies
; c
++)
989 for (i
= 0; i
< src
->nr
; )
991 dest
->iatoms
[dest
->nr
++] = src
->iatoms
[i
++];
992 for (a
= 0; a
< nral
; a
++)
994 dest
->iatoms
[dest
->nr
++] = dnum
+ src
->iatoms
[i
++];
1001 static void set_posres_params(t_idef
*idef
, gmx_molblock_t
*molb
,
1002 int i0
, int a_offset
)
1008 il
= &idef
->il
[F_POSRES
];
1010 idef
->iparams_posres_nalloc
= i1
;
1011 srenew(idef
->iparams_posres
, idef
->iparams_posres_nalloc
);
1012 for (i
= i0
; i
< i1
; i
++)
1014 ip
= &idef
->iparams_posres
[i
];
1015 /* Copy the force constants */
1016 *ip
= idef
->iparams
[il
->iatoms
[i
*2]];
1017 a_molb
= il
->iatoms
[i
*2+1] - a_offset
;
1018 if (molb
->nposres_xA
== 0)
1020 gmx_incons("Position restraint coordinates are missing");
1022 ip
->posres
.pos0A
[XX
] = molb
->posres_xA
[a_molb
][XX
];
1023 ip
->posres
.pos0A
[YY
] = molb
->posres_xA
[a_molb
][YY
];
1024 ip
->posres
.pos0A
[ZZ
] = molb
->posres_xA
[a_molb
][ZZ
];
1025 if (molb
->nposres_xB
> 0)
1027 ip
->posres
.pos0B
[XX
] = molb
->posres_xB
[a_molb
][XX
];
1028 ip
->posres
.pos0B
[YY
] = molb
->posres_xB
[a_molb
][YY
];
1029 ip
->posres
.pos0B
[ZZ
] = molb
->posres_xB
[a_molb
][ZZ
];
1033 ip
->posres
.pos0B
[XX
] = ip
->posres
.pos0A
[XX
];
1034 ip
->posres
.pos0B
[YY
] = ip
->posres
.pos0A
[YY
];
1035 ip
->posres
.pos0B
[ZZ
] = ip
->posres
.pos0A
[ZZ
];
1037 /* Set the parameter index for idef->iparams_posre */
1038 il
->iatoms
[i
*2] = i
;
1042 static void set_fbposres_params(t_idef
*idef
, gmx_molblock_t
*molb
,
1043 int i0
, int a_offset
)
1049 il
= &idef
->il
[F_FBPOSRES
];
1051 idef
->iparams_fbposres_nalloc
= i1
;
1052 srenew(idef
->iparams_fbposres
, idef
->iparams_fbposres_nalloc
);
1053 for (i
= i0
; i
< i1
; i
++)
1055 ip
= &idef
->iparams_fbposres
[i
];
1056 /* Copy the force constants */
1057 *ip
= idef
->iparams
[il
->iatoms
[i
*2]];
1058 a_molb
= il
->iatoms
[i
*2+1] - a_offset
;
1059 if (molb
->nposres_xA
== 0)
1061 gmx_incons("Position restraint coordinates are missing");
1063 /* Take flat-bottom posres reference from normal position restraints */
1064 ip
->fbposres
.pos0
[XX
] = molb
->posres_xA
[a_molb
][XX
];
1065 ip
->fbposres
.pos0
[YY
] = molb
->posres_xA
[a_molb
][YY
];
1066 ip
->fbposres
.pos0
[ZZ
] = molb
->posres_xA
[a_molb
][ZZ
];
1067 /* Note: no B-type for flat-bottom posres */
1069 /* Set the parameter index for idef->iparams_posre */
1070 il
->iatoms
[i
*2] = i
;
1074 static void gen_local_top(const gmx_mtop_t
*mtop
,
1075 bool freeEnergyInteractionsAtEnd
,
1077 gmx_localtop_t
*top
)
1079 int mb
, srcnr
, destnr
, ftype
, natoms
, mol
, nposre_old
, nfbposre_old
;
1080 gmx_molblock_t
*molb
;
1081 gmx_moltype_t
*molt
;
1082 const gmx_ffparams_t
*ffp
;
1085 gmx_mtop_atomloop_all_t aloop
;
1089 top
->atomtypes
= mtop
->atomtypes
;
1091 ffp
= &mtop
->ffparams
;
1094 idef
->ntypes
= ffp
->ntypes
;
1095 idef
->atnr
= ffp
->atnr
;
1096 idef
->functype
= ffp
->functype
;
1097 idef
->iparams
= ffp
->iparams
;
1098 idef
->iparams_posres
= NULL
;
1099 idef
->iparams_posres_nalloc
= 0;
1100 idef
->iparams_fbposres
= NULL
;
1101 idef
->iparams_fbposres_nalloc
= 0;
1102 idef
->fudgeQQ
= ffp
->fudgeQQ
;
1103 idef
->cmap_grid
= ffp
->cmap_grid
;
1104 idef
->ilsort
= ilsortUNKNOWN
;
1106 init_block(&top
->cgs
);
1107 init_blocka(&top
->excls
);
1108 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1110 idef
->il
[ftype
].nr
= 0;
1111 idef
->il
[ftype
].nalloc
= 0;
1112 idef
->il
[ftype
].iatoms
= NULL
;
1116 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
1118 molb
= &mtop
->molblock
[mb
];
1119 molt
= &mtop
->moltype
[molb
->type
];
1121 srcnr
= molt
->atoms
.nr
;
1124 blockcat(&top
->cgs
, &molt
->cgs
, molb
->nmol
);
1126 blockacat(&top
->excls
, &molt
->excls
, molb
->nmol
, destnr
, srcnr
);
1128 nposre_old
= idef
->il
[F_POSRES
].nr
;
1129 nfbposre_old
= idef
->il
[F_FBPOSRES
].nr
;
1130 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1133 ftype
== F_CONSTR
&& molt
->ilist
[F_CONSTRNC
].nr
> 0)
1135 /* Merge all constrains into one ilist.
1136 * This simplifies the constraint code.
1138 for (mol
= 0; mol
< molb
->nmol
; mol
++)
1140 ilistcat(ftype
, &idef
->il
[F_CONSTR
], &molt
->ilist
[F_CONSTR
],
1141 1, destnr
+mol
*srcnr
, srcnr
);
1142 ilistcat(ftype
, &idef
->il
[F_CONSTR
], &molt
->ilist
[F_CONSTRNC
],
1143 1, destnr
+mol
*srcnr
, srcnr
);
1146 else if (!(bMergeConstr
&& ftype
== F_CONSTRNC
))
1148 ilistcat(ftype
, &idef
->il
[ftype
], &molt
->ilist
[ftype
],
1149 molb
->nmol
, destnr
, srcnr
);
1152 if (idef
->il
[F_POSRES
].nr
> nposre_old
)
1154 /* Executing this line line stops gmxdump -sys working
1155 * correctly. I'm not aware there's an elegant fix. */
1156 set_posres_params(idef
, molb
, nposre_old
/2, natoms
);
1158 if (idef
->il
[F_FBPOSRES
].nr
> nfbposre_old
)
1160 set_fbposres_params(idef
, molb
, nfbposre_old
/2, natoms
);
1163 natoms
+= molb
->nmol
*srcnr
;
1166 if (mtop
->bIntermolecularInteractions
)
1168 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1170 ilistcat(ftype
, &idef
->il
[ftype
], &mtop
->intermolecular_ilist
[ftype
],
1171 1, 0, mtop
->natoms
);
1175 if (freeEnergyInteractionsAtEnd
&& gmx_mtop_bondeds_free_energy(mtop
))
1177 snew(qA
, mtop
->natoms
);
1178 snew(qB
, mtop
->natoms
);
1179 aloop
= gmx_mtop_atomloop_all_init(mtop
);
1180 while (gmx_mtop_atomloop_all_next(aloop
, &ag
, &atom
))
1185 gmx_sort_ilist_fe(&top
->idef
, qA
, qB
);
1191 top
->idef
.ilsort
= ilsortNO_FE
;
1196 gmx_mtop_generate_local_top(const gmx_mtop_t
*mtop
,
1197 bool freeEnergyInteractionsAtEnd
)
1199 gmx_localtop_t
*top
;
1203 gen_local_top(mtop
, freeEnergyInteractionsAtEnd
, true, top
);
1208 t_topology
gmx_mtop_t_to_t_topology(gmx_mtop_t
*mtop
)
1211 gmx_localtop_t ltop
;
1214 gen_local_top(mtop
, false, FALSE
, <op
);
1215 ltop
.idef
.ilsort
= ilsortUNKNOWN
;
1217 top
.name
= mtop
->name
;
1218 top
.idef
= ltop
.idef
;
1219 top
.atomtypes
= ltop
.atomtypes
;
1221 top
.excls
= ltop
.excls
;
1222 top
.atoms
= gmx_mtop_global_atoms(mtop
);
1223 top
.mols
= mtop
->mols
;
1224 top
.bIntermolecularInteractions
= mtop
->bIntermolecularInteractions
;
1225 top
.symtab
= mtop
->symtab
;
1227 /* We only need to free the moltype and molblock data,
1228 * all other pointers have been copied to top.
1230 * Well, except for the group data, but we can't free those, because they
1231 * are used somewhere even after a call to this function.
1233 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
1235 done_moltype(&mtop
->moltype
[mt
]);
1237 sfree(mtop
->moltype
);
1239 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
1241 done_molblock(&mtop
->molblock
[mb
]);
1243 sfree(mtop
->molblock
);
1248 std::vector
<size_t> get_atom_index(const gmx_mtop_t
*mtop
)
1251 std::vector
<size_t> atom_index
;
1252 gmx_mtop_atomloop_block_t aloopb
= gmx_mtop_atomloop_block_init(mtop
);
1255 while (gmx_mtop_atomloop_block_next(aloopb
, &atom
, &nmol
))
1257 if (atom
->ptype
== eptAtom
)
1259 atom_index
.push_back(j
);
1266 void convertAtomsToMtop(t_symtab
*symtab
,
1271 mtop
->symtab
= *symtab
;
1276 // This snew clears all entries, we should replace it by an initializer
1277 snew(mtop
->moltype
, mtop
->nmoltype
);
1278 mtop
->moltype
[0].atoms
= *atoms
;
1279 init_block(&mtop
->moltype
[0].cgs
);
1280 init_blocka(&mtop
->moltype
[0].excls
);
1282 mtop
->nmolblock
= 1;
1283 // This snew clears all entries, we should replace it by an initializer
1284 snew(mtop
->molblock
, mtop
->nmolblock
);
1285 mtop
->molblock
[0].type
= 0;
1286 mtop
->molblock
[0].nmol
= 1;
1287 mtop
->molblock
[0].natoms_mol
= atoms
->nr
;
1289 mtop
->bIntermolecularInteractions
= FALSE
;
1291 mtop
->natoms
= atoms
->nr
;
1293 gmx_mtop_finalize(mtop
);