Commit 8bc8d184 authored by peastman's avatar peastman
Browse files

Merge pull request #803 from peastman/triclinic

Python layer supports triclinic boxes
parents abadfd9d 30a53702
......@@ -572,7 +572,7 @@ with the name :file:`simulateGromacs.py`.
from sys import stdout
gro = GromacsGroFile('input.gro')
top = GromacsTopFile('input.top', unitCellDimensions=gro.getUnitCellDimensions(),
top = GromacsTopFile('input.top', periodicBoxVectors=gro.getPeriodicBoxVectors(),
includeDir='/usr/local/gromacs/share/gromacs/top')
system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer,
constraints=HBonds)
......@@ -593,10 +593,10 @@ This script is nearly identical to the previous one, just replacing
:class:`AmberInpcrdFile` and :class:`AmberPrmtopFile` with :class:`GromacsGroFile` and :class:`GromacsTopFile`.
Note that when we create the :class:`GromacsTopFile`, we specify values for two extra
options. First, we specify
:code:`unitCellDimensions=gro.getUnitCellDimensions()`\ . Unlike OpenMM and
AMBER, which can store periodic unit cell dimensions with the topology, Gromacs
only stores them with the coordinates. To let :class:`GromacsTopFile` create a :class:`Topology`
object, we therefore need to tell it the unit cell dimensions that were loaded
:code:`periodicBoxVectors=gro.getPeriodicBoxVectors()`\ . Unlike OpenMM and
AMBER, which can store periodic unit cell information with the topology, Gromacs
only stores it with the coordinates. To let :class:`GromacsTopFile` create a :class:`Topology`
object, we therefore need to tell it the periodic box vectors that were loaded
from the :file:`gro` file. You only need to do this if you are simulating a periodic
system. For implicit solvent simulations, it usually can be omitted.
......
......@@ -4,7 +4,7 @@ from simtk.unit import *
from sys import stdout
gro = GromacsGroFile('input.gro')
top = GromacsTopFile('input.top', unitCellDimensions=gro.getUnitCellDimensions(), includeDir='/usr/local/gromacs/share/gromacs/top')
top = GromacsTopFile('input.top', periodicBoxVectors=gro.getPeriodicBoxVectors(), includeDir='/usr/local/gromacs/share/gromacs/top')
system = top.createSystem(nonbondedMethod=PME, nonbondedCutoff=1*nanometer, constraints=HBonds)
integrator = LangevinIntegrator(300*kelvin, 1/picosecond, 0.002*picoseconds)
simulation = Simulation(top.topology, system, integrator)
......
......@@ -54,10 +54,12 @@ foreach(SUBDIR ${SUBDIRS})
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.sh"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.xml"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdb"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.pdbx"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prmtop"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.prm"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.inpcrd"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.crd"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.gro"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.parm7"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.rst7"
"${CMAKE_CURRENT_SOURCE_DIR}/${SUBDIR}/*.ncrst"
......
......@@ -77,14 +77,17 @@ class DCDFile(object):
header += struct.pack('<4i', 164, 4, len(list(topology.atoms())), 4)
file.write(header)
def writeModel(self, positions, unitCellDimensions=None):
def writeModel(self, positions, unitCellDimensions=None, periodicBoxVectors=None):
"""Write out a model to the DCD file.
The periodic box can be specified either by the unit cell dimensions (for a rectangular box), or the full set of box
vectors (for an arbitrary triclinic box). If neither is specified, the box vectors specified in the Topology will be
used. Regardless of the value specified, no dimensions will be written if the Topology does not represent a periodic system.
Parameters:
- positions (list) The list of atomic positions to write
- unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell. If None, the dimensions specified in
the Topology will be used. Regardless of the value specified, no dimensions will be written if the Topology does not
represent a periodic system.
- unitCellDimensions (Vec3=None) The dimensions of the crystallographic unit cell.
- periodicBoxVectors (tuple of Vec3=None) The vectors defining the periodic box.
"""
if len(list(self._topology.atoms())) != len(positions):
raise ValueError('The number of positions must match the number of atoms')
......@@ -107,12 +110,22 @@ class DCDFile(object):
# Write the data.
file.seek(0, os.SEEK_END)
boxSize = self._topology.getUnitCellDimensions()
if boxSize is not None:
if unitCellDimensions is not None:
boxSize = unitCellDimensions
size = boxSize.value_in_unit(angstroms)
file.write(struct.pack('<i6di', 48, size[0], 0, size[1], 0, 0, size[2], 48))
boxVectors = self._topology.getPeriodicBoxVectors()
if boxVectors is not None:
if getPeriodicBoxVectors is not None:
boxVectors = getPeriodicBoxVectors
elif unitCellDimensions is not None:
if is_quantity(unitCellDimensions):
unitCellDimensions = unitCellDimensions.value_in_unit(nanometers)
boxVectors = (Vec3(unitCellDimensions[0], 0, 0), Vec3(0, unitCellDimensions[1], 0), Vec3(0, 0, unitCellDimensions[2]))*nanometers
(a_length, b_length, c_length, alpha, beta, gamma) = computeLengthsAndAngles(boxVectors)
a_length = a_length.value_in_unit(angstroms)
b_length = b_length.value_in_unit(angstroms)
c_length = c_length.value_in_unit(angstroms)
angle1 = math.sin(math.pi/2-gamma)
angle2 = math.sin(math.pi/2-beta)
angle3 = math.sin(math.pi/2-alpha)
file.write(struct.pack('<i6di', 48, a_length, angle1, b_length, angle2, angle3, c_length, 48))
length = struct.pack('<i', 4*len(positions))
for i in range(3):
file.write(length)
......
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2014 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Peter Eastman, Mark Friedrichs
Contributors:
......@@ -403,9 +403,9 @@ class ForceField(object):
# Set periodic boundary conditions.
boxSize = topology.getUnitCellDimensions()
if boxSize is not None:
sys.setDefaultPeriodicBoxVectors((boxSize[0], 0, 0), (0, boxSize[1], 0), (0, 0, boxSize[2]))
boxVectors = topology.getPeriodicBoxVectors()
if boxVectors is not None:
sys.setDefaultPeriodicBoxVectors(boxVectors[0], boxVectors[1], boxVectors[2])
elif nonbondedMethod is not NoCutoff and nonbondedMethod is not CutoffNonPeriodic:
raise ValueError('Requested periodic boundary conditions for a Topology that does not specify periodic box dimensions')
......
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Lee-Ping Wang, Peter Eastman
Contributors:
......@@ -89,6 +89,18 @@ def _is_gro_box(line):
else:
return 0
def _construct_box_vectors(line):
"""Create the periodic box vectors based on the values stored in the file.
@param[in] line The line containing the description
"""
sline = line.split()
values = [float(i) for i in sline]
if len(sline) == 3:
return (Vec3(values[0], 0, 0), Vec3(0, values[1], 0), Vec3(0, 0, values[2]))*nanometers
return (Vec3(values[0], values[3], values[4]), Vec3(values[5], values[1], values[6]), Vec3(values[7], values[8], values[2]))*nanometers
class GromacsGroFile(object):
"""GromacsGroFile parses a Gromacs .gro file and constructs a set of atom positions from it.
......@@ -140,7 +152,7 @@ class GromacsGroFile(object):
xyz.append(Vec3(pos[0], pos[1], pos[2]))
elif _is_gro_box(line) and ln == na + 2:
sline = line.split()
boxes.append(tuple([float(i) for i in sline])*nanometers)
boxes.append(_construct_box_vectors(line))
xyzs.append(xyz*nanometers)
xyz = []
ln = -1
......@@ -160,7 +172,7 @@ class GromacsGroFile(object):
## A list containing the name of the residue that each atom belongs to
self.residueNames = resname
self._positions = xyzs
self._unitCellDimensions = boxes
self._periodicBoxVectors = boxes
self._numpyPositions = None
def getNumFrames(self):
......@@ -182,10 +194,21 @@ class GromacsGroFile(object):
return self._numpyPositions[frame]
return self._positions[frame]
def getPeriodicBoxVectors(self, frame=0):
"""Get the vectors defining the periodic box.
Parameters:
- frame (int=0) the index of the frame for which to get the box vectors
"""
return self._periodicBoxVectors[frame]
def getUnitCellDimensions(self, frame=0):
"""Get the dimensions of the crystallographic unit cell.
Parameters:
- frame (int=0) the index of the frame for which to get the unit cell dimensions
"""
return self._unitCellDimensions[frame]
xsize = self._periodicBoxVectors[frame][0][0].value_in_unit(nanometers)
ysize = self._periodicBoxVectors[frame][1][1].value_in_unit(nanometers)
zsize = self._periodicBoxVectors[frame][2][2].value_in_unit(nanometers)
return Vec3(xsize, ysize, zsize)*nanometers
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2014 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -358,12 +358,14 @@ class GromacsTopFile(object):
raise ValueError('Unsupported function type in [ cmaptypes ] line: '+line);
self._cmapTypes[tuple(fields[:5])] = fields
def __init__(self, file, unitCellDimensions=None, includeDir=None, defines=None):
def __init__(self, file, periodicBoxVectors=None, unitCellDimensions=None, includeDir=None, defines=None):
"""Load a top file.
Parameters:
- file (string) the name of the file to load
- unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell
- periodicBoxVectors (tuple of Vec3=None) the vectors defining the periodic box
- unitCellDimensions (Vec3=None) the dimensions of the crystallographic unit cell. For
non-rectangular unit cells, specify periodicBoxVectors instead.
- includeDir (string=None) A directory in which to look for other files
included from the top file. If not specified, we will attempt to locate a gromacs
installation on your system. When gromacs is installed in /usr/local, this will resolve
......@@ -403,6 +405,11 @@ class GromacsTopFile(object):
top = Topology()
## The Topology read from the prmtop file
self.topology = top
if periodicBoxVectors is not None:
if unitCellDimensions is not None:
raise ValueError("specify either periodicBoxVectors or unitCellDimensions, but not both")
top.setPeriodicBoxVectors(periodicBoxVectors)
else:
top.setUnitCellDimensions(unitCellDimensions)
PDBFile._loadNameReplacementTables()
for moleculeName, moleculeCount in self._molecules:
......
......@@ -8,7 +8,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2013 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Christopher M. Bruns
Contributors: Peter Eastman
......@@ -37,8 +37,10 @@ __version__ = "1.0"
from simtk.openmm.vec3 import Vec3
import simtk.unit as unit
from .. import element
from unitcell import computePeriodicBoxVectors
import warnings
import sys
import math
class PdbStructure(object):
"""
......@@ -137,7 +139,7 @@ class PdbStructure(object):
self._current_model = None
self.default_model = None
self.models_by_number = {}
self._unit_cell_dimensions = None
self._periodic_box_vectors = None
self.sequences = []
self.modified_residues = []
# read file
......@@ -170,7 +172,13 @@ class PdbStructure(object):
self._current_model._current_chain._add_ter_record()
self._reset_residue_numbers()
elif (pdb_line.find("CRYST1") == 0):
self._unit_cell_dimensions = Vec3(float(pdb_line[6:15]), float(pdb_line[15:24]), float(pdb_line[24:33]))*unit.angstroms
a_length = float(pdb_line[6:15])*0.1
b_length = float(pdb_line[15:24])*0.1
c_length = float(pdb_line[24:33])*0.1
alpha = float(pdb_line[33:40])*math.pi/180.0
beta = float(pdb_line[40:47])*math.pi/180.0
gamma = float(pdb_line[47:54])*math.pi/180.0
self._periodic_box_vectors = computePeriodicBoxVectors(a_length, b_length, c_length, alpha, beta, gamma)
elif (pdb_line.find("CONECT") == 0):
atoms = [int(pdb_line[6:11])]
for pos in (11,16,21,26):
......@@ -283,9 +291,9 @@ class PdbStructure(object):
for model in self.models:
model._finalize()
def get_unit_cell_dimensions(self):
"""Get the dimensions of the crystallographic unit cell (may be None)."""
return self._unit_cell_dimensions
def get_periodic_box_vectors(self):
"""Get the vectors defining the crystallographic unit cell (may be None)."""
return self._periodic_box_vectors
class Sequence(object):
......
"""
unitcell.py: Routines for converting between different representations of the periodic unit cell.
This is part of the OpenMM molecular simulation toolkit originating from
Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2015 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
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,
and/or sell copies of the Software, and to permit persons to whom the
Software is furnished to do so, subject to the following conditions:
The above copyright notice and this permission notice shall be included in
all copies or substantial portions of the Software.
THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,
DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR
OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE
USE OR OTHER DEALINGS IN THE SOFTWARE.
"""
__author__ = "Peter Eastman"
__version__ = "1.0"
from simtk.openmm import Vec3
from simtk.unit import nanometers, is_quantity, norm, dot
import math
def computePeriodicBoxVectors(a_length, b_length, c_length, alpha, beta, gamma):
"""Convert lengths and angles to periodic box vectors.
Lengths should be given in nanometers and angles in radians.
"""
# Compute the vectors.
a = [a_length, 0, 0]
b = [b_length*math.cos(gamma), b_length*math.sin(gamma), 0]
cx = c_length*math.cos(beta)
cy = c_length*(math.cos(alpha)-math.cos(beta)*math.cos(gamma))/math.sin(gamma)
cz = math.sqrt(c_length*c_length-cx*cx-cy*cy)
c = [cx, cy, cz]
# If any elements are very close to 0, set them to exactly 0.
for i in range(3):
if abs(a[i]) < 1e-6:
a[i] = 0.0
if abs(b[i]) < 1e-6:
b[i] = 0.0
if abs(c[i]) < 1e-6:
c[i] = 0.0
a = Vec3(*a)
b = Vec3(*b)
c = Vec3(*c)
# Make sure they're in the reduced form required by OpenMM.
c = c - b*round(c[1]/b[1])
c = c - a*round(c[0]/a[0])
b = b - a*round(b[0]/a[0])
return (a, b, c)*nanometers
def computeLengthsAndAngles(periodicBoxVectors):
"""Convert periodic box vectors to lengths and angles.
Lengths are returned in nanometers and angles in radians.
"""
(a, b, c) = vectors.value_in_unit(nanometers)
a_length = norm(a)
b_length = norm(b)
c_length = norm(c)
alpha = math.acos(dot(b, c)/(b_length*c_length))
beta = math.acos(dot(c, a)/(c_length*a_length))
gamma = math.acos(dot(a, b)/(a_length*b_length))
return (a_length, b_length, c_length, alpha, beta, gamma)
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2013 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -94,7 +94,7 @@ class Modeller(object):
# Copy over the existing model.
newTopology = Topology()
newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
newTopology.setPeriodicBoxVectors(self.topology.getPeriodicBoxVectors())
newAtoms = {}
newPositions = []*nanometer
for chain in self.topology.chains():
......@@ -140,7 +140,7 @@ class Modeller(object):
- toDelete (list) a list of Atoms, Residues, Chains, and bonds (specified as tuples of Atoms) to delete
"""
newTopology = Topology()
newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
newTopology.setPeriodicBoxVectors(self.topology.getPeriodicBoxVectors())
newAtoms = {}
newPositions = []*nanometer
deleteSet = set(toDelete)
......@@ -189,7 +189,7 @@ class Modeller(object):
else:
raise ValueError('Unknown water model: %s' % model)
newTopology = Topology()
newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
newTopology.setPeriodicBoxVectors(self.topology.getPeriodicBoxVectors())
newAtoms = {}
newPositions = []*nanometer
for chain in self.topology.chains():
......@@ -240,7 +240,7 @@ class Modeller(object):
self.topology = newTopology
self.positions = newPositions
def addSolvent(self, forcefield, model='tip3p', boxSize=None, padding=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar):
def addSolvent(self, forcefield, model='tip3p', boxSize=None, boxVectors=None, padding=None, positiveIon='Na+', negativeIon='Cl-', ionicStrength=0*molar):
"""Add solvent (both water and ions) to the model to fill a rectangular box.
The algorithm works as follows:
......@@ -250,15 +250,17 @@ class Modeller(object):
randomly selecting a water molecule and replacing it with the ion.
4. Ion pairs are added to give the requested total ionic strength.
The box size can be specified in three ways. First, you can explicitly give a box size to use. Alternatively, you can
The box size can be specified in four ways. First, you can explicitly give the vectors defining the periodic box to
use. Alternatively, for a rectangular box you can simply give the dimensions of the unit cell. Third, you can
give a padding distance. The largest dimension of the solute (along the x, y, or z axis) is determined, and a cubic
box of size (largest dimension)+2*padding is used. Finally, if neither a box size nor a padding distance is specified,
the existing Topology's unit cell dimensions are used.
box of size (largest dimension)+2*padding is used. Finally, if neither box vectors, box size, nor padding distance is specified,
the existing Topology's box vectors are used.
Parameters:
- forcefield (ForceField) the ForceField to use for determining van der Waals radii and atomic charges
- model (string='tip3p') the water model to use. Supported values are 'tip3p', 'spce', 'tip4pew', and 'tip5p'.
- boxSize (Vec3=None) the size of the box to fill with water
- boxVectors (tuple of Vec3=None) the vectors defining the periodic box to fill with water
- padding (distance=None) the padding distance to use
- positiveIon (string='Na+') the type of positive ion to add. Allowed values are 'Cs+', 'K+', 'Li+', 'Na+', and 'Rb+'
- negativeIon (string='Cl-') the type of negative ion to add. Allowed values are 'Cl-', 'Br-', 'F-', and 'I-'. Be aware
......@@ -268,18 +270,28 @@ class Modeller(object):
"""
# Pick a unit cell size.
if boxSize is not None:
if boxVectors is not None:
if is_quantity(boxVectors[0]):
boxVectors = (boxVectors[0].value_in_unit(nanometer), boxVectors[1].value_in_unit(nanometer), boxVectors[2].value_in_unit(nanometer))
box = Vec3(boxVectors[0][0], boxVectors[1][1], boxVectors[2][2])
vectors = boxVectors
elif boxSize is not None:
if is_quantity(boxSize):
boxSize = boxSize.value_in_unit(nanometer)
box = Vec3(boxSize[0], boxSize[1], boxSize[2])*nanometer
box = Vec3(boxSize[0], boxSize[1], boxSize[2])
vectors = (Vec3(boxSize[0], 0, 0), Vec3(0, boxSize[1], 0), Vec3(0, 0, boxSize[2]))
elif padding is not None:
if is_quantity(padding):
padding = padding.value_in_unit(nanometer)
maxSize = max(max((pos[i] for pos in self.positions))-min((pos[i] for pos in self.positions)) for i in range(3))
maxSize = maxSize.value_in_unit(nanometer)
box = (maxSize+2*padding)*Vec3(1, 1, 1)
vectors = (Vec3(maxSize+2*padding, 0, 0), Vec3(0, maxSize+2*padding, 0), Vec3(0, 0, maxSize+2*padding))
else:
box = self.topology.getUnitCellDimensions()
box = self.topology.getUnitCellDimensions().value_in_unit(nanometer)
vectors = self.topology.getPeriodicBoxVectors().value_in_unit(nanometer)
if box is None:
raise ValueError('Neither the box size nor padding was specified, and the Topology does not define unit cell dimensions')
box = box.value_in_unit(nanometer)
raise ValueError('Neither the box size, box vectors, nor padding was specified, and the Topology does not define unit cell dimensions')
invBox = Vec3(1.0/box[0], 1.0/box[1], 1.0/box[2])
# Identify the ion types.
......@@ -331,7 +343,7 @@ class Modeller(object):
# Copy the solute over.
newTopology = Topology()
newTopology.setUnitCellDimensions(box)
newTopology.setPeriodicBoxVectors(vectors*nanometer)
newAtoms = {}
newPositions = []*nanometer
for chain in self.topology.chains():
......@@ -378,7 +390,9 @@ class Modeller(object):
def periodicDistance(pos1, pos2):
delta = pos1-pos2
delta = [delta[i]-floor(delta[i]*invBox[i]+0.5)*box[i] for i in range(3)]
delta -= vectors[2]*floor(delta[2]*invBox[2]+0.5)
delta -= vectors[1]*floor(delta[1]*invBox[1]+0.5)
delta -= vectors[0]*floor(delta[0]*invBox[0]+0.5)
return norm(delta)
# Find the list of water molecules to add.
......@@ -472,7 +486,6 @@ class Modeller(object):
for atom2 in molAtoms:
if atom2.element == elem.hydrogen:
newTopology.addBond(atom1, atom2)
newTopology.setUnitCellDimensions(deepcopy(box)*nanometer)
self.topology = newTopology
self.positions = newPositions
......@@ -610,7 +623,7 @@ class Modeller(object):
# Loop over residues.
newTopology = Topology()
newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
newTopology.setPeriodicBoxVectors(self.topology.getPeriodicBoxVectors())
newAtoms = {}
newPositions = []*nanometer
newIndices = []
......@@ -883,7 +896,7 @@ class Modeller(object):
# Create the new Topology.
newTopology = Topology()
newTopology.setUnitCellDimensions(deepcopy(self.topology.getUnitCellDimensions()))
newTopology.setPeriodicBoxVectors(self.topology.getPeriodicBoxVectors())
newAtoms = {}
newPositions = []*nanometer
for chain in self.topology.chains():
......
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -40,7 +40,7 @@ from datetime import date
from simtk.openmm import Vec3, Platform
from simtk.openmm.app.internal.pdbstructure import PdbStructure
from simtk.openmm.app import Topology
from simtk.unit import nanometers, angstroms, is_quantity, norm, Quantity
from simtk.unit import nanometers, angstroms, is_quantity, norm, Quantity, dot
import element as elem
try:
import numpy
......@@ -142,7 +142,7 @@ class PDBFile(object):
self._positions.append(coords*nanometers)
## The atom positions read from the PDB file. If the file contains multiple frames, these are the positions in the first frame.
self.positions = self._positions[0]
self.topology.setUnitCellDimensions(pdb.get_unit_cell_dimensions())
self.topology.setPeriodicBoxVectors(pdb.get_periodic_box_vectors())
self.topology.createStandardBonds()
self.topology.createDisulfideBonds(self.positions)
self._numpyPositions = None
......@@ -250,10 +250,16 @@ class PDBFile(object):
- file (file=stdout) A file to write the file to
"""
print >>file, "REMARK 1 CREATED WITH OPENMM %s, %s" % (Platform.getOpenMMVersion(), str(date.today()))
boxSize = topology.getUnitCellDimensions()
if boxSize is not None:
size = boxSize.value_in_unit(angstroms)
print >>file, "CRYST1%9.3f%9.3f%9.3f 90.00 90.00 90.00 P 1 1 " % size
vectors = topology.getPeriodicBoxVectors()
if vectors is not None:
(a, b, c) = vectors.value_in_unit(angstroms)
a_length = norm(a)
b_length = norm(b)
c_length = norm(c)
alpha = math.acos(dot(b, c)/(b_length*c_length))*180.0/math.pi
beta = math.acos(dot(c, a)/(c_length*a_length))*180.0/math.pi
gamma = math.acos(dot(a, b)/(a_length*b_length))*180.0/math.pi
print >>file, "CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f P 1 1 " % (a_length, b_length, c_length, alpha, beta, gamma)
@staticmethod
def writeModel(topology, positions, file=sys.stdout, modelIndex=None):
......
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2014 Stanford University and the Authors.
Portions copyright (c) 2014-2015 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -33,8 +33,10 @@ __version__ = "1.0"
import os
import sys
import math
from simtk.openmm import Vec3
from simtk.openmm.app.internal.pdbx.reader.PdbxReader import PdbxReader
from simtk.openmm.app.internal.unitcell import computePeriodicBoxVectors
from simtk.openmm.app import Topology
from simtk.unit import nanometers, angstroms, is_quantity, norm, Quantity
import element as elem
......@@ -145,8 +147,9 @@ class PDBxFile(object):
cell = block.getObj('cell')
if cell is not None and cell.getRowCount() > 0:
row = cell.getRow(0)
cellSize = [float(row[cell.getAttributeIndex(attribute)]) for attribute in ('length_a', 'length_b', 'length_c')]*angstroms
self.topology.setUnitCellDimensions(cellSize)
(a, b, c) = [float(row[cell.getAttributeIndex(attribute)])*0.1 for attribute in ('length_a', 'length_b', 'length_c')]
(alpha, beta, gamma) = [float(row[cell.getAttributeIndex(attribute)])*math.pi/180.0 for attribute in ('angle_alpha', 'angle_beta', 'angle_gamma')]
self.topology.setPeriodicBoxVectors(computePeriodicBoxVectors(a, b, c, alpha, beta, gamma))
# Add bonds based on struct_conn records.
......
......@@ -6,7 +6,7 @@ Simbios, the NIH National Center for Physics-Based Simulation of
Biological Structures at Stanford, funded under the NIH Roadmap for
Medical Research, grant U54 GM072970. See https://simtk.org.
Portions copyright (c) 2012-2013 Stanford University and the Authors.
Portions copyright (c) 2012-2015 Stanford University and the Authors.
Authors: Peter Eastman
Contributors:
......@@ -33,7 +33,9 @@ __version__ = "1.0"
import os
import xml.etree.ElementTree as etree
from simtk.unit import nanometers, sqrt
from simtk.openmm.vec3 import Vec3
from simtk.unit import nanometers, sqrt, is_quantity
from copy import deepcopy
class Topology(object):
"""Topology stores the topological information about a system.
......@@ -55,7 +57,7 @@ class Topology(object):
self._numResidues = 0
self._numAtoms = 0
self._bonds = []
self._unitCellDimensions = None
self._periodicBoxVectors = None
def addChain(self):
"""Create a new Chain and add it to the Topology.
......@@ -123,16 +125,48 @@ class Topology(object):
"""Iterate over all bonds (each represented as a tuple of two Atoms) in the Topology."""
return iter(self._bonds)
def getPeriodicBoxVectors(self):
"""Get the vectors defining the periodic box.
The return value may be None if this Topology does not represent a periodic structure."""
return self._periodicBoxVectors
def setPeriodicBoxVectors(self, vectors):
"""Set the vectors defining the periodic box."""
if vectors is not None:
if not is_quantity(vectors[0][0]):
vectors = vectors*nanometers
if vectors[0][1] != 0*nanometers or vectors[0][2] != 0*nanometers:
raise ValueError("First periodic box vector must be parallel to x.");
if vectors[1][2] != 0*nanometers:
raise ValueError("Second periodic box vector must be in the x-y plane.");
if vectors[0][0] <= 0*nanometers or vectors[1][1] <= 0*nanometers or vectors[2][2] <= 0*nanometers or vectors[0][0] < 2*abs(vectors[1][0]) or vectors[0][0] < 2*abs(vectors[2][0]) or vectors[1][1] < 2*abs(vectors[2][1]):
raise ValueError("Periodic box vectors must be in reduced form.");
self._periodicBoxVectors = deepcopy(vectors)
def getUnitCellDimensions(self):
"""Get the dimensions of the crystallographic unit cell.
The return value may be None if this Topology does not represent a periodic structure.
"""
return self._unitCellDimensions
if self._periodicBoxVectors is None:
return None
xsize = self._periodicBoxVectors[0][0].value_in_unit(nanometers)
ysize = self._periodicBoxVectors[1][1].value_in_unit(nanometers)
zsize = self._periodicBoxVectors[2][2].value_in_unit(nanometers)
return Vec3(xsize, ysize, zsize)*nanometers
def setUnitCellDimensions(self, dimensions):
"""Set the dimensions of the crystallographic unit cell."""
self._unitCellDimensions = dimensions
"""Set the dimensions of the crystallographic unit cell.
This method is an alternative to setPeriodicBoxVectors() for the case of a rectangular box. It sets
the box vectors to be orthogonal to each other and to have the specified lengths."""
if dimensions is None:
self._periodicBoxVectors = None
else:
if is_quantity(dimensions):
dimensions = dimensions.value_in_unit(nanometers)
self._periodicBoxVectors = (Vec3(dimensions[0], 0, 0), Vec3(0, dimensions[1], 0), Vec3(0, 0, dimensions[2]))*nanometers
@staticmethod
def loadBondDefinitions(file):
......
......@@ -20,7 +20,7 @@ class TestForceField(unittest.TestCase):
self.topology1 = self.pdb1.topology
self.topology1.setUnitCellDimensions(Vec3(2, 2, 2))
# alalnine dipeptide with implicit water
# alanine dipeptide with implicit water
self.pdb2 = PDBFile('systems/alanine-dipeptide-implicit.pdb')
self.forcefield2 = ForceField('amber99sb.xml', 'amber99_obc.xml')
......@@ -186,6 +186,18 @@ class TestForceField(unittest.TestCase):
system2 = ff2.createSystem(modeller.topology)
self.assertEqual(XmlSerializer.serialize(system1), XmlSerializer.serialize(system2))
def test_PeriodicBoxVectors(self):
"""Test setting the periodic box vectors."""
vectors = (Vec3(5, 0, 0), Vec3(-1.5, 4.5, 0), Vec3(0.4, 0.8, 7.5))*nanometers
self.pdb1.topology.setPeriodicBoxVectors(vectors)
self.assertEqual(Vec3(5, 4.5, 7.5)*nanometers, self.pdb1.topology.getUnitCellDimensions())
system = self.forcefield1.createSystem(self.pdb1.topology)
for i in range(3):
self.assertEqual(vectors[i], self.pdb1.topology.getPeriodicBoxVectors()[i])
self.assertEqual(vectors[i], system.getDefaultPeriodicBoxVectors()[i])
class AmoebaTestForceField(unittest.TestCase):
"""Test the ForceField.createSystem() method with the AMOEBA forcefield."""
......
import unittest
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
class TestGromacsGroFile(unittest.TestCase):
"""Test the Gromacs GRO file parser"""
def test_Triclinic(self):
"""Test parsing a file that describes a triclinic box."""
gro = GromacsGroFile('systems/triclinic.gro')
self.assertEqual(len(gro.positions), 8)
expectedPositions = [
Vec3(1.744, 2.788, 3.162),
Vec3(1.048, 0.762, 2.340),
Vec3(2.489, 1.570, 2.817),
Vec3(1.027, 1.893, 3.271),
Vec3(0.937, 0.825, 0.009),
Vec3(2.290, 1.887, 3.352),
Vec3(1.266, 1.111, 2.894),
Vec3(0.933, 1.862, 3.490)]*nanometers
for (p1, p2) in zip(expectedPositions, gro.positions):
self.assertEqual(p1, p2)
expectedVectors = [
Vec3(2.5, 0, 0),
Vec3(0.5, 3.0, 0),
Vec3(0.7, 0.9, 3.5)]*nanometers
for (v1, v2) in zip(expectedVectors, gro.getPeriodicBoxVectors()):
self.assertEqual(v1, v2)
self.assertEqual(Vec3(2.5, 3.0, 3.5)*nanometers, gro.getUnitCellDimensions())
for i in range(4):
self.assertEqual(elem.chlorine, gro.elements[i])
self.assertEqual('Cl', gro.atomNames[i])
self.assertEqual('Cl', gro.residueNames[i])
for i in range(4, 8):
self.assertEqual(elem.sodium, gro.elements[i])
self.assertEqual('Na', gro.atomNames[i])
self.assertEqual('Na', gro.residueNames[i])
if __name__ == '__main__':
unittest.main()
......@@ -308,44 +308,56 @@ class TestModeller(unittest.TestCase):
self.assertTrue(len(matoms)==0 and len(m1atoms)==1 and len(m2atoms)==1)
def test_addSolventPeriodicBox(self):
""" Test the addSolvent() method; test that the three ways of passing in the periodic box all work. """
""" Test the addSolvent() method; test that the four ways of passing in the periodic box all work. """
# First way of passing in periodic box dimensions: set it in the original topology.
# First way of passing in periodic box vectors: set it in the original topology.
topology_start = self.pdb.topology
topology_start.setUnitCellDimensions(Vec3(3.5, 4.5, 5.5)*nanometers)
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield)
topology_after = modeller.getTopology()
dim3 = topology_after.getUnitCellDimensions()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertAlmostEqual(dim3[0]/nanometers, 3.5)
self.assertAlmostEqual(dim3[1]/nanometers, 4.5)
self.assertAlmostEqual(dim3[2]/nanometers, 5.5)
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.5, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.5, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.5))
# Second way of passing in the periodic box dimensions: as a parameter to addSolvent()
# Second way of passing in the periodic box vectors: with the boxSize parameter to addSolvent()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, boxSize = Vec3(3.6, 4.6, 5.6)*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getUnitCellDimensions()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertAlmostEqual(dim3[0]/nanometers, 3.6)
self.assertAlmostEqual(dim3[1]/nanometers, 4.6)
self.assertAlmostEqual(dim3[2]/nanometers, 5.6)
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.6, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 4.6, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 5.6))
# Third way of passing in the periodic box dimensions: pass a 'padding' value to addSolvent()
# Third way of passing in the periodic box vectors: with the boxVectors parameter to addSolvent()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, boxVectors = (Vec3(3.4, 0, 0), Vec3(0.5, 4.4, 0), Vec3(-1.0, -1.5, 5.4))*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(3.4, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0.5, 4.4, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(-1.0, -1.5, 5.4))
# Fourth way of passing in the periodic box vectors: pass a 'padding' value to addSolvent()
topology_start = self.pdb.topology
modeller = Modeller(topology_start, self.positions)
modeller.deleteWater()
modeller.addSolvent(self.forcefield, padding = 1.0*nanometers)
topology_after = modeller.getTopology()
dim3 = topology_after.getUnitCellDimensions()
dim3 = topology_after.getPeriodicBoxVectors()
self.assertAlmostEqual(dim3[0]/nanometers, 2.8802)
self.assertAlmostEqual(dim3[1]/nanometers, 2.8802)
self.assertAlmostEqual(dim3[2]/nanometers, 2.8802)
self.assertVecAlmostEqual(dim3[0]/nanometers, Vec3(2.8802, 0, 0))
self.assertVecAlmostEqual(dim3[1]/nanometers, Vec3(0, 2.8802, 0))
self.assertVecAlmostEqual(dim3[2]/nanometers, Vec3(0, 0, 2.8802))
def test_addSolventNeutralSolvent(self):
""" Test the addSolvent() method; test adding ions to neutral solvent. """
......@@ -874,6 +886,12 @@ class TestModeller(unittest.TestCase):
validate_equivalence(self, topology_LYN, topology_after)
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
scale = max(1.0, norm(p1),)
for i in range(3):
diff = abs(p1[i]-p2[i])/scale
self.assertTrue(diff < tol)
if __name__ == '__main__':
unittest.main()
......
import unittest
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
import cStringIO
class TestPdbFile(unittest.TestCase):
"""Test the PDB file parser"""
def test_Triclinic(self):
"""Test parsing a file that describes a triclinic box."""
pdb = PDBFile('systems/triclinic.pdb')
self.assertEqual(len(pdb.positions), 8)
expectedPositions = [
Vec3(1.744, 2.788, 3.162),
Vec3(1.048, 0.762, 2.340),
Vec3(2.489, 1.570, 2.817),
Vec3(1.027, 1.893, 3.271),
Vec3(0.937, 0.825, 0.009),
Vec3(2.290, 1.887, 3.352),
Vec3(1.266, 1.111, 2.894),
Vec3(0.933, 1.862, 3.490)]*nanometers
for (p1, p2) in zip(expectedPositions, pdb.positions):
self.assertVecAlmostEqual(p1, p2)
expectedVectors = [
Vec3(2.5, 0, 0),
Vec3(0.5, 3.0, 0),
Vec3(0.7, 0.9, 3.5)]*nanometers
for (v1, v2) in zip(expectedVectors, pdb.topology.getPeriodicBoxVectors()):
self.assertVecAlmostEqual(v1, v2, 1e-4)
self.assertVecAlmostEqual(Vec3(2.5, 3.0, 3.5)*nanometers, pdb.topology.getUnitCellDimensions(), 1e-4)
for atom in pdb.topology.atoms():
if atom.index < 4:
self.assertEqual(elem.chlorine, atom.element)
self.assertEqual('Cl', atom.name)
self.assertEqual('Cl', atom.residue.name)
else:
self.assertEqual(elem.sodium, atom.element)
self.assertEqual('Na', atom.name)
self.assertEqual('Na', atom.residue.name)
def test_WriteFile(self):
"""Write a file, read it back, and make sure it matches the original."""
pdb1 = PDBFile('systems/triclinic.pdb')
output = cStringIO.StringIO()
PDBFile.writeFile(pdb1.topology, pdb1.positions, output)
input = cStringIO.StringIO(output.getvalue())
pdb2 = PDBFile(input)
output.close();
input.close();
self.assertEqual(len(pdb2.positions), 8)
for (p1, p2) in zip(pdb1.positions, pdb2.positions):
self.assertVecAlmostEqual(p1, p2)
for (v1, v2) in zip(pdb1.topology.getPeriodicBoxVectors(), pdb2.topology.getPeriodicBoxVectors()):
self.assertVecAlmostEqual(v1, v2, 1e-4)
self.assertVecAlmostEqual(pdb1.topology.getUnitCellDimensions(), pdb2.topology.getUnitCellDimensions(), 1e-4)
for atom1, atom2 in zip(pdb1.topology.atoms(), pdb2.topology.atoms()):
self.assertEqual(atom1.element, atom2.element)
self.assertEqual(atom1.name, atom2.name)
self.assertEqual(atom1.residue.name, atom2.residue.name)
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
unit = p1.unit
p1 = p1.value_in_unit(unit)
p2 = p2.value_in_unit(unit)
scale = max(1.0, norm(p1),)
for i in range(3):
diff = abs(p1[i]-p2[i])/scale
self.assertTrue(diff < tol)
if __name__ == '__main__':
unittest.main()
import unittest
from simtk.openmm.app import *
from simtk.openmm import *
from simtk.unit import *
import simtk.openmm.app.element as elem
class TestPdbxFile(unittest.TestCase):
"""Test the PDBx/mmCIF file parser"""
def test_Triclinic(self):
"""Test parsing a file that describes a triclinic box."""
pdb = PDBxFile('systems/triclinic.pdbx')
self.assertEqual(len(pdb.positions), 8)
expectedPositions = [
Vec3(1.744, 2.788, 3.162),
Vec3(1.048, 0.762, 2.340),
Vec3(2.489, 1.570, 2.817),
Vec3(1.027, 1.893, 3.271),
Vec3(0.937, 0.825, 0.009),
Vec3(2.290, 1.887, 3.352),
Vec3(1.266, 1.111, 2.894),
Vec3(0.933, 1.862, 3.490)]*nanometers
for (p1, p2) in zip(expectedPositions, pdb.positions):
self.assertVecAlmostEqual(p1, p2)
expectedVectors = [
Vec3(2.5, 0, 0),
Vec3(0.5, 3.0, 0),
Vec3(0.7, 0.9, 3.5)]*nanometers
for (v1, v2) in zip(expectedVectors, pdb.topology.getPeriodicBoxVectors()):
self.assertVecAlmostEqual(v1, v2, 1e-6)
self.assertVecAlmostEqual(Vec3(2.5, 3.0, 3.5)*nanometers, pdb.topology.getUnitCellDimensions(), 1e-6)
for atom in pdb.topology.atoms():
if atom.index < 4:
self.assertEqual(elem.chlorine, atom.element)
self.assertEqual('Cl', atom.name)
self.assertEqual('Cl', atom.residue.name)
else:
self.assertEqual(elem.sodium, atom.element)
self.assertEqual('Na', atom.name)
self.assertEqual('Na', atom.residue.name)
def assertVecAlmostEqual(self, p1, p2, tol=1e-7):
unit = p1.unit
p1 = p1.value_in_unit(unit)
p2 = p2.value_in_unit(unit)
scale = max(1.0, norm(p1),)
for i in range(3):
diff = abs(p1[i]-p2[i])/scale
self.assertTrue(diff < tol)
if __name__ == '__main__':
unittest.main()
Triclinic test system
8
1Cl Cl 1 1.744 2.788 3.162
2Cl Cl 2 1.048 0.762 2.340
3Cl Cl 3 2.489 1.570 2.817
4Cl Cl 4 1.027 1.893 3.271
5Na Na 5 0.937 0.825 0.009
6Na Na 6 2.290 1.887 3.352
7Na Na 7 1.266 1.111 2.894
8Na Na 8 0.933 1.862 3.490
2.50000 3.00000 3.50000 0.00000 0.00000 0.50000 0.00000 0.70000 0.90000
CRYST1 25.000 30.414 36.810 74.19 79.04 80.54 P 1 1
ATOM 1 Cl Cl A 1 17.440 27.880 31.620 1.00 0.00 Cl
ATOM 2 Cl Cl A 2 10.480 7.620 23.400 1.00 0.00 Cl
ATOM 3 Cl Cl A 3 24.890 15.700 28.170 1.00 0.00 Cl
ATOM 4 Cl Cl A 4 10.270 18.930 32.710 1.00 0.00 Cl
ATOM 5 Na Na A 5 9.370 8.250 0.090 1.00 0.00 Na
ATOM 6 Na Na A 6 22.900 18.870 33.520 1.00 0.00 Na
ATOM 7 Na Na A 7 12.660 11.110 28.940 1.00 0.00 Na
ATOM 8 Na Na A 8 9.330 18.620 34.900 1.00 0.00 Na
TER 9 Na A 8
END
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