Hydrogen bond energy for TIP3 water

From: Randy J. Zauhar (r.zauhar_at_usciences.edu)
Date: Wed Feb 05 2020 - 14:34:38 CST

OK, I can’t find any thread that explicitly addresses this.

I gave students in my class a simple exercise - compute interaction energy between two water molecules in a standard orientation that would represent a ’strong’ linear hydrogen bond, all atoms in the same plane.

One water is to be translated versus the other, to get energy as a function of oxygen-oxygen separation. The idea is to estimate the stabilization of a single hydrogen bond from the curve, and to see relative contributions of Lennard-Jones and electrostatic energy.

Part of the assignment was to locate TIP3 parameters from a file I had on hand, par_all36_cgenff.parm.

The values are
OGTIP3: eps = 0.151 Rmin = 1.7682
HGTIP3: eps = 0.046 Rmin = 0.2255

As LJ is represented in CHARMm, Rmin is supposed to be separation of the same atom types at minimum energy, so these seem on the small side. They are also commented out in the file, which I did not notice at first,

Still, using those parameters with charges of -0.8 on O and +0.4 on H, I get a total minimum energy of about -20 kcal/mol! That magnitude seems really large.

Scouting around I found parameters for TIP3P water here:
https://lammps.sandia.gov/doc/Howto_tip3p.html

These express Rmin as ’sigma’, and give parameters explicitly for pairs of atom types:
O charge = -0.834
H charge = 0.417
LJ epsilon of OO = 0.1521
LJ sigma of OO = 3.1507
LJ epsilon of HH = 0.0460
LJ sigma of HH = 0.4000
LJ epsilon of OH = 0.0836
LJ sigma of OH = 1.7753

The ’sigmas’ almost correspond to treating the TIP3 Rmin as a ‘van der Waals’ radius to be added, as in some force fields, but not quite…

Anyway, using the TIP3P parameters (including the charges) I get a minimum energy of -9 kcal/mol.

That is better but is still almost double what I expect from ‘conventional wisdom’ (i.e. a strong h-bond provides about -5 kcal/mol of stabilization).

I had the students use a spreadsheet implementation, I also checked it in python (I don’t have a convenient way to get a single-point energy for this from Maestro or MOE).

Is this well understood? Something I am overlooking or misinterpreting?

Expressions I assume -

E(electro) = 332 * q1 * q2 / r12 (q in elementary charge units, r12 in Å, energy in kcal/mol)
E(LJ) = epsilon * [ (Rmin/r12)^12 - 2*(Rmin/r12)^6 ] (again, r12 in Å and E assumed in kcal/mol)

All atoms in molecule 1 see all atoms in molecule 2.

Thanks in advance,

Randy


Randy J. Zauhar, PhD

Prof. of Biochemistry

Dept. of Chemistry & Biochemistry
University of the Sciences in Philadelphia
600 S. 43rd Street
Philadelphia, PA 19104

Phone: (215)596-8691
FAX: (215)596-8543
E-mail: r.zauhar_at_usciences.edu<mailto:r.zauhar_at_usciences.edu>



“Yeah the night is gonna fall, and the vultures will surround you /
And when you’re lookin’ in the mirror what you see is gon’ astound you"

  — Death Cab for Cutie, “Monday Morning"



This archive was generated by hypermail 2.1.6 : Thu Dec 31 2020 - 23:17:12 CST