# ab initio molecular dynamics (DFT/MD)¶

Ab initio molecular dynamics uses DFT to calculate the forces between the atoms at each time step. While computationally expensive, prohibiting simulations longer than a few ps, DFT/MD can be directly applied to any system that DFT can describe.

## Transmission of hydrogen through graphene¶

Neutral atom transmission through a graphene target can be simulated with DFT/MD. However, note that TDDFT/MD is required to account for electronic stopping (Ref. ).

The following script simulates the impact of a hydrogen atom with an initial velocity corresponding to a kinetic energy of 40 keV, transmitting through the center of a hexagon in a graphene target. The impact point in this case is the center of a carbon hexagon, but this can be modified by changing the x-y position of the H atom (`projpos`).

In a realistic calculation, one should consider the choice of the XC functional, and verify the convergence of the results with respect to the grid spacing, timestep and k-points. Here, less accurate parameters are used so that the calculation can be run on 4-8 CPU cores.

```import ase.units as units
from ase import Atoms
from ase.build import graphene
from ase.io import Trajectory
from ase.md.verlet import VelocityVerlet
from gpaw import GPAW
from gpaw.utilities import h2gpts

name = 'graphene_h'

# 5 x 5 supercell of graphene
a = 2.45
gra = graphene(a=a, size=(5, 5, 1), vacuum=10)
gra.center()

# Starting position of the projectile with an impact point at the
# center of a hexagon.
# Set mass to one atomic mass unit to avoid isotope average.
atoms = gra + Atoms('H', masses=[1.0])
d = a / 3**0.5
atoms.positions[-1] = atoms.positions + (0, d, 5)
atoms.pbc = (True, True, True)

calc = GPAW(gpts=h2gpts(0.2, gra.get_cell(), idiv=8),
nbands=110,
xc='LDA',
txt=f'{name}_gs.txt')

atoms.calc = calc
atoms.get_potential_energy()

# Moving to the MD part
ekin = 100  # kinetic energy of the ion (in eV)
timestep = 0.1  # timestep in fs

# Filename for saving trajectory
ekin_str = '_ek' + str(int(ekin / 1000)) + 'keV'
strbody = name + ekin_str
traj_file = f'{name}_ek_{ekin}.traj'

# Integrator for the equations of motion, timestep depends on system
dyn = VelocityVerlet(atoms, timestep * units.fs)

# Saving the positions of all atoms after every time step
with Trajectory(traj_file, 'w', atoms) as traj:
dyn.attach(traj.write, interval=1)

# Running one timestep before impact
dyn.run(1)

# Giving the target atom a kinetic energy of ene in the -z direction
atoms[-1].momentum = -(2 * ekin * atoms[-1].mass)**0.5

# Running the simulation for 80 timesteps
dyn.run(80)
```