calculating temperature using tip4p water model in lammps

LAMMPS LAMMPS Mailing List Mirror
1

Hello,

I asked a question regarding tip4p earlier, but this is a different topic, so hopefully people don’t mind me starting a separate thread for this question. I am trying to study the lammps code, and am having issues figuring out how the fictitious m site in the tip4p model is being thermostatted in lammps, and included when the temperature is calculated. For instance, I had lammps dump the position and velocity of a 2 molecule tip4p water model after setting the velocity to 217.00 (I.E. lammps says the velocity is 217.00 K)

ITEM: ATOMS id x y z vx vy vz
1 ‘O’ 149.649 147.785 151.388 0.00267316 0.00126655 0.00451636
2 ‘H’ 149.706 148.061 152.303 0.010903 0.0256255 0.00155614
3 ‘H’ 149.776 148.59 150.886 -0.0080722 -0.0111786 -0.00203906
5 ‘O’ 149.27 150.273 149.327 -0.00464517 -0.0137007 -0.0079375
4 ‘H’ 150.01 150.065 149.897 0.00398954 -0.00106168 9.97245e-05
6 ‘H’ 150.763 150.025 149.308 0.0163819 -0.00855383 -0.0191179

If I take assume the number of constraints to be 9 (each water contributes 2 bond constraints and one angle constraint, and 3 constraints for removing velocity center of mass), then I can only calculate 217.00 for the temperature if I assume there are 8 atoms in the system, so the m atoms are included in the number of particles count, and the m atoms contribute nothing to the KE (zero mass). Is this correct, or is lammps representing the m positions as constraint in some way, since the bond distance is restrained to 0.15?

For the langevin thermostat, how is the thermostat handling the m site with no mass? The langevin thermostat will force all atoms to sample from the correct velocity distribution, but because the M site has zero mass and hence always zero KE, I feel it would lead to a temperature below the set point, if something were not adjusted compared to say a TIP3P simulation. Could someone help me out with my understanding?

2

Hello,

I asked a question regarding tip4p earlier, but this is a different topic,
so hopefully people don't mind me starting a separate thread for this
question. I am trying to study the lammps code, and am having issues
figuring out how the fictitious m site in the tip4p model is being
thermostatted in lammps, and included when the temperature is calculated.

the fictitious site does not show up anywhere, but in contributing to
forces on atoms and it cannot contribute to temperature, since it has
no mass.
another way to see this is that the virtual site is simply a
modification to the potential energy function of the water atoms that
makes the interactions non-spherical.

For instance, I had lammps dump the position and velocity of a 2 molecule
tip4p water model after setting the velocity to 217.00 (I.E. lammps says
the velocity is 217.00 K)

ITEM: ATOMS id x y z vx vy vz
1 'O' 149.649 147.785 151.388 0.00267316 0.00126655 0.00451636
2 'H' 149.706 148.061 152.303 0.010903 0.0256255 0.00155614
3 'H' 149.776 148.59 150.886 -0.0080722 -0.0111786 -0.00203906
5 'O' 149.27 150.273 149.327 -0.00464517 -0.0137007 -0.0079375
4 'H' 150.01 150.065 149.897 0.00398954 -0.00106168 9.97245e-05
6 'H' 150.763 150.025 149.308 0.0163819 -0.00855383 -0.0191179

If I take assume the number of constraints to be 9 (each water contributes 2
bond constraints and one angle constraint, and 3 constraints for removing
velocity center of mass), then I can only calculate 217.00 for the
temperature if I assume there are 8 atoms in the system, so the m atoms are
included in the number of particles count, and the m atoms contribute

no.

nothing to the KE (zero mass). Is this correct, or is lammps representing
the m positions as constraint in some way, since the bond distance is
restrained to 0.15?

no.

For the langevin thermostat, how is the thermostat handling the m site with
no mass? The langevin thermostat will force all atoms to sample from the
correct velocity distribution, but because the M site has zero mass and
hence always zero KE, I feel it would lead to a temperature below the set
point, if something were not adjusted compared to say a TIP3P simulation.
Could someone help me out with my understanding?

most likely there is something in the script you ran, that is not
entirely consistent. for example, you cannot completely trust the
velocity command for systems with SHAKE. as the assigned velocities
initially do not include corrections for velocities cancelling due to
the imposed constraints during velocity verlet timestepping. after a
few steps of MD, the temperature *should* be consistent. also, you
have to make sure that you actually do remove any center of mass
velocity in your system.

the computation of DOFs seems right. 3*6 is 18 DOFs, minus 2*3 per
water plus 3 global, gives you a remaining 9 DOFs.

axel.