[lammps-users] Computing heat flux in LAMMPS: "compute heat/flux" v/s "compute stress/atom"

LAMMPS LAMMPS Mailing List Mirror
1

Hello.

  1. With regard to your first question, I posted a similar question some time back.
    This would require unwrapped coordinates to be available as variable.

Dear Steve: Is it true that the unwrapped variables are available only for dump, but not as a variable.
Is there any way to access unwrapped coordinates as a variable?

Another problem is that the PE/atom does not include electrostatic contributions. So this problem still remains.

  1. Regarding your second point: I too would appreciate if thefix_heat_flux.cpp be modified to include the modifications made by German. The pe/atom could also be modified based on Germans additions. This decision however has to be taken by Steve and the others managing the code.

  2. Dear German,
    Attached with this mail is an input file in.lj.compute, that combines both the input files you had posted ie using variables and using compute heat/flux. For testing, I have significantly reduced the timesteps.
    Here are a few values from the flux.log file:

Step Temp Press TotEng PotEng flux[1] flux[2] flux[3] Jx Jy Jz Volume
0 40.797567 -100.90087 -63.724713 -68.275743 -1.2503741 0.35527064 -2.0326868 -1.2503741 0.35527064 -2.0326868 33823.124
1 40.860126 -102.07204 -63.724594 -68.282602 -1.2789893 0.32576403 -2.0379642 -1.267529 0.35131522 -2.0140116 33823.124
2 40.928747 -103.40437 -63.724594 -68.290258 -1.2958122 0.32209609 -2.0182231 -1.284414 0.34793316 -1.9940625 33823.124
3 41.003351 -104.87917 -63.724528 -68.298514 -1.3115363 0.318944 -1.9981188 -1.300198 0.34505271 -1.9737588 33823.124
4 41.083841 -106.50637 -63.724516 -68.30748 -1.3272088 0.31663561 -1.9776122 -1.3159291 0.34300364 -1.9530632 33823.124
5 41.170105 -108.25633 -63.724264 -68.316851 -1.3420114 0.31536598 -1.9562788 -1.3307894 0.34198112 -1.9315504 33823.124

As you can see, apart from the first step, the flux values differ slightly.
I agree that the thermal conductivity values using the two different sets agree closely.
However, I feel it would be nice to know what causes the difference in values, since both are trying to use the same formula and calculate exactly the same quantity. Also, the deviation in thermal conductivity values calculated by the two methods increases as the time increases (ie as the nuber of points over which averaging is done decreases) (in fig LJ_90K.jpg).
Do you have any clues as to what is causing this difference in values?

I haven’t been able to set up the run to calculate the thermal conductivity of water, including long range effects with the codes you have sent. I shall try to complete that soon and shall post the results.

Thanks and warm regards,
Mario

in.lj.compute (1.78 KB)

log.lammps (18.2 KB)

flux.log (4.86 KB)

LJ_90K.jpg

2

The unwrapped coords are not stored in that form, but are easy to compute.
I haven't followed all the heat/flux discussion. I'm waiting for some smart
person to summarize it succintly and tell me something that everyone
agrees on. If it's possible to get rid of the pairwise loop in
compute heat/flux
and compute heat/flux for any manybody potential, then I'd like to hear how.

Steve

3

Dear All,

I have good news regarding the debate on:
Computing heat flux in LAMMPS: “compute heat/flux” v/s “compute stress/atom”
I have been able to pin-point the problem:

I have run the same input script in two different ways (Input file in.gktc_comparison attached):

1> Run for five timesteps continuously.
2> Run in a loop five times, 1 timestep per iteration.

The concise results are as follows:
( FLUX refers to the heat flux value obtained using heat/flux command)
( Jq refers to the heat flux value obtained using stress/atom…)

Continuous run results:
Step FLUX[1] Jq[1]
0 -1.7048144e-17 -1.7048144e-17
1 -1.9299501e-17 -1.929769e-17
2 -1.8497524e-17 -1.849711e-17
3 -1.9084383e-17 -1.9084762e-17
4 -1.9916934e-17 -1.9917879e-17
5 -2.0606009e-17 -2.0606629e-17
The values given by the two methods match only for the first timestep

