Question regarding the harmonic/shift bond coefficients used in the lammps code

Dear lammps users:

I want to reproduce the results of a paper(doi:Redirecting) which utilized lammps for molecular simulation of a spin crossover material.

The potential expression given in this paper uses the harmonic/shift bond style


, and I need to convert the force constant given in the paper to the corresponding coefficients used in the lammps code. From the lammps documentation(release 19Nov2024) section 5.11.3, it states that the harmonic/shift bond potential expression is

, then it states that the spring constant is
3
, which confuses me. According to my understanding, shouldn’t the spring constant be “k = 2Umin/(r0 − rc)^2” rather than “k = Umin/[2(r0 − rc)2]”? Is it a typo or something else?

Thank you very much for your time and your help are much appreciated.

You can use the bond_write command to confirm the computed potential energy and forces used for a given set of parameters.

Thank you very much for your reply.

Unfortunately, I don’t have something like a bond table to compare the bond_write output data with, because the paper only gives the force constant and final simulation results without any supplement materials or any other kind of related data.

The reason I want to know if the k expression in the lammps documentation was a typo is that I’m on the way to figure out the cause of the discrepancy between my simulation results and the paper’s results, one of which may come from the wrong coefficients I used in the lammps code. For example, the paper gives this table below


, and I need to convert the force constants of this table to the harmonic/shift bond coefficients(Umin,r0,rc) used in the lammps code. But I don’t know what the force constants given in the table refer to. Is it Kb of the potential energy expression given in the paper? Is it k of the expression “k = Umin/[2(r0 − rc)^2]” given in the lammps documentation section 5.11.3? Is it K of the expression “K=Umin/(r0 − rc)^2”? Or is it K of the expression “K = 2Umin/(r0 − rc)^2”?

The corresponding authors haven’t replied to my asking-for-help email yet.

Thank you very much for your time.

That is not why I suggested this, but rather as a means to determine whether the parameters you use yield the desired functional form.

Given that the paper seems to provide only information about k and r0 and also since the publicly visible part of the publication only mentions using a regular harmonic potential, the inconsistency is more likely the mention of harmonic/shift. With the provided information, you can only use bond style harmonic. Apart from a constant factor in the bond energy, the result should be the same.

Since I cannot access the entire paper at the moment, could it be that the authors used the harmonic shift formulation to fit the potentials to their experimental data but then use the regular harmonic bond style with the parameters from the table for their simulations?