Dear LAMMPS users,
I have been having trouble lately trying to run simulations that will give me the heat flux between two solids and I thought I might look for some advice.
I construct 2, regular 2D grids of atoms, each with differing atom types (1 and 2), with an equal spacing in the x and y directions of 4.5 Angstroms. Atoms of the same type are harmonically bonded to their nearest neighbour in the x and y directions. Meanwhile, atoms of differing types can interact via a Lennard-Jones (12-6) potential. The two grids are separated by a small ‘vacuum gap’ of 1 nanometre (10 Angstroms in metal units). This can be seen in the picture I have attached called ‘Simulation_visual.png’, the ‘before’ section. I designate one grid to be hot (the left one) and the other to be cold. I have adjusted the cut-off of the LJ potential so that surface atoms of either grid cannot interact deeper than 1-2 layers of the other grid (so the heat transmission is mediated by surface atoms of each grid).
The objective of my simulation is quite straightforward: I want to find the heat flux between the two solids as a result of their LJ coupling. These are the two latest methods I have employed:
Method 1- Thermostatting the two solids so that the hot group remain at 300K and the cold group remain at 10K and computing the heat flux, which I try with the ‘in.ph_2D_metal_d=3.0’ script.
Method 2 - Supplying the hot group with a constant power and subtracting that same amount of power from the cold group, and waiting to see what their steady state temperatures are (and then figuring out the heat flux from the temperature difference in the steady state), which I try with the ‘in.ph_2D_metal_HEAT_d=3.0’ script.
In both cases I do not get results that I expect. In Method 1, when I print the energy tallies of the Langevin thermostats in the thermodynamical report, I find that the magnitude of the energy that the hot atoms receive from its thermostat is not equal to the magnitude of the energy lost by the cold atoms to the ‘cold’ thermostat, in general the cold atoms lose 7 times more energy than the hot atoms gain.
In method 2, when I visualise the simulation in VMD through an XYZ dump file (and I do this because the temperatures of both solids rise dramatically given enough time into the simulation run), I see that the two solids ‘squash’ themselves and almost pull themselves towards each other (see the ‘after’ section of the ‘Simulation_visual.png’ file). They don’t appear to behave like solids.
The recurring theme of my attempts is that when I visualise my simulation in VMD, the two solids do not behave like solids and instead often morph themselves through the simulation. I do not expect to see this because I would have thought the harmonic bonds between atoms in the solids would give them their ‘solidity’. I am hesitant to use the fix spring/self command since that seems quite ad-hoc, and I am unsure about whether it will affect my numerical flux results, thereby making my system dissimilar from real physical systems and not comparable (though please do tell me if I am mistaken in thinking this way). I have also tried ‘extreme’ cases, where I replace the LJ potential with harmonic bonds (i.e. a strong interaction regime), where I should expect the heat flux between the two solids to increase, but this does not happen. In fact, oddly the strong-interaction regime gives identical flux values to when I set no interaction between atom types 1 and 2 at all.
Am I perhaps performing my equilibrations incorrectly? Have I not used the right fixes that would ensure I am modelling solids (i.e. have I unintentionally written something that is actually closer to simulating a liquid or gas)? Or is my intuition wrong in that I should not, for example, expect to find equal (but opposite sign) energy tallies for the two Langevin thermostats (for hot and cold) in the case of Method 1?
If you have any pointers on how to simulate solids and how to study the heat flux between them, even if they are general and abstract, I would highly appreciate them.
Apologies for the length of this, I hope I have been clear, please let me know if I have not explained something adequately.
16-08_12.20_NORMAL_2D.pos (74.6 KB)
in.ph_2D_metal_d=3.0|attachment (3.64 KB)
in.ph_2D_metal_HEAT_d=3.0|attachment (3.79 KB)