2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/math/vecdump.h"
46 #include "gromacs/topology/atoms.h"
47 #include "gromacs/topology/idef.h"
48 #include "gromacs/topology/ifunc.h"
49 #include "gromacs/topology/symtab.h"
50 #include "gromacs/utility/compare.h"
51 #include "gromacs/utility/smalloc.h"
52 #include "gromacs/utility/stringutil.h"
53 #include "gromacs/utility/txtdump.h"
55 const char *gtypes
[egcNR
+1] = {
56 "T-Coupling", "Energy Mon.", "Acceleration", "Freeze",
57 "User1", "User2", "VCM", "Compressed X", "Or. Res. Fit", "QMMM", NULL
60 static void init_groups(gmx_groups_t
*groups
)
63 groups
->grpname
= NULL
;
64 for (int g
= 0; g
< egcNR
; g
++)
66 groups
->grps
[g
].nm_ind
= NULL
;
67 groups
->ngrpnr
[g
] = 0;
68 groups
->grpnr
[g
] = NULL
;
73 void init_mtop(gmx_mtop_t
*mtop
)
79 mtop
->molblock
= NULL
;
80 mtop
->maxres_renum
= 0;
82 init_groups(&mtop
->groups
);
83 init_block(&mtop
->mols
);
84 open_symtab(&mtop
->symtab
);
87 void init_top(t_topology
*top
)
90 init_atom(&(top
->atoms
));
91 init_atomtypes(&(top
->atomtypes
));
92 init_block(&top
->cgs
);
93 init_block(&top
->mols
);
94 init_blocka(&top
->excls
);
95 open_symtab(&top
->symtab
);
99 void done_moltype(gmx_moltype_t
*molt
)
101 done_atom(&molt
->atoms
);
102 done_block(&molt
->cgs
);
103 done_blocka(&molt
->excls
);
105 for (int f
= 0; f
< F_NRE
; f
++)
107 sfree(molt
->ilist
[f
].iatoms
);
108 molt
->ilist
[f
].nalloc
= 0;
112 void done_molblock(gmx_molblock_t
*molb
)
114 if (molb
->nposres_xA
> 0)
116 molb
->nposres_xA
= 0;
117 sfree(molb
->posres_xA
);
119 if (molb
->nposres_xB
> 0)
121 molb
->nposres_xB
= 0;
122 sfree(molb
->posres_xB
);
126 void done_gmx_groups_t(gmx_groups_t
*g
)
130 for (i
= 0; (i
< egcNR
); i
++)
132 if (NULL
!= g
->grps
[i
].nm_ind
)
134 sfree(g
->grps
[i
].nm_ind
);
135 g
->grps
[i
].nm_ind
= NULL
;
137 if (NULL
!= g
->grpnr
[i
])
143 /* The contents of this array is in symtab, don't free it here */
147 void done_mtop(gmx_mtop_t
*mtop
)
149 done_symtab(&mtop
->symtab
);
151 sfree(mtop
->ffparams
.functype
);
152 sfree(mtop
->ffparams
.iparams
);
154 for (int i
= 0; i
< mtop
->nmoltype
; i
++)
156 done_moltype(&mtop
->moltype
[i
]);
158 sfree(mtop
->moltype
);
159 for (int i
= 0; i
< mtop
->nmolblock
; i
++)
161 done_molblock(&mtop
->molblock
[i
]);
163 sfree(mtop
->molblock
);
164 done_atomtypes(&mtop
->atomtypes
);
165 done_gmx_groups_t(&mtop
->groups
);
166 done_block(&mtop
->mols
);
169 void done_top(t_topology
*top
)
171 sfree(top
->idef
.functype
);
172 sfree(top
->idef
.iparams
);
173 for (int f
= 0; f
< F_NRE
; ++f
)
175 sfree(top
->idef
.il
[f
].iatoms
);
176 top
->idef
.il
[f
].iatoms
= NULL
;
177 top
->idef
.il
[f
].nalloc
= 0;
180 done_atom(&(top
->atoms
));
183 done_atomtypes(&(top
->atomtypes
));
185 done_symtab(&(top
->symtab
));
186 done_block(&(top
->cgs
));
187 done_block(&(top
->mols
));
188 done_blocka(&(top
->excls
));
191 void done_top_mtop(t_topology
*top
, gmx_mtop_t
*mtop
)
197 for (int f
= 0; f
< F_NRE
; ++f
)
199 sfree(top
->idef
.il
[f
].iatoms
);
200 top
->idef
.il
[f
].iatoms
= NULL
;
201 top
->idef
.il
[f
].nalloc
= 0;
203 done_atom(&top
->atoms
);
204 done_block(&top
->cgs
);
205 done_blocka(&top
->excls
);
206 done_symtab(&top
->symtab
);
207 open_symtab(&mtop
->symtab
);
213 bool gmx_mtop_has_masses(const gmx_mtop_t
*mtop
)
219 return mtop
->nmoltype
== 0 || mtop
->moltype
[0].atoms
.haveMass
;
222 bool gmx_mtop_has_charges(const gmx_mtop_t
*mtop
)
228 return mtop
->nmoltype
== 0 || mtop
->moltype
[0].atoms
.haveCharge
;
231 bool gmx_mtop_has_atomtypes(const gmx_mtop_t
*mtop
)
237 return mtop
->nmoltype
== 0 || mtop
->moltype
[0].atoms
.haveType
;
240 bool gmx_mtop_has_pdbinfo(const gmx_mtop_t
*mtop
)
246 return mtop
->nmoltype
== 0 || mtop
->moltype
[0].atoms
.havePdbInfo
;
249 static void pr_grps(FILE *fp
, const char *title
, const t_grps grps
[], char **grpname
[])
253 for (i
= 0; (i
< egcNR
); i
++)
255 fprintf(fp
, "%s[%-12s] nr=%d, name=[", title
, gtypes
[i
], grps
[i
].nr
);
256 for (j
= 0; (j
< grps
[i
].nr
); j
++)
258 fprintf(fp
, " %s", *(grpname
[grps
[i
].nm_ind
[j
]]));
264 static void pr_groups(FILE *fp
, int indent
,
265 const gmx_groups_t
*groups
,
266 gmx_bool bShowNumbers
)
270 pr_grps(fp
, "grp", groups
->grps
, groups
->grpname
);
271 pr_strings(fp
, indent
, "grpname", groups
->grpname
, groups
->ngrpname
, bShowNumbers
);
273 pr_indent(fp
, indent
);
274 fprintf(fp
, "groups ");
275 for (g
= 0; g
< egcNR
; g
++)
277 printf(" %5.5s", gtypes
[g
]);
281 pr_indent(fp
, indent
);
282 fprintf(fp
, "allocated ");
284 for (g
= 0; g
< egcNR
; g
++)
286 printf(" %5d", groups
->ngrpnr
[g
]);
287 nat_max
= std::max(nat_max
, groups
->ngrpnr
[g
]);
293 pr_indent(fp
, indent
);
294 fprintf(fp
, "groupnr[%5s] =", "*");
295 for (g
= 0; g
< egcNR
; g
++)
297 fprintf(fp
, " %3d ", 0);
303 for (i
= 0; i
< nat_max
; i
++)
305 pr_indent(fp
, indent
);
306 fprintf(fp
, "groupnr[%5d] =", i
);
307 for (g
= 0; g
< egcNR
; g
++)
310 groups
->grpnr
[g
] ? groups
->grpnr
[g
][i
] : 0);
317 static void pr_moltype(FILE *fp
, int indent
, const char *title
,
318 const gmx_moltype_t
*molt
, int n
,
319 const gmx_ffparams_t
*ffparams
,
320 gmx_bool bShowNumbers
)
324 indent
= pr_title_n(fp
, indent
, title
, n
);
325 pr_indent(fp
, indent
);
326 fprintf(fp
, "name=\"%s\"\n", *(molt
->name
));
327 pr_atoms(fp
, indent
, "atoms", &(molt
->atoms
), bShowNumbers
);
328 pr_block(fp
, indent
, "cgs", &molt
->cgs
, bShowNumbers
);
329 pr_blocka(fp
, indent
, "excls", &molt
->excls
, bShowNumbers
);
330 for (j
= 0; (j
< F_NRE
); j
++)
332 pr_ilist(fp
, indent
, interaction_function
[j
].longname
,
333 ffparams
->functype
, &molt
->ilist
[j
], bShowNumbers
);
337 static void pr_molblock(FILE *fp
, int indent
, const char *title
,
338 const gmx_molblock_t
*molb
, int n
,
339 const gmx_moltype_t
*molt
)
341 indent
= pr_title_n(fp
, indent
, title
, n
);
342 pr_indent(fp
, indent
);
343 fprintf(fp
, "%-20s = %d \"%s\"\n",
344 "moltype", molb
->type
, *(molt
[molb
->type
].name
));
345 pr_int(fp
, indent
, "#molecules", molb
->nmol
);
346 pr_int(fp
, indent
, "#atoms_mol", molb
->natoms_mol
);
347 pr_int(fp
, indent
, "#posres_xA", molb
->nposres_xA
);
348 if (molb
->nposres_xA
> 0)
350 pr_rvecs(fp
, indent
, "posres_xA", molb
->posres_xA
, molb
->nposres_xA
);
352 pr_int(fp
, indent
, "#posres_xB", molb
->nposres_xB
);
353 if (molb
->nposres_xB
> 0)
355 pr_rvecs(fp
, indent
, "posres_xB", molb
->posres_xB
, molb
->nposres_xB
);
359 void pr_mtop(FILE *fp
, int indent
, const char *title
, const gmx_mtop_t
*mtop
,
360 gmx_bool bShowNumbers
)
364 if (available(fp
, mtop
, indent
, title
))
366 indent
= pr_title(fp
, indent
, title
);
367 pr_indent(fp
, indent
);
368 fprintf(fp
, "name=\"%s\"\n", *(mtop
->name
));
369 pr_int(fp
, indent
, "#atoms", mtop
->natoms
);
370 pr_int(fp
, indent
, "#molblock", mtop
->nmolblock
);
371 for (mb
= 0; mb
< mtop
->nmolblock
; mb
++)
373 pr_molblock(fp
, indent
, "molblock", &mtop
->molblock
[mb
], mb
, mtop
->moltype
);
375 pr_str(fp
, indent
, "bIntermolecularInteractions",
376 gmx::boolToString(mtop
->bIntermolecularInteractions
));
377 if (mtop
->bIntermolecularInteractions
)
379 for (j
= 0; (j
< F_NRE
); j
++)
381 pr_ilist(fp
, indent
, interaction_function
[j
].longname
,
382 mtop
->ffparams
.functype
,
383 &mtop
->intermolecular_ilist
[j
], bShowNumbers
);
386 pr_ffparams(fp
, indent
, "ffparams", &(mtop
->ffparams
), bShowNumbers
);
387 pr_atomtypes(fp
, indent
, "atomtypes", &(mtop
->atomtypes
), bShowNumbers
);
388 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
390 pr_moltype(fp
, indent
, "moltype", &mtop
->moltype
[mt
], mt
,
391 &mtop
->ffparams
, bShowNumbers
);
393 pr_groups(fp
, indent
, &mtop
->groups
, bShowNumbers
);
397 void pr_top(FILE *fp
, int indent
, const char *title
, const t_topology
*top
, gmx_bool bShowNumbers
)
399 if (available(fp
, top
, indent
, title
))
401 indent
= pr_title(fp
, indent
, title
);
402 pr_indent(fp
, indent
);
403 fprintf(fp
, "name=\"%s\"\n", *(top
->name
));
404 pr_atoms(fp
, indent
, "atoms", &(top
->atoms
), bShowNumbers
);
405 pr_atomtypes(fp
, indent
, "atomtypes", &(top
->atomtypes
), bShowNumbers
);
406 pr_block(fp
, indent
, "cgs", &top
->cgs
, bShowNumbers
);
407 pr_block(fp
, indent
, "mols", &top
->mols
, bShowNumbers
);
408 pr_str(fp
, indent
, "bIntermolecularInteractions",
409 gmx::boolToString(top
->bIntermolecularInteractions
));
410 pr_blocka(fp
, indent
, "excls", &top
->excls
, bShowNumbers
);
411 pr_idef(fp
, indent
, "idef", &top
->idef
, bShowNumbers
);
415 static void cmp_ilist(FILE *fp
, int ftype
, const t_ilist
*il1
, const t_ilist
*il2
)
420 fprintf(fp
, "comparing ilist %s\n", interaction_function
[ftype
].name
);
421 sprintf(buf
, "%s->nr", interaction_function
[ftype
].name
);
422 cmp_int(fp
, buf
, -1, il1
->nr
, il2
->nr
);
423 sprintf(buf
, "%s->iatoms", interaction_function
[ftype
].name
);
424 if (((il1
->nr
> 0) && (!il1
->iatoms
)) ||
425 ((il2
->nr
> 0) && (!il2
->iatoms
)) ||
426 ((il1
->nr
!= il2
->nr
)))
428 fprintf(fp
, "Comparing radically different topologies - %s is different\n",
433 for (i
= 0; (i
< il1
->nr
); i
++)
435 cmp_int(fp
, buf
, i
, il1
->iatoms
[i
], il2
->iatoms
[i
]);
440 static void cmp_iparm(FILE *fp
, const char *s
, t_functype ft
,
441 const t_iparams
&ip1
, const t_iparams
&ip2
, real ftol
, real abstol
)
447 for (i
= 0; i
< MAXFORCEPARAM
&& !bDiff
; i
++)
449 bDiff
= !equal_real(ip1
.generic
.buf
[i
], ip2
.generic
.buf
[i
], ftol
, abstol
);
453 fprintf(fp
, "%s1: ", s
);
454 pr_iparams(fp
, ft
, &ip1
);
455 fprintf(fp
, "%s2: ", s
);
456 pr_iparams(fp
, ft
, &ip2
);
460 static void cmp_iparm_AB(FILE *fp
, const char *s
, t_functype ft
,
461 const t_iparams
&ip1
, real ftol
, real abstol
)
463 int nrfpA
, nrfpB
, p0
, i
;
466 /* Normally the first parameter is perturbable */
468 nrfpA
= interaction_function
[ft
].nrfpA
;
469 nrfpB
= interaction_function
[ft
].nrfpB
;
474 else if (interaction_function
[ft
].flags
& IF_TABULATED
)
476 /* For tabulated interactions only the second parameter is perturbable */
481 for (i
= 0; i
< nrfpB
&& !bDiff
; i
++)
483 bDiff
= !equal_real(ip1
.generic
.buf
[p0
+i
], ip1
.generic
.buf
[nrfpA
+i
], ftol
, abstol
);
487 fprintf(fp
, "%s: ", s
);
488 pr_iparams(fp
, ft
, &ip1
);
492 static void cmp_cmap(FILE *fp
, const gmx_cmap_t
*cmap1
, const gmx_cmap_t
*cmap2
, real ftol
, real abstol
)
494 cmp_int(fp
, "cmap ngrid", -1, cmap1
->ngrid
, cmap2
->ngrid
);
495 cmp_int(fp
, "cmap grid_spacing", -1, cmap1
->grid_spacing
, cmap2
->grid_spacing
);
496 if (cmap1
->ngrid
== cmap2
->ngrid
&&
497 cmap1
->grid_spacing
== cmap2
->grid_spacing
)
501 for (g
= 0; g
< cmap1
->ngrid
; g
++)
505 fprintf(fp
, "comparing cmap %d\n", g
);
507 for (i
= 0; i
< 4*cmap1
->grid_spacing
*cmap1
->grid_spacing
; i
++)
509 cmp_real(fp
, "", i
, cmap1
->cmapdata
[g
].cmap
[i
], cmap2
->cmapdata
[g
].cmap
[i
], ftol
, abstol
);
515 static void cmp_idef(FILE *fp
, const t_idef
*id1
, const t_idef
*id2
, real ftol
, real abstol
)
518 char buf1
[64], buf2
[64];
520 fprintf(fp
, "comparing idef\n");
523 cmp_int(fp
, "idef->ntypes", -1, id1
->ntypes
, id2
->ntypes
);
524 cmp_int(fp
, "idef->atnr", -1, id1
->atnr
, id2
->atnr
);
525 for (i
= 0; (i
< std::min(id1
->ntypes
, id2
->ntypes
)); i
++)
527 sprintf(buf1
, "idef->functype[%d]", i
);
528 sprintf(buf2
, "idef->iparam[%d]", i
);
529 cmp_int(fp
, buf1
, i
, (int)id1
->functype
[i
], (int)id2
->functype
[i
]);
530 cmp_iparm(fp
, buf2
, id1
->functype
[i
],
531 id1
->iparams
[i
], id2
->iparams
[i
], ftol
, abstol
);
533 cmp_real(fp
, "fudgeQQ", -1, id1
->fudgeQQ
, id2
->fudgeQQ
, ftol
, abstol
);
534 cmp_cmap(fp
, &id1
->cmap_grid
, &id2
->cmap_grid
, ftol
, abstol
);
535 for (i
= 0; (i
< F_NRE
); i
++)
537 cmp_ilist(fp
, i
, &(id1
->il
[i
]), &(id2
->il
[i
]));
542 for (i
= 0; (i
< id1
->ntypes
); i
++)
544 cmp_iparm_AB(fp
, "idef->iparam", id1
->functype
[i
], id1
->iparams
[i
], ftol
, abstol
);
549 static void cmp_block(FILE *fp
, const t_block
*b1
, const t_block
*b2
, const char *s
)
553 fprintf(fp
, "comparing block %s\n", s
);
554 sprintf(buf
, "%s.nr", s
);
555 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
558 static void cmp_blocka(FILE *fp
, const t_blocka
*b1
, const t_blocka
*b2
, const char *s
)
562 fprintf(fp
, "comparing blocka %s\n", s
);
563 sprintf(buf
, "%s.nr", s
);
564 cmp_int(fp
, buf
, -1, b1
->nr
, b2
->nr
);
565 sprintf(buf
, "%s.nra", s
);
566 cmp_int(fp
, buf
, -1, b1
->nra
, b2
->nra
);
569 void cmp_top(FILE *fp
, const t_topology
*t1
, const t_topology
*t2
, real ftol
, real abstol
)
571 fprintf(fp
, "comparing top\n");
574 cmp_idef(fp
, &(t1
->idef
), &(t2
->idef
), ftol
, abstol
);
575 cmp_atoms(fp
, &(t1
->atoms
), &(t2
->atoms
), ftol
, abstol
);
576 cmp_block(fp
, &t1
->cgs
, &t2
->cgs
, "cgs");
577 cmp_block(fp
, &t1
->mols
, &t2
->mols
, "mols");
578 cmp_bool(fp
, "bIntermolecularInteractions", -1, t1
->bIntermolecularInteractions
, t2
->bIntermolecularInteractions
);
579 cmp_blocka(fp
, &t1
->excls
, &t2
->excls
, "excls");
583 cmp_idef(fp
, &(t1
->idef
), NULL
, ftol
, abstol
);
584 cmp_atoms(fp
, &(t1
->atoms
), NULL
, ftol
, abstol
);
588 void cmp_groups(FILE *fp
, const gmx_groups_t
*g0
, const gmx_groups_t
*g1
,
589 int natoms0
, int natoms1
)
594 fprintf(fp
, "comparing groups\n");
596 for (i
= 0; i
< egcNR
; i
++)
598 sprintf(buf
, "grps[%d].nr", i
);
599 cmp_int(fp
, buf
, -1, g0
->grps
[i
].nr
, g1
->grps
[i
].nr
);
600 if (g0
->grps
[i
].nr
== g1
->grps
[i
].nr
)
602 for (j
= 0; j
< g0
->grps
[i
].nr
; j
++)
604 sprintf(buf
, "grps[%d].name[%d]", i
, j
);
606 *g0
->grpname
[g0
->grps
[i
].nm_ind
[j
]],
607 *g1
->grpname
[g1
->grps
[i
].nm_ind
[j
]]);
610 cmp_int(fp
, "ngrpnr", i
, g0
->ngrpnr
[i
], g1
->ngrpnr
[i
]);
611 if (g0
->ngrpnr
[i
] == g1
->ngrpnr
[i
] && natoms0
== natoms1
&&
612 (g0
->grpnr
[i
] != NULL
|| g1
->grpnr
[i
] != NULL
))
614 for (j
= 0; j
< natoms0
; j
++)
616 cmp_int(fp
, gtypes
[i
], j
, ggrpnr(g0
, i
, j
), ggrpnr(g1
, i
, j
));
620 /* We have compared the names in the groups lists,
621 * so we can skip the grpname list comparison.