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