Energy Conservation in fix propel/self command

Hello!

I am simulating dense self-propelled molecules in two dimensions using LAMMPS (29 Aug 2024 version). The input scripts are attached. However, I have observed that the total energy of the system exhibits a drift over long simulations when the self-propelled force is applied. In contrast, no energy drift occurs when the self-propelled force is absent.

Could this drift be due to the self-propelled force being an external force to the system, with the potential energy not accounting for its contribution? If so, it seems the minimize command might not consider the contribution of the active force during energy minimization. I have seen that there is no fix_modify energy option available for the self-propelled force, unlike the fix addforce command.

Could you please help to resolve this energy drift and suggest any solutions?

active_dipole.in.in (946 Bytes)
Dumbbell.txt (142 Bytes)

It is only sensible to define a potential energy when a force is the gradient of that energy. Indeed, that gradient relationship is how one proves that energy is conserved (which is why a force must be the gradient of some energy function to be called a “conservative” force).

Thus fix propel/self does not have a conserved energy associated with it. This also makes it impossible to “count” such a contribution during minimisation.

1 Like

There is nothing to resolve. If you drive a car (also an object that propels itself), you have to put gas in the tank (or electricity into the battery) so that you can drive. That energy is released into the system. The same happens here, only that your cars have an infinitely large tank.

As for the minimization. That is a computation without dynamics. In my example comparing your self-propelled particles with cars, that would mean the cars are turned off, so there is nothing to add.

See @srtee’s response for an explanation using proper physics arguments instead of common sense.

Dear @srtee and Dear @akohlmey

Thank you for your helpful comments. In the active system, the virial (pressure) is defined as f_P∑<ei.ri>/(dV). For this, the default pressure calculation in lammps includes the use of fix_modify virial to account for the contribution from active stresses. Similarly, potential energy can be defined as f_P∑<ei.ri>, which represents its contribution to the potential energy. Am I correct in my understanding?