How to Modify LAMMPS: From the Prospective of a Particle Method Researcher

LAMMPS is a powerful simulator originally developed for molecular dynamics that, today, also accounts for other particle-based algorithms such as DEM, SPH, or Peridynamics. The versatility of this software is further enhanced by the fact that it is open-source and modifiable by users. This property suits particularly well Discrete Multiphysics and hybrid models that combine multiple particle methods in the same simulation. Modifying LAMMPS can be challenging for researchers with little coding experience. The available material explaining how to modify LAMMPS is either too basic or too advanced for the average researcher. In this work, we provide several examples, with increasing level of complexity, suitable for researchers and practitioners in physics and engineering, who are familiar with coding without been experts. For each feature, step by step instructions for implementing them in LAMMPS are shown to allow researchers to easily follow the procedure and compile a new version of the code. The aim is to fill a gap in the literature with particular reference to the scientific community that uses particle methods for (discrete) multiphysics.


Introduction
LAMMPS, acronym for Large-scale Atomic/Molecular Massively Parallel Simulator, was originally written in F77 by Steve Plimpton [1] in 1993 with the goal of having a large-scale parallel classical Molecular Dynamic (MD) code. The project was a Cooperative Research and Development Agreement (CRADA) between two DOE labs (Sandia and LLNL) and three companies (Cray, Bristol Myers Squibb, and Dupont). Since the initial release LAMMPS has been improved and expanded by many researchers who implemented many mesh-free computational methods such as Perydynamics, Smoothed particle hydrodynamics (SPH), Discrete Elemet Method (DEM) and many more [2][3][4][5][6][7][8][9][10][11].
Such a large number of computational methods within the same simulator allows researchers to easily combine them for the simulation of complex phenomena. In particular, our research group has used during the years LAMMPS in a variety of settings that go from • The private members, defined before the keyword public, cannot be accessed from outside the class. They can only be accessed by class or "friend" functions, which are declared as having access to class members, without themselves being members. All the class members are private by default. • The public members can be accessed from outside the class anywhere within the scope of the class object. • The protected members are similar to private members but they can be accessed by derived classes or child classes while private members cannot.
The type of inheritance is specified by the access-specifier, one of public, protected, or private. If the access-specifier is not used, then it is private by default, but public inheritance is commonly used: public members of the base class become public members of the derived class and protected members of the base class become protected members of the derived class. A base class's private members are never accessible directly from a derived class, but can be accessed through calls to the public and protected members of the base class.

Virtual Function
The signature of a function f must be declared with a virtual keyword in a base class C to allow its definition (implementation), or redefinition, in a derived class D. Then, when a derived class D object is used as an element of the base class C, and f is called, the derived class's implementation of the function is executed.
There is nothing wrong with putting the virtual in front of functions inside of the derived classes, but it is not required, unless it is known for sure that the class will not have any children who would need to override the functions of the base class. A class that declares or inherits a virtual function is called a polymorphic class.

LAMMPS Inheritance and Class Syntax
A schematic representation of the LAMMPS inheritance tree is shown in Figure 1: LAMMPS is the top-level class for the entire code, then all the core classes, highlighted in blue, inherit all the constructors, destructors, assignment operator members, friends and private members declared and defined in LAMMPS. The core classes perform LAMMPS fundamental actions. For instance, the Atom class collects and stores all the per-atom, or per-particle, data while Neighbor class builds the neighbor lists [41].
The style classes, highlighted in reds, inherit all the constructors, destructors, assignment operator members, friends and private members declared and defined in LAMMPS and in the corresponding core class. The style classes are also virtual parents class of many child classes that implement the interface defined by the parent class. For example, the fix style has around 100 child classes.
Each style is composed of a pair of files: • namestyle.h The header of the style, where the class style is defined and all the objects, methods and constructors are declared. • namestyle.cpp Where all the objects, methods and constructors declared in the class of style are defined.
When a new style is written both namestyle.h and namestyle.cpp files need to be created.
Each "family" style has its own set of methods, declared in the header and defined in the cpp file, in order to define the scope of the style. For example, the pair style are classes that set the formula(s) LAMMPS uses to compute pairwise interactions while bond style set the formula(s) to compute bond interactions between pairs of atoms [41].
Each pair style has some recurrent functions such as compute, allocate and coeff. Although the final scope of those functions can differ for different styles, they all share a similar role within the classes. When a new style is written both namestyle.h and namestyle.cpp files need to be created.
Each "family" style has its own set of methods, declared in the header and defined in the cpp file, in order to define the scope of the style. For example, the pair style are classes that set the formula(s) LAMMPS uses to compute pairwise interactions while bond style set the formula(s) to compute bond interactions between pairs of atoms [41]. Each pair style has some recurrent functions such as compute, allocate and coeff. Although the final scope of those functions can differ for different styles, they all share a similar role An example of a pair style, sph/taitwater, header in LAMMPS is shown in Listing 2.
Listing 2: Header file of sph/taitwater pair style (pair_sph_taitwater.h) All the class members are defined in the cpp file. Taking sph/taitwater pair style as reference, each method declared in Listing 2 will be defined and commented in the next sections. Although this can be style-specific, the aim is to give an overview of how the methods are defined in the cpp in LAMMPS. Albeit different style has different methods, the understanding gained can be transferred into others style, as shown in Sections 3 and 6.

