My blog has moved!

My blog has moved to here

Replaying a Lammps Trajectory

EDIT: Turns out I just didn't google the correct words, LAMMPS does have one. It's the command rerun. Oops.

LAMMPS, the molecular dynamics engine, doesn't have a replay function allowing one to go back and calculate properties on a trajectory, especially forces and energies. This came up in my work and I got started on writing some new LAMMPS code to accomplish this. I wanted to be able to calculate forces with a different potential on a trajectory. Anyway, I realized you can accomplish this just by using the LAMMPS looping variables. No new code necessary. It can be done by using the read_dump command to load particular trajectory frame. A sample input script is below

Predicting the Future of Simulations

In my day job, molecular modeling, I simulate the position and velocity of every atom in biological molecules. As you can imagine, doing this yields nearly complete information about whatever biological system I'm studying and allows me to predict the behavior of proteins. Well, that's the theory. It's more complex than that and often fails spectacularly in practice. Ignore that for now. Let's focus on the systems that can be simulated. In 1982 CHARMM, one the first molecular dynamics tool capable of simulating tiny proteins, was released. It was possible to simulate a small protein for one tenth of one billionth of one second. Now we can simulate a protein for one thousandth of a second. That's an increase of 10^7 in the speed of molecular simulations in about 30 years. In 1990, we were able to simulate lysozyme, a small protein, for few billionths of a seconds. Now we can simulate a virus for a few billionths of a second. That's an increase in about 10^4 mass of molecular simulations in 30 years. Assuming continued exponential growth (see plots below), here are my upcoming interesting milestones in computer simulation

Year	Event
2023	The first simulation of every atom in a ribosome for 10 milliseconds
2038	The first simulation of every atom in a virus for 1 second
2131	The first simulation of every atom of a cell for 1 hour
2191	We can simulate every atom of every cell in a heart for 1 heartbeat
2227	Every atom of every cell in a human brain can be simulated as fast they move in a human brain

Radial Distribution Functions of Water from Experiment

Site-site radial distribution functions, sometimes called g(r) for short, are one of the most important measured quantities in molecular simulations. As usual, wikipedia has a nice overview of them. Probably the most famous non-trivial g(r) is the water oxygen-oxygen radial distribution function. It's often taught in statistical mechanics courses that it's possible to compare a simulation g(r) with X-ray or neutron-diffraction experiments with a simple Fourier transform. In fact, this is almost never the case due to sources of error in experiments. Check out the Introduction to this article by AK Soper for more info on the kinds of error which prevent the simple Fourier transform from converting experimental structure factors to a g(r).

In any case, the point of this post is to put the raw data for an experimental g(r) for water. The AK Soper article uses the empirical potential structure refinement method for converting experimental neutron diffraction data at 298 K 1 bar into the three radial distribution functions for water (O-O, O-H, H-H). This is the best method I know of for obtaining these radial distribution functions and this is the data used to test the TIP4P/2005 water model. I've extracted the data from their paper and made it available here. Make sure to cite their paper if you use the data in a publication. Hopefully these g(r)s will help you in creating your own new special water models! See the figure below the break:

Indexing an Emacs Macro

This post is about how to increment and paste a counter while executing an emacs macro. For some reason, this is not easy to find on the Internet. So I'll write it down on my blog since I keep forgetting.

1. Set the start point of the counter
   M-: (setq i 1)
2. Start recording macro with
   C-x (
3. Insert index as needed with
   C-u M-: i
4. Increment index when finished with macro
   M-: (setq i (+ i 1))
5. Complete macro with
   C-x )
6. Execute Macro with
   C-x e

