Temperature of nitrogen is approaching zero when using a thermostat

LAMMPS LAMMPS Beginners
1

Dear LAMMPS Experts,

I hope this message finds you well.

I apologize in advance if this question may seem too basic, but I would greatly appreciate your help.

I am currently working on a simulation where I am trying to bring a liquid water droplet surrounded by nitrogen gas into equilibrium at 300K. However, I have encountered an issue where the temperature of the nitrogen particles approaches zero, and their motion appears to freeze, as seen in the trajectory files.

To set up this simulation, I based my approach on the instructions from the “Howto discussion: 8.4.4. TIP4P water model,” with the addition of nitrogen particles. I have made modifications to the code accordingly.

Simulation codes: (Using LAMMPS (29 Aug 2024 - Update 1))

units real
atom_style full
region box block -200 200 -200 200 -200 200
create_box 3 box  bond/types 2 angle/types 1 &
                extra/bond/per/atom 2 extra/angle/per/atom 1 extra/special/per/atom 2

mass 1 15.9994   #O
mass 2 1.008     #H
mass 3 14.007    #N

pair_style lj/cut/tip4p/cut 1 2 1 1 0.15 8.0
pair_coeff 1 1 0.1550 3.1536
pair_coeff 2 2 0.0    1.0
pair_coeff 3 3 0.166991 3.71  #N-N
pair_modify shift yes mix arithmetic

bond_style zero
bond_coeff 1 0.9574  # H2O
bond_coeff 2 1.10      # N2
angle_style zero
angle_coeff 1 104.52

molecule water tip3p.mol
molecule n2 n2.mol
region liquid sphere 0 0 0 35
region gas sphere 0 0 0 40 side out

create_atoms 0 random 5982 34564 liquid mol water 25367 overlap 3.0 #8930
create_atoms 0 random 1596 34564 gas mol n2 25367 overlap 2.0 #670

group liquid region liquid
group gas region gas
group water type 1 2
group n2 type 3

timestep 1.0
fix rigid all shake 0.0001 10 0 b 1 2 a 1

minimize 1E-10 1E-10 1000 10000

reset_timestep 0

velocity all create 300.0 12345

run 0

fix 1 all nvt/omp temp 300 300 100

compute water_temp water temp
compute n2_temp n2 temp

thermo_style custom step temp c_water_temp c_n2_temp

run 2000000

write_restart 1215_water_ready.restart

Molecule files:

# N2 molecule

2 atoms
1 bonds
0 angles
0 dihedrals
0 impropers

Coords

1  0.0   0.0   0.0
2  1.1   0.0   0.0

Types

1  3
2  3

Charges

1        0.0
2        0.0

Bonds

1  2  1  2


# Water molecule. TIP3P geometry

3 atoms
2 bonds
1 angles

Coords

1    0.00000  -0.06556   0.00000
2    0.75695   0.52032   0.00000
3   -0.75695   0.52032   0.00000

Types

1        1   # O
2        2   # H
3        2   # H

Charges

1       -1.040
2        0.520
3        0.520

Bonds

1   1      1      2
2   1      1      3

Angles

1   1      2      1      3

Shake Flags

1 1
2 1
3 1

Shake Atoms

1 1 2 3
2 1 2 3
3 1 2 3

Shake Bond Types

1 1 1 1
2 1 1 1
3 1 1 1

Special Bond Counts

1 2 0 0
2 1 1 0
3 1 1 0

Special Bonds

1 2 3
2 1 3
3 1 2

With above codes, when i was trying to track temperature of nitrogen, by using c_n2_temp, i got the following thermo outputs:

image
image609×775 93.5 KB

Would anyone be able to help me identify any possible errors or suggest adjustments that could resolve this issue? I would greatly appreciate any insights or guidance.

Thank you very much for your time and support!

2

Hello,

I don’t see anything blatantly wrong with your input (the cut-off is a bit short, but that’s probably not the issue here).

Did you try observing your system using VMD/ovito? Does it resemble what you expect?

Simon

3

Well this is not quite a flying ice cube (a flying … steam drop?), but if you search “flying ice cube” you will see that, very broadly, thermal energy can “concentrate” in certain modes of a system when thermalization between parts of the system is insufficient.

In this case, thermostatting both the liquid water and nitrogen gas as one system results in thermal energy concentrating in the liquid (presumably due to higher heat capacity and phonon frequencies).

You will have to experiment with different settings (and consult the literature), but I would use an NVT integrator for the liquid water and an NVE integrator for the nitrogen gas (which would be more physical too). Additionally, I would modify the water’s NVT integrator to use a temperature compute that subtracts out centre of mass motion (fix_modify temp with compute temp/com).

4

Thank you so so much for your thoughtful and meaningful reply!

I never imagined I would encounter such an issue, but your guidance has given me a new perspective. I will continue to investigate this problem diligently and explore possible solutions.

Thank you once again for your invaluable support! :laughing: