Understanding LAMMPS Source Codes: A Study Note

This is a note about learn­ing LAMMPS source codes. This note fo­cuses on com­pute style of LAMPPS which is used to com­pute cer­tain quan­tity dur­ing the sim­u­la­tion run. Of course you can as well com­pute these quan­ti­ties in post-process, how­ever it’s usu­ally faster to do it in the sim­u­la­tion since you can take ad­van­tage of the all the dis­tance, forces, et al gen­er­ated dur­ing the sim­u­la­tion in­stead of com­put­ing them again in post-process. I will go through the LAMMPS source code compute_gyration.h and compute_gyration.cpp. I am not very fa­mil­iar with c++, so I will also ex­plain some lan­guage re­lated de­tails which is what I learn when study­ing the code. Hope this ar­ti­cle can be help­ful when some­one want to mod­ify or make their own Lammps com­pute 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 spe­cific com­pute style is de­fined. If you want to write your own com­pute style, let’s say in­ter­me­di­ate scat­ter­ing func­tion. 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 in­clude the base class and name­space. Nothing spe­cial is here, you need to have this in every spe­cific com­pute style header file.

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 struc­ture in the above code class A : public B. This ba­si­cally means that our de­rived class A will in­herit all the pub­lic and pro­tected mem­ber of class B ^[1][2][3].

Next, we de­clare two types of mem­ber of our de­rived class, public and private. public is the mem­ber we want the other code can ac­cess to and private is the mem­ber which is only used in the de­rived class scope. Now let’s look at the pub­lic class mem­ber. Note that there is no type de­c­la­ra­tion of class mem­ber ComputeGyration and ~ComputeGyration. They are called Class Constructor and Class Destructor. They are usu­ally used to set up the ini­tial val­ues for cer­tain mem­ber vari­ables as we can see later in compute_gyration.cpp. Note that for some com­pute style such as compute_msd.h, the de­struc­tor is vir­tual, that is virtual ~ComputeMSD in­stead of just ~ComputeMSD. This is be­cause class ComputeMSD is also in­her­ited by de­rive class ComputeMSDNonGauss. So you need to de­calre the base de­struc­tor as be­ing vir­tual. Look at this page for more de­tails. Now let’s move for­ward.

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

here all the func­tion init, compute_scalar and compute_vector all are the base class mem­ber which are al­ready de­fined in compute.h. However they are all vir­tual func­tions ^[4][5], which means that they can be over­rided in the de­rived class, here it is the ComputeGyration. Here is a list shown in LAMMPS doc­u­men­ta­tion of some ex­am­ples of the vir­tual func­tions you can use in your de­rived class.

In our case, gy­ra­tion com­pu­ta­tion will re­turn a scalor and a vec­tor, then we need compute_scalar() and compute_vector(). Private mem­ber masstotal is the quan­tity cal­cu­lated lo­cally 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 nec­es­sary header files are in­clude. The name of these header file is self-ex­pla­nary. For in­stance, updata.h de­clare the func­tions to up­date 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­fine the what the con­struc­tor ComputeGyration ac­tu­ally does. :: is called scope op­er­a­tor, it is used to spec­ify that the func­tion be­ing de­fined is a mem­ber (in our case which is the con­struc­tor which has the same name as the its class) of the class ComputeGyration and not a reg­u­lar non-mem­ber func­tion. The struc­ture ComputeGyration : Compute() is called a Member Initializer List. It ini­tial­izes the mem­ber Compute() with the ar­gu­ments lmp, narg, arg. narg is the num­ber of ar­gu­ments pro­vided. scalar_flag, vector_flag, size_vector, extscalar and extvector all are the flags pa­ra­me­ter de­fined in Compute.h. For in­stance, scalar_flag = 1/0 in­di­cates we will/​won’t use func­tion compute_scalar() in our de­rived class. The mean­ing of each pa­ra­me­ter is ex­plained in compute.h. This line vector = new double[6] is to dy­nam­i­cally al­lo­cate the mem­ory for ar­ray of length 6. Normally the syn­tax of new op­er­a­tor is such

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

Here the line double *vector = NULL is ac­tu­ally in compute.h and compute.cpp. Where pointer vector is de­fined in compute.h and its value is set to NULL in compute.cpp.

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

The above code spec­ify de­struc­tor that is what will be ex­cuted when class ComputeGyration goes out of scope or is deleted. In this case, it delete the gy­ra­tion ten­sor vec­tor de­fined above. The syn­tax of delete op­er­a­tor for ar­ray is delete [] vector. See de­tails of new and delete.

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

This part per­form one time setup like ini­tial­iza­tion. Operator -> is just a syn­tax sugar, class->member is equiv­a­lent with (*class).member. What group->mass(igroup) does is to call the mem­ber mass() func­tion of class group, pro­vided the group-ID, and re­turn the to­tal mass of this group. How value of igroup is set can be ex­am­ined in compute.cpp. It’s the sec­ond ar­gu­ment of com­pute 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­fined in base class Compute. The value is the last timestep on which compute_scalar() was in­voked. ntimestep is the mem­ber of class update which is the cur­rent timestep. xcm func­tion of class group cal­cu­late the cen­ter-of-mass co­ords. The re­sult will be stored in xcm. gyration func­tion cal­cu­late the gy­ra­tion of a group given the to­tal mass and cen­ter of mass of the group. The to­tal mass is cal­cu­lated in init(). And in or­der for it to be ac­cessed here, it is de­fined as pri­vate in compute_gyration.h. Notice that here there is no ex­plicit code to cal­culte the gy­ra­tion scalor be­cause the mem­ber func­tion which does this job is al­ready de­fined in class group. So we just need to call it. However we also want to cal­cu­late the gy­ra­tion ten­sor, we need to write a func­tion to cal­cu­late 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 ac­tual com­pu­ta­tion of gy­ra­tion ten­sor.

Here is the list of mean­ing of each vari­able

• x: 2D ar­ray of the po­si­tion of atoms.
• mask: ar­ray of group in­for­ma­tion of each atom. if (mask[i] & groupbit) check whether the atom is in the group on which we want to per­form cal­cu­la­tion.
• type: type of atom.
• image: im­age flags of atoms. For in­stance a value of 2 means add 2 box lengths to get the un­wrapped co­or­di­nate.
• mass: mass of atoms.
• rmass: mass of atoms with fi­nite-size (meaning that it can have ro­ta­tional mo­tion). Notice that mass of such par­ti­cle is set by den­sity and di­am­e­ter, not di­rectly by the mass. That’s why they set two vari­ables rmass and mass. To ex­tract mass of atom i, use rmass[i] or mass[type[i]].
• nlocal: num­ber of atoms in one proces­sor.

Look at this line domain->unmap(x[i],image[i],unwrap), domain.cpp tells that func­tion unmap re­turn the un­wrapped co­or­di­na­tion of atoms in unwrap. The fol­low­ing sev­eral lines cal­cu­late the gy­ra­tion ten­sor. The MPI code MPI_Allreduce(rg,vector,6,MPI_DOUBLE,MPI_SUM,world) sums all the six com­po­nents of rg cal­cu­lated by each proces­sor, store the value in vector and then dis­trib­ute vector to all the proces­sors. Refer to this ar­ti­cle for de­tails.

Unfortunately, there are only a few re­sources on LAMMPS source code on the Internet. Here are two good ar­ti­cles about un­der­stand­ing and hack­ing LAMMPS code.

• A dis­sec­tion of LAMMPS classes - A Guide for Beginner LAMMPS Hackers
• How to Write Fix for LAMMPS