To execute the macro multiple times, use ALT-(# of times) C-x e.

Extending Radial Distributions

Radial distribution functions are one of the most important quantities to calculate when running a molecular dynamics simulation. They give good insight into the phase and pair potentials between molecules. Practically, they are calculated according to the following equation:

Cube-Sphere Intersection Volume

This week's math problem was particularly difficult. I wanted to know the intersection volume between a cube and a sphere. Imagine a sphere contained in a cube. The sphere will first touch the cube when its radius is equal to half the width of the cube. Until that point, its volume is just the formula for the volume of a sphere. Once the sphere's radius is 2 sqrt(3) times the width of the cube, then it completely engulfs the cube and the intersection of the 2 objects is just the volume of a cube. In between, I've found it looks like this:

Bacteria in Blender

I made this while drinking my morning coffee. I just stretched a UV Sphere, stuck two layers of hair on it (flagella, cilia) and added a few modifiers to get it bumpy. The texture is just some glass, Voronoi textures, and diffuse shaders. Meh, it's OK.

Generating n-Dimensional Lattice Coordinates

It comes up sometimes in my research that I need to generate coordinates along a lattice. For example, when generating an initial structure for a molecular dynamics simulation. I thought I share my code for this. The code works by recursively enumerating each dimension in a lattice and calling a given function each time it reaches a point. So, it could be used to execute any function along a lattice. Here's what it looks like for a 2D and 3D example. In the 3D example I didn't choose a cubic number of points on the lattice, so that it sort of stops in the middle.

Category Colors in R

I thought I'd write down the color palette I use for categories. I've seen many discussions about gradient color choices, but category colors are less often discussed. I have a function given below which generates colors where there is not meant to be an ordering to the color. It is important when data points are near, that they are easy to distinguish with the eye. It makes more sense to see it:

ggplot Style in Classic R Plots

ggplot is a highly acclaimed R package for plotting. I have had very little experience with the library because I've mostly memorized all the quirks of normal R plots. I like the ggplot default style though, so I thought write down how to replicate it in R plots. The distinctive feature of ggplot is the gray rectangle and white grid lines. This may be replicated like so:



Here's an example

DNA-Polymer Nanoparticle

DNA-Polymer nanoparticle are used to deliver DNA to cell nuclei for curing diseases at the genetic level. In my research group, we work with carboxy-betaine acrylamide polymers. I modeled an exceptionally large DNA-polymer nanoparticle below. The polymer is 80 kDa and the DNA is a single strand of 65 base pairs. After condensing it (in blender, not a true physical simulation), the nanoparticle comes out to about 15nm in diameter. 

The Accuracy of Ghose-Crippen ALogP for Peptides

I was using the Ghose-Crippen ALogP predictor implemented in CDK as a QSAR for looking at peptide solubilities. Unfortunately, this partition coefficient predictor is fit on small organic molecules and does a pretty bad job at predicting peptide partition coefficients. Here's what it looks like on single amino acids compared with experimental data.

Visualizing Size in Biochemical Systems - Update

After some helpful suggestions, I've revised the previous image I made. This time I've used the atom representation for the proteins and peptides, so atoms are used everywhere for visualization and they have the correct radii. Doesn't look as cool, but it is more consistent.

Visualizing Size in Biochemical Systems

I decided to waste a few hours and make a Blender image based on my post on the relative size of common elements in biochemical systems. Using 1 blender unit = 1 nm, I modeled a gold nanoparticle, an antibody, albumin, an iron oxide nanoparticle, a 10 amino acid length peptide, and salicylic acid.

Relative Size of Proteins, Nanoparticles, Antibodies, Peptides, and Drugs

Needing to know the relative size of these elements in biochemical systems often comes up in my research, so I've decided to write them.

Component	Size [nm]	Reference
Gold Nanoparticle	~20	Can vary, most common is 20. Frens (1973)
Iron Oxide Nanoparticle	5-15	Can be much larger with emulsion methods. This size is most common when magentic properties are desired. Called partial reduction coprecipitation method.
Antibody	15.5 x 11	PDB: 1igt measured in VMD
Human Serum Albumin	7.6 x 7	PDB: 1e7i measured in VMD
10 Amino Acid Peptide	~1.5 (coiled globule sphere)	Paper I wrote on Peptide Conformations
Salicylic Acid (drug)	0.6 x 0.5	PubChem structure measured in VMD

Approximate Restricted Integer Partition With Exact Part Number

I'm interested in a certain type of mathematical structure this week: integer partitions. An integer partition, p(n), is the set of sets of integers that add up to n. For example,

• p(3) = 1 + 1 + 1, 1 + 2, 3
• p(4) = 1 + 1 + 1 + 1, 1 + 1 + 2, 2 + 2, ...

Speeding up R

I've been programming in R for a while now. Whenever I've had performance problems, it almost always is due to data.frame usage. To check what's slowing down your R code, just use the Rprof command like so:

Changing Image Size in Tachyon in VMD

I always forget how to change the image size for rending molecular images in VMD using the (external) Tachyon rendering engine, so I thought I'd write down the steps here:

In the TCL console, get the aspect ratio correct for the image. So, for example to make a nice 1920x1080 image, we'll make the view at half of that:

 
display resize 960 540

Now, add simply add -res 1920 1080 to the end of the render command in the file render controls window. Here's an example of one of my render commands:

 
"$HOME/VMD/tachyon_LINUXAMD64" -aasamples 12 %s -format TARGA -o %s.tga  -res 1920 1080

Here's the result of one my renders:

Making nanoparticles in Blender

Today I made some nanoparticles in Blender. I used the node compositor to create the fluorophore effect. Making science in Blender is fun!