Modeling Carbon Nanotubes with GROMACS

Online there are many documents about GROMACS and Carbon Nanotube (CNTs) (1, 2, 3, ...). Even after lots of reading I didn't find what I was looking for: a clear and simple guide to model single- or multi-wall carbon nanotubes (SWNT, MWNT), either periodic and non-periodic.
It taken me a while to find a procedure that does the things as I wanted: build a periodic CNT and create the topology for GROMACS.

Build and prepare a CNT

The first step, obviously, consists in building the Carbon Nanotube, for example with YASC buildCstruct. With this script it is possible to build singlewall and multiwalls CNT, both armchair or zigzag, finite or periodic.
The output file is a TINKER xyz file that can be edited to XYZ and then converted in pdb, for example using the ChemAxon Marvin Beans utility molconvert or openbabel or saved directly in pdb using vmd. Since version 1.1, buildCstruct is able to save the structure directly in gromacs GRO format.

In order to use x2top, if you have the structure in pdb format, you must convert it into a .gro file with editconf. Also, the size of the box should be specified (mind that in gromacs the length for the vectors a, b and c is in nm, not angstrom):

editconf -f file.pdb -o file.gro -box a b c -angles 90 90 90

The gro file will then contain the CNT centered in the periodic box that has been specified. Note that if you use YASC buildCstruct to build the nanotube, it will be periodic along the b vector (Y-axis).

Be careful to have a unit cell big enough to contain all the CNT, or you won't be able to generate the topology. You can check this with gromacs ngmx or molden or using pbc in vmd. Also, mind to choose a box bigger enough to contain also the other molecules you may want to put in contact with the CNT, e.g. polymer chains.

Files for x2top

x2top requires a bunch of files in order to properly create the topology file. I wanted to use the oplsaa force field, but I didn't want to mess up the original ffoplsaa files in the share/gromacs/top directory, therefore I have created from scratch the files I needed.

.n2t itp

This file is used to translate name in topology. It reads the name of an atom in your structure and, based on its connectivity (i.e., the number of bonds, the bond lengths and the atoms bonded), tries to guess the proper atom type among those defining in the .n2t file.
My .n2t for graphene and CNT based on the oplsaa force field is shown below:

; Oplsaa-based n2t for carbon-based structures such as CNTs and graphenes
; Andrea Minoia
H    HJ    0.00       1.008  1    C 0.109                              ;Hydrogen
C    CJ    0.00      12.011  3    C 0.142   H 0.109   H 0.109 ;Periferic C
C    CJ    0.00      12.011  3    C 0.142   C 0.142   H 0.108 ;Periferic C
C    CJ    0.00      12.011  1    C 0.142                          ;Internal/periodic C
C    CJ    0.00      12.011  2    C 0.142   C 0.142            ;Internal/periodic C
C    CJ    0.00      12.011  3    C 0.142   C 0.142   C 0.142 ;Internal/periodic C

I have decided to put all the charges to zero, but you may want to do differently… it's up to you :)

.itp file

I prefer not to let x2top put the force field parameters in the topology files, and I have defined my own .itp file:

; Oplsaa-based force field for carbon-based structures such as CNTs and graphenes
; Andrea Minoia

[ defaults ]
; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1               3               yes             0.5     0.5
; parameters are taken from the OPLS force field

[ atomtypes ]
; The charges here will be overwritten by those in the rtp file
; name       mass      charge    ptype      sigma      eps
  CJ   6     12.01100     0.000       A    3.55000e-01  2.92880e-01 ;opls_147 naftalene fusion C9
  HJ   1      1.00800     0.000       A    2.42000e-01  1.25520e-01 ;opls_146 HA hydrogen benzene. I have set the charges zero

[ bondtypes ]
; i    j func        b0          kb
  CJ    CJ      1    0.14000   392459.2   ; TRP,TYR,PHE
  CJ    HJ      1    0.10800   307105.6   ; PHE, etc.

[ angletypes ]
  CJ     CJ     CJ      1   120.000    527.184   ; PHE(OL)
  CJ     CJ     HJ      1   120.000    292.880   ;
  HJ     CJ     HJ      1   117.000    292.880   ; wlj from HC-CM-HC

