src/gromacs/topology/idef.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,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 "idef.h"
  40
  41 #include <cstdio>
  42
  43 #include "gromacs/topology/ifunc.h"
  44 #include "gromacs/utility/fatalerror.h"
  45 #include "gromacs/utility/smalloc.h"
  46 #include "gromacs/utility/txtdump.h"
  47
  48 static void pr_harm(FILE *fp, const t_iparams *iparams, const char *r, const char *kr)
  49 {
  50     fprintf(fp, "%sA=%12.5e, %sA=%12.5e, %sB=%12.5e, %sB=%12.5e\n",
  51             r, iparams->harmonic.rA, kr, iparams->harmonic.krA,
  52             r, iparams->harmonic.rB, kr, iparams->harmonic.krB);
  53 }
  54
  55 void pr_iparams(FILE *fp, t_functype ftype, const t_iparams *iparams)
  56 {
  57     switch (ftype)
  58     {
  59         case F_ANGLES:
  60         case F_G96ANGLES:
  61             pr_harm(fp, iparams, "th", "ct");
  62             break;
  63         case F_CROSS_BOND_BONDS:
  64             fprintf(fp, "r1e=%15.8e, r2e=%15.8e, krr=%15.8e\n",
  65                     iparams->cross_bb.r1e, iparams->cross_bb.r2e,
  66                     iparams->cross_bb.krr);
  67             break;
  68         case F_CROSS_BOND_ANGLES:
  69             fprintf(fp, "r1e=%15.8e, r1e=%15.8e, r3e=%15.8e, krt=%15.8e\n",
  70                     iparams->cross_ba.r1e, iparams->cross_ba.r2e,
  71                     iparams->cross_ba.r3e, iparams->cross_ba.krt);
  72             break;
  73         case F_LINEAR_ANGLES:
  74             fprintf(fp, "klinA=%15.8e, aA=%15.8e, klinB=%15.8e, aB=%15.8e\n",
  75                     iparams->linangle.klinA, iparams->linangle.aA,
  76                     iparams->linangle.klinB, iparams->linangle.aB);
  77             break;
  78         case F_UREY_BRADLEY:
  79             fprintf(fp, "thetaA=%15.8e, kthetaA=%15.8e, r13A=%15.8e, kUBA=%15.8e, thetaB=%15.8e, kthetaB=%15.8e, r13B=%15.8e, kUBB=%15.8e\n", iparams->u_b.thetaA, iparams->u_b.kthetaA, iparams->u_b.r13A, iparams->u_b.kUBA, iparams->u_b.thetaB, iparams->u_b.kthetaB, iparams->u_b.r13B, iparams->u_b.kUBB);
  80             break;
  81         case F_QUARTIC_ANGLES:
  82             fprintf(fp, "theta=%15.8e", iparams->qangle.theta);
  83             for (int i = 0; i < 5; i++)
  84             {
  85                 fprintf(fp, ", c%c=%15.8e", '0'+i, iparams->qangle.c[i]);
  86             }
  87             fprintf(fp, "\n");
  88             break;
  89         case F_BHAM:
  90             fprintf(fp, "a=%15.8e, b=%15.8e, c=%15.8e\n",
  91                     iparams->bham.a, iparams->bham.b, iparams->bham.c);
  92             break;
  93         case F_BONDS:
  94         case F_G96BONDS:
  95         case F_HARMONIC:
  96             pr_harm(fp, iparams, "b0", "cb");
  97             break;
  98         case F_IDIHS:
  99             pr_harm(fp, iparams, "xi", "cx");
 100             break;
 101         case F_MORSE:
 102             fprintf(fp, "b0A=%15.8e, cbA=%15.8e, betaA=%15.8e, b0B=%15.8e, cbB=%15.8e, betaB=%15.8e\n",
 103                     iparams->morse.b0A, iparams->morse.cbA, iparams->morse.betaA,
 104                     iparams->morse.b0B, iparams->morse.cbB, iparams->morse.betaB);
 105             break;
 106         case F_CUBICBONDS:
 107             fprintf(fp, "b0=%15.8e, kb=%15.8e, kcub=%15.8e\n",
 108                     iparams->cubic.b0, iparams->cubic.kb, iparams->cubic.kcub);
 109             break;
 110         case F_CONNBONDS:
 111             fprintf(fp, "\n");
 112             break;
 113         case F_FENEBONDS:
 114             fprintf(fp, "bm=%15.8e, kb=%15.8e\n", iparams->fene.bm, iparams->fene.kb);
 115             break;
 116         case F_RESTRBONDS:
 117             fprintf(fp, "lowA=%15.8e, up1A=%15.8e, up2A=%15.8e, kA=%15.8e, lowB=%15.8e, up1B=%15.8e, up2B=%15.8e, kB=%15.8e,\n",
 118                     iparams->restraint.lowA, iparams->restraint.up1A,
 119                     iparams->restraint.up2A, iparams->restraint.kA,
 120                     iparams->restraint.lowB, iparams->restraint.up1B,
 121                     iparams->restraint.up2B, iparams->restraint.kB);
 122             break;
 123         case F_TABBONDS:
 124         case F_TABBONDSNC:
 125         case F_TABANGLES:
 126         case F_TABDIHS:
 127             fprintf(fp, "tab=%d, kA=%15.8e, kB=%15.8e\n",
 128                     iparams->tab.table, iparams->tab.kA, iparams->tab.kB);
 129             break;
 130         case F_POLARIZATION:
 131             fprintf(fp, "alpha=%15.8e\n", iparams->polarize.alpha);
 132             break;
 133         case F_ANHARM_POL:
 134             fprintf(fp, "alpha=%15.8e drcut=%15.8e khyp=%15.8e\n",
 135                     iparams->anharm_polarize.alpha,
 136                     iparams->anharm_polarize.drcut,
 137                     iparams->anharm_polarize.khyp);
 138             break;
 139         case F_THOLE_POL:
 140             fprintf(fp, "a=%15.8e, alpha1=%15.8e, alpha2=%15.8e, rfac=%15.8e\n",
 141                     iparams->thole.a, iparams->thole.alpha1, iparams->thole.alpha2,
 142                     iparams->thole.rfac);
 143             break;
 144         case F_WATER_POL:
 145             fprintf(fp, "al_x=%15.8e, al_y=%15.8e, al_z=%15.8e, rOH=%9.6f, rHH=%9.6f, rOD=%9.6f\n",
 146                     iparams->wpol.al_x, iparams->wpol.al_y, iparams->wpol.al_z,
 147                     iparams->wpol.rOH, iparams->wpol.rHH, iparams->wpol.rOD);
 148             break;
 149         case F_LJ:
 150             fprintf(fp, "c6=%15.8e, c12=%15.8e\n", iparams->lj.c6, iparams->lj.c12);
 151             break;
 152         case F_LJ14:
 153             fprintf(fp, "c6A=%15.8e, c12A=%15.8e, c6B=%15.8e, c12B=%15.8e\n",
 154                     iparams->lj14.c6A, iparams->lj14.c12A,
 155                     iparams->lj14.c6B, iparams->lj14.c12B);
 156             break;
 157         case F_LJC14_Q:
 158             fprintf(fp, "fqq=%15.8e, qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e\n",
 159                     iparams->ljc14.fqq,
 160                     iparams->ljc14.qi, iparams->ljc14.qj,
 161                     iparams->ljc14.c6, iparams->ljc14.c12);
 162             break;
 163         case F_LJC_PAIRS_NB:
 164             fprintf(fp, "qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e\n",
 165                     iparams->ljcnb.qi, iparams->ljcnb.qj,
 166                     iparams->ljcnb.c6, iparams->ljcnb.c12);
 167             break;
 168         case F_PDIHS:
 169         case F_PIDIHS:
 170         case F_ANGRES:
 171         case F_ANGRESZ:
 172             fprintf(fp, "phiA=%15.8e, cpA=%15.8e, phiB=%15.8e, cpB=%15.8e, mult=%d\n",
 173                     iparams->pdihs.phiA, iparams->pdihs.cpA,
 174                     iparams->pdihs.phiB, iparams->pdihs.cpB,
 175                     iparams->pdihs.mult);
 176             break;
 177         case F_DISRES:
 178             fprintf(fp, "label=%4d, type=%1d, low=%15.8e, up1=%15.8e, up2=%15.8e, fac=%15.8e)\n",
 179                     iparams->disres.label, iparams->disres.type,
 180                     iparams->disres.low, iparams->disres.up1,
 181                     iparams->disres.up2, iparams->disres.kfac);
 182             break;
 183         case F_ORIRES:
 184             fprintf(fp, "ex=%4d, label=%d, power=%4d, c=%15.8e, obs=%15.8e, kfac=%15.8e)\n",
 185                     iparams->orires.ex, iparams->orires.label, iparams->orires.power,
 186                     iparams->orires.c, iparams->orires.obs, iparams->orires.kfac);
 187             break;
 188         case F_DIHRES:
 189             fprintf(fp, "phiA=%15.8e, dphiA=%15.8e, kfacA=%15.8e, phiB=%15.8e, dphiB=%15.8e, kfacB=%15.8e\n",
 190                     iparams->dihres.phiA, iparams->dihres.dphiA, iparams->dihres.kfacA,
 191                     iparams->dihres.phiB, iparams->dihres.dphiB, iparams->dihres.kfacB);
 192             break;
 193         case F_POSRES:
 194             fprintf(fp, "pos0A=(%15.8e,%15.8e,%15.8e), fcA=(%15.8e,%15.8e,%15.8e), pos0B=(%15.8e,%15.8e,%15.8e), fcB=(%15.8e,%15.8e,%15.8e)\n",
 195                     iparams->posres.pos0A[XX], iparams->posres.pos0A[YY],
 196                     iparams->posres.pos0A[ZZ], iparams->posres.fcA[XX],
 197                     iparams->posres.fcA[YY], iparams->posres.fcA[ZZ],
 198                     iparams->posres.pos0B[XX], iparams->posres.pos0B[YY],
 199                     iparams->posres.pos0B[ZZ], iparams->posres.fcB[XX],
 200                     iparams->posres.fcB[YY], iparams->posres.fcB[ZZ]);
 201             break;
 202         case F_FBPOSRES:
 203             fprintf(fp, "pos0=(%15.8e,%15.8e,%15.8e), geometry=%d, r=%15.8e, k=%15.8e\n",
 204                     iparams->fbposres.pos0[XX], iparams->fbposres.pos0[YY],
 205                     iparams->fbposres.pos0[ZZ], iparams->fbposres.geom,
 206                     iparams->fbposres.r,        iparams->fbposres.k);
 207             break;
 208         case F_RBDIHS:
 209             for (int i = 0; i < NR_RBDIHS; i++)
 210             {
 211                 fprintf(fp, "%srbcA[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams->rbdihs.rbcA[i]);
 212             }
 213             fprintf(fp, "\n");
 214             for (int i = 0; i < NR_RBDIHS; i++)
 215             {
 216                 fprintf(fp, "%srbcB[%d]=%15.8e", i == 0 ? "" : ", ", i, iparams->rbdihs.rbcB[i]);
 217             }
 218             fprintf(fp, "\n");
 219             break;
 220         case F_FOURDIHS:
 221         {
 222             /* Use the OPLS -> Ryckaert-Bellemans formula backwards to get
 223              * the OPLS potential constants back.
 224              */
 225             const real *rbcA = iparams->rbdihs.rbcA;
 226             const real *rbcB = iparams->rbdihs.rbcB;
 227             real        VA[4], VB[4];
 228
 229             VA[3] = -0.25*rbcA[4];
 230             VA[2] = -0.5*rbcA[3];
 231             VA[1] = 4.0*VA[3]-rbcA[2];
 232             VA[0] = 3.0*VA[2]-2.0*rbcA[1];
 233
 234             VB[3] = -0.25*rbcB[4];
 235             VB[2] = -0.5*rbcB[3];
 236             VB[1] = 4.0*VB[3]-rbcB[2];
 237             VB[0] = 3.0*VB[2]-2.0*rbcB[1];
 238
 239             for (int i = 0; i < NR_FOURDIHS; i++)
 240             {
 241                 fprintf(fp, "%sFourA[%d]=%15.8e", i == 0 ? "" : ", ", i, VA[i]);
 242             }
 243             fprintf(fp, "\n");
 244             for (int i = 0; i < NR_FOURDIHS; i++)
 245             {
 246                 fprintf(fp, "%sFourB[%d]=%15.8e", i == 0 ? "" : ", ", i, VB[i]);
 247             }
 248             fprintf(fp, "\n");
 249             break;
 250         }
 251
 252         case F_CONSTR:
 253         case F_CONSTRNC:
 254             fprintf(fp, "dA=%15.8e, dB=%15.8e\n", iparams->constr.dA, iparams->constr.dB);
 255             break;
 256         case F_SETTLE:
 257             fprintf(fp, "doh=%15.8e, dhh=%15.8e\n", iparams->settle.doh,
 258                     iparams->settle.dhh);
 259             break;
 260         case F_VSITE2:
 261             fprintf(fp, "a=%15.8e\n", iparams->vsite.a);
 262             break;
 263         case F_VSITE3:
 264         case F_VSITE3FD:
 265         case F_VSITE3FAD:
 266             fprintf(fp, "a=%15.8e, b=%15.8e\n", iparams->vsite.a, iparams->vsite.b);
 267             break;
 268         case F_VSITE3OUT:
 269         case F_VSITE4FD:
 270         case F_VSITE4FDN:
 271             fprintf(fp, "a=%15.8e, b=%15.8e, c=%15.8e\n",
 272                     iparams->vsite.a, iparams->vsite.b, iparams->vsite.c);
 273             break;
 274         case F_VSITEN:
 275             fprintf(fp, "n=%2d, a=%15.8e\n", iparams->vsiten.n, iparams->vsiten.a);
 276             break;
 277         case F_GB12_NOLONGERUSED:
 278         case F_GB13_NOLONGERUSED:
 279         case F_GB14_NOLONGERUSED:
 280             // These could only be generated by grompp, not written in
 281             // a .top file. Now that implicit solvent is not
 282             // supported, they can't be generated, and the values are
 283             // ignored if read from an old .tpr file. So there is
 284             // nothing to print.
 285             break;
 286         case F_CMAP:
 287             fprintf(fp, "cmapA=%1d, cmapB=%1d\n", iparams->cmap.cmapA, iparams->cmap.cmapB);
 288             break;
 289         case  F_RESTRANGLES:
 290             pr_harm(fp, iparams, "ktheta", "costheta0");
 291             break;
 292         case  F_RESTRDIHS:
 293             fprintf(fp, "phiA=%15.8e, cpA=%15.8e",
 294                     iparams->pdihs.phiA, iparams->pdihs.cpA);
 295             break;
 296         case  F_CBTDIHS:
 297             fprintf(fp, "kphi=%15.8e", iparams->cbtdihs.cbtcA[0]);
 298             for (int i = 1; i < NR_CBTDIHS; i++)
 299             {
 300                 fprintf(fp, ", cbtcA[%d]=%15.8e", i-1, iparams->cbtdihs.cbtcA[i]);
 301             }
 302             fprintf(fp, "\n");
 303             break;
 304         default:
 305             gmx_fatal(FARGS, "unknown function type %d (%s) in %s line %d",
 306                       ftype, interaction_function[ftype].name, __FILE__, __LINE__);
 307     }
 308 }
 309
 310 template <typename T>
 311 static void
 312 printIlist(FILE *fp, int indent, const char *title,
 313            const t_functype *functype, const T &ilist,
 314            gmx_bool bShowNumbers,
 315            gmx_bool bShowParameters, const t_iparams *iparams)
 316 {
 317     int      i, j, k, type, ftype;
 318
 319     indent = pr_title(fp, indent, title);
 320     pr_indent(fp, indent);
 321     fprintf(fp, "nr: %d\n", ilist.size());
 322     if (ilist.size() > 0)
 323     {
 324         pr_indent(fp, indent);
 325         fprintf(fp, "iatoms:\n");
 326         for (i = j = 0; i < ilist.size(); )
 327         {
 328             pr_indent(fp, indent+INDENT);
 329             type  = ilist.iatoms[i];
 330             ftype = functype[type];
 331             if (bShowNumbers)
 332             {
 333                 fprintf(fp, "%d type=%d ", j, type);
 334             }
 335             j++;
 336             printf("(%s)", interaction_function[ftype].name);
 337             for (k = 0; k < interaction_function[ftype].nratoms; k++)
 338             {
 339                 fprintf(fp, " %3d", ilist.iatoms[i + 1 + k]);
 340             }
 341             if (bShowParameters)
 342             {
 343                 fprintf(fp, "  ");
 344                 pr_iparams(fp, ftype,  &iparams[type]);
 345             }
 346             fprintf(fp, "\n");
 347             i += 1+interaction_function[ftype].nratoms;
 348         }
 349     }
 350 }
 351
 352 void pr_ilist(FILE *fp, int indent, const char *title,
 353               const t_functype *functype, const InteractionList &ilist,
 354               gmx_bool bShowNumbers,
 355               gmx_bool bShowParameters, const t_iparams *iparams)
 356 {
 357     printIlist(fp, indent, title, functype, ilist,
 358                bShowNumbers, bShowParameters, iparams);
 359 }
 360
 361 static void pr_cmap(FILE *fp, int indent, const char *title,
 362                     const gmx_cmap_t *cmap_grid, gmx_bool bShowNumbers)
 363 {
 364     int  i, j, nelem;
 365     real dx, idx;
 366
 367     dx    = 360.0 / cmap_grid->grid_spacing;
 368     nelem = cmap_grid->grid_spacing*cmap_grid->grid_spacing;
 369
 370     if (available(fp, cmap_grid, indent, title))
 371     {
 372         fprintf(fp, "%s\n", title);
 373
 374         for (i = 0; i < static_cast<int>(cmap_grid->cmapdata.size()); i++)
 375         {
 376             idx = -180.0;
 377             fprintf(fp, "%8s %8s %8s %8s\n", "V", "dVdx", "dVdy", "d2dV");
 378
 379             fprintf(fp, "grid[%3d]={\n", bShowNumbers ? i : -1);
 380
 381             for (j = 0; j < nelem; j++)
 382             {
 383                 if ( (j%cmap_grid->grid_spacing) == 0)
 384                 {
 385                     fprintf(fp, "%8.1f\n", idx);
 386                     idx += dx;
 387                 }
 388
 389                 fprintf(fp, "%8.3f ", cmap_grid->cmapdata[i].cmap[j*4]);
 390                 fprintf(fp, "%8.3f ", cmap_grid->cmapdata[i].cmap[j*4+1]);
 391                 fprintf(fp, "%8.3f ", cmap_grid->cmapdata[i].cmap[j*4+2]);
 392                 fprintf(fp, "%8.3f\n", cmap_grid->cmapdata[i].cmap[j*4+3]);
 393             }
 394             fprintf(fp, "\n");
 395         }
 396     }
 397
 398 }
 399
 400 void pr_ffparams(FILE *fp, int indent, const char *title,
 401                  const gmx_ffparams_t *ffparams,
 402                  gmx_bool bShowNumbers)
 403 {
 404     int i;
 405
 406     indent = pr_title(fp, indent, title);
 407     pr_indent(fp, indent);
 408     fprintf(fp, "atnr=%d\n", ffparams->atnr);
 409     pr_indent(fp, indent);
 410     fprintf(fp, "ntypes=%d\n", ffparams->numTypes());
 411     for (i = 0; i < ffparams->numTypes(); i++)
 412     {
 413         pr_indent(fp, indent+INDENT);
 414         fprintf(fp, "functype[%d]=%s, ",
 415                 bShowNumbers ? i : -1,
 416                 interaction_function[ffparams->functype[i]].name);
 417         pr_iparams(fp, ffparams->functype[i], &ffparams->iparams[i]);
 418     }
 419     pr_double(fp, indent, "reppow", ffparams->reppow);
 420     pr_real(fp, indent, "fudgeQQ", ffparams->fudgeQQ);
 421     pr_cmap(fp, indent, "cmap", &ffparams->cmap_grid, bShowNumbers);
 422 }
 423
 424 void pr_idef(FILE *fp, int indent, const char *title, const t_idef *idef,
 425              gmx_bool bShowNumbers, gmx_bool bShowParameters)
 426 {
 427     int i, j;
 428
 429     if (available(fp, idef, indent, title))
 430     {
 431         indent = pr_title(fp, indent, title);
 432         pr_indent(fp, indent);
 433         fprintf(fp, "atnr=%d\n", idef->atnr);
 434         pr_indent(fp, indent);
 435         fprintf(fp, "ntypes=%d\n", idef->ntypes);
 436         for (i = 0; i < idef->ntypes; i++)
 437         {
 438             pr_indent(fp, indent+INDENT);
 439             fprintf(fp, "functype[%d]=%s, ",
 440                     bShowNumbers ? i : -1,
 441                     interaction_function[idef->functype[i]].name);
 442             pr_iparams(fp, idef->functype[i], &idef->iparams[i]);
 443         }
 444         pr_real(fp, indent, "fudgeQQ", idef->fudgeQQ);
 445
 446         for (j = 0; (j < F_NRE); j++)
 447         {
 448             printIlist(fp, indent, interaction_function[j].longname,
 449                        idef->functype, idef->il[j], bShowNumbers,
 450                        bShowParameters, idef->iparams);
 451         }
 452     }
 453 }
 454
 455 void init_idef(t_idef *idef)
 456 {
 457     idef->ntypes           = 0;
 458     idef->atnr             = 0;
 459     idef->functype         = nullptr;
 460     idef->iparams          = nullptr;
 461     idef->fudgeQQ          = 0.0;
 462     idef->iparams_posres   = nullptr;
 463     idef->iparams_fbposres = nullptr;
 464     for (int f = 0; f < F_NRE; ++f)
 465     {
 466         idef->il[f].iatoms          = nullptr;
 467         idef->il[f].nalloc          = 0;
 468         idef->il[f].nr              = 0;
 469         idef->il[f].nr_nonperturbed = 0;
 470     }
 471     idef->cmap_grid               = nullptr;
 472     idef->iparams_posres_nalloc   = 0;
 473     idef->iparams_fbposres_nalloc = 0;
 474     idef->ilsort                  = 0;
 475 }
 476
 477 void done_idef(t_idef *idef)
 478 {
 479     sfree(idef->functype);
 480     sfree(idef->iparams);
 481     sfree(idef->iparams_posres);
 482     sfree(idef->iparams_fbposres);
 483     for (int f = 0; f < F_NRE; ++f)
 484     {
 485         sfree(idef->il[f].iatoms);
 486     }
 487
 488     delete idef->cmap_grid;
 489     init_idef(idef);
 490 }