How to keep a rigid body stationary

LAMMPS LAMMPS General Discussion
1

Hi LAMMPS-users,

There is a kaolinite slab in simulation box, which contains atoms Al Si O H. And the box is filled with water molecules. In simulation, I want to keep the heavy atoms Al Si O stationary, and the H atoms thermal motion is allowed. I have tried some method but all failed. For example, “fix 1 heavyatoms setforce 0 0 0, velocity heavyatoms set 0 0 0”, and “velocity heavyatoms zero linear rigid heavyatoms, velocity heavyatoms zero angular rigid heavyatoms”. So, how can I keep the Al Si O atoms in the slab stationary, and H atoms thermal motion?

Thanks

Best regards.

2

To give meaningful advice, we have to see your input file.

3

OK, I am a new user, the files can not be uploaded. The content of input is copied as below. Thanks for you advice.

variable temp equal 298.15
units real
dimension 3
boundary p p p
atom_style full
pair_style lj/cut/coul/long 10.0 10.0
pair_modify mix arithmetic tail no
bond_style harmonic
angle_style harmonic
kspace_style pppm 1.0e-4
kspace_modify order 7 mesh 0 0 0

read_data kao_test.data
molecule h2omol H2O.txt #toff 5 boff 5 aoff 12
variable xlo equal -0.803258291
variable xhi equal 79.196741709
variable lx equal $(v_xhi-v_xlo)
variable ylo equal 0.198101712
variable yhi equal 80.198101712
variable ly equal $(v_yhi-v_ylo)
variable zlo equal -79.170501514
variable zhi equal 0.829498486
variable lz equal $(v_zhi-v_zlo)

region kaoplate block 10 69 10 69 -43 -34 side in units box
region simbox block ${xlo} ${xhi} ${ylo} ${yhi} ${zlo} ${zhi} side in units box
region kaoout block 10 69 10 69 -43 -34 side out units box
region water intersect 2 simbox kaoout units box

variable kaov equal $((69-10)(69-10)(-34+42))
variable boxv equal $(v_lxv_lyv_lz)
variable waterv equal $(v_boxv-v_kaov)
variable watern equal $(v_waterv1E-246.02E23/18)
variable waterN equal ceil(v_watern)

create_atoms 0 random ${waterN} 666666 water mol h2omol 666666 units box

group kao type 1 2 3 4 5
group heavyatoms type 1 2 3 4
group water type 6 7
group hwater type 5 6 7
group h type 5

neigh_modify delay 0 every 1 check yes exclude group heavyatoms heavyatoms

minimize 1e-6 1e-8 10000 50000

fix 1 heavyatoms rigid single
fix 2 heavyatoms setforce 0 0 0
velocity water create ${temp} 666666 dist gaussian loop geom units box
velocity h create ${temp} 666666 dist gaussian loop geom units box
velocity heavyatoms set 0 0 0 rigid 1 units box
#nve/limit
fix 0 hwater nve/limit 0.5
fix 0temp hwater langevin ${temp} ${temp} 100 666666

thermo_style custom step spcpu atoms temp pe ke etotal
thermo 1000

dump dump1 all custom 1000 test1.dump id type element x y z
dump_modify dump1 sort id element Al Si O O H O H
timestep 1
run 50000
write_data nve0.data
unfix 0
unfix 0temp
undump dump1

variable deltat equal 1
timestep ${deltat}

compute 3temp hwater temp
compute_modify 3temp dynamic/dof yes
fix 3 hwater nvt temp ${temp} ${temp} $(v_deltat*100)
fix_modify 3 temp 3temp

thermo_style custom step spcpu atoms temp pe ke etotal
thermo_modify lost ignore flush yes temp 3temp
thermo 1000

variable Nsteps equal 100000

dump dump2 all custom $(v_Nsteps/100) test2.dump id type element x y z
dump_modify dump2 sort id element Al Si O O H O H

run ${Nsteps}
write_data test.data

4

Please see the forum guidelines for suggestions about how to correctly quote input files.

That post will also explain that it is crucial to report which LAMMPS version you are using and what platform you are running on.

Since you define fix setforce after the minimize command, the minimization will move all atoms.
You have excluded the interaction between the “heavyatoms” group atoms, but the forces due to all other atoms still apply. Please also note that

You would likely want to use the “overlap” keyword here to avoid having close contacts.

This is superfluous and possibly counterproductive. You don’t need to use a rigid fix to immobilize atoms. Just just don’t include them in any time integration and they won’t move; simple as that.

These two commands would not be needed when dropping fix rigid.

Is your water model rigid (TIP3P, or SPC/E, or TIP4P, or similar)? Then there seems to be a “fix shake” command missing, otherwise, your timestep is too large for the movement of the bonds involving hydrogen atoms. Please note that with recent LAMMPS versions it is also possible to use fix shake during minimization (it will use an approximation for SHAKE with strong harmonic bonds).

I don’t see you adding or removing atoms. Why this command?

Do you want atoms to be removed from the simulation? If yes, that does not make much sense since your have periodic boundaries in all directions. Lost atoms would be a sign of some serious problems with the simulation and thus they should not be ignored.

Since you are simulating a slab geometry, you should consider using boundary p p f and kspace_modify slab 3.0 to avoid a) water molecules approaching the slab from the “wrong” side and b) decouple the long-range Coulomb interactions across the periodic boundary.