Commit cca8b902 authored by Jens Jørgen Mortensen's avatar Jens Jørgen Mortensen

Merge branch 'fix-lcaotddft-doc' into 'master'

Update LCAOTDDFT documentation

See merge request !611
parents 3d1f1bbd 7ba9bc39
Pipeline #103759899 passed with stage
in 3 minutes and 43 seconds
......@@ -3,7 +3,4 @@ from myqueue.task import task
def create_tasks():
return [
task('[email protected]:10m'),
task('[email protected]:2h', deps='lcaotddft_basis.py'),
task('lcaotddft_fig1.py', deps='lcaotddft_ag55.py'),
task('[email protected]:40m')]
......@@ -14,7 +14,7 @@ from gpaw import GPAW
calc = GPAW(mode='lcao', h=0.3, basis='dzp',
setups={'Na': '1'},
poissonsolver=poissonsolver,
convergence={'density': 1e-8})
convergence={'density': 1e-12})
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
......
......@@ -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
==============================
......@@ -177,16 +179,24 @@ benchmark their suitability for your application**.
Parallelization
===============
LCAO-TDDFT is parallelized using ScaLAPACK. It runs without ScaLAPACK,
but in this case only a single core is used for linear alrebra.
For maximum performance on large systems, it is advicable to use
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 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).
Instead of ``sl_auto``, the ScaLAPACK settings can be set by hand
as ``sl_default=(m, n, block)`` (see :ref:`manual_ScaLAPACK`,
in which case it is important that ``m * n``` equals
the total number of cores used by the calculator
and that ``max(m, n) * block < nbands``.
* Use ``parallel={'sl_default':(N, M, 64)}``; See :ref:`manual_parallel`.
* It is necessary that N*M equals the total number of cores used
by the calculator, and ``max(N,M)*64 < nbands``, where ``64`` is the used
block size. The block size can be changed to, e.g., 16 if necessary.
* Apart from parallelization of linear algrebra, normal domain and
band parallelizations can be used. As in ground-state LCAO calculations,
use band parallelization to reduce memory consumption.
It is also possible to run the code without ScaLAPACK, but it is
very inefficient for large systems as in that case only a single core
is used for linear algebra.
.. TODO
......@@ -394,33 +404,68 @@ Advanced tutorials
Plasmon resonance of silver cluster
-----------------------------------
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
`5s \rightarrow 5p` like. We require 5p orbitals in the basis set, and thus,
we must generate a custom basis set.
In this tutorial, we demonstrate the use of
:ref:`efficient parallelization settings <parallelization>` and
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.
**When calculating other systems, remember to check the convergence
with respect to the used basis sets.**
Recall :ref:`hints here <note basis sets>`.
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]_.
User-generated basis sets
-------------------------
Here is how to generate a double-zeta basis set with 5p orbital in valence
for silver for GLLB-SC potential. Note that the extra 5p valence state
effectively improves on the ordinary polarization function, so this basis set
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.
The :ref:`pvalence basis sets` distributed with GPAW and used in
the above tutorial have been generated from atomic PBE orbitals.
Similar basis sets can be generated based on atomic GLLB-SC orbitals:
.. literalinclude:: lcaotddft_basis.py
.. literalinclude:: lcaotddft_Ag55/mybasis/basis.py
We calculate the icosahedral Ag55 cluster: :download:`ag55.xyz`
The Ag55 cluster can be calculated as in the above tutorial, once
the ground-state calculation has been modified to use
the generated setup and basis set:
This code uses ScaLAPACK parallelization with 48 cores.
.. literalinclude:: lcaotddft_Ag55/mybasis/gs.py
:diff: lcaotddft_Ag55/gs.py
.. literalinclude:: lcaotddft_ag55.py
The calculation with this generated "my" p-valence basis set results only in
small differences in the spectrum in comparison to
the distributed :ref:`pvalence basis sets`:
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/Ag55_spec_basis.png
.. image:: fig1.png
The spectrum with the default dzp basis sets is also shown for reference,
resulting in unconverged spectrum due to the lack of diffuse functions.
**This demonstrates the importance of checking convergence
with respect to the used basis sets.**
Recall :ref:`hints here <note basis sets>` and
see [#Kuisma2015]_ and [#Rossi2015]_ for further discussion on the basis sets.
.. TODO
......
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(25e-3)
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': 'dzp', 'default': 'dpz'},
convergence=convergence, poissonsolver=poissonsolver,
occupations=occupations, mixer=mixer, parallel=parallel,
maxiter=1000,
txt='gs.out')
atoms.set_calculator(calc)
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
../spec.py
\ No newline at end of file
../td.py
\ No newline at end of file
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(25e-3)
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,
maxiter=1000,
txt='gs.out')
atoms.set_calculator(calc)
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
from gpaw.atom.generator import Generator
from gpaw.atom.basis import BasisMaker
from gpaw.atom.configurations import parameters, parameters_extra
xc = 'GLLBSC'
name = 'my'
if 'Ag' in parameters_extra:
args = parameters_extra['Ag'] # Choose the 11-electron setup
else:
args = parameters['Ag'] # Choose the 17-electron setup
args.update(dict(name=name, use_restart_file=False, exx=True))
# Generate setup
generator = Generator('Ag', xc, scalarrel=True)
generator.run(write_xml=True, **args)
# Generate basis
bm = BasisMaker('Ag', name='{}.{}'.format(name, xc), xc=xc, run=False)
bm.generator.run(write_xml=False, **args)
basis = bm.generate(zetacount=2, polarizationcount=0,
jvalues=[0, 1, 2], # include d, s and p
)
basis.write_xml()
from ase.io import read
from gpaw import GPAW, FermiDirac, Mixer, PoissonSolver
from gpaw import setup_paths
# Insert the path to the created basis set
setup_paths.insert(0, '.')
# 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(25e-3)
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': 'my'},
basis={'Ag': 'GLLBSC.dz', 'default': 'dpz'},
convergence=convergence, poissonsolver=poissonsolver,
occupations=occupations, mixer=mixer, parallel=parallel,
maxiter=1000,
txt='gs.out')
atoms.set_calculator(calc)
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
../spec.py
\ No newline at end of file
../td.py
\ No newline at end of file
# 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')
# Creates: Ag55_spec_basis.png
# -*- coding: utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
data_ej = np.loadtxt('spec.dat')
data_my_ej = np.loadtxt('mybasis/spec.dat')
data_dzp_ej = np.loadtxt('dzp/spec.dat')
plt.figure(figsize=(6, 6 / 1.62))
ax = plt.subplot(1, 1, 1)
ax.plot(data_ej[:, 0], data_ej[:, 1], 'k', label='p-valence')
ax.plot(data_my_ej[:, 0], data_my_ej[:, 1], 'k--', label='p-valence (my)')
ax.plot(data_dzp_ej[:, 0], data_dzp_ej[:, 1], 'k:', label='dzp')
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'Effect of LCAO basis on the spectrum')
plt.xlabel('Energy (eV)')
plt.legend(loc='upper left', title='LCAO basis')
plt.ylabel('Photoabsorption (eV$^{-1}$)')
plt.xlim(0, 6)
plt.ylim(ymin=0)
plt.tight_layout()
plt.savefig('Ag55_spec_basis.png')
from myqueue.task import task
def create_tasks():
gs = task('[email protected]:30m'),
td = task('[email protected]:4h', deps=[gs]),
spec = task('[email protected]:1m', deps=[td]),
plot = task('[email protected]:1m', deps=[spec]),
basis_my = task('[email protected]:10m', folder='mybasis'),
gs_my = task('[email protected]:30m', folder='mybasis', deps=[basis_my]),
td_my = task('[email protected]:4h', folder='mybasis', deps=[gs_my]),
spec_my = task('[email protected]:1m', folder='mybasis', deps=[td_my]),
gs_dzp = task('[email protected]:30m', folder='dzp'),
td_dzp = task('[email protected]:4h', folder='dzp', deps=[gs_dzp]),
spec_dzp = task('[email protected]:1m', folder='dzp', deps=[td_dzp]),
basis_plot = task('[email protected]:1m',
deps=[spec, spec_my, spec_dzp]),
return [gs, td, spec, plot,
basis_my, gs_my, td_my, spec_my,
gs_dzp, td_dzp, spec_dzp,
basis_plot]
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)
......@@ -20,7 +20,7 @@ ps = ExtraVacuumPoissonSolver(gpts=(512, 256, 256),
calc = GPAW(mode='lcao', h=0.3, basis='pvalence.dz', xc='LDA', nbands=6,
setups={'Na': '1'},
poissonsolver=ps,
convergence={'density': 1e-8},
convergence={'density': 1e-12},
txt='gs.out')
atoms.set_calculator(calc)
energy = atoms.get_potential_energy()
......
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)
import matplotlib.pyplot as plt
from ase.io import read
from gpaw.lcaotddft.tddfpt import transform_local_operator
transform_local_operator(gpw_file='Na8_gs.gpw',
tdop_file='Na8.TdDen',
fqop_file='Na8.FqDen',
omega=1.8,
eta=0.23)
dct = read('Na8.FqDen.imag.cube', full_output=True)
data = dct['data'][:, :, 16]
atoms = dct['atoms']
extent = [0, atoms.cell[0][0], 0, atoms.cell[1][1]]
plt.imshow(data.T, origin='lower', extent=extent)
for atom in atoms:
circle = plt.Circle((atom.position[0], atom.position[1]),
0.3, color='r', clip_on=False)
fig = plt.gcf()
fig.gca().add_artist(circle)
plt.title('Induced density of $Na_{8}$')
plt.xlabel('$\\AA$')
plt.ylabel('$\\AA$')
plt.savefig('Na8_imag.png')
from gpaw.atom.generator import Generator
from gpaw.atom.basis import BasisMaker
args = {'core': '[Kr]', 'rcut': 2.45}
generator = Generator('Ag', 'GLLBSC', scalarrel=True)
generator.N *= 2 # Increase grid resolution
generator.run(**args)
bm = BasisMaker(generator, name='GLLBSC', run=False)
basis = bm.generate(zetacount=2, polarizationcount=0,
energysplit=0.07,
jvalues=[0, 1, 2], # include d, s and p
rcutmax=12.0)
basis.write_xml()
# 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 import Atoms
from gpaw.tddft import photoabsorption_spectrum
from gpaw import PoissonSolver
from gpaw.lcaotddft.tddfpt import DensityCollector
from gpaw.lcaotddft import LCAOTDDFT
atoms = Atoms('Na8', positions=[[i * 3.0, 0, 0] for i in range(8)])
atoms.center(vacuum=5.0)
# Calculate all bands
td_calc = LCAOTDDFT(
basis='dzp', setups={'Na': '1'}, xc='LDA', h=0.3, nbands=4,
convergence={'density': 1e-7},
poissonsolver=PoissonSolver(eps=1e-14, remove_moment=1 + 3 + 5))
atoms.set_calculator(td_calc)
atoms.get_potential_energy()
td_calc.write('Na8_gs.gpw', mode='all')
td_calc = LCAOTDDFT('Na8_gs.gpw')
td_calc.attach(DensityCollector('Na8.TdDen', td_calc))
td_calc.absorption_kick([1e-4, 0.0, 0.0])
td_calc.propagate(20, 1000, 'Na8.dm')
photoabsorption_spectrum('Na8.dm', 'Na8.spec', width=0.15)
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