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