N-body simulation of ions in a harmonic trap with ion-ion Coulombic interactions


I am new to LAMMPS. I wish to simulate N=10000 singly charged ions (same mass, m) in a 3D harmonic trap with frequencies fx ≠ fy ≠ fz for about 100 trap oscillation cycles. These particles experience Coulomb repulsion from each other. I am testing this out on a test case of N=2, m=199amu, fx=fy=fz=1.6kHz with the following code, but it seems like the system is not evolving in time at all from what I see in the dump.

This seems like an embarrassingly simple problem, but I have been stuck for quite a while. I believe I am missing something critical in the code. Please enlighten, thank you!

Point particles in 3D

variable m equal 199*1.66e-27

units si # SI units
dimension 3
newton on # N3L to save computation time at the cost of communication
boundary f f f
atom_style charge

region ion_cloud_region sphere 0 0 0 0.01 # center of x y z and radius
region sim_region sphere 0 0 0 0.02
create_box 1 sim_region # types of atoms, region
create_atoms 1 random 2 1903 ion_cloud_region # types of atoms, region

mass 1 $m
set atom * charge 1.6e-19

pair_style coul/cut 0.1
pair_coeff * *
pair_coeff 1 1 0.02

#neighbor 0.3 bin # extra distance after force cutoff for pairwise list
velocity all create 50 1903 dist gaussian

variable pi equal 3.14159
variable wx equal 2*{pi}*1.6e3 variable wy equal 2*{pi}1.6e3
variable wz equal 2
variable trapXs atom -m*{wx}^2
variable trapYs atom -m*{wy}^2y
variable trapZs atom -m*{wz}^2
variable trapEnergy atom 0.5*{m}*(({wx}*v_trapXs)^2+({wy}*v_trapYs)^2+({wz}*v_trapZs)^2)

Harmonic trap

fix harmonictrap all addforce v_trapXs v_trapYs v_trapZs energy v_trapEnergy

timestep 1e-5

#dump positions all atom 1 position.bp
dump positions all custom 1 position.bp id x y z vx vy vz fx fy fz q mass

run 10000


Please check the output of the simulation carefully. There does not seem to be any time integration and LAMMPS should warn about it.