Problem with 3-site hydrogen molecules

LAMMPS LAMMPS Mailing List Mirror
1

Dear all lammps users

I’m doing simulation with TIP4P water, OPLS-UA methane, and 3-site hydrogen molecules.

The 3-site hydrogen model I’m using is from J. Chem. Phys. 123, 024507 (2005), which have massless atom at the center of hydrogen molecule with charges -0.9864e.

From paper, they used rigid hydrogen models. Also from searching, i found that lammps cannot support zero mass atom.

So I added fix rigid command to hydrogen molecules with center mass = 1e-10 and -0.9864e charge.

But simulation have some errors, so I found that hydrogen molecules move so fast and system explode.

There is no atom overlapping, I had seriously checked my .dat files. So maybe I think there is a problem with pair or bond interaction.

Also I’ve tried fix rigid/npt, but there is some errors as well.

I just add my input, and data files.

If you guys can give me some advices, I would very appreciate that.

Thank you

Dong woo Kang

<<in.test>>

Initial section

units real
dimension 3
boundary p p p
atom_style full
pair_style lj/cut/tip4p/long 5 4 1 1 0.1250 10
kspace_style pppm/tip4p 1.0e-4
bond_style harmonic
angle_style harmonic

#Atom definition section
read_data test.dat

#Setting section
pair_coeff 1 1 0.294 3.73
pair_coeff 1 2 0 0
pair_coeff 1 3 0.141564001 3.366264993
pair_coeff 1 4 0 0
pair_coeff 1 5 0.248971806 3.43688289
pair_coeff 2 2 0 0
pair_coeff 2 3 0 0
pair_coeff 2 4 0 0
pair_coeff 2 5 0 0
pair_coeff 3 3 0.068164511 3.038
pair_coeff 3 4 0 0
pair_coeff 3 5 0.119882466 3.101731516
pair_coeff 4 4 0 0
pair_coeff 4 5 0 0
pair_coeff 5 5 0.21084 3.1668

group Water type 4 5

group hydro type 2 3
fix 1 Water shake 0.0001 100 0 a 1 b 1
fix 2 hydro rigid/small molecule

thermo 1
thermo_style custom step time etotal ke temp pe ebond ecoul elong press density vol lx ly lz
dump 1 all custom 1 trajectory.lammpstrj.* id mol mass element x y z vx vy vz fx fy fz q
dump_modify 1 element M H1 MSH H2 O
neigh_modify delay 0 every 1 check yes

fix 3 all npt temp 220.0 220.0 100 iso 500 500 100
timestep 1
run 50000000

<<test.dat>>

1746 atoms
1042 bonds
521 angles

5 atom types
2 bond types
2 angle types

20 70 xlo xhi
-10 35 ylo yhi
-10 35 zlo zhi

Masses

1 16.042
2 1.008
3 0.0000000001
4 1.008
5 15.999

Bond Coeffs

1 517.6 0.95719
2 500 0.7414

Angle Coeffs

1 37.95 104.52
2 40 180

2

Try a smaller time step. 1fs is likely too large. Axel.

3

Dear Axel,

Thank you for your advice.

I’ve set timestep to 0.0001, but still there is some errors with hydrogen molecules.

Forces on massless hydrogen center are very large, so system is soon exploding.

Are there anything that I can try?

Dong Woo Kang

2019년 1월 25일 (금) 오후 4:18, Axel Kohlmeyer <akohlmey@…24…>님이 작성:

4

Dear Axel,

Thank you for your advice.

I've set timestep to 0.0001, but still there is some errors with hydrogen molecules.

Forces on massless hydrogen center are very large, so system is soon exploding.

Are there anything that I can try?

i don't know. but what you describe sounds like there could be a
problem with your force field parameters or definition.

axel.

5

Dear axel,

During trial and error, I found that there is no error when I set the ghost center atom mass = 1, but when setting below 0.1 (including 0.1), ghost atoms are separated from their molecules and missing.

And, when simulating only with hydrogen, (no water or methane) same error occurs like separating.

Is it problem with initial configuration? when I set simulation box tightly, system goes to NaN.

But when setting box size loosely, ghost atom missing.

I kept using fix rigid/small molecule command.

Cheers

Dong Woo Kang

2019년 1월 25일 (금) 오후 5:32, Axel Kohlmeyer <[email protected]>님이 작성:

6

Dear axel,

During trial and error, I found that there is no error when I set the ghost center atom mass = 1, but when setting below 0.1 (including 0.1), ghost atoms are separated from their molecules and missing.

And, when simulating only with hydrogen, (no water or methane) same error occurs like separating.

Is it problem with initial configuration? when I set simulation box tightly, system goes to NaN.

But when setting box size loosely, ghost atom missing.

I kept using fix rigid/small molecule command.

difficult to say what is happening from remote without knowing more
details. possible causes, that are consistent with your descriptions
are:

1) incorrect molecule ids in your data file, so that your molecules
are not detected correctly as rigid bodies
2) bad force field parameters, so you have point charges that can get too close
3) incorrect/inconsistent (e.g. overlapping) time integration. this is
the most likely option. have a look at your output, there should be
warnings about atoms being time integrated multiple times.

it is *definitely* possible to use very small masses of 1.0e-10 amu
for correct rigid body time integration. lots of people - including me
- have done this.

axel.

7

Dear axel,

Thank you for your advices.

I checked my input files thoroughly and find causes with time integration, and it works now.

cheers

Dong Woo Kang

2019년 1월 28일 (월) 오후 5:46, Axel Kohlmeyer <[email protected]>님이 작성: