Skip to content

GitLab

Unverified Commit 6e073b48 authored by peastman's avatar peastman Committed by GitHub
Browse files

Merge pull request #2317 from jing-huang/master 

Support CHARMM psf that contains colinear type of lonepair.
parents 325fead1 75bedabe
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 3528 additions and 777 deletions
wrappers/python/CMakeLists.txt
View file @ 6e073b48
...@@ -59,6 +59,7 @@ foreach(SUBDIR ${SUBDIRS}) ...@@ -59,6 +59,7 @@ foreach(SUBDIR ${SUBDIRS})
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.cif" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.cif"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prm" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prm"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rtf"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.crd" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.crd"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.gro" "${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.gro"
......
wrappers/python/simtk/openmm/app/charmmpsffile.py
View file @ 6e073b48
...@@ -366,25 +366,30 @@ class CharmmPsfFile(object): ...@@ -366,25 +366,30 @@ class CharmmPsfFile(object):
holder = psfsections['NUMLP NUMLPH'][1] holder = psfsections['NUMLP NUMLPH'][1]
lonepair_list = TrackedList() lonepair_list = TrackedList()
if numlp != 0 or numlph != 0: if numlp != 0 or numlph != 0:
lp_hostnum_list=[]
lp_distance_list=[] lp_distance_list=[]
lp_angle_list=[] lp_angle_list=[]
lp_dihe_list=[] lp_dihe_list=[]
for i in range(numlp): for i in range(numlp):
lpline = holder[i].split() lpline = holder[i].split()
if len(lpline)!=6 or lpline[0] != '3' or lpline[2] != 'F' or int(lpline[1]) != 4*i+1 : if len(lpline)!=6 or lpline[2] != 'F' :
raise CharmmPSFError('Lonepair format error') raise CharmmPSFError('Lonepair format error')
else: else:
lp_hostnum_list.append(int(lpline[0]))
lp_distance_list.append(float(lpline[3])) lp_distance_list.append(float(lpline[3]))
lp_angle_list.append(float(lpline[4])) lp_angle_list.append(float(lpline[4]))
lp_dihe_list.append(float(lpline[5])) lp_dihe_list.append(float(lpline[5]))
lp_atom_counter=0
for i in range(numlp): for i in range(numlp):
lpline = holder[int(i/2)+numlp].split() idall=[]
icolumn = (i%2) * 4 for j in range(lp_hostnum_list[i]+1):
id1=int(lpline[icolumn]) -1 iline = int((lp_atom_counter+j)/8)+numlp
id2=int(lpline[icolumn+1])-1 icolumn = (lp_atom_counter+j)%8
id3=int(lpline[icolumn+2])-1 idall.append(int(holder[iline].split()[icolumn])-1)
id4=int(lpline[icolumn+3])-1 if len(idall)==3:
lonepair_list.append([id1, id2, id3, id4, lp_distance_list[i], lp_angle_list[i], lp_dihe_list[i]]) idall.append(-1) # use id4=-1 to mark colinear
lonepair_list.append([idall[0], idall[1], idall[2], idall[3], lp_distance_list[i], lp_angle_list[i], lp_dihe_list[i]])
lp_atom_counter += lp_hostnum_list[i]+1
# In Drude psf, here comes anisotropic section # In Drude psf, here comes anisotropic section
if IsDrudePSF: if IsDrudePSF:
numaniso = psfsections['NUMANISO'][0] numaniso = psfsections['NUMANISO'][0]
...@@ -709,7 +714,8 @@ class CharmmPsfFile(object): ...@@ -709,7 +714,8 @@ class CharmmPsfFile(object):
if atom.type is not None: if atom.type is not None:
# This is the most reliable way of determining the element # This is the most reliable way of determining the element
atomic_num = atom.type.atomic_number atomic_num = atom.type.atomic_number
elem = element.Element.getByAtomicNumber(atomic_num) if atomic_num != 0:
elem = element.Element.getByAtomicNumber(atomic_num)
else: else:
# Figure it out from the mass # Figure it out from the mass
elem = element.Element.getByMass(atom.mass) elem = element.Element.getByMass(atom.mass)
...@@ -906,23 +912,34 @@ class CharmmPsfFile(object): ...@@ -906,23 +912,34 @@ class CharmmPsfFile(object):
bond.bond_type.req*length_conv) bond.bond_type.req*length_conv)
# Add virtual sites # Add virtual sites
if verbose: print('Adding lonepairs...') if hasattr(self, 'lonepair_list'):
for lpsite in self.lonepair_list: if verbose: print('Adding lonepairs...')
index=lpsite[0] for lpsite in self.lonepair_list:
atom1=lpsite[1] index=lpsite[0]
atom2=lpsite[2] atom1=lpsite[1]
atom3=lpsite[3] atom2=lpsite[2]
if lpsite[4] > 0 : #relative atom3=lpsite[3]
r = lpsite[4] /10.0 # in nanometer if atom3 > 0:
xweights = [-1.0, 0.0, 1.0] if lpsite[4] > 0 : # relative lonepair type
elif lpsite[4] < 0: # bisector r = lpsite[4] /10.0 # in nanometer
r = lpsite[4] / (-10.0) xweights = [-1.0, 0.0, 1.0]
xweights = [-1.0, 0.5, 0.5] elif lpsite[4] < 0: # bisector lonepair type
theta = lpsite[5] * pi / 180.0 r = lpsite[4] / (-10.0)
phi = (180.0 - lpsite[6]) * pi / 180.0 xweights = [-1.0, 0.5, 0.5]
p = [r*cos(theta), r*sin(theta)*cos(phi), r*sin(theta)*sin(phi)] theta = lpsite[5] * pi / 180.0
p = [x if abs(x) > 1e-10 else 0 for x in p] # Avoid tiny numbers caused by roundoff error phi = (180.0 - lpsite[6]) * pi / 180.0
system.setVirtualSite(index, mm.LocalCoordinatesSite(atom1, atom3, atom2, mm.Vec3(1.0, 0.0, 0.0), mm.Vec3(xweights[0],xweights[1],xweights[2]), mm.Vec3(0.0, -1.0, 1.0), mm.Vec3(p[0],p[1],p[2]))) p = [r*cos(theta), r*sin(theta)*cos(phi), r*sin(theta)*sin(phi)]
p = [x if abs(x) > 1e-10 else 0 for x in p] # Avoid tiny numbers caused by roundoff error
system.setVirtualSite(index, mm.LocalCoordinatesSite(atom1, atom3, atom2, mm.Vec3(1.0, 0.0, 0.0), mm.Vec3(xweights[0],xweights[1],xweights[2]), mm.Vec3(0.0, -1.0, 1.0), mm.Vec3(p[0],p[1],p[2])))
else: # colinear lonepair type
# find a real atom to be the third one for LocalCoordinatesSite
for bond in self.bond_list:
if (bond.atom1.idx == atom2 and bond.atom2.idx != atom1):
a3=bond.atom2.idx
elif (bond.atom2.idx == atom2 and bond.atom1.idx != atom1):
a3=bond.atom1.idx
r = lpsite[4] / 10.0 # in nanometer
system.setVirtualSite(index, mm.LocalCoordinatesSite(atom1, atom2, a3, mm.Vec3(1.0, 0.0, 0.0), mm.Vec3(1.0,-1.0, 0.0), mm.Vec3(0.0, -1.0, 1.0), mm.Vec3(r,0.0,0.0)))
# Add Bond forces # Add Bond forces
if verbose: print('Adding bonds...') if verbose: print('Adding bonds...')
force = mm.HarmonicBondForce() force = mm.HarmonicBondForce()
......
wrappers/python/tests/TestCharmmFiles.py
View file @ 6e073b48
...@@ -116,8 +116,8 @@ class TestCharmmFiles(unittest.TestCase): ...@@ -116,8 +116,8 @@ class TestCharmmFiles(unittest.TestCase):
ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole) ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
self.assertAlmostEqual(ene, 15490.0033559, delta=0.05) self.assertAlmostEqual(ene, 15490.0033559, delta=0.05)
def test_drude(self): def test_Drude(self):
"""Tests CHARMM systems with Drude force field""" """Test CHARMM systems with Drude force field"""
warnings.filterwarnings('ignore', category=CharmmPSFWarning) warnings.filterwarnings('ignore', category=CharmmPSFWarning)
psf = CharmmPsfFile('systems/ala3_solv_drude.psf') psf = CharmmPsfFile('systems/ala3_solv_drude.psf')
crd = CharmmCrdFile('systems/ala3_solv_drude.crd') crd = CharmmCrdFile('systems/ala3_solv_drude.crd')
...@@ -136,6 +136,29 @@ class TestCharmmFiles(unittest.TestCase): ...@@ -136,6 +136,29 @@ class TestCharmmFiles(unittest.TestCase):
ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole) ene = state.getPotentialEnergy().value_in_unit(kilocalories_per_mole)
self.assertAlmostEqual(ene, -1831.54, delta=0.5) self.assertAlmostEqual(ene, -1831.54, delta=0.5)
def test_Lonepair(self):
"""Test the lonepair facilities, in particular the colinear type of lonepairs"""
warnings.filterwarnings('ignore', category=CharmmPSFWarning)
psf = CharmmPsfFile('systems/chlb_cgenff.psf')
crd = CharmmCrdFile('systems/chlb_cgenff.crd')
params = CharmmParameterSet('systems/top_all36_cgenff.rtf',
'systems/par_all36_cgenff.prm')
plat = Platform.getPlatformByName('Reference')
system = psf.createSystem(params)
con = Context(system, VerletIntegrator(2*femtoseconds), plat)
con.setPositions(crd.positions)
init_coor = con.getState(getPositions=True).getPositions()
# move the position of the lonepair and recompute its coordinates
plp=12
crd.positions[plp] = Vec3(0.5, 1.0, 1.5) * angstrom
con.setPositions(crd.positions)
con.computeVirtualSites()
new_coor = con.getState(getPositions=True).getPositions()
self.assertAlmostEqual(init_coor[plp][0]/nanometers, new_coor[plp][0]/nanometers)
self.assertAlmostEqual(init_coor[plp][1]/nanometers, new_coor[plp][1]/nanometers)
self.assertAlmostEqual(init_coor[plp][2]/nanometers, new_coor[plp][2]/nanometers)
def test_InsCode(self): def test_InsCode(self):
""" Test the parsing of PSF files that contain insertion codes in their residue numbers """ """ Test the parsing of PSF files that contain insertion codes in their residue numbers """
psf = CharmmPsfFile('systems/4TVP-dmj_wat-ion.psf') psf = CharmmPsfFile('systems/4TVP-dmj_wat-ion.psf')
......
wrappers/python/tests/systems/chlb_cgenff.crd 0 → 100644
View file @ 6e073b48
* GENERATE CPD ACCORDING TO IC TABLES IN FF
* USAGE: RUNCHM READ RESNAME= SEED= FF= ARHL=
* DATE: 6/14/19 14:22:39 CREATED BY USER: xu
*
13 EXT
1 1 CHLB C1 1.3878616207 -0.0212925645 0.0000000000 CHLB 1 -0.1000000000
2 1 CHLB H1 1.9158187876 -0.9626017904 0.0000000000 CHLB 1 0.1500000000
3 1 CHLB C2 2.0909297432 1.1920897951 -0.0000000000 CHLB 1 -0.1150000000
4 1 CHLB H2 3.1713458379 1.1910776735 -0.0000000000 CHLB 1 0.1150000000
5 1 CHLB C3 1.3891041046 2.4061558141 -0.0000000000 CHLB 1 -0.1150000000
6 1 CHLB H3 1.9293927824 3.3420080347 -0.0000000000 CHLB 1 0.1150000000
7 1 CHLB C4 -0.0132207406 2.4068634907 -0.0000000000 CHLB 1 -0.1150000000
8 1 CHLB H4 -0.5543440865 3.3420028049 -0.0000000000 CHLB 1 0.1150000000
9 1 CHLB C5 -0.7124562630 1.1912684987 0.0000000000 CHLB 1 -0.1000000000
10 1 CHLB H5 -1.7916319741 1.1777934579 0.0000000000 CHLB 1 0.1500000000
11 1 CHLB C6 -0.0118030037 -0.0204026447 0.0000000000 CHLB 1 0.0600000000
12 1 CHLB CL -0.8809968380 -1.5259625634 0.0000000000 CHLB 1 -0.2100000000
13 1 CHLB LP -1.7009674118 -2.9462612145 0.0000000000 CHLB 1 0.0500000000
wrappers/python/tests/systems/chlb_cgenff.psf 0 → 100644
View file @ 6e073b48
PSF EXT CMAP CHEQ XPLOR
3 !NTITLE
* GENERATE CPD ACCORDING TO IC TABLES IN FF
* USAGE: RUNCHM READ RESNAME= SEED= FF= ARHL=
* DATE: 6/14/19 14:22:39 CREATED BY USER: xu
13 !NATOM
1 CHLB 1 CHLB C1 CG2R61 -0.100000 12.0110 0 0.00000 -0.301140E-02
2 CHLB 1 CHLB H1 HGR62 0.150000 1.00800 0 0.00000 -0.301140E-02
3 CHLB 1 CHLB C2 CG2R61 -0.115000 12.0110 0 0.00000 -0.301140E-02
4 CHLB 1 CHLB H2 HGR61 0.115000 1.00800 0 0.00000 -0.301140E-02
5 CHLB 1 CHLB C3 CG2R61 -0.115000 12.0110 0 0.00000 -0.301140E-02
6 CHLB 1 CHLB H3 HGR61 0.115000 1.00800 0 0.00000 -0.301140E-02
7 CHLB 1 CHLB C4 CG2R61 -0.115000 12.0110 0 0.00000 -0.301140E-02
8 CHLB 1 CHLB H4 HGR61 0.115000 1.00800 0 0.00000 -0.301140E-02
9 CHLB 1 CHLB C5 CG2R61 -0.100000 12.0110 0 0.00000 -0.301140E-02
10 CHLB 1 CHLB H5 HGR62 0.150000 1.00800 0 0.00000 -0.301140E-02
11 CHLB 1 CHLB C6 CG2R61 0.600000E-01 12.0110 0 0.00000 -0.301140E-02
12 CHLB 1 CHLB CL CLGR1 -0.210000 35.4530 0 0.00000 -0.301140E-02
13 CHLB 1 CHLB LP LPH 0.500000E-01 0.00000 -1 0.00000 -0.301140E-02
13 !NBOND: bonds
1 2 1 3 3 4 3 5
5 6 5 7 7 8 7 9
9 10 9 11 11 1 11 12
12 13
18 !NTHETA: angles
2 1 3 2 1 11 3 1 11
1 3 4 1 3 5 4 3 5
3 5 6 3 5 7 6 5 7
5 7 8 5 7 9 8 7 9
7 9 10 7 9 11 10 9 11
1 11 9 1 11 12 9 11 12
24 !NPHI: dihedrals
1 3 5 6 1 3 5 7
1 11 9 7 1 11 9 10
2 1 3 4 2 1 3 5
2 1 11 9 2 1 11 12
3 1 11 9 3 1 11 12
3 5 7 8 3 5 7 9
4 3 1 11 4 3 5 6
4 3 5 7 5 3 1 11
5 7 9 10 5 7 9 11
6 5 7 8 6 5 7 9
7 9 11 12 8 7 9 10
8 7 9 11 10 9 11 12
0 !NIMPHI: impropers
0 !NDON: donors
0 !NACC: acceptors
0 !NNB
0 0 0 0 0 0 0 0
0 0 0 0 0
1 0 !NGRP NST2
0 1 0
1 !MOLNT
1 1 1 1 1 1 1 1
1 1 1 1 1
1 3 !NUMLP NUMLPH
2 1 F 1.64000 0.00000 0.00000
13 12 11
0 !NCRTERM: cross-terms
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