Temperature equilibrating at wrong value

Dear all,

I am trying to simulate a crystal system of water (hydrates) using TIP4P model with SHAKE along with some gas molecules present inside the crystal. My system is 3D periodic and I am using long range coulomb as well as LJ potentials. The water and gas molecules interact via repulsive part of LJ potential only, by truncating the potential at sigma.

I am using Noose-Hoover thermostat and barostat while doing the NPT simulation at different set of temperature and pressure. However, for any pressure, there are two things happening. In some cases temperature of the system is equilibrating at the specified value, but in some cases the temperature is equilibrating at exactly 10 degrees higher than whatever temperature specified.

One major difference in the two cases is that the simulation where equilibrium temperature is incorrect has more number of gas molecules in the crystal than the simulations where equilibrium temperature is correct.

Can someone help me understand what may be wrong or where I should look?

Thanks in advance
Chaitanya

If you post your input here, we can try to take a look at it (no promises; you know more about your system than any of us). Please remember to specify the LAMMPS version that you are using (the first line of your log file will tell you), and enclose the text of your input in triple backticks

```
like this
```

so that the forum software will not misinterpret any special characters and garble the formatting.

Thank you for responding @srtee .

I am using Lammps version of September-2021

Here is the input script I am using:

# created by fftool

units real
boundary p p p

atom_style full
bond_style harmonic
angle_style harmonic

# remove hybrid if not necessary
pair_style hybrid lj/long/tip4p/long long long 1 2 1 1 0.1577 8.5 lj/cut 3.70
kspace_style pppm/disp/tip4p 1.0e-5
kspace_modify force/disp/real 0.0001
kspace_modify force/disp/kspace 0.002

read_data data.lmp

pair_coeff    1    1 lj/long/tip4p/long     0.210831     3.166800  # O O
pair_coeff    1    2 lj/long/tip4p/long     0.000000     1.000000  # O H
pair_coeff    1    3 lj/cut                 30.00000     3.700000  # O Ar
pair_coeff    2    2 lj/long/tip4p/long     0.000000     1.000000  # H H
pair_coeff    2    3 lj/cut                 0.000000     1.000000  # H Ar
pair_coeff    3    3 lj/cut                 0.000000     3.700000  # Ar Ar

fix SHAKE all shake 0.0001 20 0 b 1 a 1

neighbor 2.0 bin
neigh_modify delay 0 every 1 check yes

timestep 1.0

variable TK equal &&&&
variable PATM equal @@@@

velocity all create ${TK} 12345

fix TPSTAT all npt temp ${TK} ${TK} 100 iso ${PATM} ${PATM} 1000

thermo 1000
thermo_style custom step etotal ke pe evdwl ecoul elong temp press density

variable T equal temp
fix TAVG all ave/time 50 1000 50000 v_T file AvgT

run 200000

unfix TAVG

reset_timestep 0

dump TRAJ all custom 300 &&&&K@@@@A.lammpstrj mol type q x y z
dump_modify TRAJ element O H Ar sort id

run 900000

write_data data.eq.lmp

As I have mentioned earlier, there is no difference in the input files, but only in the number of Argon molecules. In the incorrectly equilibrated system, the number of Argon molecules is twice as that in correctly equilibrated system (64 and 128 respectively). Here are some details relating to that:

### Data file for system where T is correctly equilibrated.
created by fftool

3328 atoms  
2176 bonds
1088 angles
3 atom types
1 bond types
1 angle types
0.000000 34.620000 xlo xhi
0.000000 34.620000 ylo yhi
0.000000 34.620000 zlo zhi

Masses

   1   15.999  # O
   2    1.008  # H
   3   39.948  # Ar

Bond Coeffs

   1   517.630258     0.957200  # O-H

Angle Coeffs

   1    37.950526   104.520000  # H-O-H

AND

### Data file for system where T is INcorrectly equilibrated.

created by fftool

3392 atoms  
2176 bonds
1088 angles
3 atom types
1 bond types
1 angle types
0.000000 34.620000 xlo xhi
0.000000 34.620000 ylo yhi
0.000000 34.620000 zlo zhi

Masses

   1   15.999  # O
   2    1.008  # H
   3   39.948  # Ar

Bond Coeffs

   1   517.630258     0.957200  # O-H

Angle Coeffs

   1    37.950526   104.520000  # H-O-H

I hope you can find some headway.

Thanks

My most immediate comment is that the Ar-O potential is both very “hard” (a much larger energy scale than the O-O potential) and sharply force-discontinuous at the cutoff you have defined.

(I am not sure if you intended to use a WCA potential, but bear in mind LAMMPS’s \sigma is defined as the radius where repulsive and attractive energy cancel out – not the radius of minimum energy, which is 2^{1/6} \sigma – so a cutoff r_c = \sigma gives force discontinuity at the cutoff.)

With that in mind you may simply have a far-too-large timestep, where 1fs works when the argon-oxygen collisions are fairly rare but not when they become more frequent. Check for conservation of energy (use thermo keyword econserve to account for thermostat / barostat energy transfer).

Thanks @srtee . I will get back to you upon further analysis.

Chaitanya