Parallel Implementation of Nonadditive Gaussian Process Potentials for Monte Carlo Simulations

A strategy is presented to implement Gaussian process potentials in molecular simulations through parallel programming. Attention is focused on the three-body nonadditive energy, though all algorithms extend straightforwardly to the additive energy. The method to distribute pairs and triplets between processes is general to all potentials. Results are presented for a simulation box of argon, including full box and atom displacement calculations, which are relevant to Monte Carlo simulation. Data on speed-up are presented for up to 120 processes across four nodes. A 4-fold speed-up is observed over five processes, extending to 20-fold over 40 processes and 30-fold over 120 processes.

r 12 < r 13 < r 23 and denoting the cut-off as r c , there are three types of triplets that require non-additive LRCs: • r 12 < r 13 < r c < r 23 , denoted Type 1.
Because all interacting species are atomic they can always be rotated into the same plane, though the following method could be extended to molecules regardless. S1 r 13 > r c r 23 > r c r 12 < r c θ x 1 2 3 Figure 1: A sketch of a Type 2 triplet that shows the angle θ between r 12 and the distance x.
The Type 1 scenario gives rise to a area of integration that is not trivial to define. As such, no explicit solution to this scenario is given here. Instead it is suggested that points where r max 23 = r max 13 + r 12 are included in the training data, to permit use of the GP potential for prediction.

Type 2
Under the Type 2 scenario, the LRC accounts for the effect of 3 on the 1-2 interaction when the former is at long-range. Thus the Type 2 LRC U LRC T2 can be expressed as an integral over the position of atom 3 only. In terms of interatomic distances this therefore requires integrating over both r 13 and r 23 . As r 13 < r 23 , atom 3 must be within the hemisphere around atom 1 that leaves it closer to 1 than 2. r 13 must therefore be between the cut-off and infinity, whilst r 23 must exceed r 13 and be less than r 12 + r 13 to obey the triangle rule. Figure 1 shows that these distances can be re-written in terms of the angle θ and distance x. This leaves the integral as Here, x min is the value of x at which r 23 = r c and ρ is the density of atoms in the simulation box. The dependence of the Type 2 correction on r 12 stems from these two atoms being S2 within r c of each other. Thus this integral could be evaluated for different values of r 12 and stored in a look-up table, which would not have to undergo an update following a volume change or exchange move.

Type 3
The Type 3 scenario is similar to Type 2 as illustrated in figure 1, with the single difference being that r c < r 12 in the former. An equivalent to equation 1 can therefore be written down for the Type 3 correction. However, because r 12 is now greater than the cut-off, this distance must also be integrated over. Therefore, using equivalent notation, This integral is a purely mean field correction. Consequently, it must only be evaluated at the outset of the simulation and whenever the volume changes. The integrals in equations 1 and 2 represent a general approach to solving the long-range correction for three-body interactions as they can employ any function for E (3) .

Calculating forces
Forces are a pre-requisite for molecular dynamics simulations. It is a useful property of GP potentials that forces can be obtained by differentiating the kernel function without the need for additional training. This process is outlined for a squared exponential kernel here.
For an atom δ, forces must be found along each of the x, y and z axes. Letting d denote an inverse interatomic distance for clarity, the force on atom δ along axis η is As the forces are determined from the additive energy, N d = N perm = 1. The sum over atoms S3 represents the force being an accumulation of the forces exerted by other atoms along axis η.
As d is an inverse distance, the substitution can be used. Thus, writing k = exp − v 2 2l 2 , the derivative in equation (3) becomes Solving equation (6), evaluating ∂k ∂v and substituting the results into equation (5) gives Thus the force along axis η is The form of equation (8) implies that the same strategy used to calculate the energy can be applied to forces. Thus it is anticipated that the parallel speed-up will be similar for the force and energy calculations, but this has not yet been tested.
Direct training on quantum-chemical force data is also possible, with the above method suggested to reduce the number of expensive calculations necessary to train a potential for simulation. The method of differentiating the kernel function is likely to give accurate forces as the correlation length is a hyper-parameter of the model, which can be used to verify S4 sufficient local smoothness of the model to give confidence in the derivatives. For all of GP PES discussed in the paper the correlation lengths are consistent with derivatives being sufficiently accurate. However, quantitative testing against ab initio force data would be required to confirm this expectation. S5