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.
|
Table of Contents
|
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:
1
ffcnt_oplsaa
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 topol.top -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:
periodic_molecules=yes
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:
Hi,
Can you please answer this question.
How Gromacs differentiate same atom pairs and bond length having different charges from info in Pdb files? For example?
For example OH of carboxyl group and Phenolic OH group has same atoms and bond length but slightly different charges.
Thanks for providing the procedure of making topology file of CNT. Its a great help. By the way, Can I use this procedure to make topology file for ZnO slab of 4 layers which contains at around 1500 atoms?
I have tried but getting error again and again. Please give me some suggestions.
Regards
Sasthi
Post preview:
Close preview