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.
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"
57 #include "gpp_nextnb.h"
58 #include "hackblock.h"
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"
70 char atp[7] = "HCNOSX";
71 #define NATP (asize(atp)-1)
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 }
82 #define MARGIN_FAC 1.1
84 static gmx_bool is_bond(int nnm, t_nm2type nmt[], char *ai, char *aj, real blen)
88 for (i = 0; (i < nnm); i++)
90 for (j = 0; (j < nmt[i].nbonds); j++)
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]))
105 static int get_atype(char *nm)
109 for (i = 0; (i < NATP-1); i++)
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)
130 for (i = 0; (i < MAXATOMLIST); i++)
134 for (i = 0; (i < MAXFORCEPARAM); i++)
141 set_pbc(&pbc, -1, box);
143 for (i = 0; (i < atoms->nr); i++)
147 fprintf(stderr, "\ratom %d", i);
149 for (j = i+1; (j < atoms->nr); j++)
153 pbc_dx(&pbc, x[i], x[j], dx);
157 rvec_sub(x[i], x[j], dx);
161 if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
167 add_param_to_list (bond, &b);
172 fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
173 *atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
178 fprintf(stderr, "\ratom %d\n", i);
181 int *set_cgnr(t_atoms *atoms, gmx_bool bUsePDBcharge, real *qtot, real *mtot)
185 double qt = 0, mt = 0;
188 snew(cgnr, atoms->nr);
189 for (i = 0; (i < atoms->nr); i++)
191 if (atoms->pdbinfo && bUsePDBcharge)
193 atoms->atom[i].q = atoms->pdbinfo[i].bfac;
195 qt += atoms->atom[i].q;
196 *qtot += atoms->atom[i].q;
197 *mtot += atoms->atom[i].m;
208 gpp_atomtype_t set_atom_type(t_symtab *tab, t_atoms *atoms, t_params *bonds,
209 int *nbonds, int nnm, t_nm2type nm2t[])
211 gpp_atomtype_t atype;
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)
219 gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms",
220 nresolved, atoms->nr);
223 fprintf(stderr, "There are %d different atom types in your sample\n",
224 get_atomtype_ntypes(atype));
229 void lo_set_force_const(t_params *plist, real c[], int nrfp, gmx_bool bRound,
230 gmx_bool bDih, gmx_bool bParam)
236 for (i = 0; (i < plist->nr); i++)
240 for (j = 0; j < nrfp; j++)
249 sprintf(buf, "%.2e", plist->param[i].c[0]);
250 sscanf(buf, "%lf", &cc);
255 c[0] = plist->param[i].c[0];
260 c[0] = ((int)(c[0] + 3600)) % 360;
265 /* To put the minimum at the current angle rather than the maximum */
269 assert(nrfp <= MAXFORCEPARAM/2);
270 for (j = 0; (j < nrfp); j++)
272 plist->param[i].c[j] = c[j];
273 plist->param[i].c[nrfp+j] = c[j];
275 set_p_string(&(plist->param[i]), "");
279 void set_force_const(t_params plist[], real kb, real kt, real kp, gmx_bool bRound,
283 real c[MAXFORCEPARAM];
287 lo_set_force_const(&plist[F_BONDS], c, 2, bRound, FALSE, bParam);
289 lo_set_force_const(&plist[F_ANGLES], c, 2, bRound, FALSE, bParam);
292 lo_set_force_const(&plist[F_PDIHS], c, 3, bRound, TRUE, bParam);
295 void calc_angles_dihs(t_params *ang, t_params *dih, rvec x[], gmx_bool bPBC,
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;
305 set_pbc(&pbc, epbcXYZ, box);
309 pr_rvecs(debug, 0, "X2TOP", box, DIM);
311 for (i = 0; (i < ang->nr); i++)
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);
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);
323 ang->param[i].C0 = th;
325 for (i = 0; (i < dih->nr); i++)
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);
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);
338 dih->param[i].C0 = ph;
342 static void dump_hybridization(FILE *fp, t_atoms *atoms, int nbonds[])
346 for (i = 0; (i < atoms->nr); i++)
348 fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
352 static void print_pl(FILE *fp, t_params plist[], int ftp, const char *name,
355 int i, j, nral, nrfp;
357 if (plist[ftp].nr > 0)
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++)
365 for (j = 0; (j < nral); j++)
367 fprintf(fp, " %5s", *atomname[plist[ftp].param[i].a[j]]);
369 for (j = 0; (j < nrfp); j++)
371 if (plist[ftp].param[i].c[j] != NOTSET)
373 fprintf(fp, " %10.3e", plist[ftp].param[i].c[j]);
381 static void print_rtp(const char *filenm, const char *title, t_atoms *atoms,
382 t_params plist[], gpp_atomtype_t atype, int cgnr[])
388 fp = gmx_fio_fopen(filenm, "w");
389 fprintf(fp, "; %s\n", title);
391 fprintf(fp, "[ %s ]\n", *atoms->resinfo[0].name);
393 fprintf(fp, "[ atoms ]\n");
394 for (i = 0; (i < atoms->nr); i++)
396 tp = atoms->atom[i].type;
397 if ((tpnm = get_atomtype_name(tp, atype)) == NULL)
399 gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i);
401 fprintf(fp, "%-8s %12s %8.4f %5d\n",
402 *atoms->atomname[i], tpnm,
403 atoms->atom[i].q, cgnr[i]);
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);
413 int gmx_x2top(int argc, char *argv[])
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]",
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."
443 t_params plist[F_NRE];
445 t_atoms *atoms; /* list with all atoms */
446 gpp_atomtype_t atype;
451 char title[STRLEN], forcefield[32], ffdir[STRLEN];
452 rvec *x; /* coordinates? */
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;
460 gmx_bool bRTP, bTOP, bOPLS;
462 real cutoff, qtot, mtot;
467 { efSTX, "-f", "conf", ffREAD },
468 { efTOP, "-o", "out", ffOPTWR },
469 { efRTP, "-r", "out", ffOPTWR }
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";
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)" }
516 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
517 asize(desc), desc, asize(bugs), bugs, &oenv))
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;
532 gmx_fatal(FARGS, "Specify at least one output file");
535 /* Force field selection, interactive or direct */
536 choose_ff(strcmp(ff, "select") == 0 ? NULL : ff,
537 forcefield, sizeof(forcefield),
538 ffdir, sizeof(ffdir));
540 bOPLS = (strcmp(forcefield, "oplsaa") == 0);
543 mymol.name = gmx_strdup(molnm);
546 /* Init parameter lists */
549 /* Read coordinates */
550 get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
553 /* make space for all the atoms */
554 init_t_atoms(atoms, natoms, TRUE);
557 read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, NULL, &epbc, box);
559 sprintf(n2t, "%s", ffdir);
560 nm2t = rd_nm2type(n2t, &nnm);
563 gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)",
568 printf("There are %d name to type translations in file %s\n", nnm, n2t);
572 dump_nm2type(debug, nnm, nm2t);
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);
578 open_symtab(&symtab);
579 atype = set_atom_type(&symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
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);
592 plist[F_LJ14].nr = 0;
595 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
596 " %4d pairs, %4d bonds and %4d atoms\n",
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);
602 calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
604 set_force_const(plist, kb, kt, kp, bRound, bParam);
606 cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
607 printf("Total charge is %g, total mass is %g\n", qtot, mtot);
616 fp = ftp2FILE(efTOP, NFILE, fnm, "w");
617 print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
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);
627 print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top",
628 atoms, plist, atype, cgnr);
633 dump_hybridization(debug, atoms, nbonds);
635 close_symtab(&symtab);
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");