Hi,
I am trying to simulate a material in Gromacs. however, I am using parameters used in lammps.
There is a confusing part which is related to the bonds.
based on what I read from the article, the bond between Ti and O is at the length of 2.1 A (0.21 nm) and the force constant is 400 Kcal/mol (1673.6 kj/mol).
Now the question I have is that what type of bond it is.
I searched here as well but so far haven’t been successful to understand it.
I would be very appreciative if someone could help me understand this.
To me, it looks like a standard type-1 harmonic bond, emulating a fixed covalent bond between two atoms.
Note that the unit of a force constant is not the unit of energy, but the unit of force per distance, or energy per distance², so you need to factor in the distance unit conversion factor (1nm/1A = 10) squared (=100).
I was reading this, https://www.ks.uiuc.edu/Research/vmd/current/ug/node237.html, last night.
it looks like the force constant reported is scaled one. Meaning if the length of the bond is 2.1 A then the spring constant used is 400/(2.1^2) kcal/mol/A^2.
But this is the documentation of the Colvars module, which is used for defining and biasing collective variables. The scaling coefficient is equal to 1 by default, and is only useful for multidimensional restraints with a unique force constant.
If you’re transferring parameters from LAMMPS, you might want to check their definition of the harmonic bond, in particular whether they put the ½ factor in front if the k(x0-x)² - for example, AFAIK Gromacs and NAMD do but Amber doesn’t. Otherwise it’s usually just keeping track of units.
Thanks,
so you mean, I just have to divide 400 kcal/mol to 1 A^2 so that K = 400 kcal/mol/A^2.
In case of lammps, regarding the documentation, bond_style harmonic command — LAMMPS documentation, 1/2 is included in the K. So that k = 2 * 400 kcal/mol/A^2.
There’s no need to divide anything in principle, this is a formula not related to bonded terms but arbitrary harmonic restraints. But still, that would not change anything.
So yes, seems you need to do 2 * 400 * 4.184 * 100 to get the force constant in Gromacs units.
It’s also good to double-check that with an example parameter set that was officially published for both Gromacs and LAMMPS.
@ahosseini9473 For context, the Colvars definition of a harmonic potential includes the scaling parameter width because it is defined on collective variables of very different kinds, which may have very different units.
In a force field, the harmonic potential is always defined on a distance, and so the unit is always known: no need to scale anything there.
BTW the VMD doc page you linked is really old (VMD 1.9.3 is 8 years old). For other forum readers, I’m linking here the relevant Colvars doc page for GROMACS.
Giacomo