Commit 18aeaca3 authored by Tuomas Rossi's avatar Tuomas Rossi

Rewrite Ag55 tutorial

Simplify, use pvalence basis sets, use realistic parameters, and
split the workflow.
parent dc800ae2
......@@ -107,6 +107,8 @@ for calculating many properties during the propagation.
See :ref:`analysis` for tutorial how to extend the analysis capabilities.
.. _note basis sets:
General notes about basis sets
==============================
......@@ -182,7 +184,7 @@ ScaLAPACK and large band parallelization with ``augment_grids`` enabled.
This can be achieved with parallelization settings like
``parallel={'sl_auto': True, 'domain': 8, 'augment_grids': True}``
(see :ref:`manual_parallel`),
which would use 8 tasks for domain parallelization and the rest for
which will use groups of 8 tasks for domain parallelization and the rest for
band parallelization (for example, with a total of 144 cores this would mean
domain and band parallelizations of 8 and 18, respectively).
......@@ -402,6 +404,37 @@ Advanced tutorials
Plasmon resonance of silver cluster
-----------------------------------
This tutorial demonstrates
:ref:`the efficient parallelization options <parallelization>` and
the importance of :ref:`proper basis sets <note basis sets>`.
We calculate the photoabsorption spectrum of
:download:`an icosahedral Ag55 cluster <lcaotddft_Ag55/Ag55.xyz>`.
We use GLLB-SC potential to significantly improve the description of d states,
:ref:`pvalence basis sets` to improve the description of unoccupied states, and
11-electron Ag setup to reduce computational cost.
The workflow is the same as in the previous examples.
First, we calculate ground state (takes around 20 minutes with 36 cores):
.. literalinclude:: lcaotddft_Ag55/gs.py
Then, we carry the time propagation for 30 femtoseconds in steps of
10 attoseconds (takes around 3.5 hours with 36 cores):
.. literalinclude:: lcaotddft_Ag55/td.py
Finally, we calculate the spectrum (takes a few seconds):
.. literalinclude:: lcaotddft_Ag55/spec.py
The resulting spectrum shows an emerging plasmon resonance at around 4 eV:
.. image:: lcaotddft_Ag55/Ag55_spec.png
For further insight on plasmon resonance in metal nanoparticles,
see [#Kuisma2015]_ and [#Rossi2017]_.
One should think about what type of transitions of interest are present,
and make sure that the basis set can represent such Kohn-Sham electron and
hole wave functions. The first transitions in silver clusters will be
......@@ -415,20 +448,7 @@ is **better** than the default double-zeta polarized one.
We will use the 11-electron Ag setup, since the semi-core p states included
in the default setup are not relevant here.
.. literalinclude:: lcaotddft_Ag55/basis.py
We calculate :download:`an icosahedral Ag55 cluster <lcaotddft_Ag55/Ag55.xyz>`.
This code uses ScaLAPACK parallelization with 48 cores.
.. literalinclude:: lcaotddft_Ag55/ag55.py
Code runs for approximately two hours wall-clock time.
The resulting spectrum shows already emerging plasmonic excitation
around 4 eV.
For more details, see [#Kuisma2015]_.
.. image:: lcaotddft_Ag55/fig1.png
.. TODO
......
from ase.io import read
from gpaw import GPAW
from gpaw import PoissonSolver
from gpaw.occupations import FermiDirac
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
from gpaw.tddft.spectrum import photoabsorption_spectrum
from gpaw import setup_paths
# Set the path for the created basis set
setup_paths.insert(0, '.')
# Read the clusterm from the xyz file
atoms = read('ag55.xyz')
atoms.center(vacuum=5.0)
# Increase the accuracy of density for ground state
convergence = {'density': 1e-8}
# Increase the accuracy of PoissonSolver and
# apply multipole corrections for monopole and dipoles
poissonsolver = PoissonSolver(eps=1e-16, remove_moment=1 + 3)
# Calculate ground state in LCAO mode
calc = GPAW(xc='GLLBSC', basis='GLLBSC.dz', h=0.3, nbands=352, mode='lcao',
convergence=convergence, poissonsolver=poissonsolver,
occupations=FermiDirac(0.1),
parallel={'sl_default': (8, 6, 32), 'band': 2})
atoms.set_calculator(calc)
# Relax the ground state
atoms.get_potential_energy()
# Save the intermediate ground state result to a file
calc.write('ag55_gs.gpw', mode='all')
# Time propagation
td_calc = LCAOTDDFT('ag55_gs.gpw',
parallel={'sl_default': (8, 6, 32), 'band': 2})
DipoleMomentWriter(td_calc, 'ag55.dm')
td_calc.absorption_kick([1e-5, 0.0, 0.0])
td_calc.propagate(20, 500)
photoabsorption_spectrum('ag55.dm', 'ag55.spec', width=0.2)
# Creates: fig1.png
import numpy as np
import matplotlib.pyplot as plt
plt.figure(figsize=(6, 6 / 2 ** 0.5))
data = np.loadtxt('ag55.spec')
plt.plot(data[:, 0], data[:, 1], 'k')
plt.title(r'Absorption spectrum of $Ag_{55}$ with GLLB-SC potential')
plt.xlabel('Energy (eV)')
plt.ylabel('Absorption (arbitrary units)')
plt.ylim(ymin=0)
plt.xlim(0, 10)
plt.savefig('fig1.png')
from ase.io import read
from gpaw import GPAW, FermiDirac, Mixer, PoissonSolver
# Read the structure from the xyz file
atoms = read('Ag55.xyz')
atoms.center(vacuum=6.0)
# Increase the accuracy of density for ground state
convergence = {'density': 1e-12}
# Use occupation smearing and weak mixer to facilitate convergence
occupations = FermiDirac(0.02)
mixer = Mixer(0.02, 5, 1.0)
# Parallelzation settings
parallel = {'sl_auto': True, 'domain': 2, 'augment_grids': True}
# Increase the accuracy of PoissonSolver and
# apply multipole corrections for monopole and dipoles
poissonsolver = PoissonSolver(eps=1e-16, remove_moment=1 + 3)
# Ground-state calculation
calc = GPAW(mode='lcao', xc='GLLBSC', h=0.3, nbands=360,
setups={'Ag': '11'},
basis={'Ag': 'pvalence.dz', 'default': 'dpz'},
convergence=convergence, poissonsolver=poissonsolver,
occupations=occupations, mixer=mixer, parallel=parallel,
txt='gs.out')
atoms.set_calculator(calc)
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
# Creates: Ag55_spec.png
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
data_ej = np.loadtxt('spec.dat')
plt.figure(figsize=(6, 6 / 1.62))
ax = plt.subplot(1, 1, 1)
ax.plot(data_ej[:, 0], data_ej[:, 1], 'k')
ax.spines['right'].set_visible(False)
ax.spines['top'].set_visible(False)
ax.yaxis.set_ticks_position('left')
ax.xaxis.set_ticks_position('bottom')
plt.title(r'Absorption spectrum of Ag$_{55}$ with GLLB-SC potential')
plt.xlabel('Energy (eV)')
plt.ylabel('Photoabsorption (eV$^{-1}$)')
plt.xlim(0, 6)
plt.ylim(ymin=0)
plt.tight_layout()
plt.savefig('Ag55_spec.png')
......@@ -4,5 +4,8 @@ from myqueue.task import task
def create_tasks():
return [
task('[email protected]:10m'),
task('[email protected]:2h', deps='basis.py'),
task('fig1.py', deps='ag55.py'),
task('[email protected]:30m'),
task('[email protected]:4h', deps='gs.py'),
task('[email protected]:1m', deps='td.py'),
task('[email protected]:1m', deps='spec.py'),
]
from gpaw.tddft.spectrum import photoabsorption_spectrum
# Calculate spectrum
photoabsorption_spectrum('dm.dat', 'spec.dat', width=0.1, delta_e=0.01)
from gpaw.lcaotddft import LCAOTDDFT
from gpaw.lcaotddft.dipolemomentwriter import DipoleMomentWriter
# Parallelzation settings
parallel = {'sl_auto': True, 'domain': 2, 'augment_grids': True}
# Time propagation
td_calc = LCAOTDDFT('gs.gpw', parallel=parallel, txt='td.out')
DipoleMomentWriter(td_calc, 'dm.dat')
td_calc.absorption_kick([1e-5, 0.0, 0.0])
td_calc.propagate(10, 3000)
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