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.
46 #include "gromacs/math/utilities.h"
49 #include "gromacs/fileio/gmxfio.h"
50 #include "gromacs/utility/cstringutil.h"
51 #include "gromacs/utility/smalloc.h"
53 #include "gromacs/fileio/confio.h"
55 #include "gromacs/commandline/pargs.h"
57 #include "gromacs/math/3dview.h"
64 #include "gpp_nextnb.h"
67 #include "hackblock.h"
70 #include "gmx_fatal.h"
72 char atp[7] = "HCNOSX";
73 #define NATP (asize(atp)-1)
75 real blen[NATP][NATP] = {
76 { 0.00, 0.108, 0.105, 0.10, 0.10, 0.10 },
77 { 0.108, 0.15, 0.14, 0.14, 0.16, 0.14 },
78 { 0.105, 0.14, 0.14, 0.14, 0.16, 0.14 },
79 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.14 },
80 { 0.10, 0.16, 0.16, 0.17, 0.20, 0.17 },
81 { 0.10, 0.14, 0.14, 0.14, 0.17, 0.17 }
84 #define MARGIN_FAC 1.1
86 static gmx_bool is_bond(int nnm, t_nm2type nmt[], char *ai, char *aj, real blen)
90 for (i = 0; (i < nnm); i++)
92 for (j = 0; (j < nmt[i].nbonds); j++)
94 if ((((gmx_strncasecmp(ai, nmt[i].elem, 1) == 0) &&
95 (gmx_strncasecmp(aj, nmt[i].bond[j], 1) == 0)) ||
96 ((gmx_strncasecmp(ai, nmt[i].bond[j], 1) == 0) &&
97 (gmx_strncasecmp(aj, nmt[i].elem, 1) == 0))) &&
98 (fabs(blen-nmt[i].blen[j]) <= 0.1*nmt[i].blen[j]))
107 static int get_atype(char *nm)
111 for (i = 0; (i < NATP-1); i++)
122 void mk_bonds(int nnm, t_nm2type nmt[],
123 t_atoms *atoms, rvec x[], t_params *bond, int nbond[],
124 gmx_bool bPBC, matrix box)
132 for (i = 0; (i < MAXATOMLIST); i++)
136 for (i = 0; (i < MAXFORCEPARAM); i++)
143 set_pbc(&pbc, -1, box);
145 for (i = 0; (i < atoms->nr); i++)
149 fprintf(stderr, "\ratom %d", i);
151 for (j = i+1; (j < atoms->nr); j++)
155 pbc_dx(&pbc, x[i], x[j], dx);
159 rvec_sub(x[i], x[j], dx);
163 if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
169 add_param_to_list (bond, &b);
174 fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
175 *atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
180 fprintf(stderr, "\ratom %d\n", i);
183 int *set_cgnr(t_atoms *atoms, gmx_bool bUsePDBcharge, real *qtot, real *mtot)
187 double qt = 0, mt = 0;
190 snew(cgnr, atoms->nr);
191 for (i = 0; (i < atoms->nr); i++)
193 if (atoms->pdbinfo && bUsePDBcharge)
195 atoms->atom[i].q = atoms->pdbinfo[i].bfac;
197 qt += atoms->atom[i].q;
198 *qtot += atoms->atom[i].q;
199 *mtot += atoms->atom[i].m;
210 gpp_atomtype_t set_atom_type(t_symtab *tab, t_atoms *atoms, t_params *bonds,
211 int *nbonds, int nnm, t_nm2type nm2t[])
213 gpp_atomtype_t atype;
216 atype = init_atomtype();
217 snew(atoms->atomtype, atoms->nr);
218 nresolved = nm2type(nnm, nm2t, tab, atoms, atype, nbonds, bonds);
219 if (nresolved != atoms->nr)
221 gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms",
222 nresolved, atoms->nr);
225 fprintf(stderr, "There are %d different atom types in your sample\n",
226 get_atomtype_ntypes(atype));
231 void lo_set_force_const(t_params *plist, real c[], int nrfp, gmx_bool bRound,
232 gmx_bool bDih, gmx_bool bParam)
238 for (i = 0; (i < plist->nr); i++)
242 for (j = 0; j < nrfp; j++)
251 sprintf(buf, "%.2e", plist->param[i].c[0]);
252 sscanf(buf, "%lf", &cc);
257 c[0] = plist->param[i].c[0];
262 c[0] = ((int)(c[0] + 3600)) % 360;
267 /* To put the minimum at the current angle rather than the maximum */
271 assert(nrfp <= MAXFORCEPARAM/2);
272 for (j = 0; (j < nrfp); j++)
274 plist->param[i].c[j] = c[j];
275 plist->param[i].c[nrfp+j] = c[j];
277 set_p_string(&(plist->param[i]), "");
281 void set_force_const(t_params plist[], real kb, real kt, real kp, gmx_bool bRound,
285 real c[MAXFORCEPARAM];
289 lo_set_force_const(&plist[F_BONDS], c, 2, bRound, FALSE, bParam);
291 lo_set_force_const(&plist[F_ANGLES], c, 2, bRound, FALSE, bParam);
294 lo_set_force_const(&plist[F_PDIHS], c, 3, bRound, TRUE, bParam);
297 void calc_angles_dihs(t_params *ang, t_params *dih, rvec x[], gmx_bool bPBC,
300 int i, ai, aj, ak, al, t1, t2, t3;
301 rvec r_ij, r_kj, r_kl, m, n;
302 real sign, th, costh, ph;
307 set_pbc(&pbc, epbcXYZ, box);
311 pr_rvecs(debug, 0, "X2TOP", box, DIM);
313 for (i = 0; (i < ang->nr); i++)
315 ai = ang->param[i].AI;
316 aj = ang->param[i].AJ;
317 ak = ang->param[i].AK;
318 th = RAD2DEG*bond_angle(x[ai], x[aj], x[ak], bPBC ? &pbc : NULL,
319 r_ij, r_kj, &costh, &t1, &t2);
322 fprintf(debug, "X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
323 ai, aj, ak, norm(r_ij), norm(r_kj), th);
325 ang->param[i].C0 = th;
327 for (i = 0; (i < dih->nr); i++)
329 ai = dih->param[i].AI;
330 aj = dih->param[i].AJ;
331 ak = dih->param[i].AK;
332 al = dih->param[i].AL;
333 ph = RAD2DEG*dih_angle(x[ai], x[aj], x[ak], x[al], bPBC ? &pbc : NULL,
334 r_ij, r_kj, r_kl, m, n, &sign, &t1, &t2, &t3);
337 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",
338 ai, aj, ak, al, norm(r_ij), norm(r_kj), norm(r_kl), ph);
340 dih->param[i].C0 = ph;
344 static void dump_hybridization(FILE *fp, t_atoms *atoms, int nbonds[])
348 for (i = 0; (i < atoms->nr); i++)
350 fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
354 static void print_pl(FILE *fp, t_params plist[], int ftp, const char *name,
357 int i, j, nral, nrfp;
359 if (plist[ftp].nr > 0)
362 fprintf(fp, "[ %s ]\n", name);
363 nral = interaction_function[ftp].nratoms;
364 nrfp = interaction_function[ftp].nrfpA;
365 for (i = 0; (i < plist[ftp].nr); i++)
367 for (j = 0; (j < nral); j++)
369 fprintf(fp, " %5s", *atomname[plist[ftp].param[i].a[j]]);
371 for (j = 0; (j < nrfp); j++)
373 if (plist[ftp].param[i].c[j] != NOTSET)
375 fprintf(fp, " %10.3e", plist[ftp].param[i].c[j]);
383 static void print_rtp(const char *filenm, const char *title, t_atoms *atoms,
384 t_params plist[], gpp_atomtype_t atype, int cgnr[])
390 fp = gmx_fio_fopen(filenm, "w");
391 fprintf(fp, "; %s\n", title);
393 fprintf(fp, "[ %s ]\n", *atoms->resinfo[0].name);
395 fprintf(fp, "[ atoms ]\n");
396 for (i = 0; (i < atoms->nr); i++)
398 tp = atoms->atom[i].type;
399 if ((tpnm = get_atomtype_name(tp, atype)) == NULL)
401 gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i);
403 fprintf(fp, "%-8s %12s %8.4f %5d\n",
404 *atoms->atomname[i], tpnm,
405 atoms->atom[i].q, cgnr[i]);
407 print_pl(fp, plist, F_BONDS, "bonds", atoms->atomname);
408 print_pl(fp, plist, F_ANGLES, "angles", atoms->atomname);
409 print_pl(fp, plist, F_PDIHS, "dihedrals", atoms->atomname);
410 print_pl(fp, plist, F_IDIHS, "impropers", atoms->atomname);
415 int gmx_x2top(int argc, char *argv[])
417 const char *desc[] = {
418 "[THISMODULE] generates a primitive topology from a coordinate file.",
419 "The program assumes all hydrogens are present when defining",
420 "the hybridization from the atom name and the number of bonds.",
421 "The program can also make an [TT].rtp[tt] entry, which you can then add",
422 "to the [TT].rtp[tt] database.[PAR]",
423 "When [TT]-param[tt] is set, equilibrium distances and angles",
424 "and force constants will be printed in the topology for all",
425 "interactions. The equilibrium distances and angles are taken",
426 "from the input coordinates, the force constant are set with",
427 "command line options.",
428 "The force fields somewhat supported currently are:[PAR]",
429 "G53a5 GROMOS96 53a5 Forcefield (official distribution)[PAR]",
430 "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
431 "The corresponding data files can be found in the library directory",
432 "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
433 "information about file formats. By default, the force field selection",
434 "is interactive, but you can use the [TT]-ff[tt] option to specify",
435 "one of the short names above on the command line instead. In that",
436 "case [THISMODULE] just looks for the corresponding file.[PAR]",
438 const char *bugs[] = {
439 "The atom type selection is primitive. Virtually no chemical knowledge is used",
440 "Periodic boundary conditions screw up the bonding",
441 "No improper dihedrals are generated",
442 "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."
445 t_params plist[F_NRE];
447 t_atoms *atoms; /* list with all atoms */
448 gpp_atomtype_t atype;
453 char title[STRLEN], forcefield[32], ffdir[STRLEN];
454 rvec *x; /* coordinates? */
456 int bts[] = { 1, 1, 1, 2 };
457 matrix box; /* box length matrix */
458 int natoms; /* number of atoms in one molecule */
459 int nres; /* number of molecules? */
460 int i, j, k, l, m, ndih;
462 gmx_bool bRTP, bTOP, bOPLS;
464 real cutoff, qtot, mtot;
469 { efSTX, "-f", "conf", ffREAD },
470 { efTOP, "-o", "out", ffOPTWR },
471 { efRTP, "-r", "out", ffOPTWR }
473 #define NFILE asize(fnm)
474 static real scale = 1.1, kb = 4e5, kt = 400, kp = 5;
475 static t_restp rtp_header_settings;
476 static gmx_bool bRemoveDihedralIfWithImproper = FALSE;
477 static gmx_bool bGenerateHH14Interactions = TRUE;
478 static gmx_bool bKeepAllGeneratedDihedrals = FALSE;
479 static int nrexcl = 3;
480 static gmx_bool bParam = TRUE, bRound = TRUE;
481 static gmx_bool bPairs = TRUE, bPBC = TRUE;
482 static gmx_bool bUsePDBcharge = FALSE, bVerbose = FALSE;
483 static const char *molnm = "ICE";
484 static const char *ff = "oplsaa";
486 { "-ff", FALSE, etSTR, {&ff},
487 "Force field for your simulation. Type \"select\" for interactive selection." },
488 { "-v", FALSE, etBOOL, {&bVerbose},
489 "Generate verbose output in the top file." },
490 { "-nexcl", FALSE, etINT, {&nrexcl},
491 "Number of exclusions" },
492 { "-H14", FALSE, etBOOL, {&bGenerateHH14Interactions},
493 "Use 3rd neighbour interactions for hydrogen atoms" },
494 { "-alldih", FALSE, etBOOL, {&bKeepAllGeneratedDihedrals},
495 "Generate all proper dihedrals" },
496 { "-remdih", FALSE, etBOOL, {&bRemoveDihedralIfWithImproper},
497 "Remove dihedrals on the same bond as an improper" },
498 { "-pairs", FALSE, etBOOL, {&bPairs},
499 "Output 1-4 interactions (pairs) in topology file" },
500 { "-name", FALSE, etSTR, {&molnm},
501 "Name of your molecule" },
502 { "-pbc", FALSE, etBOOL, {&bPBC},
503 "Use periodic boundary conditions." },
504 { "-pdbq", FALSE, etBOOL, {&bUsePDBcharge},
505 "Use the B-factor supplied in a [TT].pdb[tt] file for the atomic charges" },
506 { "-param", FALSE, etBOOL, {&bParam},
507 "Print parameters in the output" },
508 { "-round", FALSE, etBOOL, {&bRound},
509 "Round off measured values" },
510 { "-kb", FALSE, etREAL, {&kb},
511 "Bonded force constant (kJ/mol/nm^2)" },
512 { "-kt", FALSE, etREAL, {&kt},
513 "Angle force constant (kJ/mol/rad^2)" },
514 { "-kp", FALSE, etREAL, {&kp},
515 "Dihedral angle force constant (kJ/mol/rad^2)" }
518 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
519 asize(desc), desc, asize(bugs), bugs, &oenv))
523 bRTP = opt2bSet("-r", NFILE, fnm);
524 bTOP = opt2bSet("-o", NFILE, fnm);
525 /* C89 requirements mean that these struct members cannot be used in
526 * the declaration of pa. So some temporary variables are needed. */
527 rtp_header_settings.bRemoveDihedralIfWithImproper = bRemoveDihedralIfWithImproper;
528 rtp_header_settings.bGenerateHH14Interactions = bGenerateHH14Interactions;
529 rtp_header_settings.bKeepAllGeneratedDihedrals = bKeepAllGeneratedDihedrals;
530 rtp_header_settings.nrexcl = nrexcl;
534 gmx_fatal(FARGS, "Specify at least one output file");
537 /* Force field selection, interactive or direct */
538 choose_ff(strcmp(ff, "select") == 0 ? NULL : ff,
539 forcefield, sizeof(forcefield),
540 ffdir, sizeof(ffdir));
542 bOPLS = (strcmp(forcefield, "oplsaa") == 0);
545 mymol.name = strdup(molnm);
548 /* Init parameter lists */
551 /* Read coordinates */
552 get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
555 /* make space for all the atoms */
556 init_t_atoms(atoms, natoms, TRUE);
559 read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, NULL, &epbc, box);
561 sprintf(n2t, "%s", ffdir);
562 nm2t = rd_nm2type(n2t, &nnm);
565 gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)",
570 printf("There are %d name to type translations in file %s\n", nnm, n2t);
574 dump_nm2type(debug, nnm, nm2t);
576 printf("Generating bonds from distances...\n");
577 snew(nbonds, atoms->nr);
578 mk_bonds(nnm, nm2t, atoms, x, &(plist[F_BONDS]), nbonds, bPBC, box);
580 open_symtab(&symtab);
581 atype = set_atom_type(&symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
583 /* Make Angles and Dihedrals */
584 snew(excls, atoms->nr);
585 printf("Generating angles and dihedrals from bonds...\n");
586 init_nnb(&nnb, atoms->nr, 4);
587 gen_nnb(&nnb, plist);
588 print_nnb(&nnb, "NNB");
589 gen_pad(&nnb, atoms, &rtp_header_settings, plist, excls, NULL, TRUE);
594 plist[F_LJ14].nr = 0;
597 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
598 " %4d pairs, %4d bonds and %4d atoms\n",
600 bOPLS ? "Ryckaert-Bellemans" : "proper",
601 plist[F_IDIHS].nr, plist[F_ANGLES].nr,
602 plist[F_LJ14].nr, plist[F_BONDS].nr, atoms->nr);
604 calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
606 set_force_const(plist, kb, kt, kp, bRound, bParam);
608 cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
609 printf("Total charge is %g, total mass is %g\n", qtot, mtot);
618 fp = ftp2FILE(efTOP, NFILE, fnm, "w");
619 print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
621 write_top(fp, NULL, mymol.name, atoms, FALSE, bts, plist, excls, atype,
622 cgnr, rtp_header_settings.nrexcl);
623 print_top_mols(fp, mymol.name, ffdir, NULL, 0, NULL, 1, &mymol);
629 print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top",
630 atoms, plist, atype, cgnr);
635 dump_hybridization(debug, atoms, nbonds);
637 close_symtab(&symtab);
640 printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",
642 printf(" Please verify atomtypes and charges by comparison to other\n");
643 printf(" topologies.\n");