molecular gear published in Nature Materials

Step-by-step rotation of a molecule-gear mounted on an atomic-scale axis
C. Manzano1*,W.-H. Soe1*, H. S.Wong, F. Ample, A. Gourdon, N. Chandrasekhar and C. Joachim1,2

compared to the small bearing (disassembled) of old:

riding the rails

MATLAB function to calculate the energy of double wall carbon nanotubes

function [U dr dt]=DWNT(n1, m1, n2, m2)
% this function uses the chiral vectors of two carbon nanotubes as inputs
%(n1 m1)-inner & (n2 m2)-outer and returns the energy of the tube in kcal
%mol^-1 atom^-1 using an empirical equation found in
%"Double-Wall carbon nanotubes...",
%Bakalis & Zerbetto, Chemical Physics Letters 463 (2008) 139-140
%function written by Tom Moore

%define energy constants
a= 1665.51; %kcal A^6 mol^-1 atom^-1
b= -92.85; %kcal A^3 mol^-1 atom^-1
c= 2.50; %kcal A^1 mol^-1 atom^-1
d= 17.24 ; %kcal A^4 mol^-1 atom^-1
f= -58.56; %kcal A^5 mol^-1 atom^-1

%define graphene lattice constant in A
k= 2.49;

% calculate the chiral angles of both tubes in radians
t1_rad=acos((2*n1+m1)/(2*(n1^2+m1^2+n1*m1)^.5));
t2_rad=acos((2*n2+m2)/(2*(n2^2+m2^2+n2*m2)^.5));

dt=t1_rad-t2_rad %may need to do a <1 check here

% calculate the radius of both tubes
r1=(2.49*(n1^2+m1^2+n1*m1)^.5)/(2*pi);
r2=(2.49*(n2^2+m2^2+n2*m2)^.5)/(2*pi);
dr= r2-r1;

%find the energy of the tube
A= (a/dr^6)+(b/dr^3)+(c/dr);
B=sin(dt)*cos(dt);
C=(d/dr^4)+(f/dr^5);
U=A+B*C

multiwall carbon nanotube & GROMACS

Force constants:

Bonds 47890 KJ/mol/nm^2
Angle 400 KJ/mol/rad^2
Dihedral 167.36 KJ/mol/rad^2

more swnts in GROMACS (back to large image files)

The tubes look a lot better when simulated without the caps. There is none of that weird buckling as before. It is interesting to note that the buckling occurred along a line containing a pentagon defect in the lattice.

Single Wall Carbon Nanotubes in 4.0.3 GROMACS and Windows Vista Part 1

Okay, I am starting to get results that are at least approximating what I would expect, by that I mean stays together and jiggles around. I'll get down to brass tacks and lay out how I did it and then later on talk about why I went about it this way, lessons learned, what's next, all that.

1. I used the oplsaa force field for this. This force field is not included in older versions of gmx. It already contains carbon in the atom types, so we don't have to add that. However we do have to add some cnt parameters to a couple of files. These files are all located in gromacs/top. The first file to change is ffoplsaa.n2t . We need to add a couple of carbon atoms bonded to two and three other carbon atoms. Open up the file in a text editor and add them to the bottom:

C opls_240 0 12.011 3 C 0.142 C 0.142 C 0.142


C opls_239 0 12.011 2 C 0.142 C 0.142

You can probably guess what's going on here: atom type, charge, mass, # of bonds, what they're bonded to and bond length in nm.

Here is my copy of the file: ffoplsaa.n2t

2. Next we are going to add some bond parameters to ffoplsaabon.itp. I just dropped these in at the bottom of the file where I saw other people had made some prior modifications:

; Added by TEM for CNTs


[ bondtypes ]

C C 1 0.14210 478900


[ angletypes ]

C C C 1 120.000 397.480


[ dihedraltypes ]

C C 1 0.000 167.360 1


Note there is some weirdness happening on the one side of the above tube. I'm not quite sure the dihedraltypes are correct, may need to revise that. Again here is my file: ffoplsaabon.itp

Okay, those were the changes. If you had read my previous post, don't worry about messing with the .rtp file. That was a messy fix and not necessary.

Next fire up cygwin and cd your directory to where you have you cnt pdb file. If you want a copy of my nano-capsule, here it is: swnt_capped.pdb

in cygwin, type: editconf -f swnt_capped.pdb -o swnt_capped.gro -d 2.0

This will create your .gro file with a 2.0 nm box around the tube.

next type: x2top -f swnt_capped.gro -o swnt_capped.top

This will create the topology file. The documentation says you can use a pdb file as input however I never got it to work with anything but the .gro file.

Next you need a simulation parameter file, md.mdp. Here is the one I used to get you started:
md.mdp

Put this file in your directory with your pdb, .gro and .top files. If you copy the text into a text file and save it as a text file, it usually doesn't work; copy an actual .mdp file and change the contents. So far the only thing I feel comfortable messing around with in this file is the simulation length and time step, lot going on here.

Now type: grompp -f md.mdp -c swnt_capped.gro -p swnt_capped.top -maxwarn 3

The experienced user may be interested to know that I used to have to set the maximum number of warnings to 150 to get past this...

And now: mdrun -s topol.tpr

The sim takes a couple of minutes to run through. You can then load the .gro and . trr file into VMD to view the results, make movies, whatever.

There is still more to go, need to figure out what is causing that buckling, a very knowledgeable professor has told me the potentials look too soft, bunch of stuff. If you have something to add to this, let me know.

Special thanks to Brian and Damian for putting up with the play-by-play emails, hope I'm not on the spam list now guys.

step closer to GROMACS and Carbon nanotubes, or some ramblings

Okay I figured out why I was getting all the short bond errors when I ran pdb2gmx (see previous post). I came across this paragraph in the gmx wiki:

"One point that is not made clearly in the manual is how to specify the connectivity between residues. This mechanism is protein-specific, and relies on the existence in each residue of atoms with the same names, i.e. the protein backbone atoms. The [ bonds ] section in the example in section 5.6.1 of the manual defines a bond between -C and N. This refers to the C atom in the previous residue. Similarly, one might use a +N prefix to refer to the N atom in the following residue."

http://wiki.gromacs.org/index.php/.rtp_file#.rtp_Residue_Topology

This meant I never defined how the residues were connected together in my .rtp file. This explains why I got a bunch of short bond warnings and my topology file never had pairs, angles or dihedrals.

I've changed the [bonds] section of the .rtp from the cnt carbons to:

NT1 ]
[atoms]
C C 0 0

C C 0 0

[ bonds ]
C +C 0.14210 478900

However I can't be certain that adding a +C to the second C in all the bonds makes sense for all of them. I guess I need to know how you put together a large protein and then figure out how you assemble a carbon nanotube using the same basic steps...