PK îu¡BÄç l( l( chemlab-v0.2/ipython.html
There is some preliminary integration between chemlab and ipython notebook, that will be extended and generalized in future releases. To see it in action, head over the example notebook
The chemlab.graphics package is one of the most interesting aspects of chemlab, that sets him apart from similar programs.
The purpose of the package is to provide a solid library to develop 3D applications to display chemical data in an flexible way. For example it’s extremely easy to build a molecular viewer and add a bunch of custom features to it.
The typical approach when developing a graphics application is to create a QtViewer instance and add 3D features to it:
>>> from chemlab.graphics import QtViewer
>>> v = QtViewer()
now let’s define a molecule. We can use the moldb module to get a water template.
>>> from chemlab.graphics.renderers import SphereRenderer
>>> from chemlab.data.moldb import water
>>> ar = v.add_renderer(AtomRenderer, water.r_array, water.type_array)
>>> v.run()
In this way you should be able to visualize a molecule where each atom is represented as a sphere. There are also a set of viewing controls:
In a similar fashion it is possible to display other features, such as boxes, arrows, lines, etc. It is useful to notice that with Viewer.add_renderer we are not passing an instance of the renderer, but we’re passing the renderer class and its respective constructor arguments. The method Viewer.add_renderer returns the actual instance.
It is possible as well to overlay 2D elements to a scene in a similar fashion, this will display a string at the screen position 300, 300:
from chemlab.graphics.uis import TextUI
tui = v.add_ui(TextUI, 300, 300, "Hello, World!")
Anyway, I encourage you to use the powerful Qt framework to provide interaction and widgets to your application.
Renderers are simply classes used to draw 3D objects. They are tecnically required to provide just one method, draw and they must take an instance of QChemlabWidget as their first argument (check out the AbstractRenderer class). In this way they provide the maximum flexibility required to build efficient opengl routines. Renderers may be subclass other renderers as well as use other renderers.
A very useful renderer is TriangleRenderer, used to render efficiently a list of triangles, it constitutes a base for writing other renderers. TriangleRenderer works basically like this, you pass the vertices, normals and colors of the triangle and it will display a triangle in the world:
from chemlab.graphics import QtViewer
from chemlab.graphics.renderers import TriangleRenderer
from chemlab.graphics.colors import green
import numpy as np
vertices = np.array([[0.0, 0.0, 0.0], [0.0, 1.0, 0.0], [1.0, 0.0, 0.0]])
normals = np.array([[0.0, 0.0, 1.0], [0.0, 0.0, 1.0], [0.0, 0.0, 1.0]])
colors = np.array([green, green, green])
v = QtViewer()
v.add_renderer(TriangleRenderer, vertices, normals, colors)
v.run()
If you pass 6 vertices/normals/colors, he will display 2 triangles and so on. As a sidenote, he is very efficient and in fact chemlab.graphics.renderers.TriangleRenderer is used as a backend for a lot of other renderers such as SphereRenderer and CylinderRenderer. If you can reduce a shape in triangles, you can easily write a renderer for it.
In addition to that, TriangleRenderer provides also a method to update vertices, normals and colors. We can demonstrate that from the last example by defining an update function that rotates our triangle:
from chemlab.graphics.transformations import rotation_matrix
def update():
y_axis = np.array([0.0, 1.0, 0.0])
# We take the [:3,:3] part because rotation_matrix can be used to
# rotate homogeneous (4D) coordinates.
rot = rotation_matrix(3.14/32, y_axis)[:3, :3]
# This is the numpy-efficient way of applying rot to each coordinate
vertices[:] = np.dot(vertices, rot.T)
normals[:] = np.dot(vertices, rot.T)
tr.update_vertices(vertices)
tr.update_normals(normals)
v.widget.repaint()
v.schedule(update, 10)
v.run()
On this ground we can develop a TetrahedronRenderer based on our TriangleRenderer. To do that we first need to understand how a tetrahedron is made, and how can we define the vertices that make the tetrahedron.
First of all, we need to have the 4 coordinates that represents a tetrahedron. Without even trying to visualize it, just pick the values straight from Wikipedia:
import numpy as np
v1 = np.array([1.0, 0.0, -1.0/np.sqrt(2)])
v2 = np.array([-1.0, 0.0, -1.0/np.sqrt(2)])
v3 = np.array([0.0, 1.0, 1.0/np.sqrt(2)])
v4 = np.array([0.0, -1.0, 1.0/np.sqrt(2)])
We can quickly verify if this is correcty by using a PointRenderer:
from chemlab.graphics import QtViewer
from chemlab.graphics.renderers import PointRenderer
from chemlab.graphics.colors import black, green, blue, red
colors = [black, green, blue, red]
v = QtViewer()
v.add_renderer(PointRenderer, np.array([v1, v2, v3, v4]), colors)
v.run()
We’ve got 4 boring points that look like they’re at the vertices of a tetrahedron. Most importantly we learned that we can use PointRenderer to quickly test shapes.
Now let’s define the four triangles (12 vertices) that represent a solid tetrahedron. It is good practice to put the triangle vertices in a certain order to estabilish which face is pointing outside and which one is pointing inside for optimization reasons. The convention is that if we specify 3 triangle vertices in clockwise order this means that the face points outwards from the solid:
We can therefore write our vertices and colors:
vertices = np.array([
v1, v4, v3,
v3, v4, v2,
v1, v3, v2,
v2, v4, v1
])
colors = [green] * 12
All is left to do is write the normals to the surface at each vertex. This is easily done by calculating the cross product of the vectors constituting two sides of a triangle, (remember that the normals should point outward):
n1 = -np.cross(v4 - v1, v3 - v1)
n2 = -np.cross(v4 - v3, v2 - v3)
n3 = -np.cross(v3 - v1, v2 - v1)
n4 = -np.cross(v4 - v2, v1 - v2)
normals = [n1, n1, n1,
n2, n2, n2,
n3, n3, n3,
n4, n4, n4]
from chemlab.graphics.renderers import TriangleRenderer
v.add_renderer(TriangleRenderer, vertices, normals, colors)
v.run()
Now that we’ve got the basic shape in place we can code the actual Renderer class to be used directly with the viewer. We will make a renderer that, given a set of coordinates will display many tetrahedra.
We can start by defining a Renderer class, inheriting from AbstractRenderer, the main thing you should notice is that you need an additional argument widget that will be passed when you use the method QtViewer.add_renderer:
from chemlab.graphics.renderers import AbstractRenderer
class TetrahedraRenderer(AbstractRenderer):
def __init__(self, widget, positions):
super(TetrahedraRenderer, self).__init__(widget)
...
The strategy to implement a multiple-tetrahedron renderer will be like this:
You can see the code in this snippet:
class TetrahedraRenderer(AbstractRenderer):
def __init__(self, widget, positions):
super(TetrahedraRenderer, self).__init__(widget)
v1 = np.array([1.0, 0.0, -1.0/np.sqrt(2)])
v2 = np.array([-1.0, 0.0, -1.0/np.sqrt(2)])
v3 = np.array([0.0, 1.0, 1.0/np.sqrt(2)])
v4 = np.array([0.0, -1.0, 1.0/np.sqrt(2)])
positions = np.array(positions)
# Vertices of a single tetrahedra
self._th_vertices = np.array([
v1, v4, v3,
v3, v4, v2,
v1, v3, v2,
v2, v4, v1
])
self._th_normals = np.array([
n1, n1, n1,
n2, n2, n2,
n3, n3, n3,
n4, n4, n4])
self.n_tetra = len(positions)
tot_vertices = []
for pos in positions:
tot_vertices.extend(self._th_vertices + pos)
# Refer to numpy.tile, this simply repeats the elements
# of the array in an efficient manner.
tot_normals = np.tile(self._th_normals, (self.n_tetra, 1))
tot_colors = [green] * self.n_tetra * 12
# !NOTICE! that we have to pass widget as the first argument
self.tr = TriangleRenderer(widget, tot_vertices,
tot_normals, tot_colors)
def draw(self):
self.tr.draw()
To demostrate let’s draw a grid of 125 tetrahedra:
positions = []
for x in range(5):
for y in range(5):
for z in range(5):
positions.append([float(x)*2, float(y)*2, float(z)*2])
v.add_renderer(TetrahedraRenderer, positions)
v.widget.camera.position = np.array([0.0, 0.0, 20.0])
v.run()
If you had any problem with the tutorial or you want to implement other kind of renderers don’t exitate to contact me. The full code of this tutorial is in chemlab/examples/tetrahedra_tutorial.py.
Webpage: | https://chemlab.github.com/chemlab |
---|---|
Project Page: | https://github.com/chemlab/chemlab |
Mailing List: | python-chemlab.googlegroups.com |
Downloads: | https://chemlab.github.com/chemlab |
Chemlab is a library that can help the user with chemistry-relevant calculations using the flexibility and power of the python programming language. It aims to be well-designed and pythonic, taking inspiration from project such as numpy and scipy.
Chemlab long term goal is to be:
Computational and theoretical chemistry is a huge field, and providing a program that encompasses all aspect of it is an impossible task. The spirit of chemlab is to provide a common ground from where you can build specific programs. For this reason it includes an fully programmable molecular viewer.
Chemlab includes a lot of utilities to programmatically download and generate geometries. The molecular viewer is very fast (it can easily animate ~100000 spheres) and the design is simple and flexible. For more information about the newest features check out the release notes in the What’s new document.
Chemlab is developer-friendly, it provides good documentation and has an easy structure to get in. Feel free to send me anything that you may do with chemlab, like supporting a new file format, a new graphic renderer, a nice example, even if you don’think it’s perfect. Send an email to the mailing list or file an issue on the github page to discuss any idea that comes to your mind. Get involved!
Table of Contents
Packages
There are a lot of file formats used and produced by chemistry applications. Each program has his way to store geometries, trajectories, energies and properties etc. chemlab tries to encompass all of those different properties by using a lightweight way to handle such differences.
The classes responsible for the I/O are subclasses of chemlab.io.handlers.IOHandler. These handlers take a file-like object as the first argument and they work all in the same way, here is an example of GroHandler:
from chemlab.io import GroIO
fd = open('waterbox.gro')
infile = GroIO(fd)
system = infile.read('system')
# Modify system as you wish...
fd = open('waterbox_out.gro', 'w')
outfile = GroIO(fd)
outfile.write('system', system)
You first create the handler instance for a certain format and then you can read a certain feature provided by the handler. In this example we read and write the system feature.
Some file formats may have some extra data for each atom, molecule or system. For example the ”.gro” file formats have his own way to call the atoms in a water molecule: OW, HW1, HW2. To handle such issues, you can write this information in the export arrays contained in the data structures, such as Atom.export, Molecule.export, and their array-based counterparts Molecule.atom_export_array, System.mol_export and System.atom_export_array.
Those attributes are especially important where you write in some data format, since you may have to provide those attribute when you initialize your Atom, Molecule and System.
You can easily open a data file without even having to search his format handler by using the utility function chemlab.io.datafile():
from chemlab.io import datafile
sys = datafile('waterbox.gro').read('system')
t, coords = datafile('traj.xtc').read('trajectory')
See also
Implementing or improving an existing IOHandler is a great way to partecipate in chemlab development. Fortuately, it’s extremely easy to setup one of them.
It boils down to a few steps:
Here is an example of the xyz handler:
import numpy as np
from chemlab.io.handlers import IOHandler
from chemlab.core import Molecule
class XyzIO(IOHandler):
'''The XYZ format is described in this wikipedia article
http://en.wikipedia.org/wiki/XYZ_file_format.
**Features**
.. method:: read("molecule")
Read the coordinates in a :py:class:`~chemlab.core.Molecule` instance.
.. method:: write("molecule", mol)
Writes a :py:class:`~chemlab.core.Molecule` instance in the XYZ format.
'''
can_read = ['molecule']
can_write = ['molecule']
def read(self, feature):
self.check_feature(feature, "read")
lines = self.fd.readlines()
num = int(lines[0])
title = lines[1]
if feature == 'title':
return title
if feature == 'molecule':
type_array = []
r_array = []
for l in lines[2:]:
type, x, y, z = l.split()
r_array.append([float(x),float(y),float(z)])
type_array.append(type)
r_array = np.array(r_array)/10 # To nm
type_array = np.array(type_array)
return Molecule.from_arrays(r_array=r_array, type_array=type_array)
def write(self, feature, mol):
self.check_feature(feature, "write")
lines = []
if feature == 'molecule':
lines.append(str(mol.n_atoms))
lines.append('Generated by chemlab')
for t, (x, y, z) in zip(mol.type_array, mol.r_array):
lines.append(' %s %.6f %.6f %.6f' %
(t, x*10, y*10, z*10))
self.fd.write('\n'.join(lines))
A few remarks:
check_feature() before performing read/write. This will check that the feature is present in the can_read/can_write list;
Molecule.from_arrays() and System.from_arrays();
EdrIO handler does not read Molecule or System at all;
You can definitely take inspiration from the handlers included in chemlab, Supported File Formats.
chemlab is currently tested on Ubuntu 13.04 and python 2.7. First install the dependencies:
$ sudo apt-get install python-numpy python-scipy python-matplotlib python-pyside python-opengl cython
Download unpack and install chemlab from the setup.py included in the package:
$ wget https://pypi.python.org/packages/source/c/chemlab/chemlab-0.2.tar.gz
$ tar xvzf chemlab-0.2.tar.gz
$ cd chemlab-0.2
$ sudo python setup.py install
Test the newly installed package by typing:
$ chemlab view tests/data/cry.gro
The molecular viewer should display a crystal, if not, file an issue on github.
Once you’re setup, you’re ready to to dig in chemlab’s features. You may start from the User Manual.
After installing the dependencies, grab the chemlab source from git:
$ git clone --recursive https://github.com/chemlab/chemlab.git
Complile the included extensions:
$ python setup.py build_ext --inplace
Just add the chemlab directory to the PYTHONPATH in your .bashrc:
export PYTHONPATH=$PYTHONPATH:/path/to/chemlab
In chemlab, atoms can be represented using the chemlab.core.Atom data structure that contains some common information about our particles like type, mass and position. Atom instances are easily created by initializing them with data
>>> from chemlab.core import Atom
>>> ar = Atom('Ar', [0.0, 0.0, 0.0])
>>> ar.type
'Ar'
>>> ar.r
np.array([0.0, 0.0, 0.0])
Note
for the atomic coordinates you should use nanometers
A chemlab.core.Molecule is an entity composed of more atoms and most of the Molecule properties are inherited from the constituent atoms. To initialize a Molecule you can, for example pass a list of atom instances to its constructor:
>>> from chemlab.core import Molecule
>>> mol = Molecule([at1, at2, at3])
Molecules are easily and efficiently manipulated through the use of numpy arrays. One of the most useful arrays contained in Molecule is the array of coordinates Molecule.r_array. The array of coordinates is a numpy array of shape (NA,3) where NA is the number of atoms in the molecule. According to the numpy broadcasting rules, if you sum two arrays with shapes (NA,3) and (3,), each row of the first array get summed by the second array. Let’s say we have a water molecule and we want to displace it randomly in a box, this is easily accomplished by initializing a Molecule at the origin and summing its coordinates by a random displacement:
import numpy as np
wat = Molecule([Atom("H", [0.0, 0.0, 0.0]),
Atom("H", [0.0, 1.0, 0.0]),
Atom("O", [0.0, 0.0, 1.0])], bonds=[[2, 0], [2, 1]])
# Shapes (NA, 3) and (3,)
wat.r_array += np.random.rand(3)
Using the same principles you can also apply other kinds of transformations such as matrices. You can for example rotate the molecule by 90 degrees around the z-axis:
from chemlab.graphics.transformations import rotation_matrix
# The transformation module returns 4x4 matrices
M = rotation_matrix(np.pi/2, np.array([0.0, 0.0, 1.0]))[:3,:3]
# slow, readable way
for i,r in enumerate(wat.r_array):
wat.r_array[i] = np.dot(M,r)
# numpy efficient way to do the same:
# wat.r_array = np.dot(wat.r_array, M.T)
The array-based API provides a massive increase in performance and a more straightforward integration with C libraries thanks to the numpy arrays. This feature comes at a cost: the data is copied between atoms and molecules, in other words the changes in the costituents atoms are not reflected in the Molecule and viceversa. Even if it may look a bit unnatural, this approach limits side effects making the code more predictable and easy to follow.
In context such as molecular simulations it is customary to introduce a new data structure called System. A System represents a collection of molecules, and optionally (but recommended) you can pass also periodic box information:
>>> from chemlab.core import System
# molecule = a list of Molecule instances
>>> s = System(molecules, boxsize=2.0)
System do not take directly Atom instances as its constituents, therefore if you need to simulate a system made of single atoms (say, a box of liquid Ar) you need to wrap the atoms into a Molecule:
>>> ar = Atom('Ar', [0.0, 0.0, 0.0])
>>> mol = Molecule([ar])
System, similarly to Molecule can expose data by using arrays and it inherits atomic data from the constituent molecules. For instance, you can easily and efficiently access all the atomic coordinates by using the attribute System.r_array. To understand the relation between Atom.r, Molecule.r_array and System.r_array you can refer to the picture below:
You can preallocate a System by using the classmethod System.empty (pretty much like you can preallocate numpy arrays with np.empty or np.zeros) and then add the molecules one by one:
import numpy as np
from chemlab.core import Atom, Molecule, System
from chemlab.graphics import display_system
# Template molecule
wat = Molecule([Atom('O', [0.00, 0.00, 0.01]),
Atom('H', [0.00, 0.08,-0.05]),
Atom('H', [0.00,-0.08,-0.05])])
# Initialize a system with four water molecules.
s = System.empty(4, 12) # 4 molecules, 12 atoms
for i in range(4):
wat.move_to(np.random.rand(3)) # randomly displace the water molecule
s.add(wat) # data gets copied each time
display_system(s)
Since the data is copied, the wat molecule act as a template so you can move it around and keep adding it to the System.
Preallocating and adding molecules is a pretty fast way to build a System, but the fastest way (in terms of processing time) is to build the system by passing ready-made arrays, this is done by using chemlab.core.System.from_arrays().
It is possible to build boxes where atoms are placed randomly by using the chemlab.core.random_lattice_box() function. A set of template molecules are copied and translated randomly on the points of a 3d lattice. This ensure that the spacing between molecules is consistent and to avoid overlaps.
To make an example box:
from chemlab.db import ChemlabDB
from chemlab.core import random_lattice_box
# Example water molecule
water = ChemlabDB().get('molecule', 'example.water')
s = random_lattice_box([water], [1000], [4.0, 4.0, 4.0])
chemlab provides an handy way to build crystal structures from the atomic coordinates and the space group information. If you have the crystallographic data, you can easily build a crystal:
from chemlab.core import Atom, Molecule, crystal
from chemlab.graphics import display_system
# Molecule templates
na = Molecule([Atom('Na', [0.0, 0.0, 0.0])])
cl = Molecule([Atom('Cl', [0.0, 0.0, 0.0])])
s = crystal([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]], # Fractional Positions
[na, cl], # Molecules
225, # Space Group
cellpar = [.54, .54, .54, 90, 90, 90], # unit cell parameters
repetitions = [5, 5, 5]) # unit cell repetitions in each direction
display_system(s)
See also
Note
If you’d like to implement a .cif file reader, you’re welcome! Drop a patch on github.
You can manipulate systems by using some simple but flexible functions. It is really easy to generate a system by selecting a part from a bigger system, this is implemented in the functions chemlab.core.subsystem_from_atoms() and chemlab.core.subsystem_from_molecules().
Those two functions take as first argument the original System, and as the second argument a selection. A selection is either a boolean array that is True when we want to select that element and False otherwise or an integer array containing the elements that we want to select. By using those two functions we can create subsystem by building those selections.
The following example shows an easy way to take the molecules that contain atoms in the region of space x > 0.5 by employing subsystem_from_atoms():
import numpy as np
from chemlab.core import crystal, Molecule, Atom, subsystem_from_atoms
from chemlab.graphics import display_system
# Template molecule
wat = Molecule([Atom('O', [0.00, 0.00, 0.01]),
Atom('H', [0.00, 0.08,-0.05]),
Atom('H', [0.00,-0.08,-0.05])])
s = crystal([[0.0, 0.0, 0.0]], [wat], 225,
cellpar = [.54, .54, .54, 90, 90, 90], # unit cell parameters
repetitions = [5, 5, 5]) # unit cell repetitions in each direction
selection = s.r_array[:, 0] > 0.5
sub_s = subsystem_from_atoms(s, selection)
display_system(sub_s)
It is also possible to select a subsystem by selecting specific molecules, in the following example we select the first 10 water molecules by using subsystem_from_molecules():
from chemlab.core import subsystem_from_molecules
selection = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])
sub_s = subsystem_from_molecules(s, selection)
Note
chemlab will provide other selection utilities in the future, if you have a specific request, file an issue on github
You can also create a system by merging two different systems. In the following example we will see how to make a NaCl/H2O interface by using chemlab.core.merge_systems():
import numpy as np
from chemlab.core import Atom, Molecule, crystal
from chemlab.core import subsystem_from_atoms, merge_systems
from chemlab.graphics import display_system
# Make water crystal
wat = Molecule([Atom('O', [0.00, 0.00, 0.01]),
Atom('H', [0.00, 0.08,-0.05]),
Atom('H', [0.00,-0.08,-0.05])])
water_crystal = crystal([[0.0, 0.0, 0.0]], [wat], 225,
cellpar = [.54, .54, .54, 90, 90, 90], # unit cell parameters
repetitions = [5, 5, 5]) # unit cell repetitions in each direction
# Make nacl crystal
na = Molecule([Atom('Na', [0.0, 0.0, 0.0])])
cl = Molecule([Atom('Cl', [0.0, 0.0, 0.0])])
nacl_crystal = crystal([[0.0, 0.0, 0.0], [0.5, 0.5, 0.5]], [na, cl], 225,
cellpar = [.54, .54, .54, 90, 90, 90],
repetitions = [5, 5, 5])
water_half = subsystem_from_atoms(water_crystal,
water_crystal.r_array[:,0] > 1.2)
nacl_half = subsystem_from_atoms(nacl_crystal,
nacl_crystal.r_array[:,0] < 1.2)
interface = merge_systems(water_half, nacl_half)
display_system(interface)
At the present time, the merging will avoid overlapping by creating a bounding box around the two systems and removing the molecules of the first system that are inside the second system bounding box. In the future there will be more clever ways to handle this overlaps.
There are two methods used to remove specific atoms and molecules from a system. chemlab.core.System.remove_molecules() and chemlab.core.System.remove_atoms(). Taking from the previous NaCl example, you may need to remove some excess ions to met the electroneutrality condition:
# n_na and n_cl are the number of Na and Cl molecules
toremove = 'Na' if n_na > n_cl else 'Cl'
nremove = abs(n_na - n_cl) # Number of indices to be removed
remove_indices = (s.type_array == toremove).nonzero()[0][:nremove]
s.remove_atoms(rem_indices)
It is possible to reorder the molecules in a System by using the method chemlab.core.System.reorder_molecules() that takes the new order as the first argument. Reordering can be useful for example to sort the molecules against a certain key.
If you use chemlab in conjunction with GROMACS, you may use the chemlab.core.System.sort() to sort the molecules according to their molecular formulas before exporting. The topology file expect to have a file with the same molecule type ordererd.
Warning
This part of chemlab is still in draft. This first, very brief implementation serves as a specification document. As we collect more feedback and feature requests there will be an expansion and a refinement of the extension functionalities.
Differents applications of chemistry may require additional data attached to each atom, molecule or system. For example you may need the velocity of the system, atomic charges or number of electrons. Chemlab should be able to provide a way to simply attach this data while retaining the selection and sorting functionalities.
The management of the atomic and molecular properties within a System is done through specific handlers. Those handlers are called attributes and fields. In the following example we may see how it’s possible to add a new field “v” to the Atom class, and transmitting this field as a “v_array” in the Molecule and System class. In those cases they basically take as their argument the attribute/field name, the type, and a function that return the default value for the field/attribute:
from chemlab.core.attributes import MArrayAttr, NDArrayAttr
from chemlab.core.fields import AtomicField
class MyAtom(Atom):
fields = Atom.fields + [AtomicField("v",
default=lambda at: np.zeros(3, np.float))]
class MyMolecule(Molecule):
attributes = Molecule.attributes + [MArrayAttr("v_array", "v", np.float,
default=lambda mol: np.zeros((mol.n_atoms, 3), np.float))]
class MySystem(System):
attributes = System.attributes + [NDArrayAttr("v_array", "v_array", np.float, 3)]
Those class are ready to use. You may want to create new instances with the Atom.from_fields, Molecule.from_arrays and System.from_arrays.
Once you’ve done your field-specific job with MyAtom/MyMolecule/MySystem you can convert back to a chemlab default class class by using the astype methods:
at = myat.astype(Atom)
mol = mymol.astype(Molecule)
sys = mysys.astype(System)
GROMACS is one of the most used packages for molecular simulations, chemlab can provide a modern and intuitive interface to generate input and analyze the output of GROMACS calculations. To illustrate the concepts we’ll perform a very simple simulation of liquid water.
This depends on the system you’re using but I believe that GROMACS is already packaged for most linux distributions and also for other operating systems.
In Ubuntu:
$ sudo apt-get install gromacs
In order to run a minimum simulation GROMACS requires to know some basic properties of the system we intend to simulate. This boils down to basically 3 ingredients:
chemlab can help us to build any system that we want and we’ll use it to write a ”.gro” file. Then we will use chemlab to visualize and analyze the result of the GROMACS simulation.
There are many ways to generate a box of water, in our example we will place 512 water molecules in a cubic grid. The advantages of doing that is the simplicity of the approach and the fact that we are naturally avoid any overlap between adiacent molecules.
To generate such a box we will:
import numpy as np
from chemlab.core import Atom, Molecule, System
from chemlab.graphics import display_system
# Spacing between two grid points
spacing = 0.3
# an 8x8x8 grid, for a total of 512 points
grid_size = (8, 8, 8)
# Preallocate the system
# 512 molecules, and 512*3 atoms
s = System.empty(512, 512*3)
# Water template, it contains export informations for gromacs
# more about export later...
water_tmp = Molecule([Atom('O', [0.0, 0.0, 0.0], export={'grotype': 'OW'}),
Atom('H', [0.1, 0.0, 0.0], export={'grotype': 'HW1'}),
Atom('H', [-0.03333, 0.09428, 0.0], export={'grotype':'HW2'})],
export={'groname': 'SOL'})
for a in range(grid_size[0]):
for b in range(grid_size[1]):
for c in range(grid_size[2]):
grid_point = np.array([a,b,c]) * spacing # array operation
water_tmp.move_to(grid_point)
s.add(water_tmp)
# Adjust boxsize for periodic boundary conditions
s.box_vectors = np.eye(3) * (8 * spacing)
# Visualize to verify that the system was setup correctly
display_system(s)
If you run this, it will display the following window:
Awesome! Now we can write the ”.gro” file. Notice that when we defined our water molecule we had to pass an export dictionary to the atoms and molecules. The export mechanism is the way used by chemlab to handle all the variety of different file formats.
In this specific case, gromacs defines its own atom and molecule names in the ”.top” file and then matches those to the ”.gro” file to infer the bonds and interactions.
TODO Add picture of the export dictionary
How do we write the .gro file? Since we’ve already setup our export information, this is an one-liner:
from chemlab.io import datafile
datafile("start.gro").write("system", s)
I’ll give you directly the gromacs input files to do an NPT simulation of water, just create those files in your working directory:
topol.top
; We simply import ready-made definitions for the molecule type
; SOL and the atom types OW, HW1 and HW2
#include "ffoplsaa.itp"
#include "spce.itp"
[ system ]
Simple box of water
[ molecules ]
SOL 512
run.mdp
integrator = md
dt = 0.001
nsteps = 200000
nstxtcout = 100
rlist = 0.9
coulombtype = pme
rcoulomb = 0.9
rvdw = 0.9
dispcorr = enerpres
tcoupl = v-rescale
tc-grps = System
ref_t = 300
tau_t = 0.1
pcoupl = berendsen
compressibility = 4.5e-5
ref_p = 1.0
gen_vel = yes
gen_temp = 300
constraints = all-bonds
To run the simulation with gromacs we have to do two steps:
Generate a parameter input, this will check that our input make sense before running the simulation:
grompp_d -f run.mdp -c start.gro -p topol.top
This will generate a bunch of files in your working directory.
Now we run the simulation, in the meantime, go grab coffee:
mdrun_d -v
This will take a while depending on your machine. If you are not a coffee drinker, don’t worry, you can stop the simulation by pressing Ctrl-C. The good news is that chemlab can read files from partial runs!
To quickly preview trajectories and system energies you can use the script chemlab included in the distribution in scripts/chemlab.
GROMACS can store the trajectory (in the form of atomic coordinates) in the .xtc file. To play the trajectory you can use the command:
$ chemlab view start.gro --traj traj.xtc
Note
the nstxtcout = 100 option in the mdp file sets the output frequency in the xtc file
You may also be interested to look at some other properties, such as the potential energy, pressure, temperature and density. This information is written by GROMACS in the ”.edr” file. You can use the chemlab script to view that:
$ chemlab gromacs energy ener.edr -e Pressure
$ chemlab gromacs energy ener.edr -e Temperature
$ chemlab gromacs energy ener.edr -e Potential
$ chemlab gromacs energy ener.edr -e Density
Warning
The chemlab gromacs command is a work in progress, the syntax may change in the future.
It is also possible to view and get the results by directly reading the files and have direct access to the xtc coordinates and the energy stored in the edr files. Take a look at the reference for chemlab.io.handlers.XtcIO and chemlab.io.handlers.EdrIO.
The tutorial is over, if you have any problem or want to know more, just drop an email on the mailing list python-chemlab@googlegroups.com or file an issue on github https://github.com/chemlab/chemlab/issues
Please activate JavaScript to enable the search functionality.
From here you can search these documents. Enter your search words into the box below and click "search". Note that the search function will automatically search for all of the words. Pages containing fewer words won't appear in the result list.
chemlab.core:
- Serialization through json with from_json and tojson for Atom, System and Molecule;
- Removing atoms and molecules from System. System.remove_atoms, System.remove_molecules;
- Experimental support for customized Atom/Molecule/System types.
- Some indexing routines: System.atom_to_molecules_indices and System.mol_to_atom_indices;
- Custom sorting of systems throught System.reorder_molecules;
- Support for bonds in molecules and experimental support for bonds in Systems throught Molecule.bonds and System.get_bonds_array
- System.merge_systems has a better overlap handling.
- Removed boxsize attribute, now you have to always specify box_vectors.
- Implemented random_lattice_box to do random solvent boxes.
chemlab.graphics:
- New Renderers: - BallAndStickRenderer - BondRenderer - WireframeRenderer
- Implemented Camera.autozoom for autoscaling
- Reimplemented BondRenderer in cython.
chemlab.io:
- New Handlers:
- MDL molfile (.mol extension)
- Chemical Markup Language (.cml extension)
chemlab.db:
- New package to handle databases
- CirDB to retrieve molecules from chemical identifier resolver
- ChemlabDB to retrieve built-in data
- LocalDB to make personal databases
chemlab.ipython:
- Preliminary ipython notebook integration. It can display Molecule and System instances by using out-of-screen rendering.
chemlab.utils:
- Implemented some (periodic/non-periodic) distance finding routines.