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 gen_local_top(const gmx_mtop_t
*mtop
,const t_inputrec
*ir
,
935 gmx_bool bMergeConstr
,
938 int mb
,srcnr
,destnr
,ftype
,ftype_dest
,mt
,natoms
,mol
,nposre_old
;
939 gmx_molblock_t
*molb
;
941 const gmx_ffparams_t
*ffp
;
944 gmx_mtop_atomloop_all_t aloop
;
948 top
->atomtypes
= mtop
->atomtypes
;
950 ffp
= &mtop
->ffparams
;
953 idef
->ntypes
= ffp
->ntypes
;
954 idef
->atnr
= ffp
->atnr
;
955 idef
->functype
= ffp
->functype
;
956 idef
->iparams
= ffp
->iparams
;
957 idef
->iparams_posres
= NULL
;
958 idef
->iparams_posres_nalloc
= 0;
959 idef
->fudgeQQ
= ffp
->fudgeQQ
;
960 idef
->cmap_grid
= ffp
->cmap_grid
;
961 idef
->ilsort
= ilsortUNKNOWN
;
963 init_block(&top
->cgs
);
964 init_blocka(&top
->excls
);
965 for(ftype
=0; ftype
<F_NRE
; ftype
++)
967 idef
->il
[ftype
].nr
= 0;
968 idef
->il
[ftype
].nalloc
= 0;
969 idef
->il
[ftype
].iatoms
= NULL
;
973 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
975 molb
= &mtop
->molblock
[mb
];
976 molt
= &mtop
->moltype
[molb
->type
];
978 srcnr
= molt
->atoms
.nr
;
981 blockcat(&top
->cgs
,&molt
->cgs
,molb
->nmol
,destnr
,srcnr
);
983 blockacat(&top
->excls
,&molt
->excls
,molb
->nmol
,destnr
,srcnr
);
985 nposre_old
= idef
->il
[F_POSRES
].nr
;
986 for(ftype
=0; ftype
<F_NRE
; ftype
++)
989 ftype
== F_CONSTR
&& molt
->ilist
[F_CONSTRNC
].nr
> 0)
991 /* Merge all constrains into one ilist.
992 * This simplifies the constraint code.
994 for(mol
=0; mol
<molb
->nmol
; mol
++) {
995 ilistcat(ftype
,&idef
->il
[F_CONSTR
],&molt
->ilist
[F_CONSTR
],
996 1,destnr
+mol
*srcnr
,srcnr
);
997 ilistcat(ftype
,&idef
->il
[F_CONSTR
],&molt
->ilist
[F_CONSTRNC
],
998 1,destnr
+mol
*srcnr
,srcnr
);
1001 else if (!(bMergeConstr
&& ftype
== F_CONSTRNC
))
1003 ilistcat(ftype
,&idef
->il
[ftype
],&molt
->ilist
[ftype
],
1004 molb
->nmol
,destnr
,srcnr
);
1007 if (idef
->il
[F_POSRES
].nr
> nposre_old
)
1009 set_posres_params(idef
,molb
,nposre_old
/2,natoms
);
1012 natoms
+= molb
->nmol
*srcnr
;
1017 top
->idef
.ilsort
= ilsortUNKNOWN
;
1021 if (ir
->efep
!= efepNO
&& gmx_mtop_bondeds_free_energy(mtop
))
1023 snew(qA
,mtop
->natoms
);
1024 snew(qB
,mtop
->natoms
);
1025 aloop
= gmx_mtop_atomloop_all_init(mtop
);
1026 while (gmx_mtop_atomloop_all_next(aloop
,&ag
,&atom
))
1031 gmx_sort_ilist_fe(&top
->idef
,qA
,qB
);
1037 top
->idef
.ilsort
= ilsortNO_FE
;
1042 gmx_localtop_t
*gmx_mtop_generate_local_top(const gmx_mtop_t
*mtop
,
1043 const t_inputrec
*ir
)
1045 gmx_localtop_t
*top
;
1049 gen_local_top(mtop
,ir
,TRUE
,top
);
1054 t_topology
gmx_mtop_t_to_t_topology(gmx_mtop_t
*mtop
)
1057 gmx_localtop_t ltop
;
1060 gen_local_top(mtop
,NULL
,FALSE
,<op
);
1062 top
.name
= mtop
->name
;
1063 top
.idef
= ltop
.idef
;
1064 top
.atomtypes
= ltop
.atomtypes
;
1066 top
.excls
= ltop
.excls
;
1067 top
.atoms
= gmx_mtop_global_atoms(mtop
);
1068 top
.mols
= mtop
->mols
;
1069 top
.symtab
= mtop
->symtab
;
1071 /* We only need to free the moltype and molblock data,
1072 * all other pointers have been copied to top.
1074 * Well, except for the group data, but we can't free those, because they
1075 * are used somewhere even after a call to this function.
1077 for(mt
=0; mt
<mtop
->nmoltype
; mt
++)
1079 done_moltype(&mtop
->moltype
[mt
]);
1081 sfree(mtop
->moltype
);
1083 for(mb
=0; mb
<mtop
->nmolblock
; mb
++)
1085 done_molblock(&mtop
->molblock
[mb
]);
1087 sfree(mtop
->molblock
);