Tidy: modernize-use-nullptr
[alexxy/gromacs.git] / src / gromacs / topology / topology.cpp
1 /*
2  * This file is part of the GROMACS molecular simulation package.
3  *
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, 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.
10  *
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.
15  *
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.
20  *
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.
25  *
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.
33  *
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.
36  */
37 #include "gmxpre.h"
38
39 #include "topology.h"
40
41 #include <cstdio>
42
43 #include <algorithm>
44
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"
54
55 const char *gtypes[egcNR+1] = {
56     "T-Coupling", "Energy Mon.", "Acceleration", "Freeze",
57     "User1", "User2", "VCM", "Compressed X", "Or. Res. Fit", "QMMM", nullptr
58 };
59
60 static void init_groups(gmx_groups_t *groups)
61 {
62     groups->ngrpname = 0;
63     groups->grpname  = nullptr;
64     for (int g = 0; g < egcNR; g++)
65     {
66         groups->grps[g].nm_ind = nullptr;
67         groups->ngrpnr[g]      = 0;
68         groups->grpnr[g]       = nullptr;
69     }
70
71 }
72
73 void init_mtop(gmx_mtop_t *mtop)
74 {
75     mtop->name         = nullptr;
76     mtop->nmoltype     = 0;
77     mtop->moltype      = nullptr;
78     mtop->nmolblock    = 0;
79     mtop->molblock     = nullptr;
80     mtop->maxres_renum = 0;
81     mtop->maxresnr     = -1;
82     init_groups(&mtop->groups);
83     init_block(&mtop->mols);
84     open_symtab(&mtop->symtab);
85 }
86
87 void init_top(t_topology *top)
88 {
89     top->name = nullptr;
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);
96 }
97
98
99 void done_moltype(gmx_moltype_t *molt)
100 {
101     done_atom(&molt->atoms);
102     done_block(&molt->cgs);
103     done_blocka(&molt->excls);
104
105     for (int f = 0; f < F_NRE; f++)
106     {
107         sfree(molt->ilist[f].iatoms);
108         molt->ilist[f].nalloc = 0;
109     }
110 }
111
112 void done_molblock(gmx_molblock_t *molb)
113 {
114     if (molb->nposres_xA > 0)
115     {
116         molb->nposres_xA = 0;
117         sfree(molb->posres_xA);
118     }
119     if (molb->nposres_xB > 0)
120     {
121         molb->nposres_xB = 0;
122         sfree(molb->posres_xB);
123     }
124 }
125
126 void done_gmx_groups_t(gmx_groups_t *g)
127 {
128     int i;
129
130     for (i = 0; (i < egcNR); i++)
131     {
132         if (nullptr != g->grps[i].nm_ind)
133         {
134             sfree(g->grps[i].nm_ind);
135             g->grps[i].nm_ind = nullptr;
136         }
137         if (nullptr != g->grpnr[i])
138         {
139             sfree(g->grpnr[i]);
140             g->grpnr[i] = nullptr;
141         }
142     }
143     /* The contents of this array is in symtab, don't free it here */
144     sfree(g->grpname);
145 }
146
147 void done_mtop(gmx_mtop_t *mtop)
148 {
149     done_symtab(&mtop->symtab);
150
151     sfree(mtop->ffparams.functype);
152     sfree(mtop->ffparams.iparams);
153
154     for (int i = 0; i < mtop->nmoltype; i++)
155     {
156         done_moltype(&mtop->moltype[i]);
157     }
158     sfree(mtop->moltype);
159     for (int i = 0; i < mtop->nmolblock; i++)
160     {
161         done_molblock(&mtop->molblock[i]);
162     }
163     sfree(mtop->molblock);
164     done_atomtypes(&mtop->atomtypes);
165     done_gmx_groups_t(&mtop->groups);
166     done_block(&mtop->mols);
167 }
168
169 void done_top(t_topology *top)
170 {
171     sfree(top->idef.functype);
172     sfree(top->idef.iparams);
173     for (int f = 0; f < F_NRE; ++f)
174     {
175         sfree(top->idef.il[f].iatoms);
176         top->idef.il[f].iatoms = nullptr;
177         top->idef.il[f].nalloc = 0;
178     }
179
180     done_atom(&(top->atoms));
181
182     /* For GB */
183     done_atomtypes(&(top->atomtypes));
184
185     done_symtab(&(top->symtab));
186     done_block(&(top->cgs));
187     done_block(&(top->mols));
188     done_blocka(&(top->excls));
189 }
190
191 void done_top_mtop(t_topology *top, gmx_mtop_t *mtop)
192 {
193     if (mtop != nullptr)
194     {
195         if (top != nullptr)
196         {
197             for (int f = 0; f < F_NRE; ++f)
198             {
199                 sfree(top->idef.il[f].iatoms);
200                 top->idef.il[f].iatoms = nullptr;
201                 top->idef.il[f].nalloc = 0;
202             }
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);
208         }
209         done_mtop(mtop);
210     }
211 }
212
213 bool gmx_mtop_has_masses(const gmx_mtop_t *mtop)
214 {
215     if (mtop == nullptr)
216     {
217         return false;
218     }
219     return mtop->nmoltype == 0 || mtop->moltype[0].atoms.haveMass;
220 }
221
222 bool gmx_mtop_has_charges(const gmx_mtop_t *mtop)
223 {
224     if (mtop == nullptr)
225     {
226         return false;
227     }
228     return mtop->nmoltype == 0 || mtop->moltype[0].atoms.haveCharge;
229 }
230
231 bool gmx_mtop_has_atomtypes(const gmx_mtop_t *mtop)
232 {
233     if (mtop == nullptr)
234     {
235         return false;
236     }
237     return mtop->nmoltype == 0 || mtop->moltype[0].atoms.haveType;
238 }
239
240 bool gmx_mtop_has_pdbinfo(const gmx_mtop_t *mtop)
241 {
242     if (mtop == nullptr)
243     {
244         return false;
245     }
246     return mtop->nmoltype == 0 || mtop->moltype[0].atoms.havePdbInfo;
247 }
248
249 static void pr_grps(FILE *fp, const char *title, const t_grps grps[], char **grpname[])
250 {
251     int i, j;
252
253     for (i = 0; (i < egcNR); i++)
254     {
255         fprintf(fp, "%s[%-12s] nr=%d, name=[", title, gtypes[i], grps[i].nr);
256         for (j = 0; (j < grps[i].nr); j++)
257         {
258             fprintf(fp, " %s", *(grpname[grps[i].nm_ind[j]]));
259         }
260         fprintf(fp, "]\n");
261     }
262 }
263
264 static void pr_groups(FILE *fp, int indent,
265                       const gmx_groups_t *groups,
266                       gmx_bool bShowNumbers)
267 {
268     int nat_max, i, g;
269
270     pr_grps(fp, "grp", groups->grps, groups->grpname);
271     pr_strings(fp, indent, "grpname", groups->grpname, groups->ngrpname, bShowNumbers);
272
273     pr_indent(fp, indent);
274     fprintf(fp, "groups          ");
275     for (g = 0; g < egcNR; g++)
276     {
277         printf(" %5.5s", gtypes[g]);
278     }
279     printf("\n");
280
281     pr_indent(fp, indent);
282     fprintf(fp, "allocated       ");
283     nat_max = 0;
284     for (g = 0; g < egcNR; g++)
285     {
286         printf(" %5d", groups->ngrpnr[g]);
287         nat_max = std::max(nat_max, groups->ngrpnr[g]);
288     }
289     printf("\n");
290
291     if (nat_max == 0)
292     {
293         pr_indent(fp, indent);
294         fprintf(fp, "groupnr[%5s] =", "*");
295         for (g = 0; g < egcNR; g++)
296         {
297             fprintf(fp, "  %3d ", 0);
298         }
299         fprintf(fp, "\n");
300     }
301     else
302     {
303         for (i = 0; i < nat_max; i++)
304         {
305             pr_indent(fp, indent);
306             fprintf(fp, "groupnr[%5d] =", i);
307             for (g = 0; g < egcNR; g++)
308             {
309                 fprintf(fp, "  %3d ",
310                         groups->grpnr[g] ? groups->grpnr[g][i] : 0);
311             }
312             fprintf(fp, "\n");
313         }
314     }
315 }
316
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, gmx_bool bShowParameters)
321 {
322     int j;
323
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++)
331     {
332         pr_ilist(fp, indent, interaction_function[j].longname,
333                  ffparams->functype, &molt->ilist[j],
334                  bShowNumbers, bShowParameters, ffparams->iparams);
335     }
336 }
337
338 static void pr_molblock(FILE *fp, int indent, const char *title,
339                         const gmx_molblock_t *molb, int n,
340                         const gmx_moltype_t *molt)
341 {
342     indent = pr_title_n(fp, indent, title, n);
343     pr_indent(fp, indent);
344     fprintf(fp, "%-20s = %d \"%s\"\n",
345             "moltype", molb->type, *(molt[molb->type].name));
346     pr_int(fp, indent, "#molecules", molb->nmol);
347     pr_int(fp, indent, "#atoms_mol", molb->natoms_mol);
348     pr_int(fp, indent, "#posres_xA", molb->nposres_xA);
349     if (molb->nposres_xA > 0)
350     {
351         pr_rvecs(fp, indent, "posres_xA", molb->posres_xA, molb->nposres_xA);
352     }
353     pr_int(fp, indent, "#posres_xB", molb->nposres_xB);
354     if (molb->nposres_xB > 0)
355     {
356         pr_rvecs(fp, indent, "posres_xB", molb->posres_xB, molb->nposres_xB);
357     }
358 }
359
360 void pr_mtop(FILE *fp, int indent, const char *title, const gmx_mtop_t *mtop,
361              gmx_bool bShowNumbers, gmx_bool bShowParameters)
362 {
363     int mt, mb, j;
364
365     if (available(fp, mtop, indent, title))
366     {
367         indent = pr_title(fp, indent, title);
368         pr_indent(fp, indent);
369         fprintf(fp, "name=\"%s\"\n", *(mtop->name));
370         pr_int(fp, indent, "#atoms", mtop->natoms);
371         pr_int(fp, indent, "#molblock", mtop->nmolblock);
372         for (mb = 0; mb < mtop->nmolblock; mb++)
373         {
374             pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
375         }
376         pr_str(fp, indent, "bIntermolecularInteractions",
377                gmx::boolToString(mtop->bIntermolecularInteractions));
378         if (mtop->bIntermolecularInteractions)
379         {
380             for (j = 0; (j < F_NRE); j++)
381             {
382                 pr_ilist(fp, indent, interaction_function[j].longname,
383                          mtop->ffparams.functype,
384                          &mtop->intermolecular_ilist[j],
385                          bShowNumbers, bShowParameters, mtop->ffparams.iparams);
386             }
387         }
388         pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
389         pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
390         for (mt = 0; mt < mtop->nmoltype; mt++)
391         {
392             pr_moltype(fp, indent, "moltype", &mtop->moltype[mt], mt,
393                        &mtop->ffparams, bShowNumbers, bShowParameters);
394         }
395         pr_groups(fp, indent, &mtop->groups, bShowNumbers);
396     }
397 }
398
399 void pr_top(FILE *fp, int indent, const char *title, const t_topology *top,
400             gmx_bool bShowNumbers, gmx_bool bShowParameters)
401 {
402     if (available(fp, top, indent, title))
403     {
404         indent = pr_title(fp, indent, title);
405         pr_indent(fp, indent);
406         fprintf(fp, "name=\"%s\"\n", *(top->name));
407         pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers);
408         pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
409         pr_block(fp, indent, "cgs", &top->cgs, bShowNumbers);
410         pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
411         pr_str(fp, indent, "bIntermolecularInteractions",
412                gmx::boolToString(top->bIntermolecularInteractions));
413         pr_blocka(fp, indent, "excls", &top->excls, bShowNumbers);
414         pr_idef(fp, indent, "idef", &top->idef, bShowNumbers, bShowParameters);
415     }
416 }
417
418 static void cmp_ilist(FILE *fp, int ftype, const t_ilist *il1, const t_ilist *il2)
419 {
420     int  i;
421     char buf[256];
422
423     fprintf(fp, "comparing ilist %s\n", interaction_function[ftype].name);
424     sprintf(buf, "%s->nr", interaction_function[ftype].name);
425     cmp_int(fp, buf, -1, il1->nr, il2->nr);
426     sprintf(buf, "%s->iatoms", interaction_function[ftype].name);
427     if (((il1->nr > 0) && (!il1->iatoms)) ||
428         ((il2->nr > 0) && (!il2->iatoms)) ||
429         ((il1->nr != il2->nr)))
430     {
431         fprintf(fp, "Comparing radically different topologies - %s is different\n",
432                 buf);
433     }
434     else
435     {
436         for (i = 0; (i < il1->nr); i++)
437         {
438             cmp_int(fp, buf, i, il1->iatoms[i], il2->iatoms[i]);
439         }
440     }
441 }
442
443 static void cmp_iparm(FILE *fp, const char *s, t_functype ft,
444                       const t_iparams &ip1, const t_iparams &ip2, real ftol, real abstol)
445 {
446     int      i;
447     gmx_bool bDiff;
448
449     bDiff = FALSE;
450     for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
451     {
452         bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], ftol, abstol);
453     }
454     if (bDiff)
455     {
456         fprintf(fp, "%s1: ", s);
457         pr_iparams(fp, ft, &ip1);
458         fprintf(fp, "%s2: ", s);
459         pr_iparams(fp, ft, &ip2);
460     }
461 }
462
463 static void cmp_iparm_AB(FILE *fp, const char *s, t_functype ft,
464                          const t_iparams &ip1, real ftol, real abstol)
465 {
466     int      nrfpA, nrfpB, p0, i;
467     gmx_bool bDiff;
468
469     /* Normally the first parameter is perturbable */
470     p0    = 0;
471     nrfpA = interaction_function[ft].nrfpA;
472     nrfpB = interaction_function[ft].nrfpB;
473     if (ft == F_PDIHS)
474     {
475         nrfpB = 2;
476     }
477     else if (interaction_function[ft].flags & IF_TABULATED)
478     {
479         /* For tabulated interactions only the second parameter is perturbable */
480         p0    = 1;
481         nrfpB = 1;
482     }
483     bDiff = FALSE;
484     for (i = 0; i < nrfpB && !bDiff; i++)
485     {
486         bDiff = !equal_real(ip1.generic.buf[p0+i], ip1.generic.buf[nrfpA+i], ftol, abstol);
487     }
488     if (bDiff)
489     {
490         fprintf(fp, "%s: ", s);
491         pr_iparams(fp, ft, &ip1);
492     }
493 }
494
495 static void cmp_cmap(FILE *fp, const gmx_cmap_t *cmap1, const gmx_cmap_t *cmap2, real ftol, real abstol)
496 {
497     cmp_int(fp, "cmap ngrid", -1, cmap1->ngrid, cmap2->ngrid);
498     cmp_int(fp, "cmap grid_spacing", -1, cmap1->grid_spacing, cmap2->grid_spacing);
499     if (cmap1->ngrid == cmap2->ngrid &&
500         cmap1->grid_spacing == cmap2->grid_spacing)
501     {
502         int g;
503
504         for (g = 0; g < cmap1->ngrid; g++)
505         {
506             int i;
507
508             fprintf(fp, "comparing cmap %d\n", g);
509
510             for (i = 0; i < 4*cmap1->grid_spacing*cmap1->grid_spacing; i++)
511             {
512                 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i], ftol, abstol);
513             }
514         }
515     }
516 }
517
518 static void cmp_idef(FILE *fp, const t_idef *id1, const t_idef *id2, real ftol, real abstol)
519 {
520     int  i;
521     char buf1[64], buf2[64];
522
523     fprintf(fp, "comparing idef\n");
524     if (id2)
525     {
526         cmp_int(fp, "idef->ntypes", -1, id1->ntypes, id2->ntypes);
527         cmp_int(fp, "idef->atnr",  -1, id1->atnr, id2->atnr);
528         for (i = 0; (i < std::min(id1->ntypes, id2->ntypes)); i++)
529         {
530             sprintf(buf1, "idef->functype[%d]", i);
531             sprintf(buf2, "idef->iparam[%d]", i);
532             cmp_int(fp, buf1, i, (int)id1->functype[i], (int)id2->functype[i]);
533             cmp_iparm(fp, buf2, id1->functype[i],
534                       id1->iparams[i], id2->iparams[i], ftol, abstol);
535         }
536         cmp_real(fp, "fudgeQQ", -1, id1->fudgeQQ, id2->fudgeQQ, ftol, abstol);
537         cmp_cmap(fp, &id1->cmap_grid, &id2->cmap_grid, ftol, abstol);
538         for (i = 0; (i < F_NRE); i++)
539         {
540             cmp_ilist(fp, i, &(id1->il[i]), &(id2->il[i]));
541         }
542     }
543     else
544     {
545         for (i = 0; (i < id1->ntypes); i++)
546         {
547             cmp_iparm_AB(fp, "idef->iparam", id1->functype[i], id1->iparams[i], ftol, abstol);
548         }
549     }
550 }
551
552 static void cmp_block(FILE *fp, const t_block *b1, const t_block *b2, const char *s)
553 {
554     char buf[32];
555
556     fprintf(fp, "comparing block %s\n", s);
557     sprintf(buf, "%s.nr", s);
558     cmp_int(fp, buf, -1, b1->nr, b2->nr);
559 }
560
561 static void cmp_blocka(FILE *fp, const t_blocka *b1, const t_blocka *b2, const char *s)
562 {
563     char buf[32];
564
565     fprintf(fp, "comparing blocka %s\n", s);
566     sprintf(buf, "%s.nr", s);
567     cmp_int(fp, buf, -1, b1->nr, b2->nr);
568     sprintf(buf, "%s.nra", s);
569     cmp_int(fp, buf, -1, b1->nra, b2->nra);
570 }
571
572 void cmp_top(FILE *fp, const t_topology *t1, const t_topology *t2, real ftol, real abstol)
573 {
574     fprintf(fp, "comparing top\n");
575     if (t2)
576     {
577         cmp_idef(fp, &(t1->idef), &(t2->idef), ftol, abstol);
578         cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
579         cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
580         cmp_block(fp, &t1->mols, &t2->mols, "mols");
581         cmp_bool(fp, "bIntermolecularInteractions", -1, t1->bIntermolecularInteractions, t2->bIntermolecularInteractions);
582         cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
583     }
584     else
585     {
586         cmp_idef(fp, &(t1->idef), nullptr, ftol, abstol);
587         cmp_atoms(fp, &(t1->atoms), nullptr, ftol, abstol);
588     }
589 }
590
591 void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
592                 int natoms0, int natoms1)
593 {
594     int  i, j;
595     char buf[32];
596
597     fprintf(fp, "comparing groups\n");
598
599     for (i = 0; i < egcNR; i++)
600     {
601         sprintf(buf, "grps[%d].nr", i);
602         cmp_int(fp, buf, -1, g0->grps[i].nr, g1->grps[i].nr);
603         if (g0->grps[i].nr == g1->grps[i].nr)
604         {
605             for (j = 0; j < g0->grps[i].nr; j++)
606             {
607                 sprintf(buf, "grps[%d].name[%d]", i, j);
608                 cmp_str(fp, buf, -1,
609                         *g0->grpname[g0->grps[i].nm_ind[j]],
610                         *g1->grpname[g1->grps[i].nm_ind[j]]);
611             }
612         }
613         cmp_int(fp, "ngrpnr", i, g0->ngrpnr[i], g1->ngrpnr[i]);
614         if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
615             (g0->grpnr[i] != nullptr || g1->grpnr[i] != nullptr))
616         {
617             for (j = 0; j < natoms0; j++)
618             {
619                 cmp_int(fp, gtypes[i], j, ggrpnr(g0, i, j), ggrpnr(g1, i, j));
620             }
621         }
622     }
623     /* We have compared the names in the groups lists,
624      * so we can skip the grpname list comparison.
625      */
626 }