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, 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/gmxassert.h"
  69 #include "gromacs/utility/smalloc.h"
  70
  71 #include "hackblock.h"
  72
  73 #define MARGIN_FAC 1.1
  74
  75 static bool is_bond(int nnm, t_nm2type nmt[], char* ai, char* aj, real blen)
  76 {
  77     int i, j;
  78
  79     for (i = 0; (i < nnm); i++)
  80     {
  81         for (j = 0; (j < nmt[i].nbonds); j++)
  82         {
  83             if ((((gmx::equalCaseInsensitive(ai, nmt[i].elem, 1))
  84                   && (gmx::equalCaseInsensitive(aj, nmt[i].bond[j], 1)))
  85                  || ((gmx::equalCaseInsensitive(ai, nmt[i].bond[j], 1))
  86                      && (gmx::equalCaseInsensitive(aj, nmt[i].elem, 1))))
  87                 && (fabs(blen - nmt[i].blen[j]) <= 0.1 * nmt[i].blen[j]))
  88             {
  89                 return TRUE;
  90             }
  91         }
  92     }
  93     return FALSE;
  94 }
  95
  96 static void mk_bonds(int                 nnm,
  97                      t_nm2type           nmt[],
  98                      t_atoms*            atoms,
  99                      const rvec          x[],
 100                      InteractionsOfType* bond,
 101                      int                 nbond[],
 102                      bool                bPBC,
 103                      matrix              box)
 104 {
 105     int   i, j;
 106     t_pbc pbc;
 107     rvec  dx;
 108     real  dx2;
 109
 110     std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
 111     if (bPBC)
 112     {
 113         set_pbc(&pbc, PbcType::Unset, box);
 114     }
 115     for (i = 0; (i < atoms->nr); i++)
 116     {
 117         if ((i % 10) == 0)
 118         {
 119             fprintf(stderr, "\ratom %d", i);
 120             fflush(stderr);
 121         }
 122         for (j = i + 1; (j < atoms->nr); j++)
 123         {
 124             if (bPBC)
 125             {
 126                 pbc_dx(&pbc, x[i], x[j], dx);
 127             }
 128             else
 129             {
 130                 rvec_sub(x[i], x[j], dx);
 131             }
 132
 133             dx2 = iprod(dx, dx);
 134             if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j], std::sqrt(dx2)))
 135             {
 136                 forceParam[0]          = std::sqrt(dx2);
 137                 std::vector<int> atoms = { i, j };
 138                 add_param_to_list(bond, InteractionOfType(atoms, forceParam));
 139                 nbond[i]++;
 140                 nbond[j]++;
 141             }
 142         }
 143     }
 144     fprintf(stderr, "\ratom %d\n", i);
 145     fflush(stderr);
 146 }
 147
 148 static int* set_cgnr(t_atoms* atoms, bool bUsePDBcharge, real* qtot, real* mtot)
 149 {
 150     int    i, n = 1;
 151     int*   cgnr;
 152     double qt = 0;
 153
 154     *qtot = *mtot = 0;
 155     snew(cgnr, atoms->nr);
 156     for (i = 0; (i < atoms->nr); i++)
 157     {
 158         if (atoms->pdbinfo && bUsePDBcharge)
 159         {
 160             atoms->atom[i].q = atoms->pdbinfo[i].bfac;
 161         }
 162         qt += atoms->atom[i].q;
 163         *qtot += atoms->atom[i].q;
 164         *mtot += atoms->atom[i].m;
 165         cgnr[i] = n;
 166         if (is_int(qt))
 167         {
 168             n++;
 169             qt = 0;
 170         }
 171     }
 172     return cgnr;
 173 }
 174
 175 static void set_atom_type(PreprocessingAtomTypes* atypes,
 176                           t_symtab*               tab,
 177                           t_atoms*                atoms,
 178                           InteractionsOfType*     bonds,
 179                           int*                    nbonds,
 180                           int                     nnm,
 181                           t_nm2type               nm2t[])
 182 {
 183     int nresolved;
 184
 185     snew(atoms->atomtype, atoms->nr);
 186     nresolved = nm2type(nnm, nm2t, tab, atoms, atypes, nbonds, bonds);
 187     if (nresolved != atoms->nr)
 188     {
 189         gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms", nresolved, atoms->nr);
 190     }
 191
 192     fprintf(stderr, "There are %zu different atom types in your sample\n", atypes->size());
 193 }
 194
 195 static void lo_set_force_const(InteractionsOfType* plist, real c[], int nrfp, bool bRound, bool bDih, bool bParam)
 196 {
 197     double cc;
 198     char   buf[32];
 199
 200     for (auto& param : plist->interactionTypes)
 201     {
 202         if (!bParam)
 203         {
 204             for (int j = 0; j < nrfp; j++)
 205             {
 206                 c[j] = NOTSET;
 207             }
 208         }
 209         else
 210         {
 211             if (bRound)
 212             {
 213                 sprintf(buf, "%.2e", param.c0());
 214                 sscanf(buf, "%lf", &cc);
 215                 c[0] = cc;
 216             }
 217             else
 218             {
 219                 c[0] = param.c0();
 220             }
 221             if (bDih)
 222             {
 223                 c[0] *= c[2];
 224                 c[0] = (static_cast<int>(c[0] + 3600)) % 360;
 225                 if (c[0] > 180)
 226                 {
 227                     c[0] -= 360;
 228                 }
 229                 /* To put the minimum at the current angle rather than the maximum */
 230                 c[0] += 180;
 231             }
 232         }
 233         GMX_ASSERT(nrfp <= MAXFORCEPARAM / 2, "Only 6 parameters may be used for an interaction");
 234         std::array<real, MAXFORCEPARAM> forceParam;
 235         for (int j = 0; (j < nrfp); j++)
 236         {
 237             forceParam[j]        = c[j];
 238             forceParam[nrfp + j] = c[j];
 239         }
 240         param = InteractionOfType(param.atoms(), forceParam);
 241     }
 242 }
 243
 244 static void set_force_const(gmx::ArrayRef<InteractionsOfType> plist, real kb, real kt, real kp, bool bRound, bool bParam)
 245 {
 246     real c[MAXFORCEPARAM];
 247
 248     c[0] = 0;
 249     c[1] = kb;
 250     lo_set_force_const(&plist[F_BONDS], c, 2, bRound, FALSE, bParam);
 251     c[1] = kt;
 252     lo_set_force_const(&plist[F_ANGLES], c, 2, bRound, FALSE, bParam);
 253     c[1] = kp;
 254     c[2] = 3;
 255     lo_set_force_const(&plist[F_PDIHS], c, 3, bRound, TRUE, bParam);
 256 }
 257
 258 static void calc_angles_dihs(InteractionsOfType* ang, InteractionsOfType* dih, const rvec x[], bool bPBC, matrix box)
 259 {
 260     int   t1, t2, t3;
 261     rvec  r_ij, r_kj, r_kl, m, n;
 262     real  costh;
 263     t_pbc pbc;
 264
 265     if (bPBC)
 266     {
 267         set_pbc(&pbc, PbcType::Xyz, box);
 268     }
 269     for (auto& angle : ang->interactionTypes)
 270     {
 271         int  ai = angle.ai();
 272         int  aj = angle.aj();
 273         int  ak = angle.ak();
 274         real th = RAD2DEG
 275                   * bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : nullptr, r_ij, r_kj, &costh, &t1, &t2);
 276         angle.setForceParameter(0, th);
 277     }
 278     for (auto dihedral : dih->interactionTypes)
 279     {
 280         int  ai = dihedral.ai();
 281         int  aj = dihedral.aj();
 282         int  ak = dihedral.ak();
 283         int  al = dihedral.al();
 284         real ph = RAD2DEG
 285                   * dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : nullptr, r_ij, r_kj, r_kl,
 286                               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     const char* tpnm;
 337
 338     fp = gmx_fio_fopen(filenm, "w");
 339     fprintf(fp, "; %s\n", title);
 340     fprintf(fp, "\n");
 341     fprintf(fp, "[ %s ]\n", *atoms->resinfo[0].name);
 342     fprintf(fp, "\n");
 343     fprintf(fp, "[ atoms ]\n");
 344     for (i = 0; (i < atoms->nr); i++)
 345     {
 346         tp = atoms->atom[i].type;
 347         if ((tpnm = atypes->atomNameFromAtomType(tp)) == nullptr)
 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", "No improper dihedrals are generated",
 387         "The atoms to atomtype translation table is incomplete ([TT]atomname2type.n2t[tt] file in "
 388         "the data directory). Please extend it and send the results back to the GROMACS crew."
 389     };
 390     FILE*                                 fp;
 391     std::array<InteractionsOfType, F_NRE> plist;
 392     t_excls*                              excls;
 393     t_nm2type*                            nm2t;
 394     t_mols                                mymol;
 395     int                                   nnm;
 396     char                                  forcefield[32], ffdir[STRLEN];
 397     rvec*                                 x; /* coordinates? */
 398     int *                                 nbonds, *cgnr;
 399     int                                   bts[] = { 1, 1, 1, 2 };
 400     matrix                                box;    /* box length matrix */
 401     int                                   natoms; /* number of atoms in one molecule  */
 402     PbcType                               pbcType;
 403     bool                                  bRTP, bTOP, bOPLS;
 404     t_symtab                              symtab;
 405     real                                  qtot, mtot;
 406     char                                  n2t[STRLEN];
 407     gmx_output_env_t*                     oenv;
 408
 409     t_filenm fnm[] = { { efSTX, "-f", "conf", ffREAD },
 410                        { efTOP, "-o", "out", ffOPTWR },
 411                        { efRTP, "-r", "out", ffOPTWR } };
 412 #define NFILE asize(fnm)
 413     real              kb = 4e5, kt = 400, kp = 5;
 414     PreprocessResidue rtp_header_settings;
 415     bool              bRemoveDihedralIfWithImproper = FALSE;
 416     bool              bGenerateHH14Interactions     = TRUE;
 417     bool              bKeepAllGeneratedDihedrals    = FALSE;
 418     int               nrexcl                        = 3;
 419     bool              bParam = TRUE, bRound = TRUE;
 420     bool              bPairs = TRUE, bPBC = TRUE;
 421     bool              bUsePDBcharge = FALSE, bVerbose = FALSE;
 422     const char*       molnm = "ICE";
 423     const char*       ff    = "oplsaa";
 424     t_pargs           pa[]  = {
 425         { "-ff",
 426           FALSE,
 427           etSTR,
 428           { &ff },
 429           "Force field for your simulation. Type \"select\" for interactive selection." },
 430         { "-v", FALSE, etBOOL, { &bVerbose }, "Generate verbose output in the top file." },
 431         { "-nexcl", FALSE, etINT, { &nrexcl }, "Number of exclusions" },
 432         { "-H14",
 433           FALSE,
 434           etBOOL,
 435           { &bGenerateHH14Interactions },
 436           "Use 3rd neighbour interactions for hydrogen atoms" },
 437         { "-alldih",
 438           FALSE,
 439           etBOOL,
 440           { &bKeepAllGeneratedDihedrals },
 441           "Generate all proper dihedrals" },
 442         { "-remdih",
 443           FALSE,
 444           etBOOL,
 445           { &bRemoveDihedralIfWithImproper },
 446           "Remove dihedrals on the same bond as an improper" },
 447         { "-pairs", FALSE, etBOOL, { &bPairs }, "Output 1-4 interactions (pairs) in topology file" },
 448         { "-name", FALSE, etSTR, { &molnm }, "Name of your molecule" },
 449         { "-pbc", FALSE, etBOOL, { &bPBC }, "Use periodic boundary conditions." },
 450         { "-pdbq",
 451           FALSE,
 452           etBOOL,
 453           { &bUsePDBcharge },
 454           "Use the B-factor supplied in a [REF].pdb[ref] file for the atomic charges" },
 455         { "-param", FALSE, etBOOL, { &bParam }, "Print parameters in the output" },
 456         { "-round", FALSE, etBOOL, { &bRound }, "Round off measured values" },
 457         { "-kb", FALSE, etREAL, { &kb }, "Bonded force constant (kJ/mol/nm^2)" },
 458         { "-kt", FALSE, etREAL, { &kt }, "Angle force constant (kJ/mol/rad^2)" },
 459         { "-kp", FALSE, etREAL, { &kp }, "Dihedral angle force constant (kJ/mol/rad^2)" }
 460     };
 461
 462     if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc,
 463                            asize(bugs), bugs, &oenv))
 464     {
 465         return 0;
 466     }
 467     bRTP = opt2bSet("-r", NFILE, fnm);
 468     bTOP = opt2bSet("-o", NFILE, fnm);
 469     /* C89 requirements mean that these struct members cannot be used in
 470      * the declaration of pa. So some temporary variables are needed. */
 471     rtp_header_settings.bRemoveDihedralIfWithImproper = bRemoveDihedralIfWithImproper;
 472     rtp_header_settings.bGenerateHH14Interactions     = bGenerateHH14Interactions;
 473     rtp_header_settings.bKeepAllGeneratedDihedrals    = bKeepAllGeneratedDihedrals;
 474     rtp_header_settings.nrexcl                        = nrexcl;
 475
 476     if (!bRTP && !bTOP)
 477     {
 478         gmx_fatal(FARGS, "Specify at least one output file");
 479     }
 480
 481     /* Force field selection, interactive or direct */
 482     choose_ff(strcmp(ff, "select") == 0 ? nullptr : ff, forcefield, sizeof(forcefield), ffdir,
 483               sizeof(ffdir));
 484
 485     bOPLS = (strcmp(forcefield, "oplsaa") == 0);
 486
 487
 488     mymol.name = gmx_strdup(molnm);
 489     mymol.nr   = 1;
 490
 491     /* Read coordinates */
 492     t_topology* top;
 493     snew(top, 1);
 494     read_tps_conf(opt2fn("-f", NFILE, fnm), top, &pbcType, &x, nullptr, box, FALSE);
 495     t_atoms* atoms = &top->atoms;
 496     natoms         = atoms->nr;
 497     if (atoms->pdbinfo == nullptr)
 498     {
 499         snew(atoms->pdbinfo, natoms);
 500     }
 501
 502     sprintf(n2t, "%s", ffdir);
 503     nm2t = rd_nm2type(n2t, &nnm);
 504     if (nnm == 0)
 505     {
 506         gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)", n2t);
 507     }
 508     else
 509     {
 510         printf("There are %d name to type translations in file %s\n", nnm, n2t);
 511     }
 512     if (debug)
 513     {
 514         dump_nm2type(debug, nnm, nm2t);
 515     }
 516     printf("Generating bonds from distances...\n");
 517     snew(nbonds, atoms->nr);
 518     mk_bonds(nnm, nm2t, atoms, x, &(plist[F_BONDS]), nbonds, bPBC, box);
 519
 520     open_symtab(&symtab);
 521     PreprocessingAtomTypes atypes;
 522     set_atom_type(&atypes, &symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
 523
 524     /* Make Angles and Dihedrals */
 525     snew(excls, atoms->nr);
 526     printf("Generating angles and dihedrals from bonds...\n");
 527     gen_pad(atoms, gmx::arrayRefFromArray(&rtp_header_settings, 1), plist, excls, {}, TRUE);
 528
 529     if (!bPairs)
 530     {
 531         plist[F_LJ14].interactionTypes.clear();
 532     }
 533     fprintf(stderr,
 534             "There are %4zu %s dihedrals, %4zu impropers, %4zu angles\n"
 535             "          %4zu pairs,     %4zu bonds and  %4d atoms\n",
 536             plist[F_PDIHS].size(), bOPLS ? "Ryckaert-Bellemans" : "proper", plist[F_IDIHS].size(),
 537             plist[F_ANGLES].size(), plist[F_LJ14].size(), plist[F_BONDS].size(), atoms->nr);
 538
 539     calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
 540
 541     set_force_const(plist, kb, kt, kp, bRound, bParam);
 542
 543     cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
 544     printf("Total charge is %g, total mass is %g\n", qtot, mtot);
 545     if (bOPLS)
 546     {
 547         bts[2] = 3;
 548         bts[3] = 1;
 549     }
 550
 551     if (bTOP)
 552     {
 553         fp = ftp2FILE(efTOP, NFILE, fnm, "w");
 554         print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
 555
 556         write_top(fp, nullptr, mymol.name.c_str(), atoms, FALSE, bts, plist, excls, &atypes, cgnr,
 557                   rtp_header_settings.nrexcl);
 558         print_top_mols(fp, mymol.name.c_str(), ffdir, nullptr, {}, gmx::arrayRefFromArray(&mymol, 1));
 559
 560         gmx_ffclose(fp);
 561     }
 562     if (bRTP)
 563     {
 564         print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top", atoms, plist, &atypes, cgnr);
 565     }
 566
 567     if (debug)
 568     {
 569         dump_hybridization(debug, atoms, nbonds);
 570     }
 571     close_symtab(&symtab);
 572
 573     printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",
 574            output_env_get_program_display_name(oenv));
 575     printf("         Please verify atomtypes and charges by comparison to other\n");
 576     printf("         topologies.\n");
 577
 578     return 0;
 579 }