GROMACS version:
GROMACS modification: Yes/No
Hi, I have a system where water molecules (1008 water molecules) are confined between two walls, graphene sheets. Also, I am doing my simulations in NVT. I want to fix one wall and move the other wall along the Z-axis with a constant force applied onto the atoms of the sheet. For this purpose, I searched a lot regarding how I should setup my simulation. I noticed two approaches, one was pressure scaling and the other one was com pulli
eql.mdp (1.5 KB)
ng.
There are a couple of things I am working on to understand, like how can I define the maximum distance the wall can move, how to convert the units of pull-coord1-k to bar or MPa.
however, after applying what I understood from reading others’ problems and answers, I got my simulation and noticed my sheet did not even move. Here the mdp file is attached.
I would be appreciative if you could give me some hints to solve this setup.
Hi,
If you desire to ‘push down’ the piston wall with a constant force in order to effectively pressurize the confined water, you can simply use accelerate and prescribe the desired acceleration per atom.
Hi Micheal,
Just saw the reply. thanks a lot.
Micheal thanks for the reply.
So if it is the case, then I assume piston will stop moving when there is net pressure, applied pressure - pressure of the fluid = 0.
am I right ?
Basically yes. The vertical position will still fluctuate a bit, but there shouldn’t be any ‘net’ motion when the confined liquid has been pressurized.
Thanks.
This is what I figured out first.
The thing is that I was reading here and some posts misdirected me because they were saying acceleration is not the way, instead we should employ pulling. But when I checked and worked with LAMMPS, I noticed they are using fix force command which is somewhat equivalent to acceleration in GROMACS.
However, I have still one problem. I know for sure that when a group of atoms are frozen, the force corresponding to freezing affects the fluctuations contributing to pressure calculation.
So if it is the case, then surface tension also gets the impact as well, thereby hydrophilicity may not be as accurate as is expected.
I have read GROMACS manual and other forums about this matter. One solution that @Masrul mentioned is to use energy group exclusion implemented in GROMACS which is not supported in GROMACS 2020 and am not sure if it is supported in the latest version of GROMACS.
So, if I want to exclude should I just exclude the non-bonded interactions between frozen atoms? Or between frozen atoms and solvent/solute? The last one seems a bit non-sense to me because I suspect if I exclude the interactions between frozen atoms and solvent/solute, then the effect of the nanosheet I have there will disappear.
I would be appreciative if you could give me your opinions regarding this.
Hi, is it necessary to freeze the graphene sheets? Could you instead restraint the positions of atoms in all directions but the one along which you are pulling/forcing the piston?
Hi, thanks for your suggestion. I can apply restraint position on the sheet, but I have to use a high force. Anyways, I am gonna give it a shot.
Thanks again.
Hi Michele,
I have a specific question regarding acceleration.
I have read through the posts here about the acceleration in old gromacs versions.
Here what I am facing now is that which version of gromacs really gives the accurate results using accelerate mdp option?
I am asking this because I simulated the same thing using GROMACS 2016.4 and GROMACS 2023.5
Guess what, there is a noticeable change in filtered water number, it seems like ~350 water molecules pass the membrane using GROMACS 2023.5, while the same simulation using GROMACS 2016.4 results in ~90 !
I am just wondering which one should I trust? besides I know that velocity component in gromacs 2016.4 that is availble in gmx energy pkg is no longer seenable in newer versions of gromacs such as 2023.5.
Hi,
As far as I remember, the only difference between the new and old acceleration code is a I/O thing (which may be related to “velocity component in gromacs 2016.4 that is availble in gmx energy pkg is no longer seenable in newer versions of gromacs such as 2023.5”), but the actual effect of acceleration groups on the simulation should be the same.
You are comparing results from two Gromacs versions that are 7 years apart from each other, are you sure the inconsistency is not caused by something else (force field, barostat/thermostat, default .mdp options)?
Hi,
Thank you so much. Every simulation parameter is the same. Nothing is different.
I am gonna run using a gromacs 2020.7, I will share my results.
But so far, what I am seeing is big inconsistency. I mean there should be inconsistency because I am running nonequilibrium MD.
Besides, I am applying an external force using acceleration. The way I calculated force was F=ma and m was Nc * Mw. Subsequently, P is ma/(LyLz). Note there is a unit conversion as well which I did not mention but overall it is a simple calculation.
For better clarification, here is a brief on system geometry:
Ly = 3.1613 nm
Lz = 5.05434 nm
Nc = 624
Mw = 12.0107 g/mol
I get the pressure for an acceleration equaling to 1 nm/ps^2, and then I use it to get to my desired pressure.
In my simulation setup, I used these parameters as below:
freezegrps = P2 P1
freezedim = N Y Y N Y Y
acc-grps = P1 P2
accelerate = 0.06419455 0 0 -0.00013009 0 0
0.06419455 and 0.00013009 are equivalent to 50 MPa and 0.1 MPa, respectively.
I checked my simulation log file, and noticed both of them have one thing in common:
acc: 0.0617287 0 0 -0.00259593 0 0 -0.00246584 0 0
nfreeze:N Y Y N Y Y N N N
First thing that is questionable here is that the first component of the first two acc sets is not the same as what I put in the mdp file. Could you please help me with this to understand better?
The next mysterious thing is that what is this -0.00246584 0 0? Is it for the other groups which are not under acceleration and freezing? If yes, why it is negative? I mean it should not be negative because the right hand of the system is under 50 MPa, while the other side is under 0.1 MPa.
Furthermore, in newer versions I get this warning:
WARNING 1 [file mdp/eql.mdp]:
There are 1248 atoms that are frozen along less then 3 dimensions and
part of COMM removal group(s), due to limitations in the code these still
contribute to the mass of the COM along frozen dimensions and therefore
the COMM correction will be too small.
To solve this warning, position restraint should be applied.
which I think it could be the potential cause for that big inconsistency ( I could be wrong ).
Update on my results from gromacs 2020.7.
It gives the same results as 2023.5 gave me.
besides, putting restraint on graphene pistons does not change anything related to the number of passed water molecules. So, it is safe to assume whether they are frozen or put restraint on, they contribute the same in what I want.