Programming the Ideal Gas

You've heard about it before. Many times. It's a model system used a lot in chemistry. The Ideal Gas law: A system where the particles can be treated as non-interacting. How do you program such a system? What do you need to describe it? Given 1) the basic introduction covered by Prof. Jan Jensen in another lecture, and 2) the following assignment (notice it is part 2), here is what the students are supposed to do once this exercise is complete:

1. Take a box with dimensions $(-1<x<1$ and $-1<y<1)$.
2. Make npart particles placed randomly in the right side of a box ($0<x<1$) but evenly from the top through the bottom
3. Give the npart particles random velocities in both directions and make it so they can move in both the positive and negative direction.
4. Move the particles around via the equation: $X_i^{n+1} = X_i^{n}+dt*dX_i^{n}$
5. Keep them inside the box at all times.
6. Plot the result.

The students start out with 1, 2 and 6 presented to them, both for inspiration and because they contain some subtle hints. The code is

One of the tricky parts is knowing "what is the i'th particle" because we need this information to update the system according to the equation above in point 4. Once you realize that you can access this information in the X, Y, dX and dY arrays as X[i] under a loop, then there is really no problem with that point either, for instance, to update all particles according to the equation in 4 in one go, we could write

Assigning a loop to iterate over the particles is not really troublesome, but understanding why is hard - even to explain it to yourself. A simple for-loop might not be the smartest thing to do, but I think it is the easiest one to go about and more importantly I often feel that many things can get lost in the whole numpy magic, especially for newcomers.

The next part is, that we should keep the particles inside the box at all times. The exercise has a little discussion about which version is the better, but in my opinion the particles should always be inside at all times. To keep them in, we ask whether a step (equation in point 4) in X and Y will result in the particle being outside the box and if that is the case, we reverse velocity that potentially brought the particle outside. We could use the following code to do that

Inserting the above two snippets of code, under a loop over the number of steps we want to take, we finally end up with the code that simulates non-interacting particles in a container.

Try and code it for yourself or use it in a course you teach. The students do find it interesting, both in the form of "what do I need to tell the computer to make it do stuff?" and "oh - the particles do move!".

Disclaimer: This code is part of a course I help teach at the University of Copenhagen - just imagine you never saw python before.

Obtaining Contact Numbers using SMARTS

Because I'm on a roll with blogging right now and I am looking for excuses to not write the next paper, here is a post about how to utilize SMARTS[1] for something cheminformatically-vant such as the contact number[2]. We will base the new code on what I did in the last post I wrote. The contact number is basically "Select an $\alpha$-Carbon and count the number $\alpha$-Carbons within a sphere of a certain radius $R$ from it". This give rise to a 1D measurement (see below) of how buried residues are. If you massage your data enough (we will get to this in a later post) it can be used to validate whether a protein is folded (somewhat) correctly or not. There is much more to this story than I have room for on my blog, but I do recommend you follow the blog of Anders Christensen, a fellow student who folds more proteins than I.

Our work today shall start from Crambine (PDB: 1CRN) which I've prepared for your viewing pleasure in pymol

From last week, the boiler-plate code we shall start from is

import openbabel
filename = "1CRN.pdb"
pattern = "[$(CC(=O)[N,O])]"
obmol = openbabel.OBMol()
obpat = openbabel.OBSmartsPattern()
obconv = openbabel.OBConversion()
obconv.SetInFormat(filename[-3:])
obconv.ReadFile(obmol, filename)
obpat.Init(pattern)
obpat.Match(obmol)
alpha_carbons = [m[0] for m in obpat.GetUMapList()]

I'm interested in $\alpha$-Carbons so the atomic primitive pattern we use to match with will be

[$(CC(=O)[N,O])]

The last part, [N,O] is a match of either a Nitrogen or an Oxygen. This Oxygen match is needed for the C-terminal. Because we are using SMARTS, this is the only thing that will change in our code once it is done. Should we be interested in locating the carboxyl groups of the backbone for instance, simply change one pattern. Of course, if you want to match $\alpha$-Carbons to carboxyl groups, this requires an extra pattern match.

Running this code will return a list of atoms. To calculate the contact number, the following function iterates over each match, forms a pair with all other matches and calculates the contact number for one residue at a time

def ContactNumbersForCutoff(matches, rcut):
  contact_numbers = []
  for i in matches:
    contact_number = 0 
    ai = obmol.GetAtom(i)
    for j in matches:
      if i == j: continue
      aj = obmol.GetAtom(j)
      if ai.GetDistance(aj) < rcut:
        contact_number +=1

    contact_numbers.append(contact_number)
  return contact_numbers

which we will invoke to generate the list of contact numbers

cn12 = ContactNumbersForCutoff(alpha_carbons, 12.0)

Since I am fond of the matplotlib[3] library for python when dealing with plots and I don't feel like doing the hard work myself anyways, I'll just plot[4] the list of numbers using the following code

