While trying to learn what happens when fixes affect energies and virials, I noticed that fix addforce doesn’t have per-atom energies enabled. Is this a conscious design choice, or simply “nobody’s gotten around to it”?
It seems to me that (1) fix addforce is commonly used, (2) without the per-atom energies, the documentation of compute pe/atom is incorrect, when it says that summing the compute pe/atom vector (e.g. with compute reduce) should give the same results as compute pe.
I’m happy to write in the per-atom energy tally if it would be welcome.
My understanding is that the potential energy depends only on the positions (and, for finite particles, orientations) of the particles. Therefore, only the kinetic energy should change upon using fix addforce.
Neither, it is more a historic choice that was changed when fixes were used to do things that they didn’t do before. You can track this down by looking at the history of features. Support for tallying energy was only added in commit c5aece5b435ff0402a711e88ea2e019d86d5a0b3 in 2016 in order to prepare LAMMPS for the addition of fix cmap (which is a 5-body force term that doesn’t fit to any of the force generating styles like pair, bond, angle, dihedral, or improper).
If you look through the source code, there is only one other fix that does support providing per-atom energies and that is fix external (which would use a callback mechanism to provide force and energy information from a coupled external code).
The challenge for fixes like fix addforce is that they are invoked in the Fix::post_force() step which has only “vflag” as an argument and not “eflag”. You can see in fix cmap how to work around it, but it is a bit ugly. The clean solution would be to change the API, but that would require changing a lot of code (for a feature used by so few fixes) and would break many external packages.
So before you go down the route of making changes, I suggest you create a feature request issue on GitHub and discuss with @sjplimp what would be the best way to proceed and whether there should be an API change (which would also probably require changes to the fix_modify behavior etc).
Force is the derivative of the potential energy, so you should be able to provide a potential energy if you have a force. In fact, fix addforce provides a global energy contribution (as indicated by energy_global_flag = 1;) so it could also provide a per-atom tally.
Kinetic energy is only changed during time integration or thermostatting which would normally be during Fix::initial_integrate() or Fix::final_integrate() or Fix::end_of_step() while fix addforce and similar are acting during the Fix::post_force() step which is after the re-computation of the forces after the position update from Fix::initial_integrate() and before the (second, half-)velocity update from Fix::final_integrate() in the velocity verlet scheme.
Ugliness is in the eye of the beholder – I did see the code in fix cmap by grepping for thermo_virial and the like, and the approach there seemed pretty intuitive to me. I liked that approach and it seems like it will be maintainable. If anything, it cost me way too much time to figure out that the code for ev_init calling ev_setup was in the header file fix.h.
In fact fix addforcealready calculates its eventual energy change per-atom and accumulates it. It should be pretty simple to shunt the result of each calculation into eatom as it goes along. But I will open a feature request as suggested and discuss more there.
For context, I’m looking into this for my work on the ELECTRODE package. Currently every energy change from fix electrode gets recorded through a force->pair->ev_tally call, which works fine in theory but makes it really difficult to disentangle and inspect all the different things we’re doing.