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!