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