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,2015,2016,2017, 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.
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/notset.h"
52 #include "gromacs/gmxpreprocess/pdb2top.h"
53 #include "gromacs/gmxpreprocess/toppush.h"
54 #include "gromacs/listed-forces/bonded.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/utilities.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/math/vecdump.h"
59 #include "gromacs/mdtypes/md_enums.h"
60 #include "gromacs/pbcutil/pbc.h"
61 #include "gromacs/topology/symtab.h"
62 #include "gromacs/topology/topology.h"
63 #include "gromacs/utility/arraysize.h"
64 #include "gromacs/utility/cstringutil.h"
65 #include "gromacs/utility/fatalerror.h"
66 #include "gromacs/utility/gmxassert.h"
67 #include "gromacs/utility/smalloc.h"
69 char atp[7] = "HCNOSX";
70 #define NATP (asize(atp)-1)
72 double blen[NATP][NATP] = {
73 { 0.00, 0.108, 0.105, 0.10, 0.10, 0.10 },
74 { 0.108, 0.15, 0.14, 0.14, 0.16, 0.14 },
75 { 0.105, 0.14, 0.14, 0.14, 0.16, 0.14 },
76 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.14 },
77 { 0.10, 0.16, 0.16, 0.17, 0.20, 0.17 },
78 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.17 }
81 #define MARGIN_FAC 1.1
83 static gmx_bool is_bond(int nnm, t_nm2type nmt[], char *ai, char *aj, real blen)
87 for (i = 0; (i < nnm); i++)
89 for (j = 0; (j < nmt[i].nbonds); j++)
91 if ((((gmx_strncasecmp(ai, nmt[i].elem, 1) == 0) &&
92 (gmx_strncasecmp(aj, nmt[i].bond[j], 1) == 0)) ||
93 ((gmx_strncasecmp(ai, nmt[i].bond[j], 1) == 0) &&
94 (gmx_strncasecmp(aj, nmt[i].elem, 1) == 0))) &&
95 (fabs(blen-nmt[i].blen[j]) <= 0.1*nmt[i].blen[j]))
104 void mk_bonds(int nnm, t_nm2type nmt[],
105 t_atoms *atoms, const rvec x[], t_params *bond, int nbond[],
106 gmx_bool bPBC, matrix box)
114 for (i = 0; (i < MAXATOMLIST); i++)
118 for (i = 0; (i < MAXFORCEPARAM); i++)
125 set_pbc(&pbc, -1, box);
127 for (i = 0; (i < atoms->nr); i++)
131 fprintf(stderr, "\ratom %d", i);
134 for (j = i+1; (j < atoms->nr); j++)
138 pbc_dx(&pbc, x[i], x[j], dx);
142 rvec_sub(x[i], x[j], dx);
146 if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
151 b.c0() = std::sqrt(dx2);
152 add_param_to_list (bond, &b);
157 fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
158 *atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
163 fprintf(stderr, "\ratom %d\n", i);
167 int *set_cgnr(t_atoms *atoms, gmx_bool bUsePDBcharge, real *qtot, real *mtot)
174 snew(cgnr, atoms->nr);
175 for (i = 0; (i < atoms->nr); i++)
177 if (atoms->pdbinfo && bUsePDBcharge)
179 atoms->atom[i].q = atoms->pdbinfo[i].bfac;
181 qt += atoms->atom[i].q;
182 *qtot += atoms->atom[i].q;
183 *mtot += atoms->atom[i].m;
194 gpp_atomtype_t set_atom_type(t_symtab *tab, t_atoms *atoms, t_params *bonds,
195 int *nbonds, int nnm, t_nm2type nm2t[])
197 gpp_atomtype_t atype;
200 atype = init_atomtype();
201 snew(atoms->atomtype, atoms->nr);
202 nresolved = nm2type(nnm, nm2t, tab, atoms, atype, nbonds, bonds);
203 if (nresolved != atoms->nr)
205 gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms",
206 nresolved, atoms->nr);
209 fprintf(stderr, "There are %d different atom types in your sample\n",
210 get_atomtype_ntypes(atype));
215 void lo_set_force_const(t_params *plist, real c[], int nrfp, gmx_bool bRound,
216 gmx_bool bDih, gmx_bool bParam)
222 for (i = 0; (i < plist->nr); i++)
226 for (j = 0; j < nrfp; j++)
235 sprintf(buf, "%.2e", plist->param[i].c[0]);
236 sscanf(buf, "%lf", &cc);
241 c[0] = plist->param[i].c[0];
246 c[0] = ((int)(c[0] + 3600)) % 360;
251 /* To put the minimum at the current angle rather than the maximum */
255 GMX_ASSERT(nrfp <= MAXFORCEPARAM/2, "Only 6 parameters may be used for an interaction");
256 for (j = 0; (j < nrfp); j++)
258 plist->param[i].c[j] = c[j];
259 plist->param[i].c[nrfp+j] = c[j];
261 set_p_string(&(plist->param[i]), "");
265 void set_force_const(t_params plist[], real kb, real kt, real kp, gmx_bool bRound,
268 real c[MAXFORCEPARAM];
272 lo_set_force_const(&plist[F_BONDS], c, 2, bRound, FALSE, bParam);
274 lo_set_force_const(&plist[F_ANGLES], c, 2, bRound, FALSE, bParam);
277 lo_set_force_const(&plist[F_PDIHS], c, 3, bRound, TRUE, bParam);
280 void calc_angles_dihs(t_params *ang, t_params *dih, const rvec x[], gmx_bool bPBC,
283 int i, ai, aj, ak, al, t1, t2, t3;
284 rvec r_ij, r_kj, r_kl, m, n;
285 real sign, th, costh, ph;
290 set_pbc(&pbc, epbcXYZ, box);
294 pr_rvecs(debug, 0, "X2TOP", box, DIM);
296 for (i = 0; (i < ang->nr); i++)
298 ai = ang->param[i].ai();
299 aj = ang->param[i].aj();
300 ak = ang->param[i].ak();
301 th = RAD2DEG*bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : nullptr,
302 r_ij, r_kj, &costh, &t1, &t2);
305 fprintf(debug, "X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
306 ai, aj, ak, norm(r_ij), norm(r_kj), th);
308 ang->param[i].c0() = th;
310 for (i = 0; (i < dih->nr); i++)
312 ai = dih->param[i].ai();
313 aj = dih->param[i].aj();
314 ak = dih->param[i].ak();
315 al = dih->param[i].al();
316 ph = RAD2DEG*dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : nullptr,
317 r_ij, r_kj, r_kl, m, n, &sign, &t1, &t2, &t3);
320 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",
321 ai, aj, ak, al, norm(r_ij), norm(r_kj), norm(r_kl), ph);
323 dih->param[i].c0() = ph;
327 static void dump_hybridization(FILE *fp, t_atoms *atoms, int nbonds[])
331 for (i = 0; (i < atoms->nr); i++)
333 fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
337 static void print_pl(FILE *fp, t_params plist[], int ftp, const char *name,
340 int i, j, nral, nrfp;
342 if (plist[ftp].nr > 0)
345 fprintf(fp, "[ %s ]\n", name);
346 nral = interaction_function[ftp].nratoms;
347 nrfp = interaction_function[ftp].nrfpA;
348 for (i = 0; (i < plist[ftp].nr); i++)
350 for (j = 0; (j < nral); j++)
352 fprintf(fp, " %5s", *atomname[plist[ftp].param[i].a[j]]);
354 for (j = 0; (j < nrfp); j++)
356 if (plist[ftp].param[i].c[j] != NOTSET)
358 fprintf(fp, " %10.3e", plist[ftp].param[i].c[j]);
366 static void print_rtp(const char *filenm, const char *title, t_atoms *atoms,
367 t_params plist[], gpp_atomtype_t atype, int cgnr[])
373 fp = gmx_fio_fopen(filenm, "w");
374 fprintf(fp, "; %s\n", title);
376 fprintf(fp, "[ %s ]\n", *atoms->resinfo[0].name);
378 fprintf(fp, "[ atoms ]\n");
379 for (i = 0; (i < atoms->nr); i++)
381 tp = atoms->atom[i].type;
382 if ((tpnm = get_atomtype_name(tp, atype)) == nullptr)
384 gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i);
386 fprintf(fp, "%-8s %12s %8.4f %5d\n",
387 *atoms->atomname[i], tpnm,
388 atoms->atom[i].q, cgnr[i]);
390 print_pl(fp, plist, F_BONDS, "bonds", atoms->atomname);
391 print_pl(fp, plist, F_ANGLES, "angles", atoms->atomname);
392 print_pl(fp, plist, F_PDIHS, "dihedrals", atoms->atomname);
393 print_pl(fp, plist, F_IDIHS, "impropers", atoms->atomname);
398 int gmx_x2top(int argc, char *argv[])
400 const char *desc[] = {
401 "[THISMODULE] generates a primitive topology from a coordinate file.",
402 "The program assumes all hydrogens are present when defining",
403 "the hybridization from the atom name and the number of bonds.",
404 "The program can also make an [REF].rtp[ref] entry, which you can then add",
405 "to the [REF].rtp[ref] database.[PAR]",
406 "When [TT]-param[tt] is set, equilibrium distances and angles",
407 "and force constants will be printed in the topology for all",
408 "interactions. The equilibrium distances and angles are taken",
409 "from the input coordinates, the force constant are set with",
410 "command line options.",
411 "The force fields somewhat supported currently are:[PAR]",
412 "G53a5 GROMOS96 53a5 Forcefield (official distribution)[PAR]",
413 "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
414 "The corresponding data files can be found in the library directory",
415 "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
416 "information about file formats. By default, the force field selection",
417 "is interactive, but you can use the [TT]-ff[tt] option to specify",
418 "one of the short names above on the command line instead. In that",
419 "case [THISMODULE] just looks for the corresponding file.[PAR]",
421 const char *bugs[] = {
422 "The atom type selection is primitive. Virtually no chemical knowledge is used",
423 "Periodic boundary conditions screw up the bonding",
424 "No improper dihedrals are generated",
425 "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."
428 t_params plist[F_NRE];
430 gpp_atomtype_t atype;
435 char forcefield[32], ffdir[STRLEN];
436 rvec *x; /* coordinates? */
438 int bts[] = { 1, 1, 1, 2 };
439 matrix box; /* box length matrix */
440 int natoms; /* number of atoms in one molecule */
442 gmx_bool bRTP, bTOP, bOPLS;
446 gmx_output_env_t *oenv;
449 { efSTX, "-f", "conf", ffREAD },
450 { efTOP, "-o", "out", ffOPTWR },
451 { efRTP, "-r", "out", ffOPTWR }
453 #define NFILE asize(fnm)
454 real kb = 4e5, kt = 400, kp = 5;
455 t_restp rtp_header_settings = { 0 };
456 gmx_bool bRemoveDihedralIfWithImproper = FALSE;
457 gmx_bool bGenerateHH14Interactions = TRUE;
458 gmx_bool bKeepAllGeneratedDihedrals = FALSE;
460 gmx_bool bParam = TRUE, bRound = TRUE;
461 gmx_bool bPairs = TRUE, bPBC = TRUE;
462 gmx_bool bUsePDBcharge = FALSE, bVerbose = FALSE;
463 const char *molnm = "ICE";
464 const char *ff = "oplsaa";
466 { "-ff", FALSE, etSTR, {&ff},
467 "Force field for your simulation. Type \"select\" for interactive selection." },
468 { "-v", FALSE, etBOOL, {&bVerbose},
469 "Generate verbose output in the top file." },
470 { "-nexcl", FALSE, etINT, {&nrexcl},
471 "Number of exclusions" },
472 { "-H14", FALSE, etBOOL, {&bGenerateHH14Interactions},
473 "Use 3rd neighbour interactions for hydrogen atoms" },
474 { "-alldih", FALSE, etBOOL, {&bKeepAllGeneratedDihedrals},
475 "Generate all proper dihedrals" },
476 { "-remdih", FALSE, etBOOL, {&bRemoveDihedralIfWithImproper},
477 "Remove dihedrals on the same bond as an improper" },
478 { "-pairs", FALSE, etBOOL, {&bPairs},
479 "Output 1-4 interactions (pairs) in topology file" },
480 { "-name", FALSE, etSTR, {&molnm},
481 "Name of your molecule" },
482 { "-pbc", FALSE, etBOOL, {&bPBC},
483 "Use periodic boundary conditions." },
484 { "-pdbq", FALSE, etBOOL, {&bUsePDBcharge},
485 "Use the B-factor supplied in a [REF].pdb[ref] file for the atomic charges" },
486 { "-param", FALSE, etBOOL, {&bParam},
487 "Print parameters in the output" },
488 { "-round", FALSE, etBOOL, {&bRound},
489 "Round off measured values" },
490 { "-kb", FALSE, etREAL, {&kb},
491 "Bonded force constant (kJ/mol/nm^2)" },
492 { "-kt", FALSE, etREAL, {&kt},
493 "Angle force constant (kJ/mol/rad^2)" },
494 { "-kp", FALSE, etREAL, {&kp},
495 "Dihedral angle force constant (kJ/mol/rad^2)" }
498 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
499 asize(desc), desc, asize(bugs), bugs, &oenv))
503 bRTP = opt2bSet("-r", NFILE, fnm);
504 bTOP = opt2bSet("-o", NFILE, fnm);
505 /* C89 requirements mean that these struct members cannot be used in
506 * the declaration of pa. So some temporary variables are needed. */
507 rtp_header_settings.bRemoveDihedralIfWithImproper = bRemoveDihedralIfWithImproper;
508 rtp_header_settings.bGenerateHH14Interactions = bGenerateHH14Interactions;
509 rtp_header_settings.bKeepAllGeneratedDihedrals = bKeepAllGeneratedDihedrals;
510 rtp_header_settings.nrexcl = nrexcl;
514 gmx_fatal(FARGS, "Specify at least one output file");
517 /* Force field selection, interactive or direct */
518 choose_ff(strcmp(ff, "select") == 0 ? nullptr : ff,
519 forcefield, sizeof(forcefield),
520 ffdir, sizeof(ffdir));
522 bOPLS = (strcmp(forcefield, "oplsaa") == 0);
525 mymol.name = gmx_strdup(molnm);
528 /* Init parameter lists */
531 /* Read coordinates */
534 read_tps_conf(opt2fn("-f", NFILE, fnm), top, &epbc, &x, nullptr, box, FALSE);
535 t_atoms *atoms = &top->atoms;
537 if (atoms->pdbinfo == nullptr)
539 snew(atoms->pdbinfo, natoms);
542 sprintf(n2t, "%s", ffdir);
543 nm2t = rd_nm2type(n2t, &nnm);
546 gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)",
551 printf("There are %d name to type translations in file %s\n", nnm, n2t);
555 dump_nm2type(debug, nnm, nm2t);
557 printf("Generating bonds from distances...\n");
558 snew(nbonds, atoms->nr);
559 mk_bonds(nnm, nm2t, atoms, x, &(plist[F_BONDS]), nbonds, bPBC, box);
561 open_symtab(&symtab);
562 atype = set_atom_type(&symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
564 /* Make Angles and Dihedrals */
565 snew(excls, atoms->nr);
566 printf("Generating angles and dihedrals from bonds...\n");
567 init_nnb(&nnb, atoms->nr, 4);
568 gen_nnb(&nnb, plist);
569 print_nnb(&nnb, "NNB");
570 gen_pad(&nnb, atoms, &rtp_header_settings, plist, excls, nullptr, TRUE);
575 plist[F_LJ14].nr = 0;
578 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
579 " %4d pairs, %4d bonds and %4d atoms\n",
581 bOPLS ? "Ryckaert-Bellemans" : "proper",
582 plist[F_IDIHS].nr, plist[F_ANGLES].nr,
583 plist[F_LJ14].nr, plist[F_BONDS].nr, atoms->nr);
585 calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
587 set_force_const(plist, kb, kt, kp, bRound, bParam);
589 cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
590 printf("Total charge is %g, total mass is %g\n", qtot, mtot);
599 fp = ftp2FILE(efTOP, NFILE, fnm, "w");
600 print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
602 write_top(fp, nullptr, mymol.name, atoms, FALSE, bts, plist, excls, atype,
603 cgnr, rtp_header_settings.nrexcl);
604 print_top_mols(fp, mymol.name, ffdir, nullptr, 0, nullptr, 1, &mymol);
610 print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top",
611 atoms, plist, atype, cgnr);
616 dump_hybridization(debug, atoms, nbonds);
618 close_symtab(&symtab);
621 printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",
622 output_env_get_program_display_name(oenv));
623 printf(" Please verify atomtypes and charges by comparison to other\n");
624 printf(" topologies.\n");