I am trying to hands on the slab + water interactions from the Silicon and water interaction calculation.
The dump file only shows silicon atoms in the simulation box.
- It doesn’t generate water molecule in the box, only silicon atoms are filled all the box.
- Even the generate silicon doesn’t look like a slab as it contract within the box through steps, not bottom fixed.
I am guessing the region setting or defining the water molecule (tip3p.mol) is wrong, but with my own beginner’s view I cannot notice which part makes the calculation doesn’t work as I intended.
the version: 02Aug23
units real
atom_style full
boundary p p p
# Define the simulation box (ensure it covers the water and silicon regions)
region box block -5 5 -5 5 0 50
create_box 3 box bond/types 1 angle/types 1 &
extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
# Define atomic masses
mass 1 28.0855 # Silicon
mass 2 15.9994 # Oxygen (Water)
mass 3 1.008 # Hydrogen (Water)
# Silicon slab setup (using a diamond lattice for silicon)
lattice diamond 5.431
region silicon block -5 5 -5 5 0 10
create_atoms 1 region silicon
# Define hybrid pair style for both Lennard-Jones and Stillinger-Weber potentials
pair_style hybrid sw lj/cut/coul/cut 8.0
# Define Stillinger-Weber potential for silicon
pair_coeff * * sw Si.sw Si NULL NULL # Silicon interaction via SW potential
# Define Lennard-Jones for water and silicon-water interactions
pair_coeff 1 2 lj/cut/coul/cut 0.1521 3.1507 # Silicon-Oxygen interaction
pair_coeff 2 2 lj/cut/coul/cut 0.0 1.0 # Oxygen-Oxygen interaction in TIP3P
pair_coeff 2 3 lj/cut/coul/cut 0.0 1.0 # Oxygen-Hydrogen interaction in TIP3P
pair_coeff 1 3 lj/cut/coul/cut 0.0 1.0 # Silicon-Hydrogen interaction
pair_coeff 3 3 lj/cut/coul/cut 0.0 1.0 # Hydrogen-Hydrogen interaction
# Define the water molecule template from a file
molecule water tip3p.mol
# Define the water region above the silicon slab (expanded for more space)
region water block -8 8 -8 8 10 50 units box
# Create water molecules randomly in the 'water' region
create_atoms 0 random 100 34564 water mol water 25367 overlap 5.0
# Print the number of atoms after water molecule creation
variable natoms equal "count(all)"
print "Number of atoms after water molecule creation: ${natoms}"
# Group atoms by type
group silicon type 1
group oxygen type 2
group hydrogen type 3
# Print the number of atoms for each type
variable nSi equal count(silicon)
variable nO equal count(oxygen)
variable nH equal count(hydrogen)
print "Number of Silicon atoms: ${nSi}"
print "Number of Oxygen atoms: ${nO}"
print "Number of Hydrogen atoms: ${nH}"
# Apply a fix for time integration (NVT ensemble)
fix 1 all nvt temp 300.0 300.0 100.0
# Dump atomic trajectory every 1000 steps, include id and type
dump 1 all custom 1000 dump.silicon_water.lammpstrj id type xs ys zs
dump_modify 1 sort id
# Run the simulation for 20000 steps
run 20000
# Write final configuration
write_data final_silicon_water.data nocoeff
Part of the printed output
...
38 Number of atoms after water molecule creation: 230
39 230 atoms in group silicon
40 0 atoms in group oxygen
41 0 atoms in group hydrogen
42 Number of Silicon atoms: 230
43 Number of Oxygen atoms: 0
44 Number of Hydrogen atoms: 0
45 Neighbor list info ...
...
tip3p.mol:
Types
1 1 # O
2 2 # H
3 2 # H
Charges
1 -0.834
2 0.417
3 0.417
Bonds
1 1 1 2
2 1 1 3
Angles
1 1 2 1 3
Shake Flags
1 1
2 1
3 1
Shake Atoms
1 1 2 3
2 1 2 3
3 1 2 3
Shake Bond Types
1 1 1 1
2 1 1 1
3 1 1 1
Special Bond Counts
1 2 0 0
2 1 1 0
3 1 1 0
Special Bonds
1 2 3
2 1 3
3 1 2
This is not correct. With this setting the water molecules only interact via charges. Please see, e.g. 8.4.3. TIP3P water model — LAMMPS documentation or the original publication describing TIP3P water.
Check the log after this command. I would suspect that overlap 5.0 is too large so all proposed insertions are rejected. Something like 0.5 should suffice.
This is probably not what you want because of your lattice command and the region command with units box below. So you likely want to use “units box” here as well.
The following things are missing:
- a bond style and corresponding parameters
- an angle style and corresponding parameters
- fix shake to constrain the water molecules to have a rigid structure.
Managed the dump file to have slab of silicon and other atoms above the silicon slab. However, the molecules on the silicon is not H2O they are Si and O gases.
Problem: water molecule should be found above the Si slab. However, instead of H2O I found Si and O atoms.
Is this because of the label in tip3p.mol file or other error in tip3p.mol file?
units real
atom_style full
boundary p p p
# Define the simulation box (ensure it covers the water and silicon regions)
region box block -5 5 -5 5 0 50
create_box 3 box bond/types 1 angle/types 1 &
extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2
# Define atomic masses
mass 1 28.0855 # Silicon
mass 2 15.9994 # Oxygen (Water)
mass 3 1.008 # Hydrogen (Water)
# Silicon slab setup (using a diamond lattice for silicon)
lattice diamond 5.431
#region silicon block -5 5 -5 5 0 10
region silicon block -5 5 -5 5 0 10 units box
create_atoms 1 region silicon
# Define hybrid pair style for both Lennard-Jones and Stillinger-Weber potentials
pair_style hybrid sw lj/cut/coul/cut 8.0
# Define Stillinger-Weber potential for silicon
pair_coeff * * sw Si.sw Si NULL NULL # Silicon interaction via SW potential
# Define Lennard-Jones for water and silicon-water interactions
pair_coeff 1 2 lj/cut/coul/cut 0.1521 3.1507 # Silicon-Oxygen interaction
pair_coeff 2 2 lj/cut/coul/cut 0.1521 3.1507 # Oxygen-Oxygen interaction in TIP3P
pair_coeff 2 3 lj/cut/coul/cut 0.0 1.0 # Oxygen-Hydrogen interaction in TIP3P
pair_coeff 1 3 lj/cut/coul/cut 0.0 1.0 # Silicon-Hydrogen interaction
pair_coeff 3 3 lj/cut/coul/cut 0.0 1.0 # Hydrogen-Hydrogen interaction
bond_style zero
bond_coeff 1 0.9574
angle_style zero
angle_coeff 1 104.52
# Define the water molecule template from a file
molecule water tip3p.mol
# Apply SHAKE algorithm to keep water molecules rigid
fix water_shake all shake 0.0001 20 0 b 1 a 1
# Define the water region above the silicon slab (expanded for more space)
region water block -8 8 -8 8 10 50 units box
# Create water molecules randomly in the 'water' region
create_atoms 0 random 100 34564 water mol water 25367 overlap 0.5
# Print the number of atoms after water molecule creation
variable natoms equal "count(all)"
print "Number of atoms after water molecule creation: ${natoms}"
# Group atoms by type
group silicon type 1
group oxygen type 2
group hydrogen type 3
# Print the number of atoms for each type
variable nSi equal count(silicon)
variable nO equal count(oxygen)
variable nH equal count(hydrogen)
print "Number of Silicon atoms: ${nSi}"
print "Number of Oxygen atoms: ${nO}"
print "Number of Hydrogen atoms: ${nH}"
# Set the skin distance to 3.0
neighbor 1.0 bin
# Neighbour list modifications to prevent missing atoms
neigh_modify every 1 delay 0 check yes
# Minimise the system before starting the MD run
minimize 1.0e-4 1.0e-6 100 1000
# Main NVT run directly at 300 K
fix 1 all nvt temp 300.0 300.0 100.0
timestep 0.1 # Reduced timestep to avoid missing bond atoms
# Dump atomic trajectory every 1000 steps, include id and type
dump 1 all custom 1000 dump.silicon_water.lammpstrj id type xs ys zs
dump_modify 1 sort id
# Run the simulation for 20000 steps
run 20000
# Write final configuration
write_data final_silicon_water.data nocoeff
partial of log:
...
44 Number of atoms after water molecule creation: 167
45 89 atoms in group silicon
46 78 atoms in group oxygen
47 0 atoms in group hydrogen
48 Number of Silicon atoms: 89
49 Number of Oxygen atoms: 78
50 Number of Hydrogen atoms: 0
51 WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (../pair.cpp:242)
52 WARNING: Using fix shake with minimization.
53 Substituting constraints with harmonic restraint forces using kbond=1987 (../fix_shake.cpp:360)
...
129 Per MPI rank memory allocation (min/avg/max) = 9.45 | 9.487 | 9.718 Mbytes
130 Step Temp E_pair E_mol TotEng Press
131 13 0 -19930.904 0 -19930.904 -97502.957
132 ERROR on proc 3: Bond atoms 132 133 missing on proc 3 at step 6975 (../ntopo_bond_partial.cpp:60)
133 Last command: run 20000
134 Abort(1) on node 3 (rank 3 in comm 0): application called MPI_Abort(MPI_COMM_WORLD, 1) - process 3
It is because you don’t know what you are doing and are not paying proper attention to advice been given. I don’t like to interact with such people, you have to find someone else to sort out your issues.
Thank you very Axel much for your guidance.
*As I noticed the tip3p.mol’s atom types and the input’s atom types (id/label) were miss matched.