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