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