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