Use molecular dynamics to calculate the surface tension of argon liquid

[Jin You Lu]

Dear Sir

I am wondering if it is possible to obtain surface tension (or local pressure) directly from md_trajectory file?
The local pressure can be obtained by using Eq. (7) in the reference.
http://marge.uochb.cas.cz/~jungwirt/paper74.pdf
The first term can be obtained by local temperature times local density.
However, the second term seems to be more complicated.

The attached file is to get a stable vapor liquid interface of liquid argon at 85K by using NVT Berendsen ensemble.

[Julian Schneider]

The surface tension is commonly calculated from MD simulations as gamma = L_z/2 * <P_zz - 1/2 * (P_xx + P_zz) > You can easily get that quantity from the md_trajectory via

Code

# Get the pressures tensors of all snapshots.
# Possibly discard the first snapshots for equilibration.
p_tensors = md_trajectory.pressureTensors()

# Strip the units to be able to use numpy.mean().
p_tensors = p_tensors.inUnitsOf(eV/Ang**3)

# Get the height L_z
cell_height = md_trajectory.cells()[0, 2, 2]

# Average the pressure components over the snapshots. 
surface_tension = numpy.mean(p_tensors[:, 2, 2] - 0.5*(p_tensors[:, 0, 0] + p_tensors[:, 1, 1]))

# Apply the units again and normalize with respect to the two surfaces and the height.
surface_tension = surface_tension*eV/Ang**3*cell_height/2.0
print surface_tension

[Jin You Lu]

Dear Julian Schneider

Thanks for your reply.
May  i know how to get local pressure, such as P_zz(z) ,P_xx(z) , and P_yy(z)? 

Best
JY

[Julian Schneider]

That's  a bit more complicated. Probably the easiest way is to compute the local pressures via the LocalStress Tensor.
Note, that in ATK the Stress and LocalStress analysis objects give the internal stress, which has the opposite sign as the (external) pressure.

I have attached a script that does that. Please check if that works for your purposes.

[Jin You Lu]

Thanks for your reply.
Actually, i did several tests to extract the surface tension and local pressure with different number of argon and simulation box.
I always get a stable density aorund 1.2 (from MD_analyzer). But when i try to get local pressure and surface tension. The values are not comparable with the results in the literature.

One of examples is as follows:
Firstly, I place 2254 argon atoms at the center of a rectangular simulation box with 384x384x2500(A^3).
and I equilibrated the system at 100k by using Berendsen algorithm. This step can be found in the "Berendsen_thermostat.py" (my attached file).

Then, I wrap the atoms into the simulation box and used NVE verlet to extract surface tension and local pressure.
For P_zz(z) normal direction of local pressure, it should be constant for any z if the equilbrium between liquid and vapor is achieved.
But coomparing to the simulation result, it is not stable near and inside the argon liquid. (Figure 1)
For surface tension, it should be around 10*(10^-3)N/M.  But i got 23.3294240614 Ang*bar( 23*10^-5N/M).

Please help to check it. Thanks for your continuous help.

 

[Julian Schneider]

The numbers you compare your results to, is that experimental results or simulation results with the same potential?

[Jin You Lu]

It is simulation results with the same potential.  You can help to check these two references.

[1]  Molecular Physics 48(6):1357-1368 · April 1983
DOI: 10.1080/00268978300100971
I did similar simulation with the same potential in the ref. [1].
We can check the normal direction of pressure in a vapor/liquid equilibrium condition.
Figure 3 and 4 show normal and tangential direction of local pressure P(z), respectively.
The normal direction is almost constant for any z position. and for Pt, there is a transition between liquid and vapor.
https://www.researchgate.net/publication/239737668_The_Pressure_Tensor_at_the_Planar_Surface_of_a_Liquid

[2] International Journal of Chemical, Molecular, Nuclear, Materials and Metallurgical Engineering Vol:3, No:9, 2009
http://waset.org/publications/3746/molecular-dynamics-simulation-of-liquid-vapor-interface-on-the-solid-surface-using-the-gear-s-algorithm

The liquid and vapor density versus temperature can be found in ref. [2].
The calculated density of my simulation by using quantumwise agrees with the Fig .2.
But for surface tension, (Fig.3) the calculated value is very different. ( the surface tension value is also confirmed by appl 107, 143105 (2015), which show surface tension is 8.16 mJ/m^2 with the same order as ref. [2].

I attached figures for the reference.
Thanks for your help.

[Julian Schneider]

I have used somewhat different expressions for the local pressure and surface tensions than normally used in the literature, beause the latter require to re-calculate all interactions explicitly.
However, with the common formulas from the literature for P_N/T (z)
http://web.mst.edu/~hale/argonMC/References/JCP.1999.Trokhymchuk.LJ.liq.sim.Tc.pots..pdf
and gamma
http://marge.uochb.cas.cz/~jungwirt/paper74.pdf

I seem to get surface tension which is in the range that you expect.
The normal pressure seems to be more constant across the interface, but it is difficult to tell, because the profile is very npoisy in the liuid region.
I have attached the script I used. Feel free to double check the equations and test other formulas.

[Jin You Lu]

Dear  Julian Schneider
Thanks! I will check it and let you know my results.
Thanks again for your continuous help. i am very appreciated.

Best
JY

[Jin You Lu]

I am back!
The normal direction of local pressure can be constant after averaging 3000 snaps. (100 time step snap once, and time step i set is 5 fs)
The lateral direction of local pressure is close to the literature.
The surface tension of Argon liquid vapor interface at 90k is 0.015 (which is close to 0.0145 N/m in the previous literature)

Thanks
Best
JY