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.?

[vdascut]

Yes, the trajectories are from Gromacs 5.04. I will run some of them using 4.5 to check.

[vdascut]

I rerun some of my trajectories with gromacs 4.5.7 (latest 4.5 version) and then used the 4.5-LS local pressure code from http://www.lacan.upc.edu/LocalStressFromMD to calculate local pressures. I used several methods for the original runs (with 4.5.7) to treat electrostatics (plain cut-off at 1.4, 1.7, 2.0nm, RF, or PME). The results are identical to the picture before. Large pressure (-15.10^3) inside the sphere and going to zero after the sphere radius. The strange thing is that the Laplace radius is very close to the real radius of the sphere for many trajectories (including of some systems with organics, salts and waters), but the surface tension seems to be low (around half the expected value). I attach again the rad_press.xvg file of the tip4p water sphere with cut-off at 2.0nm, PME run with 4.5.7 & 4.5-LS.

[ohsOllila]

Sorry for the slow response, I have been extremely busy. I have a stupid quick question: Can you check that you have velocities included in your trajectories?

[vdascut]

Hello again, yes of course velocities are included in the trajectories, sampled at the same pace as coordinates. I had run also a trajectory without velocities (out of curiosity) and apart from the warnings I got when running the LS code (missing velocities), I got completely different results, as expected.

[ohsOllila]

This is what I thought but I wanted to check it to make sure.

Anyway, I cannot understand what is going on. I was able to get a reasonable pressure distributions by running the system with Gromacs 4.5.5. I was testing only over ten frames, though. However, if I take ten frames from your trajectory I see the difference.

I got technical problems in sharing the files, but I will do it ASAP.

[ohsOllila]

I think that I found the problem now. If *tpr file for pressure calculation is done with integrator md-vv, the pressure calculation does not work properly. I changed the line in mdp for pressure calculation to: integrator = md (https://github.com/ohsOllila/pressure-tools/blob/master/test/water_drp/prescalc.mdp) and got reasonable results for all trajecotries, see e.g.
https://github.com/ohsOllila/pressure-tools/blob/master/test/water_drp/rad_pres.xvg compared to the one calculated with md-vv: https://github.com/ohsOllila/pressure-tools/blob/master/test/water_drp/rad_presVV.xvg

[vdascut]

Thank you very much. That can be it. I am currently re-running some trajectories with the pressure LS code using the md integrator, instead of md-vv. I will let you know. Thank you again.

[vdascut]

The change of integrator (from md-vv to md) indeed produced the correct surface tension for a variety of runs. Thank you. That was it. I think there was no mention of a specific integrator in the manual of LS code.

[vdascut]

Hello again, One more question for the pressure-tools. Is there a way to estimate the error on the converged values for the surface tension out of the tools? Thank you,

[ohsOllila]

There is no error estimate protocol in the tools. In the publications I have estimated error by varying the integration limits of interfacial tension calculation since the pressure fluctuations in bulk region have given the biggest uncertainty to the final results in my case.

[vdascut]

Ok, thank you. I have calculated an error based on another method. I have split the frames in groups, over the long trajectory and I calculated the surface tension in all these groups. I managed to get a stdev out of the averaging. For the bulk region I also notice some fluctuations but they cancel out usually in the integration.

[ohsOllila]

Yes, the fluctuations cancel but some of the droplets I had were so small that this issue actually affected the results, so it became the biggest uncertainty. I think that you do not have this issue.

For this reason I did not what you described. In my case that error would had been negligible compared to real error bars. If possible, it would be very useful if you could share your scripts in this repository. I would be also very interested to know how large error bars you got with your approach.

[vdascut]

Actually, I am using the same tools as published, but I rerun the trajectories for different parts of the production trajectory, then the resulting surface tension is calculated using integration by Igor or Origin software. The averages and stdev of the surface tension is then calculated based on the different values from different time windows throughout the production trajectory. The stdev in Laplace Radii is quite low (up to +-0.04 nm), while for the Surface tension the stdev is up to +-1.2 mJ/m^2. This is for several different windows of 5ns each out of a long trajectory, so convergence is quite good.

[ohsOllila]

Thanks for the information. In the end the error bars seem to be quite similar to the ones in Table 2 in here: http://dx.doi.org/10.1016/j.bpj.2012.08.023