I have a simple membrane bilayer system, where POPC is the membrane lipid and is equilibrated with 150mM KCL, and the net charge of the whole system is zero.
I am trying to get the electrostatic potential of the system across the z-axis using the tool gmx potential.
When I use the gmx potential command without -correct flag (gmx potential -f prd.0.conv.xtc -n index_2.ndx -s prd.0.tpr -o potenial_sys_2.xvg -oc charge_sys_2.xvg -of field_sys_2.xvg -sl 100 -center), I get following plot as shown in left side of the figure.
When I use the gmx potential command with -correct flag (gmx potential -f prd.0.conv.xtc -n index_2.ndx -s prd.0.tpr -o potenial_sys_2.xvg -oc charge_sys_2.xvg -of field_sys_2.xvg -sl 100 -center -correct), I get the plot shown in right of the figure.
Since my system’s net charge is already zero, I don’t understand why I get two different results in these two cases (I think both cases should give the same results.). Can anyone please explain the reason behind this?
Also, if there is any reference towards the workings of the command and related flag, please redirect me.
OK. So the commands you showed above were not exactly correct. Are there any other differences in the commands? At least they’re more similar with -center. What would it look like with -symm? I’m not sure why -correct would make a difference, but I haven’t used gmx potential much.
From the output, my impression is that you might need to simulate (and/or equilibrate) longer to get rid of the overall tilt of the potential curve. How long is your current simulation and for how long was it equilibrated?
I have done 15 ns of equilibration followed by 100 ns of production run. The system size is around ~45K atoms, and the box length is 6.7nmX6.7nmX9.6nm.
Even if it has not equilibrated properly and charges are distributed such that I get a potential plot like in figure-2, I am not sure why I am getting a different plot for -correct flag, for a system with a zero net charge.
Thank you for all the informations !
Have you found an answer to your problem ?
I am actually also working on a single bilayer system with a 150 mM of NaCl and wanted to get the potential with gmx potential. I am confronted to the same problem as you as the potential on each side of the membrane is not the same (with the right side being lower than the left one), and I don’t understand why.
The -correct option seems to get the potential how it should look like but I also have another system where there is an ionic imbalance on each side of the bilayer and if I use this option, it puts the potential of each side of the bilayer at 0 (which is not correct).
If you have found out why it does so or have any new informations, I would be very glad to read about it as I am getting quite stuck on this issue.
Hi Maya,
Unfortunately, I still haven’t found the reason for this. I will update if I do; if you know any alternative ways(or tools) of calculating this potential, please do share.
There is this link that explains the Poisson equation used to calculate the potential : Calculating 1D electrostatic potential profile from MD simulations – Alta Fang
The author here has written a python script for the main equation but it can be easily adapted if you don’t use LAMMPS as I do.
The problem is that I still get some surprising results on my end but maybe this can help you. At least, it lets me know that it is not an issue of gmx potential but from something else which I’m still trying to resolve.
I will let you know if I find an answer to my problem, so that it can maybe resolve your too !
There is a publication (https://pubs.aip.org/aip/jcp/article/130/21/215107/626544/Calculation-of-the-electrostatic-potential-of) that describes a similar observation with the symmetric POPC membrane that you have simulated. The authors recommend post-processing the trajectory by centering all atoms in respect to the center of the membrane. They also discuss several equations necessary for calculating the one-dimensional potential, which may help you resolve or better understand your current issue.
I don’t understand how the -correct option is a valid approach. It sets the net charge of each slab to zero. That is changing the physics.
I don’t see how one can improve the potential calculation. This is straightforward physics where there is only a single solution. The -center option is useful (and doesn’t interfere with physics). Symmetrizing can be useful to increase statistics for cases where the system studied is actually expected to be symmetric.
