src/programs/gmx/g_x2top.c

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