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