Troubles with coulomb potential

LAMMPS LAMMPS Beginners
1

Hello everybody. It’s my first week using LAMMPS, so I am still wrapping my mind on it. I am writing a pretty simple script to simulate a bunch of particles interacting through Coulomb potential, but it doesn’t matter which particular combination of coefficient I use, I am unable to get a useful result. We are talking about a very linear setting: a box, few particles together, pair_style could/cut, and let them interact! But nope, what happens here is that after a few steps (about 200), 99% of the particles freeze and a few of them run wildly through the space. I am posting here the script, hoping that someone would be kind to knock me over it pointing at the very dumb mistake I am doing.

3D Coulomb melt

units lj
atom_style charge
dimension 3

boundary p p p
lattice fcc 0.2

region box block 0 50 0 50 0 50
create_box 2 box
region droplet sphere 25 25 25 10

create_atoms 1 random 20 12345 droplet # anion
create_atoms 2 random 20 54321 droplet # cation

Mass to atoms

mass 1 1.0
mass 2 1.05

Charges to atoms

group g1 type 1
group g2 type 2

set group g1 charge -0.8
set group g2 charge 0.8

velocity all create 0.5 87287

pair_style lj/cut/coul/cut 5.0 8.0
pair_coeff * * 5.0 3.0

neighbor 0.2 bin
neigh_modify every 20 delay 0 check yes
fix 1 all nvt temp 0.05 0.05 1.0

dump id all atom 50 dump.lammpstrj
thermo 100
run 100000
dump m1 all movie 1000 movie.avi type type size 640 480

2

The effect you are seeing is due to using a thermostat. That will try to curb the total kinetic energy, so if some particles accelerate a lot all the other particles will be slowed down.

Consequently, the next question would be: why are these atoms so fast? That requires some speculation, but the most likely cause is that either your Lennard-Jones interaction or your Coulomb interaction is out of whack or your time step is too large or a combination of those. You can either play with the parameters and see what kind of epsilon you have to use to give meaningful results in combination with your Coulomb, or you could just switch to real (or metal) units (i.e. non-reduced units) and use some parameters from a real simulation. You can find some reasonable parameter sets for “real” units in the “rhodo” benchmark input or the “peptide” example.

3

Thank you so much, I will follow your advice! I was thinking about the energy, yes, given the strange behaviour, but I didn’t really know how to curb the effect. Now I have a starting point for the investigation. Thank you again!