Dear all,
Does temp/partial always divide by 3 when calculating the temperature or does it divide by the number of dimensions you are applying the thermostat to? The documentation suggests that it does, but this not clear.
I think I may have also found a bug. When using temp/partial for only two dimensions (doesn't matter which two), the temperature calculated by this compute is 4/3 times the temperature calculated from temp/com or temp (when there is no com motion). I don't see where this factor comes from. If temp/partial always divides by 3, I would expect this factor to be 2/3 instead of 4/3. This temperature is only off during md, not energy minimization. I set a temperature initially, and the temp/partial temperature is near temp or temp/com during minimization, but different by a factor of 4/3 during md. When using temp/partial for only 1 dimension, the compute reports a temperature of 0 during md independent of which dimension is used. The non-zero terms in the temp/partial KE tensor match those from temp/com, so the problem must be in the calculation of the temperature from the kinetic energy.
I have tried this with pure water and in a slit with and without flow and with and without a thermostat with the same result. I get the same results with the 10Nov09 and 10May10 versions.
Thanks,
Brian
The code is pretty simple - see src/compute_temp_partial.cpp.
The DOF it uses per atom is M*N where M is the number of
dims you include in the command. Can you post a simple
input script for dynamics which shows the 4/3 effect?
There should be no difference in temperature (other than statistical)
between temp/partial and temp.
Steve
Steve,
Attached is the input script and data file that I used. A little output I get with 1,2, and 3 dimensions thermostatted is also attached. The problem seems to be with SHAKE. Without it, there is no problem. I think the problem is that the number of degrees of freedom being subtracted for SHAKE is not being reduced when the number of dimensions for temp/partial is reduced. For rigid water, the number of degrees of freedom should be 3N-N = 2N for thermostatting in 3D, 2N-2N/3 = 4N/3 for thermostatting in 2D, and N-N/3 = 2N/3 for thermostatting in 1D. Instead, it looks like they are 3N-N = 2N, 2N-N=N, and N-N=0, respectively. Using the compute_modify command:
compute_modify fluidtpart extra -826
or
compute_modify fluidtpart extra -1652
for 2D or 1D thermostatting, respectively gives the correct temperatures with SHAKE turned on. I have 826 water molecules. This explains the 4/3 and 0 factors (assuming the temperature is set to 0 when the degrees of freedom are 0).
Brian
1Dlog (5.32 KB)
2Dlog (5.32 KB)
3Dlog (5.32 KB)
solvate.data (304 KB)
water.in (1.07 KB)
You're right - compute temp/partial wasn't applying the fractional weighting
for DOF to the shake constraints as well. I posted a 14May patch for this -
please see if it now works as you expect.
I believe the right formula for DOF of a body of N particles with S additional
constraints (e.g. S = 3 for rigid water due to SHAKE) is
DOF = nper/dim (dim*N - S) where dim = 2 or 3 for 2d/3d and nper = how
many flags you set in temp/partial
This is what the code now implements - see if you agree.
Thanks,
Steve
Steve,
I agree with that formula, and it works now.
Thanks,
Brian