include force function in angle table example to show the need for correct unit conversion (force is energy per radian squared)
This commit is contained in:
@ -5,8 +5,17 @@ from tabulate import AngleTabulate
|
||||
################################################################################
|
||||
import math
|
||||
|
||||
def harmonic_force(theta):
|
||||
k = 50.0
|
||||
thetazero = 120.0
|
||||
# the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians
|
||||
deg2rad = math.pi / 180.0
|
||||
t = (theta - thetazero) * deg2rad
|
||||
f = -2.0 * k * t * deg2rad
|
||||
return f
|
||||
|
||||
def harmonic_energy(theta):
|
||||
k = 100.0
|
||||
k = 50.0
|
||||
thetazero = 120.0
|
||||
# the force constant in LAMMPS is in energy per radians^2 so convert from degrees to radians
|
||||
deg2rad = math.pi / 180.0
|
||||
@ -17,5 +26,5 @@ def harmonic_energy(theta):
|
||||
################################################################################
|
||||
|
||||
if __name__ == "__main__":
|
||||
atable = AngleTabulate(harmonic_energy, units='real')
|
||||
atable = AngleTabulate(harmonic_energy, harmonic_force, units='real')
|
||||
atable.run('HARM')
|
||||
|
||||
Reference in New Issue
Block a user