pylab.plot(cn12)
pylab.savefig("contact.png")

which (after some additional touches - see the github page for the final code) results in

Here we see that residues 20, 38, 39 and 40 are clearly not buried when compared to residues 30 to 36. In a later posts, I will add a few lines of code to calculate the above information in a slightly more useful way. I'd also like to explore a quick implementation of the Half-Sphere Exposure[5] method instead as it provides even more clue to the local environment.

check out the source code on my github page.

Using OpenBabel to Make Input Files for GAMESS

One thing that people tend to hate to do, is to make repetitive tasks. One such thing could be to prepare 100 different geometry optimizations using DFT[1] and a happy mixture of basis sets[2] and functionals[3] for GAMESS[4]. In this blog post, I'll show you how you can get around manual labor and prepare 10 files with just a little bit of coding. Extending it to a 100 is an easy task after that. Lets start.

An input geometry of ethanol in the very common XYZ file format might look something very similar to this, when made in Avogadro[5]

9
ethanol.xyz
C  0.00004  0.00001  0.00002
C  0.00035 -0.00106  1.52091
H  0.02794  1.04324 -0.38022
H -0.92039 -0.49540 -0.37524
H  0.88619 -0.55034 -0.38151
H -0.89666  0.55109  1.87908
O  1.16332  0.62487  1.98684
H -0.03205 -1.05436  1.87773
H  1.11207  0.59582  2.97727

Using Open Babel[6] we can quickly make a GAMESS input file using the following command

babel -ixyz ethanol.xyz -ogamin ethanol.inp

resulting in

 $CONTRL COORD=CART UNITS=ANGS $END
 $DATA
ethanol.xyz
C1
C 6.0  0.0000400000  0.0000100000  0.0000200000 
C 6.0  0.0003500000 -0.0010600000  1.5209100000 
H 1.0  0.0279400000  1.0432400000 -0.3802200000 
H 1.0 -0.9203900000 -0.4954000000 -0.3752400000 
H 1.0  0.8861900000 -0.5503400000 -0.3815100000 
H 1.0 -0.8966600000  0.5510900000  1.8790800000 
O 8.0  1.1633200000  0.6248700000  1.9868400000 
H 1.0 -0.0320500000 -1.0543600000  1.8777300000 
H 1.0  1.1120700000  0.5958200000  2.9772700000 
 $END

but we can't run this file right away since it is obviously lacking information about the basis set, the DFT-functional and even the runtype (GAMESS assumes runtyp=energy if nothing is specified).

Ideally what we want is to generate a header for the input file and then insert it into the input file, but this header should be generated automatically for us.

Given a set of parameters (functional and basis set), the following python script (header.py) prints a correct header for GAMESS to the screen

import sys

basissets = {'sto3g' :'STO NGAUSS=3',
             'pc0' :'PC0',
             'pc1' :'PC1',
             'apc0' :'APC0',
             'apc1' :'APC1'}

functionals = {'blyp' :'BLYP',
               'b3lyp': 'B3LYP'}

basis = basissets[sys.argv[1]]
functional = functionals[sys.argv[2]]
print " $SYSTEM MWORDS=125 $END"
print " $CONTRL RUNTYP=OPTIMIZE DFTTYP=%s ISPHER=1 $END" % (functional)
print " $BASIS GBASIS=%s $END" % (basis)

Notice how we use a dictionary to store the basis set information so it can be easily recovered and printed correctly. The GAMESS basis group can be quite hideous with Pople basis sets[7].

All we require now is a small bash script, which makes loops over the various basis sets and functionals, invoke the python script we just made with the correct arguments and dump that printout information to a file (header.inp) and then use some Open Babel magic (to add a custom header with keywords, use the very secret -xf filename flag when converting to GAMESS input files) to make it work. The following script does exactly that

#!/usr/bin/env bash
for basis in sto3g pc0 pc1 apc0 apc1
do
  for functional in blyp b3lyp
  do
    filename=$basis"_"$functional".inp"
    python header.py $basis $functional > header.inp
    babel -xf header.inp -ixyz ethanol.xyz -ogamin $filename
  done
done

For information on how to use GAMESS, I recommend Professor Jan Jensens blog (http://molecularmodelingbasics.blogspot.com/) as a great resource and for general usage notes of GAMESS and the Google user group (http://groups.google.com/group/gamess) is very helpful if you have any problems or questions related to running GAMESS.

NB! To discover what options the babel file conversion tool supports for your file format (gamin in this case), simply invoke the following command

babel -H gamin

Happy Conversion!

First Post!

Alright,

I've tried it numerous times before. This time, I'll try it again.

This blog will contain blog posts on using python in chemistry related fields such as:

• Theoretical And Computational Chemistry
• Molecular Modeling
• Cheminformatics
• Teaching Chemistry

I'll also divulge into areas related to the above

Hopefully, there will be more posts soon ...