Dear all,
I am running a simulation to see streaming velocity profiles on Drude polarizable molecules with rigid SHAKE C-H bonds in a nano-channel. Things seem to be working (amazingly), including a temp/profile compute which I have set to thermostat my system at 300K via an NVT Nose-Hoover integrator for the Drude cores (with a separate NVT Nose-Hoover integrator for the Drude particles), but outputting the values of the temp/profile thermostat gives temperatures around 260K. (I have attached the relevant parts of my input script below.)
What is odd is that my system’s dynamics appear to be correctly 300K. A temp/drude compute during the simulation shows the cores’ temperature at 300K, as well as a temp/com compute rerun after the simulation on velocity dump files. I think this is related to SHAKE: when I rerun on the velocity dump files with fix shake my temperature computes are around 300K, but when I rerun on the velocity dump files without fix shake all the computes show ~ 260K temperatures, which is the same range as my in-simulation temp/profile values and also theoretically sensible (if the temp computes think there are more DOFs than there actually were during the simulation, it will calculate a lower temperature).
Is this correct / incorrect behavior, and what should I do about it (if anything)? My system shows similar dynamics under a temp/com thermostat, so I don’t think the temp/profile thermostat is doing anything wrong internally.
Cheers,
Shern
[Input script snippets:
group sol type 2:29
group MOVATOMS type 2:16
group DRUDES type 17:29
group SHAKECH type 4 5 6 7 15 16
fix SHAKE SHAKECH shake 0.0001 20 2000 t 15 16
compute TATOMS MOVATOMS temp/profile 1 0 0 z 5 out bin # various values of xflag, yflag, zflag, and nbin don’t change the behavior
fix DIRECT all drude/transform/direct
fix NVT1 MOVATOMS nvt temp 300. 300. 100
fix_modify NVT1 temp TATOMS
fix NVT2 DRUDES nvt temp 1. 1. 20
fix OUTATOMS sol ave/time 1000 1 1000 c_TATOMS[*] mode vector ave one file out_tstat.txt
fix INVERSE all drude/transform/inverse
]