Ice not melting - TIP4P/Ice

LAMMPS LAMMPS Mailing List Mirror
1

Hi everyone

I’m trying to build a simulation of ice Ih with the ice optimised parameters of the TIP4P potential (https://doi.org/10.1063/1.1931662 and https://lammps.sandia.gov/doc/Howto_tip4p.html), however, I get melting temperatures around 360 K at 1 atm instead of the close to 273 K described in the paper. I’ve tried various tweaks and twists to get closer to the 273 K but without success. I’ve among others tried varying the time step between 0.1 and 2 fs and run up to 2000000 runs.

Can you pinpoint what I’ve missed or done wrong? Or explain why the results from the TIP4P/Ice paper cannot be reproduced?

Thanks!

I’ve attached the input structure (ice_Ih.lmp) and a screenshot of it (ice_Ih.png).

Here is my input script:

units real

atom_style full

read_data ice_Ih.lmp

include tip4p_ice.potential

neighbor 2.0 bin

neigh_modify every 1 check yes

timestep 0.5

variable Nrun equal 1000000

variable Ndump equal 500

variable Text_start equal 300

variable Text_end equal 300

variable Pext_start equal 1.0

variable Pext_end equal 1.0

velocity all create ${Text_start} 122334 dist gaussian mom yes rot yes

fix contrain all shake 0.0001 100 0 b 1 a 1

thermo 100

thermo_style one

thermo_modify lost warn norm yes flush yes

dump trj all atom ${Ndump} ice_ih.lammpstrj

fix integrate all npt temp {Text_start} {Text_end} 100.0 iso {Pext_start} {Pext_end} 1000.0

run ${Nrun}

write_restart restarted.end

######### END ########

And the tip4p_ice.potential file:

TIP4P/ice for LAMMPS. Parameters set according to

https://lammps.sandia.gov/doc/Howto_tip4p.html

group hydrogen type 1

group oxygen type 2

set group oxygen charge -1.1794

set group hydrogen charge 0.5897

mass 1 1.008 # H

mass 2 15.9994 # O

pair_style lj/cut/tip4p/long 2 1 1 1 0.1577 8.5

pair_modify tail yes

kspace_style pppm/tip4p 1.0e-5

pair_coeff 1 1 0.0 0.0

pair_coeff 1 2 0.0 0.0

pair_coeff 2 2 0.21084 3.1668

bond_style harmonic

bond_coeff 1 0.0 0.9572

angle_style harmonic

angle_coeff 1 0.0 104.52

Yours sincerely,
Tor Blaabjerg Sørensen

ice_Ih.png

ice_Ih.lmp (590 KB)

2

Can you please elaborate a bit on how you are determining the melting point.

Axel.

3

Hi Alex

So far I’ve determined the melting point by visual inspection of the dump-file to see whether the hexagonal crystal structure is preserved or ‘collapsed’ into an unordered liquid phase.

A typical trial would be to raise the temperature from 150 K to, say, 320 K over 2 000 000 steps of 0,5 fs, 1 atm pressure, and then an additional run of 2 000 000 times 0,5 fs at 1 atm pressure and constant temperature of 320 K.

In a case like that, at the end of the simulation, the initial structure is essentially unaltered except for increased vibrations and some molecules briefly jumping out of their lattice points and back again. I am aware that the transition would be slow just around the melting point, but it concerns me that I cannot see any clear signs of melt after a simulation as described above.

Thank you for taking your time,

Tor Blaabjerg Sørensen
University College London

4

you are missing an important point here. melting (or freezing) is an activated process and requires nucleation. this can cause a significant hysteresis (you can have supercooled water that is liquid significantly below the melting point and also raise the temperature beyond the melting point and the system won’t melt for a very long time). this happens at the macroscopic scale (ever put a bottle of soda water into the freezer and see what happens when you open it? DON’T use a glass bottle for that kind of experiment!!), but especially with very small system sizes as in MD simulations, where you have a significant finite size effect.

to avoid this, you have to do a so-called coexistence simulation, where you have a half liquid, half solid system and then monitor which part grows and which part shrinks and the melting point is when they remain the same. LAMMPS has several structural computes categorizing atoms according to their immediate environment that may be applicable to tell the difference between liquid and solid.

so to start such a coexistence simulation, you have to set up and equilibrate a solid system at a temperature just under the expected melting point.
then apply a thermostat and time integration to only half of your system with a temperature way above the melting point, to bring the system to the liquid phase and then carefully cool it down need the melting point. then you can continue with a set of simulations starting from that restart at slightly different temperatures around the melting point to determine the temperature where the two parts of the system remain of the same size.

this is all straightfoward and has been discussed on this mailing list multiple times. it should also be discussed in the relevant published literature.

axel.

5

Thank you for the explanation and inspiration for home office experiments, I’ll try your suggested approach.

Tor Blaabjerg Sørensen