A question about variable loop

Hi everyone.
I build a silica-water system to study the hydrolysis reaction with REAXFF in a constant strain rate load. I use “fix setforce” and “velocity set” command to freeze two sides of silica bulk, then I give the left side a constant velocity and zero velocity alternately to load a quasi-static tension.
But it doesn’t work as I said. I wonder if my command is wrong or something else.
THANK YOU advanced.

Here is my code:
units real
dimension 3
boundary f p f
atom_style full

package gpu 0 device_type nvidiagpu
suffix gpu

pair_style reaxff NULL safezone 3.0 mincap 10000
read_data …/SiO2_water_surface/SiO2_water_surface.data
pair_coeff * * …/ffield_CaSiAlO_Pitman Si O O O H

neighbor 2.0 bin
neigh_modify delay 10
change_box all x delta -8 0
group water type 4 5
fix 1 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff maxiter 10000
fix 2 water wall/reflect zhi EDGE xlo 0 xhi 21
reset_atoms id
thermo_modify lost ignore

group silica_atoms type 1 2 3
variable xmax equal bound(silica_atoms,xmax)+0.5
variable xmin equal bound(silica_atoms,xmin)-0.5
variable ymax equal bound(silica_atoms,ymax)+0.5
variable ymin equal bound(silica_atoms,ymin)-0.5
variable zmax equal bound(silica_atoms,zmax)+0.5
variable zmin equal bound(silica_atoms,zmin)-0.5
variable x_left equal {xmin}+3.5 variable x_right equal {xmax}-3.5
region left block {xmin} {x_left} {ymin} {ymax} {zmin} {zmax}
region right block {x_right} {xmax} {ymin} {ymax} {zmin} {zmax}

group silica_left region left
group silica_right region right
group rest subtract all silica_left silica_right

fix left_force silica_left setforce 0 0 0
fix right_force silica_right setforce 0 0 0
velocity silica_left set 0 0 0 units box
velocity silica_right set 0 0 0 units box

dump 1 all atom 1000 dump.lammpstrj
fix 4 rest reaxff/bonds 2000 bonds.reaxff
thermo 1000
timestep 0.2

variable n loop 100
label loop
velocity silica_left set -0.00002 0 0 units box
fix 5 rest nvt temp 300.0 300.0 100.0
run 10000
velocity silica_left set 0 0 0 units box
fix 7 rest nvt temp 300.0 300.0 100.0
run 5000
next n
jump input.lammps loop


It is better to be more specific than that when you want help on a forum. Not working can mean a lot of things.

What equation of motion integrator is applied to your silica_left group? I don’t see any. Also you have a loop with 2 fixes nvt with different name applying on the same group, which will lead to the system explode.


Hi,Simon. I’m sorry that I wasn’t clear. Accurately the bulk which I give a constant velocity doesn’t move entirely. I think the command “fix setforce 0 0 0” and “velocity set 0 0 0” work because it’s freezed in ovito. But command “velocity set -0.00002 0 0” doesn’t move the bulk.

Also you have a loop with 2 fixes nvt with different name applying on the same group, which will lead to the system explode.

And BTW, I observed an incorrect dispersion of atoms on the second fix command, is there a relationship?

And there is no equation of motion integrator in silica_left group, because this group is used to load the strain.

I suppose yes

then it canno’t move

I suppose yes

So I should use the same name in "fix _ rest nvt"command?

then it canno’t move

Sorry I’m a new lammps user, I still don’t know why it cannot move even I assign a velocity to this group.

Then you should try simpler simulations first, and probably follow tutorials too. What you are trying to achieve is a very advanced molecular dynamics simulation, and you absolutely need to understand the basics of LAMMPS to make it work. Otherwise, you will keep facing new issues.

Thank you Simon, I already solve the problem of incorrect dispersion of atoms. And I will learn more about LAMMPS.

Thank you for your reply again!

My suggestion would be to avoid all possible uses of next loop or jump in LAMMPS, only demons know how to use those commands