Calculate bin-wise heat flux in water confined between solids

LAMMPS LAMMPS General Discussion
1

Hi everyone,

I’m simulating heat transfer in confined water (using the TIP4P/2005 model with fix shake) sandwiched between two solid surfaces serving as hot and cold reservoirs in an NEMD setup. My aim is to compute bin-wise (slab-wise) contributions to the convective (Jk) and virial (Jv) heat fluxes along the water slab to gain insights into local heat transfer mechanisms, such as interfacial resistances.

The compute heat/flux/tally seemed promising at first, as it tallies heat flux between two user-defined groups of atoms. However, it doesn’t include kspace contributions or bonded interactions in the tally, limiting its use for my system.

As a workaround, I’m dumping positions and velocities from the trajectory and using LAMMPS rerun commands in two separate scripts:

  • Script 1 (short-range contributions): Switch off kspace, use pair_style lj/cut/tip4p/cut (approximating short-range Coulomb with undamped 1/r up to cutoff), apply fix shake, and compute per-atom kinetic energy (ke/atom) plus centroid/stress/atom (virial only) to handle SHAKE constraints properly via the centroid approach. This gives the short-range virial and convective terms.
  • Script 2 (long-range and full potential): Use the full pair_style lj/cut/tip4p/long with pppm/tip4p kspace, compute pe/atom (full potential energy, including kspace), and stress/atom (kspace virial only).

I then bin the per-atom heat currents (e * v - S . v) using chunk/atom and reduce/chunk, and combine the binned results from both scripts for total flux profiles.

This split leverages centroid/stress/atom for better consistency in rigid molecules (vs. equal virial distribution in stress/atom), but a major drawback is the approximation in Script 1: lj/cut/tip4p/long requires kspace, so using lj/cut/tip4p/cut introduces a mismatch in short-range electrostatics (undamped 1/r vs. damped erfc(r)/r in /long), potentially causing errors in the calculated fluxes. Are there better ways to achieve accurate bin-wise fluxes using LAMMPS? Has anyone tried similar splits or code modifications to extend centroid/stress/atom for kspace?

2

I did something similar several years ago. Regarding centroid/stress/atom, indeed the code includes support for kspace but you need to modify the source code to enable it (Revised implementation of a new atomic stress definition for correct computation of heat flux with many-body interactions by donatas-surblys · Pull Request #1704 · lammps/lammps · GitHub ; note that as the LAMMPS code base is evolving you need to modify at somewhere else with newer versions). However, as the original post says, this is not verified so at least you need to do the following things:

  1. verify the results. When I was working on it I compared the results with another independent code;
  2. consider how the virtual site in TIP4P affects the heat flux (the code in LAMMPS does not consider this); you can find papers discussing this problem. The magnitude of error may or may not be important for your usage.

If you just want to match the short range force, then LAMMPS supports tabulated force fields (pair_style table); though I have no idea if it works with TIP4P or how physical the results will be.

1 Like
3

Thanks for your input.

From the discussion in the post, it seems that if kspace is used in conjunction with compute centroid/stress/atom (with the relevant flags enabled), it would still compute the kspace contribution as per stress/atom. However, this kspace contribution may not be meaningful in the centroid context. To circumvent this, my idea was to use two separate reruns to isolate the individual contributions to the heat flux.

If I may ask, what kind of system were you studying?

4

However, this kspace contribution may not be meaningful in the centroid context.

My understanding is that, the original stress/atom is correct for two-body potentials, and the centroid correction is only needed for >=3 body interactions (see e.g. J. Chem. Theory Comput. 2019, 15, 5579−5587). Since electrostatic forces are two-body in nature so using the old formula for them is fine. Again, this is only my own understanding. And just like that post said, it’s unclear how meaningful the local contributions from the kspace part are (which is essentially long-range interactions); I don’t have good suggestions on this.

If I may ask, what kind of system were you studying?

It was a crystal of small molecules, and I was using the Green-Kubo method.

1 Like
5

Thanks a lot for your time. I can test and see how accurate these calculations are and come back to this post.