Simulating ion bombardment on graphene sheets

The properties of graphene sheets can be tuned by intentionally introducing defects to the material. A promising method to obtain defects in a controlled way is by bombarding the graphene sheet with high-energy ions. Molecular dynamics simulations can be used to elucidate the mechanisms involved in such processes and to improve our understanding of how external parameters, such as the kinetic energy of the incoming ion, affect the formation of defects. In this tutorial, you simulate ion bombardment of graphene, following the protocol proposed in Ref. [BS12].

In particular, you will:

  • Learn the basic steps required to set up the calculation using Virtual NanoLab (VNL).
  • Learn how to manually introduce modifications to the Python script generated by VNL.
  • Run the calculation using the Atomistix ToolKit (ATK).

Setting up the Graphene Sheet:

Open the builder_icon Builder and do the following:

  • Use the Nanosheet plugin via Add ‣ From plugin ‣ Nanosheet to create a graphene unit cell.
  • Use the Bulk Tools ‣ Swap Axes tool to rotate the structure and the cell vectors by clicking C<->A and then Z<->X. Now the surface is oriented along the Z-axis.
  • Use Bulk Tools ‣ Repeat to create a suitable supercell by giving a desired number of repetitions in the A- and B-directions, e.g. 10x4.
  • Use the Bulk Tools ‣ Lattice parameters widget to increase the z-component of the C-vector to 50.0 Å. Make sure to keep the cartesian coordinates when changing the lattice.

Adding a Bombardment Atom

Bombardment atoms can in principle be added at any location above the surface. In the following, you can shoose to add a bombardment atom either above a graphene carbon atom (top site), or above the center of a graphene ring (hollow site).

  • To add an atom at the top site:

    • Select an atom in the center of the graphene sheet.
    • Use Coordinate Tools ‣ Translate to place the bombardment atom at the desired starting height above the graphene sheet (e.g. 30 Å), by entering a suitable z-component for the translation vector.
    • Check Copy and click Apply.
  • To add an atom in a hollow site,:

    • Select the atoms belonging to a carbon ring in the center of the sheet.
    • Click the icon AlignmentPoint02a_icon in the tool bar of the Builder.
    • Translate the atom to the desired starting position above the surface by selecting the created atom and use the Coordinate Tools ‣ Translate tool without checking the Copy box.


If you want to modify the exact position of the bombardment atom, you can do this easily, either by using the Translate tool, or by editing the coordinates of the last atom (using the Coordinate Tools ‣ Coordinate List tool).

Setting up the Simulation

You will now set up the calculator for the simulation using the script_generator_icon Script Generator. It is recommended to optimize the lattice of the graphene sheet with the classical potential before you start the simulations. Details about the calculator and the classical potential are explained below.

  • First, send the structure from the Builder to the Script Generator by clicking the sendto_icon icon .

In the Script Generator:

  • Add a calculator_icon New Calculator block.
  • Add two optimization_icon MolecularDynamics blocks.
  • Open the New Calculator widget. Select the ATK-ForceField calculator with the Tersoff_SiC_1998 potential.
  • Choose a suitable output filename in the IO section, and then click OK.


The Tersoff_SiC_1998 potential is very similar, but not identical, to the potential used in [BS12]. If you feel like using the exact same potential, you can find it in this script:


For ATK-versions older than 2017, the ATK-ForceField calculator can be found under the name ATK-Classical.

  • In the first optimization_icon MolecularDynamics block, set up the equilibration simulation with the following parameters.

The position of the bombardment atom should be fixed during the equilibration calculation. Therefore:

  • Click Add Constraints to open the Constraints widget.
  • Select the bombardment atom above the surface and click Add tag from Selection.
  • Set the constraints field to Fixed and click OK.

In the second optimization_icon MolecularDynamics block, you should set the basic parameters of the bombardment simulation:

  • Select the NVE Velocity Verlet type.
  • Set 50000 steps and a log interval of 100.
  • Enter a suitable filename for the trajectory.
  • Set the initial velocities to Configuration velocities to use the velocities after equilibration.
  • Uncheck the Remove center-of-mass momentum box.
  • Decrease the time step to 0.05 fs.


The kinetic energy of the atom is expected to be very high in this simulation. That is why a relatively small time step is necessary to guarantee the convergence. For very high incident energies (>100 eV), the time step may have to be decreased even further to 0.01 fs (cf. [BS12]).

  • Send the script to the editor_icon Editor, using the sendto_icon icon.

Modifying the Script

To let the bombardment atom approach the surface at the desired kinetic energy, a few modifications have to be made in the script.

  • Right before the second optimization_icon MolecularDynamics block, insert the following lines.
# Get the velocities from the bulk_configuration
velocities = bulk_configuration.velocities()
# Define the incident energy of the bombardment atom.
incident_energy = 10.0*eV
# Calculate the corresponding velocity and apply it to the last atom.
incident_velocity = ((2.0*incident_energy/Carbon.atomicMass())**0.5).inUnitsOf(Ang/fs)
velocities[-1, :] = [0.0, 0.0, -incident_velocity]*Ang/fs


The velocity vector is set to be perpendicular to the surface in this calculation. Arbitrary incident direction can be considered by modifying the velocities vectors in the script lines.

  • Save the script under a suitable filename, and run the simulation via the job_manager_icon JobManager or in a terminal.
  • After the simulation has completed you can visualize the trajectory using the MovieTool or the viewer_icon Viewer to inspect the damage.


[BS12](1, 2, 3) Edson P. Bellido and Jorge M. Seminario. Molecular dynamics simulations of ion-bombarded graphene. J. Phys. Chem. C, 116(6):4044–4049, 2012. doi:10.1021/jp208049t.