src/gromacs/gmxpreprocess/x2top.cpp

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