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.
   7  * Copyright (c) 2019,2020, by the GROMACS development team, led by
   8  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   9  * and including many others, as listed in the AUTHORS file in the
  10  * top-level source directory and at http://www.gromacs.org.
  11  *
  12  * GROMACS is free software; you can redistribute it and/or
  13  * modify it under the terms of the GNU Lesser General Public License
  14  * as published by the Free Software Foundation; either version 2.1
  15  * of the License, or (at your option) any later version.
  16  *
  17  * GROMACS is distributed in the hope that it will be useful,
  18  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  19  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  20  * Lesser General Public License for more details.
  21  *
  22  * You should have received a copy of the GNU Lesser General Public
  23  * License along with GROMACS; if not, see
  24  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  25  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  26  *
  27  * If you want to redistribute modifications to GROMACS, please
  28  * consider that scientific software is very special. Version
  29  * control is crucial - bugs must be traceable. We will be happy to
  30  * consider code for inclusion in the official distribution, but
  31  * derived work must not be called official GROMACS. Details are found
  32  * in the README & COPYING files - if they are missing, get the
  33  * official version at http://www.gromacs.org.
  34  *
  35  * To help us fund GROMACS development, we humbly ask that you cite
  36  * the research papers on the package. Check out http://www.gromacs.org.
  37  */
  38 #include "gmxpre.h"
  39
  40 #include "idef.h"
  41
  42 #include <cstdio>
  43
  44 #include "gromacs/topology/ifunc.h"
  45 #include "gromacs/utility/fatalerror.h"
  46 #include "gromacs/utility/smalloc.h"
  47 #include "gromacs/utility/stringstream.h"
  48 #include "gromacs/utility/textwriter.h"
  49 #include "gromacs/utility/txtdump.h"
  50
  51 static void printHarmonicInteraction(gmx::TextWriter* writer,
  52                                      const t_iparams& iparams,
  53                                      const char*      r,
  54                                      const char*      kr)
  55 {
  56     writer->writeLineFormatted("%sA=%12.5e, %sA=%12.5e, %sB=%12.5e, %sB=%12.5e", r,
  57                                iparams.harmonic.rA, kr, iparams.harmonic.krA, r,
  58                                iparams.harmonic.rB, kr, iparams.harmonic.krB);
  59 }
  60
  61 void pr_iparams(FILE* fp, t_functype ftype, const t_iparams& iparams)
  62 {
  63     gmx::StringOutputStream stream;
  64     {
  65         gmx::TextWriter writer(&stream);
  66         printInteractionParameters(&writer, ftype, iparams);
  67     }
  68     fputs(stream.toString().c_str(), fp);
  69 }
  70
  71 void printInteractionParameters(gmx::TextWriter* writer, t_functype ftype, const t_iparams& iparams)
  72 {
  73     switch (ftype)
  74     {
  75         case F_ANGLES:
  76         case F_G96ANGLES: printHarmonicInteraction(writer, iparams, "th", "ct"); break;
  77         case F_CROSS_BOND_BONDS:
  78             writer->writeLineFormatted("r1e=%15.8e, r2e=%15.8e, krr=%15.8e", iparams.cross_bb.r1e,
  79                                        iparams.cross_bb.r2e, iparams.cross_bb.krr);
  80             break;
  81         case F_CROSS_BOND_ANGLES:
  82             writer->writeLineFormatted("r1e=%15.8e, r1e=%15.8e, r3e=%15.8e, krt=%15.8e",
  83                                        iparams.cross_ba.r1e, iparams.cross_ba.r2e,
  84                                        iparams.cross_ba.r3e, iparams.cross_ba.krt);
  85             break;
  86         case F_LINEAR_ANGLES:
  87             writer->writeLineFormatted("klinA=%15.8e, aA=%15.8e, klinB=%15.8e, aB=%15.8e",
  88                                        iparams.linangle.klinA, iparams.linangle.aA,
  89                                        iparams.linangle.klinB, iparams.linangle.aB);
  90             break;
  91         case F_UREY_BRADLEY:
  92             writer->writeLineFormatted(
  93                     "thetaA=%15.8e, kthetaA=%15.8e, r13A=%15.8e, kUBA=%15.8e, thetaB=%15.8e, "
  94                     "kthetaB=%15.8e, r13B=%15.8e, kUBB=%15.8e",
  95                     iparams.u_b.thetaA, iparams.u_b.kthetaA, iparams.u_b.r13A, iparams.u_b.kUBA,
  96                     iparams.u_b.thetaB, iparams.u_b.kthetaB, iparams.u_b.r13B, iparams.u_b.kUBB);
  97             break;
  98         case F_QUARTIC_ANGLES:
  99             writer->writeStringFormatted("theta=%15.8e", iparams.qangle.theta);
 100             for (int i = 0; i < 5; i++)
 101             {
 102                 writer->writeStringFormatted(", c%c=%15.8e", '0' + i, iparams.qangle.c[i]);
 103             }
 104             writer->ensureLineBreak();
 105             break;
 106         case F_BHAM:
 107             writer->writeLineFormatted("a=%15.8e, b=%15.8e, c=%15.8e", iparams.bham.a,
 108                                        iparams.bham.b, iparams.bham.c);
 109             break;
 110         case F_BONDS:
 111         case F_G96BONDS:
 112         case F_HARMONIC: printHarmonicInteraction(writer, iparams, "b0", "cb"); break;
 113         case F_IDIHS: printHarmonicInteraction(writer, iparams, "xi", "cx"); break;
 114         case F_MORSE:
 115             writer->writeLineFormatted(
 116                     "b0A=%15.8e, cbA=%15.8e, betaA=%15.8e, b0B=%15.8e, cbB=%15.8e, betaB=%15.8e",
 117                     iparams.morse.b0A, iparams.morse.cbA, iparams.morse.betaA, iparams.morse.b0B,
 118                     iparams.morse.cbB, iparams.morse.betaB);
 119             break;
 120         case F_CUBICBONDS:
 121             writer->writeLineFormatted("b0=%15.8e, kb=%15.8e, kcub=%15.8e", iparams.cubic.b0,
 122                                        iparams.cubic.kb, iparams.cubic.kcub);
 123             break;
 124         case F_CONNBONDS: writer->ensureEmptyLine(); break;
 125         case F_FENEBONDS:
 126             writer->writeLineFormatted("bm=%15.8e, kb=%15.8e", iparams.fene.bm, iparams.fene.kb);
 127             break;
 128         case F_RESTRBONDS:
 129             writer->writeLineFormatted(
 130                     "lowA=%15.8e, up1A=%15.8e, up2A=%15.8e, kA=%15.8e, lowB=%15.8e, up1B=%15.8e, "
 131                     "up2B=%15.8e, kB=%15.8e,",
 132                     iparams.restraint.lowA, iparams.restraint.up1A, iparams.restraint.up2A,
 133                     iparams.restraint.kA, iparams.restraint.lowB, iparams.restraint.up1B,
 134                     iparams.restraint.up2B, iparams.restraint.kB);
 135             break;
 136         case F_TABBONDS:
 137         case F_TABBONDSNC:
 138         case F_TABANGLES:
 139         case F_TABDIHS:
 140             writer->writeLineFormatted("tab=%d, kA=%15.8e, kB=%15.8e", iparams.tab.table,
 141                                        iparams.tab.kA, iparams.tab.kB);
 142             break;
 143         case F_POLARIZATION:
 144             writer->writeLineFormatted("alpha=%15.8e", iparams.polarize.alpha);
 145             break;
 146         case F_ANHARM_POL:
 147             writer->writeLineFormatted("alpha=%15.8e drcut=%15.8e khyp=%15.8e",
 148                                        iparams.anharm_polarize.alpha, iparams.anharm_polarize.drcut,
 149                                        iparams.anharm_polarize.khyp);
 150             break;
 151         case F_THOLE_POL:
 152             writer->writeLineFormatted("a=%15.8e, alpha1=%15.8e, alpha2=%15.8e, rfac=%15.8e",
 153                                        iparams.thole.a, iparams.thole.alpha1, iparams.thole.alpha2,
 154                                        iparams.thole.rfac);
 155             break;
 156         case F_WATER_POL:
 157             writer->writeLineFormatted(
 158                     "al_x=%15.8e, al_y=%15.8e, al_z=%15.8e, rOH=%9.6f, rHH=%9.6f, rOD=%9.6f",
 159                     iparams.wpol.al_x, iparams.wpol.al_y, iparams.wpol.al_z, iparams.wpol.rOH,
 160                     iparams.wpol.rHH, iparams.wpol.rOD);
 161             break;
 162         case F_LJ:
 163             writer->writeLineFormatted("c6=%15.8e, c12=%15.8e", iparams.lj.c6, iparams.lj.c12);
 164             break;
 165         case F_LJ14:
 166             writer->writeLineFormatted("c6A=%15.8e, c12A=%15.8e, c6B=%15.8e, c12B=%15.8e",
 167                                        iparams.lj14.c6A, iparams.lj14.c12A, iparams.lj14.c6B,
 168                                        iparams.lj14.c12B);
 169             break;
 170         case F_LJC14_Q:
 171             writer->writeLineFormatted("fqq=%15.8e, qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e",
 172                                        iparams.ljc14.fqq, iparams.ljc14.qi, iparams.ljc14.qj,
 173                                        iparams.ljc14.c6, iparams.ljc14.c12);
 174             break;
 175         case F_LJC_PAIRS_NB:
 176             writer->writeLineFormatted("qi=%15.8e, qj=%15.8e, c6=%15.8e, c12=%15.8e", iparams.ljcnb.qi,
 177                                        iparams.ljcnb.qj, iparams.ljcnb.c6, iparams.ljcnb.c12);
 178             break;
 179         case F_PDIHS:
 180         case F_PIDIHS:
 181         case F_ANGRES:
 182         case F_ANGRESZ:
 183             writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e, phiB=%15.8e, cpB=%15.8e, mult=%d",
 184                                        iparams.pdihs.phiA, iparams.pdihs.cpA, iparams.pdihs.phiB,
 185                                        iparams.pdihs.cpB, iparams.pdihs.mult);
 186             break;
 187         case F_DISRES:
 188             writer->writeLineFormatted(
 189                     "label=%4d, type=%1d, low=%15.8e, up1=%15.8e, up2=%15.8e, fac=%15.8e)",
 190                     iparams.disres.label, iparams.disres.type, iparams.disres.low,
 191                     iparams.disres.up1, iparams.disres.up2, iparams.disres.kfac);
 192             break;
 193         case F_ORIRES:
 194             writer->writeLineFormatted(
 195                     "ex=%4d, label=%d, power=%4d, c=%15.8e, obs=%15.8e, kfac=%15.8e)",
 196                     iparams.orires.ex, iparams.orires.label, iparams.orires.power, iparams.orires.c,
 197                     iparams.orires.obs, iparams.orires.kfac);
 198             break;
 199         case F_DIHRES:
 200             writer->writeLineFormatted(
 201                     "phiA=%15.8e, dphiA=%15.8e, kfacA=%15.8e, phiB=%15.8e, dphiB=%15.8e, "
 202                     "kfacB=%15.8e",
 203                     iparams.dihres.phiA, iparams.dihres.dphiA, iparams.dihres.kfacA,
 204                     iparams.dihres.phiB, iparams.dihres.dphiB, iparams.dihres.kfacB);
 205             break;
 206         case F_POSRES:
 207             writer->writeLineFormatted(
 208                     "pos0A=(%15.8e,%15.8e,%15.8e), fcA=(%15.8e,%15.8e,%15.8e), "
 209                     "pos0B=(%15.8e,%15.8e,%15.8e), fcB=(%15.8e,%15.8e,%15.8e)",
 210                     iparams.posres.pos0A[XX], iparams.posres.pos0A[YY], iparams.posres.pos0A[ZZ],
 211                     iparams.posres.fcA[XX], iparams.posres.fcA[YY], iparams.posres.fcA[ZZ],
 212                     iparams.posres.pos0B[XX], iparams.posres.pos0B[YY], iparams.posres.pos0B[ZZ],
 213                     iparams.posres.fcB[XX], iparams.posres.fcB[YY], iparams.posres.fcB[ZZ]);
 214             break;
 215         case F_FBPOSRES:
 216             writer->writeLineFormatted(
 217                     "pos0=(%15.8e,%15.8e,%15.8e), geometry=%d, r=%15.8e, k=%15.8e",
 218                     iparams.fbposres.pos0[XX], iparams.fbposres.pos0[YY], iparams.fbposres.pos0[ZZ],
 219                     iparams.fbposres.geom, iparams.fbposres.r, iparams.fbposres.k);
 220             break;
 221         case F_RBDIHS:
 222             for (int i = 0; i < NR_RBDIHS; i++)
 223             {
 224                 writer->writeStringFormatted("%srbcA[%d]=%15.8e", i == 0 ? "" : ", ", i,
 225                                              iparams.rbdihs.rbcA[i]);
 226             }
 227             writer->ensureLineBreak();
 228             for (int i = 0; i < NR_RBDIHS; i++)
 229             {
 230                 writer->writeStringFormatted("%srbcB[%d]=%15.8e", i == 0 ? "" : ", ", i,
 231                                              iparams.rbdihs.rbcB[i]);
 232             }
 233             writer->ensureLineBreak();
 234             break;
 235         case F_FOURDIHS:
 236         {
 237             /* Use the OPLS -> Ryckaert-Bellemans formula backwards to get
 238              * the OPLS potential constants back.
 239              */
 240             const real* rbcA = iparams.rbdihs.rbcA;
 241             const real* rbcB = iparams.rbdihs.rbcB;
 242             real        VA[4], VB[4];
 243
 244             VA[3] = -0.25 * rbcA[4];
 245             VA[2] = -0.5 * rbcA[3];
 246             VA[1] = 4.0 * VA[3] - rbcA[2];
 247             VA[0] = 3.0 * VA[2] - 2.0 * rbcA[1];
 248
 249             VB[3] = -0.25 * rbcB[4];
 250             VB[2] = -0.5 * rbcB[3];
 251             VB[1] = 4.0 * VB[3] - rbcB[2];
 252             VB[0] = 3.0 * VB[2] - 2.0 * rbcB[1];
 253
 254             for (int i = 0; i < NR_FOURDIHS; i++)
 255             {
 256                 writer->writeStringFormatted("%sFourA[%d]=%15.8e", i == 0 ? "" : ", ", i, VA[i]);
 257             }
 258             writer->ensureLineBreak();
 259             for (int i = 0; i < NR_FOURDIHS; i++)
 260             {
 261                 writer->writeStringFormatted("%sFourB[%d]=%15.8e", i == 0 ? "" : ", ", i, VB[i]);
 262             }
 263             writer->ensureLineBreak();
 264             break;
 265         }
 266
 267         case F_CONSTR:
 268         case F_CONSTRNC:
 269             writer->writeLineFormatted("dA=%15.8e, dB=%15.8e", iparams.constr.dA, iparams.constr.dB);
 270             break;
 271         case F_SETTLE:
 272             writer->writeLineFormatted("doh=%15.8e, dhh=%15.8e", iparams.settle.doh, iparams.settle.dhh);
 273             break;
 274         case F_VSITE2: writer->writeLineFormatted("a=%15.8e", iparams.vsite.a); break;
 275         case F_VSITE3:
 276         case F_VSITE3FD:
 277         case F_VSITE3FAD:
 278             writer->writeLineFormatted("a=%15.8e, b=%15.8e", iparams.vsite.a, iparams.vsite.b);
 279             break;
 280         case F_VSITE3OUT:
 281         case F_VSITE4FD:
 282         case F_VSITE4FDN:
 283             writer->writeLineFormatted("a=%15.8e, b=%15.8e, c=%15.8e", iparams.vsite.a,
 284                                        iparams.vsite.b, iparams.vsite.c);
 285             break;
 286         case F_VSITEN:
 287             writer->writeLineFormatted("n=%2d, a=%15.8e", iparams.vsiten.n, iparams.vsiten.a);
 288             break;
 289         case F_GB12_NOLONGERUSED:
 290         case F_GB13_NOLONGERUSED:
 291         case F_GB14_NOLONGERUSED:
 292             // These could only be generated by grompp, not written in
 293             // a .top file. Now that implicit solvent is not
 294             // supported, they can't be generated, and the values are
 295             // ignored if read from an old .tpr file. So there is
 296             // nothing to print.
 297             break;
 298         case F_CMAP:
 299             writer->writeLineFormatted("cmapA=%1d, cmapB=%1d", iparams.cmap.cmapA, iparams.cmap.cmapB);
 300             break;
 301         case F_RESTRANGLES: printHarmonicInteraction(writer, iparams, "ktheta", "costheta0"); break;
 302         case F_RESTRDIHS:
 303             writer->writeLineFormatted("phiA=%15.8e, cpA=%15.8e", iparams.pdihs.phiA, iparams.pdihs.cpA);
 304             break;
 305         case F_CBTDIHS:
 306             writer->writeLineFormatted("kphi=%15.8e", iparams.cbtdihs.cbtcA[0]);
 307             for (int i = 1; i < NR_CBTDIHS; i++)
 308             {
 309                 writer->writeStringFormatted(", cbtcA[%d]=%15.8e", i - 1, iparams.cbtdihs.cbtcA[i]);
 310             }
 311             writer->ensureLineBreak();
 312             break;
 313         default:
 314             gmx_fatal(FARGS, "unknown function type %d (%s) in %s line %d", ftype,
 315                       interaction_function[ftype].name, __FILE__, __LINE__);
 316     }
 317 }
 318
 319 template<typename T>
 320 static void printIlist(FILE*             fp,
 321                        int               indent,
 322                        const char*       title,
 323                        const t_functype* functype,
 324                        const T&          ilist,
 325                        gmx_bool          bShowNumbers,
 326                        gmx_bool          bShowParameters,
 327                        const t_iparams*  iparams)
 328 {
 329     int i, j, k, type, ftype;
 330
 331     indent = pr_title(fp, indent, title);
 332     pr_indent(fp, indent);
 333     fprintf(fp, "nr: %d\n", ilist.size());
 334     if (ilist.size() > 0)
 335     {
 336         pr_indent(fp, indent);
 337         fprintf(fp, "iatoms:\n");
 338         for (i = j = 0; i < ilist.size();)
 339         {
 340             pr_indent(fp, indent + INDENT);
 341             type  = ilist.iatoms[i];
 342             ftype = functype[type];
 343             if (bShowNumbers)
 344             {
 345                 fprintf(fp, "%d type=%d ", j, type);
 346             }
 347             j++;
 348             printf("(%s)", interaction_function[ftype].name);
 349             for (k = 0; k < interaction_function[ftype].nratoms; k++)
 350             {
 351                 fprintf(fp, " %3d", ilist.iatoms[i + 1 + k]);
 352             }
 353             if (bShowParameters)
 354             {
 355                 fprintf(fp, "  ");
 356                 pr_iparams(fp, ftype, iparams[type]);
 357             }
 358             fprintf(fp, "\n");
 359             i += 1 + interaction_function[ftype].nratoms;
 360         }
 361     }
 362 }
 363
 364 void pr_ilist(FILE*                  fp,
 365               int                    indent,
 366               const char*            title,
 367               const t_functype*      functype,
 368               const InteractionList& ilist,
 369               gmx_bool               bShowNumbers,
 370               gmx_bool               bShowParameters,
 371               const t_iparams*       iparams)
 372 {
 373     printIlist(fp, indent, title, functype, ilist, bShowNumbers, bShowParameters, iparams);
 374 }
 375
 376 void pr_idef(FILE* fp, int indent, const char* title, const t_idef* idef, gmx_bool bShowNumbers, gmx_bool bShowParameters)
 377 {
 378     int i, j;
 379
 380     if (available(fp, idef, indent, title))
 381     {
 382         indent = pr_title(fp, indent, title);
 383         pr_indent(fp, indent);
 384         fprintf(fp, "atnr=%d\n", idef->atnr);
 385         pr_indent(fp, indent);
 386         fprintf(fp, "ntypes=%d\n", idef->ntypes);
 387         for (i = 0; i < idef->ntypes; i++)
 388         {
 389             pr_indent(fp, indent + INDENT);
 390             fprintf(fp, "functype[%d]=%s, ", bShowNumbers ? i : -1,
 391                     interaction_function[idef->functype[i]].name);
 392             pr_iparams(fp, idef->functype[i], idef->iparams[i]);
 393         }
 394         pr_real(fp, indent, "fudgeQQ", idef->fudgeQQ);
 395
 396         for (j = 0; (j < F_NRE); j++)
 397         {
 398             printIlist(fp, indent, interaction_function[j].longname, idef->functype, idef->il[j],
 399                        bShowNumbers, bShowParameters, idef->iparams);
 400         }
 401     }
 402 }
 403
 404 void init_idef(t_idef* idef)
 405 {
 406     idef->ntypes           = 0;
 407     idef->atnr             = 0;
 408     idef->functype         = nullptr;
 409     idef->iparams          = nullptr;
 410     idef->fudgeQQ          = 0.0;
 411     idef->iparams_posres   = nullptr;
 412     idef->iparams_fbposres = nullptr;
 413     for (int f = 0; f < F_NRE; ++f)
 414     {
 415         idef->il[f].iatoms                   = nullptr;
 416         idef->il[f].nalloc                   = 0;
 417         idef->il[f].nr                       = 0;
 418         idef->numNonperturbedInteractions[f] = 0;
 419     }
 420     idef->cmap_grid               = nullptr;
 421     idef->iparams_posres_nalloc   = 0;
 422     idef->iparams_fbposres_nalloc = 0;
 423     idef->ilsort                  = ilsortUNKNOWN;
 424 }
 425
 426 void done_idef(t_idef* idef)
 427 {
 428     sfree(idef->functype);
 429     sfree(idef->iparams);
 430     sfree(idef->iparams_posres);
 431     sfree(idef->iparams_fbposres);
 432     for (int f = 0; f < F_NRE; ++f)
 433     {
 434         sfree(idef->il[f].iatoms);
 435     }
 436
 437     delete idef->cmap_grid;
 438     init_idef(idef);
 439 }
 440
 441 void copy_ilist(const t_ilist* src, t_ilist* dst)
 442 {
 443     dst->nr     = src->nr;
 444     dst->nalloc = src->nr;
 445
 446     snew(dst->iatoms, dst->nalloc);
 447     for (int i = 0; i < dst->nr; ++i)
 448     {
 449         dst->iatoms[i] = src->iatoms[i];
 450     }
 451 }