Potential term support

[ahy3nz]

Different engine formats and external packages have slightly different functions, which will affect the resultant parameters for I/O:

Bond potentials

engine/package	functional form	comments
internal	V = (1/2) * k * (r-r_eq)**2	Default form in BondType
gromacs	V = (1/2) * k * (r-r0)**2	gmx bondtype 1
lammps	V = k * (r-r0)**2	bond_style harmonic
hoomd	V = (1/2) * k * (r-r0)**2	hoomd.md.bond.harmonic
parmed	V = k * (r-r0)**2	This code
openmm xml	V = (1/2) * k * (r-r0)**2	OpenMM 7.0 Theory guide

Angle potentials

• Internally, k in units of kj/deg**2 and theta_eq in units of deg

engine/package	functional form	comments
internal	V = (1/2) * k * (theta-theta_eq)**2	Default form in AngleType
gromacs	V = (1/2) * k * (theta -theta0)**2	gmx angletype 1
parmed	V = k * (theta-theta0)**2	This code

Dihedral potentials

Potential Templates

name	functional form	parameters	comments
LennardJonesPotential	V = 4 * epsilon * ((sigma/r)**12 - (sigma/r)**6	epsilon, sigma
MiePotential	V = (n/(n-m)) * (n/m)**(m/(n-m)) * epsilon * ((sigma/r)**n - (sigma/r)**m)	n, m, epsilon, sigma
BuckinghamPotential	V = a * exp(-b * r) - c * r**-6	a, b, c
HarmonicBondPotential	V = 0.5 * k * (r-r_eq)**2	k, r_eq
HarmonicAnglePotential	V = 0.5 * k * (theta - theta_eq)**2	k, theta_eq
HarmonicTorsionPotential	V = 0.5 * k * (phi - phi_eq)**2	k, phi_eq	phi_cis = 0
PeriodicTorsionPotential	V = k * (1 + cos(n * phi - phi_eq))**2	k, n, phi_eq	phi_cis = 0
OPLSTorsionPotential	V = k0 + 0.5 * k1 * (1 + cos(phi)) + 0.5 * k2 * (1 - cos(2 * phi)) + 0.5 * k3 * (1 + cos(3 * phi)) + 0.5 * k4 * (1 - cos(4 * phi))	k0, k1, k2, k3, k4	phi_cis = 0
RyckaertBellemansTorsionPotential	V = c0 * cos(phi)**0 + c1 * cos(phi)**1 + c2 * cos(phi)**2 + c3 * cos(phi)**3 + c4 * cos(phi)**4 + c5 * cos(phi)**5	c0, c1, c2, c3, c4, c5	phi_trans = 0

[mattwthompson]

I wonder if we can get sympy to intelligently compare functions that are similar but not quite exact, like the case of harmonic bond potentials sometimes having a 1/2 coefficient and sometimes having it wrapped into the force constant.

[ahy3nz]

>>> func1 = sympy.sympify('0.5 * k * (r-r0)**2')
>>> func2 = sympy.sympify('k * (r-r0)**2')
>>> func1
0.5*k*(r - r0)**2
>>> func2
k*(r - r0)**2
>>> small_k = 500
>>> big_k = 1000
>>> func1.subs({'k': small_k})
250.0*(r - r0)**2
>>> func1.subs({'k': small_k})
250.0*(r - r0)**2
>>> func1.subs({'k': big_k})
500.0*(r - r0)**2
>>> func2.subs({'k': small_k})
500*(r - r0)**2
>>> subs_1 = func1.subs({'k': big_k})
>>> subs_2 = func2.subs({'k': small_k})
>>> subs_1
500.0*(r - r0)**2
>>> subs_2
500*(r - r0)**2
>>> subs_1 == subs_2
True
>>> func1 == func2
False

It looks like if you substitute 'k' into the sympy.functions, it looks like equality can be compared. but if you compare the non-substituted sympy.functions, you won't get equality