Non-symmetric pair style - is it possible?

Hello everyone,
In the past years I have used lammps to perform some simulations on coarse grained model, so I’m not a complete noob - though I’m no expert.
I was wondering if there was a “simple” way to define non symmetric pair potentials. I have in mind something like this: we consider 2 atom types. The pair interaction is e.g. lj for 1-1 and 2-2 interactions. But is like lj with some parameters for 1-2 and different parameters, or 0 for 2-1 (e.g. the first particle is repelled or attracted by the second, but the second feels a different force - violating Newton’s principle). I understand that simply using pair_coeff in the form 2 1 does not work. Do you have any suggestions? Thank you for your time,

1 Like

I have a rather fundamental question: If you have a pair of two atoms and one is of type 1 and the other of type 2, how will you be able to tell whether it is a 1-2 or 2-1 pair?

Violating Newton’s third law will require writing a custom pair style in C++.
You will basically have to implement a LJ pair style that uses a full neighbor list instedd of the (default) half neighbor list. As you already noted, you will also have to be “creative” in how you enter potential parameters, you may have to read them from a file like it is done with manybody potentials.

I don’t even want to go into the specifics of what violating Newton’s third law means to the thermodynamics and statistical mechanics of the simulated system…

The most straightforward way to achieve this is to write a custom fix, which can request a neighbor list and iterate over it during post_force to calculate and add on any force you like, no matter how gloriously non-physical. I am deeply skeptical it is a good idea (how will you conserve momentum or energy or use a thermostat) but if I had to do it, that’s how I would.

Thank you very much for your feedback.
Indeed I was thinking of “hijacking” lammps to perform a simulation of a dynamical system with no real thermodynamics backing it up, à la “particle life”.
See e.g.:

I was hoping to be able to use lammps because it has a very convenient scripting language and is very fast - instead of having to write most code from scratch. The idea of using a homebrew version of lj potential is good, though being a lazybum I hoped I could get away with some dangerous/unlawful scripting hack :wink:
Thanks again,

I have found a full neighbor list version of pair style lj/cut called lj/cut/full, that I had written some time ago for testing purposes, lying around on my hard drive.

I am attaching it here. So to implement your model, you would just need to “massage” the code in the coeff() and init_one() functions to remove the symmetrization of the LJ parameters and disable the reordering in the Input::pair_coeff() function in input.cpp

or set one_coeff = 1; in the pair style constructor and write a reader for all parameters from a file in one go, like it is done for many-body potentials like tersoff or sw.

pair_lj_cut_full.h (1.3 KB)
pair_lj_cut_full.cpp (7.2 KB)


Thank you very much :slight_smile:

I believe this can be done in vanilla LAMMPS, albeit it is quite hacky. It requires using multiple fix nve, separately for each atom type, and pair_style lj/cut/coul/cut. You alternate between nve fixes so only one is active during a single step. Between steps, you change the charge of a chosen type of atoms. You also need fix viscous as this setup constantly creates kinetic energy, see the animation below:


The input script for the animation is below.

dimension       2
units           real
atom_style      charge

region          u block -20 20 -20 20 -0.5 0.5
create_box      2 u
mass            * 100

create_atoms    1 random 5    144 u 
create_atoms    2 random 100 14234 u 

pair_style      lj/cut/coul/cut 6.0
pair_coeff      * * 1 2

group           1 type 1
group           2 type 2

set             group 1 charge  2
set             group 2 charge -2

minimize        1e-10 1e-10 10000 10000

dump            dump all custom 10 dump.lammpstrj id element x y z
dump_modify     dump element Co Fe

thermo_style    custom step pe ke temp cpuremain
thermo          100

run             0
timestep        0.5

fix             viscous all viscous 10000.0

variable        loopa loop 20000

label           loopa
fix             nve1 1 nve
run             1 
unfix           nve1 
set             group 2 charge 2
fix             nve2 2 nve
run             1 
set             group 2 charge -2
unfix           nve2
next            loopa
jump            SELF loopa

Thanks, I’ll check this too