Solid-liquid interface

Abdul_Kareem · February 28, 2023, 5:16pm

I am trying to simulate a solid silicate and a water droplet on its surface I have attached the input files. The silicate and water structures are equilibrated separately. When trying to merge the structures and equilibrate the lammps running stops at the first step.
I have attached the log below

lmp -in input.lammps
Invalid MIT-MAGIC-COOKIE-1 keyLAMMPS (23 Jun 2022 - Update 2)
OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98)
using 1 OpenMP thread(s) per MPI task
Reading data file …
orthogonal box = (-0.58 -0.68 -0.66) to (39.32 39.22 39.24)
1 by 1 by 1 MPI processor grid
reading atoms …
6495 atoms
scanning bonds …
2 = max bonds/atom
scanning angles …
1 = max angles/atom
reading bonds …
4330 bonds
reading angles …
2165 angles
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.001 seconds
read_data CPU = 0.019 seconds
Reading data file …
orthogonal box = (-174.67 -174.67 -4.0125) to (39.32 39.22 39.24)
1 by 1 by 1 MPI processor grid
reading atoms …
10800 atoms
reading velocities …
10800 velocities
Finding 1-2 1-3 1-4 neighbors …
special bond factors lj: 0 0 0
special bond factors coul: 0 0 0
2 = max # of 1-2 neighbors
1 = max # of 1-3 neighbors
1 = max # of 1-4 neighbors
2 = max # of special neighbors
special bonds CPU = 0.002 seconds
read_data CPU = 0.041 seconds
Reading tersoff potential file SiO.tersoff with DATE: 2010-08-16
WARNING: Converting tersoff potential in metal units to real units (src/potential_file_reader.cpp:289)

Info-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info
Printed on Tue Feb 28 22:46:34 2023

Coeff status information:

Pair Coeffs:
1 1: is set
1 2: is set
1 3: is set
1 4: is set
2 2: is set
2 3: is set
2 4: is set
3 3: is set
3 4: is set
4 4: is set

Bond Coeffs:
1: is set

Angle Coeffs:
1: is set

Info-Info-Info-Info-Info-Info-Info-Info-Info-Info-Info

6495 atoms in group H2O
System init for delete_atoms …
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:239)
Neighbor list info …
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 36 36 8
4 neighbor lists, perpetual/occasional/extra = 3 1 0
(1) command delete_atoms, occasional, copy from (4)
attributes: full, newton on
pair build: copy
stencil: none
bin: none
(2) pair lj/cut/coul/cut, perpetual, half/full from (4)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(3) pair tersoff, perpetual, skip from (4)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(4) neighbor class addition, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Neighbor list info …
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 36 36 8
0 neighbor lists, perpetual/occasional/extra = 0 0 0
WARNING: Ignoring ‘compress yes’ for molecular system (src/delete_atoms.cpp:140)
Deleted 4497 atoms, new total = 12798
Deleted 2998 bonds, new total = 1332
Deleted 1499 angles, new total = 666
10800 atoms in group slab
Finding SHAKE clusters …
0 = # of size 2 clusters
0 = # of size 3 clusters
0 = # of size 4 clusters
666 = # of frozen angles
find clusters CPU = 0.001 seconds
WARNING: Using a manybody potential with bonds/angles/dihedrals and special_bond exclusions (src/pair.cpp:239)
Neighbor list info …
update every 1 steps, delay 10 steps, check yes
max neighbors/atom: 2000, page size: 100000
master list distance cutoff = 12
ghost atom cutoff = 12
binsize = 6, bins = 36 36 8
3 neighbor lists, perpetual/occasional/extra = 3 0 0
(1) pair lj/cut/coul/cut, perpetual, half/full from (3)
attributes: half, newton on
pair build: halffull/newton
stencil: none
bin: none
(2) pair tersoff, perpetual, skip from (3)
attributes: full, newton on
pair build: skip
stencil: none
bin: none
(3) neighbor class addition, perpetual
attributes: full, newton on
pair build: full/bin
stencil: full/bin/3d
bin: standard
Setting up Verlet run …
Unit style : real
Current step : 0
Time step : 1

