Commit e59c3b9e authored by Toma Susi's avatar Toma Susi

Added neutral atom option for electronic stopping and increased propagation...

Added neutral atom option for electronic stopping and increased propagation timestep to reach stopping condition within niter.
parent 67dfad7b
......@@ -61,15 +61,17 @@ Electronic stopping in graphene
-------------------------------
A more complex use for Ehrenfest dynamics is to simulate the irradiation of
materials with either chared ions or neutral atoms (see Refs. [#Ojanpera2014]_
and [#Brand2019]_).
materials with either chared ions (Ref. [#Ojanpera2014]_) or neutral atoms
(Ref. [#Brand2019]_).
The following script calculates the ground state of the projectile + target
system. An external potential is used at the hydrogen ion the converge the
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. Here, slightly less strict convergence
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
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:
......@@ -77,12 +79,13 @@ Projectile + target example:
.. literalinclude:: graphene_h_gs.py
Finally, the following script can be used for performing an electronic
stopping calculation for a hydrogen ion impacting graphene with the initial
velocity being 40 keV. 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 3 Å. (Note: this is a fairly
demanding calculation even with 32 cores.)
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 8 cores and requires close to 50 GB of
memory.)
Electronic stopping example:
......
......@@ -2,18 +2,15 @@ import ase.io as io
import numpy as np
from ase.lattice.hexagonal import Graphene
from gpaw.mixer import Mixer, MixerSum
from gpaw import GPAW
from ase import Atom, Atoms
from ase.parallel import parprint
from gpaw.utilities import h2gpts
from ase.units import Bohr
from gpaw.occupations import FermiDirac
from gpaw.external import ConstantPotential
from gpaw.mpi import world
def gaussian(x, x0, A):
E = np.linalg.norm(x-x0)
return A*np.exp(-E**2)
from gpaw import RMMDIIS
name = 'graphene_h'
......@@ -38,40 +35,58 @@ H = Atoms('H', cell=gra.cell, positions=projpos)
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}
const_pot = ConstantPotential(1.0)
mixer= Mixer(0.1,5,weight=100.0)
calc = GPAW(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)
atoms.set_calculator(calc)
atoms.get_potential_energy()
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]):
if charge == 0:
calc = GPAW(gpts=h2gpts(0.2, gra.get_cell(), idiv=8),
nbands = 110, xc='LDA', txt=name + '_gs.txt',
convergence=conv_par)
atoms.set_calculator(calc)
atoms.get_potential_energy()
calc.write(name + '.gpw', mode='all')
if charge == 1:
const_pot = ConstantPotential(1.0)
calc = GPAW(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)
atoms.set_calculator(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] += gaussian(X,X0,A)
vext.vext_g[i,j,k] += A*np.exp(-np.linalg.norm(X-X0)**2)
calc.set(convergence=conv_par, eigensolver=RMMDIIS(5), external=vext)
calc.set(convergence=conv_par, eigensolver=RMMDIIS(5), external=vext)
atoms.get_potential_energy()
calc.write(name + '.gpw', mode='all')
atoms.get_potential_energy()
calc.write(name + '.gpw', mode='all')
else:
parprint('Charge should be 0 or 1!')
......@@ -8,30 +8,27 @@ from gpaw.mpi import world
import numpy as np
name = 'graphene_h'
Ekin = 40e3
timestep = 1.0 * np.sqrt(10e3/Ekin)
Ekin = 40e3 # Kinetic energy of the ion (in eV)
timestep = 8.0 * np.sqrt(10e3/Ekin) # Adapted to the ion energy; here 0.5 as (may be too large!)
ekin_str = '_ek' + str(int(Ekin/1000)) + 'k'
amu_to_aumass = _amu/_me
strbody = name + ekin_str
traj_file = strbody + '.traj'
# The parallelization options should match the number of cores, here 32.
# The parallelization options should match the number of cores, here 8.
p_bands = 2 # Number of bands to parallelise over
dom_dc = (4,4,1) # Domain decomposition for parallelization
p_bands = 1 # Number of bands to parallelise over
dom_dc = (2,2,1) # 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
v = np.zeros((proj_idx+1,3))
delta_stop = 3.0 / Bohr
Mproj = tdcalc.atoms.get_masses()[proj_idx]
Ekin *= Mproj
Ekin = Ekin / Hartree
proj_idx = 50 # Atomic index of the projectile
delta_stop = 5.0 / Bohr # Stop condition when ion is within 5 A of cell boundary.
Mproj *= amu_to_aumass
# Setting the initial velocity according to the kinetic energy.
Mproj = tdcalc.atoms.get_masses()[proj_idx] * amu_to_aumass
Ekin *= Mproj / Hartree
v = np.zeros((proj_idx+1,3))
v[proj_idx,2] = -np.sqrt((2*Ekin)/Mproj) * Bohr / AUT
tdcalc.atoms.set_velocities(v)
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment