Issue with TIP4P Flexible Potential for water

[David_Blanton]

Hi,

I am trying to run a simulation to equilibrate H2O at 1 K in a box that is larger than the periodic boundaries of the system via the “change_box” command so that the box will have enough space for additional atoms. When I extend the box dimensions the system is equilibrating properly with energy but when uploading the coordinates into a visual tool (OVITO) the system is seen to melt. Why is this happening? My input file is below…My ultimate goal is to equilibrate H2O with SiC to calculate the equilibrium properties of the interface between the two structures.

INPUT FILE:

units real

dimension 3

boundary p p p

atom_style full

read_data data.H20.dat

#read_data H2O_648N_Coord.txt #648

#variables-------------------------------------------------------------------------

variable Text equal 1.0

variable Pext equal 0.0

variable Nrun equal 1000000

variable ts equal 0.0016

variable Tdamp equal ${ts}*100

variable Pdamp equal ${ts}*1000

variable Nf equal ${Nrun}/100

variable Ne equal 10

variable Nr equal {Nf}/{Ne}

variable Ndump equal ${Nrun}/2

variable Nr_rdf equal 0.5*{Nrun}/{Ne}

variable watMoleMass equal 18.0153 # /(g/mol)

variable nAvog equal 6.0221415e23 # Avogadro’s number

variable watMoleculeMass equal ({watMoleMass}/{nAvog}) # /(g/molecule)

variable A3_in_cm3 equal 1e-24 # Angstrom^3 in cm^3

variable nAtoms equal atoms

variable nMolecules equal v_nAtoms/3

#create groups ###---------------------------------------------------------

mass 1 1.00794 # H

mass 2 15.9994 # O

#include forcefield.TIP4P-2005.txt

TIP4P/2005 potential parameters-------------------------

group O type 2

group H type 1

#group H2O type 1 2

TIP4P/2005 flexible potenrial parameters---------------------------------

pair_style lj/cut/tip4p/long 2 1 1 1 0.1546 12.0

kspace_style pppm/tip4p 1e-5

pair_coeff 1 1 0.0 0.0

pair_coeff 1 2 1 0.0 1.5795

pair_coeff 2 2 0.1852573718 3.1644

bond_style morse

bond_coeff 1 103.38934 2.287 0.9419

angle_style harmonic

angle_coeff 1 43.95435 107.4

replicate 6 6 6

change_box all z delta -10.0 10.0 x delta -10.0 10.0 y delta -10.0 10.0 units box

neighbor-----------------------------------------------------------------

neighbor 2.0 bin

neigh_modify delay 0 every 1 check yes

#Run one minimization to settle the lattice

fix 1 all nve

fix 2 all box/relax iso 0.0 vmax 0.001

min_style cg

minimize 1.0e-12 1.0e-12 100000 1000000

#print the entire lattice (used for SIA and Vac search) -#but we only wanted printed out now

run 1

#undump startdump

unfix 1

#dumps and restart-----------------------------------------------------------

initial condition------------------------------------------------------

velocity all create ${Text} 1234546 dist gaussian mom yes rot yes

#dump waterdump all custom 1000 lammps.coord_vel id type x y z vx vy vz

dump trj all atom ${Nf} WaterTrj.data

dump_modify trj scale no sort id

dump deqlb all atom ${Nrun} FinalCoord.data

dump_modify deqlb scale no sort id

dump coord all atom ${Ne} DumpCoord.txt

dump_modify coord scale no sort id

timestep ${ts}

fix 4 all npt iso {Pext} {Pext} {Pdamp} temp {Text} {Text} {Tdamp}

variable Dens equal v_nMolecules*{watMoleculeMass}/(vol*{A3_in_cm3})

fix DensAve all ave/time {Ne} {Nr} ${Nf} v_Dens file wat.dens.data

output-------------------------------------------------------------------

thermo_style custom step temp pe ke etotal press vol

thermo_modify flush yes norm yes

thermo ${Nf}

run ${Nrun}

write_restart restart.H2O_16p

Thanks

David Blanton

[Ray_Shan]

Do you mean that your water melted? So you started with some sort of ice structures? I believe it is quite difficult to capture ice, and it is even more difficult when you have some vacuum space. When the atoms have more spaces to relax to, they may not want to stay in their original positions.

Ray

[Andrew_Jewett]

I may be missing the point of your question, but, as Ray pointed out, if your goal was to start with ice, you do need to pick the initial coordinates carefully. One way to do that is to create a DATA file. A DATA file gives you a great deal of flexibility to move the atoms to the correct locations.

If you are comfortable with the “moltemplate” tool, then there is an example for generating SPCE ice located in this directory:
moltemplate/examples/all_atom_examples/force_field_explicit/ice_crystal/
It contains the positions of the oxygen and hydrogen atoms only (not the 4th particle, but I assume you don’t need that). If you want to use one of the TIP4P water models, it should not be difficult to change the pair_style and atom charge (by editing the “spce.lt” file, which the other .lt files refer to, and renaming it to “tip4p.lt”. You can get the the parameters from http://lammps.sandia.gov/doc/Section_howto.html#howto_8). If you get this working, please send me your “tip4p.lt” file and I’ll post it for others.

It is difficult to define an exact crystal structure for ice, because 1) it has multiple phases (I assume you want “Ih” ice), AND 2) because there are additional degrees of freedom due to the positions of the hydrogen atoms.
For this reason, there are several variants of the unit cell in the moltemplate/common/ directory:
spce_ice_rect8.lt, spce_ice_rect16.lt, spce_ice_rect32.lt (…with 8, 16, and 32 water molecules, respectively). They differ only in the (amount of disorder in the) direction of the hydrogen bonds within the crystal. The 8-molecule unit cell is shown below.
The lattice constants for the ice crystal were taken from this paper:
Rottger et al.,
“Lattice constants and thermal expansion of H2O and D2O ice Ih between 10 and 265K”,
Acta Crystallogr. B, 50 (1994) 644-648

…as explained in the example, you make an array of unit cells to define larger ice crystals

cells = new SpceIceRect8 [3].move(4.521, 0.0, 0.0)
[2].move( 0.0, 7.832, 0.0)
[2].move( 0.0, 0.0, 7.362)

(see the “system.lt” file in the ice_crystal/moltemplate_files/ directory) This generates:

​

The lower-left corner of the box will be at the origin. (0,0,0) If you want to move the water molecules elsewhere, one way is to append a command like this:

cells[0-2][0-1][0-1].move(0,0,20.0)

Of course, the crystal you create need not fill the volume of the simulation box.

Non-rectangular crystals:
You can use the “delete” command (see chapter 7.10 of the moltemplate manual), and/or use non-orthogonal unit cell to create a non-rectangular ice crystal.

Note: This particular example uses the data from this paper, taken at 265K & 100Torr. Since you want a more realistic ice structure at 1K, you may want to adjust the distances between the molecules. In that case you don’t have to start from scratch. To make everything 10% smaller (for example), use this command instead:

cells = new SpceIceRect8.scale(0.9,0.9,0.9)
[3].move(4.0689, 0.0, 0.0)
[2].move( 0.0, 7.0488, 0.0)
[2].move( 0.0, 0.0, 6.6258)

(Notice that the distances in the move() commands were shrunk as well)

My apologies for posting yet another advertisement for moltemplate.

Andrew