Constructor
Any class usually include a member function called constructors. The constructor is mechanically invoked when an object of the class is created. This allows the class to initialise members or allocate storage. Unlike the other member of the class, the constructor name must match the name of the class and it does not have a return type.

Destructor
The role of destructors is to de-allocate the allocated dynamic memory, see Section 2.3.8, being mechanically invoked just before the end of the class lifetime. Similarly to constructors, destructors does not have a return type and have the same name as the class name with a tilde (∼) prefix.  int i, j, ii, jj, inum, jnum, itype, jtype; 5 double xtmp, ytmp, ztmp, delx, dely, delz, fpair; 6 7 int *ilist, *jlist, *numneigh, **firstneigh; 8 double vxtmp, vytmp, vztmp, imass, jmass, fi, fj, fvisc, h, ih, ihsq; 9 double rsq, tmp, wfd, delVdotDelR, mu, deltaE; 10 // end 11 12 if (eflag || vflag) ChemEngineering 2021, 5, 30 6 of 57 13 ev_setup(eflag, vflag); 14 else 15 evflag = vflag_fdotr = 0; 16 17 /// others variables and pointers declaration and initialisation 18 double **v = atom->vest; // pass the value of the pointer that points to a pointers 19 // pointing to the first element of velocity vector of the particles 20 double **x = atom->x; // pass the value of the pointer that points to a pointers 21 // pointing to the first element of position vector of the particles 22 double **f = atom->f; // pass the value of the pointer that points to a pointers 23 // pointing to the first element of force vector of the particles 24 double *rho = atom->rho; // pass the value of the pointer that points 25 // to the density vector of the particles 26 double *mass = atom->mass; // pass the value of the pointer that points 27 // to the mass vector of the particles 28 double *de = atom->de; // pass the value of the pointer that points 29 // to the change of internal energy of the particles 30 double *drho = atom->drho; // pass the value of the pointer that points 31 // to the change of density of the particles 32 int *type = atom->type; // pass the value of the pointer that points to the type of the particles

coeff
Similar to setting, coeff is a public void function that reads and set the coefficients used in by compute of the pair style. For each i j pair is possible to set different coefficients. The coefficients are passed in the input file with the command line pair coeff, see Listing 10. As before, examples for different pair coeff input script and the corresponding coeff are listed below: • sph/taitwater As described in the SPH for LAMMPS manual [6], the command line to invoke sph/taitwater pair coeff is shown in Listing 10. • sph/rhosum As described in the SPH for LAMMPS manual [6], the syntax to invoke the command is shown in Listing 12.

single
In single the force and energy of a single pairwise interaction, or single bond or angle (in case of bond or angle style), between two atoms is evalutated. The method is specifically invoked by the command line compute pair/local (or compute bond/local) to calculate properties of individual pair, or bond, interactions [41]. allocate is a protected void function that allocates dynamic memory. The dynamic memory allocation is used when the amount of memory needed depends on user input. As explained before, at the end of the lifetime of the class, the destructors will de-allocate the memory the memory used by allocate.   memory->create(cutsq, n + 1, n + 1, "pair:cutsq"); 12 memory->create(rho0, n + 1, "pair:rho0"); 13 memory->create(soundspeed, n + 1, "pair:soundspeed"); 14 memory->create(B, n + 1, "pair:B"); 15 memory->create(cut, n + 1, n + 1, "pair:cut"); 16 memory->create(viscosity, n + 1, n + 1, "pair:viscosity"); 17 }

