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