3D pressure field - units

[vdascut]

I am using your latest code (v. 4.5) based on Phys. Rev. Lett. 2009 to calculate the 3D pressure field of a TIP4P water droplet (sphere) in vacuum (large NVT box, 15x15x15 nm^3). I used the associated pressure-tools (spherical_av) to generate the Ptheta, Pphi (=PT the tangential component) & Prr (=PN the normal component) values for the radials, as I understand. I could not find any manual though, so I would like to ask you a couple of questions.

(1) What are the units for the Ptheta, Pphi (=PT) & Prr (=PN) output? The pressure values are calculated from the centre of the box (7.5, 7.5, 7.5) coinciding with the sphere centre (r=0) till around r=129 (I assume 12.9nm or 129A). Why as far as r=12.9nm? why not from r=0 (box, sphere centre) to r=75A or 7.5nm, given that half the box length is 7.5nm?

(2) I assume I can use the above output directly in the relation: -(1/R^2)*integral([r^2 (Ptheta+Pphi)/2-Prr]dr), with the integral from 0 to 129, οr use the relation –integral([(r/Rs)^2 (Ptheta+Pphi)/2-Prr]dr), where Rs is the Laplace radius, and the integral from 0 to 129, to calculate the surface tension. Is this correct?

[ohsOllila]

First a side note: there is recent development in the actual pressure code which I am not involved anymore:
http://redmine.gromacs.org/issues/1324
http://www.lacan.upc.edu/LocalStressFromMD
You are probably using either one of these versions. I recommend the second. As written in the redmine, I do not either understand the publication of the first or cannot reproduce it. The second one I was able to reproduce and it has improvement in the normal component issue. However, I get non-intuitive fluctuations with CFD when I calculate some dihedral potentials (I should report this better). For this reason I am not completely convinced on the CFD. Anyway, this repository contains these codes for spherical averaging which other people have not shared, as far as I can remember. And the issue has questions on these.