how to solve it. the running stops here. I have attached the input script and the data files.
data.final_sio2 (1.5 MB)
data.water (471.4 KB)
input.lammps (1.5 KB)
SiO.tersoff (1.9 KB)

patherlkd · February 28, 2023, 5:50pm

Abdul:

When I run your input.lammps on my machine with specs:

(base) ldavis@dps112:/local/ldavis$ lscpu
Architecture: x86_64
CPU op-mode(s): 32-bit, 64-bit
Address sizes: 46 bits physical, 48 bits virtual
Byte Order: Little Endian
CPU(s): 64
On-line CPU(s) list: 0-63
Vendor ID: GenuineIntel
Model name: Intel(R) Xeon(R) Gold 6242 CPU @ 2.80GHz

using LAMMPs 8 Feb 2023
I get the following:

Per MPI rank memory allocation (min/avg/max) = 36.07 | 36.07 | 36.07 Mbytes
Step Temp E_pair E_mol TotEng Press
0 2.7140548e+08 -2007893.3 0 9.8120703e+09 2.2660601e+08
1000 -nan -nan 0 -nan -nan
2000 -nan -nan 0 -nan -nan
3000 -nan -nan 0 -nan -nan
4000 -nan -nan 0 -nan -nan
5000 -nan -nan 0 -nan -nan
6000 -nan -nan 0 -nan -nan
7000 -nan -nan 0 -nan -nan
8000 -nan -nan 0 -nan -nan
9000 -nan -nan 0 -nan -nan
10000 -nan -nan 0 -nan -nan
11000 -nan -nan 0 -nan -nan
12000 -nan -nan 0 -nan -nan
13000 -nan -nan 0 -nan -nan
14000 -nan -nan 0 -nan -nan
15000 -nan -nan 0 -nan -nan
16000 -nan -nan 0 -nan -nan
17000 -nan -nan 0 -nan -nan
18000 -nan -nan 0 -nan -nan
19000 -nan -nan 0 -nan -nan
20000 -nan -nan 0 -nan -nan
21000 -nan -nan 0 -nan -nan
22000 -nan -nan 0 -nan -nan
23000 -nan -nan 0 -nan -nan
24000 -nan -nan 0 -nan -nan
25000 -nan -nan 0 -nan -nan

It keeps running. But I ctrl-c’ed to stop it because I am seeing nans in the thermo output.
My guess is your initial conditions are wrong.

Best,

L

patherlkd · February 28, 2023, 5:53pm

Ways to debug:

Simplify! You have many things going on in your simulation. A good starting point is to remove complexity either by simulating 1, 2, and then a few particles AND/OR simplifying the interactions and debug each interaction.
Visualize your initial conditions and verify it is what you expect.
Check your fixes are not conflicting and/or wrong in of themselves…

patherlkd · February 28, 2023, 5:58pm

A big hint is that your Temperature at t=0 is 2.7140548e+08 Kelvin. For reference the surface of the sun is 5,778 K.

So your temperature is definitely wrong and you should check all your simulation parameters in REAL units.

Best,

L

Abdul_Kareem · February 28, 2023, 6:09pm

Thank you

akohlmey · March 1, 2023, 12:12pm

Most likely, you equilibrated your SiO2 system with metal units but are now reading the data file in real units. That will scale all velocities by a factor of 1000, because the time unit changes from picoseconds to femtoseconds.

Moreover, it is not likely that your setup will work. The SiO.tersoff potential you are using is only suitable for bulk systems and your force field setup with hybrid/overlay may be double counting interactions and does not seem very suitable for SiO2-water interactions. Treating the SiO2 oxygen atoms like water oxygen atoms is most certainly not suitable, specifically since in the SIO.tersoff model charges are implicit and in lj/cut/coul/cut they are explicit.

When you cut SiO2 you usually have “dangling bonds” (i.e. open valences or oxygen atoms) and those typically would be saturated with hydrogen atoms. I don’t see any provisions for that here.

