Skip to content

GitLab

Commit 0c646ace authored by Zheng Gong's avatar Zheng Gong
Browse files

bug fix for constraints and rigidWater in CharmmPsfFile

parent dd85d9e0
Show whitespace changes
Inline Side-by-side
Showing with 162 additions and 45 deletions
wrappers/python/simtk/openmm/app/charmmpsffile.py
View file @ 0c646ace
...@@ -497,7 +497,7 @@ class CharmmPsfFile(object): ...@@ -497,7 +497,7 @@ class CharmmPsfFile(object):
for a1 in a2.bond_partners: for a1 in a2.bond_partners:
for a4 in a3.bond_partners: for a4 in a3.bond_partners:
pair = (min(a1.idx, a4.idx), max(a1.idx, a4.idx),) 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) pair_14_set.add(pair)
# in case there are 3,4,5-member rings # in case there are 3,4,5-member rings
...@@ -856,7 +856,7 @@ class CharmmPsfFile(object): ...@@ -856,7 +856,7 @@ class CharmmPsfFile(object):
ewaldErrorTolerance : float=0.0005 ewaldErrorTolerance : float=0.0005
The error tolerance to use if the nonbonded method is Ewald, PME, or LJPME. The error tolerance to use if the nonbonded method is Ewald, PME, or LJPME.
flexibleConstraints : bool=True 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 verbose : bool=False
Optionally prints out a running progress report Optionally prints out a running progress report
gbsaModel : str=None gbsaModel : str=None
...@@ -939,31 +939,29 @@ class CharmmPsfFile(object): ...@@ -939,31 +939,29 @@ class CharmmPsfFile(object):
pass pass
# Set up the constraints # 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...') print('Adding constraints...')
if constraints in (ff.HBonds, ff.AllBonds, ff.HAngles):
for bond in self.bond_list: for bond in self.bond_list:
if (bond.atom1.type.atomic_number != 1 and if constraints in (ff.AllBonds, ff.HAngles):
bond.atom2.type.atomic_number != 1):
continue
system.addConstraint(bond.atom1.idx, bond.atom2.idx, system.addConstraint(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv) bond.bond_type.req*length_conv)
if constraints in (ff.AllBonds, ff.HAngles): n_cons_bond += 1
for bond in self.bond_list: elif constraints is ff.HBonds:
if (bond.atom1.type.atomic_number == 1 or if bond.atom1.type.atomic_number == 1 or bond.atom2.type.atomic_number == 1:
bond.atom2.type.atomic_number == 1):
continue
system.addConstraint(bond.atom1.idx, bond.atom2.idx, system.addConstraint(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv) bond.bond_type.req*length_conv)
if rigidWater and constraints is None: n_cons_bond += 1
for bond in self.bond_list: elif rigidWater:
if (bond.atom1.type.atomic_number != 1 and if _is_bond_in_water(bond):
bond.atom2.type.atomic_number != 1):
continue
if (bond.atom1.residue.resname in WATNAMES and
bond.atom2.residue.resname in WATNAMES):
system.addConstraint(bond.atom1.idx, bond.atom2.idx, system.addConstraint(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv) bond.bond_type.req*length_conv)
n_cons_bond += 1
# Add virtual sites # Add virtual sites
if hasattr(self, 'lonepair_list'): if hasattr(self, 'lonepair_list'):
...@@ -999,14 +997,13 @@ class CharmmPsfFile(object): ...@@ -999,14 +997,13 @@ class CharmmPsfFile(object):
force = mm.HarmonicBondForce() force = mm.HarmonicBondForce()
force.setForceGroup(self.BOND_FORCE_GROUP) force.setForceGroup(self.BOND_FORCE_GROUP)
# See which, if any, energy terms we omit # See which, if any, energy terms we omit
omitall = not flexibleConstraints and constraints is ff.AllBonds omit_all = not flexibleConstraints and constraints in (ff.AllBonds, ff.HAngles)
omith = omitall or (flexibleConstraints and constraints in omit_h = not flexibleConstraints and constraints is not None
(ff.HBonds, ff.AllBonds, ff.HAngles)) omit_h_in_water = not flexibleConstraints and (constraints is not None or rigidWater)
for bond in self.bond_list: for bond in self.bond_list:
if omitall: continue if omit_all: continue
if omith and (bond.atom1.type.atomic_number == 1 or if omit_h and (bond.atom1.type.atomic_number == 1 or bond.atom2.type.atomic_number == 1): continue
bond.atom2.type.atomic_number == 1): if omit_h_in_water and _is_bond_in_water(bond): continue
continue
force.addBond(bond.atom1.idx, bond.atom2.idx, force.addBond(bond.atom1.idx, bond.atom2.idx,
bond.bond_type.req*length_conv, bond.bond_type.req*length_conv,
2*bond.bond_type.k*bond_frc_conv) 2*bond.bond_type.k*bond_frc_conv)
...@@ -1025,16 +1022,16 @@ class CharmmPsfFile(object): ...@@ -1025,16 +1022,16 @@ class CharmmPsfFile(object):
atom_constraints[c[1]].append((c[0], dist)) atom_constraints[c[1]].append((c[0], dist))
for angle in self.angle_list: for angle in self.angle_list:
# Only constrain angles including hydrogen here # Only constrain angles including hydrogen here
if (angle.atom1.type.atomic_number != 1 and if (angle.atom1.type.atomic_number != 1 and angle.atom3.type.atomic_number != 1):
angle.atom2.type.atomic_number != 1 and
angle.atom3.type.atomic_number != 1):
continue 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: if constraints is ff.HAngles:
a1 = angle.atom1.atomic_number constrained = (nh == 2 or (nh == 1 and a2 == 8))
a2 = angle.atom2.atomic_number elif rigidWater:
a3 = angle.atom3.atomic_number constrained = (nh == 2 and a2 == 8 and angle.atom1.residue.resname in WATNAMES)
nh = int(a1==1) + int(a2==1) + int(a3==1)
constrained = (nh >= 2 or (nh == 1 and a2 == 8))
else: else:
constrained = False # no constraints constrained = False # no constraints
if constrained: if constrained:
...@@ -1046,17 +1043,19 @@ class CharmmPsfFile(object): ...@@ -1046,17 +1043,19 @@ class CharmmPsfFile(object):
l2 = bond.bond_type.req * length_conv l2 = bond.bond_type.req * length_conv
# Compute the distance between the atoms and add a constraint # Compute the distance between the atoms and add a constraint
length = sqrt(l1*l1 + l2*l2 - 2*l1*l2* length = sqrt(l1*l1 + l2*l2 - 2*l1*l2*
cos(angle.angle_type.theteq)) cos(angle.angle_type.theteq*pi/180))
system.addConstraint(bond.atom1.idx, bond.atom2.idx, length) system.addConstraint(angle.atom1.idx, angle.atom3.idx, length)
n_cons_angle += 1
if flexibleConstraints or not constrained: if flexibleConstraints or not constrained:
force.addAngle(angle.atom1.idx, angle.atom2.idx, force.addAngle(angle.atom1.idx, angle.atom2.idx,
angle.atom3.idx, angle.angle_type.theteq*pi/180, angle.atom3.idx, angle.angle_type.theteq*pi/180,
2*angle.angle_type.k*angle_frc_conv) 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: for angle in self.angle_list:
# Already did the angles with hydrogen above. So skip those here # Already did the angles with hydrogen above. So skip those here
if (angle.atom1.type.atomic_number == 1 or if (angle.atom1.type.atomic_number == 1 or angle.atom3.type.atomic_number == 1):
angle.atom2.type.atomic_number == 1 or
angle.atom3.type.atomic_number == 1):
continue continue
force.addAngle(angle.atom1.idx, angle.atom2.idx, force.addAngle(angle.atom1.idx, angle.atom2.idx,
angle.atom3.idx, angle.angle_type.theteq*pi/180, angle.atom3.idx, angle.angle_type.theteq*pi/180,
...@@ -1364,9 +1363,9 @@ class CharmmPsfFile(object): ...@@ -1364,9 +1363,9 @@ class CharmmPsfFile(object):
print('Build exclusion list...') print('Build exclusion list...')
self._build_exclusion_list() self._build_exclusion_list()
if verbose: if verbose:
print('\tNumber of 1-2 pairs: %i' % len(self.pair_12_list)) print(' Number of 1-2 pairs: %i' % len(self.pair_12_list))
print('\tNumber of 1-3 pairs: %i' % len(self.pair_13_list)) print(' Number 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-4 pairs: %i' % len(self.pair_14_list))
# Add 1-4 interactions # Add 1-4 interactions
sigma_scale = 2**(-1/6) sigma_scale = 2**(-1/6)
......
wrappers/python/tests/TestCharmmFiles.py
View file @ 0c646ace
...@@ -361,6 +361,47 @@ class TestCharmmFiles(unittest.TestCase): ...@@ -361,6 +361,47 @@ class TestCharmmFiles(unittest.TestCase):
# Compare with value computed with NAMD. # Compare with value computed with NAMD.
self.assertAlmostEqual(energy, -2154.5539, delta=1e-3*abs(energy)) 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
bH_water = [[0, 1], [1, 2]]
aH_water = [[0, 2]]
bH_methanol = [[3, 4], [3, 5], [3, 6], [7, 8]]
bCO_methanol = [[3, 7]]
aH_methanol = [[4, 5], [4, 6], [5, 6], [3, 8]]
system = psf.createSystem(params, constraints=None, rigidWater=False)
assert system.getNumConstraints() == 0
system = psf.createSystem(params, constraints=None, rigidWater=True)
assert system.getNumConstraints() == 3
for i in range(3):
assert system.getConstraintParameters(i)[:2] in bH_water + aH_water
system = psf.createSystem(params, constraints=HBonds, rigidWater=False)
assert system.getNumConstraints() == 6
for i in range(6):
assert system.getConstraintParameters(i)[:2] in bH_water + bH_methanol
system = psf.createSystem(params, constraints=HBonds, rigidWater=True)
assert system.getNumConstraints() == 7
for i in range(7):
assert system.getConstraintParameters(i)[:2] in bH_water + aH_water + bH_methanol
system = psf.createSystem(params, constraints=AllBonds, rigidWater=False)
assert system.getNumConstraints() == 7
for i in range(7):
assert system.getConstraintParameters(i)[:2] in bH_water + bH_methanol + bCO_methanol
system = psf.createSystem(params, constraints=AllBonds, rigidWater=True)
assert system.getNumConstraints() == 8
for i in range(8):
assert system.getConstraintParameters(i)[:2] in bH_water + aH_water + bH_methanol + bCO_methanol
system = psf.createSystem(params, constraints=HAngles, rigidWater=False)
assert system.getNumConstraints() == 12
for i in range(12):
assert system.getConstraintParameters(i)[:2] in bH_water + aH_water + bH_methanol + bCO_methanol + aH_methanol
system = psf.createSystem(params, constraints=HAngles, rigidWater=True)
assert system.getNumConstraints() == 12
for i in range(12):
assert system.getConstraintParameters(i)[:2] in bH_water + aH_water + bH_methanol + bCO_methanol + aH_methanol
if __name__ == '__main__': if __name__ == '__main__':
unittest.main() unittest.main()
......
wrappers/python/tests/systems/water_methanol.prm 0 → 100644
View file @ 0c646ace
* 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 @ 0c646ace
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