Modeling the movement of salt ions in aqueous solution under the influence of electric fields

Brad_Roan · September 23, 2023, 5:46pm

Hello! I’m trying to model the behavior of aqueous ions in an ion trap, and the results that I’m getting seem somewhat erratic. I can’t seem, to arrive at a correct value for for the damp factor in the fix ID group-ID langevin Tstart Tstop damp. In a field of 100 volt / meter, a chloride ion of 5.88e-26 kg, charge of 1.602e-19 C, stokes radius of 1.24e-10 m, in water with a dynamic viscosity of 8.9e-4 kg/(ms), should have a ‘damp’ factor of 2.83e-14 as I understand it. Which I arrive at via [5.88e-26 kg]/(6PI*[8.9e-4 kg/(ms)][1.24e-10 m]). Given the drift velocity of sodium in water of 7.7e-8 (m^2)/(V*S), the velocity of a chloride ion in the above mentioned field would be a paltry 7.7e-6 m/s. However, when subjected to these conditions, the particle rapidly accelerates out of the bounds of the simulation, despite my 1e-6s time step.

Could someone please help me understand the error in my code?
Am I using the ‘damp’ factor properly, and messing up somewhere else?

I admit, I’m not using LAMMPS directly, but a this python wrapper but I am digging into the LAMMPS documentation to try to understand the LAMMPS file that it is outputting:

units si
atom_style charge

# Creating simulation box...
boundary m m m
region simulationDomain block -0.005  0.005 -0.005 0.005 -0.005 0.005 units box
create_box 1 simulationDomain

# Configure neighbour list...
neighbor 1 nsq
neigh_modify once yes

# Configure pairwise interactions for long range Coulombics only...
pair_style coul/cut 10
pair_coeff * *

# Placing individual ions...
create_atoms 1 single 0.0045 0.0045 0 units box

# Species...
mass 1 5.894913449999999e-26
set type 1 charge 1.6021764599999998e-19
group 1 type 1

timestep 1e-06

# Configuring additional output to flush buffer during simulation...
thermo 10000
thermo_style custom step cpu
thermo_modify flush yes

# Time integration...
group nonRigidBody union all
fix timeIntegrator nonRigidBody nve


# Creating a Linear Paul Trap... (fixID=3903137266868)
variable endCapV3903137266868		equal 1.000000e+00
variable radius3903137266868		equal 5.000000e-03
variable zLength3903137266868		equal 5.000000e-03
variable geomC3903137266868		equal 1.000000e+00

# Define frequency components.
variable oscVx39031372668680		equal 1.000000e+00
variable oscVy39031372668680		equal 1.000000e+00
variable phase39031372668680		equal "0.000000e+00 * step*dt"
variable oscConstx39031372668680		equal "v_oscVx39031372668680/(v_radius3903137266868*v_radius3903137266868)"
variable oscConsty39031372668680		equal "v_oscVy39031372668680/(v_radius3903137266868*v_radius3903137266868)"
variable statConst3903137266868		equal "v_geomC3903137266868 * v_endCapV3903137266868 / (v_zLength3903137266868 * v_zLength3903137266868)"

variable oscEX3903137266868 atom "v_oscConstx39031372668680 * cos(v_phase39031372668680) * x + v_statConst3903137266868 * x"
variable oscEY3903137266868 atom "v_oscConsty39031372668680 * cos(v_phase39031372668680) * -y + v_statConst3903137266868 * y"
variable statEZ3903137266868 atom "v_statConst3903137266868 * 2 * -z"
fix 3903137266868 all efield v_oscEX3903137266868 v_oscEY3903137266868 v_statEZ3903137266868


# Adding a langevin bath...
fix 5853822124365 all langevin 0.000000e+00 0.000000e+00 2.830000e-14 1337

dump 3902694972670 all custom 10 positions.txt id x y z

fix 3902689020878 all ave/atom 1 20 20 vx vy vz 

dump 5854270558203 all custom 200 secv.txt id f_3902689020878[1] f_3902689020878[2] f_3902689020878[3]


# Run simulation
run 10000

Thank you for any help that you can provide!

akohlmey · September 23, 2023, 11:04pm

Two comments:

your post is very difficult to read since you didn’t quote your input correctly. see Please Read This First: Guidelines and Suggestions for posting LAMMPS questions
since you are using a “meta” package to create your LAMMPS simulation input, your best changes for on-topic advice is probably to contact the author of that package or post in a suitable forum for that software.

srtee · September 25, 2023, 6:39am

One microsecond is a very implausible integration timestep for simulating the translational dynamics of atoms.

In your situation I would (1) calculate the velocity autocorrelation time for the Langevin dynamics involved, (2) choose an initial timestep of 1/10th to 1/20th of that time, and (3) work my way both downwards and upwards to confirm no major disasters. I think that should land you in the picosecond range; usually MD is done with femtosecond timesteps.

This is based purely on the heuristic that molecular modellers can frequently test the correctness of their Langevin dynamics by calculating the velocity autocorrelation, which is certainly not possible if their timestep far exceeds the velocity autocorrelation time.