1 /* -*- mode: c; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4; c-file-style: "stroustrup"; -*-
4 * This file is part of Gromacs Copyright (c) 1991-2008
5 * David van der Spoel, Erik Lindahl, Berk Hess, University of Groningen.
7 * This program is free software; you can redistribute it and/or
8 * modify it under the terms of the GNU General Public License
9 * as published by the Free Software Foundation; either version 2
10 * of the License, or (at your option) any later version.
12 * To help us fund GROMACS development, we humbly ask that you cite
13 * the research papers on the package. Check out http://www.gromacs.org
16 * Gnomes, ROck Monsters And Chili Sauce
26 #include "mtop_util.h"
29 #include "gmx_fatal.h"
31 static int gmx_mtop_maxresnr(const gmx_mtop_t *mtop, int maxres_renum)
38 for (mt = 0; mt < mtop->nmoltype; mt++)
40 atoms = &mtop->moltype[mt].atoms;
41 if (atoms->nres > maxres_renum)
43 for (r = 0; r < atoms->nres; r++)
45 if (atoms->resinfo[r].nr > maxresnr)
47 maxresnr = atoms->resinfo[r].nr;
56 void gmx_mtop_finalize(gmx_mtop_t *mtop)
60 mtop->maxres_renum = 1;
62 env = getenv("GMX_MAXRESRENUM");
65 sscanf(env, "%d", &mtop->maxres_renum);
67 if (mtop->maxres_renum == -1)
69 /* -1 signals renumber residues in all molecules */
70 mtop->maxres_renum = INT_MAX;
73 mtop->maxresnr = gmx_mtop_maxresnr(mtop, mtop->maxres_renum);
76 int ncg_mtop(const gmx_mtop_t *mtop)
82 for (mb = 0; mb < mtop->nmolblock; mb++)
85 mtop->molblock[mb].nmol*
86 mtop->moltype[mtop->molblock[mb].type].cgs.nr;
92 void gmx_mtop_remove_chargegroups(gmx_mtop_t *mtop)
98 for (mt = 0; mt < mtop->nmoltype; mt++)
100 cgs = &mtop->moltype[mt].cgs;
101 if (cgs->nr < mtop->moltype[mt].atoms.nr)
103 cgs->nr = mtop->moltype[mt].atoms.nr;
104 srenew(cgs->index, cgs->nr+1);
105 for (i = 0; i < cgs->nr+1; i++)
121 typedef struct gmx_mtop_atomlookup
123 const gmx_mtop_t *mtop;
127 } t_gmx_mtop_atomlookup;
130 gmx_mtop_atomlookup_t
131 gmx_mtop_atomlookup_init(const gmx_mtop_t *mtop)
133 t_gmx_mtop_atomlookup *alook;
135 int a_start, a_end, na, na_start = -1;
140 alook->nmb = mtop->nmolblock;
142 snew(alook->mba, alook->nmb);
145 for (mb = 0; mb < mtop->nmolblock; mb++)
147 na = mtop->molblock[mb].nmol*mtop->molblock[mb].natoms_mol;
148 a_end = a_start + na;
150 alook->mba[mb].a_start = a_start;
151 alook->mba[mb].a_end = a_end;
152 alook->mba[mb].na_mol = mtop->molblock[mb].natoms_mol;
154 /* We start the binary search with the largest block */
155 if (mb == 0 || na > na_start)
157 alook->mb_start = mb;
167 gmx_mtop_atomlookup_t
168 gmx_mtop_atomlookup_settle_init(const gmx_mtop_t *mtop)
170 t_gmx_mtop_atomlookup *alook;
172 int na, na_start = -1;
174 alook = gmx_mtop_atomlookup_init(mtop);
176 /* Check if the starting molblock has settle */
177 if (mtop->moltype[mtop->molblock[alook->mb_start].type].ilist[F_SETTLE].nr == 0)
179 /* Search the largest molblock with settle */
180 alook->mb_start = -1;
181 for (mb = 0; mb < mtop->nmolblock; mb++)
183 if (mtop->moltype[mtop->molblock[mb].type].ilist[F_SETTLE].nr > 0)
185 na = alook->mba[mb].a_end - alook->mba[mb].a_start;
186 if (alook->mb_start == -1 || na > na_start)
188 alook->mb_start = mb;
194 if (alook->mb_start == -1)
196 gmx_incons("gmx_mtop_atomlookup_settle_init called without settles");
204 gmx_mtop_atomlookup_destroy(gmx_mtop_atomlookup_t alook)
210 void gmx_mtop_atomnr_to_atom(const gmx_mtop_atomlookup_t alook,
215 int a_start, atnr_mol;
218 if (atnr_global < 0 || atnr_global >= mtop->natoms)
220 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)",
221 atnr_global, 0, mtop->natoms-1);
227 mb = alook->mb_start;
231 a_start = alook->mba[mb].a_start;
232 if (atnr_global < a_start)
236 else if (atnr_global >= alook->mba[mb].a_end)
244 mb = ((mb0 + mb1 + 1)>>1);
247 atnr_mol = (atnr_global - a_start) % alook->mba[mb].na_mol;
249 *atom = &alook->mtop->moltype[alook->mtop->molblock[mb].type].atoms.atom[atnr_mol];
252 void gmx_mtop_atomnr_to_ilist(const gmx_mtop_atomlookup_t alook,
254 t_ilist **ilist_mol, int *atnr_offset)
257 int a_start, atnr_local;
260 if (atnr_global < 0 || atnr_global >= mtop->natoms)
262 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)",
263 atnr_global, 0, mtop->natoms-1);
269 mb = alook->mb_start;
273 a_start = alook->mba[mb].a_start;
274 if (atnr_global < a_start)
278 else if (atnr_global >= alook->mba[mb].a_end)
286 mb = ((mb0 + mb1 + 1)>>1);
289 *ilist_mol = alook->mtop->moltype[alook->mtop->molblock[mb].type].ilist;
291 atnr_local = (atnr_global - a_start) % alook->mba[mb].na_mol;
293 *atnr_offset = atnr_global - atnr_local;
296 void gmx_mtop_atomnr_to_molblock_ind(const gmx_mtop_atomlookup_t alook,
298 int *molb, int *molnr, int *atnr_mol)
304 if (atnr_global < 0 || atnr_global >= mtop->natoms)
306 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)",
307 atnr_global, 0, mtop->natoms-1);
313 mb = alook->mb_start;
317 a_start = alook->mba[mb].a_start;
318 if (atnr_global < a_start)
322 else if (atnr_global >= alook->mba[mb].a_end)
330 mb = ((mb0 + mb1 + 1)>>1);
334 *molnr = (atnr_global - a_start) / alook->mba[mb].na_mol;
335 *atnr_mol = atnr_global - a_start - (*molnr)*alook->mba[mb].na_mol;
338 void gmx_mtop_atominfo_global(const gmx_mtop_t *mtop, int atnr_global,
339 char **atomname, int *resnr, char **resname)
341 int mb, a_start, a_end, maxresnr, at_loc;
342 gmx_molblock_t *molb;
343 t_atoms *atoms = NULL;
345 if (atnr_global < 0 || atnr_global >= mtop->natoms)
347 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)",
348 atnr_global, 0, mtop->natoms-1);
353 maxresnr = mtop->maxresnr;
358 if (atoms->nres <= mtop->maxres_renum)
360 /* Single residue molecule, keep counting */
361 maxresnr += mtop->molblock[mb].nmol*atoms->nres;
365 atoms = &mtop->moltype[mtop->molblock[mb].type].atoms;
367 a_end = a_start + mtop->molblock[mb].nmol*atoms->nr;
369 while (atnr_global >= a_end);
371 at_loc = (atnr_global - a_start) % atoms->nr;
372 *atomname = *(atoms->atomname[at_loc]);
373 if (atoms->nres > mtop->maxres_renum)
375 *resnr = atoms->resinfo[atoms->atom[at_loc].resind].nr;
379 /* Single residue molecule, keep counting */
380 *resnr = maxresnr + 1 + (atnr_global - a_start)/atoms->nr*atoms->nres + atoms->atom[at_loc].resind;
382 *resname = *(atoms->resinfo[atoms->atom[at_loc].resind].name);
385 typedef struct gmx_mtop_atomloop_all
387 const gmx_mtop_t *mtop;
394 } t_gmx_mtop_atomloop_all;
396 gmx_mtop_atomloop_all_t
397 gmx_mtop_atomloop_all_init(const gmx_mtop_t *mtop)
399 struct gmx_mtop_atomloop_all *aloop;
406 &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
408 aloop->maxresnr = mtop->maxresnr;
409 aloop->at_local = -1;
410 aloop->at_global = -1;
415 static void gmx_mtop_atomloop_all_destroy(gmx_mtop_atomloop_all_t aloop)
420 gmx_bool gmx_mtop_atomloop_all_next(gmx_mtop_atomloop_all_t aloop,
421 int *at_global, t_atom **atom)
425 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
431 if (aloop->at_local >= aloop->atoms->nr)
433 if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
435 /* Single residue molecule, increase the count with one */
436 aloop->maxresnr += aloop->atoms->nres;
440 if (aloop->mol >= aloop->mtop->molblock[aloop->mblock].nmol)
443 if (aloop->mblock >= aloop->mtop->nmolblock)
445 gmx_mtop_atomloop_all_destroy(aloop);
448 aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
453 *at_global = aloop->at_global;
454 *atom = &aloop->atoms->atom[aloop->at_local];
459 void gmx_mtop_atomloop_all_names(gmx_mtop_atomloop_all_t aloop,
460 char **atomname, int *resnr, char **resname)
464 *atomname = *(aloop->atoms->atomname[aloop->at_local]);
465 resind_mol = aloop->atoms->atom[aloop->at_local].resind;
466 *resnr = aloop->atoms->resinfo[resind_mol].nr;
467 if (aloop->atoms->nres <= aloop->mtop->maxres_renum)
469 *resnr = aloop->maxresnr + 1 + resind_mol;
471 *resname = *(aloop->atoms->resinfo[resind_mol].name);
474 void gmx_mtop_atomloop_all_moltype(gmx_mtop_atomloop_all_t aloop,
475 gmx_moltype_t **moltype, int *at_mol)
477 *moltype = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type];
478 *at_mol = aloop->at_local;
481 typedef struct gmx_mtop_atomloop_block
483 const gmx_mtop_t *mtop;
487 } t_gmx_mtop_atomloop_block;
489 gmx_mtop_atomloop_block_t
490 gmx_mtop_atomloop_block_init(const gmx_mtop_t *mtop)
492 struct gmx_mtop_atomloop_block *aloop;
498 aloop->atoms = &mtop->moltype[mtop->molblock[aloop->mblock].type].atoms;
499 aloop->at_local = -1;
504 static void gmx_mtop_atomloop_block_destroy(gmx_mtop_atomloop_block_t aloop)
509 gmx_bool gmx_mtop_atomloop_block_next(gmx_mtop_atomloop_block_t aloop,
510 t_atom **atom, int *nmol)
514 gmx_incons("gmx_mtop_atomloop_all_next called without calling gmx_mtop_atomloop_all_init");
519 if (aloop->at_local >= aloop->atoms->nr)
522 if (aloop->mblock >= aloop->mtop->nmolblock)
524 gmx_mtop_atomloop_block_destroy(aloop);
527 aloop->atoms = &aloop->mtop->moltype[aloop->mtop->molblock[aloop->mblock].type].atoms;
531 *atom = &aloop->atoms->atom[aloop->at_local];
532 *nmol = aloop->mtop->molblock[aloop->mblock].nmol;
537 typedef struct gmx_mtop_ilistloop
539 const gmx_mtop_t *mtop;
544 gmx_mtop_ilistloop_init(const gmx_mtop_t *mtop)
546 struct gmx_mtop_ilistloop *iloop;
556 static void gmx_mtop_ilistloop_destroy(gmx_mtop_ilistloop_t iloop)
561 gmx_bool gmx_mtop_ilistloop_next(gmx_mtop_ilistloop_t iloop,
562 t_ilist **ilist_mol, int *nmol)
566 gmx_incons("gmx_mtop_ilistloop_next called without calling gmx_mtop_ilistloop_init");
570 if (iloop->mblock == iloop->mtop->nmolblock)
572 gmx_mtop_ilistloop_destroy(iloop);
577 iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
579 *nmol = iloop->mtop->molblock[iloop->mblock].nmol;
583 typedef struct gmx_mtop_ilistloop_all
585 const gmx_mtop_t *mtop;
589 } t_gmx_mtop_ilist_all;
591 gmx_mtop_ilistloop_all_t
592 gmx_mtop_ilistloop_all_init(const gmx_mtop_t *mtop)
594 struct gmx_mtop_ilistloop_all *iloop;
606 static void gmx_mtop_ilistloop_all_destroy(gmx_mtop_ilistloop_all_t iloop)
611 gmx_bool gmx_mtop_ilistloop_all_next(gmx_mtop_ilistloop_all_t iloop,
612 t_ilist **ilist_mol, int *atnr_offset)
614 gmx_molblock_t *molb;
618 gmx_incons("gmx_mtop_ilistloop_all_next called without calling gmx_mtop_ilistloop_all_init");
623 iloop->a_offset += iloop->mtop->molblock[iloop->mblock].natoms_mol;
628 if (iloop->mol >= iloop->mtop->molblock[iloop->mblock].nmol)
632 if (iloop->mblock == iloop->mtop->nmolblock)
634 gmx_mtop_ilistloop_all_destroy(iloop);
640 iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
642 *atnr_offset = iloop->a_offset;
647 int gmx_mtop_ftype_count(const gmx_mtop_t *mtop, int ftype)
649 gmx_mtop_ilistloop_t iloop;
655 iloop = gmx_mtop_ilistloop_init(mtop);
656 while (gmx_mtop_ilistloop_next(iloop, &il, &nmol))
658 n += nmol*il[ftype].nr/(1+NRAL(ftype));
664 t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
666 t_block cgs_gl, *cgs_mol;
668 gmx_molblock_t *molb;
671 /* In most cases this is too much, but we realloc at the end */
672 snew(cgs_gl.index, mtop->natoms+1);
676 for (mb = 0; mb < mtop->nmolblock; mb++)
678 molb = &mtop->molblock[mb];
679 cgs_mol = &mtop->moltype[molb->type].cgs;
680 for (mol = 0; mol < molb->nmol; mol++)
682 for (cg = 0; cg < cgs_mol->nr; cg++)
684 cgs_gl.index[cgs_gl.nr+1] =
685 cgs_gl.index[cgs_gl.nr] +
686 cgs_mol->index[cg+1] - cgs_mol->index[cg];
691 cgs_gl.nalloc_index = cgs_gl.nr + 1;
692 srenew(cgs_gl.index, cgs_gl.nalloc_index);
697 static void atomcat(t_atoms *dest, t_atoms *src, int copies,
698 int maxres_renum, int *maxresnr)
702 int destnr = dest->nr;
706 size = destnr+copies*srcnr;
707 srenew(dest->atom, size);
708 srenew(dest->atomname, size);
709 srenew(dest->atomtype, size);
710 srenew(dest->atomtypeB, size);
714 size = dest->nres+copies*src->nres;
715 srenew(dest->resinfo, size);
718 /* residue information */
719 for (l = dest->nres, j = 0; (j < copies); j++, l += src->nres)
721 memcpy((char *) &(dest->resinfo[l]), (char *) &(src->resinfo[0]),
722 (size_t)(src->nres*sizeof(src->resinfo[0])));
725 for (l = destnr, j = 0; (j < copies); j++, l += srcnr)
727 memcpy((char *) &(dest->atomname[l]), (char *) &(src->atomname[0]),
728 (size_t)(srcnr*sizeof(src->atomname[0])));
729 memcpy((char *) &(dest->atomtype[l]), (char *) &(src->atomtype[0]),
730 (size_t)(srcnr*sizeof(src->atomtype[0])));
731 memcpy((char *) &(dest->atomtypeB[l]), (char *) &(src->atomtypeB[0]),
732 (size_t)(srcnr*sizeof(src->atomtypeB[0])));
733 memcpy((char *) &(dest->atom[l]), (char *) &(src->atom[0]),
734 (size_t)(srcnr*sizeof(src->atom[0])));
737 /* Increment residue indices */
738 for (l = destnr, j = 0; (j < copies); j++)
740 for (i = 0; (i < srcnr); i++, l++)
742 dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;
746 if (src->nres <= maxres_renum)
748 /* Single residue molecule, continue counting residues */
749 for (j = 0; (j < copies); j++)
751 for (l = 0; l < src->nres; l++)
754 dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;
759 dest->nres += copies*src->nres;
760 dest->nr += copies*src->nr;
763 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
767 gmx_molblock_t *molb;
769 init_t_atoms(&atoms, 0, FALSE);
771 maxresnr = mtop->maxresnr;
772 for (mb = 0; mb < mtop->nmolblock; mb++)
774 molb = &mtop->molblock[mb];
775 atomcat(&atoms, &mtop->moltype[molb->type].atoms, molb->nmol,
776 mtop->maxres_renum, &maxresnr);
782 void gmx_mtop_make_atomic_charge_groups(gmx_mtop_t *mtop,
783 gmx_bool bKeepSingleMolCG)
788 for (mb = 0; mb < mtop->nmolblock; mb++)
790 cgs_mol = &mtop->moltype[mtop->molblock[mb].type].cgs;
791 if (!(bKeepSingleMolCG && cgs_mol->nr == 1))
793 cgs_mol->nr = mtop->molblock[mb].natoms_mol;
794 cgs_mol->nalloc_index = cgs_mol->nr + 1;
795 srenew(cgs_mol->index, cgs_mol->nalloc_index);
796 for (cg = 0; cg < cgs_mol->nr+1; cg++)
798 cgs_mol->index[cg] = cg;
805 * The cat routines below are old code from src/kernel/topcat.c
808 static void blockcat(t_block *dest, t_block *src, int copies,
811 int i, j, l, nra, size;
815 size = (dest->nr+copies*src->nr+1);
816 srenew(dest->index, size);
819 nra = dest->index[dest->nr];
820 for (l = dest->nr, j = 0; (j < copies); j++)
822 for (i = 0; (i < src->nr); i++)
824 dest->index[l++] = nra + src->index[i];
826 nra += src->index[src->nr];
828 dest->nr += copies*src->nr;
829 dest->index[dest->nr] = nra;
832 static void blockacat(t_blocka *dest, t_blocka *src, int copies,
836 int destnr = dest->nr;
837 int destnra = dest->nra;
841 size = (dest->nr+copies*src->nr+1);
842 srenew(dest->index, size);
846 size = (dest->nra+copies*src->nra);
847 srenew(dest->a, size);
850 for (l = destnr, j = 0; (j < copies); j++)
852 for (i = 0; (i < src->nr); i++)
854 dest->index[l++] = dest->nra+src->index[i];
856 dest->nra += src->nra;
858 for (l = destnra, j = 0; (j < copies); j++)
860 for (i = 0; (i < src->nra); i++)
862 dest->a[l++] = dnum+src->a[i];
867 dest->index[dest->nr] = dest->nra;
870 static void ilistcat(int ftype, t_ilist *dest, t_ilist *src, int copies,
877 dest->nalloc = dest->nr + copies*src->nr;
878 srenew(dest->iatoms, dest->nalloc);
880 for (c = 0; c < copies; c++)
882 for (i = 0; i < src->nr; )
884 dest->iatoms[dest->nr++] = src->iatoms[i++];
885 for (a = 0; a < nral; a++)
887 dest->iatoms[dest->nr++] = dnum + src->iatoms[i++];
894 static void set_posres_params(t_idef *idef, gmx_molblock_t *molb,
895 int i0, int a_offset)
901 il = &idef->il[F_POSRES];
903 idef->iparams_posres_nalloc = i1;
904 srenew(idef->iparams_posres, idef->iparams_posres_nalloc);
905 for (i = i0; i < i1; i++)
907 ip = &idef->iparams_posres[i];
908 /* Copy the force constants */
909 *ip = idef->iparams[il->iatoms[i*2]];
910 a_molb = il->iatoms[i*2+1] - a_offset;
911 if (molb->nposres_xA == 0)
913 gmx_incons("Position restraint coordinates are missing");
915 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
916 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
917 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
918 if (molb->nposres_xB > 0)
920 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
921 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
922 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
926 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
927 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
928 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
930 /* Set the parameter index for idef->iparams_posre */
935 static void set_fbposres_params(t_idef *idef, gmx_molblock_t *molb,
936 int i0, int a_offset)
942 il = &idef->il[F_FBPOSRES];
944 idef->iparams_fbposres_nalloc = i1;
945 srenew(idef->iparams_fbposres, idef->iparams_fbposres_nalloc);
946 for (i = i0; i < i1; i++)
948 ip = &idef->iparams_fbposres[i];
949 /* Copy the force constants */
950 *ip = idef->iparams[il->iatoms[i*2]];
951 a_molb = il->iatoms[i*2+1] - a_offset;
952 if (molb->nposres_xA == 0)
954 gmx_incons("Position restraint coordinates are missing");
956 /* Take flat-bottom posres reference from normal position restraints */
957 ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
958 ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
959 ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
960 /* Note: no B-type for flat-bottom posres */
962 /* Set the parameter index for idef->iparams_posre */
967 static void gen_local_top(const gmx_mtop_t *mtop, const t_inputrec *ir,
968 gmx_bool bMergeConstr,
971 int mb, srcnr, destnr, ftype, ftype_dest, mt, natoms, mol, nposre_old, nfbposre_old;
972 gmx_molblock_t *molb;
974 const gmx_ffparams_t *ffp;
977 gmx_mtop_atomloop_all_t aloop;
981 top->atomtypes = mtop->atomtypes;
983 ffp = &mtop->ffparams;
986 idef->ntypes = ffp->ntypes;
987 idef->atnr = ffp->atnr;
988 idef->functype = ffp->functype;
989 idef->iparams = ffp->iparams;
990 idef->iparams_posres = NULL;
991 idef->iparams_posres_nalloc = 0;
992 idef->iparams_fbposres = NULL;
993 idef->iparams_fbposres_nalloc = 0;
994 idef->fudgeQQ = ffp->fudgeQQ;
995 idef->cmap_grid = ffp->cmap_grid;
996 idef->ilsort = ilsortUNKNOWN;
998 init_block(&top->cgs);
999 init_blocka(&top->excls);
1000 for (ftype = 0; ftype < F_NRE; ftype++)
1002 idef->il[ftype].nr = 0;
1003 idef->il[ftype].nalloc = 0;
1004 idef->il[ftype].iatoms = NULL;
1008 for (mb = 0; mb < mtop->nmolblock; mb++)
1010 molb = &mtop->molblock[mb];
1011 molt = &mtop->moltype[molb->type];
1013 srcnr = molt->atoms.nr;
1016 blockcat(&top->cgs, &molt->cgs, molb->nmol, destnr, srcnr);
1018 blockacat(&top->excls, &molt->excls, molb->nmol, destnr, srcnr);
1020 nposre_old = idef->il[F_POSRES].nr;
1021 nfbposre_old = idef->il[F_FBPOSRES].nr;
1022 for (ftype = 0; ftype < F_NRE; ftype++)
1025 ftype == F_CONSTR && molt->ilist[F_CONSTRNC].nr > 0)
1027 /* Merge all constrains into one ilist.
1028 * This simplifies the constraint code.
1030 for (mol = 0; mol < molb->nmol; mol++)
1032 ilistcat(ftype, &idef->il[F_CONSTR], &molt->ilist[F_CONSTR],
1033 1, destnr+mol*srcnr, srcnr);
1034 ilistcat(ftype, &idef->il[F_CONSTR], &molt->ilist[F_CONSTRNC],
1035 1, destnr+mol*srcnr, srcnr);
1038 else if (!(bMergeConstr && ftype == F_CONSTRNC))
1040 ilistcat(ftype, &idef->il[ftype], &molt->ilist[ftype],
1041 molb->nmol, destnr, srcnr);
1044 if (idef->il[F_POSRES].nr > nposre_old)
1046 /* Executing this line line stops gmxdump -sys working
1047 * correctly. I'm not aware there's an elegant fix. */
1048 set_posres_params(idef, molb, nposre_old/2, natoms);
1050 if (idef->il[F_FBPOSRES].nr > nfbposre_old)
1052 set_fbposres_params(idef, molb, nfbposre_old/2, natoms);
1055 natoms += molb->nmol*srcnr;
1060 top->idef.ilsort = ilsortUNKNOWN;
1064 if (ir->efep != efepNO && gmx_mtop_bondeds_free_energy(mtop))
1066 snew(qA, mtop->natoms);
1067 snew(qB, mtop->natoms);
1068 aloop = gmx_mtop_atomloop_all_init(mtop);
1069 while (gmx_mtop_atomloop_all_next(aloop, &ag, &atom))
1074 gmx_sort_ilist_fe(&top->idef, qA, qB);
1080 top->idef.ilsort = ilsortNO_FE;
1085 gmx_localtop_t *gmx_mtop_generate_local_top(const gmx_mtop_t *mtop,
1086 const t_inputrec *ir)
1088 gmx_localtop_t *top;
1092 gen_local_top(mtop, ir, TRUE, top);
1097 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop)
1100 gmx_localtop_t ltop;
1103 gen_local_top(mtop, NULL, FALSE, <op);
1105 top.name = mtop->name;
1106 top.idef = ltop.idef;
1107 top.atomtypes = ltop.atomtypes;
1109 top.excls = ltop.excls;
1110 top.atoms = gmx_mtop_global_atoms(mtop);
1111 top.mols = mtop->mols;
1112 top.symtab = mtop->symtab;
1114 /* We only need to free the moltype and molblock data,
1115 * all other pointers have been copied to top.
1117 * Well, except for the group data, but we can't free those, because they
1118 * are used somewhere even after a call to this function.
1120 for (mt = 0; mt < mtop->nmoltype; mt++)
1122 done_moltype(&mtop->moltype[mt]);
1124 sfree(mtop->moltype);
1126 for (mb = 0; mb < mtop->nmolblock; mb++)
1128 done_molblock(&mtop->molblock[mb]);
1130 sfree(mtop->molblock);