# * **************************************************************************
# *
# * CDTK, Chemical Dynamics Toolkit
# * A modular system for chemical dynamics applications and more
# *
# * Copyright (C) 2018
# * Ralph Welsch, DESY, <ralph.welsch@desy.de>
# * Copyright (C) 2021
# * Ludger Inhester, DESY, <ludger.inhester@desy.de>
# *
# * Copyright (C) 2020, 2021, 2022, 2023
# * Ludger Inhester, DESY, ludger.inhester@cfel.de
# *
# * This file is part of CDTK.
# *
# * CDTK is free software: you can redistribute it and/or modify
# * it under the terms of the GNU General Public License as published by
# * the Free Software Foundation, either version 3 of the License, or
# * (at your option) any later version.
# *
# * This program is distributed in the hope that it will be useful,
# * but WITHOUT ANY WARRANTY; without even the implied warranty of
# * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
# * GNU General Public License for more details.
# *
# * You should have received a copy of the GNU General Public License
# * along with this program. If not, see <http://www.gnu.org/licenses/>.
# *
# * **************************************************************************
import numpy as np
import CDTK.Tools.Conversion as cv
import CDTK.Dynamics.Constraints as cnstr
import warnings
import re
import json
def _angleDER(r1, r2, r3):
"""
for the three cartesian coordinates r1,r2,r3
returns angle theta between r1-r2 and r3-r2 and
derivative of theta with respect to r1,r2,r3
Input:
r1,r2,r3 numpy arrays dimension 3
Output: theta, dTheta_dR1, dTheta_dR2, dTheta_dR3
theta: angle
dTheta_dR1, dTheta_dR2, dTheta_dR3: derivatives of theta with respect to r1,r2,r3, respectively. 3dim. numpy arrays.
"""
r12 = r1 - r2
r32 = r3 - r2
# d12 = np.linalg.norm(r12)
d12 = np.sqrt(r12[0]**2 + r12[1]**2 + r12[2]**2)
# d32 = np.linalg.norm(r32)
d32 = np.sqrt(r32[0]**2 + r32[1]**2 + r32[2]**2)
oneOverD12d32 = 1/(d12*d32)
# cosTheta = np.dot(r12, r32) / (d12 * d32)
cosTheta = (r12[0]*r32[0] + r12[1]*r32[1] + r12[2]*r32[2]) * oneOverD12d32
theta = np.arccos(cosTheta)
if (abs(abs(cosTheta)-1.0) < 1e-8):
dTheta_dCosTheta = 0.
else:
dTheta_dCosTheta = -1. / np.sqrt(1. - cosTheta**2)
dCosTheta_dR12 = r32 * oneOverD12d32 - r12 * (cosTheta / d12**2)
dCosTheta_dR32 = r12 * oneOverD12d32 - r32 * (cosTheta / d32**2)
dTheta_dR1 = dTheta_dCosTheta * dCosTheta_dR12
dTheta_dR3 = dTheta_dCosTheta * dCosTheta_dR32
dTheta_dR2 = - dTheta_dR1 - dTheta_dR3
return theta, dTheta_dR1, dTheta_dR2, dTheta_dR3
def _cosAngleDER(r1, r2, r3):
"""
for the three cartesian coordinates r1,r2,r3
returns cos theta between r1-r2 and r3-r2 and
derivative of cos theta with respect to r1,r2,r3
Input:
r1,r2,r3 numpy arrays dimension 3
Output: costheta, dcosTheta_dR1, dcosTheta_dR2, dcosTheta_dR3
cosTheta: cos angle
dcosTheta_dR1, dcosTheta_dR2, dcosTheta_dR3: derivatives of cos theta with respect to r1,r2,r3, respectively. 3dim. numpy arrays.
"""
r12 = r1 - r2
r32 = r3 - r2
d12 = np.sqrt(r12[0]**2 + r12[1]**2 + r12[2]**2)
d32 = np.sqrt(r32[0]**2 + r32[1]**2 + r32[2]**2)
oneOverD12d32 = 1/(d12*d32)
# cosTheta = np.dot(r12, r32) / (d12 * d32)
cosTheta = (r12[0]*r32[0] + r12[1]*r32[1] + r12[2]*r32[2]) * oneOverD12d32
dCosTheta_dR12 = r32 * oneOverD12d32 - r12 * (cosTheta / d12**2)
dCosTheta_dR32 = r12 * oneOverD12d32 - r32 * (cosTheta / d32**2)
dCosTheta_dR1 = dCosTheta_dR12
dCosTheta_dR3 = dCosTheta_dR32
dCosTheta_dR2 = - dCosTheta_dR1 - dCosTheta_dR3
return cosTheta, dCosTheta_dR1, dCosTheta_dR2, dCosTheta_dR3
def _angle(r1, r2, r3):
"""
for the three cartesian coordinates r1,r2,r3
returns angle theta between r1-r2 and r3-r2
Input:
r1, r2, r3: position vectors as numpy arrays
Output:
theta: angle
"""
r12 = r1 - r2
r32 = r3 - r2
# d12 = np.linalg.norm(r12)
d12 = np.sqrt(r12[0]**2 + r12[1]**2 + r12[2]**2)
# d32 = np.linalg.norm(r32)
d32 = np.sqrt(r32[0]**2 + r32[1]**2 + r32[2]**2)
cosTheta = np.dot(r12, r32) / (d12 * d32)
theta = np.arccos(cosTheta)
return theta
def _dihedralAngle(r1, r2, r3, r4):
"""
for the four cartesian coordinates r1,r2,r3,r4
returns dihedral angle phi and derivative of phi with respect to r1,r2,r3,r4
Input:
r1,r2,r3,r4 numpy arrays dimension 3
Output: phi, dPhi_dR1, dPhi_dR2, dPhi_dR3, dPhi_dR4
phi: dihedral angle
dPhi_dR1, dPhi_dR2, dPhi_dR3, dPhi_dR4: derivatives of phi with respect to r1,r2,r3,r4, respectively. 3dim. numpy arrays.
"""
r12 = r1 - r2
r23 = r2 - r3
r34 = r3 - r4
# d12 = np.linalg.norm(r12)
d23 = np.linalg.norm(r23)
# d34 = np.linalg.norm(r34)
c1 = np.cross(r23, r34)
c2 = np.cross(r12, r23)
c3 = np.cross(r34, r12)
p1 = np.dot(r12, c1) * d23
p2 = np.dot(c1, c2)
phi = np.arctan2(p1, p2)
dPhi_dP1 = p2/(p1**2+p2**2)
dPhi_dP2 = -p1/(p1**2+p2**2)
# P1 = r12 . (r23xr34) * d23
# =r34 . (r12xr23) * d23
# =r23 . (r34xr12) * d23
dP1_dR12 = c1 * d23
dP1_dR34 = c2 * d23
dP1_dR23 = c3 * d23 + np.dot(r23, c3) * r23 / d23
# P2 = (r23 x r34 ) (r12 x r23)
# = ( r23 r12 ) (r34 r23) - (r23 r23) (r34 r12)
dP2_dR12 = r23 * np.dot(r34, r23) - d23**2 * r34
dP2_dR23 = r12 * np.dot(r34, r23) + \
np.dot(r23, r12) * r34 - 2 * r23 * np.dot(r34, r12)
dP2_dR34 = r23 * np.dot(r23, r12) - d23**2 * r12
#
dPhi_dR12 = dPhi_dP1 * dP1_dR12 + dPhi_dP2 * dP2_dR12
dPhi_dR23 = dPhi_dP1 * dP1_dR23 + dPhi_dP2 * dP2_dR23
dPhi_dR34 = dPhi_dP1 * dP1_dR34 + dPhi_dP2 * dP2_dR34
dPhi_dR1 = dPhi_dR12
dPhi_dR2 = dPhi_dR23 - dPhi_dR12
dPhi_dR3 = dPhi_dR34 - dPhi_dR23
dPhi_dR4 = -dPhi_dR34
return phi, dPhi_dR1, dPhi_dR2, dPhi_dR3, dPhi_dR4
def _EpsSigmaToC6C12(eps, sigma):
"""
converts eps, sigma to C6 C12 in kJ/mol^n
Input:
eps : epsilon in a.u.
sigma: sigma in a.u.
Output
C6: C6 in kJ/mol nm^6
C12: C12 in kJ/mol nm^12
"""
sigma /= 10. * cv.an2au # to nm
eps /= cv.kj2au # to kJ/mol
sigma6 = np.power(sigma, 6)
C6 = 4 * sigma6 * eps
C12 = 4 * sigma6 * sigma6 * eps
return C6, C12
[docs]
class forcefield(object):
"""
Template for a general force field
"""
# TODO:
# Eventually, all force fields should be specified
# by _parameter structure.
# However, few details are not implemented generally
# but only for specific force fields
# - virtual sites is currently only used for TIP4P
# and is only implemented there
# - constraints are only implemented in SPC/E
# and only implemented specifically for SPC/E
# todo: check if introduce sigma_nuc here
_parameter = {
'name': 'dummy forcefield',
'charges': {},
'sigma_nuc': {}, # Width of Gaussian nuclear charge model, bohr
'masses': {},
'C6': {}, # Lennard-Jones Parameter, unit kJ/mol nm^6
'C12': {}, # Lennard-Jones Parameter, unit kJ/mol nm^12
'bonds': [], # for each bond dictionary with parameters
'angles': [], # for each angle dictionary with parameters
'dihedrals': [], # for each dihedral dictionary with parameters
'bond_constraints': [], # bond constraints
'angle_constraints': [] # angle constraints
}
eBonds = 0.
eAngles = 0.
eDihdedral = 0.
eCoul = 0.
eNonbonbded = 0.
def __init__(self, topol=None, **opts):
self.topol = [t.split("_")[1] for t in topol]
for atom in self._parameter['charges']:
if atom not in self.topol:
raise ValueError("Cannot find atom %s in topology" % atom)
for atom in self.topol:
if atom not in self._parameter['charges']:
raise ValueError("Cannot find atom %s in parameters" % atom)
if len(topol) != len(self._parameter['charges'].keys()):
raise ValueError("Number of atoms does not match! %s != %s" %
self._parameter['charges'].keys(), topol)
self.topolMap = {k: v for v, k in enumerate(self.topol)}
if "pairs" in self._parameter:
self.pairs = self._parameter["pairs"]
else:
bondlist = self._parameter["bonds"]
if 'bond_constraints' in self._parameter:
bondlist = bondlist + self._parameter["bond_constraints"]
self.pairs = {}
for bond in bondlist:
ai = bond["ai"]
aj = bond["aj"]
if not ai in self.pairs:
self.pairs[ai] = [aj]
else:
self.pairs[ai].append(aj)
if not aj in self.pairs:
self.pairs[aj] = [ai]
else:
self.pairs[aj].append(ai)
# prepare for the evaluation of nonBonded interaction within the molecule
self.bondSteps = {}
def findDistance(a0, ai, it):
for aj in self.pairs[ai]:
if aj in self.bondSteps[a0] and self.bondSteps[a0][aj] <= it:
continue
else:
self.bondSteps[a0][aj] = it
findDistance(a0, aj, it+1)
for ai in self.topol:
self.bondSteps[ai] = {ai: 0}
for ai in self.topol:
findDistance(ai, ai, 1)
nexcl = 3
if "nexcl" in self._parameter:
nexcl = self._parameter["nexcl"]
self.nonBondedPairs = np.array([(self.topolMap[a1], self.topolMap[a2])
for a1, b in self.bondSteps.items() for a2, d in b.items() if self.topolMap[a1] < self.topolMap[a2] and d > nexcl])
self.nonBondedPairsBondSteps = np.array([d
for a1, b in self.bondSteps.items() for a2, d in b.items() if self.topolMap[a1] < self.topolMap[a2] and d > nexcl])
fudgeLJ = 1.
fudgeQQ = 1.
if "fudgeLJ" in self._parameter:
fudgeLJ = self._parameter["fudgeLJ"]
if "fudgeQQ" in self._parameter:
fudgeQQ = self._parameter["fudgeQQ"]
self.nonBondedPairsFudgeLJ = np.where(
self.nonBondedPairsBondSteps == 4, fudgeLJ, 1.)
self.nonBondedPairsFudgeQQ = np.where(
self.nonBondedPairsBondSteps == 4, fudgeQQ, 1.)
[docs]
def name(self):
"""
returns the name of the force field
"""
return self._parameter["name"]
[docs]
def getLJCombRule(self):
"""
returns the combination rule
"""
if "combRule" in self._parameter:
return self._parameter["combRule"]
else:
return None
[docs]
def hasConstraints(self):
"""
returns whether forcefield has constraints
"""
l = 0
if "bond_constraints" in self._parameter:
l += len(self._parameter["bond_constraints"])
if "angle_constraints" in self._parameter:
l += len(self._parameter["angle_constraints"])
return l > 0
[docs]
def numConstraints(self):
"""
returns number of constraints
"""
l = 0
if "bond_constraints" in self._parameter:
l += len(self._parameter["bond_constraints"])
if "angle_constraints" in self._parameter:
l += len(self._parameter["angle_constraints"])
return l
[docs]
def getConstraints(self, geom, der=False):
"""
returns constraints
Input:
geom: geometry of molecule, numpy array (natoms,3)
der: (optionial) boolen Flag for wheter to provide dg as output
Output:
if der, output is g,dg
else, output is g
g: deviation in constraints, n-dim numpy array
dg: derivative of constraints with respect to coordinates, numpy array of shape (n,na,3)
"""
# gs = []
dgs = []
nb = len(self._parameter["bond_constraints"])
na = len(self._parameter["angle_constraints"])
gs = np.zeros(na+nb)
if geom is None:
return gs
dgs = np.zeros((na+nb, geom.shape[0], geom.shape[1]))
for i, bond_constraint in enumerate(self._parameter["bond_constraints"]):
ai = bond_constraint['ai']
aj = bond_constraint['aj']
v = bond_constraint['value'] * 10 * cv.an2au
r = geom[self.topolMap[ai], :] - geom[self.topolMap[aj], :]
# d = np.linalg.norm(r)
d = np.sqrt(r[0]**2 + r[1]**2 + r[2]**2)
gs[i] = d - v
if der:
dgs[i, self.topolMap[ai], :] = r/d
dgs[i, self.topolMap[aj], :] = -r/d
for i, angle_constraint in enumerate(self._parameter["angle_constraints"]):
ai = angle_constraint['ai']
aj = angle_constraint['aj']
ak = angle_constraint['ak']
v = angle_constraint['value'] * np.pi/180
if der:
theta, dTheta_dR1, dTheta_dR2, dTheta_dR3 = _angleDER(
geom[self.topolMap[ai], :],
geom[self.topolMap[aj], :],
geom[self.topolMap[ak], :])
gs[nb + i] = theta - v
dgs[nb + i, self.topolMap[ai], :] = dTheta_dR1
dgs[nb + i, self.topolMap[aj], :] = dTheta_dR2
dgs[nb + i, self.topolMap[ak], :] = dTheta_dR3
else:
theta = _angle(
geom[self.topolMap[ai], :],
geom[self.topolMap[aj], :],
geom[self.topolMap[ak], :])
gs[nb + i] = theta - v
if der:
return gs, dgs
else:
return gs
[docs]
def getCharges(self):
"""
Returns array of charges for each atom / virtual site.
"""
return np.array([float(self._parameter['charges'][s]) for s in self.topol])
[docs]
def getSigmaNuc(self, sigma_case):
"""
Returns array with the spread (sigma) of the nuclear charge for each atom in a gaussian nuclear model. Sigma is parsed
to the QCE.
\rho(r) = Z_i * (\frac{1}{2 \pi \sigma_i^2})^(2/3) * exp(-r^2 / (2 \sigma_i^2))
Parameters
----------
sigma_case: string or float, optional
If float, it is used as sigma for all atoms.
If string, it can be either 'covrad', use the
covalent radius of the atoms; or try to convert
the string to a dictionary of the form {'atom_name': sigma_value}.
Default is 'covrad'.
Returns
-------
sigma_nuc: (n_atoms,) ndarray of floats
"""
# todo: Add equation of GNDM to documentation as implemented in libcint
# Set same sigma for all atoms
try:
sigma_case = float(sigma_case)
return np.array([sigma_case for s in self.topol])
sigma_nuc = {s:sigma_case*cv.an2au for s in self.topol}
except (TypeError, ValueError):
# Set sigma for each atom (define a dictionary)
if sigma_case == 'covrad': # Use covalent radii
sigma_nuc = {
'OW': 0.66 * cv.an2au,
'HW1': 0.31 * cv.an2au,
'HW2': 0.31 * cv.an2au,
}
else: # Try to read a dictionary from input string
try:
sigma_nuc = {}
for i in range(0, len(sigma_case), 2):
sigma_nuc[sigma_case[i]] = float(
sigma_case[i+1]) / cv.au2an
except:
raise ValueError(
'The input string should be of the form `key1 value1 key2 value2 ...`'
)
return np.array([sigma_nuc[s] for s in self.topol])
[docs]
def getMasses(self):
# ! Check units of masses
"""
Returns array of masses for each atom
"""
return np.array([float(self._parameter['masses'][s]) * cv.am2au for s in self.topol])
[docs]
def getLJsigma(self):
"""
returns array of LJ sigma parameter for each atom
"""
sigmas = []
if "C6" in self._parameter and "C12" in self._parameter:
for a in self.topol:
if a in self._parameter["C6"]:
C6 = self._parameter["C6"][a]
C12 = self._parameter["C12"][a]
sigmas.append(np.power(C12/C6, 1./6) * 10. * cv.an2au)
else:
sigmas.append(0.)
elif "sigma" in self._parameter and "epsilon" in self._parameter:
for a in self.topol:
if a in self._parameter["sigma"]:
sigma = self._parameter["sigma"][a] * 10. * cv.an2au
sigmas.append(sigma)
else:
sigmas.append(0.)
else:
raise ValueError(f"no LJ parameters defined")
return np.array(sigmas)
[docs]
def getLJeps(self):
"""
returns array of LJ epsilon parameter for each atom
"""
epsilons = []
if "C6" in self._parameter and "C12" in self._parameter:
for a in self.topol:
if a in self._parameter["C6"]:
C6 = self._parameter["C6"][a]
C12 = self._parameter["C12"][a]
epsilons.append(1/4. * C6 * C6 / C12 * cv.kj2au)
else:
epsilons.append(0.)
elif "sigma" in self._parameter and "epsilon" in self._parameter:
for a in self.topol:
if a in self._parameter["epsilon"]:
epsilon = self._parameter["epsilon"][a] * \
cv.kc2ev * cv.ev2au
epsilons.append(epsilon)
else:
epsilons.append(0.)
else:
raise ValueError(f"no LJ parameters defined")
return np.array(epsilons)
[docs]
def hasVirtualSite(self):
"""
returns True if forcefield has a virtual site logic
"""
return False
[docs]
def getVGeom(self, geom, noM=False):
"""
Input:
geom: geometry as 2 dim numpy array geom[number of atoms, 3]
noM: boolean, whether virtual sites for this molecule should be skipped
Output: vGeom, gM
vGeom: geometry with virtual sites as 2 dim numpy array geom[number of atoms, 3]
gM: Jacobi matrix elements as list ot tuples (ai,i,bj,j,g)
ai: index of real atom
i: index of atom coord (0...2)
bj: index of virtual site
j: index of virtual site coord (0...2)
g: derviative dj / di
"""
gM = [(ai, i, ai, i, 1.0) for i in range(3)
for ai in range(geom.shape[0])]
return geom.copy(), gM
[docs]
def getBonded(self, molX, topol):
"""
input:
molX: geometry of the molecule as flat numpy array
topol: topology of the molecule
Output: energy,gradient
energy: energy of bonded interaction
gradient: gradient of the energy as flat numpy array
"""
x = molX.reshape(-1, 3)
gradient = np.zeros(x.shape)
energy = 0.
self.eBonds = 0.
for bond in self._parameter["bonds"]:
ai = bond["ai"]
aj = bond["aj"]
# harmonic potential
if bond["func"] == 1:
b = bond["b"] * 10 * cv.an2au
k = bond["K"] * cv.kj2au / (10*cv.an2au)**2
e, gi, gj = self._harmonicPotential(x[self.topolMap[ai]],
x[self.topolMap[aj]],
b, k)
# Gromos-96 type quartic potential
elif bond["func"] == 2:
b = bond["b"] * 10 * cv.an2au
k = bond["K"] * cv.kj2au / (10*cv.an2au)**4
e, gi, gj = self._fourthPowerPotential(x[self.topolMap[ai]],
x[self.topolMap[aj]],
b, k)
# Morse Potential
elif bond["func"] == 3:
b = bond["b"] * 10 * cv.an2au
beta = bond["beta"] * 1/(10 * cv.an2au)
D = bond["D"] * cv.kj2au
e, gi, gj = self._morsePotential(x[self.topolMap[ai]],
x[self.topolMap[aj]],
b, beta, D)
else:
raise ValueError(
'bond stretch potential type %d not implemented' % bond["func"])
energy += e
self.eBonds += e
gradient[self.topolMap[ai]] += gi
gradient[self.topolMap[aj]] += gj
self.eAngles = 0
for angle in self._parameter["angles"]:
ai = angle["ai"]
aj = angle["aj"]
ak = angle["ak"]
if angle["func"] == 1:
# G96 type angle harmonic cos-angle potential
k = angle["K"] * cv.kj2au # kJ/(mol rad^2) -> au/rad^2
ri = x[self.topolMap[ai]]
rj = x[self.topolMap[aj]]
rk = x[self.topolMap[ak]]
# deg -> rad.
cosTheta0 = np.cos(angle["theta"] * np.pi / 180.)
e, gi, gj, gk = self._cosAnglePotential(
ri, rj, rk, cosTheta0, k)
elif angle["func"] == 2:
# harmonic angle potential
theta = angle["theta"] * np.pi / 180. # deg -> rad.
k = angle["K"] * cv.kj2au # kJ/(mol rad^2) -> au/rad^2
ri = x[self.topolMap[ai]]
rj = x[self.topolMap[aj]]
rk = x[self.topolMap[ak]]
e, gi, gj, gk = self._harmonicAnglePotential(
ri, rj, rk, theta, k)
elif angle["func"] == 2:
# harmonic angle potential
theta = angle["theta"] * np.pi / 180. # deg -> rad.
k = angle["K"] * cv.kj2au # kJ/(mol rad^2) -> au/rad^2
ri = x[self.topolMap[ai]]
rj = x[self.topolMap[aj]]
rk = x[self.topolMap[ak]]
e, gi, gj, gk = self._harmonicAnglePotential(
ri, rj, rk, theta, k)
else:
raise ValueError(
'angle potential type %d not implemented' % angle["func"])
gradient[self.topolMap[ai]] += gi
gradient[self.topolMap[aj]] += gj
gradient[self.topolMap[ak]] += gk
energy += e
self.eAngles +=e
self.eDihedrals = 0
for dihedral in self._parameter["dihedrals"]:
ai = dihedral["ai"]
aj = dihedral["aj"]
ak = dihedral["ak"]
al = dihedral["al"]
ri = x[self.topolMap[ai]]
rj = x[self.topolMap[aj]]
rk = x[self.topolMap[ak]]
rl = x[self.topolMap[al]]
if dihedral["func"] in [1, 2, 4, 9]:
phiS = dihedral["phiS"] * np.pi/180. # deg -> rad
if dihedral["func"] == 1 or dihedral["func"] == 4 or dihedral["func"] == 9:
# periodic dihedrals
n = dihedral["n"]
k = dihedral["K"] * cv.kj2au # kJ/(mol) -> a.u.
e, gi, gj, gk, gl = self._periodicDihedralPotential(
ri, rj, rk, rl, phiS, n, k)
elif dihedral["func"] == 2:
# kJ/(mol rad^2) -> a.u./rad^2
k = dihedral["K"] * cv.kj2au
e, gi, gj, gk, gl = self._harmonicDihedralPotential(
ri, rj, rk, rl, phiS, k)
elif dihedral["func"] == 3:
C = np.array([dihedral["C0"],
dihedral["C1"],
dihedral["C2"],
dihedral["C3"],
dihedral["C4"],
dihedral["C5"]]) * cv.kj2au # kJ/(mol) -> a.u.
e, gi, gj, gk, gl = self._RyckaertBellemansdihedralPotential(
ri, rj, rk, rl, C)
else:
raise ValueError(
'dihedral potential type %d not implemented' % dihedral["func"])
gradient[self.topolMap[ai]] += gi
gradient[self.topolMap[aj]] += gj
gradient[self.topolMap[ak]] += gk
gradient[self.topolMap[al]] += gl
energy += e
self.eDihedrals += e
return (energy, gradient)
[docs]
def getNonBondedInteractionPairs(self):
"""
returns the atom pairs that experience non-bonded interaction (LJ + coulomb) within the molecule
"""
return self.nonBondedPairs
[docs]
def getNonBondedInteractionPairsFudgeLJ(self):
"""
returns the fudge factor non-bonded interaction (LJ) within the molecule
"""
return self.nonBondedPairsFudgeLJ
[docs]
def getNonBondedInteractionPairsFudgeQQ(self):
"""
returns the fudge factor non-bonded interaction (Coulomb) within the molecule
"""
return self.nonBondedPairsFudgeQQ
def _harmonicPotential(self, r1, r2, b, k):
"""
Input:
r1,r2: cartesian vectors numpy array dimension 3
b: eq. bond length
K: force constant
Output:
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
"""
dist = r1 - r2
# d = np.linalg.norm(dist)
d = np.sqrt(dist[0]**2 + dist[1]**2 + dist[2]**2)
e = k/2 * (d-b)**2
g1 = k * (d-b) * dist/d
g2 = -k * (d-b) * dist/d
return e, g1, g2
def _morsePotential(self, r1, r2, b, beta, D):
"""
Morse potential
Input:
r1,r2: cartesian vectors numpy array dimension 3
b: eq. bond length
beta: steepness
D: depth
Output:
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
"""
r12 = r1-r2
# d = np.linalg.norm(r12)
d = np.sqrt(r12[0]**2 + r12[1]**2 + r12[2]**2)
exp = np.exp(-beta * (d - b))
e = D * (1. - exp)**2
dvdr = 2.0 * D * (1.-exp) * exp * beta
g1 = r12 / d * dvdr
g2 = -r12 / d * dvdr
return e, g1, g2
def _fourthPowerPotential(self, r1, r2, b, k):
"""
Input:
r1,r2: cartesian vectors numpy array dimension 3
b: eq. bond length
K: force constant
Output:
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
"""
dist = r1 - r2
d2 = (dist[0]**2 + dist[1]**2 + dist[2]**2)
# d2 = np.dot(dist, dist)
b2 = b**2
e = k/4 * (d2-b2)**2
g = k * (d2-b2) * dist
g1 = g
g2 = -g
return e, g1, g2
def _harmonicAnglePotential(self, r1, r2, r3, theta0, k):
"""
harmonic potential for bond angles
Input:
r1,r2,r3: cartesian vectors numpy array dimension 3 of atoms
theta: eq. angle
K: force constant
Output: energy,g1,g2,g3
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
g3 gradient on atom 3, numpy array dimension 3
"""
theta, dTheta_dR1, dTheta_dR2, dTheta_dR3 = _angleDER(r1, r2, r3)
e = k/2 * (theta - theta0)**2
gTheta = k * (theta - theta0)
g1 = gTheta * dTheta_dR1
g2 = gTheta * dTheta_dR2
g3 = gTheta * dTheta_dR3
return e, g1, g2, g3
def _cosAnglePotential(self, r1, r2, r3, cosTheta0, k):
"""
harmonic potential for bond angles 1/2 k (cos(phi)-cos(phi_0))**2
Input:
r1,r2,r3: cartesian vectors numpy array dimension 3 of atoms
theta: eq. angle
K: force constant
Output: energy,g1,g2,g3
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
g3 gradient on atom 3, numpy array dimension 3
"""
cosTheta, dCosTheta_dR1, dCosTheta_dR2, dCosTheta_dR3 = _cosAngleDER(
r1, r2, r3)
e = k/2 * (cosTheta - cosTheta0)**2
gCosTheta = k * (cosTheta - cosTheta0)
g1 = gCosTheta * dCosTheta_dR1
g2 = gCosTheta * dCosTheta_dR2
g3 = gCosTheta * dCosTheta_dR3
return e, g1, g2, g3
def _periodicDihedralPotential(self, r1, r2, r3, r4, phiS, n, k):
"""
periodic dihedral potetial
Input:
r1,r2,r3,r4: cartesian vectors numpy array dimension 3 of 4 atoms
phiS: eq. dihedral angle
K: force constant
n: periodicity
Output: energy,g1,g2,g3,g4
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
g3 gradient on atom 3, numpy array dimension 3
g4 gradient on atom 4, numpy array dimension 3
"""
phi, dPhi_dR1, dPhi_dR2, dPhi_dR3, dPhi_dR4 = _dihedralAngle(
r1, r2, r3, r4)
e = k * (1. + np.cos(n*phi - phiS))
gPhi = k * -n * np.sin(n*phi - phiS)
g1 = gPhi * dPhi_dR1
g2 = gPhi * dPhi_dR2
g3 = gPhi * dPhi_dR3
g4 = gPhi * dPhi_dR4
return e, g1, g2, g3, g4
def _harmonicDihedralPotential(self, r1, r2, r3, r4, phiS, k):
"""
harmonic dihedral potetial
Input:
r1,r2,r3,r4: cartesian vectors numpy array dimension 3 of 4 atoms
phiS: eq. dihedral angle
K: force constant
Output: energy,g1,g2,g3,g4
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
g3 gradient on atom 3, numpy array dimension 3
g4 gradient on atom 4, numpy array dimension 3
"""
phi, dPhi_dR1, dPhi_dR2, dPhi_dR3, dPhi_dR4 = _dihedralAngle(
r1, r2, r3, r4)
e = 0.5 * k * (phi-phiS)**2
gPhi = k * (phi-phiS)
g1 = gPhi * dPhi_dR1
g2 = gPhi * dPhi_dR2
g3 = gPhi * dPhi_dR3
g4 = gPhi * dPhi_dR4
return e, g1, g2, g3, g4
def _RyckaertBellemansdihedralPotential(self, r1, r2, r3, r4, C):
"""
RyckaertBellemansdihedra dihedral potential
Input:
r1,r2,r3,r4: cartesian vectors numpy array dimension 3 of 4 atoms
C[:]: coefficient array of shape (5)
Output: energy,g1,g2,g3,g4
energy value
g1 gradient on atom 1, numpy array dimension 3
g2 gradient on atom 2, numpy array dimension 3
g3 gradient on atom 3, numpy array dimension 3
g4 gradient on atom 4, numpy array dimension 3
"""
phi, dPhi_dR1, dPhi_dR2, dPhi_dR3, dPhi_dR4 = _dihedralAngle(
r1, r2, r3, r4)
cosPsi = np.cos(phi - np.pi)
dCosPsi_dPhi = - np.sin(phi - np.pi)
e = C[0] + C[1] * cosPsi + C[2] * cosPsi**2 + C[3] * \
cosPsi**3 + C[4] * cosPsi**4 + C[5] * cosPsi**5
gCosPsi = C[1] + 2 * C[2] * cosPsi + 3 * C[3] * \
cosPsi**2 + 4 * C[4] * cosPsi**3 + 5 * C[5] * cosPsi**4
gPhi = gCosPsi * dCosPsi_dPhi
g1 = gPhi * dPhi_dR1
g2 = gPhi * dPhi_dR2
g3 = gPhi * dPhi_dR3
g4 = gPhi * dPhi_dR4
return e, g1, g2, g3, g4
def _testHarmonicAnglePotential(self):
theta = 80.*np.pi/180.
k = 1.
thetas = np.arange(0.1*np.pi, 1.9*np.pi, 0.05*np.pi)
energies = np.zeros(thetas.shape)
eps = 1e-5
grad = np.zeros((3, 3))
nGrad = np.zeros((3, 3))
for i in range(thetas.shape[0]):
r = np.array(
[[1, 0, 0], [0, 0, 0], [np.cos(thetas[i]), np.sin(thetas[i]), 0]])
e, g1, g2, g3 = self._harmonicAnglePotential(
r[0], r[1], r[2], theta, k)
energies[i] = e
grad[0] = g1
grad[1] = g2
grad[2] = g3
for ai in range(3):
for i in range(3):
rEps = r.copy()
rEps[ai, i] += eps
eEps, ge1, ge2, ge3 = self._harmonicAnglePotential(
rEps[0], rEps[1], rEps[2], theta, k)
nGrad[ai, i] = (eEps-e)/eps
assert (np.max(np.abs((grad-nGrad).flatten())) < 1e-2)
from matplotlib import pyplot as plt
plt.plot(thetas, energies)
plt.plot(thetas, 0.5*k*(thetas-theta)**2)
plt.show()
def _testDihedralPotentials(self):
phi0 = 180.*np.pi/180.
k = 1.
n = 2
phis = np.arange(0.1*np.pi, 1.9*np.pi, 0.05*np.pi)
dihedral = np.zeros(phis.shape)
energies = np.zeros(phis.shape)
energies2 = np.zeros(phis.shape)
eps = 1e-6
nGrad = np.zeros((4, 3))
nGrad2 = np.zeros((4, 3))
for i in range(phis.shape[0]):
r = np.array(
[
[0., 0., 0.], # x-y plane
[1., 0., 0.],
[1., 1., 0.], # x-z plane
[1. + np.cos(phis[i]), 1., np.sin(phis[i])] #
])
dihedral[i] = _dihedralAngle(r[0], r[1], r[2], r[3])[0]
e, g1, g2, g3, g4 = self._periodicDihedralPotential(
r[0], r[1], r[2], r[3], phi0, n, k)
grad = np.stack([g1, g2, g3, g4])
energies[i] = e
e2, g21, g22, g23, g24 = self._harmonicDihedralPotential(
r[0], r[1], r[2], r[3], phi0, k)
energies2[i] = e2
grad2 = np.stack([g21, g22, g23, g24])
for ai in range(4):
for j in range(3):
rEps = r.copy()
rEps[ai, j] += eps
eEps, ge1, ge2, ge3, ge4 = self._periodicDihedralPotential(
rEps[0], rEps[1], rEps[2], rEps[3], phi0, n, k)
nGrad[ai, j] = (eEps-e)/eps
eEps2, ge21, ge22, ge23, eg24 = self._harmonicDihedralPotential(
rEps[0], rEps[1], rEps[2], rEps[3], phi0, k)
nGrad2[ai, j] = (eEps2-e2)/eps
assert (np.max(np.abs((grad-nGrad).flatten())) < 1e-4)
assert (np.max(np.abs((grad2-nGrad2).flatten())) < 1e-4)
from matplotlib import pyplot as plt
plt.plot(dihedral*180./np.pi, energies)
plt.show()
# test with 4 different random coordinates
r = np.random.rand(4, 3)
# analytical gradients
e, g1, g2, g3, g4 = self._periodicDihedralPotential(
r[0], r[1], r[2], r[3], phi0, n, k)
grad = np.stack([g1, g2, g3, g4])
e2, g21, g22, g23, g24 = self._harmonicDihedralPotential(
r[0], r[1], r[2], r[3], phi0, k)
grad2 = np.stack([g21, g22, g23, g24])
for ai in range(4):
for j in range(3):
rEps = r.copy()
rEps[ai, j] += eps
eEps, ge1, ge2, ge3, ge4 = self._periodicDihedralPotential(
rEps[0], rEps[1], rEps[2], rEps[3], phi0, n, k)
nGrad[ai, j] = (eEps-e)/eps
eEps2, ge21, ge22, ge23, eg24 = self._harmonicDihedralPotential(
rEps[0], rEps[1], rEps[2], rEps[3], phi0, k)
nGrad2[ai, j] = (eEps2-e2)/eps
assert (np.max(np.abs((grad-nGrad).flatten())) < 1e-4)
assert (np.max(np.abs((grad2-nGrad2).flatten())) < 1e-4)
[docs]
class ureaFF(forcefield):
# Smith et al J. Phys. Chem. B, Vol. 108, No. 3, 2004
_parameter = {
'name': 'UREA GROMOS force field Smith et al J. Phys. Chem. B, Vol. 108, No. 3, 2004',
'charges': {
'OU': -0.39000,
'CU': 0.14200,
'N1U': -0.54200,
'H11U': 0.33300,
'H12U': 0.33300,
'N2U': -0.54200,
'H21U': 0.33300,
'H22U': 0.33300
},
'masses': {
'OU': cv.periodicTable['O']['atomic_mass'],
'CU': cv.periodicTable['C']['atomic_mass'],
'N1U': cv.periodicTable['N']['atomic_mass'],
'H11U': cv.periodicTable['H']['atomic_mass'],
'H12U': cv.periodicTable['H']['atomic_mass'],
'N2U': cv.periodicTable['N']['atomic_mass'],
'H21U': cv.periodicTable['H']['atomic_mass'],
'H22U': cv.periodicTable['H']['atomic_mass']
},
'C6': { # unit kJ/mol nm^6
'CU': 0.0048868488,
'OU': 0.0023639044,
'N1U': 0.0033527574,
'N2U': 0.0033527574,
},
'C12': { # unit kJ/mol nm^12
'CU': 1.3589545e-05,
'OU': 1.5898688e-06,
'N1U': 3.9509513e-06,
'N2U': 3.9509513e-06,
},
'bonds': [
{'ai': 'CU', 'aj': 'OU', 'func': 2, 'b': 0.1265, 'K': 1.31e7},
{'ai': 'CU', 'aj': 'N1U', 'func': 2, 'b': 0.135, 'K': 1.03e7},
{'ai': 'CU', 'aj': 'N2U', 'func': 2, 'b': 0.135, 'K': 1.03e7},
{'ai': 'N1U', 'aj': 'H11U', 'func': 2, 'b': 0.1, 'K': 1.87e7},
{'ai': 'N1U', 'aj': 'H12U', 'func': 2, 'b': 0.1, 'K': 1.87e7},
{'ai': 'N2U', 'aj': 'H21U', 'func': 2, 'b': 0.1, 'K': 1.87e7},
{'ai': 'N2U', 'aj': 'H22U', 'func': 2, 'b': 0.1, 'K': 1.87e7},
],
'angles': [
{'ai': 'OU', 'aj': 'CU', 'ak': 'N1U',
'func': 2, 'theta': 121.4, 'K': 690},
{'ai': 'N1U', 'aj': 'CU', 'ak': 'N2U',
'func': 2, 'theta': 117.42, 'K': 636},
{'ai': 'CU', 'aj': 'N1U', 'ak': 'H11U',
'func': 2, 'theta': 120.0, 'K': 390},
{'ai': 'CU', 'aj': 'N1U', 'ak': 'H12U',
'func': 2, 'theta': 120.0, 'K': 390},
{'ai': 'CU', 'aj': 'N2U', 'ak': 'H21U',
'func': 2, 'theta': 120.0, 'K': 390},
{'ai': 'CU', 'aj': 'N2U', 'ak': 'H22U',
'func': 2, 'theta': 120.0, 'K': 390},
{'ai': 'H11U', 'aj': 'N1U', 'ak': 'H12U',
'func': 2, 'theta': 120.0, 'K': 445},
{'ai': 'H21U', 'aj': 'N2U', 'ak': 'H22U',
'func': 2, 'theta': 120.0, 'K': 445},
],
'dihedrals': [
{'ai': 'OU', 'aj': 'CU', 'ak': 'N1U', 'al': 'H11U',
'func': 1, 'phiS': 180., 'n': 2, 'K': 41.8},
{'ai': 'OU', 'aj': 'CU', 'ak': 'N2U', 'al': 'H22U',
'func': 1, 'phiS': 180., 'n': 2, 'K': 41.8},
{'ai': 'OU', 'aj': 'N2U', 'ak': 'N1U', 'al': 'CU',
'func': 2, 'phiS': 0., 'K': 167.42309},
{'ai': 'CU', 'aj': 'H12U', 'ak': 'H11U', 'al': 'N1U',
'func': 2, 'phiS': 0., 'K': 167.42309},
{'ai': 'CU', 'aj': 'H22U', 'ak': 'H21U', 'al': 'N2U',
'func': 2, 'phiS': 0., 'K': 167.42309},
],
}
[docs]
class furfuralFF(forcefield):
_parameter = {
'name': 'OPLSAA parameters taken from S. Rabet and G. Raabe, Fluid Phase Equilibria 554, 113331 (2022)',
'fudgeLJ': 0.5,
'fudgeQQ': 0.5,
'nbfunc': 1,
'combrule': 'OPLS',
'moltype': [],
'nexcl': 3,
'charges': {'O1': -0.45, 'C2': 0.47, 'H3': 0.122, 'C4': -0.019, 'C5': -0.154, 'H6': 0.126, 'C7': -0.154, 'H8': 0.126, 'C9': -0.019, 'H10': 0.142, 'O11': -0.19},
'masses': {'O1': 15.9994, 'C2': 12.011, 'H3': 1.008, 'C4': 12.011, 'C5': 12.011, 'H6': 1.008, 'C7': 12.011, 'H8': 1.008, 'C9': 12.011, 'H10': 1.008, 'O11': 15.9994},
'epsilon': {'O1': 0.87864, 'C2': 0.43932, 'H3': 0.06276, 'C4': 0.29288, 'C5': 0.317984, 'H6': 0.12552, 'C7': 0.317984, 'H8': 0.12552, 'C9': 0.29288, 'H10': 0.12552, 'O11': 0.58576},
'sigma': {'O1': 0.296, 'C2': 0.375, 'H3': 0.242, 'C4': 0.355, 'C5': 0.355, 'H6': 0.242, 'C7': 0.355, 'H8': 0.242, 'C9': 0.355, 'H10': 0.242, 'O11': 0.29},
'bonds': [ {'ai': 'O1', 'aj': 'C2', 'func': 1, 'b': 0.1229, 'K': 476976.0},
{'ai': 'C2', 'aj': 'H3', 'func': 1, 'b': 0.109, 'K': 284512.0},
{'ai': 'C2', 'aj': 'C4', 'func': 1, 'b': 0.149, 'K': 334720.0},
{'ai': 'C4', 'aj': 'C5', 'func': 1, 'b': 0.1367, 'K': 456892.8},
{'ai': 'C4', 'aj': 'O11', 'func': 1, 'b': 0.136, 'K': 284512.0},
{'ai': 'C5', 'aj': 'H6', 'func': 1, 'b': 0.108, 'K': 307105.6},
{'ai': 'C5', 'aj': 'C7', 'func': 1, 'b': 0.1424, 'K': 392459.2},
{'ai': 'C7', 'aj': 'H8', 'func': 1, 'b': 0.108, 'K': 307105.6},
{'ai': 'C7', 'aj': 'C9', 'func': 1, 'b': 0.1367, 'K': 456892.8},
{'ai': 'C9', 'aj': 'H10', 'func': 1, 'b': 0.108, 'K': 307105.6},
{'ai': 'C9', 'aj': 'O11', 'func': 1, 'b': 0.136, 'K': 284512.0}],
'angles': [ {'ai': 'O1', 'aj': 'C2', 'ak': 'H3', 'func': 1, 'theta': 123.0, 'K': 292.88},
{'ai': 'O1', 'aj': 'C2', 'ak': 'C4', 'func': 1, 'theta': 120.4, 'K': 669.44},
{'ai': 'C2', 'aj': 'C4', 'ak': 'C5', 'func': 1, 'theta': 120.0, 'K': 711.28},
{'ai': 'C2', 'aj': 'C4', 'ak': 'O11', 'func': 1, 'theta': 119.26, 'K': 567.94},
{'ai': 'H3', 'aj': 'C2', 'ak': 'C4', 'func': 1, 'theta': 115.0, 'K': 292.88},
{'ai': 'C4', 'aj': 'C5', 'ak': 'H6', 'func': 1, 'theta': 125.7, 'K': 292.88},
{'ai': 'C4', 'aj': 'C5', 'ak': 'C7', 'func': 1, 'theta': 107.3, 'K': 585.76},
{'ai': 'C4', 'aj': 'O11', 'ak': 'C9', 'func': 1, 'theta': 106.5, 'K': 585.76},
{'ai': 'C5', 'aj': 'C4', 'ak': 'O11', 'func': 1, 'theta': 110.6, 'K': 585.76},
{'ai': 'C5', 'aj': 'C7', 'ak': 'H8', 'func': 1, 'theta': 127.5, 'K': 292.88},
{'ai': 'C5', 'aj': 'C7', 'ak': 'C9', 'func': 1, 'theta': 107.3, 'K': 585.76},
{'ai': 'H6', 'aj': 'C5', 'ak': 'C7', 'func': 1, 'theta': 127.5, 'K': 292.88},
{'ai': 'C7', 'aj': 'C9', 'ak': 'H10', 'func': 1, 'theta': 132.1, 'K': 292.88},
{'ai': 'C7', 'aj': 'C9', 'ak': 'O11', 'func': 1, 'theta': 110.6, 'K': 585.76},
{'ai': 'H8', 'aj': 'C7', 'ak': 'C9', 'func': 1, 'theta': 125.7, 'K': 292.88},
{'ai': 'H10', 'aj': 'C9', 'ak': 'O11', 'func': 1, 'theta': 113.4, 'K': 292.88}],
'dihedrals': [ {'ai': 'O1', 'aj': 'C2', 'ak': 'C4', 'al': 'C5', 'func': 1, 'phiS': 180.0, 'K': 12.029, 'n': 2},
{'ai': 'O1', 'aj': 'C2', 'ak': 'C4', 'al': 'O11', 'func': 1, 'phiS': 180.0, 'K': 12.029, 'n': 2},
{'ai': 'C2', 'aj': 'C4', 'ak': 'C5', 'al': 'H6', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C2', 'aj': 'C4', 'ak': 'C5', 'al': 'C7', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C2', 'aj': 'C4', 'ak': 'O11', 'al': 'C9', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'H3', 'aj': 'C2', 'ak': 'C4', 'al': 'C5', 'func': 1, 'phiS': 180.0, 'K': 12.029, 'n': 2},
{'ai': 'H3', 'aj': 'C2', 'ak': 'C4', 'al': 'O11', 'func': 1, 'phiS': 180.0, 'K': 12.029, 'n': 2},
{'ai': 'C4', 'aj': 'C5', 'ak': 'C7', 'al': 'H8', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C4', 'aj': 'C5', 'ak': 'C7', 'al': 'C9', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C4', 'aj': 'O11', 'ak': 'C9', 'al': 'C7', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C4', 'aj': 'O11', 'ak': 'C9', 'al': 'H10', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C5', 'aj': 'C4', 'ak': 'O11', 'al': 'C9', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C5', 'aj': 'C7', 'ak': 'C9', 'al': 'H10', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C5', 'aj': 'C7', 'ak': 'C9', 'al': 'O11', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'H6', 'aj': 'C5', 'ak': 'C4', 'al': 'O11', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'H6', 'aj': 'C5', 'ak': 'C7', 'al': 'H8', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'H6', 'aj': 'C5', 'ak': 'C7', 'al': 'C9', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'C7', 'aj': 'C5', 'ak': 'C4', 'al': 'O11', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'H8', 'aj': 'C7', 'ak': 'C9', 'al': 'H10', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'H8', 'aj': 'C7', 'ak': 'C9', 'al': 'O11', 'func': 1, 'phiS': 180.0, 'K': 15.167, 'n': 2},
{'ai': 'O1', 'aj': 'C2', 'ak': 'H3', 'al': 'C4', 'func': 4, 'phiS': 180.0, 'K': 1.1, 'n': 2},
{'ai': 'C2', 'aj': 'C5', 'ak': 'C4', 'al': 'O11', 'func': 4, 'phiS': 180.0, 'K': 4.6024, 'n': 2},
{'ai': 'C4', 'aj': 'C7', 'ak': 'C5', 'al': 'H6', 'func': 4, 'phiS': 180.0, 'K': 4.6024, 'n': 2},
{'ai': 'C7', 'aj': 'H10', 'ak': 'C9', 'al': 'O11', 'func': 4, 'phiS': 180.0, 'K': 4.6024, 'n': 2},
{'ai': 'C9', 'aj': 'C5', 'ak': 'C7', 'al': 'H8', 'func': 4, 'phiS': 180.0, 'K': 4.6024, 'n': 2}],
}
[docs]
class waterFF_SPC(forcefield):
_parameter = {
'name': 'water SPC force field',
'charges': {
'OW': -2*0.41,
'HW1': 0.41,
'HW2': 0.41,
},
'masses': {
'OW': cv.periodicTable['O']['atomic_mass'],
'HW1': cv.periodicTable['H']['atomic_mass'],
'HW2': cv.periodicTable['H']['atomic_mass'],
},
# sigma=3.15061e-01 6.36386e-01 (amber)
# sigma 3.16557e-01 eps 6.50194e-01 (opls, charmm = literature)
# 0.0026173456 2.634129e-06 (gromos)
'C6': {'OW': 0.002618417, }, # unit kJ/mol nm^6
'C12': {'OW': 2.636964446090433e-06}, # unit kJ/mol nm^12
'bonds': [],
'angles': [],
'dihedrals': [],
'constraint_dOH': 0.1,
'constraint_angle': 109.47,
'bond_constraints': [
{'ai': 'OW', 'aj': 'HW1', 'value': 0.1},
{'ai': 'OW', 'aj': 'HW2', 'value': 0.1},
],
'angle_constraints': [
{'ai': 'HW1', 'aj': 'OW', 'ak': 'HW2',
'value': 109.47},
],
}
[docs]
class waterFF_SPCE(forcefield):
_parameter = {
'name': 'water SPC/E force field J. Phys. Chem. 91, 6269-6271 (1981)',
'charges': {
'OW': -2*0.4238,
'HW1': 0.4238,
'HW2': 0.4238,
},
'masses': {
'OW': cv.periodicTable['O']['atomic_mass'],
'HW1': cv.periodicTable['H']['atomic_mass'],
'HW2': cv.periodicTable['H']['atomic_mass'],
},
'C6': {'OW': 0.002618417, }, # unit kJ/mol nm^6
'C12': {'OW': 2.636964446090433e-06}, # unit kJ/mol nm^12
'bonds': [],
'angles': [],
'dihedrals': [],
'constraint_dOH': 0.1,
'constraint_angle': 109.47,
'bond_constraints': [
{'ai': 'OW', 'aj': 'HW1', 'value': 0.1},
{'ai': 'OW', 'aj': 'HW2', 'value': 0.1},
# {'ai': 'HW1', 'aj': 'HW2', 'value': 0.16330},
],
'angle_constraints': [
{'ai': 'HW1', 'aj': 'OW', 'ak': 'HW2',
'value': 109.47},
],
}
[docs]
class waterFF_TIP4P(forcefield):
_parameter = {
'name': 'water TIP4P force field. J. Chem. Phys. 79, 926 ',
'charges': {
'OW': -2*0.52,
'HW1': 0.52,
'HW2': 0.52,
},
'masses': {
'OW': cv.periodicTable['O']['atomic_mass'],
'HW1': cv.periodicTable['H']['atomic_mass'],
'HW2': cv.periodicTable['H']['atomic_mass'],
},
# sigma = 3.15365e-01 nm
# eps = 78.0 K * k = 0.648527.8 kJ/mol
'C6': {'OW':
_EpsSigmaToC6C12(78.0 * cv.FPC_K_B_EV * cv.ev2au,
3.15365e-01 * 10 * cv.an2au)[0]
}, # unit kJ/mol nm^6
'C12': {'OW':
_EpsSigmaToC6C12(78.0 * cv.FPC_K_B_EV * cv.ev2au,
3.15365e-01 * 10 * cv.an2au)[1]
}, # unit kJ/mol nm^12
'bonds': [],
'angles': [],
'dihedrals': [],
'bond_constraints': [
{'ai': 'OW', 'aj': 'HW1', 'value': 0.09572},
{'ai': 'OW', 'aj': 'HW2', 'value': 0.09572},
],
'angle_constraints': [
{'ai': 'HW1', 'aj': 'OW', 'ak': 'HW2',
'value': 104.52},
],
'drelom': 0.13194,
}
[docs]
def hasVirtualSite(self):
return True
[docs]
def getVGeom(self, geom, noM=False):
"""
Input:
geom: geometry as 2 dim numpy array geom[number of atoms, 3]
noM: boolean, whether virtual sites for this molecule should be skipped
Output: vGeom, gM
vGeom: geometry with virtual sites as 2 dim numpy array geom[number of atoms, 3]
gM: Jacobi matrix elements as list ot tuples (ai,i,bj,j,g)
ai: index of real atom
i: index of atom coord (0...2)
bj: index of virtual site
j: index of virtual site coord (0...2)
g: derviative dj / di
For the TIP4P water molecule, the oxygen geometry is swapped with the geometry
of the virtual site. Hydrogen positions in vGeom are identical to geom
"""
if noM:
gM = [(ai, i, ai, i, 1.0) for i in range(3)
for ai in range(geom.shape[0])]
return geom.copy(), gM
vGeom = geom.copy()
m, gM = self._getM(geom)
vGeom[0, :] = m
gM = [(ai, i, 0, j, g)
for (ai, i, j), g in np.ndenumerate(gM.reshape(3, 3, 3)) if g != 0.]
gM += [(1, i, 1, i, 1.0) for i in range(3)]
gM += [(2, i, 2, i, 1.0) for i in range(3)]
return vGeom, gM
def _getM(self, geom):
"""
Obtain virtual site M as defined in the TIP4P/2005f model. Furthermore, the derivative of the
virtual site coordinates with respect to the cartesian coordinates of the associated water molecule
are calculated.
Input:
geom: The cartesian coordinates of a water molecule. 2d numpy array
Ordering: geom[0,:] =Oxygen (x, y, z),
geom[1,:] =Hydrogen (x, y, z),
geom[2,:] =Hydrogen (x, y, z)
Units: Bohr
Output:
m: Virtual site coordinates M; Units: Bohr
g: Derivative of the virtual site coordinates with respect to the cartesian coordinates of
the associated water molecule as 2D numpy array with dimensions (9,3).
"""
drelom = self._parameter['drelom']
rOH1 = geom[1, :] - geom[0, :]
rOH2 = geom[2, :] - geom[0, :]
dOH1 = np.linalg.norm(rOH1)
dOH2 = np.linalg.norm(rOH2)
theta, dTheta_dRH1, dTheta_dRO, dTheta_dRH2 = _angleDER(
geom[1, :], geom[0, :], geom[2, :])
dDOH1_dRO = -rOH1 / dOH1
dDOH2_dRO = -rOH2 / dOH2
dDOH1_dRH1 = rOH1 / dOH1
dDOH2_dRH2 = rOH2 / dOH2
sinHalfTheta = np.sin(0.5*theta)
cosHalfTheta = np.cos(0.5*theta)
dCosHalfTheta_dRO = - 0.5 * sinHalfTheta * dTheta_dRO
dCosHalfTheta_dRH1 = - 0.5 * sinHalfTheta * dTheta_dRH1
dCosHalfTheta_dRH2 = - 0.5 * sinHalfTheta * dTheta_dRH2
dR1R2 = np.linalg.norm(rOH1+rOH2)
z = (rOH1+rOH2)/dR1R2
dZ_dRO = - 2. * np.eye(3, dtype=float) / dR1R2 +\
2./(dR1R2) * np.outer(z, z)
dZ_dRH1 = 1.0 * np.eye(3, dtype=float) / dR1R2 +\
-1./(dR1R2) * np.outer(z, z)
dZ_dRH2 = 1.0 * np.eye(3, dtype=float) / dR1R2 +\
-1./(dR1R2) * np.outer(z, z)
m = geom[0, :] + \
drelom * (dOH1+dOH2) * np.cos(0.5*theta) * z
dM_dRO = np.eye(3, dtype=float)
dM_dRO += drelom * cosHalfTheta * np.outer(z, dDOH1_dRO + dDOH2_dRO)
dM_dRO += drelom * (dOH1+dOH2) * np.outer(z, dCosHalfTheta_dRO)
dM_dRO += drelom * (dOH1+dOH2) * cosHalfTheta * dZ_dRO
dM_dRH1 = drelom * cosHalfTheta * np.outer(z, dDOH1_dRH1)
dM_dRH1 += drelom * (dOH1+dOH2) * np.outer(z, dCosHalfTheta_dRH1)
dM_dRH1 += drelom * (dOH1+dOH2) * cosHalfTheta * dZ_dRH1
dM_dRH2 = drelom * cosHalfTheta * np.outer(z, dDOH2_dRH2)
dM_dRH2 += drelom * (dOH1+dOH2) * np.outer(z, dCosHalfTheta_dRH2)
dM_dRH2 += drelom * (dOH1+dOH2) * cosHalfTheta * dZ_dRH2
return m, np.concatenate([dM_dRO.T, dM_dRH1.T, dM_dRH2.T])
[docs]
class waterFF_TIP4P2005(waterFF_TIP4P):
_parameter = {
'name': 'water TIP4P2005 force field. J. Chem. Phys. 123, 234505 ',
'charges': {
'OW': -2*0.5564,
'HW1': 0.5564,
'HW2': 0.5564,
},
'masses': {
'OW': cv.periodicTable['O']['atomic_mass'],
'HW1': cv.periodicTable['H']['atomic_mass'],
'HW2': cv.periodicTable['H']['atomic_mass'],
},
# sigma = 3.1589e-01 nm
# eps = 93.2 K
'C6': {'OW':
_EpsSigmaToC6C12(93.2 * cv.FPC_K_B_EV * cv.ev2au,
3.1589e-01 * 10 * cv.an2au)[0]
}, # unit kJ/mol nm^6
'C12': {'OW':
_EpsSigmaToC6C12(93.2 * cv.FPC_K_B_EV * cv.ev2au,
3.1589e-01 * 10 * cv.an2au)[1]
}, # unit kJ/mol nm^12
'bonds': [],
'angles': [],
'dihedrals': [],
'bond_constraints': [
{'ai': 'OW', 'aj': 'HW1', 'value': 0.09572},
{'ai': 'OW', 'aj': 'HW2', 'value': 0.09572},
],
'angle_constraints': [
{'ai': 'HW1', 'aj': 'OW', 'ak': 'HW2',
'value': 104.52},
],
'drelom': 0.13194,
}
[docs]
class waterFF_TIP4P2005f(waterFF_TIP4P):
_parameter = {
'name': 'water TIP4P2005f force field. J. Chem. Phys 135, 224516 (2011) ',
'charges': {
'OW': -1.1128,
'HW1': 0.5564,
'HW2': 0.5564,
},
'masses': {
'OW': cv.periodicTable['O']['atomic_mass'],
'HW1': cv.periodicTable['H']['atomic_mass'],
'HW2': cv.periodicTable['H']['atomic_mass'],
},
# sigma = 3.1644e-01 nm
# eps = 93.2 K
'C6': {'OW':
_EpsSigmaToC6C12(93.2 * cv.FPC_K_B_EV * cv.ev2au,
3.1644e-01 * 10 * cv.an2au)[0]
}, # unit kJ/mol nm^6
'C12': {'OW':
_EpsSigmaToC6C12(93.2 * cv.FPC_K_B_EV * cv.ev2au,
3.1644e-01 * 10 * cv.an2au)[1]
}, # unit kJ/mol nm^12
'bonds': [
{'ai': 'OW', 'aj': 'HW1',
'func': 3, 'b': 0.09419, 'D': 432.581, 'beta': 22.87},
{'ai': 'OW', 'aj': 'HW2',
'func': 3, 'b': 0.09419, 'D': 432.581, 'beta': 22.87},
],
'angles': [
{'ai': 'HW1', 'aj': 'OW', 'ak': 'HW2',
'func': 2, 'theta': 107.4, 'K': 367.810}],
'dihedrals': [],
'drelom': 0.13194,
}
[docs]
class jointForceFields(object):
# TODO: what forcefield should be used for which molecule (and potential parameters)
# should be read in from the input file
molecule2FF = {
"SOL": waterFF_TIP4P2005f,
"UREA": ureaFF,
}
def __init__(self, topol, LJCombRule="OPLS", waterForcefield=None, molecule2ffstr=None, **opts):
"""
initializes the jointForceFields object:
prepares for calculation of energies and gradients,
and initialized force field for each molecule individually
Input:
topol: list of strings where for each atom the molecule and the atom type is defined
this list must have the same order as the geometry
the strings consists of a molecule name and an atom type within the molecule
separated by '_',
E.g. "SOL1_OW" for solvent molecule 1, oxygen atom in water
LJCombRule: defines how LJ parameters are combined.
currently implemented are 'Lorentz-Berthelot' and 'OPLS'
default is 'OPLS'
"""
# override default forcefield for water
if waterForcefield is not None:
if waterForcefield == 'SPC/E':
self.molecule2FF['SOL'] = waterFF_SPCE
elif waterForcefield == 'SPC':
self.molecule2FF['SOL'] = waterFF_SPC
elif waterForcefield == 'TIP4P/2005f':
self.molecule2FF['SOL'] = waterFF_TIP4P2005f
elif waterForcefield == 'TP4P/2005':
self.molecule2FF['SOL'] = waterFF_TIP4P2005
elif waterForcefield == 'TIP4P':
self.molecule2FF['SOL'] = waterFF_TIP4P
else:
raise ValueError(
'Water force field \"%s\" not implemented' % waterForcefield)
# override general forcefield-molecule association
if molecule2ffstr is not None:
for key, value in json.loads(molecule2ffstr).items():
self.molecule2FF[key] = globals()[value]
self.topol = np.array(topol)
self.molecules = set([s.split("_")[0] for s in self.topol])
self.topolMs = {m: [s for s in topol if s.split(
'_')[0] == m] for m in self.molecules}
self.atomsOfMs = {
m: [i for i, t in enumerate(topol) if t.split('_')[0] == m]
for m in self.molecules}
self.mOfAtom = [s.split("_")[0] for s in self.topol]
self.moleculeType = [re.sub(r'[0-9]+', '', m) for m in self.molecules]
# the partial charges of all atoms
self.masses = np.zeros(self.topol.shape[0])
# the partial charges of all atoms
self.charges = np.zeros(self.topol.shape[0])
# the nuclear spread of all atoms
self.sigma_nuc = np.zeros(self.topol.shape[0])
# Define logical variable for the different active interactions
self.bonded = bool(opts.get('bonded', True))
self.coulomb = bool(opts.get('coulomb', True))
self.lennard_jones = bool(opts.get('lennard_jones', True))
self.ff = dict()
# initialize the force fields for each molecule, read out LJ parameters
eps = np.zeros((self.topol.shape[0]))
sigma = np.zeros((self.topol.shape[0]))
for m, t in zip(self.molecules, self.moleculeType):
self.ff[m] = self.molecule2FF[t](topol=self.topolMs[m])
self.masses[self.atomsOfMs[m]] = self.ff[m].getMasses()
self.charges[self.atomsOfMs[m]] = self.ff[m].getCharges()
sigma[self.atomsOfMs[m]] = self.ff[m].getLJsigma()
eps[self.atomsOfMs[m]] = self.ff[m].getLJeps()
if self.ff[m].getLJCombRule() is not None:
if self.ff[m].getLJCombRule() != LJCombRule:
raise ValueError(
f"The force field has combination rule {self.ff[m].getLJCombRule()} and you requested combination rule {LJCombRule}")
# Get nuclear spread for each atom if applicable
self.sigma_case = opts.get('sigma_case', None)
if self.sigma_case is None:
self.sigma_nuc = None
else:
for m in self.molecules:
self.sigma_nuc[self.atomsOfMs[m]] = self.ff[m].getSigmaNuc(self.sigma_case)
self.nonBondedPairs = np.array(
[(i, j) for i in range(self.charges.shape[0])
for j in range(i+1, self.charges.shape[0])
if self.mOfAtom[i] != self.mOfAtom[j]
])
fudgeFactorLJ = np.ones(self.nonBondedPairs.shape[0])
fudgeFactorQQ = np.ones(self.nonBondedPairs.shape[0])
# nonbonded pairs within a molecule:
if len(self.molecules)>0:
nonBondedPairsInMol = np.array([(self.atomsOfMs[m][ai], self.atomsOfMs[m][aj])
for ai, aj in self.ff[m].getNonBondedInteractionPairs() for m in self.molecules])
fudgeFactorLJinMol = np.array(
[fLJ for fLJ in self.ff[m].getNonBondedInteractionPairsFudgeLJ() for m in self.molecules])
fudgeFactorQQinMol = np.array(
[fQQ for fQQ in self.ff[m].getNonBondedInteractionPairsFudgeQQ() for m in self.molecules])
else:
nonBondedPairsInMol = np.array([])
fudgeFactorLJinMol = np.array([])
fudgeFactorQQinMol = np.array([])
# concatenate both parts
if self.nonBondedPairs.shape[0] == 0:
self.nonBondedPairs = nonBondedPairsInMol
elif nonBondedPairsInMol.shape[0] > 0:
self.nonBondedPairs = np.concatenate(
[self.nonBondedPairs, nonBondedPairsInMol])
self.fudgeFactorLJ = np.concatenate(
[fudgeFactorLJ, fudgeFactorLJinMol])
self.fudgeFactorQQ = np.concatenate(
[fudgeFactorQQ, fudgeFactorQQinMol])
if self.nonBondedPairs.shape[0] > 0:
self._nonBondedPairInd1 = [
np.where(self.nonBondedPairs[:, 0] == i)[0]
for i in range(self.charges.shape[0])]
self._nonBondedPairInd2 = [
np.where(self.nonBondedPairs[:, 1] == i)[0]
for i in range(self.charges.shape[0])]
else:
self._nonBondedPairInd1 = []
self._nonBondedPairInd2 = []
self.hasVirtualSites = np.any(
[self.ff[m].hasVirtualSite() for m in self.molecules])
# LJ epsilon for all atom pairs
self.LJeps = np.zeros(self.nonBondedPairs.shape[0])
# LJ sigma for all atom pairs
self.LJsigma = np.zeros(self.nonBondedPairs.shape[0])
# combine the LJ parameters according to LJCombRules
for i, (a1, a2) in enumerate(self.nonBondedPairs):
if LJCombRule == "OPLS":
self.LJsigma[i] = np.sqrt(sigma[a1]*sigma[a2])
self.LJeps[i] = np.sqrt(eps[a1]*eps[a2])
elif LJCombRule == "Lorentz-Berthelot":
self.LJsigma[i] = (sigma[a1] + sigma[a2])/2.
self.LJeps[i] = np.sqrt(eps[a1]*eps[a2])
elif LJCombRule == "geometric":
C6_1, C12_1 = _EpsSigmaToC6C12(eps[a1], sigma[a1])
C6_2, C12_2 = _EpsSigmaToC6C12(eps[a2], sigma[a2])
C6 = np.sqrt(C6_1 * C6_2)
C12 = np.sqrt(C12_1 * C12_2)
self.LJsigma[i] = np.power(C12/C6, 1./6)
self.LJeps[i] = 1/4. * C6 * C6 / C12
else:
raise ValueError("%s not implemented" % LJCombRule)
# Check if any forcefield has constraints
self._numConstraints = 0
for m in sorted(self.molecules):
if self.ff[m].hasConstraints():
self._numConstraints = self._numConstraints + \
self.ff[m].numConstraints()
[docs]
def enforceConstraints(self, x):
x = x.reshape((-1, 3))
for m in sorted(self.molecules):
if self.ff[m].hasConstraints():
constraints = cnstr.Constraints(self.ff[m].getConstraints)
x1 = constraints.enforceConstraints(x[self.atomsOfMs[m],:], self.masses[self.atomsOfMs[m]])
x[self.atomsOfMs[m],:] = x1
return x.flatten()
[docs]
def getRattleCorrection(self, vec0, vec1, dt, case="pos"):
"""
This method returns the correction to the position and velocity
of a constrained geometry needed for the RATTLE algorithm.
If case = "pos", vec0 and vec1 represent r(t) and r(t+dt) without
constraints, respectively. In this case, the correction to v(t+dt/2)
is returned as it is needed to calculate r(t+dt) with constraints.
If case = "vel" vec0 and vec1 represent r(t+dt) with constraints
and v(t+dt) without constraints. In this case the correction to
v(t+dt) with constraints is returned directly.
Parameters
----------
vec0 : (3 * n_atoms) ndarray of float.
r(t) ["pos"] or r(t+dt) ["vel"] without constraints applied.
vec1 : (3 * n_atoms) ndarray of float.
r(t+dt) ["pos"] or v(t+dt) ["vel"] with and without constraints
applied, respectively.
dt : float
Timestep
case : str, optional
Flag to calculate correction to position ["pos"] or velocity ["vel"],
by default "pos"
Returns
-------
correction: (n_atoms, 3) ndarray of float.
Correction to v(t+dt/2) ["pos"] or v(t+dt) ["vel"] for geometries
with constraints.
"""
# todo: Option for different topologies?
vec0 = vec0.reshape((-1, 3))
vec1 = vec1.reshape((-1, 3))
# Array with corrections to position according Rattle
correction = np.zeros(vec0.shape, dtype=float)
if self._numConstraints == 0:
warnings.warn("Correction is 0 if there are no constraints.")
return correction
if case != "pos" and case != "vel":
raise ImportError(
"Wrong case for RATTLE implementation, please check documentation.")
for m in sorted(self.molecules):
if self.ff[m].hasConstraints():
# ! maybe too complicated
vec0_mol = vec0[self.atomsOfMs[m], :]
vec1_mol = vec1[self.atomsOfMs[m], :]
mass = self.masses[self.atomsOfMs[m]]
constraints = cnstr.Constraints(self.ff[m].getConstraints)
if case == 'pos':
correction[self.atomsOfMs[m]] = constraints.get_correction_r(
vec0_mol, vec1_mol, dt, mass)
elif case == 'vel':
correction[self.atomsOfMs[m]] = constraints.get_correction_v(
vec0_mol, vec1_mol, dt, mass)
return correction.flatten()
[docs]
def getEnGrad(self, x, listNoM=[]):
"""
This subroutine calculates the energy and the gradient along all nuclear coordinates.
At the moment is able to calculate the energy and the gradient of the following
interactions:
* Bonded interactions
* Lennard-Jones interaction
* Coulomb interaction (deactivated for electrostatic embedding)
Parameters
----------
x : (3 * n_atoms) ndarray of float.
Coordinates of the nuclear centers.
listNoM : list[int], optional
Indices labeling the atoms for which virtual sites are disabled, by default []
Returns
-------
energy: float.
Energy given by the force field. Units are Hartree!
gradient: (3 * n_atoms) ndarray of float.
Gradient of the energy along all nuclear coordinates
as flat numpy array. The order is as in the input.
Units are Hartree/Bohr.
"""
geom = x.reshape((-1, 3))
gradient = np.zeros(geom.shape)
energy = 0.
eBonds = 0.
eAngles = 0.
eDihedrals = 0.
# Bonded interactions
if self.bonded:
for m in self.molecules:
geomOfM = geom[self.atomsOfMs[m], :]
e, g = self.ff[m].getBonded(
geomOfM.flatten(),
self.topolMs[m])
energy += e
eBonds += self.ff[m].eBonds
eAngles += self.ff[m].eAngles
eDihedrals += self.ff[m].eDihedrals
gradient[self.atomsOfMs[m], :] += g.reshape(-1, 3)
#print(f"ebond: {eBonds} eAngles:{eAngles} eDihedrals:{eDihedrals}")
if self.nonBondedPairs.shape[0] == 0:
# Return if no non-bonded interactions
return energy, gradient.flatten()
eNonBonded = 0
eCoul = 0
## Calculate non-bonded interaction energy and gradient ##
# !: Avoid recalculating distances between atoms as much as possible.
dist = geom[self.nonBondedPairs[:, 0]] - \
geom[self.nonBondedPairs[:, 1]]
d = np.linalg.norm(dist, axis=1)
distOverD = dist * (1/d)[:, np.newaxis]
# Return warning if any value of d is below 0.4 bohr (~0.2 angstrom)
if np.any(d < 0.4):
warnings.warn("Some distances are smaller than 0.2 angstrom.")
# Lennard-Jones interaction
if self.lennard_jones:
ener_tmp, grad_tmp = self.lennard_jones_egrad(
d, distOverD, gradient.shape, self.fudgeFactorLJ)
energy += ener_tmp
gradient += grad_tmp
eNonBonded += ener_tmp
# Coulomb interaction
if self.coulomb:
if self.hasVirtualSites: # calculation with virtual sites
vGeom, vGeomD = self.geomWithVirtualSites(geom, listNoM)
dist = vGeom[self.nonBondedPairs[:, 0]] - \
vGeom[self.nonBondedPairs[:, 1]]
d = np.linalg.norm(dist, axis=1)
distOverD = dist * (1/d)[:, np.newaxis]
ener_tmp, grad_tmp = self.coulomb_egrad(
d, distOverD, gradient.shape, self.fudgeFactorQQ)
grad_tmp = self.gradient_virt2orig(grad_tmp, vGeomD)
else:
ener_tmp, grad_tmp = self.coulomb_egrad(
d, distOverD, gradient.shape, self.fudgeFactorQQ)
energy += ener_tmp
gradient += grad_tmp
eCoul += ener_tmp
# Convert virtual geometry to the original coordinates (if necessary).
#print(f"EvdW: {eNonBonded} Ecoul {eCoul}")
return energy, gradient.flatten()
# TODO: Find a way to avoid parsing grad_shape (also for LJ)?
[docs]
def coulomb_egrad(self, dist_norm, dist_uvec, grad_shape, fudgeFactor):
"""
Calculate the Coulomb interaction energy and gradient.
Parameters
----------
dist_norm : (n_nonbonded) ndarray of floats.
Norm of the distance between all pairs of atoms.
dist_uvec : (n_nonbonded, 3)
Unit vector between all pairs of atoms.
grad_shape : tuple of int.
Shape of the gradient (n_nuclei, 3)
fudgeFactor : (n_nonbonded) ndarray of floats.
Scaling factor for each non-bonded pair
Returns
-------
energy_coul : float.
Coulomb interaction energy.
gradient_coul : (n_nuclei, 3) ndarray of float.
Gradient of the Coulomb interaction energy.
"""
q1q2 = self.charges[self.nonBondedPairs[:, 0]] * \
self.charges[self.nonBondedPairs[:, 1]]
q1q2OverD = q1q2 / dist_norm * fudgeFactor
q1q2OverD2 = -q1q2OverD / dist_norm
gradPairs = dist_uvec * q1q2OverD2[:, np.newaxis] * fudgeFactor[:, np.newaxis]
gradient = np.zeros(grad_shape)
for a1 in range(grad_shape[0]):
# ?: This += might not be necessary.
gradient[a1, :] += (gradPairs[self._nonBondedPairInd1[a1], :]).sum(axis=0)
gradient[a1, :] += - \
(gradPairs[self._nonBondedPairInd2[a1], :]).sum(axis=0)
return q1q2OverD.sum(), gradient
[docs]
def lennard_jones_egrad(self, dist_norm, dist_uvec, grad_shape, fudgeFactor):
"""
Calculates the energy and gradient of the Lennard-Jones interaction.
Parameters
----------
dist_norm : (n_nonbonded) ndarray of floats.
Norm of the distance between all pairs of atoms.
dist_uvec : (n_nonbonded, 3)
Unit vector between all pairs of atoms.
grad_shape : tuple of int.
Shape of the gradient (n_nuclei, 3)
fudgeFactor : (n_nonbonded) ndarray of floats.
Scaling factor for each non-bonded pair
Returns
-------
energy_LJ : float.
Lennard-Jones interaction energy.
gradient_LJ : (n_nuclei, 3) ndarray of float.
Gradient of the Lennard-Jones interaction energy.
"""
LJ12 = (self.LJsigma[:] / dist_norm)**12
LJ6 = (self.LJsigma[:] / dist_norm)**6
LJ = 4.0 * self.LJeps[:] * (LJ12[:] - LJ6[:]) * fudgeFactor[:]
dLJdd = 4.0 * \
self.LJeps[:] * (-12.0/dist_norm[:] * LJ12[:] +
6.0/dist_norm[:] * LJ6[:])
gradPairs = dist_uvec * dLJdd[:, np.newaxis] * fudgeFactor[:, np.newaxis]
gradient = np.zeros(grad_shape)
for a1 in range(grad_shape[0]):
gradient[a1,
:] += (gradPairs[self._nonBondedPairInd1[a1], :]).sum(axis=0)
gradient[a1, :] += - \
(gradPairs[self._nonBondedPairInd2[a1], :]).sum(axis=0)
return LJ.sum(), gradient
[docs]
def geomWithVirtualSites(self, geom, listNoM=[]):
"""
Function translates geometry into geometry with virtual sites that can be used for
Coulomb calculation.
Parameters
----------
geom : (n_atoms, 3) ndarray of floats.
Positions of each atomic nuclei.
listNoM : list, optional
List of atom indices for which virtual sites should be skipped, by default []
Returns
-------
vGeom : (n_atoms, 3) ndarray of floats.
Geometry with virtual sites.
vGeomD: list.
List of elements of the Jacobian matrix between real and virtual geometry
each element is a tuple (ai, i, aj, j, g), where:
* `ai`: real geometry atomic index.
* `i`: real geometry coordinate index 0...2.
* `aj`: virtual geometry site index.
* `j`: virtual geometry coordinate index 0...2.
* `g`: derivative d r_j / d r_i.
"""
vGeom = np.zeros(geom.shape)
vGeomD = []
for m in self.molecules:
if self.ff[m].hasVirtualSite():
vGeomMol, gM = self.ff[m].getVGeom(
geom[self.atomsOfMs[m]],
noM=any([i in listNoM for i in self.atomsOfMs[m]]))
vGeom[self.atomsOfMs[m]] = vGeomMol
vGeomD += [(self.atomsOfMs[m][ai], i, self.atomsOfMs[m][aj], j, float(g))
for (ai, i, aj, j, g) in gM]
else:
# if not virtual sites for this molecule
# copy geometry and add unit derivative values into vGeomD
vGeom[self.atomsOfMs[m]] = geom[self.atomsOfMs[m]]
vGeomD += [(ai, i, ai, i, 1.0)
for ai in self.atomsOfMs[m] for i in range(3)]
return vGeom, vGeomD
[docs]
def velWithVirtualSites(self, velocity, vGeomD):
"""
This function translates velocities of original coordinates to velocities
of virtual sites.
.. math::
\\frac{d M_k}{d t} = \\sum_{x_i} \\frac{\\partial M_k}{\\partial x_i} \\frac{d x_i}{d}
Parameters
----------
velocity : (n_atoms, 3) ndarray of floats.
Valocity of each virtual site position (some virtual positions are already real ones).
vGeomD : list
List of elements of the Jacobian matrix between real and virtual geometry
each element is a tuple (ai, i, aj, j, g), where:
* `ai`: real geometry atomic index.
* `i`: real geometry coordinate index 0...2.
* `aj`: virtual geometry site index.
* `j`: virtual geometry coordinate index 0...2.
* `g`: derivative d r_j / d r_i.
Returns
-------
vel_res : (n_atoms, 3) ndarray of floats.
Velocity of virtual sites.
"""
vel_res = np.zeros(velocity.shape)
for (ai, i, aj, j, g) in vGeomD:
vel_res[aj, j] += velocity[ai, i] * g
return vel_res
[docs]
def gradient_virt2orig(self, gradient, vGeomD):
"""
This function translates gradient of virtual sites into gradient of
original geometry.
.. math::
\\frac{\\partial V}{\\partial x_i} = \\sum_{M_j} \\frac{\\partial V}{\\partial M_j} \\frac{\\partial M_j}{\\partial x_i}
Parameters
----------
gradient : (n_atoms, 3) ndarray of floats.
Gradient at each virtual site position (some virtual positions are already real ones).
vGeomD : list
List of elements of the Jacobian matrix between real and virtual geometry
each element is a tuple (ai, i, aj, j, g), where:
* `ai`: real geometry atomic index.
* `i`: real geometry coordinate index 0...2.
* `aj`: virtual geometry site index.
* `j`: virtual geometry coordinate index 0...2.
* `g`: derivative d r_j / d r_i.
Returns
-------
grad_res : (n_atoms, 3) ndarray of floats.
Gradient at each atom `real` position.
"""
grad_res = np.zeros(gradient.shape)
for (ai, i, aj, j, g) in vGeomD:
grad_res[ai, i] += g * gradient[aj, j]
return grad_res