calculation of the surface tension

Dear all,

In order to check if I understand what LAMMPS does, I have calculated the surface tension of my system (periodic liquid with two surfaces along z-direction) in two different ways.

First, I have done it directly in the input script using the stess tensor as:

fix 2 all temp/rescale 100 235.0 235.0 10.0 1.0
compute my_pres all pressure 2_temp
variable my_gamma equal lz*(c_my_pres[3]-(c_my_pres[1]+c_my_pres[2])/2.0)/3.0/2.0 # [bars * Angstroems] = [10^{-5} N/m]

This has yielded at the MD-step #10 the value 2068.2295x10^{-5} N/m.

Second, I have output the particle coordinates, velocities and forces using the dump command,

and calculated the surface tension as

gamma = SUM_over_all_particles ( m*(vzvz - (vxvx + vyvy)/2.) + (zfz - (xfx + yfy)/2.) ) / (32lx*ly)

and the result was -0.099977427, which is different from the first version and looks strange because it is negative.

An elementary c++ program that read the dump file and computes the gamma is attached below together with the compiled executable and the trial dump file. The executable is called surface_tension and needs two command line argument as an input (the first argument is the name of the dump file and the second argument the time step for which the surface tension should be computed).

I have check very carefully many time all the units conversion and my elementary c++ code and have not found an error.

This is why I assume I am doing something wrong with LAMMPS.

I would appreciate very much any clue regarding why this two analogical way of computation of the surface tension yield different results.

The version of LAMMPS that I use is dated by April 2, 2014.

Thank you!
Anton

surface_tension.cpp (1.77 KB)

surface_tension (15.3 KB)

dumpfile (164 KB)

in.nodd (903 Bytes)

data.Hg.dp.N1000.lxlz.d13.69.relaxed (49.3 KB)

You are going to have to improve substantially your pitch to get any developer to help you with this request. Lammps does what you are asking it to do. Your claim about the two approaches having to render the same result is nothing but just that: “your claim”. This post seems to look for answers to a personal inquire that does not seem to point out to a bug in the code or a miss-implementation of a certain method.
I guess you have placed yourself at the mercy of the elements…
Carlos

Dear all,

In order to check if I understand what LAMMPS does, I have calculated the surface tension of my system (periodic liquid with two surfaces along z-direction) in two different ways.
First, I have done it directly in the input script using the stess tensor as:

fix 2 all temp/rescale 100 235.0 235.0 10.0 1.0

I most certainly would not look seriously at any input for a production run that is using temp/rescale. That is just a bad idea and will mess up your dynamics in unpredictable ways.

Axel.

Dear Axel, I am currently not interested in dynamics. I just want to understand how does lammps work. That is why I took just a couple of lammps steps to see what it does.

Dear Carlos, thank you for your critics. But if you would carefully read my email you would notice that I am NOT “claiming” anything and that there was any offensive from my onto lammps. instead I said that most probably I am doing something wrong with lammps. the lammps documentation sometimes does not give fully description of some lammps features. I have provided all the scripts, which are actually working and well-commented and quite small. I would also appreciate some practical remarks from your side on how to improve more my pitch so that it would be more understandable? Thank you!

Dear Axel, the temperature compute for the aforementioned problem plays no role. if to calculate the stress tensor without the contribution from the kinetic energy and compare it with the virial contribution computed from the dump the result still unfornunately stays different (in the case if to use my lammps scripts where i have confused something most probably). The kinetic energy contribution to the surface tension should vanish in the equilibrium, and equilibrium properties do not depend on the type of thermostat.

Dear Axel, the temperature compute for the aforementioned problem plays no role. if to calculate the stress tensor without the contribution from the kinetic energy and compare it with the virial contribution computed from the dump the result still unfornunately stays different (in the case if to use my lammps scripts where i have confused something most probably). The kinetic energy contribution to the surface tension should vanish in the equilibrium, and equilibrium properties do not depend on the type of thermostat.

I don’t care. That doesn’t make it any less stupid. Why do a thing that is bad, when it cost you nothing to do it right from the beginning?

Of course you have the freedom to choose differently, but so do I have the freedom to ignore who does things that I consider unreasonable.

Dear Axel, thank you for the reasonable discussion. I will correct the script and check what it will produce.

Dear Carlos, thank you for your critics.

Wan't criticizing but advising.

But if you would carefully read my email you would notice that I am NOT
"claiming" anything and that there was any offensive from my onto lammps.

Not sure why you felt I was deeming your email upsetting or offensive. All
I said was that you provide two formulas that don't look identical (one
with pressure computes and the other one with velocities and forces) and
two codes (lammps and yours) and then you tell the world the outcomes
should be the same. I call that a personal claim. Developers don't know in
many cases what surface tension is or how should be calculated thus, your
email proves confusing by adding extra lingo that does not help with your
final goal of getting the issue resolved.

instead I said that most probably I am doing something wrong with lammps.
the lammps documentation sometimes does not give fully description of some
lammps features. I have provided all the scripts, which are actually
working and well-commented and quite small. I would also appreciate some
practical remarks from your side on how to improve more my pitch so that it
would be more understandable? Thank you!

Again, my point here is that you did not phrase the Q properly, in my
opinion. If your two formulas are the same then you should explicitly state
something like:" please note that the second one represents the same
calculation as done within lammps with the pressure compute but instead
using velocities and forces in an external script". If that is not the case
then realize you are posting a physics Q not a lammps Q. Do you have any
idea how long people spend reading at email requests? Extra words, lingo
and unclear statements can result in the developers overseeing some real
issue because the person posting the request could not explain his/her
ideas the right way.

Finally, to improve your pitch my advice is:
1-) realize what this list is about and phrase your request accordingly
2-) be brief and to the point to maximize attention
3-) do not assume the developers know all topics in chemistry/physics
4-) do your homework by browsing the archives first for similar Qs to yours
5-) try to avoid at all cost sounding like you are "begging with a shotgun"

Carlos
PS: I am having a really good day so I hope my message doesn't come to you
as cranky but as "mature" :wink:

you don’t say what units you are running your simulation in,
but for real or metal units, LAMMPS has internal unit conversion
constants defined. In your 2nd expression for gamma, m*v^2 will
not yield energy in real or metal units, unless one of those
conversion factors is used.

Steve

In addition to the above comments, I believe you also need to check the formula you’re asking LAMMPS to use—I’m not sure where the extra factor of 3 in the divisor is coming from. It should just be (lz/2) * [pzz - (pxx + pyy)/2].

—AEI