3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Groningen Machine for Chemical Simulation
42 #include "gromacs/fileio/gmxfio.h"
47 #include "gromacs/fileio/confio.h"
59 #include "gpp_nextnb.h"
62 #include "hackblock.h"
66 char atp[7] = "HCNOSX";
67 #define NATP (asize(atp)-1)
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 }
78 #define MARGIN_FAC 1.1
80 static gmx_bool is_bond(int nnm, t_nm2type nmt[], char *ai, char *aj, real blen)
84 for (i = 0; (i < nnm); i++)
86 for (j = 0; (j < nmt[i].nbonds); j++)
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]))
101 static int get_atype(char *nm)
105 for (i = 0; (i < NATP-1); i++)
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)
126 for (i = 0; (i < MAXATOMLIST); i++)
130 for (i = 0; (i < MAXFORCEPARAM); i++)
137 set_pbc(&pbc, -1, box);
139 for (i = 0; (i < atoms->nr); i++)
143 fprintf(stderr, "\ratom %d", i);
145 for (j = i+1; (j < atoms->nr); j++)
149 pbc_dx(&pbc, x[i], x[j], dx);
153 rvec_sub(x[i], x[j], dx);
157 if (is_bond(nnm, nmt, *atoms->atomname[i], *atoms->atomname[j],
163 add_param_to_list (bond, &b);
168 fprintf(debug, "Bonding atoms %s-%d and %s-%d\n",
169 *atoms->atomname[i], i+1, *atoms->atomname[j], j+1);
174 fprintf(stderr, "\ratom %d\n", i);
177 int *set_cgnr(t_atoms *atoms, gmx_bool bUsePDBcharge, real *qtot, real *mtot)
181 double qt = 0, mt = 0;
184 snew(cgnr, atoms->nr);
185 for (i = 0; (i < atoms->nr); i++)
187 if (atoms->pdbinfo && bUsePDBcharge)
189 atoms->atom[i].q = atoms->pdbinfo[i].bfac;
191 qt += atoms->atom[i].q;
192 *qtot += atoms->atom[i].q;
193 *mtot += atoms->atom[i].m;
204 gpp_atomtype_t set_atom_type(t_symtab *tab, t_atoms *atoms, t_params *bonds,
205 int *nbonds, int nnm, t_nm2type nm2t[])
207 gpp_atomtype_t atype;
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)
215 gmx_fatal(FARGS, "Could only find a forcefield type for %d out of %d atoms",
216 nresolved, atoms->nr);
219 fprintf(stderr, "There are %d different atom types in your sample\n",
220 get_atomtype_ntypes(atype));
225 void lo_set_force_const(t_params *plist, real c[], int nrfp, gmx_bool bRound,
226 gmx_bool bDih, gmx_bool bParam)
232 for (i = 0; (i < plist->nr); i++)
236 for (j = 0; j < nrfp; j++)
245 sprintf(buf, "%.2e", plist->param[i].c[0]);
246 sscanf(buf, "%lf", &cc);
251 c[0] = plist->param[i].c[0];
256 c[0] = ((int)(c[0] + 3600)) % 360;
261 /* To put the minimum at the current angle rather than the maximum */
265 for (j = 0; (j < nrfp); j++)
267 plist->param[i].c[j] = c[j];
268 plist->param[i].c[nrfp+j] = c[j];
270 set_p_string(&(plist->param[i]), "");
274 void set_force_const(t_params plist[], real kb, real kt, real kp, gmx_bool bRound,
278 real c[MAXFORCEPARAM];
282 lo_set_force_const(&plist[F_BONDS], c, 2, bRound, FALSE, bParam);
284 lo_set_force_const(&plist[F_ANGLES], c, 2, bRound, FALSE, bParam);
287 lo_set_force_const(&plist[F_PDIHS], c, 3, bRound, TRUE, bParam);
290 void calc_angles_dihs(t_params *ang, t_params *dih, rvec x[], gmx_bool bPBC,
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;
300 set_pbc(&pbc, epbcXYZ, box);
304 pr_rvecs(debug, 0, "X2TOP", box, DIM);
306 for (i = 0; (i < ang->nr); i++)
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);
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);
318 ang->param[i].C0 = th;
320 for (i = 0; (i < dih->nr); i++)
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);
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);
333 dih->param[i].C0 = ph;
337 static void dump_hybridization(FILE *fp, t_atoms *atoms, int nbonds[])
341 for (i = 0; (i < atoms->nr); i++)
343 fprintf(fp, "Atom %5s has %1d bonds\n", *atoms->atomname[i], nbonds[i]);
347 static void print_pl(FILE *fp, t_params plist[], int ftp, const char *name,
350 int i, j, nral, nrfp;
352 if (plist[ftp].nr > 0)
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++)
360 for (j = 0; (j < nral); j++)
362 fprintf(fp, " %5s", *atomname[plist[ftp].param[i].a[j]]);
364 for (j = 0; (j < nrfp); j++)
366 if (plist[ftp].param[i].c[j] != NOTSET)
368 fprintf(fp, " %10.3e", plist[ftp].param[i].c[j]);
376 static void print_rtp(const char *filenm, const char *title, t_atoms *atoms,
377 t_params plist[], gpp_atomtype_t atype, int cgnr[])
383 fp = gmx_fio_fopen(filenm, "w");
384 fprintf(fp, "; %s\n", title);
386 fprintf(fp, "[ %s ]\n", *atoms->resinfo[0].name);
388 fprintf(fp, "[ atoms ]\n");
389 for (i = 0; (i < atoms->nr); i++)
391 tp = atoms->atom[i].type;
392 if ((tpnm = get_atomtype_name(tp, atype)) == NULL)
394 gmx_fatal(FARGS, "tp = %d, i = %d in print_rtp", tp, i);
396 fprintf(fp, "%-8s %12s %8.4f %5d\n",
397 *atoms->atomname[i], tpnm,
398 atoms->atom[i].q, cgnr[i]);
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);
408 int gmx_x2top(int argc, char *argv[])
410 const char *desc[] = {
411 "[THISMODULE] 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 [THISMODULE] just looks for the corresponding file.[PAR]",
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."
438 t_params plist[F_NRE];
440 t_atoms *atoms; /* list with all atoms */
441 gpp_atomtype_t atype;
446 char title[STRLEN], forcefield[32], ffdir[STRLEN];
447 rvec *x; /* coordinates? */
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;
455 gmx_bool bRTP, bTOP, bOPLS;
457 real cutoff, qtot, mtot;
462 { efSTX, "-f", "conf", ffREAD },
463 { efTOP, "-o", "out", ffOPTWR },
464 { efRTP, "-r", "out", ffOPTWR }
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";
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)" }
511 if (!parse_common_args(&argc, argv, 0, NFILE, fnm, asize(pa), pa,
512 asize(desc), desc, asize(bugs), bugs, &oenv))
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;
527 gmx_fatal(FARGS, "Specify at least one output file");
530 /* Force field selection, interactive or direct */
531 choose_ff(strcmp(ff, "select") == 0 ? NULL : ff,
532 forcefield, sizeof(forcefield),
533 ffdir, sizeof(ffdir));
535 bOPLS = (strcmp(forcefield, "oplsaa") == 0);
538 mymol.name = strdup(molnm);
541 /* Init parameter lists */
544 /* Read coordinates */
545 get_stx_coordnum(opt2fn("-f", NFILE, fnm), &natoms);
548 /* make space for all the atoms */
549 init_t_atoms(atoms, natoms, TRUE);
552 read_stx_conf(opt2fn("-f", NFILE, fnm), title, atoms, x, NULL, &epbc, box);
554 sprintf(n2t, "%s", ffdir);
555 nm2t = rd_nm2type(n2t, &nnm);
558 gmx_fatal(FARGS, "No or incorrect atomname2type.n2t file found (looking for %s)",
563 printf("There are %d name to type translations in file %s\n", nnm, n2t);
567 dump_nm2type(debug, nnm, nm2t);
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);
573 open_symtab(&symtab);
574 atype = set_atom_type(&symtab, atoms, &(plist[F_BONDS]), nbonds, nnm, nm2t);
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);
587 plist[F_LJ14].nr = 0;
590 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
591 " %4d pairs, %4d bonds and %4d atoms\n",
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);
597 calc_angles_dihs(&plist[F_ANGLES], &plist[F_PDIHS], x, bPBC, box);
599 set_force_const(plist, kb, kt, kp, bRound, bParam);
601 cgnr = set_cgnr(atoms, bUsePDBcharge, &qtot, &mtot);
602 printf("Total charge is %g, total mass is %g\n", qtot, mtot);
611 fp = ftp2FILE(efTOP, NFILE, fnm, "w");
612 print_top_header(fp, ftp2fn(efTOP, NFILE, fnm), TRUE, ffdir, 1.0);
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);
622 print_rtp(ftp2fn(efRTP, NFILE, fnm), "Generated by x2top",
623 atoms, plist, atype, cgnr);
628 dump_hybridization(debug, atoms, nbonds);
630 close_symtab(&symtab);
633 printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",
635 printf(" Please verify atomtypes and charges by comparison to other\n");
636 printf(" topologies.\n");