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