Dear All,
I am trying to simulate contact angle of water on a quartz surface. I have created and equilibrated a quartz slab using tersoff potential. I have also created a water cluster and equilibrated it using the SPC model. I combined both these geometries using a Python script. This combined geometry is the input file for my NVT simulations to calculate contact angle. My input file is as follows:
units real
dimension 3
processors * * *
boundary p p f
read data
atom_style full
read_data contact.dat
group water type 4 5
#potentials
pair_style hybrid lj/cut/coul/long 10.0 lj/cut/coul/long 10.0 lj/cut/coul/long 10.0 lj/cut/coul/long 10 lj/cut/coul/long 10 lj/cut/c$
pair_coeff * * tersoff SiO2.tersoff Si O O NULL NULL #Si-O in slab
pair_coeff 4 5 lj/cut/coul/long 1 0.0 0.0 #H-O(water)
pair_coeff 4 4 lj/cut/coul/long 2 0.0 0.0 #H-H(water)
pair_coeff 5 5 lj/cut/coul/long 3 0.15535 3.166 #O-O(water)
pair_coeff 1 5 lj/cut/coul/long 4 0.34844479 2.633 #Si-O interface
pair_coeff 3 5 lj/cut/coul/long 5 8.763 1.8 #O-O interface
pair_coeff 1 4 lj/cut/coul/long 6 0.0 0.0
pair_coeff 2 4 lj/cut/coul/long 7 0.0 0.0
pair_coeff 2 5 lj/cut/coul/long 8 0.0 0.0
pair_coeff 3 4 lj/cut/coul/long 9 0.0 0.0
kspace_style pppm 1.0e-5
kspace_modify slab 3.0
bond_style harmonic
bond_coeff 1 0.0 1.149
angle_style harmonic
angle_coeff 1 0 109.47
neighbor 2.0 bin
neigh_modify delay 10 check yes
timestep 2.0
#fixes and computes
fix constrain water shake 1.0e-4 100 0 b 1 a 1
fix zwalls all wall/reflect zlo EDGE zhi 5.8
compute cont water chunk/atom bin/cylinder y center 25 0 0 5.6 14 3 discard yes
fix 1 water ave/chunk 100 1000 100000 cont density/mass ave running file density.profile”
fix integrate all nvt temp 300.0 300.0 100.0
thermo 10000
thermo_style custom step temp lx ly lz
dump 2 all custom 10000 contact.lammpstrj element x y z
dump_modify 2 element Si O O H O
run 2000000
My data file contains the coordinates of all the atoms. It also has the bonds and angles data for the water molecule as given by VMD ‘writelammpsdata’ command. However, I get the following error:
ERROR on proc 0: Shake atoms 374 366 382 missing on proc 0 at step 0 (…/fix_shake.cpp:546)
Last command: run 2000000
ERROR on proc 1: Shake atoms 756 748 764 missing on proc 1 at step 0 (…/fix_shake.cpp:546)
Last command: run 2000000
ERROR on proc 2: Shake atoms 757 749 765 missing on proc 2 at step 0 (…/fix_shake.cpp:546)
difficult to say anything without knowing more details. most likely your data file is inconsistent.
Here are the steps that I followed to create the input file:
-
Created an unit cell of quartz, replicated (using replicate command) it in 3 dimensions, minimised energy and equilibrated at 300K, 1atm using tersoff potential (tersoff potential does not require bond/angle information).
-
Created a water molecule, replicated (using replicate command) it in 3 dimensions, minimised energy and equilibrated at 300K, 1atm using SPC model.
-
To get the bond and angle information of the replicated water cluster, I used writelammpsdata command in VMD
-
Now, my aim is to place the water cluster on the quartz slab and run a NVT simulation of the composite geometry. For that I wrote a Python script that modifies the output trajectory files by bringing the COMs of both the geometries to (0,0,0) and then puts the water cluster on the slab by adjusting the z coordinates. I simply copied the bond and angle information of the water cluster generated by VMD into this input file (I think the spatial shifting of the cluster will not affect the bond and angle data). All other interactions in the system are either through tersoff potential or non-bonded LJ potential that does not require bond/angle data.
Hope this helps.
Here are the steps that I followed to create the input file:
-
Created an unit cell of quartz, replicated (using replicate command) it in 3 dimensions, minimised energy and equilibrated at 300K, 1atm using tersoff potential (tersoff potential does not require bond/angle information).
-
Created a water molecule, replicated (using replicate command) it in 3 dimensions, minimised energy and equilibrated at 300K, 1atm using SPC model.
-
To get the bond and angle information of the replicated water cluster, I used writelammpsdata command in VMD
-
Now, my aim is to place the water cluster on the quartz slab and run a NVT simulation of the composite geometry. For that I wrote a Python script that modifies the output trajectory files by bringing the COMs of both the geometries to (0,0,0) and then puts the water cluster on the slab by adjusting the z coordinates. I simply copied the bond and angle information of the water cluster generated by VMD into this input file (I think the spatial shifting of the cluster will not affect the bond and angle data). All other interactions in the system are either through tersoff potential or non-bonded LJ potential that does not require bond/angle data.
there are two potential issues with your step 4:
depending on how you translate the coordinates, you may have incorrect atom positions due to inconsistent image flags (check your output)
and depending on how you handle atom IDs when merging the two subsystems, the atom IDs describing the bonds may point to the wrong atoms without adjustments.
So you need to verify that your data file is correct. I suggest repeating the procedure with a very small number of atoms, so you can check it manually.
You may have more luck by just reversing the order of which subsystem you add to which.
An alternative would be to use change_box and displace_atoms in LAMMPS at the end of the preparation for each of the subsystems so that each part of the system is at its final position and the box dimensions are the same. then you can use read_data twice and use proper reserved types on the first read_data and suitable offsets on the second read.
there really isn’t much else that can be done from remote.
you just have to make sure that your data file is consistent and correct.
axel.
Sir, I thank you for your valuable suggestions. I went through the data file created by VMD and did find some erroneous bond and angle connections. I used the method suggested by you to create the composite geometry (used create_box and then used read_data twice with proper offsets). In order to get a correct topology for the water cluster, I used the data.SPCE file provided in the LAMMPS examples. I shifted the slab instead of the water cluster so as to not change the H,O coordinates in water. But unfortunately, I get the same error. I even tried increasing the cutoff distance for neighbour list creation and also tried creating the neighbour list more frequently using the neigh_modify command. But nothing seems to work. Can you suggest any other plausible solution to this problem? Are the non periodic boundary conditions or the rigid wall that is being used to simulate the interface creating this issue ? Here is the input file:
units real
dimension 3
processors * * *
boundary p p f
read data
atom_style full
region 1 block -32 65 -32 65 -15 37
create_box 5 1 bond/types 1 angle/types 1 dihedral/types 0 improper/types 0 extra/bond/per/atom 4 extra/angle/per/atom 4
read_data water_1.dat add merge
read_data slab.dat add append offset 2 0 0 0 0 shift -30 -30 -13
group water type 1 2
#potentials
pair_style hybrid lj/cut/coul/long 10.0 lj/cut/coul/long 10.0 lj/cut/coul/long 10.0 lj/cut/coul/long 10 lj/cut/coul/long 10 tersoff
pair_coeff * * tersoff SiO2.tersoff NULL NULL Si O O
pair_coeff 1 2 lj/cut/coul/long 1 0.0 0.0 #H-O(water)
pair_coeff 2 2 lj/cut/coul/long 2 0.0 0.0 #H-H(water)
pair_coeff 1 1 lj/cut/coul/long 3 0.15535 3.166 #O-O(water)
pair_coeff 1 3 lj/cut/coul/long 4 0.34844479 2.633 #Si-O interface
pair_coeff 1 5 lj/cut/coul/long 5 8.763 1.8 #O-O interface
pair_coeff 2 3 lj/cut/coul/long 6 0.0 0.0
pair_coeff 2 4 lj/cut/coul/long 7 0.0 0.0
pair_coeff 1 4 lj/cut/coul/long 8 0.0 0.0
pair_coeff 2 5 lj/cut/coul/long 9 0.0 0.0
kspace_style pppm 1.0e-5
kspace_modify slab 3.0
bond_style harmonic
bond_coeff 1 0.0 1.149
angle_style harmonic
angle_coeff 1 0 109.47
neighbor 5.0 bin
neigh_modify delay 10 check yes
timestep 1.0
dump 2 all custom 10000 contact.lammpstrj element x y z
dump_modify 2 element Si O O H O
#fixes and computes
fix constrain water shake 1.0e-4 200 0 b 1 a 1
fix zwalls all wall/reflect zlo EDGE zhi -1
compute cont water chunk/atom bin/cylinder y center 25 0 0 12 20 5 discard yes
fix 1 water ave/chunk 100 1000 100000 cont density/mass ave running file density.profile”
fix integrate all nvt temp 300.0 300.0 100.0
thermo 10000
thermo_style custom step temp lx ly lz
run 2000000
one problem when using the data file from the LAMMPS examples is that it was equilibrated for fully periodic boundaries and thus has atoms with image flag values != 0.
that would be a problem when loading such a data file into a system with (partially) non-periodic boundaries, as it would read atoms that are outside the box and thus should be ignored, or if they get included by accident, they will have incorrect positions.
this problem is made worse by using read_data twice, because that will change the box (and thus the offset used for the image flags). so you will have to process the data file(s) before you can load them to have box dimensions where all image flags for non-periodic dimensions are zero. it may be easier in general to just have the box and positions adjusted in a separate input file (using change_box and displace_atoms and then write out a new data file with write data) so that the box dimensions for all read in data files are the same.
I would also make adjustments to atom/bond/angle/dihedral/etc types and just write out data files with “nocoeff” to skip writing forcefield data.
axel.
Sir,
I did the following after your last email:
- Unwrapped the coordinates (using the image flags) of all the water atoms in data.spce file provided within LAMMPS examples. I found out the maximum and minimum value of the z-coordinates (because z direction will finally have the non-periodic boundary) and created a box large enough to contain all the atoms(both the water as well as the slab atoms). For example, if the maximum value of unwrapped z coordinate comes out as 36, I set zhi as 37. xhi,xlo,yhi,ylo,zlo are determined from the position of the slab. I then made the z image flags of all water atoms 0 with the set image command. I did not alter the x,y image flags as those directions will have periodic boundaries. The input script is :
units real
atom_style full
read_data water.dat
change_box all x final -30.199999999999999 63 y final -30.199999999999999 63 z final -20 37 units box
set type 1*2 image NULL NULL 0
write_data data_water.dat
-
Created a simulation box with the same dimensions for the slab. Shifted the slab atoms with displace_atoms by sufficient amounts so that the water cluster sits on it in the composite geometry. I ensured that there is no overlap and there is a gap of 1ang between the slab and water cluster. Again ensured that the image flags in the z dimension are 0 with the set image command. Input script:
units real
atom_style full
read_data slab.dat
change_box all x delta -30 -30 y delta -30 -30 z final -20 37 units box
displace_atoms all move -30 -30 -13 units box
set type 1*3 image NULL NULL 0
write_data data_slab.dat
-
Used the data files generated by these write_data commands in my final simulation. As a big enough simulation box is already created and the slab atoms have already been moved to their required positions, I did not use the create_box command in the final simulation. shift keyword was also not required with the read_data commands. Here is the geometry creation part of the script:
read_data data_water.dat extra/atom/types 3 extra/bond/types 0 extra/angle/types 0 extra/bond/per/atom 0 extra/angle/per/atom 0
read_data data_slab_1.dat add append offset 2 0 0 0 0
But unfortunately the same error of missing atoms crop up.
Could I interpret your suggestions correctly ? Can manual removal of the atoms that are creating trouble be of any help ? Also, is there a way I can visualise the composite geometry before the simulation starts ? I think creating a restart file will not be of any help as the problem occurs in the 0th step. I am a complete newbie in LAMMPS and I am banking on this conversation to get this technical problem solved. I thank you once again for your valuable suggestions.
sorry, but this is going far too far now. I am not your adviser and what you (very obviously) need is advice, tutoring, and supervision in person. this is the job of your adviser/supervisor and not the purpose of this mailing list. I can give some pointers, but not spoon-feed and individually tutor you to solve your problems.
besides, if you are completely new to the topic of simulations, you should not even be trying to do a simulation of such a complex compound system, but rather develop and hone your skills with simpler problems. again, it is the job of your adviser/supervisor to organize this and arrange for proper training and tutoring before you can be asked to take on a more complex problem.
I have given as much and as detailed advice and spent as much time as I can afford on helping you, but there is a limit which you have now reached.
And I have already taken into account your status as a beginner and gone much farther and provided more detail and specific instructions than I would have done for somebody more experienced. As for solving technical problems, most can be resolved on your own by paying attention to detail, carefully studying the documentation, observing and checking everything at every step, and using common sense.
Axel.
I completely understand. Thank you for giving so much time to my problems.