NVT simulation of rigid co2 molecules

LAMMPS LAMMPS Beginners
1

Hello,

I am simulating 8 molecules of CO2 using fix rigid/nvt/small at T=338K.
The mean temperatures calculated from thermo_style output and
from the total kinetic energy ( compute ke ) agree with each other, however
are 8% greater than the imposed temperature.
I have used 37 degrees of freedom to convert the kinetic energy to temperature
(=5 dof/molecule *8 molecules - 3).
I have included my input files ( CO2.txt and simulation.in) and the log file (simProg.log).

Where am I going wrong with the temperature calculation ?

Cheers,
Thejas

simProg.log (325.4 KB)
simulation.in (3.2 KB)
CO2.txt (408 Bytes)

2

This looks reasonable, but there are two issues and one concern:

  • Your system is tiny. Temperature is not well defined for atomic scale systems, but the smaller the system the more it will fluctuate and be subject to finite size effects.
  • You are not doing any time averaging. You either need to split your system in two parts where the first part is for equilibration and the second for collecting data after equilibrium is reached (note that this also applies to rotation versus translation of the rigid molecules, which is not so easy to see if you only look at the total temperature based on atoms), or you use a sufficiently large sliding window time averaging and wait until the numbers are properly stabilized.
  • Your timestep of 1fs may be a bit on the large side for your system. The rotational part of rigid body integrators is far more sensitive to divergence by too large a timestep with numerical integration that the translational part. So, if this was a water molecule where you can use fix shake, you might even have a 2fs timestep, but with fix rigid (also the water) would need 0.25fs or perhaps 0.5fs for a stable and energy conserving time integration.
3

Hi Axel,

many thanks for your reply.

I tried system sizes of 8, 27 and 125 molecules.
Indeed, the calculated temperature, using time averaging over a time window,
approaches the imposed value with increasing system size.
As suggested, I introduced a longer initialization phase
before the production phase, and the timestep was always below 0.25 fs.

I have attached the updated input file and progress file.
Thank you.

Cheers,
Thejas

meanTemp.dat (7.3 KB)
simProg.log (188.1 KB)
simulation.in (3.4 KB)

4

8 molecules is definitely to small. Due to periodic boundaries the molecule to the left and the right of a given molecule is the same molecule. So there are intrinsic restrictions from that because the motions are not independent.
I would expect that anything below 1000 molecules will exhibit such finite size effects (more indirect of course for the larger systems).

I did some systematic studies of finite size effects with water some 20+ years ago and found that you could see finite size effects on different system properties for up to 1024 molecules. At the time using 512 water molecules for simulating a pure bulk water system was considered state of the art. Time has moved on and computers have become much more powerful. So by running simulations with so few molecules you are setting yourself up to be in serious trouble once you want to publish any results from those runs.

5

Many thanks for your comments, Axel.

Point noted. The tiny system size only served as a starting point
to check the correctness of the input script.

Cheers,
Thejas