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, 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/gmxlib/readinp.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/legacyheaders/copyrite.h"
55 #include "gromacs/legacyheaders/names.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/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 real 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 static int get_atype(char *nm
)
108 for (i
= 0; (i
< NATP
-1); i
++)
119 void mk_bonds(int nnm
, t_nm2type nmt
[],
120 t_atoms
*atoms
, rvec x
[], t_params
*bond
, int nbond
[],
121 gmx_bool bPBC
, matrix box
)
129 for (i
= 0; (i
< MAXATOMLIST
); i
++)
133 for (i
= 0; (i
< MAXFORCEPARAM
); i
++)
140 set_pbc(&pbc
, -1, box
);
142 for (i
= 0; (i
< atoms
->nr
); i
++)
146 fprintf(stderr
, "\ratom %d", i
);
148 for (j
= i
+1; (j
< atoms
->nr
); j
++)
152 pbc_dx(&pbc
, x
[i
], x
[j
], dx
);
156 rvec_sub(x
[i
], x
[j
], dx
);
160 if (is_bond(nnm
, nmt
, *atoms
->atomname
[i
], *atoms
->atomname
[j
],
165 b
.c0() = std::sqrt(dx2
);
166 add_param_to_list (bond
, &b
);
171 fprintf(debug
, "Bonding atoms %s-%d and %s-%d\n",
172 *atoms
->atomname
[i
], i
+1, *atoms
->atomname
[j
], j
+1);
177 fprintf(stderr
, "\ratom %d\n", i
);
180 int *set_cgnr(t_atoms
*atoms
, gmx_bool bUsePDBcharge
, real
*qtot
, real
*mtot
)
187 snew(cgnr
, atoms
->nr
);
188 for (i
= 0; (i
< atoms
->nr
); i
++)
190 if (atoms
->pdbinfo
&& bUsePDBcharge
)
192 atoms
->atom
[i
].q
= atoms
->pdbinfo
[i
].bfac
;
194 qt
+= atoms
->atom
[i
].q
;
195 *qtot
+= atoms
->atom
[i
].q
;
196 *mtot
+= atoms
->atom
[i
].m
;
207 gpp_atomtype_t
set_atom_type(t_symtab
*tab
, t_atoms
*atoms
, t_params
*bonds
,
208 int *nbonds
, int nnm
, t_nm2type nm2t
[])
210 gpp_atomtype_t atype
;
213 atype
= init_atomtype();
214 snew(atoms
->atomtype
, atoms
->nr
);
215 nresolved
= nm2type(nnm
, nm2t
, tab
, atoms
, atype
, nbonds
, bonds
);
216 if (nresolved
!= atoms
->nr
)
218 gmx_fatal(FARGS
, "Could only find a forcefield type for %d out of %d atoms",
219 nresolved
, atoms
->nr
);
222 fprintf(stderr
, "There are %d different atom types in your sample\n",
223 get_atomtype_ntypes(atype
));
228 void lo_set_force_const(t_params
*plist
, real c
[], int nrfp
, gmx_bool bRound
,
229 gmx_bool bDih
, gmx_bool bParam
)
235 for (i
= 0; (i
< plist
->nr
); i
++)
239 for (j
= 0; j
< nrfp
; j
++)
248 sprintf(buf
, "%.2e", plist
->param
[i
].c
[0]);
249 sscanf(buf
, "%lf", &cc
);
254 c
[0] = plist
->param
[i
].c
[0];
259 c
[0] = ((int)(c
[0] + 3600)) % 360;
264 /* To put the minimum at the current angle rather than the maximum */
268 GMX_ASSERT(nrfp
<= MAXFORCEPARAM
/2, "Only 6 parameters may be used for an interaction");
269 for (j
= 0; (j
< nrfp
); j
++)
271 plist
->param
[i
].c
[j
] = c
[j
];
272 plist
->param
[i
].c
[nrfp
+j
] = c
[j
];
274 set_p_string(&(plist
->param
[i
]), "");
278 void set_force_const(t_params plist
[], real kb
, real kt
, real kp
, gmx_bool bRound
,
281 real c
[MAXFORCEPARAM
];
285 lo_set_force_const(&plist
[F_BONDS
], c
, 2, bRound
, FALSE
, bParam
);
287 lo_set_force_const(&plist
[F_ANGLES
], c
, 2, bRound
, FALSE
, bParam
);
290 lo_set_force_const(&plist
[F_PDIHS
], c
, 3, bRound
, TRUE
, bParam
);
293 void calc_angles_dihs(t_params
*ang
, t_params
*dih
, rvec x
[], gmx_bool bPBC
,
296 int i
, ai
, aj
, ak
, al
, t1
, t2
, t3
;
297 rvec r_ij
, r_kj
, r_kl
, m
, n
;
298 real sign
, th
, costh
, ph
;
303 set_pbc(&pbc
, epbcXYZ
, box
);
307 pr_rvecs(debug
, 0, "X2TOP", box
, DIM
);
309 for (i
= 0; (i
< ang
->nr
); i
++)
311 ai
= ang
->param
[i
].ai();
312 aj
= ang
->param
[i
].aj();
313 ak
= ang
->param
[i
].ak();
314 th
= RAD2DEG
*bond_angle(x
[ai
], x
[aj
], x
[ak
], bPBC
? &pbc
: NULL
,
315 r_ij
, r_kj
, &costh
, &t1
, &t2
);
318 fprintf(debug
, "X2TOP: ai=%3d aj=%3d ak=%3d r_ij=%8.3f r_kj=%8.3f th=%8.3f\n",
319 ai
, aj
, ak
, norm(r_ij
), norm(r_kj
), th
);
321 ang
->param
[i
].c0() = th
;
323 for (i
= 0; (i
< dih
->nr
); i
++)
325 ai
= dih
->param
[i
].ai();
326 aj
= dih
->param
[i
].aj();
327 ak
= dih
->param
[i
].ak();
328 al
= dih
->param
[i
].al();
329 ph
= RAD2DEG
*dih_angle(x
[ai
], x
[aj
], x
[ak
], x
[al
], bPBC
? &pbc
: NULL
,
330 r_ij
, r_kj
, r_kl
, m
, n
, &sign
, &t1
, &t2
, &t3
);
333 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",
334 ai
, aj
, ak
, al
, norm(r_ij
), norm(r_kj
), norm(r_kl
), ph
);
336 dih
->param
[i
].c0() = ph
;
340 static void dump_hybridization(FILE *fp
, t_atoms
*atoms
, int nbonds
[])
344 for (i
= 0; (i
< atoms
->nr
); i
++)
346 fprintf(fp
, "Atom %5s has %1d bonds\n", *atoms
->atomname
[i
], nbonds
[i
]);
350 static void print_pl(FILE *fp
, t_params plist
[], int ftp
, const char *name
,
353 int i
, j
, nral
, nrfp
;
355 if (plist
[ftp
].nr
> 0)
358 fprintf(fp
, "[ %s ]\n", name
);
359 nral
= interaction_function
[ftp
].nratoms
;
360 nrfp
= interaction_function
[ftp
].nrfpA
;
361 for (i
= 0; (i
< plist
[ftp
].nr
); i
++)
363 for (j
= 0; (j
< nral
); j
++)
365 fprintf(fp
, " %5s", *atomname
[plist
[ftp
].param
[i
].a
[j
]]);
367 for (j
= 0; (j
< nrfp
); j
++)
369 if (plist
[ftp
].param
[i
].c
[j
] != NOTSET
)
371 fprintf(fp
, " %10.3e", plist
[ftp
].param
[i
].c
[j
]);
379 static void print_rtp(const char *filenm
, const char *title
, t_atoms
*atoms
,
380 t_params plist
[], gpp_atomtype_t atype
, int cgnr
[])
386 fp
= gmx_fio_fopen(filenm
, "w");
387 fprintf(fp
, "; %s\n", title
);
389 fprintf(fp
, "[ %s ]\n", *atoms
->resinfo
[0].name
);
391 fprintf(fp
, "[ atoms ]\n");
392 for (i
= 0; (i
< atoms
->nr
); i
++)
394 tp
= atoms
->atom
[i
].type
;
395 if ((tpnm
= get_atomtype_name(tp
, atype
)) == NULL
)
397 gmx_fatal(FARGS
, "tp = %d, i = %d in print_rtp", tp
, i
);
399 fprintf(fp
, "%-8s %12s %8.4f %5d\n",
400 *atoms
->atomname
[i
], tpnm
,
401 atoms
->atom
[i
].q
, cgnr
[i
]);
403 print_pl(fp
, plist
, F_BONDS
, "bonds", atoms
->atomname
);
404 print_pl(fp
, plist
, F_ANGLES
, "angles", atoms
->atomname
);
405 print_pl(fp
, plist
, F_PDIHS
, "dihedrals", atoms
->atomname
);
406 print_pl(fp
, plist
, F_IDIHS
, "impropers", atoms
->atomname
);
411 int gmx_x2top(int argc
, char *argv
[])
413 const char *desc
[] = {
414 "[THISMODULE] generates a primitive topology from a coordinate file.",
415 "The program assumes all hydrogens are present when defining",
416 "the hybridization from the atom name and the number of bonds.",
417 "The program can also make an [REF].rtp[ref] entry, which you can then add",
418 "to the [REF].rtp[ref] database.[PAR]",
419 "When [TT]-param[tt] is set, equilibrium distances and angles",
420 "and force constants will be printed in the topology for all",
421 "interactions. The equilibrium distances and angles are taken",
422 "from the input coordinates, the force constant are set with",
423 "command line options.",
424 "The force fields somewhat supported currently are:[PAR]",
425 "G53a5 GROMOS96 53a5 Forcefield (official distribution)[PAR]",
426 "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
427 "The corresponding data files can be found in the library directory",
428 "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
429 "information about file formats. By default, the force field selection",
430 "is interactive, but you can use the [TT]-ff[tt] option to specify",
431 "one of the short names above on the command line instead. In that",
432 "case [THISMODULE] just looks for the corresponding file.[PAR]",
434 const char *bugs
[] = {
435 "The atom type selection is primitive. Virtually no chemical knowledge is used",
436 "Periodic boundary conditions screw up the bonding",
437 "No improper dihedrals are generated",
438 "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."
441 t_params plist
[F_NRE
];
443 gpp_atomtype_t atype
;
448 char 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 */
455 gmx_bool bRTP
, bTOP
, bOPLS
;
459 gmx_output_env_t
*oenv
;
462 { efSTX
, "-f", "conf", ffREAD
},
463 { efTOP
, "-o", "out", ffOPTWR
},
464 { efRTP
, "-r", "out", ffOPTWR
}
466 #define NFILE asize(fnm)
467 static real 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 [REF].pdb[ref] 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
= gmx_strdup(molnm
);
541 /* Init parameter lists */
544 /* Read coordinates */
547 read_tps_conf(opt2fn("-f", NFILE
, fnm
), top
, &epbc
, &x
, NULL
, box
, FALSE
);
548 t_atoms
*atoms
= &top
->atoms
;
550 if (atoms
->pdbinfo
== NULL
)
552 snew(atoms
->pdbinfo
, natoms
);
555 sprintf(n2t
, "%s", ffdir
);
556 nm2t
= rd_nm2type(n2t
, &nnm
);
559 gmx_fatal(FARGS
, "No or incorrect atomname2type.n2t file found (looking for %s)",
564 printf("There are %d name to type translations in file %s\n", nnm
, n2t
);
568 dump_nm2type(debug
, nnm
, nm2t
);
570 printf("Generating bonds from distances...\n");
571 snew(nbonds
, atoms
->nr
);
572 mk_bonds(nnm
, nm2t
, atoms
, x
, &(plist
[F_BONDS
]), nbonds
, bPBC
, box
);
574 open_symtab(&symtab
);
575 atype
= set_atom_type(&symtab
, atoms
, &(plist
[F_BONDS
]), nbonds
, nnm
, nm2t
);
577 /* Make Angles and Dihedrals */
578 snew(excls
, atoms
->nr
);
579 printf("Generating angles and dihedrals from bonds...\n");
580 init_nnb(&nnb
, atoms
->nr
, 4);
581 gen_nnb(&nnb
, plist
);
582 print_nnb(&nnb
, "NNB");
583 gen_pad(&nnb
, atoms
, &rtp_header_settings
, plist
, excls
, NULL
, TRUE
);
588 plist
[F_LJ14
].nr
= 0;
591 "There are %4d %s dihedrals, %4d impropers, %4d angles\n"
592 " %4d pairs, %4d bonds and %4d atoms\n",
594 bOPLS
? "Ryckaert-Bellemans" : "proper",
595 plist
[F_IDIHS
].nr
, plist
[F_ANGLES
].nr
,
596 plist
[F_LJ14
].nr
, plist
[F_BONDS
].nr
, atoms
->nr
);
598 calc_angles_dihs(&plist
[F_ANGLES
], &plist
[F_PDIHS
], x
, bPBC
, box
);
600 set_force_const(plist
, kb
, kt
, kp
, bRound
, bParam
);
602 cgnr
= set_cgnr(atoms
, bUsePDBcharge
, &qtot
, &mtot
);
603 printf("Total charge is %g, total mass is %g\n", qtot
, mtot
);
612 fp
= ftp2FILE(efTOP
, NFILE
, fnm
, "w");
613 print_top_header(fp
, ftp2fn(efTOP
, NFILE
, fnm
), TRUE
, ffdir
, 1.0);
615 write_top(fp
, NULL
, mymol
.name
, atoms
, FALSE
, bts
, plist
, excls
, atype
,
616 cgnr
, rtp_header_settings
.nrexcl
);
617 print_top_mols(fp
, mymol
.name
, ffdir
, NULL
, 0, NULL
, 1, &mymol
);
623 print_rtp(ftp2fn(efRTP
, NFILE
, fnm
), "Generated by x2top",
624 atoms
, plist
, atype
, cgnr
);
629 dump_hybridization(debug
, atoms
, nbonds
);
631 close_symtab(&symtab
);
634 printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",
635 output_env_get_program_display_name(oenv
));
636 printf(" Please verify atomtypes and charges by comparison to other\n");
637 printf(" topologies.\n");