GROMACS version: 2022.4
GROMACS modification: No
Dear GROMACS users,
I’m trying to simulate thin liquid film on a solid substrate and observe the breakup (or dewetting) of the thin film. Because I’m only interested in the liquid, I would like to freeze the solid substrate to simplify the analysis. The interaction between the solid and liquid is still left on since is is needed to break the film. To do that, I define two groups, one named “Fluid” with all the liquid atoms, and one named “bWall” with all the solid atoms. I then add the following two lines to the system.mdp file
and try to keep the temperature of the solid substrate to zero (so it does not move).
I then use grompp to generate initial structure gmx grompp -f system.mdp -c system.pdb -p system.top -o system.tpr -n system.ndx
and run the MD simulation by gmx mdrun -deffnm system -s system.tpr -g system.log -o system.trr -c system.gro -nt 30
You can find the .top file system.top (992 Bytes)
and the .mdp file system.mdp (7.0 KB)
here.
The argon.itp file is
[ moleculetype ]
Ar 1
[ atoms ]
1 Ar 1 Ar Ar 1 0
and the platinum.itp file is
[ moleculetype ]
Pt 1
[ atoms ]
1 Pt 1 Pt Pt 1 0
The simulation runs without errors, however the solid atoms move. But they shouldn’t, right? Did I do something wrong? I’m new to Gromacs and I’m afraid I might have done something very stupid… Thanks a lot in advance!
Thanks for the idea. I tried generating index group file with gmx make_ndx -f system.pdb -o system.ndx and the problem persists. I think it’s not an index group file problem then.
Thanks for the idea. I tried and the platinum atoms still move. Not like they fly around, just oscillates. Which gets me to think maybe this is just an error in the output file? I’ll elaborate in the message below.
After looking carefully into the position data I noticed that i) the displacement of platinum atoms are all 3.92 angstrom, which is exactly the lattice parameter (side length of cubic cell) ii) the displacement only happens in x and y direction.
For you information and to answer Michele’s question, I’m using fcc lattice for the platinum plate. Each cell contains 4 atoms and each cell has side length 3.92 angstrom. So the surface is smooth.
The way I extract position is by first gmx dump -f system.trr > out.dump to turn the trajectory file into human readable form (sorry I don’t know C++), and then read the file line by line (Python). A typical line in the human readable form is
x[ 1]={ 0.00000e+00, 1.96200e-01, 6.96200e-01}
in which I assume x[ 1] tells the index of the atom.
I’m beginning to think maybe there is an error in the gmx dump functionality? Or somehow the index of platinum atoms gets mixed when GROMACS writes to the system.trr file?
Use trjconv -dump to write out two coordinate files from successive time points and compare the coordinates of the supposedly frozen atoms. This way, you don’t have to deal with internal file formats and cryptic stuff, and you can see the evolution of the coordinates. It would be very surprising to me if freezing didn’t work, the code is pretty explicit about not updating coordinates.
Thanks for the advice, I can confirm that the frozen atoms are not moving, sorry for the false alarm! I think it’s the gmx dump function not keeping track of the index of atoms that caused the problem. Is there another option to convert .trr or .xtc files into human readable form? I tried gmx trajconv -dump 0 -o system.gro -s system.tpr but it gives the following error
Fatal error:
reading tpx file (system.tpr) version 127 with version 112 program
This is not a huge problem so don’t worry if there isn’t, and thanks for your help!
That error just means you’re using inconsistent versions of GROMACS (older version of trjconv than the version of grompp used to generate the .tpr file). gmx traj -ox can also print coordinates.
After looking carefully into the position data I noticed that i) the displacement of platinum atoms are all 3.92 angstrom, which is exactly the lattice parameter (side length of cubic cell) ii) the displacement only happens in x and y direction.
I had a similar issue once, but I was providing some gro/pdb file as -r. In that case the substrate’s molecules were moving because the order of atoms in the reference system and in the initial configuration (-c) were different.
Have you run some energy minimzation steps (steep) before starting?
For you information and to answer Michele’s question, I’m using fcc lattice for the platinum plate. Each cell contains 4 atoms and each cell has side length 3.92 angstrom. So the surface is smooth.
I had a similar issue once, but I was providing some gro/pdb file as -r . In that case the substrate’s molecules were moving because the order of atoms in the reference system and in the initial configuration (-c ) were different.
Yes I do believe this is a problem with the index, not that the freeze function is broken.
Have you run some energy minimzation steps (steep ) before starting?
No, the film layer and vapour layer are equilibrated separetely first by NVE and then by NVT, and then the system is “glued” together. But do you think energy minimization is needed? The problem that we are investigating is the mean passage time for film rupture (in average how long does it take for film to rupture start from flat). We were afraid that an energy minimization would actually rupture the film for us (film rupture because surface energy + disjoining pressure energy is minimized). But we are no experts in MD simulations so any opinion is welcomed!
Thanks for the information! Yes slip length is tricky… We used frozen solid, LJ liquid and Poiseuille flow with low pressure gradient (which is the driving force of Poiseuille flow, so low shear) and found that the slip length is around 0.5 nm, which is minimal (similar to this paper https://pubs.acs.org/doi/abs/10.1021/acs.langmuir.1c00352). It seems that many things would influence the measurement (which flow to use, shear rate, particle model etc).
I think disjoining pressure is not the issue here and I don’t think it would be affected much when running energy minimization (unless you expect the film to rupture almost immediately). Also, you can run only a handful of step (like 1000) to make sure that only the ‘ugliest’ interactions are regularized (e.g. two atoms being within each others vdW radius).
And that brings me to the second observation: how could you ‘glue’ the film(s) together without making the system explode? I would imagine there could be overlaps between molecules in the two layers that make potental energy skyrocket.
Regarding the slip length, 0.5nm is indeed quite small (I am surprised that a combination of a LJ fluid and a LJ solid could yield that number).
I am afraid this doesn’t get us closer to understanding why a frozen substrate is moving, though…
Sorry I forgot to mention this: the frozen substrate isn’t moving, it was an observational mistake. Some of my solid atoms were sitting on the periodic boundaries, and as soon as I moved them a tiny bit away from the simulation cell boundaries (0.3 angstrom) they are not oscillating in my visualisation program ;)
Also, you can run only a handful of step (like 1000) to make sure that only the ‘ugliest’ interactions are regularized (e.g. two atoms being within each others vdW radius).
I’ll have a go at energy minimization. Thanks.
And that brings me to the second observation: how could you ‘glue’ the film(s) together without making the system explode? I would imagine there could be overlaps between molecules in the two layers that make potental energy skyrocket.
I don’t glue them together closely, but leave a gap of 0.17 nm. This 0.17 nm gap is measured as the first peak in the averaged density profile in separate simulations.