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 "copyrite.h"
  44 #include "gromacs/math/utilities.h"
  45 #include "macros.h"
  46 #include "bondf.h"
  47 #include "gromacs/fileio/gmxfio.h"
  48 #include "string2.h"
  49 #include "smalloc.h"
  50 #include "sysstuff.h"
  51 #include "gromacs/fileio/confio.h"
  52 #include "physics.h"
  53 #include "gromacs/commandline/pargs.h"
  54 #include "vec.h"
  55 #include "random.h"
  56 #include "gromacs/math/3dview.h"
  57 #include "txtdump.h"
  58 #include "readinp.h"
  59 #include "names.h"
  60 #include "toppush.h"
  61 #include "pdb2top.h"
  62 #include "gen_ad.h"
  63 #include "gpp_nextnb.h"
  64 #include "vec.h"
  65 #include "atomprop.h"
  66 #include "hackblock.h"
  67
  68 #include "nm2type.h"
  69
  70 char atp[7] = "HCNOSX";
  71 #define NATP (asize(atp)-1)
  72
  73 real blen[NATP][NATP] = {
  74     {  0.00,  0.108, 0.105, 0.10, 0.10, 0.10 },
  75     {  0.108, 0.15,  0.14,  0.14, 0.16, 0.14 },
  76     {  0.105, 0.14,  0.14,  0.14, 0.16, 0.14 },
  77     {  0.10,  0.14,  0.14,  0.14, 0.17, 0.14 },
  78     {  0.10,  0.16,  0.16,  0.17, 0.20, 0.17 },
  79     {  0.10,  0.14,  0.14,  0.14, 0.17, 0.17 }
  80 };
  81
  82 #define MARGIN_FAC 1.1
  83
  84 static gmx_bool is_bond(int nnm, t_nm2type nmt[], char *ai, char *aj, real blen)
  85 {
  86     int i, j;
  87
  88     for (i = 0; (i < nnm); i++)
  89     {
  90         for (j = 0; (j < nmt[i].nbonds); j++)
  91         {
  92             if ((((gmx_strncasecmp(ai, nmt[i].elem, 1) == 0) &&
  93                   (gmx_strncasecmp(aj, nmt[i].bond[j], 1) == 0)) ||
  94                  ((gmx_strncasecmp(ai, nmt[i].bond[j], 1) == 0) &&
  95                   (gmx_strncasecmp(aj, nmt[i].elem, 1) == 0))) &&
  96                 (fabs(blen-nmt[i].blen[j]) <= 0.1*nmt[i].blen[j]))
  97             {
  98                 return TRUE;
  99             }
 100         }
 101     }
 102     return FALSE;
 103 }
 104
 105 static int get_atype(char *nm)
 106 {
 107     int i, aai = NATP-1;
 108
 109     for (i = 0; (i < NATP-1); i++)
 110     {
 111         if (nm[0] == atp[i])
 112         {
 113             aai = i;
 114             break;
 115         }
 116     }
 117     return aai;
 118 }
 119
 120 void mk_bonds(int nnm, t_nm2type nmt[],
 121               t_atoms *atoms, rvec x[], t_params *bond, int nbond[],
 122               gmx_bool bPBC, matrix box)
 123 {
 124     t_param b;
 125     int     i, j;
 126     t_pbc   pbc;
 127     rvec    dx;
 128     real    dx2;
 129
 130     for (i = 0; (i < MAXATOMLIST); i++)
 131     {
 132         b.a[i] = -1;
 133     }
 134     for (i = 0; (i < MAXFORCEPARAM); i++)
 135     {
 136         b.c[i] = 0.0;
 137     }
 138
 139     if (bPBC)
 140     {
 141         set_pbc(&pbc, -1, box);
 142     }
 143     for (i = 0; (i < atoms->nr); i++)
 144     {
 145         if ((i % 10) == 0)
 146         {
 147             fprintf(stderr, "\ratom %d", i);
 148         }
 149         for (j = i+1; (j < atoms->nr); j++)
 150         {
 151             if (bPBC)
 152             {
 153                 pbc_dx(&pbc, x[i], x[j], dx);
 154             }
 155             else
 156             {
 157                 rvec_sub(x[i], x[j], dx);
 158             }
 159
 160             dx2 = iprod(dx, dx);
 161             if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
 162                         sqrt(dx2)))
 163             {
 164                 b.AI = i;
 165                 b.AJ = j;
 166                 b.C0 = sqrt(dx2);
 167                 add_param_to_list (bond, &b);
 168                 nbond[i]++;
 169                 nbond[j]++;
 170                 if (debug)
 171                 {
 172                     fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
 173                             *atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
 174                 }
 175             }
 176         }
 177     }
 178     fprintf(stderr, "\ratom %d\n", i);
 179 }
 180
 181 int *set_cgnr(t_atoms *atoms, gmx_bool bUsePDBcharge, real *qtot, real *mtot)
 182 {
 183     int     i, n = 1;
 184     int    *cgnr;
 185     double  qt = 0, mt = 0;
 186
 187     *qtot = *mtot = 0;
 188     snew(cgnr, atoms->nr);
 189     for (i = 0; (i < atoms->nr); i++)
 190     {
 191         if (atoms->pdbinfo && bUsePDBcharge)
 192         {
 193             atoms->atom[i].q = atoms->pdbinfo[i].bfac;
 194         }
 195         qt     += atoms->atom[i].q;
 196         *qtot  += atoms->atom[i].q;
 197         *mtot  += atoms->atom[i].m;
 198         cgnr[i] = n;
 199         if (is_int(qt))
 200         {
 201             n++;
 202             qt = 0;
 203         }
 204     }
 205     return cgnr;
 206 }
 207
 208 gpp_atomtype_t set_atom_type(t_symtab *tab, t_atoms *atoms, t_params *bonds,
 209                              int *nbonds, int nnm, t_nm2type nm2t[])
 210 {
 211     gpp_atomtype_t atype;
 212     int            nresolved;
 213
 214     atype = init_atomtype();
 215     snew(atoms->atomtype, atoms->nr);
 216     nresolved = nm2type(nnm, nm2t, tab, atoms, atype, nbonds, bonds);
 217     if (nresolved != atoms->nr)
 218     {
 219         gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms",
 220                   nresolved, atoms->nr);
 221     }
 222
 223     fprintf(stderr, "There are %d different atom types in your sample\n",
 224             get_atomtype_ntypes(atype));
 225
 226     return atype;
 227 }
 228
 229 void lo_set_force_const(t_params *plist, real c[], int nrfp, gmx_bool bRound,
 230                         gmx_bool bDih, gmx_bool bParam)
 231 {
 232     int    i, j;
 233     double cc;
 234     char   buf[32];
 235
 236     for (i = 0; (i < plist->nr); i++)
 237     {
 238         if (!bParam)
 239         {
 240             for (j = 0; j < nrfp; j++)
 241             {
 242                 c[j] = NOTSET;
 243             }
 244         }
 245         else
 246         {
 247             if (bRound)
 248             {
 249                 sprintf(buf, "%.2e", plist->param[i].c[0]);
 250                 sscanf(buf, "%lf", &cc);
 251                 c[0] = cc;
 252             }
 253             else
 254             {
 255                 c[0] = plist->param[i].c[0];
 256             }
 257             if (bDih)
 258             {
 259                 c[0] *= c[2];
 260                 c[0]  = ((int)(c[0] + 3600)) % 360;
 261                 if (c[0] > 180)
 262                 {
 263                     c[0] -= 360;
 264                 }
 265                 /* To put the minimum at the current angle rather than the maximum */
 266                 c[0] += 180;
 267             }
 268         }
 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         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     free(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 }