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:¶
- Use the Nanosheet plugin via to create a graphene unit cell.
- Use the 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 to create a suitable supercell by giving a desired number of repetitions in the A- and B-directions, e.g. 10x4.
- Use the 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 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,:
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 thetool).
Setting up the Simulation¶
You will now set up the calculator for the simulation using the 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.
In the Script Generator:
- Add a New Calculator block.
- Add two 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 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.
- 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]).
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.
1 2 3 4 5 6 7 8
# 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 bulk_configuration.setVelocities(velocities)
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 JobManager or in a terminal.
- After the simulation has completed you can visualize the trajectory using the MovieTool or the Viewer to inspect the damage.