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