Kelvin-Voigt Bond Style
We can use what we learned in the previous section to generate a new dissipative bond potential that can be used to model viscoelastic materials. The Kelvin-Voigt model [44] is used to model viscoelastic material as a purely viscous damper and purely elastic spring connected in parallel as shown in Figure 2.
Since the two components of the model are arranged in parallel, the strain in each component is identical: On the other hand, the total stress σ tot will be split into σ spring and σ damper to have ε spring = ε damper . Thus we have Combining Equations (1) and (2) with the constitutive relation for both the spring and the dumper, σ spring = kε and σ damper = bε , is possible to write that where k is the elastic modulus and b is the coefficient of viscosity. Equation (3) relates stress to strain and strain rate for a Kelvin-Voigt material [44]. Similarly to bond test to write a new pair style called bond kv we take the bond harmonic pair style as reference. The new pair style is declared and initialised in bond_kv.h and bond_kv.cpp saved in the /src/MOLECULE directory and its hierarchy is shown in Figure 3.
where k is the elastic modulus and b is the coefficient of viscosity. Equation (3) relates stress to strain and strain rate for a Kelvin-Voigt material [44].

Validation
The bond kv pair style has been validated by Sahputra et al. [45] in their Discrete Multiphysics model for encapsulate particles with a soft outer shell.

bond_kv.cpp
All the functions will be the same as in the reference bond harmonic. However, in our new bond kv, we need to substitute the "BondHarmonic" text by a new "BondKv" text, as can be seen in Listing 17 and 18. Form now on, when we show a side-by-side comparison between the reference and the modified file, we highlight in yellow the modified lines and in red the deleted lines.

Validation
The bond kv pair style has been validated by Sahputra et al. [45] in their Discrete Multiphysics model for encapsulate particles with a soft outer shell.

bond_kv.cpp
All the functions will be the same as in the reference bond harmonic. However, in our new bond kv, we need to substitute the "BondHarmonic" text by a new "BondKv" text, as can be seen in Listings 17 and 18. Form now on, when we show a side-by-side comparison between the reference and the modified file, we highlight in yellow the modified lines and in red the deleted lines.

bond_kv.h
In the header of the new pair style we need to substitute the "BondHarmonic" text by a new "BondKv" text as well as declare a new protected member in the class, the pointer to b.

