simulating water using fix shake

I did not look at your files but one question immediately posed up: How did you cut your box/water molecules and during the cutting process did you cut though any bonds with water molecules?

Ray

I cut it by vmd.

First, I read my output file .dcd in vmd. I use “set sel [atomselect top x<0 ]” to choose the water molecules that I do not want. Then I use “$sel get resid” to get their resid, and type "sel set [atomselect top “not resid #”] to choose the molecules I want. This is because each molecule has the same resid, in the last command I can avoid cutting through bonds. Then I output the file into data file for LAMMPS to read. However, this is what I figure out by myself and I am not sure is this method totally correct. What I am sure is that I will not cut through bonds.
Thank you.

Miao-Chun

I cut it by vmd.

First, I read my output file .dcd in vmd. I use "set sel [atomselect top
x<0 ]" to choose the water molecules that I do not want.

​to extract the atoms that you want directly and to make sure to have whole
molecules, you can ​try a selection like the following:

set sel [atomselect top "not (same fragment as (x < 0))"]

​this one may also work​, if molecules are properly labeled with molecule
ids, but fragment should be safer.

​set sel [atomselect top "not (same residue as (x < 0))"]​

using "resid" instead is not a good idea, since "resid" is not required to
be unique and this can lead to unwanted inclusion of stray atoms that just
by chance have the same residue id, without being in the same molecule.

​axel.​

Dear lammps users,

Thank you Axel for correcting the cutting water method.
I recut my water, but step 0 still happens nan in pressure and the simulation crush. Does anyone has any idea? I attach my simulation below, this simulation is after cutting the water into desire shape, I merged the water and pyrophyllite together. I make two pyrophyllite clay particles which sandwich the water.
Thank you.

Miao-Chun

units real
dimension 3
boundary p p p
atom_style full

read_data system.data
neighbor 2.0 bin
neigh_modify delay 4 every 2 check yes

set potential function and parameters

pair_style lj/cut/coul/long 13.0

#pair_coeff i j \epslon \sigma
#water
pair_coeff 1 1 0.152 3.53 #O-O
pair_coeff 2 2 0.046 0.44 #H-H

#pyrophyllite
pair_coeff 3 3 0.001 7.4 # Si-Si
pair_coeff 4 4 0.150 6.3 # Al-Al
pair_coeff 5 5 0.0001 2.4 # H1-H1
pair_coeff 6 6 1.0 3.0 # O1-O1
pair_coeff 7 7 6.0 2.8 # O2-O2
pair_coeff 8 8 6.0 2.8 # O3-O3

pair_modify mix arithmetic tail no

bond_style harmonic
bond_coeff 1 450.0 0.9572 #water

angle_style harmonic
angle_coeff 1 34.933 104.52 #water

kspace_style pppm 1.0e-5
#kspace_modify slab 3.0e-5

group slab type 3 4 5 6 7 8
group water type 1 2

compute Tliq water temp

thermo 10
thermo_style custom step pe etotal press temp c_Tliq
thermo_modify norm no lost ignore lost/bond ignore

timestep 0.25

velocity water create 298.0 482748 dist gaussian
fix 1 water nvt temp 298.0 298.0 100
fix shake water shake 1.0e-4 100 0 b 1 a 1
#fix 2 water temp/rescale 10 298.0 298.0 0.02 1.0
dump dcd_1 all dcd 10 nvt_rescale_298.dcd
run 2000000

unfix 2
undump dcd_1
dump dcd_2 all dcd 100 nvt_298.dcd
run 10000

Dear lammps users,

Thank you Axel for correcting the cutting water method.
I recut my water, but step 0 still happens nan in pressure and the
simulation crush. Does anyone has any idea? I attach my simulation below,
this simulation is after cutting the water into desire shape, I merged the
water and pyrophyllite together. I make two pyrophyllite clay particles
which sandwich the water.

could be you have overlapping atoms, possibly due to a miscalculation when
merging structure, or due to periodic boundaries.

try running the command:

delete_atoms overlap 0.5 all all

and see if any atoms are being deleted. if yes, you have atoms that are
very close. you can find out which ones by writing out a data file and
comparing with the original data file.

axel.

Yes, there are atoms that are too close and it happens in the water. The simulation deletes 1845 atoms from the water molecules, though I calculate the number of water angles(615) and water bonds(1230) that has been cut, and the number of bonds and angles do meet the water molecules number, I am still concerned that whether the simulation has delete incomplete water molecules and let HO or only O or H working in my simulation. Besides, does this phenomenon happen because my water isn’t right?
I happen to see that my water crushes during the simulation (cannot find pppm error) (I remember Axel does tell me about this error before about the aggressive set up in my command, it was solved when I am using SPCE water model, but it occurs in here while this water model is tip3p), but I do not put it in mind because the system has been minimized. I think I should not ignore this problem, so below is the input script of my water.

Thank you!

Miao-Chun

units real
dimension 3
boundary p p p
atom_style full
read_data system.data

alternate way to set charges - besides manually in data file

group ox type 1
group hy type 2
set group ox charge -0.830
set group hy charge 0.415

Flexible SPC/E Potential Parameters

Zhang et al., Fluid Phase Equilibria, 262 (2007) 210-216

pair_style lj/cut/coul/long 13.0
pair_modify mix arithmetic tail no
pair_coeff 1 1 0.152 3.15
pair_coeff 2 2 0.046 0.4

bond_style harmonic
bond_coeff 1 450.0 0.9572
angle_style harmonic
angle_coeff 1 34.933 104.52
kspace_style pppm 1.0e-5

velocity all create 1.0 12345689 dist uniform
fix shake all shake 1.0e-4 100 0 b 1 a 1
fix 2 all npt temp 298.0 298.0 100.0 iso 0.0 0.0 100.0

compute 3 all property/atom fx fy fz
compute 4 all reduce sum c_3[1] c_3[2] c_3[3]

neighbor 2.0 bin
neigh_modify delay 2 every 1 check yes
thermo 10
thermo_style custom step temp press vol pe etotal c_4[1] c_4[2] c_4[3]
thermo_modify norm no flush yes

#dumps
dump waterdump all dcd 100 npt_water.dcd

#run variables
timestep 0.25
run 100000

Yes, there are atoms that are too close and it happens in the water. The
simulation deletes 1845 atoms from the water molecules, though I calculate
the number of water angles(615) and water bonds(1230) that has been cut,
and the number of bonds and angles do meet the water molecules number, I am
still concerned that whether the simulation has delete incomplete water
molecules and let HO or only O or H working in my simulation. Besides, does
this phenomenon happen because my water isn't right?

​i don't know. i don't have a crystal ball. but extrapolating from previous
experience, i'd say there is a high probability that something went wrong
when merging your two systems. this can often be found out with
visualization. remember to visualize the simulation cell as well.​

​figuring out how to properly set up a complex system can be challenging.
it is commendable that you are trying to be careful. however, thes​e things
are extremely difficult to debug from remote. so you really have to check
every detail by yourself. perhaps first practice with something simpler and
smaller.

axel.

I happen to see that my water crushes during the simulation (cannot find

Sorry I made a mistake here. I was too tired and I simulate the wrong file. Please forgive me of not paying attention.
Actually, there are two atoms that overlapped with each other. I write out the data file, but I don’t have an idea how to check these two files together.
Am I going to look over it line by line or is there a faster way to work with?
Thank you.

Miao-Chun

Dear users,

I solve the problem of how to find out the delete atoms.
However, I still occurs the error at step 0 with nan in pressure and the simulation crush. Is there any other probabilities?
Thank you!

Miao-Chun