How to set the COMPASS class2 potential parameter for water?

Dear lammps-users,

I’m now calculating the dielectric constant of water via lammps, which the flexible water potentials are better suited. The attached paper(http://journals.aps.org/prb/abstract/10.1103/PhysRevB.31.2643) by K. Toukan and A. Rahman is really a good flexible spc water model, which introducing an anharmonic Morse potential for the O-H stretch vibrational part, so that I want to use this model.

According to the equ(3) from the above mentioned paper, I checked the lammps website and found that the bond_style morse and angle_style class2 can meet my requests for using a flexible water model. My input file for this coefficient part is following, under the real unit.

pair_coeff * * 0.0000 1.000

pair_coeff 1 1 0.1554253 3.165492 #OW OW

pair_coeff 1 2 0.0000 0.00000 #OW HW

pair_coeff 2 2 0.0000 0.00000 #HW HW

bond_style morse # D(dissociation energy), alpha(inverse distance), r0

bond_coeff 1 109.6797323 0.196314751 1.0 #OW-HW

angle_style class2 # Ea(angle term), Ebb(bond-bond term), Eba(bond-angle term)

angle_coeff 1 109.5 0.3805 0 0 # Ea: theta0, K2, K3, K4

angle_coeff 1 bb 2.4139579349904396e-23__/NA__ 1.0 1.0 # Ebb: bb, M, r1 r2

angle_coeff 1 ba 0.228 0.228 1.0 1.0 # Eba: ba N1 N2 r1 r2

Because of the units of the parameters are not right, I always get the error that “Bond atom missing …”. I checked some related papers but still not got the right unit. So does anyone can help me to fix the angle_coeff or bond_coeff? The aforementioned NA is the Avogadro Constant, if it’s right to divide NA to get the right coefficient? Or anyone could tell me if I did it right?

Any suggestions are appreciated very much. Thanks for your time.

Best,

Zhu Liu

PhysRevB.31.2643.pdf (289 KB)

Real units for energy are Kcal/mole. Which means how much
energy in Kcals there would be if you had Avogadro’s # of bonds.

That should be enough info to convert any energy value they
give you in the paper to “real” units for LAMMPS.

Steve