special_bonds command issue

Hi, Lammps users,

I simulate SPC/E water behavior on a rigid wall, but experienced weired results with respect to “special_bonds” command. (sorry for long explanation below)

I think: since H-O-H are bonded, it is unnecessary to compute the lj/coulombic forces between bonded atoms in each atoms. I delete the Lennard-Jones & Coulombic force between H and O bonded in each water molecules, and also the H-H interactions angled in the same water molecules. I only consider the interactions between atoms that belong to different molecules.

However i got some weird results: the water did not move at all after they achieve equilibrium at T=300k; the temperature is almost 300K; no water molecule escape.

I also tried simulation without special_bonds command, and it seems the water molecules are moving randomly as imagined at equilibrium at T=300K, which is more close to common sense.

I checked the mannual for details: The rationale for using these weighting factors is that the interaction between a pair of bonded atoms is all (or mostly) specified by the bond, angle, dihedral potentials, and thus the non-bonded Lennard-Jones or Coulombic interaction between the pair of atoms should be excluded (or reduced by a weighting factor). The 1st of the 3 coefficients (LJ or Coulombic) is the weighting factor on 1-2 atom pairs, which are pairs of atoms directly bonded to each other. The 2nd coefficient is the weighting factor on 1-3 atom pairs which are those separated by 2 bonds (e.g. the two H atoms in a water molecule). The 3rd coefficient is the weighting factor on 1-4 atom pairs which are those separated by 3 bonds (e.g. the 1st and 4th atoms in a dihedral interaction). Thus if the 1-2 coefficient is set to 0.0, then the pairwise interaction is effectively turned off for all pairs of atoms bonded to each other. If it is set to 1.0, then that interaction will be at full strength.

My thinking follows the above command explanation. But it seems this command did not show what i want.

My question is how to set the weighting factor for H-O-H which did not influence the results? Thanks!

The details of input are as following:

group hy type 1 #H of H2O

group ox type 2 #O of H2O

group Mg type 3

group O type 4

group water type 1 2

group boundary type 3 4

pair_style lj/cut/coul/long 10 10

kspace_style pppm 0.0001

pair_coeff 1 1 0 0 #hy-hy

pair_coeff 2 2 0.1554 3.166 #ox-ox

pair_coeff 3 3 9.03E-07 5.909 #Mg-Mg

pair_coeff 4 4 0.1554 3.166 #O2-O2

pair_modify mix arithmetic #shift yes

special_bonds lj/coul 0.0 0.0 0.5

velocity boundary set 0 0 0 units box

velocity water create 300 49875678 mom yes rot yes dist gaussian units box

fix 1 boundary setforce 0 0 0

fix 2 boundary rigid single

fix 3 water shake 1.00E-04 20 0 b 1 a 1 #H-O and H-O-H

fix NVT water nvt temp 300 300 100

Paul can probably look at this …

Steve

I don’t see anything obviously wrong with your input script, but you’ve probably got an error in there (or in your data file) somewhere. The special_bonds command you’ve got there seems fine. I wouldn’t expect it to work properly without the special_bonds command. I’ll paste below some sample input commands I’ve used to simulate water. This might be helpful for you in debugging your input.

Paul

units real

neigh_modify delay 4 every 1

atom_style full

bond_style harmonic

angle_style charmm

pair_style lj/charmm/coul/long 8 10

pair_modify mix arithmetic

kspace_style pppm 1e-4

read_data water_equil.data

special_bonds charmm

velocity all create 400.0 12345678 dist uniform

fix 1 all shake 1e-6 500 0 m 1.0 a 1

fix 2 all npt 308.0 308.0 100.0 xyz 1.0 1.0 1000.0

thermo 100

timestep 2.0

#dump 1 all custom 100 water.dump tag type x y z fx fy fz

restart 1000 water.res

run 1000