Skip to content

GitLab

Unverified Commit 4442c095 authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2598 from z-gong/psf 

bug fix for handling constraints and rigidWater in CharmmPsfFile
parents dd85d9e0 d47c51b9
Hide whitespace changes
Inline Side-by-side
Showing with 155 additions and 45 deletions
wrappers/python/simtk/openmm/app/charmmpsffile.py
View file @ 4442c095
......@@ -497,7 +497,7 @@ class CharmmPsfFile(object):
for a1 in a2.bond_partners:
for a4 in a3.bond_partners:
pair = (min(a1.idx, a4.idx), max(a1.idx, a4.idx),)
if a1 != a4:
if a1 != a3 and a2 != a4 and a1 != a4:
pair_14_set.add(pair)
# in case there are 3,4,5-member rings
......@@ -856,7 +856,7 @@ class CharmmPsfFile(object):
ewaldErrorTolerance : float=0.0005
The error tolerance to use if the nonbonded method is Ewald, PME, or LJPME.
flexibleConstraints : bool=True
Are our constraints flexible or not?
If True, parameters for constrained degrees of freedom will be added to the System
verbose : bool=False
Optionally prints out a running progress report
gbsaModel : str=None
......@@ -939,31 +939,29 @@ class CharmmPsfFile(object):
pass
# Set up the constraints
if verbose and (constraints is not None and not rigidWater):
def _is_bond_in_water(bond):
return bond.atom1.residue.resname in WATNAMES and \
tuple(sorted([bond.atom1.type.atomic_number, bond.atom2.type.atomic_number])) == (1, 8)
n_cons_bond = n_cons_angle = 0
if verbose and (constraints is not None or rigidWater):
print('Adding constraints...')
if constraints in (ff.HBonds, ff.AllBonds, ff.HAngles):
for bond in self.bond_list:
if (bond.atom1.type.atomic_number != 1 and
bond.atom2.type.atomic_number != 1):
continue
system.addConstraint(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv)
if constraints in (ff.AllBonds, ff.HAngles):
for bond in self.bond_list:
if (bond.atom1.type.atomic_number == 1 or
bond.atom2.type.atomic_number == 1):
continue
for bond in self.bond_list:
if constraints in (ff.AllBonds, ff.HAngles):
system.addConstraint(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv)
if rigidWater and constraints is None:
for bond in self.bond_list:
if (bond.atom1.type.atomic_number != 1 and
bond.atom2.type.atomic_number != 1):
continue
if (bond.atom1.residue.resname in WATNAMES and
bond.atom2.residue.resname in WATNAMES):
n_cons_bond += 1
elif constraints is ff.HBonds:
if bond.atom1.type.atomic_number == 1 or bond.atom2.type.atomic_number == 1:
system.addConstraint(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv)
n_cons_bond += 1
elif rigidWater:
if _is_bond_in_water(bond):
system.addConstraint(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv)
n_cons_bond += 1
# Add virtual sites
if hasattr(self, 'lonepair_list'):
......@@ -999,14 +997,13 @@ class CharmmPsfFile(object):
force = mm.HarmonicBondForce()
force.setForceGroup(self.BOND_FORCE_GROUP)
# See which, if any, energy terms we omit
omitall = not flexibleConstraints and constraints is ff.AllBonds
omith = omitall or (flexibleConstraints and constraints in
(ff.HBonds, ff.AllBonds, ff.HAngles))
omit_all = not flexibleConstraints and constraints in (ff.AllBonds, ff.HAngles)
omit_h = not flexibleConstraints and constraints is not None
omit_h_in_water = not flexibleConstraints and (constraints is not None or rigidWater)
for bond in self.bond_list:
if omitall: continue
if omith and (bond.atom1.type.atomic_number == 1 or
bond.atom2.type.atomic_number == 1):
continue
if omit_all: continue
if omit_h and (bond.atom1.type.atomic_number == 1 or bond.atom2.type.atomic_number == 1): continue
if omit_h_in_water and _is_bond_in_water(bond): continue
force.addBond(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv,
2*bond.bond_type.k*bond_frc_conv)
......@@ -1025,16 +1022,16 @@ class CharmmPsfFile(object):
atom_constraints[c[1]].append((c[0], dist))
for angle in self.angle_list:
# Only constrain angles including hydrogen here
if (angle.atom1.type.atomic_number != 1 and
angle.atom2.type.atomic_number != 1 and
angle.atom3.type.atomic_number != 1):
if (angle.atom1.type.atomic_number != 1 and angle.atom3.type.atomic_number != 1):
continue
a1 = angle.atom1.type.atomic_number
a2 = angle.atom2.type.atomic_number
a3 = angle.atom3.type.atomic_number
nh = int(a1==1) + int(a3==1)
if constraints is ff.HAngles:
a1 = angle.atom1.atomic_number
a2 = angle.atom2.atomic_number
a3 = angle.atom3.atomic_number
nh = int(a1==1) + int(a2==1) + int(a3==1)
constrained = (nh >= 2 or (nh == 1 and a2 == 8))
constrained = (nh == 2 or (nh == 1 and a2 == 8))
elif rigidWater:
constrained = (nh == 2 and a2 == 8 and angle.atom1.residue.resname in WATNAMES)
else:
constrained = False # no constraints
if constrained:
......@@ -1046,17 +1043,19 @@ class CharmmPsfFile(object):
l2 = bond.bond_type.req * length_conv
# Compute the distance between the atoms and add a constraint
length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*
cos(angle.angle_type.theteq))
system.addConstraint(bond.atom1.idx, bond.atom2.idx, length)
cos(angle.angle_type.theteq*pi/180))
system.addConstraint(angle.atom1.idx, angle.atom3.idx, length)
n_cons_angle += 1
if flexibleConstraints or not constrained:
force.addAngle(angle.atom1.idx, angle.atom2.idx,
angle.atom3.idx, angle.angle_type.theteq*pi/180,
2*angle.angle_type.k*angle_frc_conv)
if verbose and (constraints is not None or rigidWater):
print(' Number of bond constraints:', n_cons_bond)
print(' Number of angle constraints:', n_cons_angle)
for angle in self.angle_list:
# Already did the angles with hydrogen above. So skip those here
if (angle.atom1.type.atomic_number == 1 or
angle.atom2.type.atomic_number == 1 or
angle.atom3.type.atomic_number == 1):
if (angle.atom1.type.atomic_number == 1 or angle.atom3.type.atomic_number == 1):
continue
force.addAngle(angle.atom1.idx, angle.atom2.idx,
angle.atom3.idx, angle.angle_type.theteq*pi/180,
......@@ -1364,9 +1363,9 @@ class CharmmPsfFile(object):
print('Build exclusion list...')
self._build_exclusion_list()
if verbose:
print('\tNumber of 1-2 pairs: %i' % len(self.pair_12_list))
print('\tNumber of 1-3 pairs: %i' % len(self.pair_13_list))
print('\tNumber of 1-4 pairs: %i' % len(self.pair_14_list))
print(' Number of 1-2 pairs: %i' % len(self.pair_12_list))
print(' Number of 1-3 pairs: %i' % len(self.pair_13_list))
print(' Number of 1-4 pairs: %i' % len(self.pair_14_list))
# Add 1-4 interactions
sigma_scale = 2**(-1/6)
......
wrappers/python/tests/TestCharmmFiles.py
View file @ 4442c095
......@@ -361,6 +361,40 @@ class TestCharmmFiles(unittest.TestCase):
# Compare with value computed with NAMD.
self.assertAlmostEqual(energy, -2154.5539, delta=1e-3*abs(energy))
def test_Constraints(self):
"""Test that bond and angles constraints are correctly added into the system"""
psf = CharmmPsfFile('systems/water_methanol.psf')
params = CharmmParameterSet('systems/water_methanol.prm')
# the system is made of one water molecule and one methanol molecule
hBonds_water = [[0, 1], [1, 2]]
hAngles_water = [[0, 2]]
hBonds_methanol = [[3, 4], [3, 5], [3, 6], [7, 8]]
allBonds_methanol = hBonds_methanol + [[3, 7]]
hAngles_methanol = [[4, 5], [4, 6], [5, 6], [3, 8]]
system = psf.createSystem(params, constraints=None, rigidWater=False)
self.assertEqual(system.getNumConstraints(), 0)
system = psf.createSystem(params, constraints=None, rigidWater=True)
self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
sorted(hBonds_water + hAngles_water))
system = psf.createSystem(params, constraints=HBonds, rigidWater=False)
self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
sorted(hBonds_water + hBonds_methanol))
system = psf.createSystem(params, constraints=HBonds, rigidWater=True)
self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
sorted(hBonds_water + hAngles_water + hBonds_methanol))
system = psf.createSystem(params, constraints=AllBonds, rigidWater=False)
self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
sorted(hBonds_water + allBonds_methanol))
system = psf.createSystem(params, constraints=AllBonds, rigidWater=True)
self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
sorted(hBonds_water + hAngles_water + allBonds_methanol))
system = psf.createSystem(params, constraints=HAngles, rigidWater=False)
self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
sorted(hBonds_water + hAngles_water + allBonds_methanol + hAngles_methanol))
system = psf.createSystem(params, constraints=HAngles, rigidWater=True)
self.assertEqual(sorted(system.getConstraintParameters(i)[:2] for i in range(system.getNumConstraints())),
sorted(hBonds_water + hAngles_water + allBonds_methanol + hAngles_methanol))
if __name__ == '__main__':
unittest.main()
......
wrappers/python/tests/systems/water_methanol.prm 0 → 100644
View file @ 4442c095
* CHARMM FORCE FIELD GENERATED BY fftool
*
ATOMS
MASS 1 Hw 1.0080
MASS 2 Ow 15.9990
MASS 3 CTO 12.0110
MASS 4 H1O 1.0080
MASS 5 OH 15.9990
MASS 6 HO 1.0080
BONDS
Ow Hw 517.630258 1.000000
H1O CTO 339.997610 1.090000
OH CTO 320.004780 1.410000
HO OH 552.999522 0.945000
ANGLES
Hw Ow Hw 37.950526 109.470000
H1O CTO H1O 32.994742 107.800000
H1O CTO OH 35.002390 109.500000
CTO OH HO 54.995220 108.500000
DIHEDRALS
HO OH CTO H1O 0.225000 3 0
NONBONDED
!VLJ = Eps,i,j[(Rmin,i,j/ri,j)**12 - 2(Rmin,i,j/ri,j)**6]
!epsilon: kcal/mole, Eps,i,j = sqrt(eps,i * eps,j)
!Rmin/2: A, Rmin,i,j = Rmin/2,i + Rmin/2,j
!atom ignored epsilon Rmin/2 ignored eps,1-4 Rmin/2,1-4
Hw 0.000000 -0.000000 0.000000 0.000000 -0.000000 0.000000
Ow 0.000000 -0.155425 1.776577 0.000000 -0.077713 1.776577
CTO 0.000000 -0.065999 1.964309 0.000000 -0.033000 1.964309
H1O 0.000000 -0.030000 1.403078 0.000000 -0.015000 1.403078
OH 0.000000 -0.170000 1.751041 0.000000 -0.085000 1.751041
HO 0.000000 -0.000000 0.000000 0.000000 -0.000000 0.000000
END
wrappers/python/tests/systems/water_methanol.psf 0 → 100644
View file @ 4442c095
PSF
1 !NTITLE
REMARKS Created by fftool
9 !NATOM
1 S 1 SPCE H Hw 0.423800 1.0080 0 0.0000 0.0000
2 S 1 SPCE O Ow -0.847600 15.9990 0 0.0000 0.0000
3 S 1 SPCE H Hw 0.423800 1.0080 0 0.0000 0.0000
4 S 2 meth C CTO 0.145000 12.0110 0 0.0000 0.0000
5 S 2 meth H H1O 0.040000 1.0080 0 0.0000 0.0000
6 S 2 meth H H1O 0.040000 1.0080 0 0.0000 0.0000
7 S 2 meth H H1O 0.040000 1.0080 0 0.0000 0.0000
8 S 2 meth O OH -0.683000 15.9990 0 0.0000 0.0000
9 S 2 meth H HO 0.418000 1.0080 0 0.0000 0.0000
7 !NBOND: bonds
1 2 2 3 4 5 4 6
4 7 4 8 8 9
8 !NTHETA: angles
1 2 3 5 4 6 5 4 7
5 4 8 6 4 7 6 4 8
7 4 8 4 8 9
3 !NPHI: dihedrals
9 8 4 5 9 8 4 6
9 8 4 7
0 !NIMPHI: impropers
0 !NDON: donors
0 !NACC: acceptors
0 !NNB
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment