# Water model functions
# Caroline Arnold 2018
import numpy as np
import CDTK.Tools.Conversion as cv
import CDTK.Tools.Mathematics as mat
[docs]
def water_to_point_charge( r_xyz, types, mode='SPC' ):
"""
Parametrize a set of water molecules by point charges
Parameters:
r_xyz - Water XYZ geometries
types - Atom types in r_xyz
mode - Parametrization to be applied [default: SPC]
Returns:
point charge xyz / c
"""
# 3 site water model
pc_xyz = r_xyz
pc_charges = np.zeros(len(types), dtype='float')
for i in range(pc_charges.size):
if types[i] == 'O':
pc_charges[i] = waterModel[mode]['q_O']
elif types[i] == 'H':
pc_charges[i] = waterModel[mode]['q_H']
else:
raise ValueError("Unrecognized atom type")
# 4 site water model
# Sort water molecules together
# Define M position
# Assign xyz positions and point charges
return pc_xyz, pc_charges
[docs]
def lennard_jones_potential( r, mode='SPC' ):
"""
Construct Lennard-Jones potential
V_LJ = 4 * eps * ( (sigma / r)**12 - (sigma / r)**6)
Parameters:
r - distance between atoms
mode - Water model [default: SPC]
Returns:
Lennart-Jones potential
"""
eps = waterModel[mode]['eps']
sigma = waterModel[mode]['sigma']
return 4. * eps * ( (sigma / r)**12 - (sigma / r)**6)
[docs]
def lennard_jones_gradient( r_vec, mode='SPC' ):
"""
Calculate analytical gradient from Lennard-Jones type interaction
F_LJ = 4 * eps * ( 12*(sigma / r)**12 - 6*(sigma / r)**6) * e_r / r
Parameters:
r_vec - vector between atoms
mode - water model [default: SPC]
Returns:
Analytical gradient
"""
eps = waterModel[mode]['eps']
sigma = waterModel[mode]['sigma']
r = np.linalg.norm(r_vec)
e_r = r_vec / r # force direction
return 4. * eps * ( 12. * ( sigma / r )**12 - 6. * ( sigma / r )**6 ) * e_r / r
[docs]
def coulomb_potential( r, q_i, q_j ):
"""
Construct Coulomb potential
r - distance [a.u.]
q_i,q_j - charges [a.u.]
return
Coulomb potential
"""
# SI: k = 1. / 4. / cv.PI / cv.FPC_EPSILON_0
return q_i * q_j / r
[docs]
def morse_potential(r, mode='TIP4P/2005f'):
"""
Construct Morse potential for intramolecular bond stretching
"""
D_r = waterModel[mode]['D_r']
beta = waterModel[mode]['beta']
r_eq = waterModel[mode]['r_eq']
e = D_r * (1 - np.exp( -beta*(r - r_eq )))**2
g = 0.0
return e, g
[docs]
def angle_bending_potential(theta, mode='TIP4P/2005f'):
"""
Construct angle bending potential for intramolecular angle bending
"""
K_theta = waterModel[mode]['K_theta']
theta_eq = waterModel[mode]['theta_eq']
e = 0.5 * K_theta * (theta - theta_eq)**2
g = 0.0
return e, g
[docs]
def get_m_site_xyz(r_xyz, mode='TIP4P/2005f'):
"""
Return M-site coordinates for given water molecule
The M site is located along the bisector of the dihedral angle, coplanar with the molecule
Parameters:
r_xyz - H2O xyz coordinates, ordered as OHH
mode - water model to be applied [default: TIP4P/2005f]
"""
oh1 = r_xyz[1] - r_xyz[0]
oh2 = r_xyz[2] - r_xyz[0]
rot_axis = np.cross( oh1, oh2 )
rot_axis /= np.linalg.norm(rot_axis)
alpha = 0.5 * np.arccos( np.dot(oh1, oh2) / np.linalg.norm(oh1) / np.linalg.norm(oh2) )
RM = mat.rodriguesMatrix( rot_axis, alpha )
om = np.dot( RM, oh1 )
om /= np.linalg.norm(om)
if 'd_OM_rel' in waterModel[mode]:
m_xyz = r_xyz[0] + om * waterModel[mode]['d_OM_rel']
else:
raise ValueError('M-site requires at least 4-site water model')
return m_xyz
# =====================================================
# Dictionary of water model parameters
#
# Format
# waterModel[mode]
#
# Contains
# q_O : partial charge on oxygen site [e]
# q_H : partial charge on hydrogen site [e]
# eps : Lennard-Jones potential well depth [eV]
# sigma : Lennard-Jones vanishing distance [an]
# r_eq : Distance oxygen-hydrogen [an]
# theta_eq : Angle hydrogen-oxygen-hydrogen [deg]
# d_OM_rel : 4-site model distance
# D_r : bond strength [ kJ / mol]
# beta : bond width [nm-1]
# K_theta : angle bending strength constant [kJ / mol / rad^2]
# TODO Units: a.u. throughout
# TIP4P/2005 and TIP4P/2005f: J Chem Phys 135, 224516 (2011)
waterModel = {
'SPC' : { 'q_O' : -0.82, # TODO citation needed
'q_H' : +0.41,
'eps' : 0.1544*cv.ev2au,
'sigma' : 3.16556*cv.an2au,
'r_eq' : 1.0*cv.an2au,
'theta_eq' : 109.47
},
'TIP4P/2005' : { 'eps' : 93.2*cv.FPC_K_B_EV,
'sigma' : 3.1589,
'q_H' : 0.5564,
'd_OM_rel' : 0.13194,
'r_eq' : 0.9572,
'theta_eq' : 104.52
},
'TIP4P/2005f' : { 'eps' : 93.2 * cv.FPC_K_B_EV * cv.ev2au,
'sigma' : 3.1644 * cv.an2au,
'q_H' : 0.5564,
'q_O' : -1.1128,
'd_OM_rel' : 0.13194 * cv.an2au,
'r_eq' : 0.9419 * cv.an2au,
'theta_eq' : 107.4 * np.pi / 180.0,
'D_r' : 432.581 * cv.kj2ev * cv.ev2au,
'beta' : 22.87 / (cv.an2au * 10.0),
'K_theta' : 367.810 * cv.kj2ev * cv.ev2au
}
}