# Understanding LAMMPS Source Codes: A Study Note

This is a note about learning `LAMMPS`

source codes. This note focuses on compute style of `Lammps`

which is used to compute certain quantity during the simulation run. Of course you can as well compute these quantities in post-process, however it’s usually faster to do it in the simulation since you can take advantage of the all the distance, forces, et al generated during the simulation instead of computing them again in post-process. I will go through the `LAMMPS`

source code `compute_gyration.h`

and `compute_gyration.cpp`

. I am not very familiar with `c++`

, so I will also explain some language related details which is what I learn when studying the code. Hope this article can be helpful when someone want to modify or make their own `Lammps`

compute style.

`compute_gyration.h`

```
#ifdef COMPUTE_CLASS
ComputeStyle(gyration,ComputeGyration)
#else
#ifndef LMP_COMPUTE_GYRATION_H
#define LMP_COMPUTE_GYRATION_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeGyration : public Compute {
public:
ComputeGyration(class LAMMPS *, int, char **);
~ComputeGyration();
void init();
double compute_scalar();
void compute_vector();
private:
double masstotal;
};
}
#endif
#endif
```

First part of this code

```
#ifdef COMPUTE_CLASS
ComputeStyle(gyration,ComputeGyration)
#else
```

is where this speciﬁc compute style is deﬁned. If you want to write your own compute style, let’s say **intermediate scattering function**. Then we write like this

```
#ifdef COMPUTE_CLASS
ComputeStyle(isf,ComputeISF) // ISF stands for intermediate scattering function
#else
```

Move to the rest part. `#include "compute.h"`

and `namespace LAMMPS_NS`

is to include the base class and namespace. Nothing special is here, you need to have this in every speciﬁc compute style header ﬁle.

```
class ComputeGyration : public Compute {
public:
ComputeGyration(class LAMMPS *, int, char **);
~ComputeGyration();
void init();
double compute_scalar();
void compute_vector();
private:
double masstotal;
```

You can see there is a overal structure in the above code `class A : public B`

. This basically means that our derived class A will inherit all the public and protected member of class B. More details can be found here, here and here

Next, we declare two types of member of our derived class, `public`

and `private`

. `public`

is the member we want the other code can access to and `private`

is the member which is only used in the derived class scope. Now let’s look at the public class member. Note that there is no type declaration of class member `ComputeGyration`

and `~ComputeGyration`

. They are called *Class Constructor* and *Class Destructor*. They are usually used to set up the initial values for certain member variables as we can see later in `compute_gyration.cpp`

. Note that for some compute style such as `compute_msd.h`

, the destructor is virtual, that is `virtual ~ComputeMSD`

instead of just `~ComputeMSD`

. This is because class `ComputeMSD`

is also inherited by derive class `ComputeMSDNonGauss`

. So you need to decalre the base destructor as being virtual. Look at this page for more details. Now let’s move forward.

```
void init();
double compute_scalar();
void compute_vector();
```

here all the function `init`

, `compute_scalar`

and `compute_vector`

all are the base class member which are already deﬁned in `compute.h`

. However they are all virtual functions, which means that they can be overrided in the derived class, here it is the `ComputeGyration`

. This and this pages provide some basic explanations for the use of virtual functions. Here is a list shown in LAMMPS documentation of **some examples** of the virtual functions you can use in your derived class.

In our case, gyration computation will return a scalor and a vector, then we need `compute_scalar()`

and `compute_vector()`

. Private member `masstotal`

is the quantity calculated locally which is only used within the class and not needed for the rest of the codes.

`compute_gyration.cpp`

Now let’s look at the `compute_gyration.cpp`

.

```
#include <math.h>
#include "compute_gyration.h"
#include "update.h"
#include "atom.h"
#include "group.h"
#include "domain.h"
#include "error.h"
```

Here the necessary header ﬁles are include. The name of these header ﬁle is self-explanary. For instance, `updata.h`

declare the functions to update the timestep, et al.

```
ComputeGyration::ComputeGyration(LAMMPS *lmp, int narg, char **arg) :
Compute(lmp, narg, arg)
{
if (narg != 3) error->all(FLERR,"Illegal compute gyration command");
scalar_flag = vector_flag = 1;
size_vector = 6;
extscalar = 0;
extvector = 0;
vector = new double[6];
}
```

The above code deﬁne the what the constructor `ComputeGyration`

actually does. `::`

is called *scope operator*, it is used to specify that the function being deﬁned is a member (in our case which is the constructor which has the same name as the its class) of the class `ComputeGyration`

and not a regular non-member function. The structure `ComputeGyration : Compute()`

is called a **Member Initializer List**. It initializes the member `Compute()`

with the arguments `lmp, narg, arg`

. `narg`

is the number of arguments provided. `scalar_flag`

, `vector_flag`

, `size_vector`

, `extscalar`

and `extvector`

all are the ﬂags parameter deﬁned in `Compute.h`

. For instance, `scalar_flag = 1/0`

indicates we will/won’t use function `compute_scalar()`

in our derived class. The meaning of each parameter is explained in `compute.h`

. This line `vector = new double[6]`

is to dynamically allocate the memory for array of length 6. Normally the syntax of **new** operator is such

```
double *vector = NULL;
vector = new double[6];
```

Here the line `double *vector = NULL`

is actually in `compute.h`

and `compute.cpp`

. Where pointer `vector`

is deﬁned in `compute.h`

and its value is set to `NULL`

in `compute.cpp`

.

```
ComputeGyration::~ComputeGyration()
{
delete [] vector;
}
```

The above code speﬁcy destructor that is what will be excuted when class `ComputeGyration`

goes out of scope or is deleted. In this case, it delete the gyration tensor vector deﬁned above. The syntax of **delete** operator for array is `delete [] vector`

. For details of **new** and **delete** can be found here.

```
void ComputeGyration::init()
{
masstotal = group->mass(igroup);
}
```

This part perform one time setup like initialization. Operator **->** is just a syntax sugar, `class->member`

is equivalent with `(*class).member`

. What `group->mass(igroup)`

does is to call the member `mass()`

function of class `group`

, provided the group-ID, and return the total mass of this group. How value of `igroup`

is set can be examined in `compute.cpp`

. It’s the second argument of compute style.

```
double ComputeGyration::compute_scalar()
{
invoked_scalar = update->ntimestep;
double xcm[3];
group->xcm(igroup,masstotal,xcm);
scalar = group->gyration(igroup,masstotal,xcm);
return scalar;
}
```

`invoked_scalar`

is deﬁned in base class `Compute`

. The value is the last timestep on which `compute_scalar()`

was invoked. `ntimestep`

is the member of class `update`

which is the current timestep. `xcm`

function of class `group`

calculate the center-of-mass coords. The result will be stored in `xcm`

. `gyration`

function calculate the gyration of a group given the total mass and center of mass of the group. The total mass is calculated in `init()`

. And in order for it to be accessed here, it is deﬁned as private in `compute_gyration.h`

. Notice that here there is no explicit code to calculte the gyration scalor because the member function which does this job is already deﬁned in class `group`

. So we just need to call it. However we also want to calculate the gyration tensor, we need to write a function to calculate it.

```
void ComputeGyration::compute_vector()
{
invoked_vector = update->ntimestep;
double xcm[3];
group->xcm(igroup,masstotal,xcm);
double **x = atom->x;
int *mask = atom->mask;
int *type = atom->type;
imageint *image = atom->image;
double *mass = atom->mass;
double *rmass = atom->rmass;
int nlocal = atom->nlocal;
double dx,dy,dz,massone;
double unwrap[3];
double rg[6];
rg[0] = rg[1] = rg[2] = rg[3] = rg[4] = rg[5] = 0.0;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
if (rmass) massone = rmass[i];
else massone = mass[type[i]];
domain->unmap(x[i],image[i],unwrap);
dx = unwrap[0] - xcm[0];
dy = unwrap[1] - xcm[1];
dz = unwrap[2] - xcm[2];
rg[0] += dx*dx * massone;
rg[1] += dy*dy * massone;
rg[2] += dz*dz * massone;
rg[3] += dx*dy * massone;
rg[4] += dx*dz * massone;
rg[5] += dy*dz * massone;
}
MPI_Allreduce(rg,vector,6,MPI_DOUBLE,MPI_SUM,world);
if (masstotal == 0.0) return;
for (int i = 0; i < 6; i++) vector[i] = vector[i]/masstotal;
}
```

The above code do the actual computation of gyration tensor.

Here is the list of meaning of each variable

`x`

: 2D array of the position of atoms.`mask`

: array of group information of each atom.`if (mask[i] & groupbit)`

check whether the atom is in the group on which we want to perform calculation.`type`

: type of atom.`image`

: image ﬂags of atoms. For instance a value of 2 means add 2 box lengths to get the unwrapped coordinate.`mass`

: mass of atoms.`rmass`

: mass of atoms with ﬁnite-size (meaning that it can have rotational motion). Notice that mass of such particle is set by density and diameter, not directly by the mass. That’s why they set two variables`rmass`

and`mass`

. To extract mass of atom`i`

, use`rmass[i]`

or`mass[type[i]]`

.`nlocal`

: number of atoms in one processor.

Look at this line `domain->unmap(x[i],image[i],unwrap)`

, `domain.cpp`

tells that function `unmap`

return the unwrapped coordination of atoms in `unwrap`

. The following several lines calculate the gyration tensor. The MPI code `MPI_Allreduce(rg,vector,6,MPI_DOUBLE,MPI_SUM,world)`

sums all the six components of `rg`

calculated by each processor, store the value in `vector`

and then distribute `vector`

to all the processors. Refer to this article for details.

Here are two good articles about understanding and hacking LAMMPS code.