1Writing new styles 2------------------ 3 4The :doc:`Modify` section of the manual gives an overview of how LAMMPS can 5be extended by writing new classes that derive from existing 6parent classes in LAMMPS. Here, some specific coding 7details are provided for writing code for LAMMPS. 8 9Writing a new fix style 10^^^^^^^^^^^^^^^^^^^^^^^ 11 12Writing fixes is a flexible way of extending LAMMPS. Users can 13implement many things using fixes: 14 15- changing particles attributes (positions, velocities, forces, etc.). Examples: FixNVE, FixFreeze. 16- reading/writing data. Example: FixRestart. 17- adding or modifying properties due to geometry. Example: FixWall. 18- interacting with other subsystems or external code: Examples: FixTTM, FixExternal, FixLATTE 19- saving information for analysis or future use (previous positions, 20 for instance). Examples: Fix AveTime, FixStoreState. 21 22 23All fixes are derived from the Fix base class and must have a 24constructor with the signature: ``FixPrintVel(class LAMMPS *, int, char **)``. 25 26Every fix must be registered in LAMMPS by writing the following lines 27of code in the header before include guards: 28 29.. code-block:: c 30 31 #ifdef FIX_CLASS 32 FixStyle(print/vel,FixPrintVel) 33 #else 34 /* the definition of the FixPrintVel class comes here */ 35 ... 36 #endif 37 38Where ``print/vel`` is the style name of your fix in the input script and 39``FixPrintVel`` is the name of the class. The header file would be called 40``fix_print_vel.h`` and the implementation file ``fix_print_vel.cpp``. 41These conventions allow LAMMPS to automatically integrate it into the 42executable when compiling and associate your new fix class with the designated 43keyword when it parses the input script. 44 45Let's write a simple fix which will print the average velocity at the end 46of each timestep. First of all, implement a constructor: 47 48.. code-block:: C++ 49 50 FixPrintVel::FixPrintVel(LAMMPS *lmp, int narg, char **arg) 51 : Fix(lmp, narg, arg) 52 { 53 if (narg < 4) 54 error->all(FLERR,"Illegal fix print/vel command"); 55 56 nevery = force->inumeric(FLERR,arg[3]); 57 if (nevery <= 0) 58 error->all(FLERR,"Illegal fix print/vel command"); 59 } 60 61In the constructor you should parse your fix arguments which are 62specified in the script. All fixes have pretty much the same syntax: 63``fix <fix-ID> <fix group> <fix name> <fix arguments ...>``. The 64first 3 parameters are parsed by Fix base class constructor, while 65``<fix arguments>`` should be parsed by you. In our case, we need to 66specify how often we want to print an average velocity. For instance, 67once in 50 timesteps: ``fix 1 print/vel 50``. There is a special variable 68in the Fix class called ``nevery`` which specifies how often the method 69``end_of_step()`` is called. Thus all we need to do is just set it up. 70 71The next method we need to implement is ``setmask()``: 72 73.. code-block:: C++ 74 75 int FixPrintVel::setmask() 76 { 77 int mask = 0; 78 mask |= FixConst::END_OF_STEP; 79 return mask; 80 } 81 82Here the user specifies which methods of your fix should be called 83during execution. The constant ``END_OF_STEP`` corresponds to the 84``end_of_step()`` method. The most important available methods that 85are called during a timestep and the order in which they are called 86are shown in the previous section. 87 88.. code-block:: C++ 89 90 void FixPrintVel::end_of_step() 91 { 92 // for add3, scale3 93 using namespace MathExtra; 94 95 double** v = atom->v; 96 int nlocal = atom->nlocal; 97 double localAvgVel[4]; // 4th element for particles count 98 memset(localAvgVel, 0, 4 * sizeof(double)); 99 for (int particleInd = 0; particleInd < nlocal; ++particleInd) { 100 add3(localAvgVel, v[particleInd], localAvgVel); 101 } 102 localAvgVel[3] = nlocal; 103 double globalAvgVel[4]; 104 memset(globalAvgVel, 0, 4 * sizeof(double)); 105 MPI_Allreduce(localAvgVel, globalAvgVel, 4, MPI_DOUBLE, MPI_SUM, world); 106 scale3(1.0 / globalAvgVel[3], globalAvgVel); 107 if ((comm->me == 0) && screen) { 108 fmt::print(screen,"{}, {}, {}\n", 109 globalAvgVel[0], globalAvgVel[1], globalAvgVel[2]); 110 } 111 } 112 113In the code above, we use MathExtra routines defined in 114``math_extra.h``. There are bunch of math functions to work with 115arrays of doubles as with math vectors. It is also important to note 116that LAMMPS code should always assume to be run in parallel and that 117atom data is thus distributed across the MPI ranks. Thus you can 118only process data from local atoms directly and need to use MPI library 119calls to combine or exchange data. For serial execution, LAMMPS 120comes bundled with the MPI STUBS library that contains the MPI library 121function calls in dummy versions that only work for a single MPI rank. 122 123In this code we use an instance of Atom class. This object is stored 124in the Pointers class (see ``pointers.h``) which is the base class of 125the Fix base class. This object contains references to various class 126instances (the original instances are created and held by the LAMMPS 127class) with all global information about the simulation system. 128Data from the Pointers class is available to all classes inherited from 129it using protected inheritance. Hence when you write you own class, 130which is going to use LAMMPS data, don't forget to inherit from Pointers 131or pass an Pointer to it to all functions that need access. When writing 132fixes we inherit from class Fix which is inherited from Pointers so 133there is no need to inherit from it directly. 134 135The code above computes average velocity for all particles in the 136simulation. Yet you have one unused parameter in fix call from the 137script: ``group_name``. This parameter specifies the group of atoms 138used in the fix. So we should compute average for all particles in the 139simulation only if ``group_name == "all"``, but it can be any group. 140The group membership information of an atom is contained in the *mask* 141property of and atom and the bit corresponding to a given group is 142stored in the groupbit variable which is defined in Fix base class: 143 144.. code-block:: C++ 145 146 for (int i = 0; i < nlocal; ++i) { 147 if (atom->mask[i] & groupbit) { 148 // Do all job here 149 } 150 } 151 152Class Atom encapsulates atoms positions, velocities, forces, etc. User 153can access them using particle index. Note, that particle indexes are 154usually changed every few timesteps because of neighbor list rebuilds 155and spatial sorting (to improve cache efficiency). 156 157Let us consider another Fix example: We want to have a fix which stores 158atoms position from previous time step in your fix. The local atoms 159indexes may not be valid on the next iteration. In order to handle 160this situation there are several methods which should be implemented: 161 162- ``double memory_usage()``: return how much memory the fix uses (optional) 163- ``void grow_arrays(int)``: do reallocation of the per particle arrays in your fix 164- ``void copy_arrays(int i, int j, int delflag)``: copy i-th per-particle 165 information to j-th. Used when atom sorting is performed. if delflag is set 166 and atom j owns a body, move the body information to atom i. 167- ``void set_arrays(int i)``: sets i-th particle related information to zero 168 169Note, that if your class implements these methods, it must call add calls of 170add_callback and delete_callback to constructor and destructor. Since we want 171to store positions of atoms from previous timestep, we need to add 172``double** xold`` to the header file. Than add allocation code 173to the constructor: 174 175.. code-block:: C++ 176 177 FixSavePos::FixSavePos(LAMMPS *lmp, int narg, char **arg), xold(nullptr) 178 { 179 //... 180 memory->create(xold, atom->nmax, 3, "FixSavePos:x"); 181 atom->add_callback(0); 182 } 183 184 FixSavePos::~FixSavePos() { 185 atom->delete_callback(id, 0); 186 memory->destroy(xold); 187 } 188 189Implement the aforementioned methods: 190 191.. code-block:: C++ 192 193 double FixSavePos::memory_usage() 194 { 195 int nmax = atom->nmax; 196 double bytes = 0.0; 197 bytes += nmax * 3 * sizeof(double); 198 return bytes; 199 } 200 201 void FixSavePos::grow_arrays(int nmax) 202 { 203 memory->grow(xold, nmax, 3, "FixSavePos:xold"); 204 } 205 206 void FixSavePos::copy_arrays(int i, int j, int delflag) 207 { 208 memcpy(xold[j], xold[i], sizeof(double) * 3); 209 } 210 211 void FixSavePos::set_arrays(int i) 212 { 213 memset(xold[i], 0, sizeof(double) * 3); 214 } 215 216 int FixSavePos::pack_exchange(int i, double *buf) 217 { 218 int m = 0; 219 buf[m++] = xold[i][0]; 220 buf[m++] = xold[i][1]; 221 buf[m++] = xold[i][2]; 222 223 return m; 224 } 225 226 int FixSavePos::unpack_exchange(int nlocal, double *buf) 227 { 228 int m = 0; 229 xold[nlocal][0] = buf[m++]; 230 xold[nlocal][1] = buf[m++]; 231 xold[nlocal][2] = buf[m++]; 232 233 return m; 234 } 235 236Now, a little bit about memory allocation. We use the Memory class which 237is just a bunch of template functions for allocating 1D and 2D 238arrays. So you need to add include ``memory.h`` to have access to them. 239 240Finally, if you need to write/read some global information used in 241your fix to the restart file, you might do it by setting flag 242``restart_global = 1`` in the constructor and implementing methods void 243``write_restart(FILE *fp)`` and ``void restart(char *buf)``. 244If, in addition, you want to write the per-atom property to restart 245files additional settings and functions are needed: 246 247- a fix flag indicating this needs to be set ``restart_peratom = 1;`` 248- ``atom->add_callback()`` and ``atom->delete_callback()`` must be called 249 a second time with the final argument set to 1 instead of 0 (indicating 250 restart processing instead of per-atom data memory management). 251- the functions ``void pack_restart(int i, double *buf)`` and 252 ``void unpack_restart(int nlocal, int nth)`` need to be implemented 253 254