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[], InteractionTypeParameters *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 void set_atom_type(PreprocessingAtomTypes     *atypes,
 182                           t_symtab                   *tab,
 183                           t_atoms                    *atoms,
 184                           InteractionTypeParameters  *bonds,
 185                           int                        *nbonds,
 186                           int                         nnm,
 187                           t_nm2type                   nm2t[])
 188 {
 189     int            nresolved;
 190
 191     snew(atoms->atomtype, atoms->nr);
 192     nresolved = nm2type(nnm, nm2t, tab, atoms, atypes, nbonds, bonds);
 193     if (nresolved != atoms->nr)
 194     {
 195         gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms",
 196                   nresolved, atoms->nr);
 197     }
 198
 199     fprintf(stderr, "There are %zu different atom types in your sample\n",
 200             atypes->size());
 201 }
 202
 203 static void lo_set_force_const(InteractionTypeParameters *plist, real c[], int nrfp, bool bRound,
 204                                bool bDih, bool bParam)
 205 {
 206     double cc;
 207     char   buf[32];
 208
 209     for (int i = 0; (i < plist->nr); i++)
 210     {
 211         if (!bParam)
 212         {
 213             for (int 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 (int 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(gmx::ArrayRef<InteractionTypeParameters> 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(InteractionTypeParameters *ang, InteractionTypeParameters *dih, const rvec x[], bool bPBC,
 268                              matrix box)
 269 {
 270     int    t1, t2, t3;
 271     rvec   r_ij, r_kj, r_kl, m, n;
 272     real   costh;
 273     t_pbc  pbc;
 274
 275     if (bPBC)
 276     {
 277         set_pbc(&pbc, epbcXYZ, box);
 278     }
 279     for (int i = 0; (i < ang->nr); i++)
 280     {
 281         int  ai = ang->param[i].ai();
 282         int  aj = ang->param[i].aj();
 283         int  ak = ang->param[i].ak();
 284         real 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 (int i = 0; (i < dih->nr); i++)
 289     {
 290         int  ai = dih->param[i].ai();
 291         int  aj = dih->param[i].aj();
 292         int  ak = dih->param[i].ak();
 293         int  al = dih->param[i].al();
 294         real 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     for (int i = 0; (i < atoms->nr); i++)
 303     {
 304         fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
 305     }
 306 }
 307
 308 static void print_pl(FILE *fp, gmx::ArrayRef<const InteractionTypeParameters> plist, int ftp, const char *name,
 309                      char ***atomname)
 310 {
 311     if (plist[ftp].nr > 0)
 312     {
 313         fprintf(fp, "\n");
 314         fprintf(fp, "[ %s ]\n", name);
 315         int nral = interaction_function[ftp].nratoms;
 316         int nrfp = interaction_function[ftp].nrfpA;
 317         for (int i = 0; (i < plist[ftp].nr); i++)
 318         {
 319             for (int j = 0; (j < nral); j++)
 320             {
 321                 fprintf(fp, "  %5s", *atomname[plist[ftp].param[i].a[j]]);
 322             }
 323             for (int 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,
 336                       const char                                     *title,
 337                       t_atoms                                        *atoms,
 338                       gmx::ArrayRef<const InteractionTypeParameters>  plist,
 339                       PreprocessingAtomTypes                         *atypes,
 340                       int                                             cgnr[])
 341 {
 342     FILE       *fp;
 343     int         i, tp;
 344     const 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 = atypes->atomNameFromAtomType(tp)) == 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     std::array<InteractionTypeParameters, F_NRE> plist;
 402     t_excls                                     *excls;
 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     PreprocessResidue                            rtp_header_settings;
 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     /* 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     PreprocessingAtomTypes atypes;
 532     set_atom_type(&atypes, &symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
 533
 534     /* Make Angles and Dihedrals */
 535     snew(excls, atoms->nr);
 536     printf("Generating angles and dihedrals from bonds...\n");
 537     init_nnb(&nnb, atoms->nr, 4);
 538     gen_nnb(&nnb, plist);
 539     print_nnb(&nnb, "NNB");
 540     gen_pad(&nnb, atoms, gmx::arrayRefFromArray(&rtp_header_settings, 1), plist, excls, {}, TRUE);
 541     done_nnb(&nnb);
 542
 543     if (!bPairs)
 544     {
 545         plist[F_LJ14].nr = 0;
 546     }
 547     fprintf(stderr,
 548             "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
 549             "          %4d pairs,     %4d bonds and  %4d atoms\n",
 550             plist[F_PDIHS].nr,
 551             bOPLS ? "Ryckaert-Bellemans" : "proper",
 552             plist[F_IDIHS].nr, plist[F_ANGLES].nr,
 553             plist[F_LJ14].nr, plist[F_BONDS].nr, atoms->nr);
 554
 555     calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
 556
 557     set_force_const(plist, kb, kt, kp, bRound, bParam);
 558
 559     cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
 560     printf("Total charge is %g, total mass is %g\n", qtot, mtot);
 561     if (bOPLS)
 562     {
 563         bts[2] = 3;
 564         bts[3] = 1;
 565     }
 566
 567     if (bTOP)
 568     {
 569         fp = ftp2FILE(efTOP, NFILE, fnm, "w");
 570         print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
 571
 572         write_top(fp, nullptr, mymol.name.c_str(), atoms, FALSE, bts, plist, excls, &atypes,
 573                   cgnr, rtp_header_settings.nrexcl);
 574         print_top_mols(fp, mymol.name.c_str(), ffdir, nullptr, {}, gmx::arrayRefFromArray(&mymol, 1));
 575
 576         gmx_ffclose(fp);
 577     }
 578     if (bRTP)
 579     {
 580         print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top",
 581                   atoms, plist, &atypes, cgnr);
 582     }
 583
 584     if (debug)
 585     {
 586         dump_hybridization(debug, atoms, nbonds);
 587     }
 588     close_symtab(&symtab);
 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 }