Biomodeling: graphene
Intro
We will do 4_graphene.
Graphene is a promising material due to its interesting physical and chemical properties such as electrical conductivity and mechanical durability.
The classical MD simulations can handle the mechanics and chemistry of graphene.
Definitely, electrical conductivity requires QM treatments.
A typical application of graphene is the nanopore instrument.
For example, Ref 1 shows the application for osmosis (or desalination) using graphene (CNT to be precise).
In this tutorial, we will learn:
How to prepare the MD simulation setup of graphene systems.
How to calculate the ionic conductance through a nanopore.
References
Kalra, A., Garde, S. & Hummer, G. Osmotic water transport through carbon nanotube membranes. Proc Natl Acad Sci U S A 100, 10175-10180 (2003). https://doi.org:10.1073/pnas.1633354100
Hello, graphene!
As usual, we need a graphene version of
Structure file (PDB)
and, an ITP file including moleculetype matching the PDB file.
Because graphene is not officially included in conventional force fields, it is not so easy to prepare those.
To make the process convenient, we created a web server
Options – Dimension of graphene
zigzag side along the x axis
armchair side along the y axis.
Options – Periodic boundary conditions
In many cases, we assume that the graphene sheet is very large.
To emulate such large graphene sheet using nm-scale box in MD, we can use PBC.
Exactly same logic as the lipid membranes.
Options – Hole
We can put a nanopre in the graphene to induce an ionic current through the graphene sheet.
Options – Terminal groups
In most MD simulations of graphene, people use blunt-ended carbon atoms, which are chemically incorrect.
Carbon atoms at the edge should form chemical bonds with some oxygen or hydrogen atoms.
Perform MD
The web server gives you a package including most files necessary for MD.
Downloaded output from the web server includes:
amber99sb-ildn-phi-bsc0.ff: force field files that will be included in topol.top
topol.top: topology file in Gromacs format.
final.pdb: the PDB file of generated graphene.
graphene.itp: moleculetype matching final.pdb
graphene.posres.itp: position restraints file if you want to fix graphene in space.
md.???: output files from a short MD simulation on the server. final.pdb is the output of this MD run.
Load final.pdb to VMD and see how it looks like.
See periodic images.
If you chose PBC on the web server, than moleculetype definition will have covalent bonds between carbon atoms at the edges.
Solvation
Confirm the box size in final.pdb. Adjust it using text editor as you like.
CRYST1 48.854 50.777 69.904 90.00 90.00 90.00 P 1 1
Note that if you are using PBC, you cannot modify the x and y dimensions.
Add water molecules using gmx solvate
gmx solvate -cp final.pdb -o conf.pdb -cs -p topol.top
-cp: Input PDB
-o: Output PDB with water molecules.
-cs: Just tells that you want to add water
-p: The number of added water molecules will be specified at the end of topol.top file.
Generate tpr file to use as an input for gmx genion
gmx grompp -f mini.mdp -maxwarn 40 -c conf.pdb -r conf.pdb
Add ions using gmx genion
gmx genion -conc 1 -s topol.tpr -o conf.pdb -p topol.top
-conc: the concentration of NaCl in M
-s: input TPR file
-o: output PDB
-p: The number of added ions will be specified at the end of topol.top file.
Minimization
You can use mini.mdp downloaded from the web server.
gmx grompp -f mini.mdp -r conf.pdb -c conf.pdb -maxwarn 40
gmx mdrun -c conf.pdb -nt 4
-c: output PDB
md.mdp – Pressure coupling
You can use md.mdp downloaded from the web server.
Mostly same options as other simulations.
Pressure coupling should be anisotropic because graphene sheet is anisotropic.
; pressure coupling
Pcoupl = berendsen
Pcoupltype = anisotropic
; Time constant (ps), compressibility (1/bar) and reference P (bar)
tau-p = 1.0
compressibility = 4.5e-5 4.5e-5 4.5e-5
ref-p = 1 1 1
md.mdp – Electric field
We want to apply a constant electric field along the z axis such that ions move through the pore.
The height of the box in the z axis is 7 nm for this example (could be different for you).
I want to apply a potential difference of 1 V.
Then, the electric field will be 1 V / 7 nm = 0.142857 V/nm.
Add the following lines to md.mdp
; Electric fields
; Format for electric-field-x, etc. is: four real variables:
; amplitude (V/nm), frequency omega (1/ps), time for the pulse peak (ps),
; and sigma (ps) width of the pulse. Omega = 0 means static field,
; sigma = 0 means no pulse, leaving the field to be a cosine function.
electric-field-x = 0 0 0 0
electric-field-y = 0 0 0 0
electric-field-z = 0.142857 0 0 0
Perform MD
gmx grompp -c conf.pdb -f md.mdp -r conf.pdb -maxwarn 40
gmx mdrun
Run for about 10 ns.
Load the trajectory to VMD and see ions moving through the pore.
The rate of ion movements is the ionic current.
Visualizing ionic flow
Ionic flows can be visualized by using the densflux application on the server.
We need to upload a PDB file and an XTC file that includes the trajectory under an electric field.
To reduce the file size, we recommend uploading an XTC file with only ions extracted.
Write XTC file with only ions:
gmx trjconv -f traj_comp.xtc -o tmp.xtc -pbc mol -ur compact -center -boxcenter zero -s topol.tpr
-f: input XTC file
-o: output XTC file (upload this to the server).
-pbc mol -ur compact: make graphene whole when it is broken due to PBC
-center -boxcenter zero: center the coordinates such that the center is located at the origin (0,0,0). The analysis program assumes this.
For "Select group for centering", choose "Other". This is the group for graphene.
For "Select group for output", choose Ion. This includes both NA and CL.
Write PDB file with only ions:
gmx trjconv -f conf.pdb -o tmp.pdb -center -boxcenter zero -s topol.tpr
Do exactly the same as above.
Upload tmp.pdb and tmp.xtc as shown in Figure.
You need to specify the time step in your XTC file. In the downloaded file, it is set to 10 ps.
After uploading is done, you need to select group that you want to analyze.
"NA" indicates sodium ions. Put "NA" and "check". Then, "Analysis".
Attached image is under an extremely high voltage, 10 V.
Colored streamlines indicate the flux of Na ions.
Scalebar unit of flux is nm-2 ns-1.
Gray background heatmap indicates the local concentration in molar unit, M.
Homework
Try to reproduce the CNT system by the Hummer group (Ref. 1)