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