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