Re: [External] Re: Hydrogen bond energy for TIP3 water

From: Randy J. Zauhar (r.zauhar_at_usciences.edu)
Date: Wed Feb 05 2020 - 19:13:41 CST

Nice answer Kim, thank you!

Randy

On 5Feb, 2020, at 6:45 PM, Sharp, Kim <sharpk_at_pennmedicine.upenn.edu<mailto:sharpk_at_pennmedicine.upenn.edu>> wrote:

CAUTION: This email originated from outside the organization. Do not click links or open attachments unless you recognize the sender and know the content is safe.



Randy,
the original TIP3P was paramterized to reproduce the enthalpy of vaporizatoin of water:
(
Comparison of simple potential functions for simulating liquid water
J Chem. Phys. Jorgensen, W. L. Chandrasekhar, J. Madura, J. D.Impey, R. W. Klein, M. L. 1983

at 300K this is about 10.5 kcal/mol. if you assume that on average 2 bonds/water are broken (2=4/2 to
avoid double counting), then that would be about -5.2 kcal/mole/hbond, v. close to the -5 kcal/mol
you quoted. but none of these hbonds would be close to the
pair energy minimum due to thermal fluc. and bulk effects. Whether that could account for the loss of about
4 cal/mole I don’t know. Also original ’TIP3P’ parameters may have been tweaked in the CHARMM files
over the years to play nice with proteins. it would not take much change in sigma to have a large effect.
best
Kim Sharp
On Feb 5, 2020, at 4:58 PM, Victor Kwan <vkwan8_at_uwo.ca<mailto:vkwan8_at_uwo.ca>> wrote:

Dear Randy,

1. The LJ parameters in the CHARMMFF files is in the form of \epsilon (kcal/mol) and Rmin/2 (Angstrom)

2. The Lorentz-Berthelot rules is used to combine LJ parameters for a pair of heteroatom

3. \sigma is where the inter-particle potential is zero (c.f. Rmin, where the potential is the minimum). It is related to Rmin by: Rmin = 2^(1/6) \sigma

4. With 2. and 3., you get the TIP3P parameters on https://lammps.sandia.gov/doc/Howto_tip3p.html

5. Without knowing the exact detail of your calculation - I can't comment on whether the numbers are correct or not. "Stabilization energy" is always a quantity with respect to some reference state.

Kind regards,

Victor Kwan


On Wed, Feb 5, 2020 at 3:35 PM Randy J. Zauhar <r.zauhar_at_usciences.edu<mailto:r.zauhar_at_usciences.edu>> wrote:
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



“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"






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 : Fri Dec 31 2021 - 23:17:08 CST