Molecular Dynamics Simulations of Active Matter using LAMMPS
MMolecular Dynamics Simulations of Active Matter using LAMMPS
C. S. Dias
1, 2, ∗ Departamento de F´ısica, Faculdade de Ciˆencias,Universidade de Lisboa, 1749-016 Lisboa, Portugal Centro de F´ısica Te´orica e Computacional, Faculdade de Ciˆencias,Universidade de Lisboa, 1749-016 Lisboa, Portugal
LAMMPS is a widely popular classical Molecular Dynamics package. It was designed for materialsmodeling but it is well prepared for simulations in Soft Matter. The use packages like LAMMPShas advantages and disadvantages. The main advantage is the optimization of the methods, mainlyfor parallel computing. The main disadvantage is that, due to the complexity of the code, it ha longlearning curve. One purpose of these notes is to shorten that curve for researchers are starting to useLAMMPS to simulate soft active matter. In these notes, we first discuss some Molecular Dynamicsmethods implemented in LAMMPS. We present an hands-on introduction to first-time users and wefinish with an advanced hands-on section, where we implement and test Active Brownian particlessimulations.
Keywords : LAMMPS, Molecular Dynamics, Langevin Dynamics, Active Matter, Active Brownianparticles
I. INTRODUCTION
LAMMPS is an acronym for Large-scaleAtomic/Molecular Massively Parallel Simulator. Itis a classical Molecular Dynamics package optimizefor parallelization and mainly directed for materialsmodeling. Its development began in the mid 1990s, withits first versions in FORTRAN 77 and 90. It was buildunder a cooperative research & development agreement(CRADA) between two DOE labs (Sandia and LLNL)and 3 companies (Cray, Bristol Myers Squibb, andDupont). The coding effort was led by Steve Plimptonat Sandia [1–6].After the initial development in FORTRAN, the cur-rent versions are written in C++ and it was released asan open source code in 2004. The package is freely avail-able for download under GPL, and is designed to be easyto modify or extend. The package is well documented(see https://lammps.sandia.gov ).LAMMPS is a very powerful package for Molecular Dy-namics simulations, and the use of this type of packagecan have many advantages and disadvantages. The ad-vantages are clear, since it is designed to be modular,it is relatively easy to implement a new model and isoptimize for parallel computing. The use of this typeof packages gives us a confidence in the code since it istested by many. In the other hand, it has also a fewdisadvantages. For instance, also because it is modular,the implemented routines are thought-out to be generic,which lead to a lower performance than methods imple-mented for specific models. One main disadvantage isalso the complexity of the code, since it is difficult tofollow everything implemented in the code, it has a longinitial learning curve. The objective of these notes is toshorten that initial learning curve for researchers startingto simulate soft active matter with LAMMPS. ∗ [email protected] These notes are divided into three main sections. Thefirst section discusses basic Molecular Dynamics meth-ods behind LAMMPS and their implementation. Thesecond section is a hands-on introduction to LAMMPS.The third section is a more advanced hands-on sectionwhere we implement a custom method for the simulationof Active Brownian particles.
II. BASIC MOLECULAR DYNAMICS INLAMMPS
In this section we discuss the basic Molecular Dynam-ics methods implemented in LAMMPS. Since MolecularDynamics simulations consist in integrating the equa-tions of Newton for many particles [7], we start by dis-cussing the integration method. This method needs theforce computation, that is discussed in the following sec-tion. Some concepts of parallelization are discussed next.We then show how Langevin Dynamics simulations areperformed in LAMMPS, and the last section has a briefdiscussion on reduced units.
A. Velocity Verlet integration
Many different methods can be used to integrate theequations of Newton. Here we discuss the velocity Verletmethod. This method is similar to the classical Verletmethod and consists of two integration steps [8]: (cid:126)r ( t + ∆ t ) = (cid:126)r ( t ) + (cid:126)v ( t )∆ t + 12 (cid:126)a ( t )∆ t , (1)and, (cid:126)v ( t + ∆ t ) = (cid:126)v ( t ) + (cid:126)a ( t ) + (cid:126)a ( t + ∆ t )2 ∆ t, (2)where (cid:126)r is the position, (cid:126)v the velocity, and (cid:126)a the acceler-ation. ∆ t is the integration timestep. One advantage of a r X i v : . [ c ond - m a t . s o f t ] F e b the velocity Verlet method when compared to the classi-cal Verlet one is that it deals directly with the velocity,which avoid problems in the definition of the initial ve-locity,. LAMMPS implementation:
One could implementdirectly the two steps given by Eqs. 1 and 2, howeverthese steps need to save, at all time, two acceleration (orforce) vectors, which take up memory and communica-tion time between parallel jobs. The most common im-plementation, and the one used in LAMMPS, computesthe velocity twice, using an intermediate half timestepjump forward. The steps are the following:1. Compute (cid:126)v ( t + ∆ t/
2) = (cid:126)v ( t ) + (cid:126)a ( t ) ∆ t ;2. Compute (cid:126)r ( t + ∆ t ) = (cid:126)r ( t ) + (cid:126)v ( t + ∆ t/ t ;3. Compute the forces at (cid:126)r ( t + ∆ t );4. Compute (cid:126)v ( t + ∆ t ) = (cid:126)v ( t + ∆ t/
2) + (cid:126)a ( t + ∆ t )∆ t .With this implementation, only one vector of forces isused since they computed between the two steps. B. Force computation
The force calculation is usually the bottleneck for anefficient Molecular Dynamics simulation. For the sim-plest case of pair potentials, for instance, the calculationis done over all pairs of particles. One first optimizationis by considering that the potential interaction betweenparticles can be neglected above a certain distance, r cut .Taking that into account, pairs of particles at distancesabove r cut do not need to be computed. Different meth-ods can be used to take advantage of that. The mostcommon ones are the Verlet list (see Fig. 1(a)) and thecell list (see Fig. 1(b)) methods. The first consists on sav-ing a list of neighbor particles and update the list fromtime to time. The second is to discretize the system intocells and only check particles in the neighboring cells. LAMMPS implementation:
In LAMMPS, bothmethods are implemented.The calculation of the force isdone using the Verlet list (see Fig. 1(a)). The list is de-fined as all particles that are at a distance of r cut + skin ,where the r cut depends on the potential and the skin isdefined (command “neighbor” in the LAMMPS script).The cell list method (see Fig. 1(b)) is then used to updatethe particles in the Verlet list, otherwise it would be nec-essary to run all pairs of particles. The length of the cellcan be defined and it is a trade-off between the time ittakes to run all cells if the cells are too small, or the num-ber of unnecessary particles check if they are two large.LAMMPS consider the length of the cell r cut / FIG. 1. (a) Neighbor list (Verlet list) used to choose the pairsof particles to compute the force. r cut is the cutoff distanceof the potential and skin is an extra distance added to r cut as a threshold to build the list of neighbors. (b) Cell listmethod where the system is divided into cells. In LAMMPSit is not used to compute the force but to build the neighborlist without searching all particles.FIG. 2. Spacial domain decomposition by processes 1 to 4.Particles are allocated to each process and particles in thered region that belong to a process are ghost particles in theprocess of the adjacent domain. C. Parallel computing
One of the main advantages of Molecular Dynamicssimulations is its scalability when split into multiple pro-cesses. Basically, most tasks can be done in parallel,needing only information from the positions of other par-ticles during the force calculation. One of the most com-mon method for parallel computing is the spatial decom-position method [6]. As you can see in Fig. 2, the spaceis divided into domains that are assigned to different pro-cesses. Particles inside the red region that belong to acertain process are also ghost particles to processes in ad-jacent domains. This way, when the force is computed,all particles in the domain plus the ghost particles areused.
LAMMPS implementation:
In LAMMPS, thelength of the red region is given by r cut + skin (thecommand “comm modify” in the LAMMPS script canbe used to change the length of the red region). Themapping of the spacial domains can be controlled inLAMMPS (command “processors” in the LAMMPSscript), however, by default, LAMMPS optimize the do-main distribution in such a way to minimize the red re-gion volume. In the Molecular Dynamics pipeline oneneeds to add a step where the communication betweenprocesses is done and the particles that belong to eachprocess are updated. In LAMMPS, this communica-tion step occurs after the first integration of the veloc-ity Verlet method (before step 3 of force computation inSec. II A). D. Langevin dynamics
Soft matter dynamics usually occur inside viscous en-vironments (fluid), so to perform Molecular Dynamicssimulations one needs to add the interaction with thefluid. It is possible to add the particles of the fluid andperform brute force simulations, however the length andtime scales will be very short. One approximation isto coarse-grain the interaction with the fluid using theLangevin dynamics approach. Langevin dynamics fol-lows the Langevin equation where two new terms of forceare added: viscous drag and thermal noise. The Langevinequation for the translational degrees of freedom, in onedimension, is given by, m ˙ v ( t ) = −∇ x V ( x ) − γ t v ( t ) + ξ t ( t ) , (3)where m is the mass of colloids, γ t the translationaldamping coefficient related to a drag force exerted bythe medium, V the external potential, and ξ t ( t ) is thestochastic term, from the thermal noise, given from adistribution of zero mean and satisfying the relation, < ξ t ( t ) ξ t ( t (cid:48) ) > = 2 k B T γ t δ ( t − t (cid:48) ) , (4)where δ is the Dirac delta, T is the temperature, k B theBoltzmann constant. For a rigid body, one also needs theLangevin equation for the rotational degrees of freedom.The Langevin equation for the one-dimensional diffusionof a single-axis rotation is given by, I ˙ w ( t ) = −∇ θ V ( θ ) − γ r w ( t ) + ξ r ( t ) , (5)where I is the moment of inertia and γ r the rotationaldamping coefficient along that single-axis rotation, ξ r ( t )is the stochastic torque from the thermal noise, givenfrom a distribution of zero mean and satisfying the samerelation of, < ξ r ( t ) ξ r ( t (cid:48) ) > = 2 k B T γ r δ ( t − t (cid:48) ) . (6)For both translational and rotational motion, above adamping time (time for the inertia to become negligible),the Brownian diffusive regime is recovered. We obtainthe translational diffusion coefficient given by the Stokes-Einstein relation as [9], D t = k B Tγ t , (7) and the rotational diffusion coefficient as, D r = k B Tγ r . (8)The damping times for the translational and rotationalmotions are τ t = m/γ t and τ r = I/γ r , respectively.The translational and rotational diffusion coefficientspresented here are generic for passive or active particle.In the case of active particles, the rotational diffusion co-efficient sets the time scale for the self-propulsion direc-tional motion to forget its initial orientation (for activeparticles the effective diffusion coefficient depends also onthe propulsion which is discussed in Sec.IV). For passiveparticles, one can find a relation between translationaland rotational diffusion. For spherical passive particles,these coefficients take the form [9], D t = k B T πηR , (9)and, D r = k B T πηR , (10)where η is the viscosity of the fluid medium and R theradius of the particles. It is possible to obtain a relationbetween the rotational and translational diffusion coeffi-cient given by, D r D t = 34 R . (11) LAMMPS implementation:
The variables, inLAMMPS, that parameterize the Langevin dynamics arethe translational damping time, τ t and the mass m (com-mand “fix langevin” in the LAMMPS script). Since onlythe translational damping time is defined, an isotropicscale factor α = τ t /τ r , between translational and rota-tional damping, is introduced. For spherical particles,using Eq. 11, one can write the scale factor as α = 10 / α is an external parameter(“angmom” in the langevin fix), and for rigid more com-plex bodies, if they are approximately spherical you needto add the factor of 10/3 directly into the code [10], orfor more complex shapes the scale factor is one [11]. InLAMMPS, the Langevin equation for the translationalmotion, in one dimension, is given by, m ˙ v ( t ) = −∇ r V ( r ) − mτ t v ( t ) + (cid:114) mk B Tτ t ξ ( t ) , (12)and for the rotational motion, around one single axis ofrotation, is, I ˙ ω ( t ) = −∇ θ V ( θ ) − αIτ t ω ( t ) + (cid:114) αIk B Tτ t ξ ( t ) . (13)In eqs. 12 and 13, v and ω are the translational and angu-lar velocity, m and I are mass and Inertia of the colloidalparticle, and V is the external potential. In LAMMPS,for efficiency [12], ξ ( t ) is an uniform distribution from − . .
5, which forces a correction in the stochasticterm of the Langevin equation since the second momentof the uniform distribution is given by, σ = (cid:90) . − . x dx = 23 (cid:18) (cid:19) = 224 , (14)and imposes a correction of (cid:112) /
24, that can be seen inthe source code in the file fix langevin.cpp. In the samesource file, the integration timestep is included in thedenominator inside the square root because LAMMPSmultiply all forces (including the stochastic term fromthe thermal noise) by the timestep in the velocity Verletintegration (see Sec. II A). The timestep in the stochasticterm comes from the variance of the random number dis-tribution which is σ = 2 D ∆ t (which is the multiplicativefactor in the stochastic term of the Langevin equation).This way, the diffusion coefficient does not depend on thetimestep. E. Reduced Units
When modeling Soft Matter, we are usually workingwith length scales in the order of nanometers to microm-eters and energy scales in the order of 10 − to 10 − eV. Ifwe decide to use standard SI units, we will have to workwith very small values. It is, of course, more convenientto work with values near unity, which can be done byusing the size of the simulation elementary particles asthe reference. Reduced units are not only used for con-venience. Performing numerical calculations with verysmall values can lead to rounding errors due to the ma-chine precision.Using dimensionless units also facilitatesscaling, where a single model can help to understand alarge class of problems at different scales [13]. LAMMPS implementation:
A common scale forreduced units is the one of the Lennard-Jones potential[13], given by, V LJ ( r ) = 4 (cid:15) (cid:20)(cid:16) σr (cid:17) − (cid:16) σr (cid:17) (cid:21) , (15)where r is the distance between two particles, (cid:15) the depthof the potential well, and σ the width of the potential(distance at which the potential is zero). LAMMPS usesthese units and, without loss of generality, sets the funda-mental quantities m , σ , (cid:15) , and k B to unity. The masses,distances, energies are multiples of these fundamentalvalues. The formulas relating the reduced or unitlessquantity are show in Table I.The conversion to real units can be sometimes con-fusing. Nonetheless, it is always better to use reducedunits instead of SI units (which LAMMPS allows). Theconversion between real units and reduced units can bedone by using the relations in Table I. For instance, Soft TABLE I. LAMMPS reduced units.Variable ValueTemperature ( T ) (cid:15)k B Time ( τ ) (cid:113) mσ (cid:15) Boltzmann constant ( k B ) 1Energy ( (cid:15) ) 1Distance ( σ ) 1Mass ( m ) 1Translational damping time ( τ t ) (cid:113) mσ (cid:15) Matter experiments are usually performed at room tem-perature of T = 293 K . So taking Boltzmann constant, k B = 1 . × − m kgs − K − , we obtain the energyunit of (cid:15) = 4 . × − J . The same can be done with thelength units, for instance, if we set the mass of a colloidas M = 10 m and the diameter as d = 10 σ in the simula-tions, and consider the real values of colloid mass 10 − g and radius 0 . µ m, we recover the fundamental quantitiesof σ = 10 − m and m = 10 − kg respectively. The samecan be done with the time units. If we set a timestep inthe simulations of ∆ t = 0 .
01, this tells us that each sim-ulation step takes 1 . × − s . For Langevin Dynamicssimulations it is sometimes clearer to use the Browniantime [14] (time for a particle to diffuse the square of itsdiameter) to set the timescale. III. LAMMPS: BASIC HANDS-ON
In this section we will have a more practical approachto LAMMPS. We will start by downloading and com-pile LAMMPS, run a test example, and performing somemeasurements on a Langevin Dynamics simulation ofpassive particles.
A. Compiling LAMMPS
LAMMPS was primarily developed for Linux systems,mainly because most of the high-performance computersare Linux based. With this in mind, if you have a Linuxsystem it is straightforward to download the most re-cent tarball ( https://lammps.sandia.gov/download.html ), uncompress and compile performing “make mpi”or “make serial” in the source folder (src/). If you havea windows 10 or newer I recommend to use the WindowsSubsystem for Linux (WSL). Clear instructions are avail-able at https://lammps.sandia.gov/doc/Howto_wsl.html . For older versions of windows I recommend youto install some Linux virtual machine. MacOS users caninstall like in the Linux case.In LAMMPS one can find extra packages to includesome specific features. These packages were added bythe LAMMPS developers or by other users, which areconsistent with the style of the rest of LAMMPS. Tocheck what packages are installed just run the command“make ps” in the source folder. In this course we will use,for instance, the ASPHERE package and the COLLOIDpackage. You can add them by running the commands“make yes-asphere” and “make yes-colloid”. You need torecompile after adding a new package.
B. Running LAMMPS
When compiling the code with the “make mpi” or“make serial”, an executable file is created in the sourcefolder with name lmp mpi or lmp serial. To run, one justneed to use the command: ./lmp_mpi -in lj.in to run in one core, or mpirun -np 4 ./lmp_mpi -in lj.in to run in four cores (the lmp serial executable do not al-low parallelization). The lj.in file is a LAMMPS scriptwhere the simulation parameters and methods are in-cluded. An example of the script for Lennard-Jones in-teracting particles in a box with periodic boundary con-ditions is the following: units lj atom_style atomic dimension 3 boundary p p p region box block 0 30 0 30 0 30 units box create_box 1 box lattice sc 0.1 create_atoms 1 box mass 1 1.0 velocity all create 1 176217 dist gaussian pair_style lj/cut 2.5 pair_coeff 1 1 1 1 dump DUMPXYZ all xyz 5 out.xyz timestep 0.005 fix 1 all nve thermo_style custom step temp pe ke etotal thermo 100 run 5000 Lines 1 and 2 define main properties of the simula-tions like the reduced (LJ) units and the atom style“atomic”defines particles as point particles. The lines 4to 7 define and create the simulation box in three dimen-sions with periodic boundary conditions in all directions
FIG. 3. (a) Initial and (b) final configuration of the NVE sim-ulation using the test script. The visualization was performedusing Ovito. (boundary p p p). A region named “box” is defined inline 6 and the simulation box with the properties of thatregion is created. Lines 9 to 10 create the particles in theregion “box”. They are distributed on a square cube (sc)lattice with a reduced density of particles of ρ ∗ = 0 . ρ ∗ = N/V with N the number of particles and V = L x × L y × L z the volume in reduced units. Thecommand “mass” sets the mass of a particle of type 1 as1 and the command “velocity” sets their velocities witha Gaussian distribution around 0 and a variance equalto the temperature 1. The value 176217 is the randomnumber generator seed.In line 14 and 15 the Lennard-Jones potential is de-fined, with a cutoff radius of 2.5 and coefficients of σ = 1and (cid:15) = 1 (two last values). The first two values of the co-efficients indicates which pair of particles type will havethose coefficients. Different particles could have differentinteractions defined. In 17 a xyz file is created with thepositions of all the particles every 5 time steps. Thesetype of files can be read in different external software,which I recommend Ovito ( ) tobe simple to download and use, and very powerful. InFig. 3, the initial and final configuration of a simulation,using the proposed script, are shown.In line 19 we define the simulation timestep. One cantest larger timesteps to check LAMMPS error outputswhen the simulation is unstable. In line 20, the fix nveis applied to all particles and perform integration of theequations of motion in constant number of atoms (N),volume (V), and energy (E). The “fix” commands inLAMMPS are used to change something in the integra-tion pipeline. In lines 22 and 23 the output informationis defined (if no file is defined the outputs are saved inthe lammps log file). In this case: timestep, tempera-ture, potential energy, kinetic energy, and total energy.One simple test is to check the constant total energy ofthe NVE integration. Line 25 defines the total numberof timesteps in the simulation. C. Langevin Dynamics of passive particles inLAMMPS
To test the concepts introduced in Sec. II D, we per-form Langevin Dynamics simulations in LAMMPS usingthe following script: units lj atom_style sphere variable temperature equal 1 variable damping equal 1 dimension 3 boundary p p p read_data inp.data dump DUMPXYZ all xyz 50 out.xyz timestep 0.005 fix 1 all nve fix 2 all langevin ${temperature} & ${temperature} ${damping} 82763871 compute msd_col all msd thermo_style custom step c_msd_col[1] & c_msd_col[2] c_msd_col[3] c_msd_col[4] ke thermo 100 run 50000 Some commands are the same as the previous scriptand the only relevant to Langevin Dynamics is the com-mand in line 16 where the Langevin fix is added to theNVE integration. In this fix, basically, extra terms areadded to the forces and torque following Eqs. 12 and 13.We add the starting and stopping temperature (which,in this case, we defined constant over the full simulationrun), and the Damping time introduced in Sec. II D. Thelast value is the random number generator seed.In this script we introduced a new atom style in line2, which is the sphere type, where particles are consid-ered as spheres, with a specific diameter and density. Wealso introduced the concept of variables in line 4 and 5which we use in the Langevin fix. One major differencein this script is in line 10 where we read the particles po-sitions from an external file (inp.data) with the followingstructure: FIG. 4. Mean square displacement as a function of time fora single particle and for 10000 particles. The inset shows thetrajectory of a single particle. Atoms
We define the number of particles and particle typesin line 3 and 5, and the size of the simulation box inlines 7, 8, and 9. The particles positions are definedin the Atoms sections, where we can see the first (andonly) particle (column 1 gives then particle number) oftype one (column 2), with a diameter of one (column3), a particle density of 1.9099 (column 4) and positionsin x, y, and z (columns 5, 6, and 7). The order of theparticles properties is specific for the sphere atom style,the manual of the read data command has examples forother atom styles. The density was chosen such that aparticle of unit diameter has unit mass.The final new command introduced in the run scriptis in line 19. The command “compute” is used byLAMMPS to perform calculations during the simulation.Here, we compute the mean square displacement (msd)of all particles (only one in this data file). The nameof the compute is msd col, which is used in the output“thermo style” (line 21) , which is a vector with four co-ordinates (x, y, z, and total), summed and averaged overall particles.In Fig. 4 we plot the output of the proposed script fora single particle and for 10000 particles. To implementthe many particles case, just add 10000 consecutive linesin the inp.data file. Since no pair potential has been de-fined one can just define the same position for all parti-cles. The proper slope of the mean square displacement,in three dimensions, of 6 D is recovered, where D is thediffusion coefficient of spherical passive particles given bythe Stokes-Einstein relation of Eq. 9. IV. IMPLEMENTATION OF ACTIVEBROWNIAN PARTICLES IN LAMMPS
As stated before, LAMMPS is an open source codeunder GPL. This means you are free to download, use,study, and even change. At this point, more than 95%of LAMMPS is made of add-on files to the core ones.Despite being a very complete and powerful tool, somespecific cases are not implemented and a user could needto add new features.One feature, relevant for our topic of active matter, isthe inclusion of Active Brownian particles simulations,i.e., particles that have a persistent motion in a specificdirection. This feature has been recently implemented inLAMMPS (see command “fix propel/self”), and othershave used custom implementations [15]. However, in thissection we discuss how to build a new class for ActiveBrownian particles, which can be used as an example onhow to implement custom methods into LAMMPS. Theimplementation of Active Brownian particles consists onadding a force (cid:126)F A = F(cid:126)e ( t ) [16, 17] where F is the mag-nitude of the activity (a constant) and (cid:126)e ( t ) is the unitorientation vector. The orientation vector evolution overtime is given by Eq. 13.Most of the features to be added to LAMMPS requirewriting a new C++ class. The simplest way is to findsome class that does something similar to the one wewant to add. For instance, in the case of an active par-ticle, one just need to add a force to a specific preferreddirection of the particle. Since this is not a pair interac-tion force, we should not build a new pair class. In thiscase, we just add a force to the ones already computed inthe pair classes much like is done in the Langevin dynam-ics case. The simplest classes one could use as examplesare the fix addforce or the fix setforce, which basicallyadd a value or set a value to the force array in a spe-cific function of the class, that LAMMPS uses for thosecases, called post force. The only thing we need to worryis how to define the orientation of the active particles.One option is to use ellipsoid particles (which we candefine all axis of equal size if we want spheres), whichthen integrates the rotational motion through the use ofquaternions [13]. In the case of ellipsoids, LAMMPS de-fines the x-axis as the reference axis and the rotation isdefined by the quaternions values using that same axisas a reference. You just need to have a fail check on yournew build class that your atom style is ellipsoid.For a fast implementation of active particles inLAMMPS one could do a simpler task. This task canbe resumed into four steps: • Copy the fix langevin class to a fix langevin activeone (each class has two files, in this casefix langevin.cpp and fix langevin.h). You just needto change the class name inside the files, and infix langevin active.h change the call of the fix from “langevin” to “langevin/active” for instance. Thisway, we use a fix, which is a bit more complex tofully understand than the two examples used be-fore, but is far easier to implement. For starter,the use of ellipsoids is already contemplated. • In the constructor function, in the beginning of theclass, we can access the arguments of the fix. Anoption for ellipsoids is the ”angmom” (which is thescale factor α of the rotational diffusion in Eq. 13),which has a numerical value after the option. Weonly need to add a second value with the force (justcopy the line with the variable “ascale” and addyour new variable from “arg[iarg+2]”) that repre-sents the activity, which will be saved into a vari-able defined in the header file. • The final step to have a fully functioning class withactivity is to add, in the angmom thermostat()function, a few lines inside the cycle that runs overthe local particles (particles in the same process)which will add a force to the forces vector, in thedirection given by the ellipsoid quaternions. Onecan convert quaternions to the three main axis ofthe ellipsoids using a function in MathExtra calledq to exyz (use only the principal axis of the ellip-soid).After implementing this new fix, we just need to recom-pile LAMMPS with the new files in the source folder. Wethen run a very similar script as the one for the Langevindynamics, units lj atom_style ellipsoid variable t equal 0.1 variable d equal 0.1 dimension 3 boundary p p p read_data inp.data compute shape all property/atom shapex shapey & shapez compute quat all property/atom quatw quati & quatj quatk dump positions all custom 50 out.dump & id type x y z & c_quat[1] c_quat[2] c_quat[3] c_quat[4] & c_shape[1] c_shape[2] c_shape[3] timestep 0.05 fix 1 all nve/asphere fix 2 all langevin/active $t $t $d & compute msd_col all msd thermo_style custom step c_msd_col[1] & c_msd_col[2] c_msd_col[3] & c_msd_col[4] ke thermo 100 run 50000 Notice only the change in the atom style and the langevinfix (and a different dump with shape and orientation totest the ellipsoids proper dynamics). In this new fixlangevin/active, after the angmom option we add theangmom value and the activity (100 in this script). Youcan use the angmom to control the rotational diffusioneffect, which in this case is set to 10 /
3, the value forspherical passive particles. Since we are dealing with el-lipsoids, the inp.data file changes to, Atoms
Ellipsoids
Here, the third column in the Atoms section, is not thediameter, since the shape is defined in collumns 2, 3, and4 of the Ellipsoids section, but a flag indicating that thisparticles is an ellipsoid. In the Ellipsoids section, theshape is defined by the diameter in x, y, and z (columns2, 3, and 4), which is a sphere in this example, and thenthe four last values are the quaternions quatw, quatx,quaty, and quatz.One simple test to this new active particles class isto measure the velocity at zero temperature for differentvalues of the force. The particles rapidly reach the ter-minal velocity and, using Eq. 12 at zero temperature, weobtain the relation v = τ t m F . This relation is shown inFig. 5.Another test is to measure the mean square displace-ment of an active particle. In Fig. 6, we can see thetwo typical regimes of an active particle, with a ballisticmotion at short times and a diffusive motion for longer FIG. 5. Terminal velocity as a function of the force, for activeparticles at zero temperature.FIG. 6. Mean square displacement of an active particle as afunction of time., performed for 10000 non-interacting activeparticles. The two regimes of ballistic and diffusive are seen,where the diffusive regime is linear and proportional to 6 D eff . times. The diffusive motion has an effective diffusion co-efficient given by [18], D eff = D t + 16 v D r , (16)where D t is the translational diffusion coefficient, v thevelocity, which we shown in Fig. 5 to be τ t /m , and D r therotational diffusion coefficient. D r can be related with D t using Eqs. 7 and 8, and the scaling factor α = τ t /τ r ,which gives the relation, D t D r = α Im , (17)and, for spherical passive particles, recovers Eq. 11. V. FINAL DISCUSSION
The LAMMPS package is a very powerful tool.However, like when implementing Molecular Dynamicsfrom scratch, one needs to be very careful with theparametrization of our models. We always need to startwith fewer particles and small systems, without paral-lelization. Then, as done in these notes, perform simplebut effective tests to check if we recover the expectedphysics.
VI. ACKNOWLEDGMENTS
We acknowledge financial support from the PortugueseFoundation for Science and Technology (FCT) under Contracts no. PTDC/FIS-MAC/28146/2017 (LISBOA-01-0145-FEDER-028146), CEECIND/00586/2017,UIDB/00618/2020, and UIDP/00618/2020. The noteswere part of the ”Initial Training on Numerical Methodsof Active Matter” organized by MSCA-ITN Active-Matter (Grant: 812780) funded by the EU H2020program. [1] P. J. in ’t Veld, S. J. Plimpton, and G. S. Grest, Comp.Phys. Commun. , 320 (2008).[2] M. L. Parks, R. B. Lehoucq, S. J. Plimpton, and S. A.Silling, Comp. Phys. Commun. , 777 (2008).[3] W. M. Brown, P. Wang, S. J. Plimpton, and A. N. Thar-rington, Comp. Phys. Commun. , 898 (2011).[4] W. M. Brown, T. D. Nguyen, M. Fuentes-Cabrera, J. D.Fowlkes, P. D. Rack, M. Berger, and A. S. Bland, Proc.Comp. Sci. , 186 (2012).[5] M. K. Petersen, J. B. Lechman, S. J. Plimpton, G. S.Grest, P. J. In ’T Veld, and P. R. Schunk, J. Chem.Phys. , 174106 (2010).[6] S. Plimpton, J. Comp. Phys. , 1 (1995).[7] M. P. Allen and D. J. Tildesley, Computer simulation ofliquids (Oxford University Press, Oxford, 1987).[8] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R.Wilson, J. Chem. Phys. , 637 (1982).[9] M. G. Mazza, N. Giovambattista, H. E. Stanley, andF. W. Starr, Phys. Rev. E , 031203 (2007). [10] C. S. Dias, C. Braga, N. A. M. Ara´ujo, and M. M. Teloda Gama, Soft Matt. , 1550 (2016).[11] H. P. Melo, C. S. Dias, and N. A. Ara´ujo, Commun.Phys. , 154 (2020).[12] B. Dunweg and P. Wolfgang, Int. J. Mod. Phys. C , 817(1991).[13] D. C. Rapaport, The Art of Molecular Dynamics Simu-lation (Cambridge University Press, Cambridge, 1996).[14] C. S. Dias, N. A. M. Araujo, and M. M. Telo da Gama,J. Phys.: Cond. Matt. , 014001 (2018).[15] J. Stenhammar, D. Marenduzzo, R. J. Allen, and M. E.Cates, Soft Matt. , 1489 (2014).[16] M. R. Shaebani, A. Wysocki, R. G. Winkler, G. Gomp-per, and H. Rieger, Nat. Rev. Phys. , 181 (2020).[17] G. Volpe, S. Gigan, and G. Volpe, Am. J. Phys. , 659(2014).[18] C. Bechinger, R. Di Leonardo, H. L¨owen, C. Reichhardt,G. Volpe, and G. Volpe, Rev. Mod. Phys.88