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