# ******************************************************************************
# pysimm.forcefield.gaff module
# ******************************************************************************
#
# ******************************************************************************
# License
# ******************************************************************************
# The MIT License (MIT)
#
# Copyright (c) 2016 Michael E. Fortunato, Coray M. Colina
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN
# THE SOFTWARE.
import os
from itertools import permutations, combinations
import gasteiger
from pysimm.system import Angle, Dihedral, Improper
from forcefield import Forcefield
[docs]class Gaff(Forcefield):
"""pysimm.forcefield.Gaff
Forcefield object with typing rules for Gaff model.
By default reads data file in forcefields subdirectory.
Attributes:
ff_name: gaff
pair_style: lj
ff_class: 1
"""
def __init__(self, db_file=None):
if not db_file and db_file is not False:
db_file = os.path.join(os.path.dirname(os.path.realpath(__file__)),
os.pardir, os.pardir, 'dat', 'forcefields', 'gaff.json')
Forcefield.__init__(self, db_file)
self.name = 'gaff'
self.pair_style = 'lj'
self.bond_style = 'harmonic'
self.angle_style = 'harmonic'
self.dihedral_style = 'fourier'
self.improper_style = 'cvff'
self.ff_class = '1'
[docs] def assign_ptypes(self, s):
"""pysimm.forcefield.Gaff.assign_ptypes
Gaff specific particle typing rules.
Requires :class:`~pysimm.system.System` object :class:`~pysimm.system.Particle` objects have bonds defined.
*** use System.add_particle_bonding() to ensure this ***
*** Not entirely inclusive - some atom types not used ***
Args:
s: :class:`~pysimm.system.System`
Returns:
None
"""
all_types = set()
s.pair_style = self.pair_style
s.add_particle_bonding()
for p in s.particles:
p.bond_orders = [x.order for x in p.bonds]
if None in p.bond_orders:
error_print('error: bond orders are not set')
p.bond_elements = [x.a.elem if p is x.b else x.b.elem for x in
p.bonds]
p.nbonds = len(p.bond_elements)
if p.linker:
p.nbonds += 1
for p in s.particles:
if p.elem == 'H':
if 'O' in p.bond_elements:
water = False
for pb in p.bonded_to:
if pb.elem == 'O' and pb.bond_elements.count('H') == 2:
water = True
if water:
p.type_name = 'hw'
else:
p.type_name = 'ho'
elif 'N' in p.bond_elements:
p.type_name = 'hn'
elif 'P' in p.bond_elements:
p.type_name = 'hp'
elif 'S' in p.bond_elements:
p.type_name = 'hs'
elif 'C' in p.bond_elements:
for pb in p.bonded_to:
if pb.elem == 'C':
elctrwd = (pb.bond_elements.count('N') +
pb.bond_elements.count('O') +
pb.bond_elements.count('F') +
pb.bond_elements.count('P') +
pb.bond_elements.count('S') +
pb.bond_elements.count('Cl'))
if 4 in pb.bond_orders or 'A' in pb.bond_orders:
p.type_name = 'ha'
elif elctrwd == 0:
p.type_name = 'hc'
elif pb.nbonds == 4 and elctrwd == 1:
p.type_name = 'h1'
elif pb.nbonds == 4 and elctrwd == 2:
p.type_name = 'h2'
elif pb.nbonds == 4 and elctrwd == 3:
p.type_name = 'h3'
elif pb.nbonds == 3 and elctrwd == 1:
p.type_name = 'h4'
elif pb.nbonds == 3 and elctrwd == 2:
p.type_name = 'h5'
elif p.elem == 'C':
if p.nbonds == 3 and 'O' in p.bond_elements:
p.type_name = 'c'
elif 4 in p.bond_orders or 'A' in p.bond_orders:
for pb in p.bonded_to:
if pb.elem != 'C' and pb.elem != 'H':
p.type_name = 'cc'
if not p.type_name:
p.type_name = 'ca'
elif p.nbonds == 4:
p.type_name = 'c3'
elif p.nbonds == 3 and not 'O' in p.bond_elements:
p.type_name = 'c2'
elif p.nbonds == 2:
p.type_name = 'c1'
elif p.elem == 'F':
p.type_name = 'f'
elif p.elem == 'Cl':
p.type_name = 'cl'
elif p.elem == 'Br':
p.type_name = 'br'
elif p.elem == 'I':
p.type_name = 'i'
elif p.elem == 'N':
if 2 in p.bond_orders and p.nbonds == 2:
p.type_name = 'n2'
elif 2 in p.bond_orders and p.nbonds == 3:
p.type_name = 'na'
elif 3 in p.bond_orders:
p.type_name = 'n1'
elif p.bond_elements.count('O') == 2:
p.type_name = 'no'
elif p.nbonds <= 3:
amide = False
aromatic_ring = False
for pb in p.bonded_to:
if pb.elem == 'C':
if 4 in pb.bond_orders or 'A' in pb.bond_orders:
aromatic_ring = True
for b in pb.bonds:
bp = b.a if pb is b.b else b.b
if bp.elem == 'O' and b.order == 2:
amide = True
if amide:
p.type_name = 'n'
elif (4 in p.bond_orders or 'A' in p.bond_orders) and 'C' in p.bond_elements:
p.type_name = 'nc'
elif aromatic_ring:
p.type_name = 'nh'
else:
p.type_name = 'n3'
elif p.nbonds == 4:
p.type_name = 'n4'
else:
print(p.elem, p.nbonds, p.bond_elements, p.bond_orders)
elif p.elem == 'O':
if p.nbonds == 1:
p.type_name = 'o'
elif p.bond_elements.count('H') == 2:
p.type_name = 'ow'
elif p.bond_elements.count('H') == 1:
p.type_name = 'oh'
else:
p.type_name = 'os'
elif p.elem == 'P':
if 4 in p.bond_orders or 'A' in p.bond_orders:
p.type_name = 'pc'
elif p.nbonds == 2:
p.type_name = 'p2'
elif p.nbonds == 3 and p.bond_elements.count('H') == 3:
p.type_name = 'p3'
elif p.nbonds == 3:
p.type_name = 'p4'
elif p.nbonds == 4:
p.type_name = 'p5'
elif p.elem == 'S':
if p.nbonds == 1:
p.type_name = 's'
elif p.nbonds == 2 and 2 in p.bond_orders:
p.type_name = 's2'
elif p.nbonds == 3:
p.type_name = 's4'
elif p.nbonds == 4:
p.type_name = 's6'
elif len(set(p.bond_orders)) == 1 and p.bond_orders[0] == 1:
if 'H' in p.bond_elements:
p.type_name = 'sh'
elif p.nbonds == 2:
p.type_name = 'ss'
else:
print 'cant type particle %s' % p.tag
return p
type_ = self.particle_types.get(p.type_name)
if not type_:
print(p.tag, p.elem, p.type_name)
all_types.add(self.particle_types.get(p.type_name)[0])
for pt in all_types:
pt = pt.copy()
s.particle_types.add(pt)
for p in s.particles:
pt = s.particle_types.get(p.type_name)
if pt:
p.type = pt[0]
[docs] def assign_btypes(self, s):
"""pysimm.forcefield.Gaff.assign_btypes
Gaff specific bond typing rules.
Requires :class:`~pysimm.system.System` object :class:`~pysimm.system.Particle` objects have bonds, type and type.name defined.
*** use after assign_ptypes ***
Args:
s: :class:`~pysimm.system.System`
Returns:
None
"""
all_types = set()
s.bond_style = self.bond_style
for b in s.bonds:
bt = self.bond_types.get('%s,%s' % (b.a.type.name, b.b.type.name))
if bt:
b.type_name = bt[0].name
else:
print ('couldnt type this bond %s,%s'
% (b.a.type.name, b.b.type.name))
return b
all_types.add(self.bond_types.get(b.type_name)[0])
for bt in all_types:
bt = bt.copy()
s.bond_types.add(bt)
for b in s.bonds:
bt = s.bond_types.get(b.type_name)
if bt:
b.type = bt[0]
[docs] def assign_atypes(self, s):
"""pysimm.forcefield.Gaff.assign_atypes
Gaff specific boanglend typing rules.
Requires :class:`~pysimm.system.System` object :class:`~pysimm.system.Particle` objects have bonds, type and type.name defined.
*** use after assign_ptypes ***
Args:
s: :class:`~pysimm.system.System`
Returns:
None
"""
all_types = set()
s.angle_style = self.angle_style
s.add_particle_bonding()
for p in s.particles:
for p1 in p.bonded_to:
for p2 in p.bonded_to:
if p1 is not p2:
unique = True
for a in s.angles:
if ((a.a is p1 and a.b is p and a.c is p2) or
(a.a is p2 and a.b is p and a.c is p1)):
unique = False
if unique:
at = self.angle_types.get('%s,%s,%s'
% (p1.type.name,
p.type.name,
p2.type.name))
if at:
s.angles.add(Angle(type_name=at[0].name,
a=p1, b=p, c=p2))
all_types.add(at[0])
else:
print ('I cant type this angle %s,%s,%s'
% (p1.type.name,
p.type.name,
p2.type.name))
for at in all_types:
at = at.copy()
s.angle_types.add(at)
for a in s.angles:
at = s.angle_types.get(a.type_name)
if at:
a.type = at[0]
[docs] def assign_dtypes(self, s):
"""pysimm.forcefield.Gaff.assign_dtypes
Gaff specific dihedral typing rules.
Requires :class:`~pysimm.system.System` object :class:`~pysimm.system.Particle` objects have bonds, type and type.name defined.
*** use after assign_ptypes ***
Args:
s: :class:`~pysimm.system.System`
Returns:
None
"""
all_types = set()
s.dihedral_style = self.dihedral_style
for b in s.bonds:
for p1 in b.a.bonded_to:
for p2 in b.b.bonded_to:
if p1 is b.b or p2 is b.a:
continue
unique = True
for d in s.dihedrals:
if ((d.a == p1 and d.b == b.a and
d.c == b.b and d.d == p2) or
(d.a == p2 and d.b == b.b and
d.c == b.a and d.d == p1)):
unique = False
if unique:
p1_name = p1.type.name
a_name = b.a.type.name
b_name = b.b.type.name
p2_name = p2.type.name
dt = self.dihedral_types.get('%s,%s,%s,%s'
% (p1_name, a_name,
b_name, p2_name))
if dt:
if len(dt) == 1:
all_types.add(dt[0])
s.dihedrals.add(Dihedral(type_name=dt[0].name,
a=p1, b=b.a,
c=b.b, d=p2))
else:
index = 0
x = 5
for i in range(len(dt)):
if dt[i].name.count('X') < x:
index = i
x = dt[i].name.count('X')
dt = dt[index]
all_types.add(dt)
s.dihedrals.add(Dihedral(type_name=dt.name,
a=p1, b=b.a,
c=b.b, d=p2))
else:
print ('I cant type this dihedral %s,%s,%s,%s'
% (p1_name, a_name, b_name, p2_name))
for dt in all_types:
dt = dt.copy()
s.dihedral_types.add(dt)
for d in s.dihedrals:
dt = s.dihedral_types.get(d.type_name, item_wildcard=None)
if dt:
d.type = dt[0]
[docs] def assign_itypes(self, s):
"""pysimm.forcefield.Gaff.assign_itypes
Gaff specific improper typing rules.
There are none.
Args:
s: :class:`~pysimm.system.System`
Returns:
None
"""
all_types = set()
s.improper_style = self.improper_style
for p in s.particles:
if len(p.bonded_to) == 3:
for perm in permutations(p.bonded_to, 3):
p1_name = perm[0].type.eq_improper or perm[0].type.name
p2_name = perm[1].type.eq_improper or perm[1].type.name
p3_name = perm[2].type.eq_improper or perm[2].type.name
it = self.improper_types.get(','.join([p.type.name, p1_name,
p2_name, p3_name]), order=True)
if it:
all_types.add(it[0])
bonded_to = p.bonded_to.get('all')
s.impropers.add(Improper(type_name=it[0].name,
a=p, b=bonded_to[0],
c=bonded_to[1],
d=bonded_to[2]))
break
for it in all_types:
it = it.copy()
s.improper_types.add(it)
for i in s.impropers:
it = s.improper_types.get(i.type_name)
if it:
i.type = it[0]
[docs] def assign_charges(self, s, charges='gasteiger'):
"""pysimm.forcefield.Gaff.assign_charges
Charge assignment. Gasteiger is default for now.
Args:
s: :class:`~pysimm.system.System`
charges: gasteiger
Returns:
None
"""
if charges == 'gasteiger':
print('adding gasteiger charges')
gasteiger.set_charges(s)