2 * This file is part of the GROMACS molecular simulation package.
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.
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.
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.
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.
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.
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.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/fileio/confio.h"
45 #include "gromacs/fileio/gmxfio.h"
46 #include "gromacs/gmxpreprocess/gen_ad.h"
47 #include "gromacs/gmxpreprocess/gpp_nextnb.h"
48 #include "gromacs/gmxpreprocess/hackblock.h"
49 #include "gromacs/gmxpreprocess/nm2type.h"
50 #include "gromacs/gmxpreprocess/pdb2top.h"
51 #include "gromacs/gmxpreprocess/toppush.h"
52 #include "gromacs/legacyheaders/copyrite.h"
53 #include "gromacs/legacyheaders/macros.h"
54 #include "gromacs/legacyheaders/names.h"
55 #include "gromacs/legacyheaders/readinp.h"
56 #include "gromacs/legacyheaders/txtdump.h"
57 #include "gromacs/listed-forces/bonded.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"
67 char atp[7] = "HCNOSX";
68 #define NATP (asize(atp)-1)
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 }
79 #define MARGIN_FAC 1.1
81 static gmx_bool is_bond(int nnm, t_nm2type nmt[], char *ai, char *aj, real blen)
85 for (i = 0; (i < nnm); i++)
87 for (j = 0; (j < nmt[i].nbonds); j++)
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]))
102 static int get_atype(char *nm)
106 for (i = 0; (i < NATP-1); i++)
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)
127 for (i = 0; (i < MAXATOMLIST); i++)
131 for (i = 0; (i < MAXFORCEPARAM); i++)
138 set_pbc(&pbc, -1, box);
140 for (i = 0; (i < atoms->nr); i++)
144 fprintf(stderr, "\ratom %d", i);
146 for (j = i+1; (j < atoms->nr); j++)
150 pbc_dx(&pbc, x[i], x[j], dx);
154 rvec_sub(x[i], x[j], dx);
158 if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
164 add_param_to_list (bond, &b);
169 fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
170 *atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
175 fprintf(stderr, "\ratom %d\n", i);
178 int *set_cgnr(t_atoms *atoms, gmx_bool bUsePDBcharge, real *qtot, real *mtot)
182 double qt = 0, mt = 0;
185 snew(cgnr, atoms->nr);
186 for (i = 0; (i < atoms->nr); i++)
188 if (atoms->pdbinfo && bUsePDBcharge)
190 atoms->atom[i].q = atoms->pdbinfo[i].bfac;
192 qt += atoms->atom[i].q;
193 *qtot += atoms->atom[i].q;
194 *mtot += atoms->atom[i].m;
205 gpp_atomtype_t set_atom_type(t_symtab *tab, t_atoms *atoms, t_params *bonds,
206 int *nbonds, int nnm, t_nm2type nm2t[])
208 gpp_atomtype_t atype;
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)
216 gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms",
217 nresolved, atoms->nr);
220 fprintf(stderr, "There are %d different atom types in your sample\n",
221 get_atomtype_ntypes(atype));
226 void lo_set_force_const(t_params *plist, real c[], int nrfp, gmx_bool bRound,
227 gmx_bool bDih, gmx_bool bParam)
233 for (i = 0; (i < plist->nr); i++)
237 for (j = 0; j < nrfp; j++)
246 sprintf(buf, "%.2e", plist->param[i].c[0]);
247 sscanf(buf, "%lf", &cc);
252 c[0] = plist->param[i].c[0];
257 c[0] = ((int)(c[0] + 3600)) % 360;
262 /* To put the minimum at the current angle rather than the maximum */
266 assert(nrfp <= MAXFORCEPARAM/2);
267 for (j = 0; (j < nrfp); j++)
269 plist->param[i].c[j] = c[j];
270 plist->param[i].c[nrfp+j] = c[j];
272 set_p_string(&(plist->param[i]), "");
276 void set_force_const(t_params plist[], real kb, real kt, real kp, gmx_bool bRound,
280 real c[MAXFORCEPARAM];
284 lo_set_force_const(&plist[F_BONDS], c, 2, bRound, FALSE, bParam);
286 lo_set_force_const(&plist[F_ANGLES], c, 2, bRound, FALSE, bParam);
289 lo_set_force_const(&plist[F_PDIHS], c, 3, bRound, TRUE, bParam);
292 void calc_angles_dihs(t_params *ang, t_params *dih, rvec x[], gmx_bool bPBC,
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;
302 set_pbc(&pbc, epbcXYZ, box);
306 pr_rvecs(debug, 0, "X2TOP", box, DIM);
308 for (i = 0; (i < ang->nr); i++)
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);
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);
320 ang->param[i].C0 = th;
322 for (i = 0; (i < dih->nr); i++)
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);
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);
335 dih->param[i].C0 = ph;
339 static void dump_hybridization(FILE *fp, t_atoms *atoms, int nbonds[])
343 for (i = 0; (i < atoms->nr); i++)
345 fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
349 static void print_pl(FILE *fp, t_params plist[], int ftp, const char *name,
352 int i, j, nral, nrfp;
354 if (plist[ftp].nr > 0)
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++)
362 for (j = 0; (j < nral); j++)
364 fprintf(fp, " %5s", *atomname[plist[ftp].param[i].a[j]]);
366 for (j = 0; (j < nrfp); j++)
368 if (plist[ftp].param[i].c[j] != NOTSET)
370 fprintf(fp, " %10.3e", plist[ftp].param[i].c[j]);
378 static void print_rtp(const char *filenm, const char *title, t_atoms *atoms,
379 t_params plist[], gpp_atomtype_t atype, int cgnr[])
385 fp = gmx_fio_fopen(filenm, "w");
386 fprintf(fp, "; %s\n", title);
388 fprintf(fp, "[ %s ]\n", *atoms->resinfo[0].name);
390 fprintf(fp, "[ atoms ]\n");
391 for (i = 0; (i < atoms->nr); i++)
393 tp = atoms->atom[i].type;
394 if ((tpnm = get_atomtype_name(tp, atype)) == NULL)
396 gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i);
398 fprintf(fp, "%-8s %12s %8.4f %5d\n",
399 *atoms->atomname[i], tpnm,
400 atoms->atom[i].q, cgnr[i]);
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);
410 int gmx_x2top(int argc, char *argv[])
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]",
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."
440 t_params plist[F_NRE];
442 t_atoms *atoms; /* list with all atoms */
443 gpp_atomtype_t atype;
448 char title[STRLEN], forcefield[32], ffdir[STRLEN];
449 rvec *x; /* coordinates? */
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;
457 gmx_bool bRTP, bTOP, bOPLS;
459 real cutoff, qtot, mtot;
464 { efSTX, "-f", "conf", ffREAD },
465 { efTOP, "-o", "out", ffOPTWR },
466 { efRTP, "-r", "out", ffOPTWR }
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";
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)" }
513 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
514 asize(desc), desc, asize(bugs), bugs, &oenv))
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;
529 gmx_fatal(FARGS, "Specify at least one output file");
532 /* Force field selection, interactive or direct */
533 choose_ff(strcmp(ff, "select") == 0 ? NULL : ff,
534 forcefield, sizeof(forcefield),
535 ffdir, sizeof(ffdir));
537 bOPLS = (strcmp(forcefield, "oplsaa") == 0);
540 mymol.name = gmx_strdup(molnm);
543 /* Init parameter lists */
546 /* Read coordinates */
547 get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
550 /* make space for all the atoms */
551 init_t_atoms(atoms, natoms, TRUE);
554 read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, NULL, &epbc, box);
556 sprintf(n2t, "%s", ffdir);
557 nm2t = rd_nm2type(n2t, &nnm);
560 gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)",
565 printf("There are %d name to type translations in file %s\n", nnm, n2t);
569 dump_nm2type(debug, nnm, nm2t);
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);
575 open_symtab(&symtab);
576 atype = set_atom_type(&symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
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);
589 plist[F_LJ14].nr = 0;
592 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
593 " %4d pairs, %4d bonds and %4d atoms\n",
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);
599 calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
601 set_force_const(plist, kb, kt, kp, bRound, bParam);
603 cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
604 printf("Total charge is %g, total mass is %g\n", qtot, mtot);
613 fp = ftp2FILE(efTOP, NFILE, fnm, "w");
614 print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
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);
624 print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top",
625 atoms, plist, atype, cgnr);
630 dump_hybridization(debug, atoms, nbonds);
632 close_symtab(&symtab);
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");