PK vˆB˜[7;¯ ;¯ chemlab-v0.1/graphics.html
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 is in its early developement and it provides the most basic data structures. The molecular viewer has a solid ground and can actually draw and play trajectories in an efficient way. To get started be sure to check the User Manual.
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 work all in the same way, here is an example of GroHandler:
from chemlab.io import GroIO
infile = GroIO('waterbox.gro')
system = infile.read('system')
# Modify system as you wish...
outfile = GroIO('waterbox_out.gro')
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 __init__(self, filename):
self.filename = filename
def read(self, feature):
self.check_feature(feature, "read")
lines = open(self.filename).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))
open(self.filename, 'w').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 12.10 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 chemalb from the setup.py included in the package:
$ wget https://pypi.python.org/packages/source/c/chemlab/chemlab-0.1.tar.gz
$ tar xvzf chemlab-0.1.tar.gz
$ cd chemlab-0.1
$ 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 contained in 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])])
# 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().
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.
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.
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.boxsize = 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 ener.edr -e Pressure
$ chemlab gromacs ener.edr -e Temperature
$ chemlab gromacs ener.edr -e Potential
$ chemlab gromacs 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.