Analytical vs. numerical derivatives command


Is there a command in LAMMPS to automatically perform numerical differentiation, in order to test analytical derivatives vs. numerical derivatives? I couldn't find anything in the manual.


If you want to take a time derivative of a per-atom property I think you can use fix state/store and an atom-style variable to achieve this. I am not sure if there is an equivalent fix for global quantities.

Hi Stefan,

Thanks. I should have precised though that I’m interested in derivatives wrt atomic position - mostly to test numerical stability of potentials.


Right. In that case it is more difficult I guess. I have no clue apart from writing your own fix or compute as conceptually it is not that hard, but maybe someone else has a quicker solution.

Well right now I have a script which loops over atoms, x, y and z coordinates, and computing numerical derivatives from that can easily be scripted, but I was wondering if there was a LAMMPS command to directly achieve this (which is the case in a lot of MD packages I think). It would of course be more convenient!


Interesting idea. If you want the derivative wrt each atom,
how would such a command work, other than doing a 3N

length loop to displace each atom in each direction?

Which for a million atoms and a many-body or long-range
potential, would take a long time.


Yes, that is the issue - in practice, for bulk sustems, I’m just picking up a certain number of atoms randomly to do perturbation on and check the derivatives. But for small systems like oligomers for instance, the 3N loop is totally doable. This makes it probably difficult to build a general command for it, and input scripting is well suited for this kind of task.


As a midway solution you could use the library interface LAMMPS has to leverage force computations to it. I personally don’t like it because it is “too C’ish” for me so I kludged some C+±like wrappers around it, which I can share if you want. Then it really is as easy as looping over coordinates and recalculating. This way you can test your implementations in a LAMMPS-like pair_style first and then, if it suits your needs, it is already inclusion-ready.