Noble-Abel Stiffened-Gas Pair Style
In the SPH framework is possible to determine all the particles properties by solving the particle form of the continuity equation [6,26] the momentum equation [6,26] and the energy conservation equation [6,26] However, to be able to solve this set of equations an Equation of State (EOS) linking the pressure P and the density ρ is needed [46]. In the user-SPH package of LAMMPS one EOS is used for the liquid (Tait's EOS) and one for gas phase (ideal gas EOS). In this section we will implement a new EOS for the liquid phase. Note that with similar steps is also possible to implement a new gas EOS.
Le Métayer and Saurel [47] combined the "Noble-Abel" and the "Stiffened-Gas" EOS proposing a new EOS called Noble-Abel Stiffened-Gas (NASG), suitable for multiphase flow. The expression of the EOS does not change with the phase considered. For each phases, the pressure and temperature are calculated as function of density and specific internal energy, e.g., and temperature-wise where P, ρ, e, and q are, respectively, the pressure, the density, the specific internal energy, and the heat bond of the corresponding phase. γ, P ∞ , q, and b are constant coefficients that defines the thermodynamic properties of the fluid. For this new pair style, called sph/nasgliquid, we take as a reference the sph/taitwater pair style declared and initialised in pair_sph_taitwater.h and pair_sph_taitwater.cpp files in the directory /src/USER-SPH. All the files regarding sph/nasgliquid must be saved in the /src/USER-SPH directory and its hierarchy is shown in Figure 4.. we will implement a new EOS for the liquid phase. Note that with similar steps is also possible to implement a new gas EOS. Le Métayer & Saurel [47] combined the "Noble-Abel" and the "Stiffened-Gas" EOS proposing a new EOS called Noble-Abel Stiffened-Gas (NASG), suitable for multiphase flow. The expression of the EOS does not change with the phase considered. For each phases, the pressure and temperature are calculated as function of density and specific internal energy, e.g., and temperature wise where P, ρ, e, and q are, respectively, the pressure, the density, the specific internal energy, and the heat bond of the corresponding phase. γ, P ∞ , q, and b are constant coefficients that defines the thermodynamic properties of the fluid. For this new pair style, called sph/nasgliquid, we take as a reference the sph/taitwater pair style declared and initialised in pair_sph_taitwater.h and pair_sph_taitwater.cpp files in the directory /src/USER-SPH. All the files regarding sph/nasgliquid must be saved in the /src/USER-SPH directory.

validation
The sph/nasgliquid pair style has validated by Albano & Alexiadis [26] to study the Rayleigh collapse of an empty cavity.

Validation
The sph/nasgliquid pair style has validated by Albano and Alexiadis [26] to study the Rayleigh collapse of an empty cavity.

pair_sph_nasgliquid.cpp
All the functions will be the same as in the reference sph/taitwater. However, in our new sph/nasgliquid, we need to substitute the "PairSPHTaitwater" text in "PairSPHNasgliquid", as can be seen in Listings 35 and 36.

pair_sph_nasgliquid.h
In the header of the new pair style we need to substitute the "PairSPHTaitwater" text in "PairSPHNasgliquid" as well as declare new protected members in the class, the pointers to the new arguments.

Multiphase (Liquid-Gas) Heat Exchange Pair Style
In LAMMPS thermal conductivity between SPH particles is enabled using the sph/heatconduction pair style inside the user-SPH package. However, the pair style is designed only for mono phase fluid where the thermal conductivities is constant (κ i = κ). When more than one phase is present, the heat conduction at the interface can be implemented by using [6,26] In the new pair style, called sph/heatgasliquid, one phase is assumed to be liquid with an initial temperature of T l,0 and the other is assumed to be and ideal gas. Each time-step the temperature of the fluid is updated as [26].
where T l,0 is the reference temperature, E 0 the internal energy in [J], E l internal energy [J] at the current time step and C p,l is heat capacity of the fluid in [J K −1 ]. The temperature of the gas is updated following the ideal EOS [26].
where MM is the molar mass [kg kmol −1 ], e g is the specific internal energy in [J kg −1 ], γ is the heat capacity ratio and R is the ideal gas constant in [J K −1 kmol −1 ]. Generally the choice of the reference states E l,0 is arbitrary, but if the Equation of State (EOS) used for the phase is function of both density and internal energy of the reference state will be determined by the EOS.
In the sph/heatgasliquid pair style is important to check if the i-th and j-th particles are liquid or gas phase to apply either Equation (10) or Equation (11). This "phase check" is explained in Section 5.2 compute function is modified.
For the energy balance the new pair style needs T l,0 , E l,0 , C p,l and κ l for the liquid phase and κ g for the gas phase. Moreover, for the phase check, the particle types of each phases must be specified. All this informations is passed by the user in the in the input file.
The reference pair style is sph/heatconduction. It is declared and initialised in the pair_sph_heatconduction.cpp pair_sph_heatconduction.cpp files in the directory /src/USER-SPH. All the files regarding sph/heatgasliquid must be saved in the /src/USER-SPH directory and its hierarchy is shown in Figure 5.
where MM is the molar mass [kg kmol -1 ], e g is the specific internal energy in [J kg -1 ], γ is the heat capacity ratio and R is the ideal gas constant in [J K -1 kmol -1 ]. Generally the choice of the reference states E l,0 is arbitrary, but if the equation of state (EOS) used for the phase is function of both density and internal energy of the reference state will be determined by the EOS.
In the sph/heatgasliquid pair style is important to check if the i-th and j-th particles are liquid or gas phase to apply either Equation (10) or Equation (11). This "phase check" is explained in Section 5.2 compute function is modified.
For the energy balance the new pair style needs T l,0 , E l,0 , C p,l and κ l for the liquid phase and κ g for the gas phase. Moreover, for the phase check, the particle types of each phases must be specified. All this informations is passed by the user in the in the input file.
The reference pair style is sph/heatconduction. It is declared and initialised in the pair_sph_heatconduction.cpp pair_sph_heatconduction.cpp files in the directory /src/USER-SPH. All the files regarding sph/heatgasliquid must be saved in the /src/USER-SPH directory.

Validation
The sph/heatgasliquid pair style has validated by Albano & Alexiadis [26] to study the role of the heat diffusion in for a gas filled Rayleigh collapse in water.

Validation
The sph/heatgasliquid pair style has validated by Albano and Alexiadis [26] to study the role of the heat diffusion in for a gas filled Rayleigh collapse in water.

pair_sph_heatgasliquid.cpp
All the functions will be the same as in the reference sph/heatconduction. However, in our new sph/heatgasliquid, we need to substitute the "PairSPHHeatConduction" text in "PairSPHHeatgasliquid", as can be seen in Listings 52 and 53.

pair_sph_heatgasliquid.h
In the header of the new pair style we need to substitute the "PairSPHHeatConduction" text in "PairSPHHeatgasliquid" and declare new protected members in the class.

Full Stationary Fix Style
In LAMMPS a fix style is any operation that is applied to the system, usually to a group of particles, during time stepping or minimisation used to alter some property of the sytem [41]. There are hundreds of fixes defined in LAMMPS and new ones can be added. Usually fixes are used for time integration, force constraints, boundary conditions and diagnostics.
In the user-sph package in LAMMPS there is the so called meso/stationary fix used to set boundary condition. With meso/stationary is possible to fix position and velocity for a group of particles, walls as example, but internal energy and density will be updated. In some cases, it is useful to have a fully stationary conditions that maintains constant also the energy and the density. For this new fix, called meso/fullstationary, we take as a reference the meso/stationary fix declared and initialised in fix_meso_stationary.h and fix_meso_stationary.cpp files in the directory /src/USER-SPH. All the files regarding meso/fullstationary must be saved in the /src/USER-SPH directory and its hierarchy is shown in Figure 6.

Full Stationary Fix Style
In LAMMPS a fix style is any operation that is applied to the system, usually to a group of particles, during time stepping or minimisation used to alter some property of the sytem [41]. There are hundreds of fixes defined in LAMMPS and new ones can be added. Usually fixes are used for time integration, force constraints, boundary conditions and diagnostics.
In the user-sph package in LAMMPS there is the so called meso/stationary fix used to set boundary condition. With meso/stationary is possible to fix position and velocity for a group of particles, walls as example, but internal energy and density will be updated. In some cases, it is useful to have a fully stationary conditions that maintains constant also the energy and the density. For this new fix, called meso/fullstationary, we take as a reference the meso/stationary fix declared and initialised in fix_meso_stationary.h and fix_meso_stationary.cpp files in the directory /src/USER-SPH. All the files regarding meso/fullstationary must be saved in the /src/USER-SPH directory.

LAMMPS
Modify Fix FixMesoFullStationary FixMesoStationary Figure 6. Class hierarchy of the new fix style

validation
The meso/fullstationary has been used in the validation of the new viscosity class to set the boundary condition of a constant asymmetric heated walls, see Section 7.2.

fix_meso_fullstationary.cpp
All the functions will be the same as in the reference meso/stationary. However, in our new fullstationary, we need to substitute the "FixMesoStationary" text in "FixMeso-FullStationary", as can be seen in Listing 69 and 70

Validation
The meso/fullstationary has been used in the validation of the new viscosity class to set the boundary condition of a constant asymmetric heated walls, see Section 7.2.

fix_meso_fullstationary.cpp
All the functions will be the same as in the reference meso/stationary. However, in our new fullstationary, we need to substitute the "FixMesoStationary" text in "FixMeso-FullStationary", as can be seen in Listings 69 and 70.

fix_mes_fullstationary.h
In the header of the new fix we need to substitute the "FixMesoStationary" text in "FixMesoFullStationary".

Viscosity Class
Viscosity in the SPH method has been addressed with different solutions [46]. Shock waves, for example, have been a challenge to model due to the arise of numerical oscillations around the shocked region. Monaghan solved this problem with the introduction of the so-called Monaghan artificial viscosity [48]. Artificial viscosity is still used nowadays for energy dissipation and to prevent unphysical penetration for particles approaching each other [25,49]. The SPH package of LAMMPS uses the following artificial viscosity expression [6], within the sph/idealgas and sph/taitwater pair style.
where α is the dimensionless dissipation factor, c i and c j are the speed of sound of particle i and j. The dissipation factor, α, can be linked with the real viscosity in term of [6] where c is the speed of sound, ρ the density, µ the dynamic viscosity and h the smoothing length. The artificial viscosity approach performs well at a high Reynolds number but better solutions are available for laminar flow: Morris et al. [50] approximated and implemented the viscosity momentum term for SPH. The same solution can be found in the sph/taitwater/morris pair style with the expression [6].
where µ is the real dynamic viscosity.
In LAMMPS both the dissipation factor and the dynamic viscosity are treated as a constant between a pair of particles when they interact within the smoothing length. In this section we want to make the viscosity a per atom property instead of a pair property only existing within a pair style. Moreover, five temperature dependent viscosity models are added. For this example, no reference file is used; a new class, Viscosity, is implemented in LAMMPS from scratch and its hierarchy is shown in Figure 7.
solutions are available for laminar flow: Morris et al. [50] approximated and implemented the viscosity momentum term for SPH. The same solution can be found in the sph/taitwater/morris pair style with the expression [6] ∑ j m i m j (µ i + µ j )v ij where µ is the real dynamic viscosity. In LAMMPS both the dissipation factor and the dynamic viscosity are treated as a constant between a pair of particles when they interact within the smoothing length. In this section we want to make the viscosity a per atom property instead of a pair property only existing within a pair style. Moreover, five temperature dependent viscosity models are added. For this example, no reference file is used; a new class, Viscosity, is implemented in LAMMPS from scratch.

Temperature dependant viscosity
In literature multiples empirical models that correlate viscosity with temperature are available [51][52][53]. In the new viscosity class five different viscosity models have been implemented:

1.
Andrade's equation [54] where µ is the viscosity in [Kg m -1 s -1 ], T is the static temperature in Kelvin, A, B, C and D are fluid-dependent dimensional coefficients available in literature.

Temperature Dependant Viscosity
In literature multiples empirical models that correlate viscosity with temperature are available [51][52][53]. In the new viscosity class five different viscosity models have been implemented:

1.
Andrade's equation [54] where µ is the viscosity in [Kg m −1 s −1 ], T is the static temperature in Kelvin, A, B, C and D are fluid-dependent dimensional coefficients available in literature.

2.
Arrhenius viscosity by Raman [55,56] where µ is the dynamic viscosity in [Kg m −1 s −1 ], T is the temperature in Kelvin, C 1 and C 2 are fluid-dependent dimensional coefficients available in literature.

3.
Sutherland's viscosity [57,58] for gas phase Sutherland's law can be expressed as: where µ is the viscosity in [Kg m −1 s −1 ], T is the static temperature in Kelvin, C 1 and C 2 are dimensional coefficients.

4.
Power-Law viscosity law [57] for gas phase A power-law viscosity law with two coefficients has the form : where µ is the viscosity in [Kg m −1 s −1 ], T is the static temperature in Kelvin, and B is a dimensional coefficient.

Constant viscosity
With constant viscosity both dissipation factor and dynamic viscosity will be constant during the simulation.
When the artificial viscosity is used the dissipation factor of Equation (12) is defined as the arithmetic mean of the dissipation factors of i-th particle and j-th particle.
where α ij is the dissipation factor of the particles pair i and j.

Validation
In order to validate the new Viscosity class, we will study the effect of asymmetrically heating walls in a channel flow, and more specifically the effect on the velocity field of the fluid. The data obtained with our model will be compared with the analytical solution obtained by Sameen and Govindarajan [59].
The water flows between two walls in the x-direction with periodic conditions. The walls are set at different temperatures T cold and T hot , see Figure 8. Both water and walls are modelled as fluid following the tait EOS. The physical properties of the walls are set constant throughout the simulation using the full stationary conditions described in Section 6. To match the condition used by Sameen and Govindarajan we set the cold wall temperature to T cold = 295 K and the temperature dependence of the dynamic viscosity described by the Arrhenius model, Equation (16), with C 1 = 0.000183 [Ns m −2 ] and C 2 = 1879.9 K [59]. To describe the asymmetric heating Sameen and Govindarajan introduced the parameter m, defined as: where µ re f = µ hot is the viscosity at the hot wall in the case of asymmetric heating and µ cold is the viscosity at the cold wall. By combining (16) and (20), with the given T cold , is possible to express the temperature difference of the walls ∆T as function of m. Figure 9 shows the viscosity trend for different values of m and the corresponding ∆T. Sometimes, in particle methods, instantaneous data can be noisy (scattered) as can be seen from the blue circles of both Figures 9 and 10. In all the cases considered, the model is in good agreement with the work of Sameen and Govindarajan. Figure 10 shows the dimensionless velocity trend for different values of m. Again, the model is in good agreement with the analytical solution of Sameen and Govindarajan always laying within the velocity scattered points. In both our model and in the analytical solution the maximum of the velocity shifts to the right as m increases. We can conclude that our model is in good agreement with the literature, showing the typical viscosity and velocity profiles for asymmetric heating confirming the correct functionality of the new viscosity class.

New Abstract Class: Viscosity
To implement the new viscosity model a new abstract class has been created, called Viscosity. The class has no attribute, and one virtual method: compute_visc, that is used to compute the viscosity using one of the Equations (15)- (18). As usual, the Viscosity class is divided in two files, see Listings 78 and 79. As it is an abstract class, it cannot be instantiated. It is used as a base, a mold, to implement the viscosity models. All implemented viscosity classes, such as the ones implementing the Arrhenius viscosity or the Sutherland viscosity, will inherit from this class.  This type of base class is called an interface, though as the code is written in C++, there is no actual difference in the implementation. The difference is only in concepts.
This structure allows for a very simple procedure to add a new viscosity type to LAMMPS, as one doesn't have to go through all of the code everytime a new viscosity type is implemented. All that is required is to implement a new viscosity class inheriting from the Viscosity abstract class and modify the add_viscosity function. The details of the changes required for those two actions are detailed later in this section.
Another structure one might think of to implement the viscosity abstract base class would be a template. Indeed, templates are more efficient than inherited classes as inherited classes create additional virtual calls when calling the class's methods. However, the choice of which viscosity should be called is made at runtime, and not at compile time, which means the abstract base class would be a better fit. When runtime polymorphism is needed, the structure preferred is an abstract base class.
The abstract class is not the most efficient implementation, but it allows for simplicity of use, which is important considering most of LAMMPS users are not programmers. In this work, we have chosen to sacrifice a bit of efficiency to gain ease of use.

Implementing a New Viscosity Class
In this section the steps to implement one of Equations (15)- (18) are shown, using the four parameter exponential viscosity as an example.
A new class is created that inherits from the Viscosity abstract class. The new class have as much attributes as the viscosity type has parameters. In this example that means four, as shown in the header in Listing 80.   6 We want to be able to choose which type of viscosity is being used in the simulation from the input file, using a new command called viscosity. Let's discuss the implementation of this feature. First we need to define the viscosity command. This is done by modifying the execute_command method of the Input class. We then define a new function called add_viscosity, whose declaration is shown in Listing 86 and definition in Listing 87. This function will have to be modified each time one wants to create a new viscosity class. In add_viscosity, the element arg[0] is the string representing the type of viscosity. For each viscosity class, the method performs the following procedure: Listing 88: New include (atom.cpp) 1 #include <string.h> 2 #include <iostream> 3 #include "viscosity_four_parameter_exp.h" 4 #include "viscosity_sutherland_law.h" 5 #include "viscosity_power_law_gas.h" 6 #include "viscosity_arrhenius.h" 7 #include "viscosity_constant.h" The viscosity attribute is initialised to NULL in the constructor, see Listing 89. In the destructor of the Atom class, we add a line to delete the viscosity attribute, see Listing 90.

memory->destroy(viscosity);
The extract function is modified to process the viscosity attribute, see Listing 91.

Using compute_Visc in SPH Pair Styles: Tait Water Implementation
The dynamic viscosity is used to compute the artificial viscosity force, that is used in the compute function of the following SPH pair style: sph/idealgas, sph/lj, sph/taitwater and sph/taitwater/morris. In this section the steps to implement compute_visc in sph/taitwater are shown, the others required a similar procedure.

}
Inside the compute function of the sph/taitwater pair style we need to declare a new set of variables. Where e is the energy and cv the heat capacity, now needed to calculate the temperature and thus the viscosity.
The dynamic viscosities µ i and µ j are calculated for each atoms, using the formula implemented in the compute_visc method. The temperature for the i-th atom is obtained using T i = e i /cv i . It is important to note that using such expression for the energy balance prevents the reference state of the internal energy to be set at 0.
The constant viscosity matrix element is replaced by the formula defined in Equation (19) Viscosity is now a per atom property, this means that we don't have to pass its value then the pair style is invoked. For this reason we need to delete the viscosity related lines inside coeff.    force->bounds(FLERR,arg [1], atom->ntypes, jlo, jhi); 12 13 double rho0_one = force->numeric(FLERR,arg [2]); 14 double soundspeed_one = force->numeric(FLERR,arg [3]); 15 double cut_one = force->numeric(FLERR,arg [4]); 16 double B_one = soundspeed_one*soundspeed_one*rho0_one/7;  The last modification is in init_one. Again, we delete lines related to the former viscosity attribute. At this stage, the software is designed to only run in serial. Changes need to be made to make it run with Message Passing Interface (MPI). This will allow the software to run in parallel: some computations being independent from each other, they can be performed at the same time. Instead of using one processor for a long time, we will use multiple processors for a shorter period. The simulation will therefore take more computing resources but will take a lot shorter to compute. The original SPH module can already be run with MPI however as we have modified the code that is no longer true. We need to make additional changes to the software. All those changes are located in the Atom Vec Meso class of the SPH module.
In LAMMPS, the different MPI processes have to communicate with each other as the computations they perform are not completely independent from each other. They need data from other processes in order to perform their own calculations. They communicate with each other using a buffer that will contain all the necessary data. The buffer is simply an array that we will fill with the data. The different methods for packing and unpacking this buffer are defined in the Atom Vec Meso class. We need to add a new data to transmit: the calculated viscosity.
The first thing to do is to increase the size of the buffers in their initialisation so they can accept the viscosity value, an example is shown in Listings 103 and 104.  9 // we also communicate de and drho in reverse direction 10 comm_f_only = 0; 11 // 3 + rho + e + vest [3], that means we may 12 // only communicate 5 in hybrid   After making those changes for all the methods in the class, the software can be run using mpirun.

