Ehrenfest dynamics (TDDFT/MD)

For a brief introduction to the Ehrenfest dynamics theory and the details of its implementation in GPAW, see Ehrenfest theory. The original implementation by Ari Ojanpera is described in Ref. [1].

See also

• Time-propagation TDDFT

• gpaw.tddft.TDDFT

• gpaw.tddft.TDDFT.propagate()

Ground state

Similar to static TDDFT calculations, one has to start with a standard ground state simulation. In TDDFT, one can use larger grid spacing than for geometry optimization, so for example, if you use h=0.25 for geometry optimization, try h=0.3 for TDDFT to save a lot of time (the same spacing should be used for the ground state and TDDFT calculators).

Ground state example:

from ase import Atoms
from gpaw import GPAW

name = 'h2_diss'

# Create H2 molecule in the center of a box aligned along the z axis.
d_bond = 0.754  # H2 equilibrium bond length
atoms = Atoms('H2', positions=[(0, 0, 0), (0, 0, d_bond)])
atoms.set_pbc(False)
atoms.center(vacuum=4.0)

# Set groundstate calculator and get and save wavefunctions
calc = GPAW(mode='fd', h=0.3, nbands=1, basis='dzp', txt=name + '_gs.txt',
            symmetry={'point_group': False})
atoms.calc = calc
atoms.get_potential_energy()
calc.write(name + '_gs.gpw', mode='all')

Simulating H2 dissociation

We then simulate Ehrenfest dynamics of the H2 molecule in a very intense laser field to observe its dissociation.

For Ehrenfest dynamics we must use the parameter propagator='EFSICN' for the TDDFT calculator tdcalc to take into account the necessary corrections in the propagation of the time-dependent Kohn-Sham equation. The parameter solver='BiCGStab' to use the stabilized biconjugate gradient method (BiCGStab) is generally recommended for Ehrenfest dynamics calculations for any system more complex than simple dimers or dimers.

H2 dissociation example:

from gpaw.tddft import TDDFT
from gpaw.tddft.ehrenfest import EhrenfestVelocityVerlet
from gpaw.tddft.laser import CWField
from ase.units import Hartree, Bohr, AUT
from ase.io import Trajectory
from ase.parallel import parprint

name = 'h2_diss'

# Ehrenfest simulation parameters
timestep = 10.0  # timestep given in attoseconds
ndiv = 10  # write trajectory every 10 timesteps
niter = 500  # run for 500 timesteps

# TDDFT calculator with an external potential emulating an intense
# harmonic laser field aligned (CWField uses by default the z axis)
# along the H2 molecular axis.
tdcalc = TDDFT(name + '_gs.gpw',
               txt=name + '_td.txt',
               propagator='EFSICN',
               solver='BiCGStab',
               td_potential=CWField(1000 * Hartree, 1 * AUT, 10))

# For Ehrenfest dynamics, we use this object for the Velocity Verlet dynamics.
ehrenfest = EhrenfestVelocityVerlet(tdcalc)

# Trajectory to save the dynamics.
traj = Trajectory(name + '_td.traj', 'w', tdcalc.get_atoms())

# Propagates the dynamics for niter timesteps.
for i in range(1, niter + 1):
    ehrenfest.propagate(timestep)

    if tdcalc.atoms.get_distance(0, 1) > 2.0:
        # Stop simulation if H-H distance is greater than 2 A
        parprint('Dissociated!')
        break

    # Every ndiv timesteps, save an image in the trajectory file.
    if i % ndiv == 0:
        # Currently needed with Ehrenfest dynamics to save energy,
        # forces and velocitites.
        epot = tdcalc.get_td_energy() * Hartree
        F_av = ehrenfest.F * Hartree / Bohr  # forces
        v_av = ehrenfest.v * Bohr / AUT  # velocities
        atoms = tdcalc.atoms.copy()
        # Needed to save the velocities to the trajectory:
        atoms.set_velocities(v_av)

        traj.write(atoms, energy=epot, forces=F_av)

traj.close()

The distance between the H atoms at the end of the dynamics is more than 2 Å and thus their bond was broken by the intense laser field.

Electronic stopping in graphene

A more complex use for Ehrenfest dynamics is to simulate the irradiation of materials with either chared ions (Ref. [2]) or neutral atoms (Ref. [3]).

The following script calculates the ground state of the projectile + target system, with the parameter charge defining its charge state. For ionisation state +1, an external potential is used at the hydrogen ion the converge the calculation. One might also have to change the default convergence parameters depending on the projectile used, and to verify the convergence of the results with respect to the timestep and k-points. Here, slightly less strict criteria are used. 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).

Projectile + target example:

import numpy as np

from ase.lattice.hexagonal import Graphene
from gpaw import GPAW
from ase import Atoms
from ase.parallel import parprint
from gpaw.utilities import h2gpts
from ase.units import Bohr
from gpaw.external import ConstantPotential
from gpaw import RMMDIIS

name = 'graphene_h'

# 5 x 5 supercell of graphene
index1 = 5
index2 = 5
a = 2.45
c = 3.355

gra = Graphene(symbol='C',
               latticeconstant={'a': a, 'c': c},
               size=(index1, index2, 1),
               pbc=(1, 1, 0))

gra.center(vacuum=15.0, axis=2)
gra.center()

# Starting position of the projectile with an impact point at the
# center of a hexagon
projpos = [[gra[15].position[0], gra[15].position[1] + 1.41245, 25.0]]

H = Atoms('H', cell=gra.cell, positions=projpos)

# Combine target and projectile
atoms = gra + H
atoms.set_pbc(True)

# Specify the charge state of the projectile, either
# 0 (neutral) or 1 (singly charged ion)
charge = 0  # default for neutral projectile

# Two-stage convergence is needed for the charged system
conv_fast = {'energy': 1.0, 'density': 1.0, 'eigenstates': 1.0}
conv_par = {'energy': 0.001, 'density': 1e-3, 'eigenstates': 1e-7}

if charge == 0:
    calc = GPAW(mode='fd',
                gpts=h2gpts(0.2, gra.get_cell(), idiv=8),
                nbands=110,
                xc='LDA',
                txt=name + '_gs.txt',
                convergence=conv_par,
                symmetry={'point_group': False})

    atoms.calc = calc
    atoms.get_potential_energy()
    calc.write(name + '.gpw', mode='all')

if charge == 1:
    const_pot = ConstantPotential(1.0)

    calc = GPAW(mode='fd',
                gpts=h2gpts(0.2, gra.get_cell(), idiv=8),
                nbands=110,
                xc='LDA',
                charge=1,
                txt=name + '_gs.txt',
                convergence=conv_fast,
                external=const_pot,
                symmetry={'point_group': False})

    atoms.calc = calc
    atoms.get_potential_energy()

    # External potential used to prevent charge tranfer from graphene to ion.
    A = 1.0
    X0 = atoms.positions[50] / Bohr
    rcut = 3.0 / Bohr
    vext = calc.hamiltonian.vext
    gd = calc.hamiltonian.finegd
    n_c = gd.n_c
    h_c = gd.get_grid_spacings()
    b_c = gd.beg_c
    vext.vext_g.flags.writeable = True
    vext.vext_g[:] = 0.0
    for i in range(n_c[0]):
        for j in range(n_c[1]):
            for k in range(n_c[2]):
                x = h_c[0] * (b_c[0] + i)
                y = h_c[1] * (b_c[1] + j)
                z = h_c[2] * (b_c[2] + k)
                X = np.array([x, y, z])
                dist = np.linalg.norm(X - X0)
                if dist < rcut:
                    vext.vext_g[i, j, k] += A * np.exp(-dist**2)

    atoms.calc = calc.new(convergence=conv_par, eigensolver=RMMDIIS(5),
                          external=vext, txt=name + '_vext_gs.txt')

    atoms.get_potential_energy()
    atoms.calc.write(name + '.gpw', mode='all')

else:
    parprint('Charge should be 0 or 1!')

