[lammps-users] compute heat/flux

I'm coming back to our discussion about use energy/virial per atom
instead of compute heat flux,

This is a par of script for calculation of heat flux:

compute PE all pe/atom
compute KE all ke/atom

variable CX atom (vx[]*(c_KE[]+c_PE[]))
compute Jcx all reduce sum v_CX

variable CY atom (vy[]*(c_KE[]+c_PE[]))
compute Jcy all reduce sum v_CY

variable CZ atom (vz[]*(c_KE[]+c_PE[]))
compute Jcz all reduce sum v_CZ

compute SA all stress/atom virial

# J_x = S_xx*V_x + S_xy*V_y + S_xz*V_z
variable VX atom
-(c_SA[][1]*vx[]+c_SA[][4]*vy[]+c_SA[][5]*vz[])/1.6021765e6
compute Jvx all reduce sum v_VX
# J_y = S_yx*V_x + S_yy*V_y + S_yz*V_z, yx=xy
variable VY atom
-(c_SA[][4]*vx[]+c_SA[][2]*vy[]+c_SA[][6]*vz[])/1.6021765e6
compute Jvy all reduce sum v_VY
# J_z = S_zx*V_x + S_zy*V_y + S_zz*V_z, yx=xy
variable VZ atom
-(c_SA[][5]*vx[]+c_SA[][6]*vy[]+c_SA[][3]*vz[])/1.6021765e6
compute Jvz all reduce sum v_VZ

variable Jx equal (c_Jcx+c_Jvx)
variable Jy equal (c_Jcy+c_Jvy)
variable Jz equal (c_Jcz+c_Jvz)

thermo_style custom step temp press etotal pe v_Jx v_Jy v_Jz vol

See comparison of two results, one obtained with compute heat_fix and
another with this script on fig LJ_90K.jpg
on linux cluster using 8 processors.
It looks like it should work with any type of potentials.
Reference Schelling PR B65, 144306 - it's not important how you treat
three-body term.

However, and it's confusing me, both methods demonstrate dependence of
the number of processors (see Fig LJ.jpg)

LJ_90K.jpg

LJ.jpg

Dear German,
Thanks for keeping this alive. I have a quick question. What is the cutoff you are using? For the case of 8 processors that the each processor has 3X3X3 box size. I am thinking that some interactions might be getting calculating twice (perhaps) due to the size constraints. Is it possible for you to run a larger box may be 12X12X12 for 1 and 8 processors and observe the differences.

Regards,
Vikas

So you've added this diagnostic to a run to compute the flux
yourself. Are you getting the same run on any number
of procs (independent of the diagnostic)? If so, and the
diagnostic just depends on quantities like per-atom eng
and virial (which also should not change as a function
of the number of procs), then I don't see how you would
get different answers as P changes.

Steve

Dear Vikas,

It was mistype. The result are presented for 6x6x6. I'll try 12x12x12

Best,
German

So you've added this diagnostic to a run to compute the flux
yourself. Are you getting the same run on any number
of procs (independent of the diagnostic)? If so, and the
diagnostic just depends on quantities like per-atom eng
and virial (which also should not change as a function
of the number of procs),

Yes they are different for different number of processors.
I checked it on linux cluster - compiled with intel_mkl and openmpi,
and also on cray xt4 computer - compiled with "storm".
It has this dependence of number of processors.

German

Dear Vikas,
I've followed Stive's recommendation to add loop geom in velocity command and
started Equilibration in nve for 100000 iterations at the same volume.
Previously
I did Equilibration and thermalisation in NPT for 100000 iterations
before Equilibration.
Now results for different number of processors (fig LJ.jpg) and
different type of
heat flux calculations (fig LJ_90K.jpg) looks reasonable.

Best,
German

LJ.jpg

LJ_90K.jpg

Dear Vikas,
Do you know if computer heat flux will works with potential hybrid /overlay
For example

pair_style hybrid/overlay buck/coul/long 10.0 morse 10.0
pair_coeff 1 1 buck/coul/long 294.6197 0.327022 0.0
pair_coeff 1 2 buck/coul/long 693.6010218 0.327022 0.0
pair_coeff 1 2 morse 0.57745 1.65 2.369
pair_coeff 2 2 buck/coul/long 1632.892739 0.327022 3.948169

Thank you in advance,
German

I think it should - although you will not get the long range
component with or without hybrid.

Steve

Dear German,
I think it should work. You can check it this way. Do couple of test simulations.

  1. With no overlay and single pair 1 2 interaction
  2. With overlay and both pairs of 1 2 interaction but make the appropriate coefficients equal to zero in the second pair (which you did not use in case 1).
  3. With overlay and all coefficients non-zero

Case 1 and 2 should give you same answers while case 3 will be different than both. You can try to do this with both pairs seperately just to be sure. Let us know what you find if you do this.

Best Regards,
Vikas