This all looks like you are trying to teach yourself how to set up such a calculation and are way in over your head without competent supervision. This is a recipe for lots of frustration, wasted time, and bogus simulation results. You should get competent in-person help.

Abdul_Kareem · March 8, 2023, 6:08pm

Hi Akohlmey, Thank you for your valuable suggestion. Instead of SiO2, now I tried using graphene layers trying to reproduce the results from the paper " Molecular dynamics simulations of the contact angle between water droplets and graphite surfaces". I have used almost the same parameters as mentioned in this paper except using run_style rRESPA and using kspace_style pppm which is not in support with charmm forcefield. while running, I am getting the same problem which I faced in Sio2. The simulation didn’t start. I am not able to figure it out. searching for shows that it may due to the overlap of atoms. using delete_atoms command I checked whether there is any but there was none. I have attached the input files for your reference.
input.lammps (1.6 KB)
graphene.data (2.5 MB)
water_512.data (246.7 KB)

akohlmey · March 8, 2023, 7:05pm

What do you mean by that?

When I run your input, I get an error that an atom has moved too far. This is with the current version of LAMMPS (8 Feb 2023).

akohlmey · March 8, 2023, 7:16pm

You have a serious problem with how you do things.

The box sizes for your two data files are very different, but in your water box you have non-zero image flags. If you read the second data file, the box dimensions get changed to accommodate both systems, but the image flags or their corresponding coordinates are not corrected for that, which is equivalent to a massive displacement, which then creates serious problems with fix shake.

The best way to do this is to either merge the two data files externally with some scripting either written by yourself or using VMD and Tcl scripting with TopoTools.

Another option would be to use only the data file for graphene, but add the extra atom, bond, and angle types for your water and create the water with LAMMPS using a molecule file by filling a region. If you have a sufficiently recent LAMMPS version you can use the overlap option for random placement of water molecules and I have found recently empirically that this works quite well for water with a setting of 1.333 angstrom.

akohlmey · March 8, 2023, 8:03pm

As an example, here is an input based on yours showcasing how to add a water droplet to your graphene system. If you want a film on a slab, you need to input the graphene atoms differently.
For simplicity, I am also keeping those atoms immobile.

spce.mol (618 Bytes)

units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic
dihedral_style opls
improper_style harmonic
special_bonds lj/coul 0.0 0.0 0.5
pair_style    lj/charmm/coul/charmm 10 12 
#kspace_style pppm 1.0e-4 

read_data graphene.data extra/atom/types 2 extra/bond/types 1 extra/angle/types 1

mass 2 15.9994
mass 3 1.008

pair_coeff * *  0.00000  0.00000
pair_coeff 1 1  0.066047 3.4
pair_coeff 1 2  0.0478   3.581
pair_coeff 2 2  0.15533 3.16557

bond_coeff 1 469 1.4
bond_coeff 2 1000.0 1.0

angle_coeff 1 63 120
angle_coeff 2 10000.0 109.47

dihedral_coeff 1 0 7.25 0 0
improper_coeff 1 5 180

region liquid block 30.0 70.0 30.0 70.0 2.0 27.0
molecule water spce.mol toff 1 boff 1 aoff 1
create_atoms 0 random 512 23462134 liquid mol water 85342 overlap 1.33333

group H2O type 2 3 
group slab type 1

fix             walls slab setforce 0.0 0.0 0.0
fix		mynvt all nvt temp 10.0 10.0 10
fix 		myshk H2O shake 1.0e-3 100 10000 b 2 a 2

velocity H2O create 100.0 465135

thermo 100
minimize  1.0e-4 1.0e-4 10000000 10000

reset_timestep 0

write_data minimized.data pair ij

timestep       1 	

dump 		dp1 all atom 1000 dump_graphene_water_nvt.lammpstrj
thermo		1000

run		5000
write_data	data.graphene_water

Abdul_Kareem · March 11, 2023, 6:30am

Thank You Axel Kohlmeyer.