water molecules out of box boundary

LAMMPS LAMMPS Mailing List Mirror
1

Hi users,

I met a problem that the water molecules are not confined in my periodic box. At the beginning, I have about 9900 water molecules (SPCE water), after my simulation, ~ 1.5 ns, only about few hundred left. Other molecules are dispersed outside of the box. But the system never reported any error about losing atoms. Besides, I have 36 big particles in the box as well, but they are moving around but staying in the box. I am using a very small timestep 0.1 fs, and I did try smaller, the same result. Later, I also tried with a coarse-grain model water developed by Dr.Valeria Molinero, still same thing happens.

Is this problem possibly caused by my charged atoms? I saw a warning during the simulation said my system charge is not neutral.
This is my input script. Many thanks for your help.

----------------- Init Section -----------------

#include “system.in.init_m”
units real
atom_style hybrid full sphere
pair_style hybrid yukawa/colloid 2.0 2.5 lj/charmm/coul/long 9.0 10.0 10.0
bond_style hybrid harmonic fene
angle_style harmonic
kspace_style ewald 0.0001
pair_modify mix arithmetic
special_bonds fene

----------------- Atom Definition Section -----------------

read_data “systemw.data”

----------------- Settings Section -----------------

#include “system.in.settings_m”
bond_coeff 1 fene 1 2.5 1.5 1.0
bond_coeff 2 harmonic 200.0 1.0
angle_coeff 1 200.0 109.47
pair_coeff * * yukawa/colloid 100.0 2.3
pair_coeff 2 2 yukawa/colloid 100.0 2.3
pair_coeff 3 3 lj/charmm/coul/long 0.1553 3.166
pair_coeff 4 4 lj/charmm/coul/long 0.0 2.058

Load the atom coordinates:

include “systemw.in.coords”

----------------- Run Section -----------------

The lines above define the system you want to simulate.

What you do next is up to you.

Typically a user would minimize and equilibrate

the system using commands similar to the following:

---- examples ----

group particle molecule <> 1 36
group water type 3 4
fix 1 particle rigid molecule
communicate multi vel yes
fix fSHAKE water shake 0.0001 10 50 b 2 a 1

– minimize –

timestep 0.005
thermo_style custom step temp vol press etotal
thermo 100
dump min all custom 100 min.lammpstrj id mol type x y z ix iy iz
fix fxlan all langevin 298 298 100.0 48279
fix fxnve all nve
fix tmp all temp/rescale 100 280 330 1.0 1.0
run 10000

unfix fSHAKE

minimize 1.0e-4 1.0e-6 100 1000

(Note: Some fixes, for example “shake”, interfere with the minimize command,

You can use the “unfix” command to disable them before minimization.)

– declare time step for normal MD –

timestep 0.1

– run at constant pressure (Nose-Hoover)–

fix fxnpt all npt temp 300.0 300.0 100.0 iso 1.0 1.0 1000.0 drag 1.0

– ALTERNATELY, run at constant volume (Nose-Hoover) –

fix fxnvt all nvt temp 300.0 300.0 500.0 tchain 1

– ALTERNATELY, run at constant volume using Langevin dynamics. –

– (This is good for sparse CG polymers in implicit solvent.) –

#velocity all create 300 245634
#fix fxlan all langevin 300.0 300.0 100 48279
#fix fxnve all nve
#fix fxnvt all nvt temp 300 300 500
thermo_style custom step temp vol press etotal
thermo 5000
dump 1 all custom 5000 test.lammpstrj id mol type x y z ix iy iz
run 5000000
#run 5000000
#run 5000000

2

I assume you are viewing the dump file in VMD.
VMD does not "wrap" the coordinates of atoms or molecules by default.
If they look as though they are leaving the simulation box, this has
nothing to do with the way LAMMPS sees them.

To move all the water molecules into the same box (without cutting
them in half), try this:

Select the "Extensions"->"Tk Console" menu in VMD and enter:

pbc wrap -compound res -all
pbc box

If you need to, you can change the position of the box this way:

pbc wrap -compound res -all -shiftcenterrel {-0.3 -0.5 0}
pbc box -shiftcenterrel {-0.3 -0.5 0}

Cheers
Andrew

3

Try printing out the keyword “atoms” with your thermo-style.

If you’re saying that the atom count drops from 3*9900 to
some small number during the simulation and LAMMPS doesn’t
complain, then I don’t see how that is possible.

Steve

4

Additionally, go back and tune your script to make sure your simulations are meaningful.

Quick check: you are simultaneously using a langevin thermostat and fix/rescale on all the atoms of your system at two distinct temperatures., i.e.,

fix fxlan all langevin 298 298 100.0 48279
fix tmp all temp/rescale 100 280 330 1.0 1.0

This feels like complete nonsense unless you just forgot to comment one of the two in the input deck you attached.

Carlos

5

Hi Andrew,

I tried to “wrap” them in VMD and it works. Thanks so much.

Jing

6

Hi Carlos,

I didn’t actually run the fix/rescale in my simulation. But thank you for pointing it out anyway. My problem has solved. It turns out to be my VMD showing problem with the dump file.

Jing

7

Make a practice to always weed out all unused extras form the script, especially the wrong ones. If you ever have to come back to your input after six months to try reproducing or generating new data you’ll remember my words…

Carlos

8

I really appreciate your advice. I will follow :slight_smile: