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