src/gromacs/gmxpreprocess/x2top.c

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