# Heat Flux Calculation using Fix Langevin Output and Fix Heat

Hello,

I am trying to use both the Langevin fix and heat fix to calculate thermal conductivity. My question is on how to calculate the heat fluxes. Using the output generated from the fix langevin command. In the documentation example it calculates dQ as follows

dQ = 8000atoms * 0.5*(0.905+0.947) / 100/ 18.82^2 / 2

dQ = (Atoms in system)*(Average of tallied_energies)/(time_of_simulation)/(Area)/2

why are the fix outputs 0.905 and 0.947 the same sign? If one is a source and the other a sink shouldn’t the signs be opposite? Why multiply by the total number of atoms in the system?

In the simulation I am running to calculate the thermal conductivity my model appears as follows.

(Fixed End)(Source)-------------------->>>heat flow>>>---------------------(Sink)(Fixed End)

I apply the source and sink until a steady state is reached.

The total number of atoms in the system are 51200.

The total number of atoms in the source are 1850.

The total number of atoms in the sink are 1850

The units are metal

The timestep is 0.0005 picoseconds

After a steady state has been reached the time step is reset and at time step 10,000,000 of the reset run the outputs from the Langevin fixes are

Langevin Source fix output: -5564.9849 eV

Langevin Sink fix output: 5560.3983 eV

(Since the fixes tally energy the units are eV?)

The cross sectional area is 768.8*10^-20 m^2.

There is no need to divide by 2 since the heat flows only in one direction.

If I follow the example then…

dQ = 51200 atoms * 0.5(5564.9849 eV+5560.3983eV)/(0.0005ps*10,000,000steps)/Area

(I took both fix outputs here to be the same sign. As a side note averaging the fixes with one positive and one negative seems to produce arbitrary results since the values are nearly identical and the value produced is a result of small fluctuations from the outputs.)

dQ = 56961.961984 eV/ps/Area

This did not seem right so I modified the dQ equation to

dQ = (Average of tallied_energies)/(time_of_simulation)/(Area)

this produces

dQ = ((5564.9849 eV+5560.3983eV)/2) /(0.0005ps*10,000,000steps)/Area

dQ = 1.11eV/ps/Area

When I use the fix heat command

Fix Source_Fix Source_Group Heat 1 0.5

Fix Sink_Fix Sink_Group Heat 1 -0.5

It produces a very similar temperature profile in the same temperature range, but the heat fluxes do not match.

The documentation sample calculation for the fix heat is as follows.

dQ = (time in tau)(energy delta specified in fix heat)/(number_of_time_stepstime_step)/Area/2

The heat flux is using the previous time values and not dividing by 2…

dQ = ((1ps/0.0005ps time_in_tau)(0.5 eV/ps)/(10,000,000steps 0.0005ps run_time_in_tau)/Area

dQ = (2000 time_in_tau)*0.5 eV/ps/( 5000 run_time_tau)/Area

dQ = 0.2 eV/ps/Area

(For a steady state simulation should run time in tau and time in tau be considered for the fix heat calculation since run_time_in_tau is somewhat arbitrary? dQ can be made any number by making the simulation long or short, but in a steady state simulation the system is not changing in time. Does the following better reflect dQ for a steady state calculation?

dQ = (energy delta specified in fix heat)/Area

in which case

dQ = 0.5eV/ps/Area

This is nearly half the value calculated with the modified dQ for the Langevin calculation.

)

Can someone explain the way to calculate the dQ for both the fix Langevin and fix heat methods?

Thanks for any help.

From,

Marcus

I am trying to use both the Langevin fix and heat fix to calculate thermal conductivity. My question is on >how to calculate the heat fluxes. Using the output generated from the fix langevin command. In the >documentation example it calculates dQ as follows

dQ = 8000atoms * 0.5*(0.905+0.947) / 100/ 18.82^2 / 2

dQ = (Atoms in system)*(Average of tallied_energies)/(time_of_simulation)/(Area)/2

why are the fix outputs 0.905 and 0.947 the same sign? If one is a source and the other a sink shouldn’t the >signs be opposite? Why multiply by the total number of atoms in the system?

The outputs in the log file for the two thermostatted regions

are opposite signs. The flux between them is just the

ave of the 2 magnitudes which is what 0.5*() is doing.

It’s multiplying by 8000 b/c the scripts are in LJ units

and output values divided by N. The formula wants

the total flux, so it multiplies by N.