[lammps-users] Lennard-Jones + Coulomb

Hi everybody,

As a newcomer to LAMMPS, I definitely need your help.
I would like to run a simulation with only two types of pair potentials (no bond, angle,dihedral,etc): Lennard-Jones and Coulomb. As a first approach, I chose the lj/cut/coul/cut pair_style.
My input script contains the following lines.
Is it correct to use the line "atom_style charge" combined with the command "pair_style lj/cut/coul/cut 8.75 20" ?
(This script also calls the file system.dat that defines the masse of the atoms and the atoms via the following line:
1 1 -1 -331.074383386 -331.074383386 -1.32953354535 (Atom-ID, Atom-type, Charge, x, y, z))

units real
atom_style charge

pair_style lj/cut/coul/cut 8.75 20

read_data system.dat

pair_coeff * * 100 3.5 8.75 20

# Define groups
group water 1
group solid_neutral 2
group solid_charged 3
group cations 4
group anions 5

# temp controllers
temperature water_temp water full
temperature cations_temp cations full
temperature anions_temp anions full
temperature solid_neutral_temp solid_neutral full
temperature solid_charged_temp solid_charged full

velocity water create 298.0 887723 temp water_temp
velocity cations create 298.0 887723 temp cations_temp
velocity anions create 298.0 887723 temp anions_temp
velocity solid_neutral create 298.0 887723 temp solid_neutral_temp
velocity solid_charged create 298.0 887723 temp solid_charged_temp

# Fix solid atom around initial position
fix tether solid_neutral spring/self 100
fix tether solid_charged spring/self 100

# Keep solid temperature constant
fix 1 solid_neutral nvt 1.00 1.00 10
fix 2 solid_charged nvt 1.00 1.00 10
fix 3 water nve
fix 4 cations nve
fix 5 anions nve

# Enforce 2d calculation
fix 1 all enforce2d

# Save results every 100 frames (positions and thermodynamics data)
dump 1 all xyz 100 electrokinetics.out
thermo 100
run 100000

So far it crashes immmdiately. There must be something wrong with the initial position of the atoms or an inproper charge or energy (unit problem). I am looking for the problem but I just would like to know if I am going in the right direction.

I have another question. Most of the atoms are not charged, so the LJ potential is enough. Some atoms are charged, therefore I need to had the Coulombic interaction. But so far all the atoms undergo both potentials. I just put either a zero or non-zero charge whether the atom is interacting with LJ or LJ+Coulomb.
Is it possible to run the LJ potential only on some types of atoms and run both LJ and Coulombic on other types of atoms ?

I thank you ahead for your help and possible answers.

Have a good day.

Just use atom_style charge and the lj/cut/coul/cut potential
even if most atoms are not charged. If your initial atom coords
or charges produce large forces, then the system can crash.


Thank for your answer Steve,

The simulation is now, at least, starting. There were some mistakes in the input script like "group water 1" instead of "group water type 1" which lead to "Segmentation fault". The initial atom positions seems to be actually correct.


Steve Plimpton wrote:

The command "group water 1" should have thrown an error
and not crashed. It is due to a bug in how the command is
parsed. I'll post a patch this AM.


dear list

i am trying to simulate the granular flow in a vertical channel
i modified fix_gran_diag.cpp for my problem
i made x-flowing direction , gravity acts in x-direction only
y-direction will be width and walls will be at y=0.0 and y=10.0
z-spanwise direction
i have total 4069 particles and all are active, i am not having ant frozen base as i kept flat frictional walls at both ends in y-direction

my problem is that system is not reaching steady state

velocity in x-direction kept on increasing even after 2 crores of time steps

please let me know any suggestions if you have regarding this problem

my input script :

# LAMMPS benchmark of granular flow
# vertical chute flow of 4069 atoms
# x-flowing direction-gravity,y-width,z-spanwise

units lj
atom_style granular
boundary p f p
newton off

read_data chute.dat

pair_style gran/history 200000.0 200.0 0.5 1

neighbor 0.1 bin
neigh_modify every 1 delay 0

timestep 0.0001

#group bottom type 2
#group active subtract all bottom
#neigh_modify exclude group bottom bottom

fix 1 all gravity chute 90.0
#fix 2 bottom freeze
fix 3 all wall/gran yplane 0.0 10.0 200.0 0.5
fix 4 all nve/gran
fix 5 all gran/diag 1000 tmp 0.5

thermo_style granular
thermo 1000
thermo_modify lost ignore
dump 1 all atom 1000 dump.chute
dump_modify 1 scale no
restart 50000 tmp.restart

run 20000000

thanks and regards

Your gravity appears to be at 90 degrees or in the x direction. This
means the particles are free-falling towards x. So they will
accelerate forever, hardly touching the walls?

More typical is a gravity direction of 20-30 degrees, tilted towards
x, as in examples/in.pour. Too small a value and the particles don't
flow due to friction w/ the walls. Too large, and they will become