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

From: Sharp, Kim (sharpk_at_pennmedicine.upenn.edu)
Date: Wed Feb 05 2020 - 17:45:57 CST

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

This archive was generated by hypermail 2.1.6 : Fri Dec 31 2021 - 23:17:08 CST