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,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 "x2top.h"
  40
  41 #include <cmath>
  42 #include <cstring>
  43
  44 #include "gromacs/commandline/pargs.h"
  45 #include "gromacs/fileio/confio.h"
  46 #include "gromacs/fileio/gmxfio.h"
  47 #include "gromacs/gmxpreprocess/gen_ad.h"
  48 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
  49 #include "gromacs/gmxpreprocess/grompp_impl.h"
  50 #include "gromacs/gmxpreprocess/nm2type.h"
  51 #include "gromacs/gmxpreprocess/notset.h"
  52 #include "gromacs/gmxpreprocess/pdb2top.h"
  53 #include "gromacs/gmxpreprocess/toppush.h"
  54 #include "gromacs/gmxpreprocess/toputil.h"
  55 #include "gromacs/listed_forces/bonded.h"
  56 #include "gromacs/math/units.h"
  57 #include "gromacs/math/utilities.h"
  58 #include "gromacs/math/vec.h"
  59 #include "gromacs/math/vecdump.h"
  60 #include "gromacs/mdtypes/md_enums.h"
  61 #include "gromacs/pbcutil/pbc.h"
  62 #include "gromacs/topology/symtab.h"
  63 #include "gromacs/topology/topology.h"
  64 #include "gromacs/utility/arraysize.h"
  65 #include "gromacs/utility/cstringutil.h"
  66 #include "gromacs/utility/fatalerror.h"
  67 #include "gromacs/utility/gmxassert.h"
  68 #include "gromacs/utility/smalloc.h"
  69
  70 #include "hackblock.h"
  71
  72 #define MARGIN_FAC 1.1
  73
  74 static bool is_bond(int nnm, t_nm2type nmt[], char* ai, char* aj, real blen)
  75 {
  76     int i, j;
  77
  78     for (i = 0; (i < nnm); i++)
  79     {
  80         for (j = 0; (j < nmt[i].nbonds); j++)
  81         {
  82             if ((((gmx::equalCaseInsensitive(ai, nmt[i].elem, 1))
  83                   && (gmx::equalCaseInsensitive(aj, nmt[i].bond[j], 1)))
  84                  || ((gmx::equalCaseInsensitive(ai, nmt[i].bond[j], 1))
  85                      && (gmx::equalCaseInsensitive(aj, nmt[i].elem, 1))))
  86                 && (fabs(blen - nmt[i].blen[j]) <= 0.1 * nmt[i].blen[j]))
  87             {
  88                 return TRUE;
  89             }
  90         }
  91     }
  92     return FALSE;
  93 }
  94
  95 static void mk_bonds(int                 nnm,
  96                      t_nm2type           nmt[],
  97                      t_atoms*            atoms,
  98                      const rvec          x[],
  99                      InteractionsOfType* bond,
 100                      int                 nbond[],
 101                      bool                bPBC,
 102                      matrix              box)
 103 {
 104     int   i, j;
 105     t_pbc pbc;
 106     rvec  dx;
 107     real  dx2;
 108
 109     std::array<real, MAXFORCEPARAM> forceParam = { 0.0 };
 110     if (bPBC)
 111     {
 112         set_pbc(&pbc, -1, box);
 113     }
 114     for (i = 0; (i < atoms->nr); i++)
 115     {
 116         if ((i % 10) == 0)
 117         {
 118             fprintf(stderr, "\ratom %d", i);
 119             fflush(stderr);
 120         }
 121         for (j = i + 1; (j < atoms->nr); j++)
 122         {
 123             if (bPBC)
 124             {
 125                 pbc_dx(&pbc, x[i], x[j], dx);
 126             }
 127             else
 128             {
 129                 rvec_sub(x[i], x[j], dx);
 130             }
 131
 132             dx2 = iprod(dx, dx);
 133             if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j], std::sqrt(dx2)))
 134             {
 135                 forceParam[0]          = std::sqrt(dx2);
 136                 std::vector<int> atoms = { i, j };
 137                 add_param_to_list(bond, InteractionOfType(atoms, forceParam));
 138                 nbond[i]++;
 139                 nbond[j]++;
 140             }
 141         }
 142     }
 143     fprintf(stderr, "\ratom %d\n", i);
 144     fflush(stderr);
 145 }
 146
 147 static int* set_cgnr(t_atoms* atoms, bool bUsePDBcharge, real* qtot, real* mtot)
 148 {
 149     int    i, n = 1;
 150     int*   cgnr;
 151     double qt = 0;
 152
 153     *qtot = *mtot = 0;
 154     snew(cgnr, atoms->nr);
 155     for (i = 0; (i < atoms->nr); i++)
 156     {
 157         if (atoms->pdbinfo && bUsePDBcharge)
 158         {
 159             atoms->atom[i].q = atoms->pdbinfo[i].bfac;
 160         }
 161         qt += atoms->atom[i].q;
 162         *qtot += atoms->atom[i].q;
 163         *mtot += atoms->atom[i].m;
 164         cgnr[i] = n;
 165         if (is_int(qt))
 166         {
 167             n++;
 168             qt = 0;
 169         }
 170     }
 171     return cgnr;
 172 }
 173
 174 static void set_atom_type(PreprocessingAtomTypes* atypes,
 175                           t_symtab*               tab,
 176                           t_atoms*                atoms,
 177                           InteractionsOfType*     bonds,
 178                           int*                    nbonds,
 179                           int                     nnm,
 180                           t_nm2type               nm2t[])
 181 {
 182     int nresolved;
 183
 184     snew(atoms->atomtype, atoms->nr);
 185     nresolved = nm2type(nnm, nm2t, tab, atoms, atypes, nbonds, bonds);
 186     if (nresolved != atoms->nr)
 187     {
 188         gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms", nresolved, atoms->nr);
 189     }
 190
 191     fprintf(stderr, "There are %zu different atom types in your sample\n", 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, epbcXYZ, 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 = 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 = RAD2DEG
 284                   * dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : nullptr, r_ij, r_kj, r_kl,
 285                               m, n, &t1, &t2, &t3);
 286         dihedral.setForceParameter(0, ph);
 287     }
 288 }
 289
 290 static void dump_hybridization(FILE* fp, t_atoms* atoms, int nbonds[])
 291 {
 292     for (int i = 0; (i < atoms->nr); i++)
 293     {
 294         fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
 295     }
 296 }
 297
 298 static void
 299 print_pl(FILE* fp, gmx::ArrayRef<const InteractionsOfType> plist, int ftp, const char* name, char*** atomname)
 300 {
 301     if (!plist[ftp].interactionTypes.empty())
 302     {
 303         fprintf(fp, "\n");
 304         fprintf(fp, "[ %s ]\n", name);
 305         int nrfp = interaction_function[ftp].nrfpA;
 306         for (const auto& param : plist[ftp].interactionTypes)
 307         {
 308             gmx::ArrayRef<const int>  atoms      = param.atoms();
 309             gmx::ArrayRef<const real> forceParam = param.forceParam();
 310             for (const auto& atom : atoms)
 311             {
 312                 fprintf(fp, "  %5s", *atomname[atom]);
 313             }
 314             for (int j = 0; (j < nrfp); j++)
 315             {
 316                 if (forceParam[j] != NOTSET)
 317                 {
 318                     fprintf(fp, "  %10.3e", forceParam[j]);
 319                 }
 320             }
 321             fprintf(fp, "\n");
 322         }
 323     }
 324 }
 325
 326 static void print_rtp(const char*                             filenm,
 327                       const char*                             title,
 328                       t_atoms*                                atoms,
 329                       gmx::ArrayRef<const InteractionsOfType> plist,
 330                       PreprocessingAtomTypes*                 atypes,
 331                       int                                     cgnr[])
 332 {
 333     FILE*       fp;
 334     int         i, tp;
 335     const char* tpnm;
 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         if ((tpnm = atypes->atomNameFromAtomType(tp)) == nullptr)
 347         {
 348             gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i);
 349         }
 350         fprintf(fp, "%-8s  %12s  %8.4f  %5d\n", *atoms->atomname[i], tpnm, atoms->atom[i].q, cgnr[i]);
 351     }
 352     print_pl(fp, plist, F_BONDS, "bonds", atoms->atomname);
 353     print_pl(fp, plist, F_ANGLES, "angles", atoms->atomname);
 354     print_pl(fp, plist, F_PDIHS, "dihedrals", atoms->atomname);
 355     print_pl(fp, plist, F_IDIHS, "impropers", atoms->atomname);
 356
 357     gmx_fio_fclose(fp);
 358 }
 359
 360 int gmx_x2top(int argc, char* argv[])
 361 {
 362     const char* desc[] = {
 363         "[THISMODULE] generates a primitive topology from a coordinate file.",
 364         "The program assumes all hydrogens are present when defining",
 365         "the hybridization from the atom name and the number of bonds.",
 366         "The program can also make an [REF].rtp[ref] entry, which you can then add",
 367         "to the [REF].rtp[ref] database.[PAR]",
 368         "When [TT]-param[tt] is set, equilibrium distances and angles",
 369         "and force constants will be printed in the topology for all",
 370         "interactions. The equilibrium distances and angles are taken",
 371         "from the input coordinates, the force constant are set with",
 372         "command line options.",
 373         "The force fields somewhat supported currently are:[PAR]",
 374         "G53a5  GROMOS96 53a5 Forcefield (official distribution)[PAR]",
 375         "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
 376         "The corresponding data files can be found in the library directory",
 377         "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
 378         "information about file formats. By default, the force field selection",
 379         "is interactive, but you can use the [TT]-ff[tt] option to specify",
 380         "one of the short names above on the command line instead. In that",
 381         "case [THISMODULE] just looks for the corresponding file.[PAR]",
 382     };
 383     const char* bugs[] = {
 384         "The atom type selection is primitive. Virtually no chemical knowledge is used",
 385         "Periodic boundary conditions screw up the bonding", "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     int                                   epbc;
 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(&argc, argv, 0, NFILE, fnm, asize(pa), pa, asize(desc), desc,
 462                            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     /* Force field selection, interactive or direct */
 481     choose_ff(strcmp(ff, "select") == 0 ? nullptr : ff, forcefield, sizeof(forcefield), ffdir,
 482               sizeof(ffdir));
 483
 484     bOPLS = (strcmp(forcefield, "oplsaa") == 0);
 485
 486
 487     mymol.name = gmx_strdup(molnm);
 488     mymol.nr   = 1;
 489
 490     /* Read coordinates */
 491     t_topology* top;
 492     snew(top, 1);
 493     read_tps_conf(opt2fn("-f", NFILE, fnm), top, &epbc, &x, nullptr, box, FALSE);
 494     t_atoms* atoms = &top->atoms;
 495     natoms         = atoms->nr;
 496     if (atoms->pdbinfo == nullptr)
 497     {
 498         snew(atoms->pdbinfo, natoms);
 499     }
 500
 501     sprintf(n2t, "%s", ffdir);
 502     nm2t = rd_nm2type(n2t, &nnm);
 503     if (nnm == 0)
 504     {
 505         gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)", n2t);
 506     }
 507     else
 508     {
 509         printf("There are %d name to type translations in file %s\n", nnm, n2t);
 510     }
 511     if (debug)
 512     {
 513         dump_nm2type(debug, nnm, nm2t);
 514     }
 515     printf("Generating bonds from distances...\n");
 516     snew(nbonds, atoms->nr);
 517     mk_bonds(nnm, nm2t, atoms, x, &(plist[F_BONDS]), nbonds, bPBC, box);
 518
 519     open_symtab(&symtab);
 520     PreprocessingAtomTypes atypes;
 521     set_atom_type(&atypes, &symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
 522
 523     /* Make Angles and Dihedrals */
 524     snew(excls, atoms->nr);
 525     printf("Generating angles and dihedrals from bonds...\n");
 526     gen_pad(atoms, gmx::arrayRefFromArray(&rtp_header_settings, 1), plist, excls, {}, TRUE);
 527
 528     if (!bPairs)
 529     {
 530         plist[F_LJ14].interactionTypes.clear();
 531     }
 532     fprintf(stderr,
 533             "There are %4zu %s dihedrals, %4zu impropers, %4zu angles\n"
 534             "          %4zu pairs,     %4zu bonds and  %4d atoms\n",
 535             plist[F_PDIHS].size(), bOPLS ? "Ryckaert-Bellemans" : "proper", plist[F_IDIHS].size(),
 536             plist[F_ANGLES].size(), plist[F_LJ14].size(), plist[F_BONDS].size(), atoms->nr);
 537
 538     calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
 539
 540     set_force_const(plist, kb, kt, kp, bRound, bParam);
 541
 542     cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
 543     printf("Total charge is %g, total mass is %g\n", qtot, mtot);
 544     if (bOPLS)
 545     {
 546         bts[2] = 3;
 547         bts[3] = 1;
 548     }
 549
 550     if (bTOP)
 551     {
 552         fp = ftp2FILE(efTOP, NFILE, fnm, "w");
 553         print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
 554
 555         write_top(fp, nullptr, mymol.name.c_str(), atoms, FALSE, bts, plist, excls, &atypes, cgnr,
 556                   rtp_header_settings.nrexcl);
 557         print_top_mols(fp, mymol.name.c_str(), ffdir, nullptr, {}, gmx::arrayRefFromArray(&mymol, 1));
 558
 559         gmx_ffclose(fp);
 560     }
 561     if (bRTP)
 562     {
 563         print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top", atoms, plist, &atypes, cgnr);
 564     }
 565
 566     if (debug)
 567     {
 568         dump_hybridization(debug, atoms, nbonds);
 569     }
 570     close_symtab(&symtab);
 571
 572     printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",
 573            output_env_get_program_display_name(oenv));
 574     printf("         Please verify atomtypes and charges by comparison to other\n");
 575     printf("         topologies.\n");
 576
 577     return 0;
 578 }