Calculate force on atoms given the energy of configuration, trajectory of simulation and position of atoms in each frame

Question

Background

I am using a neural network to calculate the potential energy of atoms in a configuration and then adding energy of all atoms to compare it with the true energy of the configuration(label) to update the weights.

I am doing so for 10000 configurations of 100 atoms each for liquid Argon at 100K.

Note I don't have the actual energy of each atom and only the energy of the configuration.

Problem

I need to calculate the forces on each atom to see if the sum of forces for all the atoms is zero in X,Y,Z axis. In case it is not I want to find out how off is it from zero.

How do I calculate the force on atoms using their position and potential energy for a given trajectory. I know I need to do $\frac{dE}{dr}$ but what should be my E and r. Also since I don,t a formula for E(r), what should my $dE$ and $dr$ be respectively?

Answer 1

tl,dr: Forces are equivalent to the negative of the gradient, use small, arbitrary $\mathrm{d}r$ for each atom

The problem, as you mention, with many neural network methods in chemistry at the moment, is that the $E = f(x,y,z)$ is a black box and thus analytical gradients (and thus forces) are often not available.

First, as you clearly know, $F = -\nabla E(x,y,z)$ so you need to solve $\frac{\mathrm{d}E}{\mathrm{d}r}$ for your neural network. (I indicate this for anyone reading who doesn't realize that the force will inherently be the negative of the gradient of the energy function.

So what you need is a finite difference method. This is, of course, a common problem (e.g., on Math Stack Exchange or Wikipedia):

$$ \frac{\partial E}{\partial x} \approx \frac{f(x+\epsilon)-f(x)}{\epsilon}$$

Picking a good step $\epsilon$ is a bit tricky, but I'd probably pick something like 0.01 Å or 0.001 Å, as a start. In your case, it sounds like you want Cartesian gradients, so you'll take your $f(x,y,z)$ and perform the partial derivatives for all three axes (i.e., the initial point $(x,y,z)$ plus 3 extra energy evaluations per atom.

That is, you'd loop through each atom, move it a bit in the x, y, and z directions and evaluate the system to get the gradient/force on that atom. You'll, therefore, have $3N$ energy evaluations to get all the forces on all the atoms.

You can perform better (and more computationally expensive) finite difference methods if you choose.