Simulating a small molecule with custom topology

GROMACS version: 2018
GROMACS modification: No

I am trying to simulate Cobalt Hexammine in Gromacs. I have found a source that lists all of the bonded parameters (bonds, angles) as well as the nonbonded parameters (sigma, epsilon, and charge).
I created an .itp file that lists the atoms with their charge and mass (in the [atoms] section) as well as the bonds and angles (in the [bonds] and [angles] sections, respectively). (There are no dihedral constraints). However, I have three questions:

  1. Do I still need to modify my bonded force field (“ffbonded.itp”)? Or is the topology enough?
  2. Where do I include my non-bonded interactions (sigma and epsilon)?
  3. Do I need to include anything else except the CoHex pdb and this topology file?
    I’ve read through the tutorials and manual, but I’m still a bit confused. Any help or points to resources would be greatly appreciated.

What force field are you trying to use? Most of the force fields implemented in GROMACS have the atomtypes for the amine and chloride defined, but none to my knowledge have the atomtype for cobalt defined (I could be wrong on this).

Depending on what force field you use, the topology in which you generate will differ.

I’m using Charmm at the moment (I may use Amber as well in the future). However, the study that I am trying to replicate uses a very specific set of parameters for all atoms in the Cobalt Hexamine so that they reproduce observed experimental behavior.

If your simulating cobalt (III) hexamine only and using a specific set of parameters from the literature, it is probably best to generate the topology lists you need yourself. Since you are using CHARMM, I would suggest your topology look like

#include “charmm27.ff/forcefield.itp”

[ atomtypes ]
;name at.num mass charge ptype sigma (nm) epsilon (kJ/mol)
(you can complete the rest on each of these lines)

[ bondtypes ]
; atomtype 1 atomtype2 func
(If func 1 is used, then provide b0 (nm) and kb (kJ/mol) respectively on each line)

[ angletypes ]
; atomtype 1 atomtype2 atomtype3 func
(If func 1 is used, then provide theta0 (deg) and k_theta (kJ/(mol rad^2)) respectively on each line)

#include “charmm27.ff/spce.itp”
(If you want to include some water for example)

#include “cobalt_hexammine.itp” (which contains the following)
[ moleculetype ]
; resname nrexcl
“residue name of your choice” (I’ll use CHN) 3

[ atoms ]
; nr type resnr residu atom cgnr charge mass
1 Co 1 CHN Co 1 (fill) (fill)
2 N 1 CHN N1 1 (fill) (fill)
3 H 1 CHN H1 1 (fill) (fill)
(repeat this for each atom which should provide a total of 25 entries)

[ bonds ]
; ai aj func
1 2 1
etc. for each bond

[ angles ]
; ai aj ak func
1 2 3 1
etc. for each angle


[ system ]
“Name whatever you want”

[ molecules ]

So to answer your 3 questions more directly

  1. You in theory could add the 3 atomtypes to the ffnonbonded.itp file and the 2 bondtypes and 3 angletypes to the ffbonded.itp file. I just add the lists directly to the topology file because its easier to do so and unlike previous GROMACS versions, you can have multiple parameter-level directives present and order does not matter with the exception that atomtypes need to be defined before they can be referred to in a parameter-level directive.
  2. sigma and epsilon go in the [ atomtypes ] section. You can add another [ atomtypes ] list yourself or modify the one in ffnonbonded.itp.
  3. Three files are needed - topology (.top), coordinate (.gro/.pdb), and MD parameter (.mdp) files to start a simulation. The coordinate file must have residues listed in the same order as in [ molecules ] in topology.

I would suggest reading the following manual sections here on how to generate topology or understand how it works.

This is exceedingly helpful. Thank you.

Last question, I promise:
When I define tha “atom types” section, does this “overwrite” the definition for N and H for the other parts of my system? I other words, if I am simulating CoHex and DNA, will the N and H in the DNA still get the correct CHARMM parameters?

Thanks again for all of your help.

Yes, there will be overwriting. If the 2nd [ atomtypes ] directive above contains an atomtype present in the 1st [ atomtypes ] directive that is contained within ffnonbonded.itp, it will overwrite it. Therefore, I suggest naming the atomtypes in the 2nd [ atomtypes ] list differently and uniquely to prevent such overwriting from occurring as DNA uses the “H” atomtype if I recall correctly.

  • “Co” can stay as “Co” since there is no cobalt atom defined in [ atomtypes ] of ffnonbonded.itp.
  • “H” can be named as “H_Co” for example to prevent overwriting of “H” information
  • “N” can be named as “N_Co” for example to prevent overwriting of “N” information

In any case, by naming these atomtypes differently, you prevent the parameters for CoHex from clashing with that of DNA.

Perfect. Thank you, again. This has been extremely helpful.