Skip to content

GitLab

Commit e1a4e015 authored by Justin MacCallum's avatar Justin MacCallum
Browse files

Revert "Changed gb settings to use strings rather than named constants." 

This reverts commit 4572d2e3.
parent 40929077
Hide whitespace changes
Inline Side-by-side
Showing with 123 additions and 115 deletions
wrappers/python/simtk/openmm/app/amberprmtopfile.py
View file @ e1a4e015
......@@ -10,7 +10,7 @@ Portions copyright (c) 2012 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
Permission is hereby granted, free of charge, to any person obtaining a
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,
......@@ -40,23 +40,23 @@ import simtk.unit as unit
import simtk.openmm as mm
# Enumerated values for implicit solvent model
# Prefer to use strings now, but these are here for backwards compatibility
HCT = 'HCT'
OBC1 = 'OBC1'
OBC2 = 'OBC2'
GBn = 'GBn'
HCT = object()
OBC1 = object()
OBC2 = object()
GBn = object()
class AmberPrmtopFile(object):
"""AmberPrmtopFile parses an AMBER prmtop file and constructs a Topology and (optionally) an OpenMM System from it."""
def __init__(self, file):
"""Load a prmtop file."""
top = Topology()
## The Topology read from the prmtop file
self.topology = top
# Load the prmtop file
prmtop = amber_file_parser.PrmtopLoader(file)
self._prmtop = prmtop
......@@ -82,7 +82,7 @@ class AmberPrmtopFile(object):
atomName = atomReplacements[atomName]
# Try to guess the element.
upper = atomName.upper()
if upper.startswith('CL'):
element = elem.chlorine
......@@ -96,17 +96,17 @@ class AmberPrmtopFile(object):
except KeyError:
element = None
top.addAtom(atomName, element, r)
# Add bonds to the topology
atoms = list(top.atoms())
for bond in prmtop.getBondsWithH():
top.addBond(atoms[bond[0]], atoms[bond[1]])
for bond in prmtop.getBondsNoH():
top.addBond(atoms[bond[0]], atoms[bond[1]])
# Set the periodic box size.
if prmtop.getIfBox():
top.setUnitCellDimensions(tuple(x.value_in_unit(unit.nanometer) for x in prmtop.getBoxBetaAndDimensions()[1:4])*unit.nanometer)
......@@ -114,7 +114,7 @@ class AmberPrmtopFile(object):
constraints=None, rigidWater=True, implicitSolvent=None, soluteDielectric=1.0, solventDielectric=78.5, removeCMMotion=True,
ewaldErrorTolerance=0.0005):
"""Construct an OpenMM System representing the topology described by this prmtop file.
Parameters:
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
......@@ -122,8 +122,7 @@ class AmberPrmtopFile(object):
- constraints (object=None) Specifies which bonds angles should be implemented with constraints.
Allowed values are None, HBonds, AllBonds, or HAngles.
- rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
- implicitSolvent (object=None) If not None, the implicit solvent
model to use. Allowed values are 'HCT', 'OBC1', 'OBC2', or 'GBn'.
- implicitSolvent (object=None) If not None, the implicit solvent model to use. Allowed values are HCT, OBC1, OBC2, or GBn.
- soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
- solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
- removeCMMotion (boolean=True) If true, a CMMotionRemover will be added to the System
......@@ -149,15 +148,24 @@ class AmberPrmtopFile(object):
constraintString = constraintMap[constraints]
else:
raise ValueError('Illegal value for constraints')
sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString,
nonbondedCutoff=nonbondedCutoff, nonbondedMethod=methodMap[nonbondedMethod],
flexibleConstraints=False, gbmodel=implicitSolvent,
soluteDielectric=soluteDielectric, solventDielectric=solventDielectric,
rigidWater=rigidWater)
if implicitSolvent is None:
implicitString = None
elif implicitSolvent == HCT:
implicitString = 'HCT'
elif implicitSolvent == OBC1:
implicitString = 'OBC1'
elif implicitSolvent == OBC2:
implicitString = 'OBC2'
elif implicitSolvent == GBn:
implicitString = 'GBn'
else:
raise ValueError('Illegal value for implicit solvent model')
sys = amber_file_parser.readAmberSystem(prmtop_loader=self._prmtop, shake=constraintString, nonbondedCutoff=nonbondedCutoff,
nonbondedMethod=methodMap[nonbondedMethod], flexibleConstraints=False, gbmodel=implicitString,
soluteDielectric=soluteDielectric, solventDielectric=solventDielectric, rigidWater=rigidWater)
for force in sys.getForces():
if isinstance(force, mm.NonbondedForce):
force.setEwaldErrorTolerance(ewaldErrorTolerance)
if removeCMMotion:
sys.addForce(mm.CMMotionRemover())
return sys
return sys
\ No newline at end of file
wrappers/python/simtk/openmm/app/gromacstopfile.py
View file @ e1a4e015
......@@ -10,7 +10,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
Permission is hereby granted, free of charge, to any person obtaining a
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,
......@@ -45,11 +45,11 @@ HBonds = ff.HBonds
AllBonds = ff.AllBonds
HAngles = ff.HAngles
OBC2 = 'OBC2'
OBC2 = prmtop.OBC2
class GromacsTopFile(object):
"""GromacsTopFile parses a Gromacs top file and constructs a Topology and (optionally) an OpenMM System from it."""
class _MoleculeType(object):
"""Inner class to store information about a molecule type."""
def __init__(self):
......@@ -60,7 +60,7 @@ class GromacsTopFile(object):
self.exclusions = []
self.pairs = []
self.cmaps = []
def _processFile(self, file):
append = ''
for line in open(file):
......@@ -69,7 +69,7 @@ class GromacsTopFile(object):
else:
self._processLine(append+' '+line, file)
append = ''
def _processLine(self, line, file):
"""Process one line from a file."""
if ';' in line:
......@@ -79,13 +79,13 @@ class GromacsTopFile(object):
if stripped.startswith('*') or len(stripped) == 0:
# A comment or empty line.
return
elif stripped.startswith('[') and not ignore:
# The start of a category.
if not stripped.endswith(']'):
raise ValueError('Illegal line in .top file: '+line)
self._currentCategory = stripped[1:-1].strip()
elif stripped.startswith('#'):
# A preprocessor command.
fields = stripped.split()
......@@ -127,7 +127,7 @@ class GromacsTopFile(object):
if len(self._ifStack) == 0:
raise ValueError('Unexpected line in .top file: '+line)
del(self._ifStack[-1])
elif not ignore:
# A line of data for the current category
if self._currentCategory is None:
......@@ -166,7 +166,7 @@ class GromacsTopFile(object):
self._processPairType(line)
elif self._currentCategory == 'cmaptypes':
self._processCmapType(line)
def _processDefaults(self, line):
"""Process the [ defaults ] line."""
fields = line.split()
......@@ -179,7 +179,7 @@ class GromacsTopFile(object):
if fields[2].lower() == 'no':
raise ValueError('gen_pairs=no is not supported')
self._defaults = fields
def _processMoleculeType(self, line):
"""Process a line in the [ moleculetypes ] category."""
fields = line.split()
......@@ -188,14 +188,14 @@ class GromacsTopFile(object):
type = GromacsTopFile._MoleculeType()
self._moleculeTypes[fields[0]] = type
self._currentMoleculeType = type
def _processMolecule(self, line):
"""Process a line in the [ molecules ] category."""
fields = line.split()
if len(fields) < 2:
raise ValueError('Too few fields in [ molecules ] line: '+line);
self._molecules.append((fields[0], int(fields[1])))
def _processAtom(self, line):
"""Process a line in the [ atoms ] category."""
if self._currentMoleculeType is None:
......@@ -204,7 +204,7 @@ class GromacsTopFile(object):
if len(fields) < 5:
raise ValueError('Too few fields in [ atoms ] line: '+line);
self._currentMoleculeType.atoms.append(fields)
def _processBond(self, line):
"""Process a line in the [ bonds ] category."""
if self._currentMoleculeType is None:
......@@ -215,7 +215,7 @@ class GromacsTopFile(object):
if fields[2] != '1':
raise ValueError('Unsupported function type in [ bonds ] line: '+line);
self._currentMoleculeType.bonds.append(fields)
def _processAngle(self, line):
"""Process a line in the [ angles ] category."""
if self._currentMoleculeType is None:
......@@ -226,7 +226,7 @@ class GromacsTopFile(object):
if fields[3] not in ('1', '5'):
raise ValueError('Unsupported function type in [ angles ] line: '+line);
self._currentMoleculeType.angles.append(fields)
def _processDihedral(self, line):
"""Process a line in the [ dihedrals ] category."""
if self._currentMoleculeType is None:
......@@ -237,7 +237,7 @@ class GromacsTopFile(object):
if fields[4] not in ('1', '2', '3', '4', '9'):
raise ValueError('Unsupported function type in [ dihedrals ] line: '+line);
self._currentMoleculeType.dihedrals.append(fields)
def _processExclusion(self, line):
"""Process a line in the [ exclusions ] category."""
if self._currentMoleculeType is None:
......@@ -246,7 +246,7 @@ class GromacsTopFile(object):
if len(fields) < 2:
raise ValueError('Too few fields in [ exclusions ] line: '+line);
self._currentMoleculeType.exclusions.append(fields)
def _processPair(self, line):
"""Process a line in the [ pairs ] category."""
if self._currentMoleculeType is None:
......@@ -257,7 +257,7 @@ class GromacsTopFile(object):
if fields[2] != '1':
raise ValueError('Unsupported function type in [ pairs ] line: '+line);
self._currentMoleculeType.pairs.append(fields)
def _processCmap(self, line):
"""Process a line in the [ cmaps ] category."""
if self._currentMoleculeType is None:
......@@ -266,14 +266,14 @@ class GromacsTopFile(object):
if len(fields) < 6:
raise ValueError('Too few fields in [ pairs ] line: '+line);
self._currentMoleculeType.cmaps.append(fields)
def _processAtomType(self, line):
"""Process a line in the [ atomtypes ] category."""
fields = line.split()
if len(fields) < 7:
raise ValueError('Too few fields in [ atomtypes ] line: '+line);
self._atomTypes[fields[0]] = fields
def _processBondType(self, line):
"""Process a line in the [ bondtypes ] category."""
fields = line.split()
......@@ -282,7 +282,7 @@ class GromacsTopFile(object):
if fields[2] != '1':
raise ValueError('Unsupported function type in [ bondtypes ] line: '+line);
self._bondTypes[tuple(fields[:2])] = fields
def _processAngleType(self, line):
"""Process a line in the [ angletypes ] category."""
fields = line.split()
......@@ -291,7 +291,7 @@ class GromacsTopFile(object):
if fields[3] not in ('1', '5'):
raise ValueError('Unsupported function type in [ angletypes ] line: '+line);
self._angleTypes[tuple(fields[:3])] = fields
def _processDihedralType(self, line):
"""Process a line in the [ dihedraltypes ] category."""
fields = line.split()
......@@ -305,14 +305,14 @@ class GromacsTopFile(object):
self._dihedralTypes[key].append(fields)
else:
self._dihedralTypes[key] = [fields]
def _processImplicitType(self, line):
"""Process a line in the [ implicit_genborn_params ] category."""
fields = line.split()
if len(fields) < 6:
raise ValueError('Too few fields in [ implicit_genborn_params ] line: '+line);
self._implicitTypes[fields[0]] = fields
def _processPairType(self, line):
"""Process a line in the [ pairtypes ] category."""
fields = line.split()
......@@ -321,7 +321,7 @@ class GromacsTopFile(object):
if fields[2] != '1':
raise ValueError('Unsupported function type in [ pairtypes ] line: '+line);
self._pairTypes[tuple(fields[:2])] = fields
def _processCmapType(self, line):
"""Process a line in the [ cmaptypes ] category."""
fields = line.split()
......@@ -330,10 +330,10 @@ class GromacsTopFile(object):
if fields[5] != '1':
raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line);
self._cmapTypes[tuple(fields[:5])] = fields
def __init__(self, file, unitCellDimensions=None, includeDir='/usr/local/gromacs/share/gromacs/top', defines={}):
"""Load a top file.
Parameters:
- file (string) the name of the file to load
- unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell
......@@ -343,9 +343,9 @@ class GromacsTopFile(object):
"""
self._includeDirs = (os.path.dirname(file), includeDir)
self._defines = defines
# Parse the file.
self._currentCategory = None
self._ifStack = []
self._moleculeTypes = {}
......@@ -359,9 +359,9 @@ class GromacsTopFile(object):
self._pairTypes = {}
self._cmapTypes = {}
self._processFile(file)
# Create the Topology from it.
top = Topology()
## The Topology read from the prmtop file
self.topology = top
......@@ -371,9 +371,9 @@ class GromacsTopFile(object):
if moleculeName not in self._moleculeTypes:
raise ValueError("Unknown molecule type: "+moleculeName)
moleculeType = self._moleculeTypes[moleculeName]
# Create the specified number of molecules of this type.
for i in range(moleculeCount):
atoms = []
lastResidue = None
......@@ -393,9 +393,9 @@ class GromacsTopFile(object):
atomName = fields[4]
if atomName in atomReplacements:
atomName = atomReplacements[atomName]
# Try to guess the element.
upper = atomName.upper()
if upper.startswith('CL'):
element = elem.chlorine
......@@ -409,16 +409,16 @@ class GromacsTopFile(object):
except KeyError:
element = None
atoms.append(top.addAtom(atomName, element, r))
# Add bonds to the topology
for fields in moleculeType.bonds:
top.addBond(atoms[int(fields[0])-1], atoms[int(fields[1])-1])
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):
"""Construct an OpenMM System representing the topology described by this prmtop file.
Parameters:
- nonbondedMethod (object=NoCutoff) The method to use for nonbonded interactions. Allowed values are
NoCutoff, CutoffNonPeriodic, CutoffPeriodic, Ewald, or PME.
......@@ -426,7 +426,7 @@ class GromacsTopFile(object):
- constraints (object=None) Specifies which bonds and angles should be implemented with constraints.
Allowed values are None, HBonds, AllBonds, or HAngles.
- rigidWater (boolean=True) If true, water molecules will be fully rigid regardless of the value passed for the constraints argument
- implicitSolvent (string=None) If not None, the implicit solvent model to use. The only allowed value is 'OBC2'.
- implicitSolvent (object=None) If not None, the implicit solvent model to use. The only allowed value is OBC2.
- soluteDielectric (float=1.0) The solute dielectric constant to use in the implicit solvent model.
- solventDielectric (float=78.5) The solvent dielectric constant to use in the implicit solvent model.
- ewaldErrorTolerance (float=0.0005) The error tolerance to use if nonbondedMethod is Ewald or PME.
......@@ -434,7 +434,7 @@ class GromacsTopFile(object):
Returns: the newly created System
"""
# Create the System.
sys = mm.System()
boxSize = self.topology.getUnitCellDimensions()
if boxSize is not None:
......@@ -443,7 +443,7 @@ class GromacsTopFile(object):
raise ValueError('Illegal nonbonded method for a non-periodic system')
nb = mm.NonbondedForce()
sys.addForce(nb)
if implicitSolvent == 'OBC2':
if implicitSolvent is OBC2:
gb = mm.GBSAOBCForce()
gb.setSoluteDielectric(soluteDielectric)
gb.setSolventDielectric(solventDielectric)
......@@ -461,33 +461,33 @@ class GromacsTopFile(object):
topologyAtoms = list(self.topology.atoms())
exceptions = []
fudgeQQ = float(self._defaults[4])
# Loop over molecules and create the specified number of each type.
for moleculeName, moleculeCount in self._molecules:
moleculeType = self._moleculeTypes[moleculeName]
for i in range(moleculeCount):
# Record the types of all atoms.
baseAtomIndex = sys.getNumParticles()
atomTypes = [atom[1] for atom in moleculeType.atoms]
try:
[self._atomTypes[t][1] for t in atomTypes]
except KeyError as e:
raise ValueError('Unknown atom type: '+e.message)
# Add atoms.
for fields in moleculeType.atoms:
if len(fields) >= 8:
mass = float(fields[7])
else:
mass = float(self._atomTypes[fields[1]][2])
sys.addParticle(mass)
# Add bonds.
atomBonds = [{} for x in range(len(moleculeType.atoms))]
for fields in moleculeType.bonds:
atoms = [int(x)-1 for x in fields[:2]]
......@@ -522,9 +522,9 @@ class GromacsTopFile(object):
# Record information that will be needed for constraining angles.
atomBonds[atoms[0]][atoms[1]] = length
atomBonds[atoms[1]][atoms[0]] = length
# Add angles.
degToRad = math.pi/180
for fields in moleculeType.angles:
atoms = [int(x)-1 for x in fields[:3]]
......@@ -571,7 +571,7 @@ class GromacsTopFile(object):
bonds.addBond(baseAtomIndex+atoms[0], baseAtomIndex+atoms[2], float(params[2]), k)
# Add torsions.
for fields in moleculeType.dihedrals:
atoms = [int(x)-1 for x in fields[:4]]
types = tuple(atomTypes[i] for i in atoms)
......@@ -617,7 +617,7 @@ class GromacsTopFile(object):
rb = mm.RBTorsionForce()
sys.addForce(rb)
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.
for fields in moleculeType.cmaps:
......@@ -648,7 +648,7 @@ class GromacsTopFile(object):
baseAtomIndex+atoms[1], baseAtomIndex+atoms[2], baseAtomIndex+atoms[3], baseAtomIndex+atoms[4])
# Set nonbonded parameters for particles.
for fields in moleculeType.atoms:
params = self._atomTypes[fields[1]]
if len(fields) > 6:
......@@ -656,7 +656,7 @@ class GromacsTopFile(object):
else:
q = float(params[3])
nb.addParticle(q, float(params[5]), float(params[6]))
if implicitSolvent == 'OBC2':
if implicitSolvent is OBC2:
if fields[1] not in self._implicitTypes:
raise ValueError('No implicit solvent parameters specified for atom type: '+fields[1])
gbparams = self._implicitTypes[fields[1]]
......@@ -664,9 +664,9 @@ class GromacsTopFile(object):
for fields in moleculeType.bonds:
atoms = [int(x)-1 for x in fields[:2]]
bondIndices.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1]))
# Record nonbonded exceptions.
for fields in moleculeType.pairs:
atoms = [int(x)-1 for x in fields[:2]]
types = tuple(atomTypes[i] for i in atoms)
......@@ -681,15 +681,15 @@ class GromacsTopFile(object):
atom1params = nb.getParticleParameters(baseAtomIndex+atoms[0])
atom2params = nb.getParticleParameters(baseAtomIndex+atoms[1])
exceptions.append((baseAtomIndex+atoms[0], baseAtomIndex+atoms[1], atom1params[0]*atom2params[0]*fudgeQQ, params[0], params[1]))
# Create nonbonded exceptions.
nb.createExceptionsFromBonds(bondIndices, fudgeQQ, float(self._defaults[3]))
for exception in exceptions:
nb.addException(exception[0], exception[1], exception[2], float(exception[3]), float(exception[4]), True)
# Finish configuring the NonbondedForce.
methodMap = {ff.NoCutoff:mm.NonbondedForce.NoCutoff,
ff.CutoffNonPeriodic:mm.NonbondedForce.CutoffNonPeriodic,
ff.CutoffPeriodic:mm.NonbondedForce.CutoffPeriodic,
......
wrappers/python/simtk/openmm/app/internal/amber_file_parser.py
View file @ e1a4e015
......@@ -16,7 +16,7 @@ Portions copyright (c) 2012-2013 Stanford University and the Authors.
Authors: Randall J. Radmer, John D. Chodera, Peter Eastman
Contributors: Christoph Klein, Michael R. Shirts
Permission is hereby granted, free of charge, to any person obtaining a
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,
......@@ -92,7 +92,7 @@ class PrmtopLoader(object):
>>> import os, os.path
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
>>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
>>> prmtop = PrmtopLoader(prmtop_filename)
>>> prmtop = PrmtopLoader(prmtop_filename)
"""
def __init__(self, inFilename):
......@@ -101,7 +101,7 @@ class PrmtopLoader(object):
ARGUMENTS
inFilename (string) - AMBER 'new-style' prmtop file, probably generated with one of the AMBER tleap/xleap/sleap
inFilename (string) - AMBER 'new-style' prmtop file, probably generated with one of the AMBER tleap/xleap/sleap
"""
......@@ -152,7 +152,7 @@ class PrmtopLoader(object):
Parameter:
- pointerLabel: a string matching one of the following:
NATOM : total number of atoms
NATOM : total number of atoms
NTYPES : total number of distinct atom types
NBONH : number of bonds containing hydrogen
MBONA : number of bonds not containing hydrogen
......@@ -183,7 +183,7 @@ class PrmtopLoader(object):
NMXRS : number of atoms in the largest residue
IFCAP : set to 1 if the CAP option from edit was specified
"""
index = POINTER_LABEL_LIST.index(pointerLabel)
index = POINTER_LABEL_LIST.index(pointerLabel)
return float(self._raw_data['POINTERS'][index])
def getNumAtoms(self):
......@@ -350,7 +350,7 @@ class PrmtopLoader(object):
bondPointers=self._raw_data["BONDS_INC_HYDROGEN"]
self._bondListWithH = self._getBonds(bondPointers)
return self._bondListWithH
def getBondsNoH(self):
"""Return list of bonded atom pairs, K, and Rmin for each bond with no hydrogen"""
......@@ -509,11 +509,11 @@ class PrmtopLoader(object):
def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmodel=None, soluteDielectric=1.0, solventDielectric=78.5, nonbondedCutoff=None, nonbondedMethod='NoCutoff', scee=1.2, scnb=2.0, mm=None, verbose=False, EwaldErrorTolerance=None, flexibleConstraints=True, rigidWater=True):
"""
Create an OpenMM System from an Amber prmtop file.
ARGUMENTS (specify one or the other, but not both)
prmtop_filename (String) - name of Amber prmtop file (new-style only)
prmtop_loader (PrmtopLoader) - the loaded prmtop file
OPTIONAL ARGUMENTS
shake (String) - if 'h-bonds', will SHAKE all bonds to hydrogen and water; if 'all-bonds', will SHAKE all bonds and water (default: None)
gbmodel (String) - if 'OBC', OBC GBSA will be used; if 'GBVI', GB/VI will be used (default: None)
......@@ -547,10 +547,10 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
>>> prmtop_filename = os.path.join(directory, 'alanine-dipeptide.prmtop')
>>> system = readAmberSystem(prmtop_filename)
>>> system = readAmberSystem(prmtop_filename)
"""
if prmtop_filename is None and prmtop_loader is None:
raise Exception("Must specify a filename or loader")
if prmtop_filename is not None and prmtop_loader is not None:
......@@ -567,7 +567,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
if prmtop.getIfPert()>0:
raise Exception("perturbation not currently supported")
if prmtop.getIfBox()>1:
raise Exception("only standard periodic boxes are currently supported")
......@@ -596,9 +596,9 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
if isWater[iAtom] and isWater[jAtom]:
system.addConstraint(iAtom, jAtom, rMin)
# Add harmonic bonds.
if verbose: print "Adding bonds..."
if verbose: print "Adding bonds..."
force = mm.HarmonicBondForce()
if flexibleConstraints or (shake not in ('h-bonds', 'all-bonds', 'h-angles')):
for (iAtom, jAtom, k, rMin) in prmtop.getBondsWithH():
......@@ -610,7 +610,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
system.addForce(force)
# Add harmonic angles.
if verbose: print "Adding angles..."
if verbose: print "Adding angles..."
force = mm.HarmonicAngleForce()
if shake == 'h-angles':
numConstrainedBonds = system.getNumConstraints()
......@@ -638,7 +638,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
l1 = bond[1]
elif bond[0] == kAtom:
l2 = bond[1]
# Compute the distance between atoms and add a constraint
length = math.sqrt(l1*l1 + l2*l2 - 2*l1*l2*math.cos(aMin))
system.addConstraint(iAtom, kAtom, length)
......@@ -647,14 +647,14 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
system.addForce(force)
# Add torsions.
if verbose: print "Adding torsions..."
if verbose: print "Adding torsions..."
force = mm.PeriodicTorsionForce()
for (iAtom, jAtom, kAtom, lAtom, forceConstant, phase, periodicity) in prmtop.getDihedrals():
force.addTorsion(iAtom, jAtom, kAtom, lAtom, periodicity, phase, forceConstant)
system.addForce(force)
# Add nonbonded interactions.
if verbose: print "Adding nonbonded interactions..."
if verbose: print "Adding nonbonded interactions..."
force = mm.NonbondedForce()
if (prmtop.getIfBox() == 0):
# System is non-periodic.
......@@ -663,12 +663,12 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
elif nonbondedMethod == 'CutoffNonPeriodic':
if nonbondedCutoff is None:
raise Exception("No cutoff value specified")
force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
force.setNonbondedMethod(mm.NonbondedForce.CutoffNonPeriodic)
force.setCutoffDistance(nonbondedCutoff)
else:
raise Exception("Illegal nonbonded method for a non-periodic system")
else:
# System is periodic.
# System is periodic.
# Set periodic box vectors for periodic system
(boxBeta, boxX, boxY, boxZ) = prmtop.getBoxBetaAndDimensions()
d0 = units.Quantity(0.0, units.angstroms)
......@@ -676,16 +676,16 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
yVec = units.Quantity((d0, boxY, d0))
zVec = units.Quantity((d0, d0, boxZ))
system.setDefaultPeriodicBoxVectors(xVec, yVec, zVec)
# Set cutoff.
if nonbondedCutoff is None:
# Compute cutoff automatically.
min_box_width = min([boxX / units.nanometers, boxY / units.nanometers, boxZ / units.nanometers])
CLEARANCE_FACTOR = 0.97 # reduce the cutoff to be a bit smaller than 1/2 smallest box length
CLEARANCE_FACTOR = 0.97 # reduce the cutoff to be a bit smaller than 1/2 smallest box length
nonbondedCutoff = units.Quantity((min_box_width * CLEARANCE_FACTOR) / 2.0, units.nanometers)
if nonbondedMethod != 'NoCutoff':
force.setCutoffDistance(nonbondedCutoff)
# Set nonbonded method.
if nonbondedMethod == 'NoCutoff':
force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
......@@ -724,7 +724,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
excludeParams = (0.0, 0.1, 0.0)
for iAtom in range(prmtop.getNumAtoms()):
for jAtom in excludedAtoms[iAtom]:
if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue
if min((iAtom, jAtom), (jAtom, iAtom)) in excludedAtomPairs: continue
force.addException(iAtom, jAtom, excludeParams[0], excludeParams[1], excludeParams[2])
system.addForce(force)
......@@ -812,7 +812,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
if gbmodel == 'OBC2':
gb.addParticle(charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1])
else:
gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]])
gb.addParticle([charges[iAtom], gb_parms[iAtom][0], gb_parms[iAtom][1]])
system.addForce(gb)
if nonbondedMethod == 'NoCutoff':
gb.setNonbondedMethod(mm.NonbondedForce.NoCutoff)
......@@ -835,7 +835,7 @@ def readAmberSystem(prmtop_filename=None, prmtop_loader=None, shake=None, gbmode
def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbose=False, asNumpy=False):
"""
Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file.
Read atomic coordinates (and optionally, box vectors) from Amber formatted coordinate file.
ARGUMENTS
......@@ -852,13 +852,13 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
Read coordinates in vacuum.
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-gbsa')
>>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
>>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
>>> coordinates = readAmberCoordinates(crd_filename)
Read coordinates in solvent.
>>> directory = os.path.join(os.getenv('YANK_INSTALL_DIR'), 'test', 'systems', 'alanine-dipeptide-explicit')
>>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
>>> crd_filename = os.path.join(directory, 'alanine-dipeptide.inpcrd')
>>> [coordinates, box_vectors] = readAmberCoordinates(crd_filename, read_box=True)
"""
......@@ -921,7 +921,7 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
velocities = newvel
# Assign units.
velocities = units.Quantity(velocities, units.angstroms/units.picoseconds)
# Read box size if present
box_vectors = None
if (read_box):
......@@ -943,20 +943,20 @@ def readAmberCoordinates(filename, read_box=False, read_velocities=False, verbos
else:
a = units.Quantity(mm.Vec3(box_dimensions[0], 0.0, 0.0), units.angstroms)
b = units.Quantity(mm.Vec3(0.0, box_dimensions[1], 0.0), units.angstroms)
c = units.Quantity(mm.Vec3(0.0, 0.0, box_dimensions[2]), units.angstroms)
c = units.Quantity(mm.Vec3(0.0, 0.0, box_dimensions[2]), units.angstroms)
box_vectors = [a,b,c]
else:
raise Exception("Don't know what to do with box vectors: %s" % line)
# Close file
infile.close()
if box_vectors and velocities:
return (coordinates, box_vectors, velocities)
if box_vectors:
return (coordinates, box_vectors)
if velocities:
return (coordinates, velocities)
return (coordinates, velocities)
return coordinates
#=============================================================================================
......
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