Invoking, Selecting and Computing a Viscosity Object
To compute the new viscosity a new argument was added to the compute command: viscosities. This allows the user to use the compute command to output the dynamic viscosity to the dump file. This can be done by the following command: 1 compute viscosities_peratom all meso/viscosities/atom The implementation of this feature is simple, as it is very similar to other compute argument implementation. All that needs to be done is to modify another compute's implementation, such as compute_meso_rho_atom so it processes the variable viscosities instead of rho.
The viscosity used in the simulation can be invoked in the input file, using the following command: For example, to invoke the four parameter exponential viscosity, we can write in the input file: 1 viscosity FourParameterExp C1 C2 C3 C4 As stated earlier, this list can easily be extended by the user by modifying the add_viscosity function defined earlier.

Conclusions
Particle methods are very versatile and can be applied in a variety of applications, ranging from modelling of molecules to the simulation of galaxies. Their power is even amplified when they are coupled together within a discrete multiphysics framework. This versatility matches well with LAMMPS, which is a particle simulator, whose open-source code can be extended with new functionalities. However, modifying LAMMPS can be challenging for researchers with little coding experience and the available support material on how to modify LAMMPS is either too basic or too advanced for the average researcher. Moreover, most of the available material focuses on MD; while the aim of this paper is to support researchers that use other particle methods such as SPH or DEM.
In this work, we present several examples, explained step-by-step and with increasing level of complexity. We begin with simple cases and concluding with more complex ones: Section 3 shows the implementation of the Kelvin-Voigt bond style used to model encapsulate particles with a soft outer shell and validated validated by simulating spherical homogeneous linear elastic and viscoelastic particles [45]; Section 7 show how to implement a new per-atom temperature dependant viscosity property and is validated finding the same viscosity and velocity trend shown by Sameen and Govindarajan [59] in their analytical solution for a channel flow in a asymmetrical heating walls.
The work perfectly fits in the "Discrete Multiphysics: Modelling Complex Systems with Particle Methods" special issue by sharing some in dept know-how and "trick and trades" developed by our group in years of use of LAMMPS. In fact, the aim is to support, in several ways, researchers that use computational particle methods. Often researchers tend to write their own code. The advantage of this approach is that the code is well understood by the researcher and, therefore, easily extendible. However, this sometimes implies reinventing the wheel and countless hours of debugging. Familiarity with a code like LAMMPS, which has an active community of practice and is periodically enriched with new features would be beneficial to this type of researchers allowing them to save considerable time. In the long term, there is another advantage. Modules written for in-house code are hardly sharable. At the moment, the largest portion of the LAMMPS community is dedicated to MD. While this article was under review, for instance, a new book dedicated to modifying LAMMPS came out [60]. However, it focuses only on MD and it does not mention other discrete methods like SPH or DEM. Instead, the aim of this paper is to make LAMMPS more accessible for the Discrete Multiphysics community facilitating sharing reusable code among practitioners in this field.  The crucial line for DMP simulations is the hybrid keyword of the atom_style, which allows for combining different particle models. The keyword meso refers to the SPH model and bond, in the case under consideration, to the LSM. The angle keyword corresponds to angular springs, but, as it will be clear later, it is not used in this simulation.
The following section contains several variables that are going to be used later on. In particular, the initial particle distance is dL and their mass m. The resolution of the membrane is Nt times higher than the fluid. The initial distance between membrane particle is therefore db=dL/Nt and their mass mM=m/Nt.