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;
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;
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) {
631 if (iloop->mblock == iloop->mtop->nmolblock)
633 gmx_mtop_ilistloop_all_destroy(iloop);
639 iloop->mtop->moltype[iloop->mtop->molblock[iloop->mblock].type].ilist;
641 *atnr_offset = iloop->a_offset;
646 int gmx_mtop_ftype_count(const gmx_mtop_t *mtop,int ftype)
648 gmx_mtop_ilistloop_t iloop;
654 iloop = gmx_mtop_ilistloop_init(mtop);
655 while (gmx_mtop_ilistloop_next(iloop,&il,&nmol))
657 n += nmol*il[ftype].nr/(1+NRAL(ftype));
663 t_block gmx_mtop_global_cgs(const gmx_mtop_t *mtop)
665 t_block cgs_gl,*cgs_mol;
667 gmx_molblock_t *molb;
670 /* In most cases this is too much, but we realloc at the end */
671 snew(cgs_gl.index,mtop->natoms+1);
675 for(mb=0; mb<mtop->nmolblock; mb++)
677 molb = &mtop->molblock[mb];
678 cgs_mol = &mtop->moltype[molb->type].cgs;
679 for(mol=0; mol<molb->nmol; mol++)
681 for(cg=0; cg<cgs_mol->nr; cg++)
683 cgs_gl.index[cgs_gl.nr+1] =
684 cgs_gl.index[cgs_gl.nr] +
685 cgs_mol->index[cg+1] - cgs_mol->index[cg];
690 cgs_gl.nalloc_index = cgs_gl.nr + 1;
691 srenew(cgs_gl.index,cgs_gl.nalloc_index);
696 static void atomcat(t_atoms *dest, t_atoms *src, int copies,
697 int maxres_renum, int *maxresnr)
705 size=destnr+copies*srcnr;
706 srenew(dest->atom,size);
707 srenew(dest->atomname,size);
708 srenew(dest->atomtype,size);
709 srenew(dest->atomtypeB,size);
713 size=dest->nres+copies*src->nres;
714 srenew(dest->resinfo,size);
717 /* residue information */
718 for (l=dest->nres,j=0; (j<copies); j++,l+=src->nres)
720 memcpy((char *) &(dest->resinfo[l]),(char *) &(src->resinfo[0]),
721 (size_t)(src->nres*sizeof(src->resinfo[0])));
724 for (l=destnr,j=0; (j<copies); j++,l+=srcnr)
726 memcpy((char *) &(dest->atomname[l]),(char *) &(src->atomname[0]),
727 (size_t)(srcnr*sizeof(src->atomname[0])));
728 memcpy((char *) &(dest->atomtype[l]),(char *) &(src->atomtype[0]),
729 (size_t)(srcnr*sizeof(src->atomtype[0])));
730 memcpy((char *) &(dest->atomtypeB[l]),(char *) &(src->atomtypeB[0]),
731 (size_t)(srcnr*sizeof(src->atomtypeB[0])));
732 memcpy((char *) &(dest->atom[l]),(char *) &(src->atom[0]),
733 (size_t)(srcnr*sizeof(src->atom[0])));
736 /* Increment residue indices */
737 for (l=destnr,j=0; (j<copies); j++)
739 for (i=0; (i<srcnr); i++,l++)
741 dest->atom[l].resind = dest->nres+j*src->nres+src->atom[i].resind;
745 if (src->nres <= maxres_renum)
747 /* Single residue molecule, continue counting residues */
748 for (j=0; (j<copies); j++)
750 for (l=0; l<src->nres; l++)
753 dest->resinfo[dest->nres+j*src->nres+l].nr = *maxresnr;
758 dest->nres += copies*src->nres;
759 dest->nr += copies*src->nr;
762 t_atoms gmx_mtop_global_atoms(const gmx_mtop_t *mtop)
766 gmx_molblock_t *molb;
768 init_t_atoms(&atoms,0,FALSE);
770 maxresnr = mtop->maxresnr;
771 for(mb=0; mb<mtop->nmolblock; mb++)
773 molb = &mtop->molblock[mb];
774 atomcat(&atoms,&mtop->moltype[molb->type].atoms,molb->nmol,
775 mtop->maxres_renum,&maxresnr);
781 void gmx_mtop_make_atomic_charge_groups(gmx_mtop_t *mtop,
782 gmx_bool bKeepSingleMolCG)
787 for(mb=0; mb<mtop->nmolblock; mb++)
789 cgs_mol = &mtop->moltype[mtop->molblock[mb].type].cgs;
790 if (!(bKeepSingleMolCG && cgs_mol->nr == 1))
792 cgs_mol->nr = mtop->molblock[mb].natoms_mol;
793 cgs_mol->nalloc_index = cgs_mol->nr + 1;
794 srenew(cgs_mol->index,cgs_mol->nalloc_index);
795 for(cg=0; cg<cgs_mol->nr+1; cg++)
797 cgs_mol->index[cg] = cg;
804 * The cat routines below are old code from src/kernel/topcat.c
807 static void blockcat(t_block *dest,t_block *src,int copies,
814 size=(dest->nr+copies*src->nr+1);
815 srenew(dest->index,size);
818 nra = dest->index[dest->nr];
819 for (l=dest->nr,j=0; (j<copies); j++)
821 for (i=0; (i<src->nr); i++)
823 dest->index[l++] = nra + src->index[i];
825 nra += src->index[src->nr];
827 dest->nr += copies*src->nr;
828 dest->index[dest->nr] = nra;
831 static void blockacat(t_blocka *dest,t_blocka *src,int copies,
835 int destnr = dest->nr;
836 int destnra = dest->nra;
840 size=(dest->nr+copies*src->nr+1);
841 srenew(dest->index,size);
845 size=(dest->nra+copies*src->nra);
846 srenew(dest->a,size);
849 for (l=destnr,j=0; (j<copies); j++)
851 for (i=0; (i<src->nr); i++)
853 dest->index[l++] = dest->nra+src->index[i];
855 dest->nra += src->nra;
857 for (l=destnra,j=0; (j<copies); j++)
859 for (i=0; (i<src->nra); i++)
861 dest->a[l++] = dnum+src->a[i];
866 dest->index[dest->nr] = dest->nra;
869 static void ilistcat(int ftype,t_ilist *dest,t_ilist *src,int copies,
876 dest->nalloc = dest->nr + copies*src->nr;
877 srenew(dest->iatoms,dest->nalloc);
879 for(c=0; c<copies; c++)
881 for(i=0; i<src->nr; )
883 dest->iatoms[dest->nr++] = src->iatoms[i++];
884 for(a=0; a<nral; a++)
886 dest->iatoms[dest->nr++] = dnum + src->iatoms[i++];
893 static void set_posres_params(t_idef *idef,gmx_molblock_t *molb,
900 il = &idef->il[F_POSRES];
902 idef->iparams_posres_nalloc = i1;
903 srenew(idef->iparams_posres,idef->iparams_posres_nalloc);
906 ip = &idef->iparams_posres[i];
907 /* Copy the force constants */
908 *ip = idef->iparams[il->iatoms[i*2]];
909 a_molb = il->iatoms[i*2+1] - a_offset;
910 if (molb->nposres_xA == 0)
912 gmx_incons("Position restraint coordinates are missing");
914 ip->posres.pos0A[XX] = molb->posres_xA[a_molb][XX];
915 ip->posres.pos0A[YY] = molb->posres_xA[a_molb][YY];
916 ip->posres.pos0A[ZZ] = molb->posres_xA[a_molb][ZZ];
917 if (molb->nposres_xB > 0)
919 ip->posres.pos0B[XX] = molb->posres_xB[a_molb][XX];
920 ip->posres.pos0B[YY] = molb->posres_xB[a_molb][YY];
921 ip->posres.pos0B[ZZ] = molb->posres_xB[a_molb][ZZ];
925 ip->posres.pos0B[XX] = ip->posres.pos0A[XX];
926 ip->posres.pos0B[YY] = ip->posres.pos0A[YY];
927 ip->posres.pos0B[ZZ] = ip->posres.pos0A[ZZ];
929 /* Set the parameter index for idef->iparams_posre */
934 static void set_fbposres_params(t_idef *idef,gmx_molblock_t *molb,
941 il = &idef->il[F_FBPOSRES];
943 idef->iparams_fbposres_nalloc = i1;
944 srenew(idef->iparams_fbposres,idef->iparams_fbposres_nalloc);
947 ip = &idef->iparams_fbposres[i];
948 /* Copy the force constants */
949 *ip = idef->iparams[il->iatoms[i*2]];
950 a_molb = il->iatoms[i*2+1] - a_offset;
951 if (molb->nposres_xA == 0)
953 gmx_incons("Position restraint coordinates are missing");
955 /* Take flat-bottom posres reference from normal position restraints */
956 ip->fbposres.pos0[XX] = molb->posres_xA[a_molb][XX];
957 ip->fbposres.pos0[YY] = molb->posres_xA[a_molb][YY];
958 ip->fbposres.pos0[ZZ] = molb->posres_xA[a_molb][ZZ];
959 /* Note: no B-type for flat-bottom posres */
961 /* Set the parameter index for idef->iparams_posre */
966 static void gen_local_top(const gmx_mtop_t *mtop,const t_inputrec *ir,
967 gmx_bool bMergeConstr,
970 int mb,srcnr,destnr,ftype,ftype_dest,mt,natoms,mol,nposre_old,nfbposre_old;
971 gmx_molblock_t *molb;
973 const gmx_ffparams_t *ffp;
976 gmx_mtop_atomloop_all_t aloop;
980 top->atomtypes = mtop->atomtypes;
982 ffp = &mtop->ffparams;
985 idef->ntypes = ffp->ntypes;
986 idef->atnr = ffp->atnr;
987 idef->functype = ffp->functype;
988 idef->iparams = ffp->iparams;
989 idef->iparams_posres = NULL;
990 idef->iparams_posres_nalloc = 0;
991 idef->iparams_fbposres = NULL;
992 idef->iparams_fbposres_nalloc = 0;
993 idef->fudgeQQ = ffp->fudgeQQ;
994 idef->cmap_grid = ffp->cmap_grid;
995 idef->ilsort = ilsortUNKNOWN;
997 init_block(&top->cgs);
998 init_blocka(&top->excls);
999 for(ftype=0; ftype<F_NRE; ftype++)
1001 idef->il[ftype].nr = 0;
1002 idef->il[ftype].nalloc = 0;
1003 idef->il[ftype].iatoms = NULL;
1007 for(mb=0; mb<mtop->nmolblock; mb++)
1009 molb = &mtop->molblock[mb];
1010 molt = &mtop->moltype[molb->type];
1012 srcnr = molt->atoms.nr;
1015 blockcat(&top->cgs,&molt->cgs,molb->nmol,destnr,srcnr);
1017 blockacat(&top->excls,&molt->excls,molb->nmol,destnr,srcnr);
1019 nposre_old = idef->il[F_POSRES].nr;
1020 nfbposre_old = idef->il[F_FBPOSRES].nr;
1021 for(ftype=0; ftype<F_NRE; ftype++)
1024 ftype == F_CONSTR && molt->ilist[F_CONSTRNC].nr > 0)
1026 /* Merge all constrains into one ilist.
1027 * This simplifies the constraint code.
1029 for(mol=0; mol<molb->nmol; mol++) {
1030 ilistcat(ftype,&idef->il[F_CONSTR],&molt->ilist[F_CONSTR],
1031 1,destnr+mol*srcnr,srcnr);
1032 ilistcat(ftype,&idef->il[F_CONSTR],&molt->ilist[F_CONSTRNC],
1033 1,destnr+mol*srcnr,srcnr);
1036 else if (!(bMergeConstr && ftype == F_CONSTRNC))
1038 ilistcat(ftype,&idef->il[ftype],&molt->ilist[ftype],
1039 molb->nmol,destnr,srcnr);
1042 if (idef->il[F_POSRES].nr > nposre_old)
1044 /* Executing this line line stops gmxdump -sys working
1045 * correctly. I'm not aware there's an elegant fix. */
1046 set_posres_params(idef,molb,nposre_old/2,natoms);
1048 if (idef->il[F_FBPOSRES].nr > nfbposre_old)
1050 set_fbposres_params(idef,molb,nfbposre_old/2,natoms);
1053 natoms += molb->nmol*srcnr;
1058 top->idef.ilsort = ilsortUNKNOWN;
1062 if (ir->efep != efepNO && gmx_mtop_bondeds_free_energy(mtop))
1064 snew(qA,mtop->natoms);
1065 snew(qB,mtop->natoms);
1066 aloop = gmx_mtop_atomloop_all_init(mtop);
1067 while (gmx_mtop_atomloop_all_next(aloop,&ag,&atom))
1072 gmx_sort_ilist_fe(&top->idef,qA,qB);
1078 top->idef.ilsort = ilsortNO_FE;
1083 gmx_localtop_t *gmx_mtop_generate_local_top(const gmx_mtop_t *mtop,
1084 const t_inputrec *ir)
1086 gmx_localtop_t *top;
1090 gen_local_top(mtop,ir,TRUE,top);
1095 t_topology gmx_mtop_t_to_t_topology(gmx_mtop_t *mtop)
1098 gmx_localtop_t ltop;
1101 gen_local_top(mtop,NULL,FALSE,<op);
1103 top.name = mtop->name;
1104 top.idef = ltop.idef;
1105 top.atomtypes = ltop.atomtypes;
1107 top.excls = ltop.excls;
1108 top.atoms = gmx_mtop_global_atoms(mtop);
1109 top.mols = mtop->mols;
1110 top.symtab = mtop->symtab;
1112 /* We only need to free the moltype and molblock data,
1113 * all other pointers have been copied to top.
1115 * Well, except for the group data, but we can't free those, because they
1116 * are used somewhere even after a call to this function.
1118 for(mt=0; mt<mtop->nmoltype; mt++)
1120 done_moltype(&mtop->moltype[mt]);
1122 sfree(mtop->moltype);
1124 for(mb=0; mb<mtop->nmolblock; mb++)
1126 done_molblock(&mtop->molblock[mb]);
1128 sfree(mtop->molblock);