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, 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", NULL
  58 };
  59
  60 static void init_groups(gmx_groups_t *groups)
  61 {
  62     groups->ngrpname = 0;
  63     groups->grpname  = NULL;
  64     for (int g = 0; g < egcNR; g++)
  65     {
  66         groups->grps[g].nm_ind = NULL;
  67         groups->ngrpnr[g]      = 0;
  68         groups->grpnr[g]       = NULL;
  69     }
  70
  71 }
  72
  73 void init_mtop(gmx_mtop_t *mtop)
  74 {
  75     mtop->name         = NULL;
  76     mtop->nmoltype     = 0;
  77     mtop->moltype      = NULL;
  78     mtop->nmolblock    = 0;
  79     mtop->molblock     = NULL;
  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 = NULL;
  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_mtop(gmx_mtop_t *mtop, gmx_bool bDoneSymtab)
 127 {
 128     if (bDoneSymtab)
 129     {
 130         done_symtab(&mtop->symtab);
 131     }
 132
 133     sfree(mtop->ffparams.functype);
 134     sfree(mtop->ffparams.iparams);
 135
 136     for (int i = 0; i < mtop->nmoltype; i++)
 137     {
 138         done_moltype(&mtop->moltype[i]);
 139     }
 140     sfree(mtop->moltype);
 141     for (int i = 0; i < mtop->nmolblock; i++)
 142     {
 143         done_molblock(&mtop->molblock[i]);
 144     }
 145     sfree(mtop->molblock);
 146     done_block(&mtop->mols);
 147 }
 148
 149 void done_top(t_topology *top)
 150 {
 151     sfree(top->idef.functype);
 152     sfree(top->idef.iparams);
 153     for (int f = 0; f < F_NRE; ++f)
 154     {
 155         sfree(top->idef.il[f].iatoms);
 156         top->idef.il[f].iatoms = NULL;
 157         top->idef.il[f].nalloc = 0;
 158     }
 159
 160     done_atom(&(top->atoms));
 161
 162     /* For GB */
 163     done_atomtypes(&(top->atomtypes));
 164
 165     done_symtab(&(top->symtab));
 166     done_block(&(top->cgs));
 167     done_block(&(top->mols));
 168     done_blocka(&(top->excls));
 169 }
 170
 171 bool gmx_mtop_has_masses(const gmx_mtop_t *mtop)
 172 {
 173     if (mtop == nullptr)
 174     {
 175         return false;
 176     }
 177     return mtop->nmoltype == 0 || mtop->moltype[0].atoms.haveMass;
 178 }
 179
 180 bool gmx_mtop_has_charges(const gmx_mtop_t *mtop)
 181 {
 182     if (mtop == nullptr)
 183     {
 184         return false;
 185     }
 186     return mtop->nmoltype == 0 || mtop->moltype[0].atoms.haveCharge;
 187 }
 188
 189 static void pr_grps(FILE *fp, const char *title, const t_grps grps[], char **grpname[])
 190 {
 191     int i, j;
 192
 193     for (i = 0; (i < egcNR); i++)
 194     {
 195         fprintf(fp, "%s[%-12s] nr=%d, name=[", title, gtypes[i], grps[i].nr);
 196         for (j = 0; (j < grps[i].nr); j++)
 197         {
 198             fprintf(fp, " %s", *(grpname[grps[i].nm_ind[j]]));
 199         }
 200         fprintf(fp, "]\n");
 201     }
 202 }
 203
 204 static void pr_groups(FILE *fp, int indent,
 205                       const gmx_groups_t *groups,
 206                       gmx_bool bShowNumbers)
 207 {
 208     int nat_max, i, g;
 209
 210     pr_grps(fp, "grp", groups->grps, groups->grpname);
 211     pr_strings(fp, indent, "grpname", groups->grpname, groups->ngrpname, bShowNumbers);
 212
 213     pr_indent(fp, indent);
 214     fprintf(fp, "groups          ");
 215     for (g = 0; g < egcNR; g++)
 216     {
 217         printf(" %5.5s", gtypes[g]);
 218     }
 219     printf("\n");
 220
 221     pr_indent(fp, indent);
 222     fprintf(fp, "allocated       ");
 223     nat_max = 0;
 224     for (g = 0; g < egcNR; g++)
 225     {
 226         printf(" %5d", groups->ngrpnr[g]);
 227         nat_max = std::max(nat_max, groups->ngrpnr[g]);
 228     }
 229     printf("\n");
 230
 231     if (nat_max == 0)
 232     {
 233         pr_indent(fp, indent);
 234         fprintf(fp, "groupnr[%5s] =", "*");
 235         for (g = 0; g < egcNR; g++)
 236         {
 237             fprintf(fp, "  %3d ", 0);
 238         }
 239         fprintf(fp, "\n");
 240     }
 241     else
 242     {
 243         for (i = 0; i < nat_max; i++)
 244         {
 245             pr_indent(fp, indent);
 246             fprintf(fp, "groupnr[%5d] =", i);
 247             for (g = 0; g < egcNR; g++)
 248             {
 249                 fprintf(fp, "  %3d ",
 250                         groups->grpnr[g] ? groups->grpnr[g][i] : 0);
 251             }
 252             fprintf(fp, "\n");
 253         }
 254     }
 255 }
 256
 257 static void pr_moltype(FILE *fp, int indent, const char *title,
 258                        const gmx_moltype_t *molt, int n,
 259                        const gmx_ffparams_t *ffparams,
 260                        gmx_bool bShowNumbers)
 261 {
 262     int j;
 263
 264     indent = pr_title_n(fp, indent, title, n);
 265     pr_indent(fp, indent);
 266     fprintf(fp, "name=\"%s\"\n", *(molt->name));
 267     pr_atoms(fp, indent, "atoms", &(molt->atoms), bShowNumbers);
 268     pr_block(fp, indent, "cgs", &molt->cgs, bShowNumbers);
 269     pr_blocka(fp, indent, "excls", &molt->excls, bShowNumbers);
 270     for (j = 0; (j < F_NRE); j++)
 271     {
 272         pr_ilist(fp, indent, interaction_function[j].longname,
 273                  ffparams->functype, &molt->ilist[j], bShowNumbers);
 274     }
 275 }
 276
 277 static void pr_molblock(FILE *fp, int indent, const char *title,
 278                         const gmx_molblock_t *molb, int n,
 279                         const gmx_moltype_t *molt)
 280 {
 281     indent = pr_title_n(fp, indent, title, n);
 282     pr_indent(fp, indent);
 283     fprintf(fp, "%-20s = %d \"%s\"\n",
 284             "moltype", molb->type, *(molt[molb->type].name));
 285     pr_int(fp, indent, "#molecules", molb->nmol);
 286     pr_int(fp, indent, "#atoms_mol", molb->natoms_mol);
 287     pr_int(fp, indent, "#posres_xA", molb->nposres_xA);
 288     if (molb->nposres_xA > 0)
 289     {
 290         pr_rvecs(fp, indent, "posres_xA", molb->posres_xA, molb->nposres_xA);
 291     }
 292     pr_int(fp, indent, "#posres_xB", molb->nposres_xB);
 293     if (molb->nposres_xB > 0)
 294     {
 295         pr_rvecs(fp, indent, "posres_xB", molb->posres_xB, molb->nposres_xB);
 296     }
 297 }
 298
 299 void pr_mtop(FILE *fp, int indent, const char *title, const gmx_mtop_t *mtop,
 300              gmx_bool bShowNumbers)
 301 {
 302     int mt, mb, j;
 303
 304     if (available(fp, mtop, indent, title))
 305     {
 306         indent = pr_title(fp, indent, title);
 307         pr_indent(fp, indent);
 308         fprintf(fp, "name=\"%s\"\n", *(mtop->name));
 309         pr_int(fp, indent, "#atoms", mtop->natoms);
 310         pr_int(fp, indent, "#molblock", mtop->nmolblock);
 311         for (mb = 0; mb < mtop->nmolblock; mb++)
 312         {
 313             pr_molblock(fp, indent, "molblock", &mtop->molblock[mb], mb, mtop->moltype);
 314         }
 315         pr_str(fp, indent, "bIntermolecularInteractions",
 316                gmx::boolToString(mtop->bIntermolecularInteractions));
 317         if (mtop->bIntermolecularInteractions)
 318         {
 319             for (j = 0; (j < F_NRE); j++)
 320             {
 321                 pr_ilist(fp, indent, interaction_function[j].longname,
 322                          mtop->ffparams.functype,
 323                          &mtop->intermolecular_ilist[j], bShowNumbers);
 324             }
 325         }
 326         pr_ffparams(fp, indent, "ffparams", &(mtop->ffparams), bShowNumbers);
 327         pr_atomtypes(fp, indent, "atomtypes", &(mtop->atomtypes), bShowNumbers);
 328         for (mt = 0; mt < mtop->nmoltype; mt++)
 329         {
 330             pr_moltype(fp, indent, "moltype", &mtop->moltype[mt], mt,
 331                        &mtop->ffparams, bShowNumbers);
 332         }
 333         pr_groups(fp, indent, &mtop->groups, bShowNumbers);
 334     }
 335 }
 336
 337 void pr_top(FILE *fp, int indent, const char *title, const t_topology *top, gmx_bool bShowNumbers)
 338 {
 339     if (available(fp, top, indent, title))
 340     {
 341         indent = pr_title(fp, indent, title);
 342         pr_indent(fp, indent);
 343         fprintf(fp, "name=\"%s\"\n", *(top->name));
 344         pr_atoms(fp, indent, "atoms", &(top->atoms), bShowNumbers);
 345         pr_atomtypes(fp, indent, "atomtypes", &(top->atomtypes), bShowNumbers);
 346         pr_block(fp, indent, "cgs", &top->cgs, bShowNumbers);
 347         pr_block(fp, indent, "mols", &top->mols, bShowNumbers);
 348         pr_str(fp, indent, "bIntermolecularInteractions",
 349                gmx::boolToString(top->bIntermolecularInteractions));
 350         pr_blocka(fp, indent, "excls", &top->excls, bShowNumbers);
 351         pr_idef(fp, indent, "idef", &top->idef, bShowNumbers);
 352     }
 353 }
 354
 355 static void cmp_ilist(FILE *fp, int ftype, const t_ilist *il1, const t_ilist *il2)
 356 {
 357     int  i;
 358     char buf[256];
 359
 360     fprintf(fp, "comparing ilist %s\n", interaction_function[ftype].name);
 361     sprintf(buf, "%s->nr", interaction_function[ftype].name);
 362     cmp_int(fp, buf, -1, il1->nr, il2->nr);
 363     sprintf(buf, "%s->iatoms", interaction_function[ftype].name);
 364     if (((il1->nr > 0) && (!il1->iatoms)) ||
 365         ((il2->nr > 0) && (!il2->iatoms)) ||
 366         ((il1->nr != il2->nr)))
 367     {
 368         fprintf(fp, "Comparing radically different topologies - %s is different\n",
 369                 buf);
 370     }
 371     else
 372     {
 373         for (i = 0; (i < il1->nr); i++)
 374         {
 375             cmp_int(fp, buf, i, il1->iatoms[i], il2->iatoms[i]);
 376         }
 377     }
 378 }
 379
 380 static void cmp_iparm(FILE *fp, const char *s, t_functype ft,
 381                       const t_iparams &ip1, const t_iparams &ip2, real ftol, real abstol)
 382 {
 383     int      i;
 384     gmx_bool bDiff;
 385
 386     bDiff = FALSE;
 387     for (i = 0; i < MAXFORCEPARAM && !bDiff; i++)
 388     {
 389         bDiff = !equal_real(ip1.generic.buf[i], ip2.generic.buf[i], ftol, abstol);
 390     }
 391     if (bDiff)
 392     {
 393         fprintf(fp, "%s1: ", s);
 394         pr_iparams(fp, ft, &ip1);
 395         fprintf(fp, "%s2: ", s);
 396         pr_iparams(fp, ft, &ip2);
 397     }
 398 }
 399
 400 static void cmp_iparm_AB(FILE *fp, const char *s, t_functype ft,
 401                          const t_iparams &ip1, real ftol, real abstol)
 402 {
 403     int      nrfpA, nrfpB, p0, i;
 404     gmx_bool bDiff;
 405
 406     /* Normally the first parameter is perturbable */
 407     p0    = 0;
 408     nrfpA = interaction_function[ft].nrfpA;
 409     nrfpB = interaction_function[ft].nrfpB;
 410     if (ft == F_PDIHS)
 411     {
 412         nrfpB = 2;
 413     }
 414     else if (interaction_function[ft].flags & IF_TABULATED)
 415     {
 416         /* For tabulated interactions only the second parameter is perturbable */
 417         p0    = 1;
 418         nrfpB = 1;
 419     }
 420     bDiff = FALSE;
 421     for (i = 0; i < nrfpB && !bDiff; i++)
 422     {
 423         bDiff = !equal_real(ip1.generic.buf[p0+i], ip1.generic.buf[nrfpA+i], ftol, abstol);
 424     }
 425     if (bDiff)
 426     {
 427         fprintf(fp, "%s: ", s);
 428         pr_iparams(fp, ft, &ip1);
 429     }
 430 }
 431
 432 static void cmp_cmap(FILE *fp, const gmx_cmap_t *cmap1, const gmx_cmap_t *cmap2, real ftol, real abstol)
 433 {
 434     cmp_int(fp, "cmap ngrid", -1, cmap1->ngrid, cmap2->ngrid);
 435     cmp_int(fp, "cmap grid_spacing", -1, cmap1->grid_spacing, cmap2->grid_spacing);
 436     if (cmap1->ngrid == cmap2->ngrid &&
 437         cmap1->grid_spacing == cmap2->grid_spacing)
 438     {
 439         int g;
 440
 441         for (g = 0; g < cmap1->ngrid; g++)
 442         {
 443             int i;
 444
 445             fprintf(fp, "comparing cmap %d\n", g);
 446
 447             for (i = 0; i < 4*cmap1->grid_spacing*cmap1->grid_spacing; i++)
 448             {
 449                 cmp_real(fp, "", i, cmap1->cmapdata[g].cmap[i], cmap2->cmapdata[g].cmap[i], ftol, abstol);
 450             }
 451         }
 452     }
 453 }
 454
 455 static void cmp_idef(FILE *fp, const t_idef *id1, const t_idef *id2, real ftol, real abstol)
 456 {
 457     int  i;
 458     char buf1[64], buf2[64];
 459
 460     fprintf(fp, "comparing idef\n");
 461     if (id2)
 462     {
 463         cmp_int(fp, "idef->ntypes", -1, id1->ntypes, id2->ntypes);
 464         cmp_int(fp, "idef->atnr",  -1, id1->atnr, id2->atnr);
 465         for (i = 0; (i < std::min(id1->ntypes, id2->ntypes)); i++)
 466         {
 467             sprintf(buf1, "idef->functype[%d]", i);
 468             sprintf(buf2, "idef->iparam[%d]", i);
 469             cmp_int(fp, buf1, i, (int)id1->functype[i], (int)id2->functype[i]);
 470             cmp_iparm(fp, buf2, id1->functype[i],
 471                       id1->iparams[i], id2->iparams[i], ftol, abstol);
 472         }
 473         cmp_real(fp, "fudgeQQ", -1, id1->fudgeQQ, id2->fudgeQQ, ftol, abstol);
 474         cmp_cmap(fp, &id1->cmap_grid, &id2->cmap_grid, ftol, abstol);
 475         for (i = 0; (i < F_NRE); i++)
 476         {
 477             cmp_ilist(fp, i, &(id1->il[i]), &(id2->il[i]));
 478         }
 479     }
 480     else
 481     {
 482         for (i = 0; (i < id1->ntypes); i++)
 483         {
 484             cmp_iparm_AB(fp, "idef->iparam", id1->functype[i], id1->iparams[i], ftol, abstol);
 485         }
 486     }
 487 }
 488
 489 static void cmp_block(FILE *fp, const t_block *b1, const t_block *b2, const char *s)
 490 {
 491     char buf[32];
 492
 493     fprintf(fp, "comparing block %s\n", s);
 494     sprintf(buf, "%s.nr", s);
 495     cmp_int(fp, buf, -1, b1->nr, b2->nr);
 496 }
 497
 498 static void cmp_blocka(FILE *fp, const t_blocka *b1, const t_blocka *b2, const char *s)
 499 {
 500     char buf[32];
 501
 502     fprintf(fp, "comparing blocka %s\n", s);
 503     sprintf(buf, "%s.nr", s);
 504     cmp_int(fp, buf, -1, b1->nr, b2->nr);
 505     sprintf(buf, "%s.nra", s);
 506     cmp_int(fp, buf, -1, b1->nra, b2->nra);
 507 }
 508
 509 void cmp_top(FILE *fp, const t_topology *t1, const t_topology *t2, real ftol, real abstol)
 510 {
 511     fprintf(fp, "comparing top\n");
 512     if (t2)
 513     {
 514         cmp_idef(fp, &(t1->idef), &(t2->idef), ftol, abstol);
 515         cmp_atoms(fp, &(t1->atoms), &(t2->atoms), ftol, abstol);
 516         cmp_block(fp, &t1->cgs, &t2->cgs, "cgs");
 517         cmp_block(fp, &t1->mols, &t2->mols, "mols");
 518         cmp_bool(fp, "bIntermolecularInteractions", -1, t1->bIntermolecularInteractions, t2->bIntermolecularInteractions);
 519         cmp_blocka(fp, &t1->excls, &t2->excls, "excls");
 520     }
 521     else
 522     {
 523         cmp_idef(fp, &(t1->idef), NULL, ftol, abstol);
 524         cmp_atoms(fp, &(t1->atoms), NULL, ftol, abstol);
 525     }
 526 }
 527
 528 void cmp_groups(FILE *fp, const gmx_groups_t *g0, const gmx_groups_t *g1,
 529                 int natoms0, int natoms1)
 530 {
 531     int  i, j;
 532     char buf[32];
 533
 534     fprintf(fp, "comparing groups\n");
 535
 536     for (i = 0; i < egcNR; i++)
 537     {
 538         sprintf(buf, "grps[%d].nr", i);
 539         cmp_int(fp, buf, -1, g0->grps[i].nr, g1->grps[i].nr);
 540         if (g0->grps[i].nr == g1->grps[i].nr)
 541         {
 542             for (j = 0; j < g0->grps[i].nr; j++)
 543             {
 544                 sprintf(buf, "grps[%d].name[%d]", i, j);
 545                 cmp_str(fp, buf, -1,
 546                         *g0->grpname[g0->grps[i].nm_ind[j]],
 547                         *g1->grpname[g1->grps[i].nm_ind[j]]);
 548             }
 549         }
 550         cmp_int(fp, "ngrpnr", i, g0->ngrpnr[i], g1->ngrpnr[i]);
 551         if (g0->ngrpnr[i] == g1->ngrpnr[i] && natoms0 == natoms1 &&
 552             (g0->grpnr[i] != NULL || g1->grpnr[i] != NULL))
 553         {
 554             for (j = 0; j < natoms0; j++)
 555             {
 556                 cmp_int(fp, gtypes[i], j, ggrpnr(g0, i, j), ggrpnr(g1, i, j));
 557             }
 558         }
 559     }
 560     /* We have compared the names in the groups lists,
 561      * so we can skip the grpname list comparison.
 562      */
 563 }