[ dihedraltypes ]
  CJ     CJ     CJ     CJ      3     30.33400   0.00000 -30.33400   0.00000   0.00000   0.00000 ; aromatic ring
  HJ     CJ     CJ     HJ      3     30.33400   0.00000 -30.33400   0.00000   0.00000   0.00000 ; aromatic ring
  HJ     CJ     CJ     CJ      3     30.33400   0.00000 -30.33400   0.00000   0.00000   0.00000 ; aromatic ring

No dihedrals are used at the moment.

.rtp file

Finally, I needed to create a residue topology file for oplsaa, in which I did not define any residue:

; New format introduced in Gromacs 3.1.4.
; Dont use this forcefield with earlier versions.

; Oplsaa-based rtp for carbon-based structures such as CNTs and graphenes
; Andrea Minoia

; NB: OPLS chargegroups are not strictly neutral, since we mainly
; use them to optimize the neighborsearching. For accurate simulations
; you should use PME.

[ bondedtypes ]
; Col 1: Type of bond
; Col 2: Type of angles
; Col 3: Type of proper dihedrals
; Col 4: Type of improper dihedrals
; Col 5: Generate all dihedrals if 1, only heavy atoms of 0.
; Col 6: Number of excluded neighbors for nonbonded interactions
; Col 7: Generate 1,4 interactions between pairs of hydrogens if 1
; Col 8: Remove propers over the same bond as an improper if it is 1
; bonds  angles  dihedrals  impropers all_dihedrals nrexcl HH14 RemoveDih
     1       1          3          1        0         3      1     0

FF.dat file

Here I specify the name of the force field:


Generate the topology with x2top

Once all the files are ready, it is time to generate the topology using x2top:

x2top -f file.gro -o -ff cnt_oplsaa -name CNT -noparam -pbc

The option -pbc will generate all bonds, angles and dihedral that cross the periodic boundary conditions (this is why in the .gro file we need to specify the good size of the cell). If your nanotube is not periodic, then no need to go for -pbc.

Correct the topology

The last step consists in change from 1 to 3 the function used to describe the dihedrals.

.mdp file

Turns on the periodicity in your mdp file using:

pbc = xyz

That's all… now you are ready to model your periodic CNT or carbon-based structure.

Solvating CNT

To solvate a CNT we can use the gromacs utility genbox with the options -cp (the nanotube gro file) and -cs (solvent box gro file).
Usually, I previously equilibrate a solvent box in NPT and then create a supercell in order to have a box bigger than the cell of the CNT. This because the box of the solute (-cp) will be the box of the final system.

genbox will try to fill all the empty space with solvent molecules, even inside the CNT where is possible. If this is not desirable, we can write a script to create an index file that does not contain the atomic indices for the atoms in the molecules inside the CNT. Then, use editconf to generate the final system:

editconf -f starting_system.gro -n index.ndx -o final_system.gro

Remove atoms inside a carbon nanotube with YASC nosolincnt

The index file to use with editconf in order to remove the molecules of solvent from inside the CNT, can be easily written by YASC nosolincnt, a tcl script for VMD.
After loading the system into VMD, open the TCL console and execute the script:

source /path/of/the/script/nosolincnt.tcl

You will be prompted to insert the fragment number of the CNT (play with the representations if you are not sure about this) and the CNT radius in Angstrom. Those info are used to create a cylindrical selection of the given radius taken from the X-,Z-coords of the CNT center of mass. You can check the result of this selection by pressing C when the view window is active. If everything is good, then press W and the index file, without the selected atoms, will be written.

In the photo gallery below, a solvated CNT with solvent molecules included in it, the molecules to remove and the final system are displayed, respectively.

A video tutorial to illustrate the procedure is here.

Complex case: solvating CNT+polymer

I don't really know what interest could have model a solvated CNT alone, surely, if we add a polymer to the game, things will get more interesting. Below the flowchart I have found to solvate a CNT interacting with a polymer chain.


That's it.

Troubleshoots, tips and tricks

Grompp and non-maching atom names

If you had generated the gro file starting by a pdb file, depending on how you have generated the pdb file, residue names could be missing. In this case I have experienced a series of warning with grompp:

Warning: atom name X in _FILE.top_ and _FILE.gro_ does not match...

I solved the problem by adding a residue name to the .gro and .top files.

VI editor

learn it! it is your best buddy, in particular when working with gromacs and huge text files ;)

Rate this page:

rating: 0+x

Add a New Comment
or Sign in as Wikidot user
(will not be published)
- +
Unless otherwise stated, the content of this page is licensed under Creative Commons Attribution-ShareAlike 3.0 License