Finally, the following script can be used for performing an electronic stopping calculation for a hydrogen atom impacting graphene with the initial velocity being 40 keV. In the charged state, the external potential is automatically set to zero when the TDDFT object is initialized and hence does not affect the calculation. The calculation ends when the distance between the projectile and the bottom of the supercell is less than 5 Å. (Note: this is a fairly demanding calculation with 32 cores and requires close to 50 GB of memory.)

Electronic stopping example:

from ase.units import _amu, _me, Bohr, AUT, Hartree
from gpaw.tddft import TDDFT
from gpaw.tddft.ehrenfest import EhrenfestVelocityVerlet
from ase.io import Trajectory
import numpy as np

name = 'graphene_h'
Ekin = 40e3  # kinetic energy of the ion (in eV)

# Adapted to the ion energy; here 4 as (probably too large!)
timestep = 16.0 * np.sqrt(10e3 / Ekin)
ekin_str = '_ek' + str(int(Ekin / 1000)) + 'k'
strbody = name + ekin_str
traj_file = strbody + '.traj'

# The parallelization options should match the number of cores, here 24.
p_bands = 2  # number of bands to parallelise over
dom_dc = (2, 2, 3)  # domain decomposition for parallelization
parallel = {'band': p_bands, 'domain': dom_dc}

tdcalc = TDDFT(name + '.gpw',
               propagator='EFSICN',
               solver='BiCGStab',
               txt=strbody + '_td.txt',
               parallel=parallel)

proj_idx = 50  # atomic index of the projectile
delta_stop = 5.0 / Bohr  # stop condition when ion is within 5 A of boundary.

# Setting the initial velocity according to the kinetic energy.
amu_to_aumass = _amu / _me
Mproj = tdcalc.atoms.get_masses()[proj_idx] * amu_to_aumass
Ekin *= 1 / Hartree
v = np.zeros((proj_idx + 1, 3))
v[proj_idx, 2] = -np.sqrt((2 * Ekin) / Mproj) * Bohr / AUT
tdcalc.atoms.set_velocities(v)

evv = EhrenfestVelocityVerlet(tdcalc)
traj = Trajectory(traj_file, 'w', tdcalc.get_atoms())

trajdiv = 1  # number of timesteps between trajectory images
densdiv = 10  # number of timesteps between saved electron densities
niters = 100  # total number of timesteps to propagate

for i in range(niters):
    # Stopping condition when projectile z coordinate passes threshold
    if evv.x[proj_idx, 2] < delta_stop:
        tdcalc.write(strbody + '_end.gpw', mode='all')
        break

    # Saving trajectory file every trajdiv timesteps
    if i % trajdiv == 0:
        F_av = evv.F * Hartree / Bohr  # forces converted from atomic units
        v_av = evv.v * Bohr / AUT  # velocities converted from atomic units
        epot = tdcalc.get_td_energy() * Hartree  # energy
        atoms = tdcalc.get_atoms().copy()
        atoms.set_velocities(v_av)

        traj.write(atoms, energy=epot, forces=F_av)

    # Saving electron density every densdiv timesteps
    if (i != 0 and i % densdiv == 0):
        tdcalc.write(strbody + '_step' + str(i) + '.gpw')

    evv.propagate(timestep)

traj.close()

References

[1]

A. Ojanpera, V. Havu, L. Lehtovaara, M. Puska, “Nonadiabatic Ehrenfest molecular dynamics within the projector augmented-wave method”, J. Chem. Phys. 136, 144103 (2012).

[2]

A. Ojanpera, Arkady V. Krasheninnikov, M. Puska, “Electronic stopping power from first-principles calculations with account for core electron excitations and projectile ionization”, Phys. Rev. B 89, 035120 (2014).

[3]

C. Brand, M. Debiossac, T. Susi, F. Aguillon, J. Kotakoski, P. Roncin, M. Arndt, “Coherent diffraction of hydrogen through the 246 pm lattice of graphene”, New J. Phys. 21, 033004 (2019).