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.
7 * Copyright (c) 2018,2019,2020, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
45 #include "gromacs/commandline/pargs.h"
46 #include "gromacs/fileio/confio.h"
47 #include "gromacs/fileio/gmxfio.h"
48 #include "gromacs/gmxpreprocess/gen_ad.h"
49 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
50 #include "gromacs/gmxpreprocess/grompp_impl.h"
51 #include "gromacs/gmxpreprocess/nm2type.h"
52 #include "gromacs/gmxpreprocess/notset.h"
53 #include "gromacs/gmxpreprocess/pdb2top.h"
54 #include "gromacs/gmxpreprocess/toppush.h"
55 #include "gromacs/gmxpreprocess/toputil.h"
56 #include "gromacs/listed_forces/bonded.h"
57 #include "gromacs/math/units.h"
58 #include "gromacs/math/utilities.h"
59 #include "gromacs/math/vec.h"
60 #include "gromacs/math/vecdump.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/pbcutil/pbc.h"
63 #include "gromacs/topology/symtab.h"
64 #include "gromacs/topology/topology.h"
65 #include "gromacs/utility/arraysize.h"
66 #include "gromacs/utility/cstringutil.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/gmxassert.h"
69 #include "gromacs/utility/smalloc.h"
71 #include "hackblock.h"
73 #define MARGIN_FAC 1.1
75 static bool is_bond(int nnm
, t_nm2type nmt
[], char* ai
, char* aj
, real blen
)
79 for (i
= 0; (i
< nnm
); i
++)
81 for (j
= 0; (j
< nmt
[i
].nbonds
); j
++)
83 if ((((gmx::equalCaseInsensitive(ai
, nmt
[i
].elem
, 1))
84 && (gmx::equalCaseInsensitive(aj
, nmt
[i
].bond
[j
], 1)))
85 || ((gmx::equalCaseInsensitive(ai
, nmt
[i
].bond
[j
], 1))
86 && (gmx::equalCaseInsensitive(aj
, nmt
[i
].elem
, 1))))
87 && (fabs(blen
- nmt
[i
].blen
[j
]) <= 0.1 * nmt
[i
].blen
[j
]))
96 static void mk_bonds(int nnm
,
100 InteractionsOfType
* bond
,
110 std::array
<real
, MAXFORCEPARAM
> forceParam
= { 0.0 };
113 set_pbc(&pbc
, PbcType::Unset
, box
);
115 for (i
= 0; (i
< atoms
->nr
); i
++)
119 fprintf(stderr
, "\ratom %d", i
);
122 for (j
= i
+ 1; (j
< atoms
->nr
); j
++)
126 pbc_dx(&pbc
, x
[i
], x
[j
], dx
);
130 rvec_sub(x
[i
], x
[j
], dx
);
134 if (is_bond(nnm
, nmt
, *atoms
->atomname
[i
], *atoms
->atomname
[j
], std::sqrt(dx2
)))
136 forceParam
[0] = std::sqrt(dx2
);
137 std::vector
<int> atoms
= { i
, j
};
138 add_param_to_list(bond
, InteractionOfType(atoms
, forceParam
));
144 fprintf(stderr
, "\ratom %d\n", i
);
148 static int* set_cgnr(t_atoms
* atoms
, bool bUsePDBcharge
, real
* qtot
, real
* mtot
)
155 snew(cgnr
, atoms
->nr
);
156 for (i
= 0; (i
< atoms
->nr
); i
++)
158 if (atoms
->pdbinfo
&& bUsePDBcharge
)
160 atoms
->atom
[i
].q
= atoms
->pdbinfo
[i
].bfac
;
162 qt
+= atoms
->atom
[i
].q
;
163 *qtot
+= atoms
->atom
[i
].q
;
164 *mtot
+= atoms
->atom
[i
].m
;
175 static void set_atom_type(PreprocessingAtomTypes
* atypes
,
178 InteractionsOfType
* bonds
,
185 snew(atoms
->atomtype
, atoms
->nr
);
186 nresolved
= nm2type(nnm
, nm2t
, tab
, atoms
, atypes
, nbonds
, bonds
);
187 if (nresolved
!= atoms
->nr
)
189 gmx_fatal(FARGS
, "Could only find a forcefield type for %d out of %d atoms", nresolved
, atoms
->nr
);
192 fprintf(stderr
, "There are %zu different atom types in your sample\n", atypes
->size());
195 static void lo_set_force_const(InteractionsOfType
* plist
, real c
[], int nrfp
, bool bRound
, bool bDih
, bool bParam
)
200 for (auto& param
: plist
->interactionTypes
)
204 for (int j
= 0; j
< nrfp
; j
++)
213 sprintf(buf
, "%.2e", param
.c0());
214 sscanf(buf
, "%lf", &cc
);
224 c
[0] = (static_cast<int>(c
[0] + 3600)) % 360;
229 /* To put the minimum at the current angle rather than the maximum */
233 GMX_ASSERT(nrfp
<= MAXFORCEPARAM
/ 2, "Only 6 parameters may be used for an interaction");
234 std::array
<real
, MAXFORCEPARAM
> forceParam
;
235 for (int j
= 0; (j
< nrfp
); j
++)
237 forceParam
[j
] = c
[j
];
238 forceParam
[nrfp
+ j
] = c
[j
];
240 param
= InteractionOfType(param
.atoms(), forceParam
);
244 static void set_force_const(gmx::ArrayRef
<InteractionsOfType
> plist
, real kb
, real kt
, real kp
, bool bRound
, bool bParam
)
246 real c
[MAXFORCEPARAM
];
250 lo_set_force_const(&plist
[F_BONDS
], c
, 2, bRound
, FALSE
, bParam
);
252 lo_set_force_const(&plist
[F_ANGLES
], c
, 2, bRound
, FALSE
, bParam
);
255 lo_set_force_const(&plist
[F_PDIHS
], c
, 3, bRound
, TRUE
, bParam
);
258 static void calc_angles_dihs(InteractionsOfType
* ang
, InteractionsOfType
* dih
, const rvec x
[], bool bPBC
, matrix box
)
261 rvec r_ij
, r_kj
, r_kl
, m
, n
;
267 set_pbc(&pbc
, PbcType::Xyz
, box
);
269 for (auto& angle
: ang
->interactionTypes
)
275 * bond_angle(x
[ai
], x
[aj
], x
[ak
], bPBC
? &pbc
: nullptr, r_ij
, r_kj
, &costh
, &t1
, &t2
);
276 angle
.setForceParameter(0, th
);
278 for (auto dihedral
: dih
->interactionTypes
)
280 int ai
= dihedral
.ai();
281 int aj
= dihedral
.aj();
282 int ak
= dihedral
.ak();
283 int al
= dihedral
.al();
285 * dih_angle(x
[ai
], x
[aj
], x
[ak
], x
[al
], bPBC
? &pbc
: nullptr, r_ij
, r_kj
, r_kl
,
286 m
, n
, &t1
, &t2
, &t3
);
287 dihedral
.setForceParameter(0, ph
);
291 static void dump_hybridization(FILE* fp
, t_atoms
* atoms
, int nbonds
[])
293 for (int i
= 0; (i
< atoms
->nr
); i
++)
295 fprintf(fp
, "Atom %5s has %1d bonds\n", *atoms
->atomname
[i
], nbonds
[i
]);
300 print_pl(FILE* fp
, gmx::ArrayRef
<const InteractionsOfType
> plist
, int ftp
, const char* name
, char*** atomname
)
302 if (!plist
[ftp
].interactionTypes
.empty())
305 fprintf(fp
, "[ %s ]\n", name
);
306 int nrfp
= interaction_function
[ftp
].nrfpA
;
307 for (const auto& param
: plist
[ftp
].interactionTypes
)
309 gmx::ArrayRef
<const int> atoms
= param
.atoms();
310 gmx::ArrayRef
<const real
> forceParam
= param
.forceParam();
311 for (const auto& atom
: atoms
)
313 fprintf(fp
, " %5s", *atomname
[atom
]);
315 for (int j
= 0; (j
< nrfp
); j
++)
317 if (forceParam
[j
] != NOTSET
)
319 fprintf(fp
, " %10.3e", forceParam
[j
]);
327 static void print_rtp(const char* filenm
,
330 gmx::ArrayRef
<const InteractionsOfType
> plist
,
331 PreprocessingAtomTypes
* atypes
,
338 fp
= gmx_fio_fopen(filenm
, "w");
339 fprintf(fp
, "; %s\n", title
);
341 fprintf(fp
, "[ %s ]\n", *atoms
->resinfo
[0].name
);
343 fprintf(fp
, "[ atoms ]\n");
344 for (i
= 0; (i
< atoms
->nr
); i
++)
346 tp
= atoms
->atom
[i
].type
;
347 if ((tpnm
= atypes
->atomNameFromAtomType(tp
)) == nullptr)
349 gmx_fatal(FARGS
, "tp = %d, i = %d in print_rtp", tp
, i
);
351 fprintf(fp
, "%-8s %12s %8.4f %5d\n", *atoms
->atomname
[i
], tpnm
, atoms
->atom
[i
].q
, cgnr
[i
]);
353 print_pl(fp
, plist
, F_BONDS
, "bonds", atoms
->atomname
);
354 print_pl(fp
, plist
, F_ANGLES
, "angles", atoms
->atomname
);
355 print_pl(fp
, plist
, F_PDIHS
, "dihedrals", atoms
->atomname
);
356 print_pl(fp
, plist
, F_IDIHS
, "impropers", atoms
->atomname
);
361 int gmx_x2top(int argc
, char* argv
[])
363 const char* desc
[] = {
364 "[THISMODULE] generates a primitive topology from a coordinate file.",
365 "The program assumes all hydrogens are present when defining",
366 "the hybridization from the atom name and the number of bonds.",
367 "The program can also make an [REF].rtp[ref] entry, which you can then add",
368 "to the [REF].rtp[ref] database.[PAR]",
369 "When [TT]-param[tt] is set, equilibrium distances and angles",
370 "and force constants will be printed in the topology for all",
371 "interactions. The equilibrium distances and angles are taken",
372 "from the input coordinates, the force constant are set with",
373 "command line options.",
374 "The force fields somewhat supported currently are:[PAR]",
375 "G53a5 GROMOS96 53a5 Forcefield (official distribution)[PAR]",
376 "oplsaa OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)[PAR]",
377 "The corresponding data files can be found in the library directory",
378 "with name [TT]atomname2type.n2t[tt]. Check Chapter 5 of the manual for more",
379 "information about file formats. By default, the force field selection",
380 "is interactive, but you can use the [TT]-ff[tt] option to specify",
381 "one of the short names above on the command line instead. In that",
382 "case [THISMODULE] just looks for the corresponding file.[PAR]",
384 const char* bugs
[] = {
385 "The atom type selection is primitive. Virtually no chemical knowledge is used",
386 "Periodic boundary conditions screw up the bonding", "No improper dihedrals are generated",
387 "The atoms to atomtype translation table is incomplete ([TT]atomname2type.n2t[tt] file in "
388 "the data directory). Please extend it and send the results back to the GROMACS crew."
391 std::array
<InteractionsOfType
, F_NRE
> plist
;
396 char forcefield
[32], ffdir
[STRLEN
];
397 rvec
* x
; /* coordinates? */
399 int bts
[] = { 1, 1, 1, 2 };
400 matrix box
; /* box length matrix */
401 int natoms
; /* number of atoms in one molecule */
403 bool bRTP
, bTOP
, bOPLS
;
407 gmx_output_env_t
* oenv
;
409 t_filenm fnm
[] = { { efSTX
, "-f", "conf", ffREAD
},
410 { efTOP
, "-o", "out", ffOPTWR
},
411 { efRTP
, "-r", "out", ffOPTWR
} };
412 #define NFILE asize(fnm)
413 real kb
= 4e5
, kt
= 400, kp
= 5;
414 PreprocessResidue rtp_header_settings
;
415 bool bRemoveDihedralIfWithImproper
= FALSE
;
416 bool bGenerateHH14Interactions
= TRUE
;
417 bool bKeepAllGeneratedDihedrals
= FALSE
;
419 bool bParam
= TRUE
, bRound
= TRUE
;
420 bool bPairs
= TRUE
, bPBC
= TRUE
;
421 bool bUsePDBcharge
= FALSE
, bVerbose
= FALSE
;
422 const char* molnm
= "ICE";
423 const char* ff
= "oplsaa";
429 "Force field for your simulation. Type \"select\" for interactive selection." },
430 { "-v", FALSE
, etBOOL
, { &bVerbose
}, "Generate verbose output in the top file." },
431 { "-nexcl", FALSE
, etINT
, { &nrexcl
}, "Number of exclusions" },
435 { &bGenerateHH14Interactions
},
436 "Use 3rd neighbour interactions for hydrogen atoms" },
440 { &bKeepAllGeneratedDihedrals
},
441 "Generate all proper dihedrals" },
445 { &bRemoveDihedralIfWithImproper
},
446 "Remove dihedrals on the same bond as an improper" },
447 { "-pairs", FALSE
, etBOOL
, { &bPairs
}, "Output 1-4 interactions (pairs) in topology file" },
448 { "-name", FALSE
, etSTR
, { &molnm
}, "Name of your molecule" },
449 { "-pbc", FALSE
, etBOOL
, { &bPBC
}, "Use periodic boundary conditions." },
454 "Use the B-factor supplied in a [REF].pdb[ref] file for the atomic charges" },
455 { "-param", FALSE
, etBOOL
, { &bParam
}, "Print parameters in the output" },
456 { "-round", FALSE
, etBOOL
, { &bRound
}, "Round off measured values" },
457 { "-kb", FALSE
, etREAL
, { &kb
}, "Bonded force constant (kJ/mol/nm^2)" },
458 { "-kt", FALSE
, etREAL
, { &kt
}, "Angle force constant (kJ/mol/rad^2)" },
459 { "-kp", FALSE
, etREAL
, { &kp
}, "Dihedral angle force constant (kJ/mol/rad^2)" }
462 if (!parse_common_args(&argc
, argv
, 0, NFILE
, fnm
, asize(pa
), pa
, asize(desc
), desc
,
463 asize(bugs
), bugs
, &oenv
))
467 bRTP
= opt2bSet("-r", NFILE
, fnm
);
468 bTOP
= opt2bSet("-o", NFILE
, fnm
);
469 /* C89 requirements mean that these struct members cannot be used in
470 * the declaration of pa. So some temporary variables are needed. */
471 rtp_header_settings
.bRemoveDihedralIfWithImproper
= bRemoveDihedralIfWithImproper
;
472 rtp_header_settings
.bGenerateHH14Interactions
= bGenerateHH14Interactions
;
473 rtp_header_settings
.bKeepAllGeneratedDihedrals
= bKeepAllGeneratedDihedrals
;
474 rtp_header_settings
.nrexcl
= nrexcl
;
478 gmx_fatal(FARGS
, "Specify at least one output file");
481 /* Force field selection, interactive or direct */
482 choose_ff(strcmp(ff
, "select") == 0 ? nullptr : ff
, forcefield
, sizeof(forcefield
), ffdir
,
485 bOPLS
= (strcmp(forcefield
, "oplsaa") == 0);
488 mymol
.name
= gmx_strdup(molnm
);
491 /* Read coordinates */
494 read_tps_conf(opt2fn("-f", NFILE
, fnm
), top
, &pbcType
, &x
, nullptr, box
, FALSE
);
495 t_atoms
* atoms
= &top
->atoms
;
497 if (atoms
->pdbinfo
== nullptr)
499 snew(atoms
->pdbinfo
, natoms
);
502 sprintf(n2t
, "%s", ffdir
);
503 nm2t
= rd_nm2type(n2t
, &nnm
);
506 gmx_fatal(FARGS
, "No or incorrect atomname2type.n2t file found (looking for %s)", n2t
);
510 printf("There are %d name to type translations in file %s\n", nnm
, n2t
);
514 dump_nm2type(debug
, nnm
, nm2t
);
516 printf("Generating bonds from distances...\n");
517 snew(nbonds
, atoms
->nr
);
518 mk_bonds(nnm
, nm2t
, atoms
, x
, &(plist
[F_BONDS
]), nbonds
, bPBC
, box
);
520 open_symtab(&symtab
);
521 PreprocessingAtomTypes atypes
;
522 set_atom_type(&atypes
, &symtab
, atoms
, &(plist
[F_BONDS
]), nbonds
, nnm
, nm2t
);
524 /* Make Angles and Dihedrals */
525 snew(excls
, atoms
->nr
);
526 printf("Generating angles and dihedrals from bonds...\n");
527 gen_pad(atoms
, gmx::arrayRefFromArray(&rtp_header_settings
, 1), plist
, excls
, {}, TRUE
);
531 plist
[F_LJ14
].interactionTypes
.clear();
534 "There are %4zu %s dihedrals, %4zu impropers, %4zu angles\n"
535 " %4zu pairs, %4zu bonds and %4d atoms\n",
536 plist
[F_PDIHS
].size(), bOPLS
? "Ryckaert-Bellemans" : "proper", plist
[F_IDIHS
].size(),
537 plist
[F_ANGLES
].size(), plist
[F_LJ14
].size(), plist
[F_BONDS
].size(), atoms
->nr
);
539 calc_angles_dihs(&plist
[F_ANGLES
], &plist
[F_PDIHS
], x
, bPBC
, box
);
541 set_force_const(plist
, kb
, kt
, kp
, bRound
, bParam
);
543 cgnr
= set_cgnr(atoms
, bUsePDBcharge
, &qtot
, &mtot
);
544 printf("Total charge is %g, total mass is %g\n", qtot
, mtot
);
553 fp
= ftp2FILE(efTOP
, NFILE
, fnm
, "w");
554 print_top_header(fp
, ftp2fn(efTOP
, NFILE
, fnm
), TRUE
, ffdir
, 1.0);
556 write_top(fp
, nullptr, mymol
.name
.c_str(), atoms
, FALSE
, bts
, plist
, excls
, &atypes
, cgnr
,
557 rtp_header_settings
.nrexcl
);
558 print_top_mols(fp
, mymol
.name
.c_str(), ffdir
, nullptr, {}, gmx::arrayRefFromArray(&mymol
, 1));
564 print_rtp(ftp2fn(efRTP
, NFILE
, fnm
), "Generated by x2top", atoms
, plist
, &atypes
, cgnr
);
569 dump_hybridization(debug
, atoms
, nbonds
);
571 close_symtab(&symtab
);
573 printf("\nWARNING: topologies generated by %s can not be trusted at face value.\n",
574 output_env_get_program_display_name(oenv
));
575 printf(" Please verify atomtypes and charges by comparison to other\n");
576 printf(" topologies.\n");