This tutorial will describe how to implement a new pair potential and
test it. Let us say we want to implement a Gaussian core potential:

(3) |

When implementing a new pair potential there are a number of things to do:

- Implement the potential and force into the RUMD source code
- Include the new potential in the Python interface
- Compile
- Test the potential
- Use it

The `PairPotential.h`

in the directory <RUMD-HOME>/include /rumd
contains all the pair potentials, implemented as classes in C++.

We start with a review of the general 12-6 Lennard-Jones potential - it is
the first pair potential in the file (search for `Pot_IPL_12_6`

). Note
that this is not the one most typically used in simulations; that would be
one of its derived classes `Pot_LJ_12_6`

. The difference is only in how
the parameters are specified: for the `Pot_IPL_12_6`

potential, rather
than `Epsilon`

and `Sigma`

, we specify the coefficients of the IPL
terms, `A12`

and `A6`

. We will
make a copy of it and change it to the Gaussian core potential.

The implementation is split into two major parts: The first part is the
function `SetParams(...)`

which is called by the user via the
python interface. The last part is the function `ComputeInteraction`

where the actual computation (on the GPU) is taking place.

`SetParams(...)`

is called via the python interface once for
each type of interaction, with `i`

and `j`

denoting the types
of the two particles involved in the particular
interaction. The user specifies the coefficients and
and the cut-off (in other potentials the parameters might include a
characteristic length and
characteristic energy and possibly more parameters).
The parameters of the potential are communicated by storing them in
`h_params[i][j][]`

. NOTE: the correctness of the program relies on `h_params[i][j][0]`

containing the cut-off in absolute units.
In other words, interactions between particles of type `i`

and `j`

are only computed if the distance between the two particles is less than
`h_params[i][j][0]`

. The rest
of the stored parameters are chosen to facilitate fast computation.

The actual computation of the potential and the associated force
is done on the GPU by the
function `ComputeInteraction(...)`

. When this is called
`dist2`

contains the square of the distance between the two
particles in question, and `param`

contains the appropriate
parameters as set up by ` SetParams(...)`

. From this information,
the function calculates (i.e. the force, see below in equation
(5)) and the potential. In the Lennard-Jones
potential it is convenient to define two inverse lengths
`invDist2`

and `invDist6`

to minimize the number of
calculations. With these
lengths it is easy to calculate and the potential energy.
`(*my_f).w +=`

is the contribution to the potential energy.

We will now proceed to implement a Gaussian core potential. First
step is to write up the equations we will implement. We will here
implement the Gaussian core potential (repeated below)

(4) |

We have to provide the function which is used
by RUMD to calculate the force:

Now copy the entire Lennard-Jones potential and paste it somewhere in the file
`PairPotential.h`

. Replace `Pot_IPL_12_6`

with `Pot_Gaussian`

(four places).
Choose a default identity string, e.g. `SetID_String("potGauss");`

.

We will change the arguments in the ` SetParams(...)`

to match those of
the other LJ potential `Pot_LJ_12_6`

which uses and
instead of and , since
the Gaussian potential also has a characteristic length and a
characteristic energy . But we need to make corresponding
changes to the `h_params[i][j][]`

's. For the Gaussian potential they
become:

`h_params[i][j][0] = Rcut * Sigma;`

`h_params[i][j][1] = 1.f / ( Sigma * Sigma );`

`h_params[i][j][2] = Epsilon;`

`Rcut`

is to be interpreted as a scaled
cutoff, in units of the corresponding .
Now we can write the part taking place on the GPU. Since we do
not need any inverse distances delete the lines computing
`invDist2`

and `invDist6`

. We define instead a
float called `Dist2`

:
float Dist2 = dist2 * param[1];which is . To avoid computing the exponential more than once, we define:

float exponential = exp(-Dist2);With these definitions, is simply implemented as:

float s = 2.f * param[1] * param[2] * exponential;The next line is the potential:

(*my_f).w += param[2] * exponential;

This finishes the actual implementation, and we now proceed with a number of more technical changes ensuring that RUMD knows your new pair potential exists.

- In the file
`Potentials.i`

in the directory <RUMD-HOME>/Swig copy the entry (see below) relating to`Pot_IPL_12_6`

, replace`Pot_IPL_12_6`

with your potential name, change the arguments to`SetParams`

in accordance with how it was defined in`PairPotential.h`

, and replace the docstring with a similar brief description of your potential. To be 100% clear about what needs to be copied, here is the entry for`Pot_IPL_12_6`

:// docstring for Pot_IPL_12_6 %feature("autodoc","Standard 12-6 Lennard-Jones (cut and shifted potential) v(r) = ( A12/r^12 + A6/r^6 ) - v_cut s(r) = 12 A12/r^14 + 6 A6/r^8 = -r^-1*dv/dr (F_ij = s * r_ij)") Pot_IPL_12_6; class Pot_IPL_12_6 : public PairPotential { public: Pot_IPL_12_6( CutoffMethod cutoff_method ); void SetParams(unsigned int i, unsigned int j, float A12, float A6, float Rcut); };

- In the directory <RUMD-HOME> do

`make`

If everything compiled without errors the new pair potential will now be available from the python interface, and you need to test it. Run a simulation as described in the Level 1 tutorial, but choose the new potential.

Check that the simulation ran with the correct potential by plotting
it by adding
`potential.WritePotentials(sim.sample.GetSimulationBox())`

to your script. This command produces a file ''potentials_`PotID`

.dat'' where the potential is followed by the radial force in each column for all
possible particle pairs.

Make sure that the cut-off is done correctly by the program by zooming in on that distance. An example is shown in figure 5.

Next, run a NVE simulation and test that the total energy is conserved (i.e., have fluctuations that are much smaller than for the potential and kinetic energy).Have fun with RUMD!

Heine Larsen 2017-07-21