Water Temperature 2/3 desired value

Dear Lammps Users,

I am conducting a project using lammps involving evaporation of water inside nanopores. One of the measurements I’m trying to take is the temperature of the water molecules (tip4p). To do this I have written a Matlab code but the answers have been off and running a test simulation of only 1000 water molecules shows answers are exactly 2/3 of the actual temperature.

Some assumptions I’ve made:

  1. Temperature is set using ‘fix fx tip4p nvt temp 300.0 300.0 200.0’ command

  2. Lammps temperature is given by ‘thermo_style custom step temp’ command

  3. Matlab temperature is given by: 1/2 m v^2 = 3/2 N k T

  4. Kinetic energy is evaluated from both hydrogen and oxygen atom velocities

  5. N is the total number of hydrogen and oxygen atoms
    Questions I have:

  6. Is the temperature equation I’m using correct?

  7. Should KE be derived from just oxygen atom velocity and water molecule mass or by current method?

  8. What am I misunderstanding that leads to matlab producing 2/3 lammps result?

Attached are the results I’m obtaining and the matlab script I used along with the lammps run file (although the simulations appear to run without any issues).

Any advice you can offer is greatly appreciated.

Best regards,

Sam Shilliday
The University of Edinburgh is a charitable body, registered in Scotland, with registration number SC005336.

TemperaturePlot.png

run.int.nvt (1.45 KB)

water_temp.m (4.5 KB)

Dear Lammps Users,

I am conducting a project using lammps involving evaporation of water inside nanopores. One of the measurements I'm trying to take is the temperature of the water molecules (tip4p). To do this I have written a Matlab code but the answers have been off and running a test simulation of only 1000 water molecules shows answers are exactly 2/3 of the actual temperature.

Some assumptions I've made:

Temperature is set using 'fix fx tip4p nvt temp 300.0 300.0 200.0' command
Lammps temperature is given by 'thermo_style custom step temp' command
Matlab temperature is given by: 1/2 m v^2 = 3/2 N k T
Kinetic energy is evaluated from both hydrogen and oxygen atom velocities
N is the total number of hydrogen and oxygen atoms

Questions I have:

Is the temperature equation I'm using correct?

no.

Should KE be derived from just oxygen atom velocity and water molecule mass or by current method?

current.

What am I misunderstanding that leads to matlab producing 2/3 lammps result?

you are not counting degrees of freedom correctly. you have 9 DOF with
3 atoms per water, but due to using fix shake, you are removing 3
DOFs, so you are left with 6. 6/9 is the same as 2/3. also you have to
remove an additional 3 DOFs globally, since your system is invariant
to translation of the total system.

axel.