I’m trying to calculate the thermal conductivity using the Green Kubo method for an ensemble of rigid bodies (fix rigid/nve/small), for which I get very strange results. To test this, I considered an ensemble of dimers and got very different results (1-2 orders of magnitude) when using shake vs. fix rigid/nve/small.
To further investigate this, I tried a system with one dimer with the same initial conditions. I get the same PE (identically 0) and KE (which is now conserved) and similar stress/atom elements. In this case, the stress contributions come only from the constraint fixes. However, the heat flux and thermal conductivity are very different. I’ve attached the input files, and appended the salient features to the end of this email. I’m using LAMMPS version 23Jun18, and I get the same results for the single dimer using 3Mar20.
Why is there a difference between shake and fix rigid/nve/small? Is it possible to calculate the correct heat flux and thermal conductivity with fix rigid/nve/small?
Please see the attached input files and log files for 1000 steps. I’ve only attached the main lammps log file, but other data (e.g. sum of atomic stresses) are also outputted in separate files through ave/time fixes.
I was wondering whether you have had a chance to look at this? It appears that the per-atom stresses are not identical for shake and rigid/nve/small, so I suspect the difference in heat flux originates from differences in the applied internal constraint forces.