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

Dear Steve, German, Vikas and other LAMMPS users,

I have a few doubts regarding the contents of some previous posts on heat flux calculation,
as well as with the heat flux calculations in LAMMPS:

--------------------------------------[Questions for German:]-----------------------------------------------------

Dear German,

In your script, as a step to calculate flux, you have used:
variable VX atom -(c_SA[][1]*vx[]+c_SA[][4]*vy[]+c_SA[][5]*vz[])/1.6021765e6

Could you please tell me why you have divided by 1.602…?
This division is only for the Jv component and not the Jc component.

Could you please share the code you have to include electrostatic effects in heat flux calculations?

--------------------------------------[Question for Vikas/Everyone]-----------------------------------------------------

Dear Vikas,

I am not too familiar with C++. I was looking at the compute_heat_flux.cpp file:


148 xtmp = x[i][0];
149 ytmp = x[i][1];
150 ztmp = x[i][2];

167 delx = xtmp - x[j][0];
168 dely = ytmp - x[j][1];
169 delz = ztmp - x[j][2];

174 eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair);

181 double vx = 0.5*(v[i][0]+v[j][0]);
182 double vy = 0.5*(v[i][1]+v[j][1]);
183 double vz = 0.5*(v[i][2]+v[j][2]);
184 fdotv = factor * fpair * (delxvx + delyvy + delz*vz)

186 Jv[0] += fdotvdelx;
187 Jv[1] += fdotv
188 Jv[2] += fdotv*delz;

Is it that, for the value of the force between a pair of atoms,
the magnitude of the total force has been used, instead of the vector components?

The reason I am asking this is:
I was trying to compare the values of the heat flux obtained by two methods:
viz using compute heat/flux v/s compute stress/atom.

I am expecting the flux values to be exactly the same,
but the values differ.
The kinetic contributions to the flux are identical, but the potential contributions differ.
I am trying to understand why the values are different.

The input script used for this comparison (in.gktc_comparison) is attached.

--------------------------------------[Question for Everyone]-----------------------------------------------------

Has anyone tried to calculate the thermal conductivity using the Einstein formula?
(Similar to MSD being used to calculate diffusivity)

i.e. Lamda_x = 0.5 * 1/(VkBT^2) d/dt < [e_x(t) - e_x(0)]^2 >
-for long t.

This would require unwrapped coordinates to be available as a variable.
I can currently find unwrapped variables only in the dump command.
Is it possible to access these as variables?
To handle charged systems, this would also require pe/atom to include long range contributions.

Since you have obtained long range contributions for the flux calculations,
would it also be possible to obtain the contributions to PE?]

Awaiting your reply and thanks in advance.

Warm regards,

in.gktc_comparison (2.36 KB)

Dear Mario and Majid,
Attached please find tar file with modification of compute_heat_flux
including long-rage Coulomb interactions.
I hope I did not forget anything. I would really appreciate any
comments. The code works with pair potentials and
Ewald summation only.

Answers to the questions:

1) the coefficient 1.6021765e6 is coming from text of stress per atom
In the code trace of this value divided by 3*volume gives pressure in
Bars, see manual,
for units metal the coefficient is 1.6021765e6. Thus I divide stress
by 1.6... to return to units metal.
All values used to calculate J_c are in units metal so I did not
change anything in this place.

2) See attachment.

3) It uses scalar multiplication of force vector by velocity vector.
The difference in potential energy contribution looks strange to me.
According to you script you should have
difference in J_v part because you did not return stress values to
units you use.
I've compared Ar with LJ and both method gives the same results.

4) I did not. The attached ewald.cpp contains contribution to PE.

All the best,

coulomb.tar.gz (8 KB)

OK..some comments and questions for improvement.

1- COMPUTE H/F: for simplification (I hope), why not go back to the
definition of heat current and use (d/dt)*sum_i(r_i*E_i)? Ei is total
No virials and works with any multibody potential.

2- PER ATOM ENERGY COMPUTE: for generalizing, why not move your modification
in compute_h/f.cpp to compute per atom (ewald) long range interaction energy
to compute_pe_atom...so that it may be called whenever a per atom energy is