Potential energy computations in lammps

Hello.

I’m simulating a Brownian dynamics system of semiflexible, bead-spring polymer networks in LAMMPS using fix langevin and fix nve/limit, applying displacements to groups of atoms via fix move, and then computing the work done on the network in post processing. I’m noticing a large discrepancy in some of the values.

When I compute the work numerically as a sum of the force against each atom’s trajectory over time, I get values on the the order of 10^7. Looking at the kinetic and potential energies in the thermo output, they both remain in the 10^1 range for the entire simulation. Note that I’m using nondimensionalized LJ units in my simulations.

This is confusing (and if I’m missing anything here I’d love to know), but I think it could be due to the fact that ke and pe are computed continuously in LAMMPS while my work computation takes the dot products of an instantaneous force and a displacement between snapshots that are 1000s of timesteps apart.

What’s more concerning is looking at the bond data. For my purposes, I’m interested in the potential energy changes due to different types of bonds. Using compute bond/local and the engpot keyword, I can look at the potential energy in each bond at every dumped timestep, and some individually have more potential energy than the potential energy of the system at the same time step, not even considering the sum of potential energies over all bonds. What’s going on here? Shouldn’t the potential energy of the system be strictly greater than the sum of the potential energy of all the bonds given that I have pair potentials and angle potentials in my simulation?

Hello.

[...]

This is confusing (and if I'm missing anything here I'd love to know), but I think it could be due to the fact that ke and pe are computed continuously in LAMMPS while my work computation takes the dot products of an instantaneous force and a displacement between snapshots that are 1000s of timesteps apart.

sorry, but kinetic energy and potential energy are both state
functions, i.e. they are fully defined by their state, kinetic energy
by mass and velocity, and potential energy by position and potential
function. their history does not matter (except for models where the
potential function itself is history dependent). the work you talk
about on the other hand is not, as it depends on the path and is the
integral over F dot r dt along the exact path of each particle.

What's more concerning is looking at the bond data. For my purposes, I'm interested in the potential energy changes due to different types of bonds. Using compute bond/local and the engpot keyword, I can look at the potential energy in each bond at every dumped timestep, and some individually have more potential energy than the potential energy of the system at the same time step, not even considering the sum of potential energies over all bonds. What's going on here? Shouldn't the potential energy of the system be strictly greater than the sum of the potential energy of all the bonds given that I have pair potentials and angle potentials in my simulation?

why should that be? that entirely depends on the definition of the
potential functions and their absolute value is completely arbitrary.
if you use a lennard-jones potential, for example, the potential
energy can be < 0.

axel.