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