REPLIES TO QUESTIONS:
(1) The pressure unit is bar (Note that code in http://www.lacan.upc.edu/LocalStressFromMD may have different sign convention than my so you might get -bar). The calculation is reached unrealistically far for the reason which I do not really remember. I have just ignored the numbers which are obviously too far. This could be fixed in the code.

(2) I am not sure what is R in your first relation. If it is the Laplace radius, then the relation becomes the same as the second and I think that, yes, you can use it to calculate the Laplace surface tension. Laplace radius you can calculate from Eq. (4) in http://dx.doi.org/10.1016/j.bpj.2012.08.023. Regarding the integration limits, in principle you should not go to obviously too large values (beyond the box). However, the integrand should be zero inside and outside of your droplet (Ptheta=Pphi=Prr in bulk), thus the limits do not matter as long as they cover the interface.

I have never calculated interfacial tensions or these pressure components from atomistic simulations in spherical symmetry. I would be very interested to see what you get. There may be some issues related to PME interactions.

[vdascut]

Thank you for your reply. I am using the latest code from:
http://repo.or.cz/w/gromacs.git/shortlog/refs/heads/release-4-5-localpressure
So probably it's the original one and not the ones in the publications you mention.
(1) I understand. Ok.
(2) R is the real sphere radius, but yes I can use the Laplace version. I use a large cut-off with reaction field, instead of PME, so this part should be ok. I agree about the limits, that was my first thought, but I wanted to be sure that the radials (r=0-129) refer to actual r=0-12.9nm and a conversion is not needed. Indeed, I can observe the Ptheta=Pphi=Prr (approximately) in the bulk of the sphere, so I can neglect this region in the integration. Then I have a large positive peak of the [r^2 (Ptheta+Pphi)/2-Prr] around the surface of the sphere, which of course gives me a negative surface tension for some reason. That's the problem. The absolute value "seems" ok (around 70 mN/m at 278K, although it is above the 59 mN/m of TIP4P at 300K).

[ohsOllila]

The code you use is the one I referred in the first link above:
http://redmine.gromacs.org/issues/1324
I recommed rather using this:
http://www.lacan.upc.edu/LocalStressFromMD

(2)

R is the real sphere radius

The question is now that how do you define the sphere real radius? This may not be the same as the Laplace radius (although the Tolman length should be really small for water). Anyway, the equation to calculate the surface/interfacial tension depends the intended usage of the values.

I wanted to be sure that the radials (r=0-129) refer to actual r=0-12.9nm

This depends which was the voxel size you used in the pressure calculation. The conversion factor depends on this (if I recall correctly 0.1nm is the default so this conversion would be ok). You can check the reasonability by plotting the density distribution together with pressure.

gives me a negative surface tension for some reason

This is most likely due to different sign convention used in the new local pressure versions, i.e. you get values with wrong sign by using this repository with the new pressure code versions.

around 70 mN/m at 278K, although it is above the 59 mN/m of TIP4P at 300K

Could the difference come from the usage of cut-off with reaction field, instead of PME in the droplet calculation?

[vdascut]

Thank you again for your replies. I will have a look at the results from the code at:
http://www.lacan.upc.edu/LocalStressFromMD as well. I will let you know. It might be the reaction field, I have some test runs with PME, so I will rerun those as well with mdrun_LS to compare results.

[vdascut]

Hello again,
I rerun the trajectories (coordinates & velocities sampled every 10ps) with the code of http://www.lacan.upc.edu/LocalStressFromMD, while using your tensor tools for averaging in the spherical coordinates (spherical_av). Now the signs are correct (big negative peak of the PT-Prr around the surface). I get really smoother curves for the PT-Prr profiles for two TIP4P spheres with around 6900 water molecules each and with the electrostatics treated by RF & PME, but with huge rerun-time (slower by a factor of 10-15 compared to the code I had used before). The PT-Prr profile in the sphere bulk is now really closer to zero, as it is outside of the sphere. If I use the equations (3) & (4) from http://dx.doi.org/10.1016/j.bpj.2012.08.023, the Laplace radii Rs are calculated at 3.39nm (RF) & at 3.36nm (PME). I chose to integrate between 2nm and 7.5nm (half the box size, with r=0 at the center of the water droplet). This is justified by the fact that in the bulk (r<2nm) the PT-Prr approaches zero with some minor fluctuations and contributes practically nothing to integration. The equation (3) gives me a value of σs equal to 2570 bar_A = 25.7 mN/m (RF) or 2690 bar_A = 26.9 mN/m (PME). This can't be correct. Am I loosing something in the process?

[ohsOllila]

I think it should work in principle like you described. Is ~3.4nm close to the real radius of the droplet from the density profiles? Could you share the rad_pres.xvg file? I could take a look. From bulk pressure difference, physical radius and Laplace equation you can roughly estimate what should be the tension. Then you see if the problem is in the pressure calculation or in the tension calculation.

[vdascut]

The real radius is close to 3.4, yes, around 3.5nm i would say. The rad_pres.xvg (TIP4P sphere, RF) is attached with a jpg extension instead due to write limitations on the repository. Just rename extension to xvg.

[ohsOllila]

Something goes wrong now, there is like -15000 pressure inside the droplet according to the file. Could you share the localpressure.dat file, written by the pressure code and/or the run input files? You can also create a test folder into this repo and make pull request, if you find that easier in sharing files.

I also understood that you got different surface tension values (more reasonable?) with other pressure code version. Could you share the files from that as well? I am suspecting that there has been some change in the output format.

[vdascut]

Thank you again for your time. The original dat0 file is around 230MB (you can download from: https://www.dropbox.com/s/bk23to9hsj4tqa6/tip4p_droplet_LS.dat0?dl=0). I am not sure if I can upload it over here due to its size. I attach the older xvg file (from the previous version code, using again your tensor tools). The pressure profiles are also exhibiting large pressures inside the water cluster.

[ohsOllila]

I used my tools and tools by Vanegas et al. to the dat0 file you sent and got the same results, see 2D pressure profiles in test folder in this repo (2Dlpp from this repo, prof from Vanegas et al.). The full 3D tensor looked the same as well, but I did not upload due to the size.

So it seems that there is no problem in this analysis. It indicates that you actually have such weird pressures, or then the pressure code does not work for one reason or another. The local pressure code give ~-820 for the total pressure value. What do you get with g_energy?

[vdascut]

For the (RF) NVT 15x15x15 nm^3 box (6900x TIP4P waters):
Energy: Pressure
Average: 0.0246057
Err.Est.: 0.18
RMSD 10.4523
Tot-Drift: 0.473935 (bar)

                             Average         Err.Est   RMSD    Tot-Drift

Pressure 0.0246057 0.18 10.4523 0.473935 (bar)
Pres-XX 0.0587412 0.59 14.3323 0.567242 (bar)
Pres-XY -0.326447 0.51 8.74117 3.03044 (bar)
Pres-XZ -0.141184 0.45 8.42659 2.54251 (bar)
Pres-YX -0.330609 0.51 8.86218 2.95042 (bar)
Pres-YY 0.540273 0.35 14.0358 0.517445 (bar)
Pres-YZ -0.0239573 0.35 8.13787 -0.00233738 (bar)
Pres-ZX -0.069289 0.42 8.42662 2.17443 (bar)
Pres-ZY 0.0207717 0.35 8.13245 -0.27151 (bar)
Pres-ZZ -0.525197 0.5 13.9116 0.337117 (bar)

[ohsOllila]

Can you share your run input files? I try to reproduce this.

[vdascut]

I 've used a short trajectory (0-5ns total, PME) to rerun with the localpressure code (between 1-5ns) . The results are not that different from RF. You can download the 1-4ns trajectory and rerun it with the localpressure code from: https://www.dropbox.com/s/rshucfjn5l1azzb/water_drp.tgz?dl=0

[ohsOllila]

I got the same results as you.

I was thinking this a bit and playing around. I think that the reason for the results is that you have ran simulations with Gromacs 5 (right?). However, the pressure codes are written for Gromacs 4.5. Since Gromacs 4.6. the cut-off treatment has been changed. So the cut-off settings you are using in the simulation are not available in Gromacs 4.5. and older. Due to this forces used in the simulation and pressure calculation do not match and the results are wrong.

I made 10 frame trajectory with Gromacs 4.5. and tested. It seemed reasonable but with lot of fluctuations naturally. Could you test what you get if you run with Gromacs 4.5.?