[lammps-users] compute heat/flux

Dear all,

Does the "compute heat/flux" support charged pairwise potentials such as
Coulombics interaction?


Ali Rajabpour

The compute heat/flux command authored by Reese, Philip, Vikas
was released today in a 6Jul09 patch.


Hi Reese,
I'm also interested in code for heat flux.
I would really appreciate if you can send it to me.
German Samolyuk

Hi Bill,
I have working, validated code for computing the heat flux and a

python script for computing the necessary autocorrelation and its
integration. This was developed by Vikas Varshney, Phillip Howell and
myself. I am cleaning it up and giving it to Steve P. in the next week
or so.

I can send you the cleaned up version and the correlation script at

the same time I send it to Steve if you would like.


WJ Evans wrote:

Hi Steve,

There are several advantages of the Green-Kubo (GK) method compared

to the

non-equilibrium method (NEMD). 1) Because GK simulations are done

with the

system in equilibrium, there is no imposed driving force and so the


is always in the linear-response regime. 2) The GK method can compute


entire heat flux tensor with one simulation. NEMD would require a


run for each component of the heat flux vector. There the GK method


itself much easier to the study of anisotropic effects. 3) A key


to the GK method is that it is much less sensitive to finite model size
effects because a heat source and sink are not required. Hence the GK


can use much smaller models than the NEMD methods. For example, good


for diamond were obtained with only a 4000 atom model with a

dimension of

about 2.8 nm. The NEMD method would require a model comparable tot he


free path length (174 nm).

Of course there are disadvantages of the GK method compared to the NEMD
method but I have only discussed the positive attributes of the GK


A way to implement the GK in lammps would be beneficial to many


I would be interested in obtaining a copy of any GK source code

additions to

lammps if anyone has been successful in implementing it.

Bill Evans

From: "Steve Plimpton" <[email protected]>
To: <[email protected]>
Sent: Tuesday, April 14, 2009 9:24 AM
Subject: [lammps-users] Fwd: Thermal Conductivity Calculation

usingGreen-Kubo Relations

From: Steve Plimpton <[email protected]>
Date: Tue, Apr 14, 2009 at 7:23 AM
Subject: Re: [lammps-users] Thermal Conductivity Calculation using

Green-Kubo Relations

To: German Samolyuk <[email protected]...>
Cc: Mario Pinto <[email protected]...>, [email protected]...

I don't think anyone has answered my question of how
is this different than the Muller-Plathe method. For
viscosity the MP method is almost always better than
the traditional NEMD way of computing it. I'm guessing
the same is true of MP for thermal conductivity versus
the GK method. MP will converge much more quickly
than GK b/c the induced temp profile in MP is stable, whereas the

induced heat flux in GK is a quantity with large fluctuations. In
viscosity it is the off-diagonal pressure component that fluctuates

and MP avoids computing it at all.

I also disagree with the statement below that GK can easily
be implemented for many-body potentials. It would require
inserting new code in every such potential in LAMMPS (AIREBO,

Tersoff, MEAM, ReaxFF, etc). By contrast the MP method just

works with any potential without the need to alter the potential at all.
It also just works with long-range Coulombics.

Are there other downsides to MP for thermal K that I am missing?


the coefficient of thermal conductivity for the system in thermal

quilibrium, according to Green-Kubo method, equals

k=1/(3*volume*k_B*T*T) \int^{\inf}_{0}<J(t)J(0)> dt

there <> means ensemble average.
What is good in this expression - it's calculated in thermal

equilibrium, not so good - the ensemble average should be
calculated. Usually it's done by slicing results of one long enough
calculations by overlapping in time samples separated by t_{shift}

enough for the samples to be independent.


good - thermal equilibrium, bad - long calculations to reach

convergence in <> calculations.

Usually it's good to have value of thermal conductivity calculated

by two methods.

The many-body potentials was not inserted yet, I don't think it's

something very difficult.

I think Asegun Henry did it for LAMMPS.


A conceptual question on this topic. I don't know the details of

Green-Kubo for thermal conductivity. But isn't it essentially the
inverse of the Muller-Plathe method encoded in fix


