Coupling lammps through external force

Post date: May 20, 2016 6:31:29 PM

Lammps can be used as a library to an external code as well. Typically during the coupling process, the information between lammps and external code is exchanged. Particle position x, velocity v and force f are usually exchanged. Another nice feature is that those variables are public members, thus can be accessed anywhere. Thus, it is quite simple to change particle information x and v but NOT force f, which will be discussed below.

Let us look at the simple.cpp provided in COUPLE/simple in lammps folder. It provides two functions to gather and scatter atom information. We can easily change x outside of normal lammps run loop and send it back, as demonstrated in the source code.

lammps_gather_atoms(lmp,"x",1,3,x) x[0] = epsilon; lammps_scatter_atoms(lmp,"x",1,3,x)

The same is true for velocity v. Actually the gather/scatter function is not needed. You can directly change it like this:

double **x = lmp->atom->x; x[1][2] = epsilon;

However, if you want to change f, both approaches fail.

The problem is that during each time step, the force vector will be cleared once. e.g., through force_clear() in verlet.cpp. Anything we changed for f outside lammps between time steps will be erased! However, this is not the case for position and velocity, as they are preserved and continuous between time steps.

Lammps provide a "fix external" to solve this problem. Basically it saves the external force to a public member **fexternal, and then add it the **f during post_force(). You have to set a function pointer in the fix. I find it difficult to use. Then I followed "fix external" source code and created a new fix style. Instead of setting function pointer to callback, I directly change the public external force vector **fexternal in my external driver code through fix pointer static_cast conversion, as I need to access the public member of a child class.

FixMyStyle *pf; pf = static_cast(lmp->modify->fix[ifix]); double **fe = pf->fexternal; fe[i][0]=epsilon;

Then the fix will add **fexternal to **f through post_force(). Now the external force has been successfully added.