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