I.e. in one case you impose a temperature difference and tally the

heat flux that results. In the other case you exchange bits of heat
flux and monitor the temperature profile that results.

If the two methods are essentially equivalent, and the GK has the

various difficulties you mention, and what you've done doesn't look
like it would address many-body potentials, then why not just use
fix thermal/conductivity?


Hi all,
I was not sure if it's a good idea to all lammps users, because it

still in testing phase.

If you think it will be interesting to everybody I will do it.

I've attached slightly modified style.h file. I added compute_hf.h

to it

and compute_hf.cpp, computer_hf.h files. I've used
as a prototype, as soon as it has cycle over pair in it.

It should be compiled

make yes-dpd,

- it uses velocities of ghost atoms.

I've also attached input file for Ar computation with two
(in comment).
the results could be compared with results presented in chapter 2.18
of S. Yip (ed.), Handbook of Materials Modeling, 763–771. c 2005

Springer. Printed in the Netherlands.

I've deleted first 1000 iterations and calculated averaged values


20000 samples with 50 steps shift.

Known problems:

1) It's working with pair potentials only,
2) I do no know how to include Coulomb interaction.

I would really appreciate any comments and suggestions ...


Dear German,

I would be very grateful if you could share your code. Could you


it to me and also put it on the mailing list?

Mario Pinto

I have some variant of code with calculation of j(t) - heat

flux, actually energy flux.

For the calculation of lambda in the few components systems. I

tested it on a linux cluster using 32 Xenon processors for Ar -
solid and liquid.

It looks like it's working OK. So I can share it. I will appreciate
any information of
mistakes in the code.

But I do not know how to include Coulomb interactions in this code.
As I understand pair->single() returns only real space part of

coulomb interaction. How I can calculate rest of the sum which
is calculated in kspace?

German Samolyuk

LAMMPS doesn't store Fij anywhere, it is computed in the pair

potentials, summed to Fi and Fj and discarded. You could write
a fix that looped over the neighbor list, and called

on each pair, which will return Fij to use however you wish.

The single() method also computes Eij.


Dear All,

I would like to calculate the thermal conductivity of a simple


component LJ system using LAMMPS. The Green-Kubo relation to


thermal conductivity (lambda) is:

lambda = ( 1/(3*V*kB*T^2) ) * Integral( 0 to infinity){




j(t) = sum (over i) {v ei} + 0.5* Sum (over i,j; i not eq j) {


rij }

However, I have not been able to locate the variable, (if it


gives the force between a given pair of particles. This would be
required to
perform the above calculation.

Please can you point me as to whether such a variable exists /


other way
I can get a handle to this quantity.
If not, could you please provide access to this variable in

the next

of LAMMPS? Similarly, a variable that gives the potential between
specified particles would also be helpful.

Any suggestions in performing the above calculations are welcome.

Mario Pinto

PS: I have looked at the thermal/conductivity fix, but that works
particles of the same mass; However in the near future I would


out computations for systems consisting of a mixture of


Dear Steve,
I have the same question.
As I understood from text, it dos not include kspace Coulomb contribution.

It should work with pair styles that have a cutoff Coulomb (not long range).


Correct - the heat/flux compute is not including
any long-range Coulombic contribution. Not even
sure how to do that.


The expressions for long range Coulomb contribution can be found, for
example, in Phys Rev B60, p292, S. Motoyama ....
ewald.cpp calculates total contribution to energy and virial from
logn range part of Coulomb interaction.
It will be easy to add long range contribution to heat flux if this
values are calculated per atom.

Dear German,
I have calculated the per atom long range Coulumbic energy for my system (polymer networks) for thermal conductivity calculation. (Polymer 50, 2009,3378). However, I did it as a post-processing of data. If Steve is willing, we can put up some effort to incorporate that in LAMMPS.

Many Regards,

I'm willing to look at how per-atom eng/virial can be computed
for long-range Coulombics (both Ewald and PPPM). If it's a parallel
computation then we could add it. But, how can you get long-range
heat-flux contributions
from a per-atom eng/virial quantity when you can't do that for pairwise
interactions? If you could, we wouldn't have needed to write
compute heat/flux.