Heat Flux

Dear Users,

I want to check equality of these two methods to obtain J_x ( heat
flux in x direction) but each of these methods gives different answers
for J_x.
are not these two methods equal ?

Method 1:

compute Stress all stress/atom virial
compute KE all ke/atom
compute PE all pe/atom
variable Etot atom c_KE+c_PE
variable jx atom
compute 1 all reduce sum v_jx
fix 1 all ave/time 1 1000 1000 c_1 file jx.txt

Method 2:

compute Stress all stress/atom virial
compute KE all ke/atom
compute PE all pe/atom
compute Flux all heat/flux KE PE Stress
variable jx equal c_Flux[1]/vol
fix 1 all ave/time 1 1000 1000 v_jx file jx.txt

Thanks a lot,

Unless I've missed something, then those expressions should be equivalent. What sort of differences are you seeing?

I also believe that these two should be equivalent. Answer for the
first method is 25 times greater than for the second method.


Of course, this order of magnitude is not constant! it sometimes rise to ~1000.


I was just wondering shouldn’t you use [1], [2], and [3] for vx, vy and vz?


I don't think so. Am I correct in thinking that in the atom style variable command, vx gives you the x component of the velocity for each atom?


Would you mind attaching your full input deck? I've just put the following lines into one of my input scripts, and the two fluxes that I see are the same:

compute s all stress/atom virial
compute ke all ke/atom
compute pe all pe/atom
variable etot atom c_ke+c_pe
variable jx_i atom (v_etot*vx)-(c_s[1]*vx)-(c_s[4]*vy)-(c_s[5]*vz)
compute jx_m all reduce sum v_jx_i
compute jx all heat/flux ke pe s

thermo 10
thermo_style custom step c_jx_m c_jx[1]

Did you check the units? For the first to work you'll need the stress
to be in the same units as the energy.

That's a good point- the simulation I ran was in LJ units, where the conversion factor is 1.

Thanks for your replies.

Dear Niall, I have solved my issue by your valuable guiding. The
conversion factor for metal units (units for my input) ,as mentioned
in update.cpp file in src directory, is not 1 unlike your input units
(LJ). it is approximatly 10e6. Now my problem, how do I calculate this
conversion factor for metal units?

Thanks again.

I assume the one you've seen in update.cpp is nktv2p. Dividing by this constant turns stress*volume units (bar angstrom^3 for metal units) into energy (eV) units.

Something like this would do it:
variable conv equal 1.602176463e6
variable Sxx atom c_Stress[1]/conv

Your quite right. (sorry for delay!)

Thank you very much.