src/gromacs/gmxpreprocess/x2top.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-2008, The GROMACS development team.
   6  * Copyright (c) 2013,2014,2015,2016,2017 by the GROMACS development team.
   7  * Copyright (c) 2018,2019,2020,2021, 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 "x2top.h"
  41
  42 #include <cmath>
  43 #include <cstring>
  44
  45 #include "gromacs/commandline/pargs.h"
  46 #include "gromacs/fileio/confio.h"
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "gromacs/gmxpreprocess/gen_ad.h"
  49 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  50 #include "gromacs/gmxpreprocess/grompp_impl.h"
  51 #include "gromacs/gmxpreprocess/nm2type.h"
  52 #include "gromacs/gmxpreprocess/notset.h"
  53 #include "gromacs/gmxpreprocess/pdb2top.h"
  54 #include "gromacs/gmxpreprocess/toppush.h"
  55 #include "gromacs/gmxpreprocess/toputil.h"
  56 #include "gromacs/listed_forces/bonded.h"
  57 #include "gromacs/math/units.h"
  58 #include "gromacs/math/utilities.h"
  59 #include "gromacs/math/vec.h"
  60 #include "gromacs/math/vecdump.h"
  61 #include "gromacs/mdtypes/md_enums.h"
  62 #include "gromacs/pbcutil/pbc.h"
  63 #include "gromacs/topology/symtab.h"
  64 #include "gromacs/topology/topology.h"
  65 #include "gromacs/utility/arraysize.h"
  66 #include "gromacs/utility/cstringutil.h"
  67 #include "gromacs/utility/fatalerror.h"
  68 #include "gromacs/utility/filestream.h"
  69 #include "gromacs/utility/gmxassert.h"
  70 #include "gromacs/utility/logger.h"
  71 #include "gromacs/utility/loggerbuilder.h"
  72 #include "gromacs/utility/smalloc.h"
  73
  74 #include "hackblock.h"
  75
  76 #define MARGIN_FAC 1.1
  77
  78 static bool is_bond(int nnm, t_nm2type nmt[], char* ai, char* aj, real blen)
  79 {
  80     int i, j;
  81
  82     for (i = 0; (i < nnm); i++)
  83     {
  84         for (j = 0; (j < nmt[i].nbonds); j++)
  85         {
  86             if ((((gmx::equalCaseInsensitive(ai, nmt[i].elem, 1))
  87                   && (gmx::equalCaseInsensitive(aj, nmt[i].bond[j], 1)))
  88                  || ((gmx::equalCaseInsensitive(ai, nmt[i].bond[j], 1))
  89                      && (gmx::equalCaseInsensitive(aj, nmt[i].elem, 1))))
  90                 && (fabs(blen - nmt[i].blen[j]) <= 0.1 * nmt[i].blen[j]))
  91             {
  92                 return TRUE;
  93             }
  94         }
  95     }
  96     return FALSE;
  97 }
  98
  99 static void mk_bonds(int                 nnm,
 100                      t_nm2type           nmt[],
 101                      t_atoms*            atoms,
 102                      const rvec          x[],
 103                      InteractionsOfType* bond,
 104                      int                 nbond[],
 105                      bool                bPBC,
 106                      matrix              box)
 107 {
 108     int   i, j;
 109     t_pbc pbc;
 110     rvec  dx;
 111     real  dx2;
 112
 113     std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
 114     if (bPBC)
 115     {
 116         set_pbc(&pbc, PbcType::Unset, box);
 117     }
 118     for (i = 0; (i < atoms->nr); i++)
 119     {
 120         for (j = i + 1; (j < atoms->nr); j++)
 121         {
 122             if (bPBC)
 123             {
 124                 pbc_dx(&pbc, x[i], x[j], dx);
 125             }
 126             else
 127             {
 128                 rvec_sub(x[i], x[j], dx);
 129             }
 130
 131             dx2 = iprod(dx, dx);
 132             if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j], std::sqrt(dx2)))
 133             {
 134                 forceParam[0]          = std::sqrt(dx2);
 135                 std::vector<int> atoms = { i, j };
 136                 add_param_to_list(bond, InteractionOfType(atoms, forceParam));
 137                 nbond[i]++;
 138                 nbond[j]++;
 139             }
 140         }
 141     }
 142 }
 143
 144 static int* set_cgnr(t_atoms* atoms, bool bUsePDBcharge, real* qtot, real* mtot)
 145 {
 146     int    i, n = 1;
 147     int*   cgnr;
 148     double qt = 0;
 149
 150     *qtot = *mtot = 0;
 151     snew(cgnr, atoms->nr);
 152     for (i = 0; (i < atoms->nr); i++)
 153     {
 154         if (atoms->pdbinfo && bUsePDBcharge)
 155         {
 156             atoms->atom[i].q = atoms->pdbinfo[i].bfac;
 157         }
 158         qt += atoms->atom[i].q;
 159         *qtot += atoms->atom[i].q;
 160         *mtot += atoms->atom[i].m;
 161         cgnr[i] = n;
 162         if (is_int(qt))
 163         {
 164             n++;
 165             qt = 0;
 166         }
 167     }
 168     return cgnr;
 169 }
 170
 171 static void set_atom_type(PreprocessingAtomTypes* atypes,
 172                           t_symtab*               tab,
 173                           t_atoms*                atoms,
 174                           InteractionsOfType*     bonds,
 175                           int*                    nbonds,
 176                           int                     nnm,
 177                           t_nm2type               nm2t[],
 178                           const gmx::MDLogger&    logger)
 179 {
 180     int nresolved;
 181
 182     snew(atoms->atomtype, atoms->nr);
 183     nresolved = nm2type(nnm, nm2t, tab, atoms, atypes, nbonds, bonds);
 184     if (nresolved != atoms->nr)
 185     {
 186         gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms", nresolved, atoms->nr);
 187     }
 188
 189     GMX_LOG(logger.info)
 190             .asParagraph()
 191             .appendTextFormatted("There are %zu different atom types in your sample", atypes->size());
 192 }
 193
 194 static void lo_set_force_const(InteractionsOfType* plist, real c[], int nrfp, bool bRound, bool bDih, bool bParam)
 195 {
 196     double cc;
 197     char   buf[32];
 198
 199     for (auto& param : plist->interactionTypes)
 200     {
 201         if (!bParam)
 202         {
 203             for (int j = 0; j < nrfp; j++)
 204             {
 205                 c[j] = NOTSET;
 206             }
 207         }
 208         else
 209         {
 210             if (bRound)
 211             {
 212                 sprintf(buf, "%.2e", param.c0());
 213                 sscanf(buf, "%lf", &cc);
 214                 c[0] = cc;
 215             }
 216             else
 217             {
 218                 c[0] = param.c0();
 219             }
 220             if (bDih)
 221             {
 222                 c[0] *= c[2];
 223                 c[0] = (static_cast<int>(c[0] + 3600)) % 360;
 224                 if (c[0] > 180)
 225                 {
 226                     c[0] -= 360;
 227                 }
 228                 /* To put the minimum at the current angle rather than the maximum */
 229                 c[0] += 180;
 230             }
 231         }
 232         GMX_ASSERT(nrfp <= MAXFORCEPARAM / 2, "Only 6 parameters may be used for an interaction");
 233         std::array<real, MAXFORCEPARAM> forceParam;
 234         for (int j = 0; (j < nrfp); j++)
 235         {
 236             forceParam[j]        = c[j];
 237             forceParam[nrfp + j] = c[j];
 238         }
 239         param = InteractionOfType(param.atoms(), forceParam);
 240     }
 241 }
 242
 243 static void set_force_const(gmx::ArrayRef<InteractionsOfType> plist, real kb, real kt, real kp, bool bRound, bool bParam)
 244 {
 245     real c[MAXFORCEPARAM];
 246
 247     c[0] = 0;
 248     c[1] = kb;
 249     lo_set_force_const(&plist[F_BONDS], c, 2, bRound, FALSE, bParam);
 250     c[1] = kt;
 251     lo_set_force_const(&plist[F_ANGLES], c, 2, bRound, FALSE, bParam);
 252     c[1] = kp;
 253     c[2] = 3;
 254     lo_set_force_const(&plist[F_PDIHS], c, 3, bRound, TRUE, bParam);
 255 }
 256
 257 static void calc_angles_dihs(InteractionsOfType* ang, InteractionsOfType* dih, const rvec x[], bool bPBC, matrix box)
 258 {
 259     int   t1, t2, t3;
 260     rvec  r_ij, r_kj, r_kl, m, n;
 261     real  costh;
 262     t_pbc pbc;
 263
 264     if (bPBC)
 265     {
 266         set_pbc(&pbc, PbcType::Xyz, box);
 267     }
 268     for (auto& angle : ang->interactionTypes)
 269     {
 270         int  ai = angle.ai();
 271         int  aj = angle.aj();
 272         int  ak = angle.ak();
 273         real th = gmx::c_rad2Deg
 274                   * bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : nullptr, r_ij, r_kj, &costh, &t1, &t2);
 275         angle.setForceParameter(0, th);
 276     }
 277     for (auto dihedral : dih->interactionTypes)
 278     {
 279         int  ai = dihedral.ai();
 280         int  aj = dihedral.aj();
 281         int  ak = dihedral.ak();
 282         int  al = dihedral.al();
 283         real ph =
 284                 gmx::c_rad2Deg
 285                 * dih_angle(
 286                         x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : nullptr, r_ij, r_kj, r_kl, m, n, &t1, &t2, &t3);
 287         dihedral.setForceParameter(0, ph);
 288     }
 289 }
 290
 291 static void dump_hybridization(FILE* fp, t_atoms* atoms, int nbonds[])
 292 {
 293     for (int i = 0; (i < atoms->nr); i++)
 294     {
 295         fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
 296     }
 297 }
 298
 299 static void
 300 print_pl(FILE* fp, gmx::ArrayRef<const InteractionsOfType> plist, int ftp, const char* name, char*** atomname)
 301 {
 302     if (!plist[ftp].interactionTypes.empty())
 303     {
 304         fprintf(fp, "\n");
 305         fprintf(fp, "[ %s ]\n", name);
 306         int nrfp = interaction_function[ftp].nrfpA;
 307         for (const auto& param : plist[ftp].interactionTypes)
 308         {
 309             gmx::ArrayRef<const int>  atoms      = param.atoms();
 310             gmx::ArrayRef<const real> forceParam = param.forceParam();
 311             for (const auto& atom : atoms)
 312             {
 313                 fprintf(fp, "  %5s", *atomname[atom]);
 314             }
 315             for (int j = 0; (j < nrfp); j++)
 316             {
 317                 if (forceParam[j] != NOTSET)
 318                 {
 319                     fprintf(fp, "  %10.3e", forceParam[j]);
 320                 }
 321             }
 322             fprintf(fp, "\n");
 323         }
 324     }
 325 }
 326
 327 static void print_rtp(const char*                             filenm,
 328                       const char*                             title,
 329                       t_atoms*                                atoms,
 330                       gmx::ArrayRef<const InteractionsOfType> plist,
 331                       PreprocessingAtomTypes*                 atypes,
 332                       int                                     cgnr[])
 333 {
 334     FILE* fp;
 335     int   i, tp;
 336
 337     fp = gmx_fio_fopen(filenm, "w");
 338     fprintf(fp, "; %s\n", title);
 339     fprintf(fp, "\n");
 340     fprintf(fp, "[ %s ]\n", *atoms->resinfo[0].name);
 341     fprintf(fp, "\n");
 342     fprintf(fp, "[ atoms ]\n");
 343     for (i = 0; (i < atoms->nr); i++)
 344     {
 345         tp        = atoms->atom[i].type;
 346         auto tpnm = atypes->atomNameFromAtomType(tp);
 347         if (!tpnm.has_value())
 348         {
 349             gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i);
 350         }
 351         fprintf(fp, "%-8s  %12s  %8.4f  %5d\n", *atoms->atomname[i], *tpnm, atoms->atom[i].q, cgnr[i]);
 352     }
 353     print_pl(fp, plist, F_BONDS, "bonds", atoms->atomname);
 354     print_pl(fp, plist, F_ANGLES, "angles", atoms->atomname);
 355     print_pl(fp, plist, F_PDIHS, "dihedrals", atoms->atomname);
 356     print_pl(fp, plist, F_IDIHS, "impropers", atoms->atomname);
 357
 358     gmx_fio_fclose(fp);
 359 }
 360
 361 int gmx_x2top(int argc, char* argv[])
 362 {
 363     const char* desc[] = {
 364         "[THISMODULE] generates a primitive topology from a coordinate file.",
 365         "The program assumes all hydrogens are present when defining",
 366         "the hybridization from the atom name and the number of bonds.",
 367         "The program can also make an [REF].rtp[ref] entry, which you can then add",
 368         "to the [REF].rtp[ref] database.[PAR]",
 369         "When [TT]-param[tt] is set, equilibrium distances and angles",
 370         "and force constants will be printed in the topology for all",
 371         "interactions. The equilibrium distances and angles are taken",
 372         "from the input coordinates, the force constant are set with",
 373         "command line options.",
 374         "The force fields somewhat supported currently are:[PAR]",
 375         "G53a5  GROMOS96 53a5 Forcefield (official distribution)[PAR]",
 376         "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
 377         "The corresponding data files can be found in the library directory",
 378         "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
 379         "information about file formats. By default, the force field selection",
 380         "is interactive, but you can use the [TT]-ff[tt] option to specify",
 381         "one of the short names above on the command line instead. In that",
 382         "case [THISMODULE] just looks for the corresponding file.[PAR]",
 383     };
 384     const char* bugs[] = {
 385         "The atom type selection is primitive. Virtually no chemical knowledge is used",
 386         "Periodic boundary conditions screw up the bonding",
 387         "No improper dihedrals are generated",
 388         "The atoms to atomtype translation table is incomplete ([TT]atomname2type.n2t[tt] file in",
 389         "the data directory). Please extend it and send the results back to the GROMACS crew."
 390     };
 391     FILE*                                 fp;
 392     std::array<InteractionsOfType, F_NRE> plist;
 393     t_excls*                              excls;
 394     t_nm2type*                            nm2t;
 395     t_mols                                mymol;
 396     int                                   nnm;
 397     char                                  forcefield[32], ffdir[STRLEN];
 398     rvec*                                 x; /* coordinates? */
 399     int *                                 nbonds, *cgnr;
 400     int                                   bts[] = { 1, 1, 1, 2 };
 401     matrix                                box;    /* box length matrix */
 402     int                                   natoms; /* number of atoms in one molecule  */
 403     PbcType                               pbcType;
 404     bool                                  bRTP, bTOP, bOPLS;
 405     t_symtab                              symtab;
 406     real                                  qtot, mtot;
 407     char                                  n2t[STRLEN];
 408     gmx_output_env_t*                     oenv;
 409
 410     t_filenm fnm[] = { { efSTX, "-f", "conf", ffREAD },
 411                        { efTOP, "-o", "out", ffOPTWR },
 412                        { efRTP, "-r", "out", ffOPTWR } };
 413 #define NFILE asize(fnm)
 414     real              kb = 4e5, kt = 400, kp = 5;
 415     PreprocessResidue rtp_header_settings;
 416     bool              bRemoveDihedralIfWithImproper = FALSE;
 417     bool              bGenerateHH14Interactions     = TRUE;
 418     bool              bKeepAllGeneratedDihedrals    = FALSE;
 419     int               nrexcl                        = 3;
 420     bool              bParam = TRUE, bRound = TRUE;
 421     bool              bPairs = TRUE, bPBC = TRUE;
 422     bool              bUsePDBcharge = FALSE, bVerbose = FALSE;
 423     const char*       molnm = "ICE";
 424     const char*       ff    = "oplsaa";
 425     t_pargs           pa[]  = {
 426         { "-ff",
 427           FALSE,
 428           etSTR,
 429           { &ff },
 430           "Force field for your simulation. Type \"select\" for interactive selection." },
 431         { "-v", FALSE, etBOOL, { &bVerbose }, "Generate verbose output in the top file." },
 432         { "-nexcl", FALSE, etINT, { &nrexcl }, "Number of exclusions" },
 433         { "-H14",
 434           FALSE,
 435           etBOOL,
 436           { &bGenerateHH14Interactions },
 437           "Use 3rd neighbour interactions for hydrogen atoms" },
 438         { "-alldih",
 439           FALSE,
 440           etBOOL,
 441           { &bKeepAllGeneratedDihedrals },
 442           "Generate all proper dihedrals" },
 443         { "-remdih",
 444           FALSE,
 445           etBOOL,
 446           { &bRemoveDihedralIfWithImproper },
 447           "Remove dihedrals on the same bond as an improper" },
 448         { "-pairs", FALSE, etBOOL, { &bPairs }, "Output 1-4 interactions (pairs) in topology file" },
 449         { "-name", FALSE, etSTR, { &molnm }, "Name of your molecule" },
 450         { "-pbc", FALSE, etBOOL, { &bPBC }, "Use periodic boundary conditions." },
 451         { "-pdbq",
 452           FALSE,
 453           etBOOL,
 454           { &bUsePDBcharge },
 455           "Use the B-factor supplied in a [REF].pdb[ref] file for the atomic charges" },
 456         { "-param", FALSE, etBOOL, { &bParam }, "Print parameters in the output" },
 457         { "-round", FALSE, etBOOL, { &bRound }, "Round off measured values" },
 458         { "-kb", FALSE, etREAL, { &kb }, "Bonded force constant (kJ/mol/nm^2)" },
 459         { "-kt", FALSE, etREAL, { &kt }, "Angle force constant (kJ/mol/rad^2)" },
 460         { "-kp", FALSE, etREAL, { &kp }, "Dihedral angle force constant (kJ/mol/rad^2)" }
 461     };
 462
 463     if (!parse_common_args(
 464                 &argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs, &oenv))
 465     {
 466         return 0;
 467     }
 468     bRTP = opt2bSet("-r", NFILE, fnm);
 469     bTOP = opt2bSet("-o", NFILE, fnm);
 470     /* C89 requirements mean that these struct members cannot be used in
 471      * the declaration of pa. So some temporary variables are needed. */
 472     rtp_header_settings.bRemoveDihedralIfWithImproper = bRemoveDihedralIfWithImproper;
 473     rtp_header_settings.bGenerateHH14Interactions     = bGenerateHH14Interactions;
 474     rtp_header_settings.bKeepAllGeneratedDihedrals    = bKeepAllGeneratedDihedrals;
 475     rtp_header_settings.nrexcl                        = nrexcl;
 476
 477     if (!bRTP && !bTOP)
 478     {
 479         gmx_fatal(FARGS, "Specify at least one output file");
 480     }
 481
 482     gmx::LoggerBuilder builder;
 483     builder.addTargetStream(gmx::MDLogger::LogLevel::Info, &gmx::TextOutputFile::standardOutput());
 484     builder.addTargetStream(gmx::MDLogger::LogLevel::Warning, &gmx::TextOutputFile::standardError());
 485     gmx::LoggerOwner logOwner(builder.build());
 486     gmx::MDLogger    logger(logOwner.logger());
 487
 488
 489     /* Force field selection, interactive or direct */
 490     choose_ff(strcmp(ff, "select") == 0 ? nullptr : ff, forcefield, sizeof(forcefield), ffdir, sizeof(ffdir), logger);
 491
 492     bOPLS = (strcmp(forcefield, "oplsaa") == 0);
 493
 494
 495     mymol.name = gmx_strdup(molnm);
 496     mymol.nr   = 1;
 497
 498     /* Read coordinates */
 499     t_topology* top;
 500     snew(top, 1);
 501     read_tps_conf(opt2fn("-f", NFILE, fnm), top, &pbcType, &x, nullptr, box, FALSE);
 502     t_atoms* atoms = &top->atoms;
 503     natoms         = atoms->nr;
 504     if (atoms->pdbinfo == nullptr)
 505     {
 506         snew(atoms->pdbinfo, natoms);
 507     }
 508
 509     sprintf(n2t, "%s", ffdir);
 510     nm2t = rd_nm2type(n2t, &nnm);
 511     if (nnm == 0)
 512     {
 513         gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)", n2t);
 514     }
 515     else
 516     {
 517         GMX_LOG(logger.info)
 518                 .asParagraph()
 519                 .appendTextFormatted("There are %d name to type translations in file %s", nnm, n2t);
 520     }
 521     if (debug)
 522     {
 523         dump_nm2type(debug, nnm, nm2t);
 524     }
 525     GMX_LOG(logger.info).asParagraph().appendTextFormatted("Generating bonds from distances...");
 526     snew(nbonds, atoms->nr);
 527     mk_bonds(nnm, nm2t, atoms, x, &(plist[F_BONDS]), nbonds, bPBC, box);
 528
 529     open_symtab(&symtab);
 530     PreprocessingAtomTypes atypes;
 531     set_atom_type(&atypes, &symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t, logger);
 532
 533     /* Make Angles and Dihedrals */
 534     snew(excls, atoms->nr);
 535     GMX_LOG(logger.info)
 536             .asParagraph()
 537             .appendTextFormatted("Generating angles and dihedrals from bonds...");
 538     gen_pad(atoms, gmx::arrayRefFromArray(&rtp_header_settings, 1), plist, excls, {}, TRUE, {});
 539
 540     if (!bPairs)
 541     {
 542         plist[F_LJ14].interactionTypes.clear();
 543     }
 544     GMX_LOG(logger.info)
 545             .asParagraph()
 546             .appendTextFormatted(
 547                     "There are %4zu %s dihedrals, %4zu impropers, %4zu angles\n"
 548                     "          %4zu pairs,     %4zu bonds and  %4d atoms\n",
 549                     plist[F_PDIHS].size(),
 550                     bOPLS ? "Ryckaert-Bellemans" : "proper",
 551                     plist[F_IDIHS].size(),
 552                     plist[F_ANGLES].size(),
 553                     plist[F_LJ14].size(),
 554                     plist[F_BONDS].size(),
 555                     atoms->nr);
 556
 557     calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
 558
 559     set_force_const(plist, kb, kt, kp, bRound, bParam);
 560
 561     cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
 562     GMX_LOG(logger.info).asParagraph().appendTextFormatted("Total charge is %g, total mass is %g", qtot, mtot);
 563     if (bOPLS)
 564     {
 565         bts[2] = 3;
 566         bts[3] = 1;
 567     }
 568
 569     if (bTOP)
 570     {
 571         fp = ftp2FILE(efTOP, NFILE, fnm, "w");
 572         print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
 573
 574         write_top(fp,
 575                   nullptr,
 576                   mymol.name.c_str(),
 577                   atoms,
 578                   FALSE,
 579                   bts,
 580                   plist,
 581                   excls,
 582                   &atypes,
 583                   cgnr,
 584                   rtp_header_settings.nrexcl);
 585         print_top_mols(fp, mymol.name.c_str(), ffdir, nullptr, {}, gmx::arrayRefFromArray(&mymol, 1));
 586
 587         gmx_ffclose(fp);
 588     }
 589     if (bRTP)
 590     {
 591         print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top", atoms, plist, &atypes, cgnr);
 592     }
 593
 594     if (debug)
 595     {
 596         dump_hybridization(debug, atoms, nbonds);
 597     }
 598     close_symtab(&symtab);
 599
 600     GMX_LOG(logger.warning)
 601             .asParagraph()
 602             .appendTextFormatted(
 603                     "Topologies generated by %s can not be trusted at face value. "
 604                     "Please verify atomtypes and charges by comparison to other topologies.",
 605                     output_env_get_program_display_name(oenv));
 606
 607     return 0;
 608 }