1. In the MD Analyzer, the temperature profile that is shown, is that the average temperature at each point between tmin and tmax?
The cell is divided into bins along the selected axis, and the average temperature of the atoms in each bin is calculated by averaging over all snapshots.
Using the tmin and tmax parameters, you can average only over a selected time window of the simulation, e.g. if you increase tmin you will discard the initial section of the trajectory, where the system has not yet reached a steady-state. You can vary the tmin and tmax parameters to monitor the convergence of the temperature profile.
2. How exactly is the heat rate reported at the end of the log file, calculated? Is it also an average heat rate over all the exchange processes in the simulation? Or is it also averaged over all the simulation steps (so that essentially higher the simulation steps, lower the heat rate)?
I suppose you are talking about the average heat current? This is not a heating rate, as it is not really associated with the temperature. It is calculated as the energy that is externally transferred from the heat sink to the heat source per time. The average is carried out over the entire simulation, so it includes also parts before the steady-state has been reached, but normally, if the simulations is long enough, it should have converged to a constant value at the end of the simulation.
3. Do you have a suggestion for a "good" heat flux? Or could you suggest a method for determining it based on maybe material parameters?
It is difficult to say, because it depends on many things, such as the cross-sectional area of your system or the average temperature. In general, the lower the heat flux, the less non-linear effects are included, but the more difficult to extract an accurate gradient from the temperature profile. I usually use exchange intervals of 100-200, but you probably need to make some tests for the system you are interested in.
4. In a paper http://www.tandfonline.com/doi/pdf/10.1080/0026897031000068578, which also uses the same formalism, it is stated that atoms at the interface are loosely bound by weak harmonic springs to avoid dissolution. Can that be done in VNL?
This will be possible in VNL-2016, but it is not possible in VNL-2015.
5. Does this simulation consider atomic vibrational energy, or is it only translational energy that is considered?
What do you mean by "translational" and "vibrational" energy? In MD simulations there is no such distinction, both contributions are natively included in the equations of motion.
6. Any system, irrespective of whether it has an interface, would come to the same temperature throughout if allowed sufficient time to relax. So instead of calling temperature profile at equilibrium, I would prefer calling it temperature profile at steady state, where the heat exchange rate and the time of relaxation are consistently chosen. It seems, therefore, that the choice of this combination is very critical. Do you have any suggestion regarding this, or any rule of thumb to follow based on material used, cross-section area, etc.?
You are right that steady-state is a more suitable expression than equilibrium in this case. In NEMD you only have one adjustable parameter, which is the exchange interval, so I don't quite understand what you mean by relaxation time?
7. Do you have any suggestion about the length of the "Hot" and "Cold" regions? I believe this also has a bearing on the resultant heat flux, since a longer length would possibly give a broader distribution of velocities and therefore, higher flux--is that understanding correct? Apart from that, what else needs to be considered for these regions? Once again, could you give a good number of atoms that you think should be included in these regions?
The regions should be small compared to the central part, through which you want to measure the conductance, because there is always a non-linear temperature change around the sink and source region, which one would like to keep comparably narrow to avoid blurring the profile. As you correctly say, the lengths might have some influence on the temperature difference of the two regions, but I have not yet come across a situation where this would have a decisive effect on the end results. Some papers suggest 10% of the length of the central part, I typically use between 5 and 10 Angstrom, but in the end you need to investigate yourself.
8. Some papers have set the boundary condition as temperature, and have calculated the heat flux. Is there a way in VNL to do so, i.e. have thermostats connected to some of the atoms, and see the evolution?
You can apply separate NVTNoseHooverChain-thermostats only to the two sink- and source-regions, to keep them at well-defined temperatures (in VNL-2016 you will be able to do this with all thermostats) by saying
method = NVTNoseHooverChain(
time_step=1*femtoSecond,
reservoir_temperatures=[('heat_source', 330*Kelvin), ('heat_sink', 270*Kelvin)],
thermostat_timescale=100*femtoSecond,
initial_velocity=initial_velocity
)
You can retrieve the total transferred kinetic energy of each thermostat, by querying
heat_flux_hot = method.thermostats()[0].transferredKineticEnergy()/(number_of_steps*time_step)
heat_flux_cold = method.thermostats()[1].transferredKineticEnergy()/(number_of_steps*time_step)
The problem here is that both thermostats are not synchronized, i.e. the heat removed from the cold region may not be equal to the heat added to the hot region.
We are planning to extend the implementation such that you can specify the heat flux, but I don't know if this will make it into VNL-2016.
9. If I am interested in changing the exchange rate over time, which cannot be currently done directly in VNL, is it possible to carry over the results of one momentum exchange simulation to another? I know I can carry over the bulk configuration and velocities, but at this point I am a bit confused if that should be sufficient for the next momentum exchange simulation I would be looking for. Do you think I need to read some other results from the previous simulation and use that for some calculation in the next momentum exchange step?
If you restart the simulation from the last step of a previous MD simulation, using
initial_velocity=ConfigurationVelocities()
, then it works essentially as if you continue the same simulation, as all information is stored in the positions and the velocities. So use can use this approach to change parameters during the simulation.