Run in loop:
Step FLUX[1] Jq[1]
0 -1.7048144e-17 -1.7048144e-17
1 -1.9299501e-17 -1.929769e-17

1 -1.929769e-17 -1.929769e-17
2 -1.8492888e-17 -1.849711e-17

2 -1.849711e-17 -1.849711e-17
3 -1.9079635e-17 -1.9084762e-17

(more detailed results appended at the end of this mail)

As we can see, the flux value printed out, obtained using compute heat/flux, changes value for the same timestep, when run in a loop. The value printed the second time matches with the value obtained by using stress/atom!

This seems to me like some bug in the value being printed out. If the value printed the second time is the correct value, this also means that one can get rid of the pairwise loop in compute heat/flux by using a suitable combination of variables instead.

Steve, please could you look into this?

Thanks,
Mario

PS:

in.gktc_comparison (2.62 KB)

4

Dear Steve,

There is some problem in the flux value being printed out by the fix heat/flux command.
I have described the problem in the mail below.
Could you please go through it?

Sorry for reposting this issue.

Thanks.
Mario

5

yes, I'll look at it, but it may not be until next week

Steve

6

Here is the problem - if this (last) comm line is added to compute_heat_flux.cpp
then the 2 methods give the same answer:

// invoke half neighbor list (will copy or build if necessary)
  neighbor->build_one(list->index);
// invoke ghost comm to insure ghost vels are up-to-date
  comm->forward_comm();

The computation uses ghost atom velocities. They are communicated
in the middle of the timestep, but are slightly out-of-date at the end
of the timestep (when compute heat/flux is called), b/c velocity Verlet
updates them at the end of the timestep. So they need to
be re-communicated.

I'm unclear on what this means:

a) is this a bug that needs to be patched in compute heat/flux?
b) the formula method (in the attached input script) does
    not use ghost atom velocity, as it is not a pairwise computation?
c) does this mean compute heat/flux is unnecessary, b/c heat/flux
    can actually be computed via formulas (as in the script) without
    requiring its own compute and a (costly) loop over the neighbor list?
    Note that if this is correct, then since the formulas are complicated,
    we could have a new compute heat/flux that computed them without
    the need to loop over the neighbor list.
d) does this mean the script method will work for all potentials in LAMMPS,
   including manybody pair, bond, etc (which compute heat/flux will not)
e) does it mean the only thing lacking is the long-range contribution,
    which is really a separate issue, since we also don't tally per-atom
    long-range virial/energy contributions - if we did this (e.g. for Ewald),
    then the heat/flux formulas would also then be correct

Steve

7

Here is the problem - if this (last) comm line is added to compute_heat_flux.cpp
then the 2 methods give the same answer:

// invoke half neighbor list (will copy or build if necessary)
neighbor->build_one(list->index);
// invoke ghost comm to insure ghost vels are up-to-date
comm->forward_comm();

The computation uses ghost atom velocities. They are communicated
in the middle of the timestep, but are slightly out-of-date at the end
of the timestep (when compute heat/flux is called), b/c velocity Verlet
updates them at the end of the timestep. So they need to
be re-communicated.

I'm unclear on what this means:

a) is this a bug that needs to be patched in compute heat/flux?
b) the formula method (in the attached input script) does
not use ghost atom velocity, as it is not a pairwise computation?

The two expressions for heat flux are equivalent

\sum_{i<j} r_{ij} F_{ij} (v_i+v_j)/2 = \sum_{i<j} r_{ij} F_{ij} v_i/2
+ \sum_{i<j} r_{ij} F_{ij} v_j/2 =>
rename i by j and j by i in the second part => \sum_{i<j} r_{ij}
F_{ij} v_i/2 + \sum_{j<i} r_{ji} F_{ji} v_i/2 =>
use r_{ji}=-r_{ij} and F_{ji}=-F_{ij} => \sum_{i<j} r_{ij} F_{ij} v_i
The last expression is used quite often , see for example Phys Rev B65, 144306
The tensor \sum_{j} r_{ij} F_{ij} corresponds to stress per atom expression.

c) does this mean compute heat/flux is unnecessary, b/c heat/flux
can actually be computed via formulas (as in the script) without
requiring its own compute and a (costly) loop over the neighbor list?
Note that if this is correct, then since the formulas are complicated,
we could have a new compute heat/flux that computed them without
the need to loop over the neighbor list.
d) does this mean the script method will work for all potentials in LAMMPS,
including manybody pair, bond, etc (which compute heat/flux will not)

It looks like.
I tested it for Tersoff, SW and MEAM potentials for Si and MEAM for SiC.
It gives reasonable results.

e) does it mean the only thing lacking is the long-range contribution,
which is really a separate issue, since we also don't tally per-atom
long-range virial/energy contributions - if we did this (e.g. for Ewald),
then the heat/flux formulas would also then be correct

I have part of the code with Coulomb per atom for energy and stress
but as I mentioned for Ewald
summation only.

German

8

Dear Steve,

Thank you very much for looking into the matter so quickly.
Answers to your questions:

a) is this a bug that needs to be patched in compute heat/flux?

Yes. Velocities need to be up to date.

b) the formula method (in the attached input script) does
not use ghost atom velocity, as it is not a pairwise computation?

(See German’s reply)

c) does this mean compute heat/flux is unnecessary, b/c heat/flux
can actually be computed via formulas (as in the script) without
requiring its own compute and a (costly) loop over the neighbor list?
Note that if this is correct, then since the formulas are complicated,
we could have a new compute heat/flux that computed them without
the need to loop over the neighbor list.

Yes.

d) does this mean the script method will work for all potentials in LAMMPS,
including manybody pair, bond, etc (which compute heat/flux will not)

Yes. As long as the stress/atom (virial) and pe/atom are computed correctly for the potentials.

e) does it mean the only thing lacking is the long-range contribution,
which is really a separate issue, since we also don’t tally per-atom
long-range virial/energy contributions - if we did this (e.g. for Ewald),
then the heat/flux formulas would also then be correct

Yes. In fact, as German has mentioned, he has already made additions
to the appropriate LAMMPS files, so that these components can be calculated for the Ewald case.
If this is incorporated in LAMMPS, that would be great!

Thanks,
Mario

Disclaimer: I do not claim to be an expert, and would be glad to have
opinions/disagreements/criticism etc. from all on these issues.

9

ok - I want to wait for Reese to weigh in, but it
looks like I will remove the current compute heat/flux and
replace it with a simpler version that computes your scripted
formulas. All pair/bond/angle etc potentials in LAMMPS
compute the correct per-atom virial and energy, as well
as constraints like SHAKE. The only
exception is the long-range contribution. Looks like we
can add that to Ewald due to German's contribution. But I've
never seen a PPPM formulation of this.

Steve

10

ok - I want to wait for Reese to weigh in, but it
looks like I will remove the current compute heat/flux and
replace it with a simpler version that computes your scripted
formulas.

I was thinning that script was mine :slight_smile:

All pair/bond/angle etc potentials in LAMMPS
compute the correct per-atom virial and energy, as well
as constraints like SHAKE. The only
exception is the long-range contribution. Looks like we
can add that to Ewald due to German's contribution. But I've
never seen a PPPM formulation of this.

I don't know how to do it for parallel computations. For Ewald it was simple.

German

11

Well, it took me a while to get back to it, but I just released a patch
with the new formulation of the heat flux computation that German
worked out and others clarified for me. This one is actually simpler
to calculate since it leverages other computes, and more robust,
since it includes all ff contributions, except long-range Coulombics.
The previous one was just pairwise contributions.

German, Mario, Reese, please check the formulas in the new doc page.
Possibly there is an extra factor of V, in the J and kappa formulas?
Also is the I<J summation correct in both versions of the J formula?

Thanks to all who participated in this discussion,
Steve

12

Dear Steve,

Everything looks fine to me.
Thanks a lot!

Best,
Mario

13

Thank you Steve. Everything looks fine.
One more thing, I have modification of long range Coulomb interaction energy per atom and stress per atom
for Ewald summation only. If it’s interesting I can send it.
German.