src/gromacs/gmxpreprocess/x2top.cpp

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