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

Merge pull request #2101 from peastman/top

Support more features in Gromacs top files
parents 9dac72f5 c30fa8e6
...@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of ...@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org. Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2016 Stanford University and the Authors. Portions copyright (c) 2012-2018 Stanford University and the Authors.
Authors: Peter Eastman Authors: Peter Eastman
Contributors: Jason Swails Contributors: Jason Swails
...@@ -115,6 +115,7 @@ class GromacsTopFile(object): ...@@ -115,6 +115,7 @@ class GromacsTopFile(object):
self.dihedrals = [] self.dihedrals = []
self.exclusions = [] self.exclusions = []
self.pairs = [] self.pairs = []
self.constraints = []
self.cmaps = [] self.cmaps = []
self.vsites2 = [] self.vsites2 = []
self.has_virtual_sites = False self.has_virtual_sites = False
...@@ -236,6 +237,8 @@ class GromacsTopFile(object): ...@@ -236,6 +237,8 @@ class GromacsTopFile(object):
self._processExclusion(line) self._processExclusion(line)
elif self._currentCategory == 'pairs': elif self._currentCategory == 'pairs':
self._processPair(line) self._processPair(line)
elif self._currentCategory == 'constraints':
self._processConstraint(line)
elif self._currentCategory == 'cmap': elif self._currentCategory == 'cmap':
self._processCmap(line) self._processCmap(line)
elif self._currentCategory == 'atomtypes': elif self._currentCategory == 'atomtypes':
...@@ -274,7 +277,7 @@ class GromacsTopFile(object): ...@@ -274,7 +277,7 @@ class GromacsTopFile(object):
raise ValueError('Too few fields in [ defaults ] line: '+line) raise ValueError('Too few fields in [ defaults ] line: '+line)
if fields[0] != '1': if fields[0] != '1':
raise ValueError('Unsupported nonbonded type: '+fields[0]) raise ValueError('Unsupported nonbonded type: '+fields[0])
if not (fields[1] == '2' or fields[1] == '1'): if not fields[1] in ('1', '2', '3'):
raise ValueError('Unsupported combination rule: '+fields[1]) raise ValueError('Unsupported combination rule: '+fields[1])
if fields[2].lower() == 'no': if fields[2].lower() == 'no':
self._genpairs = False self._genpairs = False
...@@ -334,7 +337,7 @@ class GromacsTopFile(object): ...@@ -334,7 +337,7 @@ class GromacsTopFile(object):
fields = line.split() fields = line.split()
if len(fields) < 5: if len(fields) < 5:
raise ValueError('Too few fields in [ dihedrals ] line: '+line) raise ValueError('Too few fields in [ dihedrals ] line: '+line)
if fields[4] not in ('1', '2', '3', '4', '9'): if fields[4] not in ('1', '2', '3', '4', '5', '9'):
raise ValueError('Unsupported function type in [ dihedrals ] line: '+line) raise ValueError('Unsupported function type in [ dihedrals ] line: '+line)
self._currentMoleculeType.dihedrals.append(fields) self._currentMoleculeType.dihedrals.append(fields)
...@@ -358,6 +361,15 @@ class GromacsTopFile(object): ...@@ -358,6 +361,15 @@ class GromacsTopFile(object):
raise ValueError('Unsupported function type in [ pairs ] line: '+line) raise ValueError('Unsupported function type in [ pairs ] line: '+line)
self._currentMoleculeType.pairs.append(fields) self._currentMoleculeType.pairs.append(fields)
def _processConstraint(self, line):
"""Process a line in the [ constraints ] category."""
if self._currentMoleculeType is None:
raise ValueError('Found [ constraints ] section before [ moleculetype ]')
fields = line.split()
if len(fields) < 4:
raise ValueError('Too few fields in [ constraints ] line: '+line)
self._currentMoleculeType.constraints.append(fields)
def _processCmap(self, line): def _processCmap(self, line):
"""Process a line in the [ cmaps ] category.""" """Process a line in the [ cmaps ] category."""
if self._currentMoleculeType is None: if self._currentMoleculeType is None:
...@@ -408,7 +420,7 @@ class GromacsTopFile(object): ...@@ -408,7 +420,7 @@ class GromacsTopFile(object):
fields = line.split() fields = line.split()
if len(fields) < 7: if len(fields) < 7:
raise ValueError('Too few fields in [ dihedraltypes ] line: '+line) raise ValueError('Too few fields in [ dihedraltypes ] line: '+line)
if fields[4] not in ('1', '2', '3', '4', '9'): if fields[4] not in ('1', '2', '3', '4', '5', '9'):
raise ValueError('Unsupported function type in [ dihedraltypes ] line: '+line) raise ValueError('Unsupported function type in [ dihedraltypes ] line: '+line)
key = tuple(fields[:5]) key = tuple(fields[:5])
if fields[4] == '9' and key in self._dihedralTypes: if fields[4] == '9' and key in self._dihedralTypes:
...@@ -575,7 +587,7 @@ class GromacsTopFile(object): ...@@ -575,7 +587,7 @@ class GromacsTopFile(object):
def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer, def createSystem(self, nonbondedMethod=ff.NoCutoff, nonbondedCutoff=1.0*unit.nanometer,
constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None): constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, ewaldErrorTolerance=0.0005, removeCMMotion=True, hydrogenMass=None):
"""Construct an OpenMM System representing the topology described by this """Construct an OpenMM System representing the topology described by this
prmtop file. top file.
Parameters Parameters
---------- ----------
...@@ -587,6 +599,8 @@ class GromacsTopFile(object): ...@@ -587,6 +599,8 @@ class GromacsTopFile(object):
constraints : object=None constraints : object=None
Specifies which bonds and angles should be implemented with Specifies which bonds and angles should be implemented with
constraints. Allowed values are None, HBonds, AllBonds, or HAngles. constraints. Allowed values are None, HBonds, AllBonds, or HAngles.
Regardless of this value, constraints that are explicitly specified
in the top file will always be included.
rigidWater : boolean=True rigidWater : boolean=True
If true, water molecules will be fully rigid regardless of the value If true, water molecules will be fully rigid regardless of the value
passed for the constraints argument passed for the constraints argument
...@@ -622,8 +636,8 @@ class GromacsTopFile(object): ...@@ -622,8 +636,8 @@ class GromacsTopFile(object):
raise ValueError('Illegal nonbonded method for a non-periodic system') raise ValueError('Illegal nonbonded method for a non-periodic system')
nb = mm.NonbondedForce() nb = mm.NonbondedForce()
sys.addForce(nb) sys.addForce(nb)
if self._defaults[1] == '1': if self._defaults[1] in ('1', '3'):
lj = mm.CustomNonbondedForce('sqrt(A1*A2)/r^12-sqrt(C1*C2)/r^6') lj = mm.CustomNonbondedForce('A1*A2/r^12-C1*C2/r^6')
lj.addPerParticleParameter('C') lj.addPerParticleParameter('C')
lj.addPerParticleParameter('A') lj.addPerParticleParameter('A')
sys.addForce(lj) sys.addForce(lj)
...@@ -802,7 +816,7 @@ class GromacsTopFile(object): ...@@ -802,7 +816,7 @@ class GromacsTopFile(object):
dihedralType = fields[4] dihedralType = fields[4]
reversedTypes = types[::-1]+(dihedralType,) reversedTypes = types[::-1]+(dihedralType,)
types = types+(dihedralType,) types = types+(dihedralType,)
if (dihedralType in ('1', '4', '9') and len(fields) > 7) or (dihedralType == '3' and len(fields) > 10) or (dihedralType == '2' and len(fields) > 6): if (dihedralType in ('1', '4', '5', '9') and len(fields) > 7) or (dihedralType == '3' and len(fields) > 10) or (dihedralType == '2' and len(fields) > 6):
paramsList = [fields] paramsList = [fields]
else: else:
# Look for a matching dihedral type. # Look for a matching dihedral type.
...@@ -847,6 +861,9 @@ class GromacsTopFile(object): ...@@ -847,6 +861,9 @@ class GromacsTopFile(object):
if rb is None: if rb is None:
rb = mm.RBTorsionForce() rb = mm.RBTorsionForce()
sys.addForce(rb) sys.addForce(rb)
if dihedralType == '5':
# Convert Fourier coefficients to RB coefficients.
c = [c[1]+0.5*(c[0]+c[2]), 0.5*(-c[0]+3*c[2]), -c[1]+4*c[3], -2*c[2], -4*c[3], 0]
rb.addTorsion(baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], c[0], c[1], c[2], c[3], c[4], c[5]) rb.addTorsion(baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], c[0], c[1], c[2], c[3], c[4], c[5])
# Add CMAP terms. # Add CMAP terms.
...@@ -896,9 +913,14 @@ class GromacsTopFile(object): ...@@ -896,9 +913,14 @@ class GromacsTopFile(object):
else: else:
if self._defaults[1] == '1': if self._defaults[1] == '1':
nb.addParticle(q, 1.0, 0.0) nb.addParticle(q, 1.0, 0.0)
lj.addParticle([float(params[6]), float(params[7])]) lj.addParticle([math.sqrt(float(params[6])), math.sqrt(float(params[7]))])
elif self._defaults[1] == '2': elif self._defaults[1] == '2':
nb.addParticle(q, float(params[6]), float(params[7])) nb.addParticle(q, float(params[6]), float(params[7]))
elif self._defaults[1] == '3':
nb.addParticle(q, 1.0, 0.0)
sigma = float(params[6])
epsilon = float(params[7])
lj.addParticle([math.sqrt(4*epsilon*sigma**6), math.sqrt(4*epsilon*sigma**12)])
for fields in moleculeType.atoms: for fields in moleculeType.atoms:
if implicitSolvent is OBC2: if implicitSolvent is OBC2:
...@@ -909,6 +931,10 @@ class GromacsTopFile(object): ...@@ -909,6 +931,10 @@ class GromacsTopFile(object):
for fields in moleculeType.bonds: for fields in moleculeType.bonds:
atoms = [int(x)-1 for x in fields[:2]] atoms = [int(x)-1 for x in fields[:2]]
bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1])) bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1]))
for fields in moleculeType.constraints:
if fields[2] == '1':
atoms = [int(x)-1 for x in fields[:2]]
bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1]))
if has_nbfix_terms: if has_nbfix_terms:
for fields in moleculeType.bonds: for fields in moleculeType.bonds:
atoms = [int(x)-1 for x in fields[:2]] atoms = [int(x)-1 for x in fields[:2]]
...@@ -932,19 +958,29 @@ class GromacsTopFile(object): ...@@ -932,19 +958,29 @@ class GromacsTopFile(object):
for fields in moleculeType.pairs: for fields in moleculeType.pairs:
atoms = [int(x)-1 for x in fields[:2]] atoms = [int(x)-1 for x in fields[:2]]
types = tuple(atomTypes[i] for i in atoms) types = tuple(atomTypes[i] for i in atoms)
atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0])
atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
atom1params = [x.value_in_unit_system(unit.md_unit_system) for x in atom1params]
atom2params = [x.value_in_unit_system(unit.md_unit_system) for x in atom2params]
if len(fields) >= 5: if len(fields) >= 5:
params = fields[3:5] params = [float(x) for x in fields[3:5]]
elif types in self._pairTypes: elif types in self._pairTypes:
params = self._pairTypes[types][3:5] params = [float(x) for x in self._pairTypes[types][3:5]]
elif types[::-1] in self._pairTypes: elif types[::-1] in self._pairTypes:
params = self._pairTypes[types[::-1]][3:5] params = [float(x) for x in self._pairTypes[types[::-1]][3:5]]
elif not self._genpairs: elif not self._genpairs:
raise ValueError('No pair parameters defined for atom ' raise ValueError('No pair parameters defined for atom '
'types %s and gen-pairs is "no"' % types) 'types %s and gen-pairs is "no"' % types)
elif has_nbfix_terms:
continue
else: else:
continue # We'll use the automatically generated parameters # Generate the parameters based on the atom parameters.
atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0]) if self._defaults[1] == '2':
atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1]) params = [0.5*(atom1params[1]+atom2params[1]), fudgeLJ*math.sqrt(atom1params[2]*atom2params[2])]
else:
atom1lj = lj.getParticleParameters(baseAtomIndex+atoms[0])
atom2lj = lj.getParticleParameters(baseAtomIndex+atoms[1])
params = [fudgeLJ*atom1lj[0]*atom2lj[0], fudgeLJ*atom1lj[1]*atom2lj[1]]
pairs.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1])) pairs.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
for fields in moleculeType.exclusions: for fields in moleculeType.exclusions:
atoms = [int(x)-1 for x in fields] atoms = [int(x)-1 for x in fields]
...@@ -960,6 +996,13 @@ class GromacsTopFile(object): ...@@ -960,6 +996,13 @@ class GromacsTopFile(object):
vsite = mm.TwoParticleAverageSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], (1-c1), c1) vsite = mm.TwoParticleAverageSite(baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], (1-c1), c1)
sys.setVirtualSite(baseAtomIndex+atoms[0], vsite) sys.setVirtualSite(baseAtomIndex+atoms[0], vsite)
# Add explicitly specified constraints.
for fields in moleculeType.constraints:
atoms = [int(x)-1 for x in fields[:2]]
length = float(fields[2])
sys.addConstraint(baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], length)
# Create nonbonded exceptions. # Create nonbonded exceptions.
if not has_nbfix_terms: if not has_nbfix_terms:
...@@ -1012,24 +1055,23 @@ class GromacsTopFile(object): ...@@ -1012,24 +1055,23 @@ class GromacsTopFile(object):
for exclusion in exclusions: for exclusion in exclusions:
nb.addException(exclusion[0], exclusion[1], 0.0, 1.0, 0.0, True) nb.addException(exclusion[0], exclusion[1], 0.0, 1.0, 0.0, True)
if self._defaults[1] == '1': if self._defaults[1] in ('1', '3'):
# We're using a CustomNonbondedForce for LJ interactions, so also create a CustomBondForce # We're using a CustomNonbondedForce for LJ interactions, so also create a CustomBondForce
# to handle the exceptions. # to handle the exceptions.
pair_bond = mm.CustomBondForce('138.935456*q/r-C/r^6+A/r^12') pair_bond = mm.CustomBondForce('-C/r^6+A/r^12')
pair_bond.addPerBondParameter('q')
pair_bond.addPerBondParameter('C') pair_bond.addPerBondParameter('C')
pair_bond.addPerBondParameter('A') pair_bond.addPerBondParameter('A')
sys.addForce(pair_bond) sys.addForce(pair_bond)
lj.createExclusionsFromBonds(bondIndices, 3) lj.createExclusionsFromBonds(bondIndices, 3)
for pair in pairs: for pair in pairs:
nb.addException(pair[0], pair[1], pair[2], 1.0, 0.0, True) nb.addException(pair[0], pair[1], pair[2], 1.0, 0.0, True)
pair_bond.addBond(pair[0], pair[1], [pair[2],float(pair[3]), float(pair[4])]) pair_bond.addBond(pair[0], pair[1], [pair[3], pair[4]])
for exclusion in exclusions: for exclusion in exclusions:
lj.addExclusion(exclusion[0], exclusion[1]) lj.addExclusion(exclusion[0], exclusion[1])
elif self._defaults[1] == '2': elif self._defaults[1] == '2':
for pair in pairs: for pair in pairs:
nb.addException(pair[0], pair[1], pair[2], float(pair[3]), float(pair[4]), True) nb.addException(pair[0], pair[1], pair[2], pair[3], pair[4], True)
# Finish configuring the NonbondedForce. # Finish configuring the NonbondedForce.
...@@ -1042,7 +1084,7 @@ class GromacsTopFile(object): ...@@ -1042,7 +1084,7 @@ class GromacsTopFile(object):
nb.setNonbondedMethod(methodMap[nonbondedMethod]) nb.setNonbondedMethod(methodMap[nonbondedMethod])
nb.setCutoffDistance(nonbondedCutoff) nb.setCutoffDistance(nonbondedCutoff)
nb.setEwaldErrorTolerance(ewaldErrorTolerance) nb.setEwaldErrorTolerance(ewaldErrorTolerance)
if self._defaults[1] == '1': if self._defaults[1] in ('1', '3'):
methodMap = {ff.NoCutoff:mm.CustomNonbondedForce.NoCutoff, methodMap = {ff.NoCutoff:mm.CustomNonbondedForce.NoCutoff,
ff.CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic, ff.CutoffNonPeriodic:mm.CustomNonbondedForce.CutoffNonPeriodic,
ff.CutoffPeriodic:mm.CustomNonbondedForce.CutoffPeriodic, ff.CutoffPeriodic:mm.CustomNonbondedForce.CutoffPeriodic,
......
...@@ -63,6 +63,22 @@ class TestGromacsTopFile(unittest.TestCase): ...@@ -63,6 +63,22 @@ class TestGromacsTopFile(unittest.TestCase):
ene = context.getState(getEnergy=True).getPotentialEnergy() ene = context.getState(getEnergy=True).getPotentialEnergy()
self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), -346.940915296) self.assertAlmostEqual(ene.value_in_unit(kilojoules_per_mole), -346.940915296)
def test_ionic(self):
"""Test simulating an ionic liquid"""
gro = GromacsGroFile('systems/ionic.gro')
top = GromacsTopFile('systems/ionic.top', periodicBoxVectors=gro.getPeriodicBoxVectors())
system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1.2)
for f in system.getForces():
if isinstance(f, CustomNonbondedForce):
f.setUseLongRangeCorrection(True)
context = Context(system, VerletIntegrator(1*femtosecond),
Platform.getPlatformByName('Reference'))
context.setPositions(gro.positions)
energy = context.getState(getEnergy=True).getPotentialEnergy().value_in_unit(kilojoules_per_mole)
self.assertAlmostEqual(energy, 3135.33, delta=energy*0.005)
self.assertEqual(1400, system.getNumConstraints())
def test_Cutoff(self): def test_Cutoff(self):
"""Test to make sure the nonbondedCutoff parameter is passed correctly.""" """Test to make sure the nonbondedCutoff parameter is passed correctly."""
......
This source diff could not be displayed because it is too large. You can view the blob instead.
; created by fftool
[ defaults ]
; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ
1 3 yes 0.5 0.5
[ atomtypes ]
; name at.nr mass charge ptype sigma epsilon
CBT 6 12.0110 0.3500 A 3.50000e-01 2.76140e-01
F1 9 18.9980 -0.1600 A 3.11800e-01 2.55400e-01
SBT 16 32.0660 1.0200 A 3.55000e-01 1.04600e+00
NBT 7 14.0000 -0.6600 A 3.25000e-01 7.11280e-01
OBT 8 15.9990 -0.5300 A 3.15000e-01 8.37360e-01
NAP 7 14.0070 0.1500 A 3.25000e-01 7.11280e-01
CAPO 6 12.0110 0.0000 A 3.55000e-01 2.92880e-01
CAPM 6 12.0110 -0.0700 A 3.55000e-01 2.92880e-01
CAPP 6 12.0110 0.0200 A 3.55000e-01 2.92880e-01
C1 6 12.0110 -0.1700 A 3.50000e-01 2.76140e-01
HAP 1 1.0080 0.1500 A 2.42000e-01 1.25520e-01
C2 6 12.0110 0.0100 A 3.50000e-01 2.76140e-01
H1 1 1.0080 0.1300 A 2.50000e-01 1.25520e-01
CS 6 12.0110 -0.1200 A 3.50000e-01 2.76140e-01
HC 1 1.0080 0.0600 A 2.50000e-01 1.25520e-01
CT 6 12.0110 -0.1800 A 3.50000e-01 2.76140e-01
[ moleculetype ]
; name nrexcl
tf2N- 3
[ atoms ]
; nr type resnr residu atom cgnr charge
1 CBT 1 tf2 CBT 1 0.3500
2 F1 1 tf2 F1 1 -0.1600
3 F1 1 tf2 F1 1 -0.1600
4 F1 1 tf2 F1 1 -0.1600
5 SBT 1 tf2 SBT 1 1.0200
6 NBT 1 tf2 NBT 1 -0.6600
7 OBT 1 tf2 OBT 1 -0.5300
8 OBT 1 tf2 OBT 1 -0.5300
9 SBT 1 tf2 SBT 1 1.0200
10 OBT 1 tf2 OBT 1 -0.5300
11 OBT 1 tf2 OBT 1 -0.5300
12 CBT 1 tf2 CBT 1 0.3500
13 F1 1 tf2 F1 1 -0.1600
14 F1 1 tf2 F1 1 -0.1600
15 F1 1 tf2 F1 1 -0.1600
[ bonds ]
; ai aj func b0 kb
2 1 1 0.13230 369800.0
3 1 1 0.13230 369800.0
4 1 1 0.13230 369800.0
5 1 1 0.18180 195000.0
6 5 1 0.15700 313700.0
7 5 1 0.14370 533100.0
8 5 1 0.14370 533100.0
9 6 1 0.15700 313700.0
10 9 1 0.14370 533100.0
11 9 1 0.14370 533100.0
12 9 1 0.18180 195000.0
13 12 1 0.13230 369800.0
14 12 1 0.13230 369800.0
15 12 1 0.13230 369800.0
[ constraints ]
; ai aj func b0
[ angles ]
; ai aj ak func th0 cth
2 1 3 1 107.100 781.000
2 1 4 1 107.100 781.000
2 1 5 1 111.700 694.000
3 1 4 1 107.100 781.000
3 1 5 1 111.700 694.000
4 1 5 1 111.700 694.000
1 5 6 1 103.500 764.000
1 5 7 1 102.600 870.000
1 5 8 1 102.600 870.000
6 5 7 1 113.600 789.000
6 5 8 1 113.600 789.000
7 5 8 1 118.500 969.000
5 6 9 1 125.600 671.000
6 9 10 1 113.600 789.000
6 9 11 1 113.600 789.000
6 9 12 1 103.500 764.000
10 9 11 1 118.500 969.000
10 9 12 1 102.600 870.000
11 9 12 1 102.600 870.000
9 12 13 1 111.700 694.000
9 12 14 1 111.700 694.000
9 12 15 1 111.700 694.000
13 12 14 1 107.100 781.000
13 12 15 1 107.100 781.000
14 12 15 1 107.100 781.000
[ dihedrals ]
; ai aj ak al func coefficients
6 5 1 2 5 0.00000 0.00000 1.32200 0.00000
6 5 1 3 5 0.00000 0.00000 1.32200 0.00000
6 5 1 4 5 0.00000 0.00000 1.32200 0.00000
7 5 1 2 5 0.00000 0.00000 1.45100 0.00000
7 5 1 3 5 0.00000 0.00000 1.45100 0.00000
7 5 1 4 5 0.00000 0.00000 1.45100 0.00000
8 5 1 2 5 0.00000 0.00000 1.45100 0.00000
8 5 1 3 5 0.00000 0.00000 1.45100 0.00000
8 5 1 4 5 0.00000 0.00000 1.45100 0.00000
9 6 5 1 5 32.77300 -10.42000 -3.19500 0.00000
9 6 5 7 5 0.00000 0.00000 -0.01500 0.00000
9 6 5 8 5 0.00000 0.00000 -0.01500 0.00000
10 9 6 5 5 0.00000 0.00000 -0.01500 0.00000
11 9 6 5 5 0.00000 0.00000 -0.01500 0.00000
12 9 6 5 5 32.77300 -10.42000 -3.19500 0.00000
13 12 9 6 5 0.00000 0.00000 1.32200 0.00000
13 12 9 10 5 0.00000 0.00000 1.45100 0.00000
13 12 9 11 5 0.00000 0.00000 1.45100 0.00000
14 12 9 6 5 0.00000 0.00000 1.32200 0.00000
14 12 9 10 5 0.00000 0.00000 1.45100 0.00000
14 12 9 11 5 0.00000 0.00000 1.45100 0.00000
15 12 9 6 5 0.00000 0.00000 1.32200 0.00000
15 12 9 10 5 0.00000 0.00000 1.45100 0.00000
15 12 9 11 5 0.00000 0.00000 1.45100 0.00000
[ pairs ]
; ai aj func
6 2 1
6 3 1
6 4 1
7 2 1
7 3 1
7 4 1
8 2 1
8 3 1
8 4 1
9 1 1
9 7 1
9 8 1
10 5 1
11 5 1
12 5 1
13 6 1
13 10 1
13 11 1
14 6 1
14 10 1
14 11 1
15 6 1
15 10 1
15 11 1
[ moleculetype ]
; name nrexcl
bpyri+ 3
[ atoms ]
; nr type resnr residu atom cgnr charge
1 NAP 1 bpy NAP 1 0.1500
2 CAPO 1 bpy CAPO 1 0.0000
3 CAPM 1 bpy CAPM 1 -0.0700
4 CAPP 1 bpy CAPP 1 0.0200
5 CAPM 1 bpy CAPM 1 -0.0700
6 CAPO 1 bpy CAPO 1 0.0000
7 C1 1 bpy C1 1 -0.1700
8 HAP 1 bpy HAP 1 0.1500
9 HAP 1 bpy HAP 1 0.1500
10 HAP 1 bpy HAP 1 0.1500
11 HAP 1 bpy HAP 1 0.1500
12 HAP 1 bpy HAP 1 0.1500
13 C2 1 bpy C2 1 0.0100
14 H1 1 bpy H1 1 0.1300
15 H1 1 bpy H1 1 0.1300
16 CS 1 bpy CS 1 -0.1200
17 HC 1 bpy HC 1 0.0600
18 HC 1 bpy HC 1 0.0600
19 CT 1 bpy CT 1 -0.1800
20 HC 1 bpy HC 1 0.0600
21 HC 1 bpy HC 1 0.0600
22 HC 1 bpy HC 1 0.0600
23 HC 1 bpy HC 1 0.0600
24 HC 1 bpy HC 1 0.0600
[ bonds ]
; ai aj func b0 kb
2 1 1 0.13400 404200.0
3 2 1 0.14000 392460.0
4 3 1 0.14000 392460.0
5 4 1 0.14000 392460.0
6 5 1 0.14000 392460.0
7 1 1 0.14660 282000.0
13 7 1 0.15290 224200.0
16 13 1 0.15290 224200.0
19 16 1 0.15290 224200.0
6 1 1 0.13400 404200.0
[ constraints ]
; ai aj func b0
8 2 1 0.10800
9 3 1 0.10800
10 4 1 0.10800
11 5 1 0.10800
12 6 1 0.10800
14 7 1 0.10900
15 7 1 0.10900
17 13 1 0.10900
18 13 1 0.10900
20 16 1 0.10900
21 16 1 0.10900
22 19 1 0.10900
23 19 1 0.10900
24 19 1 0.10900
[ angles ]
; ai aj ak func th0 cth
2 1 7 1 119.800 585.800
2 1 6 1 120.400 585.800
7 1 6 1 119.800 585.800
1 2 3 1 120.000 585.800
1 2 8 1 120.000 292.900
3 2 8 1 120.000 292.900
2 3 4 1 120.000 527.200
2 3 9 1 120.000 292.900
4 3 9 1 120.000 292.900
3 4 5 1 120.000 527.200
3 4 10 1 120.000 292.900
5 4 10 1 120.000 292.900
4 5 6 1 120.000 527.200
4 5 11 1 120.000 292.900
6 5 11 1 120.000 292.900
5 6 12 1 120.000 292.900
5 6 1 1 120.000 585.800
12 6 1 1 120.000 292.900
1 7 13 1 112.700 488.300
1 7 14 1 110.700 313.800
1 7 15 1 110.700 313.800
13 7 14 1 110.700 313.800
13 7 15 1 110.700 313.800
14 7 15 1 107.800 276.100
7 13 16 1 112.700 488.300
7 13 17 1 110.700 313.800
7 13 18 1 110.700 313.800
16 13 17 1 110.700 313.800
16 13 18 1 110.700 313.800
17 13 18 1 107.800 276.100
13 16 19 1 112.700 488.300
13 16 20 1 110.700 313.800
13 16 21 1 110.700 313.800
19 16 20 1 110.700 313.800
19 16 21 1 110.700 313.800
20 16 21 1 107.800 276.100
16 19 22 1 110.700 313.800
16 19 23 1 110.700 313.800
16 19 24 1 110.700 313.800
22 19 23 1 107.800 276.100
22 19 24 1 107.800 276.100
23 19 24 1 107.800 276.100
[ dihedrals ]
; ai aj ak al func coefficients
3 2 1 7 5 0.00000 12.55200 0.00000 0.00000
3 2 1 6 5 0.00000 12.55200 0.00000 0.00000
8 2 1 7 5 0.00000 12.55200 0.00000 0.00000
8 2 1 6 5 0.00000 12.55200 0.00000 0.00000
4 3 2 1 5 0.00000 30.33400 0.00000 0.00000
4 3 2 8 5 0.00000 30.33400 0.00000 0.00000
9 3 2 1 5 0.00000 30.33400 0.00000 0.00000
9 3 2 8 5 0.00000 30.33400 0.00000 0.00000
5 4 3 2 5 0.00000 30.33400 0.00000 0.00000
5 4 3 9 5 0.00000 30.33400 0.00000 0.00000
10 4 3 2 5 0.00000 30.33400 0.00000 0.00000
10 4 3 9 5 0.00000 30.33400 0.00000 0.00000
6 5 4 3 5 0.00000 30.33400 0.00000 0.00000
6 5 4 10 5 0.00000 30.33400 0.00000 0.00000
11 5 4 3 5 0.00000 30.33400 0.00000 0.00000
11 5 4 10 5 0.00000 30.33400 0.00000 0.00000
12 6 5 4 5 0.00000 30.33400 0.00000 0.00000
12 6 5 11 5 0.00000 30.33400 0.00000 0.00000
1 6 5 4 5 0.00000 30.33400 0.00000 0.00000
1 6 5 11 5 0.00000 30.33400 0.00000 0.00000
13 7 1 2 5 0.00000 1.09200 0.00000 0.79300
13 7 1 6 5 0.00000 1.09200 0.00000 0.79300
14 7 1 2 5 0.00000 0.00000 0.00000 0.00000
14 7 1 6 5 0.00000 0.00000 0.00000 0.00000
15 7 1 2 5 0.00000 0.00000 0.00000 0.00000
15 7 1 6 5 0.00000 0.00000 0.00000 0.00000
16 13 7 1 5 -7.47970 3.16420 -1.20260 0.00000
16 13 7 14 5 0.00000 0.00000 1.25520 0.00000
16 13 7 15 5 0.00000 0.00000 1.25520 0.00000
17 13 7 1 5 0.00000 0.00000 0.36700 0.00000
17 13 7 14 5 0.00000 0.00000 1.25520 0.00000
17 13 7 15 5 0.00000 0.00000 1.25520 0.00000
18 13 7 1 5 0.00000 0.00000 0.36700 0.00000
18 13 7 14 5 0.00000 0.00000 1.25520 0.00000
18 13 7 15 5 0.00000 0.00000 1.25520 0.00000
19 16 13 7 5 5.43920 -0.20920 0.83680 0.00000
19 16 13 17 5 0.00000 0.00000 1.25520 0.00000
19 16 13 18 5 0.00000 0.00000 1.25520 0.00000
20 16 13 7 5 0.00000 0.00000 1.25520 0.00000
20 16 13 17 5 0.00000 0.00000 1.25520 0.00000
20 16 13 18 5 0.00000 0.00000 1.25520 0.00000
21 16 13 7 5 0.00000 0.00000 1.25520 0.00000
21 16 13 17 5 0.00000 0.00000 1.25520 0.00000
21 16 13 18 5 0.00000 0.00000 1.25520 0.00000
22 19 16 13 5 0.00000 0.00000 1.25520 0.00000
22 19 16 20 5 0.00000 0.00000 1.25520 0.00000
22 19 16 21 5 0.00000 0.00000 1.25520 0.00000
23 19 16 13 5 0.00000 0.00000 1.25520 0.00000
23 19 16 20 5 0.00000 0.00000 1.25520 0.00000
23 19 16 21 5 0.00000 0.00000 1.25520 0.00000
24 19 16 13 5 0.00000 0.00000 1.25520 0.00000
24 19 16 20 5 0.00000 0.00000 1.25520 0.00000
24 19 16 21 5 0.00000 0.00000 1.25520 0.00000
5 6 1 2 5 0.00000 12.55200 0.00000 0.00000
5 6 1 7 5 0.00000 12.55200 0.00000 0.00000
12 6 1 2 5 0.00000 12.55200 0.00000 0.00000
12 6 1 7 5 0.00000 12.55200 0.00000 0.00000
2 6 1 7 5 0.00000 9.20480 0.00000 0.00000
3 1 2 8 5 0.00000 9.20480 0.00000 0.00000
2 4 3 9 5 0.00000 9.20480 0.00000 0.00000
3 5 4 10 5 0.00000 9.20480 0.00000 0.00000
4 6 5 11 5 0.00000 9.20480 0.00000 0.00000
5 1 6 12 5 0.00000 9.20480 0.00000 0.00000
[ pairs ]
; ai aj func
3 7 1
3 6 1
8 7 1
8 6 1
4 1 1
4 8 1
9 1 1
9 8 1
5 2 1
5 9 1
10 2 1
10 9 1
6 10 1
11 3 1
11 10 1
12 4 1
12 11 1
1 11 1
13 2 1
13 6 1
14 2 1
14 6 1
15 2 1
15 6 1
16 1 1
16 14 1
16 15 1
17 1 1
17 14 1
17 15 1
18 1 1
18 14 1
18 15 1
19 7 1
19 17 1
19 18 1
20 7 1
20 17 1
20 18 1
21 7 1
21 17 1
21 18 1
22 13 1
22 20 1
22 21 1
23 13 1
23 20 1
23 21 1
24 13 1
24 20 1
24 21 1
5 7 1
12 2 1
12 7 1
[ system ]
simbox
[ molecules ]
tf2N- 100
bpyri+ 100
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