# Source code for moldesign.forcefields.terms

# Copyright 2016 Autodesk Inc.
#
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#
# Unless required by applicable law or agreed to in writing, software
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and

import math

import numpy as np

from moldesign import units as u
from moldesign.forcefields.forcefield import FFTerm
from moldesign.geom import coords as geo

[docs]class HarmonicBondTerm(FFTerm):
def __init__(self, a1, a2, k, d0):
self.atoms = [a1, a2]
self.k = k
self.d0 = d0
self.bond = a1.bond_graph.get(a2, None)

def __str__(self):
return 'Harmonic bond ({atoms}), k={k.magnitude:6.3f}{k.units}, equil={d0.magnitude:6.2}{d0.units}'.format(
atoms=','.join(map(str, self.atoms)),
k=self.k,
d0=self.d0)

[docs]    def coord(self):
return geo.distance(*self.atoms)

[docs]    def energy(self):
return u.default.convert(self.k * (self.coord() - self.d0)**2)

[docs]class HarmonicAngleTerm(FFTerm):
def __init__(self, a1, a2, a3, k, theta0):
self.atoms = [a1, a2, a3]
self.k = k
self.theta0 = theta0
self.bonds = [a1.bond_graph.get(a2,None),
a2.bond_graph.get(a3,None)]

def __str__(self):
return 'Harmonic angle ({atoms}), k={k.magnitude:6.3f}{k.units}, equil={d0.magnitude:6.2}{d0.units}'.format(
atoms=','.join(map(str, self.atoms)),
k=self.k.defunits(),
d0=self.theta0.defunits())

[docs]    def coord(self):
return geo.angle(*self.atoms)

[docs]    def energy(self):
return u.default.convert(self.k * (self.coord() - self.theta0)**2)

[docs]class PeriodicTorsionTerm(FFTerm):
def __init__(self, a1, a2, a3, a4, n, v_n, gamma):
self.atoms = [a1, a2, a3, a4]
self.n = n
self.v_n = v_n
self.gamma = gamma
self.bonds = [a1.bond_graph.get(a2, None),
a2.bond_graph.get(a3, None),
a3.bond_graph.get(a4, None)]

def __str__(self):
return 'Periodic torsion ({atoms}),' \
' n={n}, k={k.magnitude:6.3f}{k.units}, equil={d0.magnitude:6.2}{d0.units}'.format(
atoms=','.join(map(str, self.atoms)),
k=self.v_n.defunits(),
n=self.n,
d0=self.gamma)

[docs]    def coord(self):
return geo.dihedral(*self.atoms)

[docs]    def energy(self):
return u.default.convert(0.5 * self.v_n * (1.0 + np.cos(self.n * self.coord() - self.gamma)))

[docs]class LennardJonesSigmaEps(FFTerm):
def __init__(self, atom, sigma, epsilon):
self.atom = atom
self.sigma = sigma
self.epsilon = epsilon

[docs]    def get_sigma(self, other):
return (self.sigma + other.sigma) / 2.0

[docs]    def get_epsilon(self, other):
return math.sqrt(self.sigma * other.sigma)