Convert gmx_mtop_t to C++
[gromacs.git] / src / gromacs / topology / topology.cpp
blob198e1a5af2f18d38a71c31982d57f85fc1958fdc
1 /*
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,2017,2018, 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.
37 #include "gmxpre.h"
39 #include "topology.h"
41 #include <cstdio>
43 #include <algorithm>
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/strconvert.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", nullptr
60 static void init_groups(gmx_groups_t *groups)
62 groups->ngrpname = 0;
63 groups->grpname = nullptr;
64 for (int g = 0; g < egcNR; g++)
66 groups->grps[g].nr = 0;
67 groups->grps[g].nm_ind = nullptr;
68 groups->ngrpnr[g] = 0;
69 groups->grpnr[g] = nullptr;
74 void init_mtop(gmx_mtop_t *mtop)
76 mtop->name = nullptr;
78 // TODO: Move to ffparams when that is converted to C++
79 mtop->ffparams.functype = nullptr;
80 mtop->ffparams.iparams = nullptr;
81 mtop->ffparams.cmap_grid.ngrid = 0;
82 mtop->ffparams.cmap_grid.grid_spacing = 0;
83 mtop->ffparams.cmap_grid.cmapdata = nullptr;
85 mtop->moltype.clear();
86 mtop->molblock.clear();
87 mtop->bIntermolecularInteractions = FALSE;
88 mtop->intermolecular_ilist = nullptr;
90 mtop->natoms = 0;
91 mtop->maxres_renum = 0;
92 mtop->maxresnr = -1;
93 init_atomtypes(&mtop->atomtypes);
94 init_groups(&mtop->groups);
95 open_symtab(&mtop->symtab);
98 void init_top(t_topology *top)
100 top->name = nullptr;
101 init_atom(&(top->atoms));
102 init_atomtypes(&(top->atomtypes));
103 init_block(&top->cgs);
104 init_block(&top->mols);
105 init_blocka(&top->excls);
106 open_symtab(&top->symtab);
110 gmx_moltype_t::gmx_moltype_t() :
111 name(nullptr),
112 cgs(),
113 excls()
115 init_t_atoms(&atoms, 0, FALSE);
117 for (int ftype = 0; ftype < F_NRE; ftype++)
119 ilist[ftype].nr = 0;
120 ilist[ftype].nr_nonperturbed = 0;
121 ilist[ftype].iatoms = nullptr;
122 ilist[ftype].nalloc = 0;
126 gmx_moltype_t::~gmx_moltype_t()
128 done_atom(&atoms);
129 done_block(&cgs);
130 done_blocka(&excls);
132 for (int f = 0; f < F_NRE; f++)
134 sfree(ilist[f].iatoms);
135 ilist[f].nalloc = 0;
139 gmx_molblock_t::gmx_molblock_t()
141 type = -1;
142 nmol = 0;
143 nposres_xA = 0;
144 posres_xA = nullptr;
145 nposres_xB = 0;
146 posres_xB = nullptr;
148 natoms_mol = 0;
149 globalAtomStart = -1;
150 globalAtomEnd = -1;
151 globalResidueStart = -1;
152 residueNumberStart = -1;
153 moleculeIndexStart = -1;
156 gmx_molblock_t::~gmx_molblock_t()
158 if (nposres_xA > 0)
160 nposres_xA = 0;
161 sfree(posres_xA);
163 if (nposres_xB > 0)
165 nposres_xB = 0;
166 sfree(posres_xB);
170 void done_gmx_groups_t(gmx_groups_t *g)
172 int i;
174 for (i = 0; (i < egcNR); i++)
176 if (nullptr != g->grps[i].nm_ind)
178 sfree(g->grps[i].nm_ind);
179 g->grps[i].nm_ind = nullptr;
181 if (nullptr != g->grpnr[i])
183 sfree(g->grpnr[i]);
184 g->grpnr[i] = nullptr;
187 /* The contents of this array is in symtab, don't free it here */
188 sfree(g->grpname);
191 gmx_mtop_t::gmx_mtop_t()
193 init_mtop(this);
196 gmx_mtop_t::~gmx_mtop_t()
198 done_symtab(&symtab);
200 sfree(ffparams.functype);
201 sfree(ffparams.iparams);
202 for (int i = 0; i < ffparams.cmap_grid.ngrid; i++)
204 sfree(ffparams.cmap_grid.cmapdata[i].cmap);
206 sfree(ffparams.cmap_grid.cmapdata);
208 moltype.clear();
209 molblock.clear();
210 done_atomtypes(&atomtypes);
211 done_gmx_groups_t(&groups);
214 void done_top(t_topology *top)
216 sfree(top->idef.functype);
217 sfree(top->idef.iparams);
218 for (int f = 0; f < F_NRE; ++f)
220 sfree(top->idef.il[f].iatoms);
221 top->idef.il[f].iatoms = nullptr;
222 top->idef.il[f].nalloc = 0;
225 done_atom(&(top->atoms));
227 /* For GB */
228 done_atomtypes(&(top->atomtypes));
230 done_symtab(&(top->symtab));
231 done_block(&(top->cgs));
232 done_block(&(top->mols));
233 done_blocka(&(top->excls));
236 void done_top_mtop(t_topology *top, gmx_mtop_t *mtop)
238 if (mtop != nullptr)
240 if (top != nullptr)
242 for (int f = 0; f < F_NRE; ++f)
244 sfree(top->idef.il[f].iatoms);
245 top->idef.il[f].iatoms = nullptr;
246 top->idef.il[f].nalloc = 0;
248 done_atom(&top->atoms);
249 done_block(&top->cgs);
250 done_blocka(&top->excls);
251 done_block(&top->mols);
252 done_symtab(&top->symtab);
253 open_symtab(&mtop->symtab);
256 // Note that the rest of mtop will be freed by the destructor
260 bool gmx_mtop_has_masses(const gmx_mtop_t *mtop)
262 if (mtop == nullptr)
264 return false;
266 return mtop->moltype.empty() || mtop->moltype[0].atoms.haveMass;
269 bool gmx_mtop_has_charges(const gmx_mtop_t *mtop)
271 if (mtop == nullptr)
273 return false;
275 return mtop->moltype.empty() || mtop->moltype[0].atoms.haveCharge;
278 bool gmx_mtop_has_atomtypes(const gmx_mtop_t *mtop)
280 if (mtop == nullptr)
282 return false;
284 return mtop->moltype.empty() || mtop->moltype[0].atoms.haveType;
287 bool gmx_mtop_has_pdbinfo(const gmx_mtop_t *mtop)
289 if (mtop == nullptr)
291 return false;
293 return mtop->moltype.empty() || mtop->moltype[0].atoms.havePdbInfo;
296 static void pr_grps(FILE *fp, const char *title, const t_grps grps[], char **grpname[])
298 int i, j;
300 for (i = 0; (i < egcNR); i++)
302 fprintf(fp, "%s[%-12s] nr=%d, name=[", title, gtypes[i], grps[i].nr);
303 for (j = 0; (j < grps[i].nr); j++)
305 fprintf(fp, " %s", *(grpname[grps[i].nm_ind[j]]));
307 fprintf(fp, "]\n");
311 static void pr_groups(FILE *fp, int indent,
312 const gmx_groups_t *groups,
313 gmx_bool bShowNumbers)
315 int nat_max, i, g;
317 pr_grps(fp, "grp", groups->grps, groups->grpname);
318 pr_strings(fp, indent, "grpname", groups->grpname, groups->ngrpname, bShowNumbers);
320 pr_indent(fp, indent);
321 fprintf(fp, "groups ");
322 for (g = 0; g < egcNR; g++)
324 printf(" %5.5s", gtypes[g]);
326 printf("\n");
328 pr_indent(fp, indent);
329 fprintf(fp, "allocated ");
330 nat_max = 0;
331 for (g = 0; g < egcNR; g++)
333 printf(" %5d", groups->ngrpnr[g]);
334 nat_max = std::max(nat_max, groups->ngrpnr[g]);
336 printf("\n");
338 if (nat_max == 0)
340 pr_indent(fp, indent);
341 fprintf(fp, "groupnr[%5s] =", "*");
342 for (g = 0; g < egcNR; g++)
344 fprintf(fp, " %3d ", 0);
346 fprintf(fp, "\n");
348 else
350 for (i = 0; i < nat_max; i++)
352 pr_indent(fp, indent);
353 fprintf(fp, "groupnr[%5d] =", i);
354 for (g = 0; g < egcNR; g++)
356 fprintf(fp, " %3d ",
357 groups->grpnr[g] ? groups->grpnr[g][i] : 0);
359 fprintf(fp, "\n");
364 static void pr_moltype(FILE *fp, int indent, const char *title,
365 const gmx_moltype_t *molt, int n,
366 const gmx_ffparams_t *ffparams,
367 gmx_bool bShowNumbers, gmx_bool bShowParameters)
369 int j;
371 indent = pr_title_n(fp, indent, title, n);
372 pr_indent(fp, indent);
373 fprintf(fp, "name=\"%s\"\n", *(molt->name));
374 pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers);
375 pr_block(fp, indent, "cgs", &molt->cgs, bShowNumbers);
376 pr_blocka(fp, indent, "excls", &molt->excls, bShowNumbers);
377 for (j = 0; (j < F_NRE); j++)
379 pr_ilist(fp, indent, interaction_function[j].longname,
380 ffparams->functype, &molt->ilist[j],
381 bShowNumbers, bShowParameters, ffparams->iparams);
385 static void pr_molblock(FILE *fp, int indent, const char *title,
386 const gmx_molblock_t *molb, int n,
387 const std::vector<gmx_moltype_t> &molt)
389 indent = pr_title_n(fp, indent, title, n);
390 pr_indent(fp, indent);
391 fprintf(fp, "%-20s = %d \"%s\"\n",
392 "moltype", molb->type, *(molt[molb->type].name));
393 pr_int(fp, indent, "#molecules", molb->nmol);
394 pr_int(fp, indent, "#atoms_mol", molb->natoms_mol);
395 pr_int(fp, indent, "#posres_xA", molb->nposres_xA);
396 if (molb->nposres_xA > 0)
398 pr_rvecs(fp, indent, "posres_xA", molb->posres_xA, molb->nposres_xA);
400 pr_int(fp, indent, "#posres_xB", molb->nposres_xB);
401 if (molb->nposres_xB > 0)
403 pr_rvecs(fp, indent, "posres_xB", molb->posres_xB, molb->nposres_xB);
407 void pr_mtop(FILE *fp, int indent, const char *title, const gmx_mtop_t *mtop,
408 gmx_bool bShowNumbers, gmx_bool bShowParameters)
410 if (available(fp, mtop, indent, title))
412 indent = pr_title(fp, indent, title);
413 pr_indent(fp, indent);
414 fprintf(fp, "name=\"%s\"\n", *(mtop->name));
415 pr_int(fp, indent, "#atoms", mtop->natoms);
416 pr_int(fp, indent, "#molblock", mtop->molblock.size());
417 for (size_t mb = 0; mb < mtop->molblock.size(); mb++)
419 pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
421 pr_str(fp, indent, "bIntermolecularInteractions",
422 gmx::boolToString(mtop->bIntermolecularInteractions));
423 if (mtop->bIntermolecularInteractions)
425 for (int j = 0; j < F_NRE; j++)
427 pr_ilist(fp, indent, interaction_function[j].longname,
428 mtop->ffparams.functype,
429 &mtop->intermolecular_ilist[j],
430 bShowNumbers, bShowParameters, mtop->ffparams.iparams);
433 pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
434 pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
435 for (size_t mt = 0; mt < mtop->moltype.size(); mt++)
437 pr_moltype(fp, indent, "moltype", &mtop->moltype[mt], mt,
438 &mtop->ffparams, bShowNumbers, bShowParameters);
440 pr_groups(fp, indent, &mtop->groups, bShowNumbers);
444 void pr_top(FILE *fp, int indent, const char *title, const t_topology *top,
445 gmx_bool bShowNumbers, gmx_bool bShowParameters)
447 if (available(fp, top, indent, title))
449 indent = pr_title(fp, indent, title);
450 pr_indent(fp, indent);
451 fprintf(fp, "name=\"%s\"\n", *(top->name));
452 pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers);
453 pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
454 pr_block(fp, indent, "cgs", &top->cgs, bShowNumbers);
455 pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
456 pr_str(fp, indent, "bIntermolecularInteractions",
457 gmx::boolToString(top->bIntermolecularInteractions));
458 pr_blocka(fp, indent, "excls", &top->excls, bShowNumbers);
459 pr_idef(fp, indent, "idef", &top->idef, bShowNumbers, bShowParameters);
463 static void cmp_ilist(FILE *fp, int ftype, const t_ilist *il1, const t_ilist *il2)
465 int i;
466 char buf[256];
468 fprintf(fp, "comparing ilist %s\n", interaction_function[ftype].name);
469 sprintf(buf, "%s->nr", interaction_function[ftype].name);
470 cmp_int(fp, buf, -1, il1->nr, il2->nr);
471 sprintf(buf, "%s->iatoms", interaction_function[ftype].name);
472 if (((il1->nr > 0) && (!il1->iatoms)) ||
473 ((il2->nr > 0) && (!il2->iatoms)) ||
474 ((il1->nr != il2->nr)))
476 fprintf(fp, "Comparing radically different topologies - %s is different\n",
477 buf);
479 else
481 for (i = 0; (i < il1->nr); i++)
483 cmp_int(fp, buf, i, il1->iatoms[i], il2->iatoms[i]);
488 static void cmp_iparm(FILE *fp, const char *s, t_functype ft,
489 const t_iparams &ip1, const t_iparams &ip2, real ftol, real abstol)
491 int i;
492 gmx_bool bDiff;
494 bDiff = FALSE;
495 for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
497 bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], ftol, abstol);
499 if (bDiff)
501 fprintf(fp, "%s1: ", s);
502 pr_iparams(fp, ft, &ip1);
503 fprintf(fp, "%s2: ", s);
504 pr_iparams(fp, ft, &ip2);
508 static void cmp_iparm_AB(FILE *fp, const char *s, t_functype ft,
509 const t_iparams &ip1, real ftol, real abstol)
511 int nrfpA, nrfpB, p0, i;
512 gmx_bool bDiff;
514 /* Normally the first parameter is perturbable */
515 p0 = 0;
516 nrfpA = interaction_function[ft].nrfpA;
517 nrfpB = interaction_function[ft].nrfpB;
518 if (ft == F_PDIHS)
520 nrfpB = 2;
522 else if (interaction_function[ft].flags & IF_TABULATED)
524 /* For tabulated interactions only the second parameter is perturbable */
525 p0 = 1;
526 nrfpB = 1;
528 bDiff = FALSE;
529 for (i = 0; i < nrfpB && !bDiff; i++)
531 bDiff = !equal_real(ip1.generic.buf[p0+i], ip1.generic.buf[nrfpA+i], ftol, abstol);
533 if (bDiff)
535 fprintf(fp, "%s: ", s);
536 pr_iparams(fp, ft, &ip1);
540 static void cmp_cmap(FILE *fp, const gmx_cmap_t *cmap1, const gmx_cmap_t *cmap2, real ftol, real abstol)
542 cmp_int(fp, "cmap ngrid", -1, cmap1->ngrid, cmap2->ngrid);
543 cmp_int(fp, "cmap grid_spacing", -1, cmap1->grid_spacing, cmap2->grid_spacing);
544 if (cmap1->ngrid == cmap2->ngrid &&
545 cmap1->grid_spacing == cmap2->grid_spacing)
547 int g;
549 for (g = 0; g < cmap1->ngrid; g++)
551 int i;
553 fprintf(fp, "comparing cmap %d\n", g);
555 for (i = 0; i < 4*cmap1->grid_spacing*cmap1->grid_spacing; i++)
557 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i], ftol, abstol);
563 static void cmp_idef(FILE *fp, const t_idef *id1, const t_idef *id2, real ftol, real abstol)
565 int i;
566 char buf1[64], buf2[64];
568 fprintf(fp, "comparing idef\n");
569 if (id2)
571 cmp_int(fp, "idef->ntypes", -1, id1->ntypes, id2->ntypes);
572 cmp_int(fp, "idef->atnr", -1, id1->atnr, id2->atnr);
573 for (i = 0; (i < std::min(id1->ntypes, id2->ntypes)); i++)
575 sprintf(buf1, "idef->functype[%d]", i);
576 sprintf(buf2, "idef->iparam[%d]", i);
577 cmp_int(fp, buf1, i, (int)id1->functype[i], (int)id2->functype[i]);
578 cmp_iparm(fp, buf2, id1->functype[i],
579 id1->iparams[i], id2->iparams[i], ftol, abstol);
581 cmp_real(fp, "fudgeQQ", -1, id1->fudgeQQ, id2->fudgeQQ, ftol, abstol);
582 cmp_cmap(fp, &id1->cmap_grid, &id2->cmap_grid, ftol, abstol);
583 for (i = 0; (i < F_NRE); i++)
585 cmp_ilist(fp, i, &(id1->il[i]), &(id2->il[i]));
588 else
590 for (i = 0; (i < id1->ntypes); i++)
592 cmp_iparm_AB(fp, "idef->iparam", id1->functype[i], id1->iparams[i], ftol, abstol);
597 static void cmp_block(FILE *fp, const t_block *b1, const t_block *b2, const char *s)
599 char buf[32];
601 fprintf(fp, "comparing block %s\n", s);
602 sprintf(buf, "%s.nr", s);
603 cmp_int(fp, buf, -1, b1->nr, b2->nr);
606 static void cmp_blocka(FILE *fp, const t_blocka *b1, const t_blocka *b2, const char *s)
608 char buf[32];
610 fprintf(fp, "comparing blocka %s\n", s);
611 sprintf(buf, "%s.nr", s);
612 cmp_int(fp, buf, -1, b1->nr, b2->nr);
613 sprintf(buf, "%s.nra", s);
614 cmp_int(fp, buf, -1, b1->nra, b2->nra);
617 void cmp_top(FILE *fp, const t_topology *t1, const t_topology *t2, real ftol, real abstol)
619 fprintf(fp, "comparing top\n");
620 if (t2)
622 cmp_idef(fp, &(t1->idef), &(t2->idef), ftol, abstol);
623 cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
624 cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
625 cmp_block(fp, &t1->mols, &t2->mols, "mols");
626 cmp_bool(fp, "bIntermolecularInteractions", -1, t1->bIntermolecularInteractions, t2->bIntermolecularInteractions);
627 cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
629 else
631 cmp_idef(fp, &(t1->idef), nullptr, ftol, abstol);
632 cmp_atoms(fp, &(t1->atoms), nullptr, ftol, abstol);
636 void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
637 int natoms0, int natoms1)
639 int i, j;
640 char buf[32];
642 fprintf(fp, "comparing groups\n");
644 for (i = 0; i < egcNR; i++)
646 sprintf(buf, "grps[%d].nr", i);
647 cmp_int(fp, buf, -1, g0->grps[i].nr, g1->grps[i].nr);
648 if (g0->grps[i].nr == g1->grps[i].nr)
650 for (j = 0; j < g0->grps[i].nr; j++)
652 sprintf(buf, "grps[%d].name[%d]", i, j);
653 cmp_str(fp, buf, -1,
654 *g0->grpname[g0->grps[i].nm_ind[j]],
655 *g1->grpname[g1->grps[i].nm_ind[j]]);
658 cmp_int(fp, "ngrpnr", i, g0->ngrpnr[i], g1->ngrpnr[i]);
659 if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
660 (g0->grpnr[i] != nullptr || g1->grpnr[i] != nullptr))
662 for (j = 0; j < natoms0; j++)
664 cmp_int(fp, gtypes[i], j, ggrpnr(g0, i, j), ggrpnr(g1, i, j));
668 /* We have compared the names in the groups lists,
669 * so we can skip the grpname list comparison.