Linear response TDDFT 2 - indexed version¶
Ground state¶
The linear response TDDFT calculation needs a converged ground state calculation with a set of unoccupied states.
We demonstrate the code usage for
(R)-methyloxirane molecule.
First, we calculate the ground state:
from ase.io import read
from gpaw import GPAW
atoms = read('r-methyloxirane.xyz')
atoms.center(vacuum=8)
calc = GPAW(h=0.2,
nbands=14,
xc='LDA',
poissonsolver={'remove_moment': 1 + 3 + 5},
convergence={'bands': 'occupied'},
txt='gs.out')
atoms.calc = calc
atoms.get_potential_energy()
calc.write('gs.gpw', mode='all')
Text output: gs.out.
Then, we converge unoccupied states:
from gpaw import GPAW
calc = GPAW('gs.gpw')
calc = calc.fixed_density(nbands=42,
convergence={'bands': 40},
maxiter=1000,
txt='unocc.out')
calc.write('unocc.gpw', mode='all')
Text output: unocc.out.
Let’s have a look on the states in the output file:
Band Eigenvalues Occupancy
0 -27.01684 2.00000
1 -18.44358 2.00000
2 -16.21613 2.00000
3 -14.70939 2.00000
4 -12.66138 2.00000
5 -11.68439 2.00000
6 -10.60330 2.00000
7 -9.93510 2.00000
8 -9.26324 2.00000
9 -8.94170 2.00000
10 -7.44985 2.00000
11 -6.29526 2.00000
12 -0.43775 0.00000
13 0.01845 0.00000
14 0.03453 0.00000
15 0.14067 0.00000
16 0.63631 0.00000
17 0.67821 0.00000
18 0.71377 0.00000
19 0.82988 0.00000
20 0.86534 0.00000
21 0.90683 0.00000
22 1.00434 0.00000
23 1.02512 0.00000
24 1.06549 0.00000
25 1.19626 0.00000
26 1.27747 0.00000
27 1.34996 0.00000
28 1.37359 0.00000
29 1.41483 0.00000
30 1.43999 0.00000
31 1.47051 0.00000
32 1.53497 0.00000
33 1.58404 0.00000
34 1.63794 0.00000
35 1.66528 0.00000
36 1.71123 0.00000
37 1.79165 0.00000
38 1.83837 0.00000
39 1.86016 0.00000
40 1.89781 0.00000
41 1.99248 0.00000
We see that the Kohn-Sham eigenvalue difference between HOMO and the highest converged unoccupied state is about 8.2 eV. Thus, all Kohn-Sham single-particle transitions up to this energy difference can be calculated from these unoccupied states. If more is needed, then more unoccupied states would need to be converged.
Note
Converging unoccupied states in some systems may require tuning the eigensolver. See the possible options in the manual.
Calculating response matrix and spectrum¶
The next step is to calculate the response matrix with
LrTDDFT2.
A very important convergence parameter is the number of Kohn-Sham
single-particle transitions used to calculate the response matrix.
This can be set through state indices
(see the parameters of LrTDDFT2),
or as demonstrated here,
through an energy cutoff parameter max_energy_diff.
This parameter defines the maximum energy difference of
the Kohn-Sham transitions included in the calculation.
Note! If the used gpw file does not contain enough unoccupied states so
that all single-particle transitions defined by max_energy_diff
can be included, then the calculation does not usually make sense.
Thus, check carefully the states in the unoccupied states calculation
(see the example above).
Note also! The max_energy_diff parameter does not mean that
the TDDFT excitations would be converged up to this energy.
Typically, the max_energy_diff needs to be much larger than the smallest
excitation energy of interest to obtain well converged results.
Checking the convergence with respect to the number of states
included in the calculation is crucial.
In this script, we set max_energy_diff to 7 eV.
We also show how to parallelize calculation
over Kohn-Sham electron-hole (eh) pairs with
LrCommunicators
(8 tasks are used for each GPAW calculator):
from gpaw.mpi import world
from gpaw import GPAW
from gpaw.lrtddft2 import LrTDDFT2
from gpaw.lrtddft2.lr_communicators import LrCommunicators
# Maximum energy difference for Kohn-Sham transitions
# included in the calculation
max_energy_diff = 7.0 # eV
# Set parallelization
dd_size = 2 * 2 * 2
eh_size = world.size // dd_size
assert eh_size * dd_size == world.size
lr_comms = LrCommunicators(world, dd_size, eh_size)
calc = GPAW('unocc.gpw',
communicator=lr_comms.dd_comm)
lr = LrTDDFT2('lr2', calc,
fxc='LDA',
max_energy_diff=max_energy_diff,
recalculate=None, # Change this to force recalculation
lr_communicators=lr_comms,
txt='lr2_with_%05.2feV.out' % max_energy_diff)
# This is the expensive part triggering the calculation!
lr.calculate()
# Get and write spectrum
spec = lr.get_spectrum('spectrum_with_%05.2feV.dat'
% max_energy_diff,
min_energy=0,
max_energy=10,
energy_step=0.01,
width=0.1)
# Get and write transitions
trans = lr.get_transitions('transitions_with_%05.2feV.dat'
% max_energy_diff,
min_energy=0.0,
max_energy=10.0)
Text output (lr2_with_07.00eV.out)
shows the number of Kohn-Sham transitions within the set 7 eV limit:
Number of electron-hole pairs: 6
Maximum energy difference: 6.973
Calculating spectrum (2021-04-07 20:51:40.541549).
Spectrum calculated (2021-04-07 20:51:40.598142).
Calculating transitions (2021-04-07 20:51:40.600258).
Transitions calculated (2021-04-07 20:51:40.601826).
The TDDFT excitations are
(transitions_with_07.00eV.dat):
# energy osc str rot str osc str x osc str y osc str z units: eVcgs
5.90014859 0.01415298 -22.29668296 0.00454107 0.00651019 0.03140768
6.31736848 0.03467841 14.60038494 0.02842660 0.00920665 0.06640197
6.36238701 0.01720349 3.39478662 0.02010759 0.01895798 0.01254491
6.48395718 0.00716565 2.79790582 0.00023826 0.00011389 0.02114480
6.94316199 0.00763759 3.06129417 0.00220299 0.01634170 0.00436809
7.00422991 0.00557088 -1.47961764 0.01088486 0.00165336 0.00417441
Restarting and recalculating¶
The calculation can be restarted with the same scipt. As an example, here we increase the energy cutoff to 8 eV. The matrix elements calculated earlier up to 7 eV are reused, and only the missing matrix elements are calculated:
from gpaw.mpi import world
from gpaw import GPAW
from gpaw.lrtddft2 import LrTDDFT2
from gpaw.lrtddft2.lr_communicators import LrCommunicators
# Maximum energy difference for Kohn-Sham transitions
# included in the calculation
max_energy_diff = 8.0 # eV
# Set parallelization
dd_size = 2 * 2 * 2
eh_size = world.size // dd_size
assert eh_size * dd_size == world.size
lr_comms = LrCommunicators(world, dd_size, eh_size)
calc = GPAW('unocc.gpw',
communicator=lr_comms.dd_comm)
lr = LrTDDFT2('lr2', calc,
fxc='LDA',
max_energy_diff=max_energy_diff,
recalculate=None, # Change this to force recalculation
lr_communicators=lr_comms,
txt='lr2_with_%05.2feV.out' % max_energy_diff)
# This is the expensive part triggering the calculation!
lr.calculate()
# Get and write spectrum
spec = lr.get_spectrum('spectrum_with_%05.2feV.dat'
% max_energy_diff,
min_energy=0,
max_energy=10,
energy_step=0.01,
width=0.1)
# Get and write transitions
trans = lr.get_transitions('transitions_with_%05.2feV.dat'
% max_energy_diff,
min_energy=0.0,
max_energy=10.0)
Text output (lr2_with_08.00eV.out)
shows the number of Kohn-Sham transitions within the set limit:
Number of electron-hole pairs: 28
Maximum energy difference: 7.961
Calculating spectrum (2021-04-07 20:52:59.518653).
Spectrum calculated (2021-04-07 20:52:59.539593).
Calculating transitions (2021-04-07 20:52:59.544254).
Transitions calculated (2021-04-07 20:52:59.546225).
The TDDFT excitations are
(transitions_with_08.00eV.dat):
# energy osc str rot str osc str x osc str y osc str z units: eVcgs
5.89434225 0.01286016 -22.38610116 0.00467060 0.00433082 0.02957905
6.31647544 0.03382317 14.26698519 0.02727170 0.00916365 0.06503416
6.35484909 0.01351514 1.51504417 0.01876752 0.01279369 0.00898422
6.47796623 0.00636743 2.73711797 0.00037130 0.00014909 0.01858190
6.93880614 0.00685675 1.21782395 0.00109041 0.01398186 0.00549798
6.99529312 0.00416800 -1.06396322 0.00787600 0.00276062 0.00186738
7.00273395 0.00197321 -1.12516199 0.00020712 0.00560001 0.00011251
7.05405681 0.00784394 7.18685093 0.01310312 0.00970403 0.00072468
7.13526003 0.00154719 -0.79666419 0.00084143 0.00023426 0.00356589
7.16399561 0.00312080 -1.29655193 0.00469785 0.00406846 0.00059610
7.22806831 0.00171861 -0.15080148 0.00000630 0.00000010 0.00514944
7.32175725 0.01331729 9.18716103 0.00011708 0.00621724 0.03361755
7.35667273 0.00338142 -3.72164113 0.00596230 0.00002769 0.00415426
7.38311944 0.00988269 9.30214077 0.02726894 0.00020590 0.00217322
7.48875052 0.00646825 -1.75487586 0.00004130 0.01487458 0.00448887
7.49725391 0.00010539 -0.01049960 0.00013052 0.00006307 0.00012257
7.52261324 0.01096163 -7.77494030 0.00756206 0.00082819 0.02449464
7.61452934 0.00005488 -0.88697000 0.00006544 0.00000277 0.00009643
7.64617018 0.00446563 0.94471763 0.00562240 0.00001166 0.00776282
7.68148265 0.00069125 -0.59263534 0.00137820 0.00006153 0.00063400
7.72419854 0.01576408 -4.16106780 0.01281219 0.03447988 0.00000016
7.75473167 0.01454875 -3.37307638 0.01581073 0.02741790 0.00041763
7.77425369 0.01075815 11.02321627 0.00394681 0.02830611 0.00002151
7.78502206 0.01614515 -6.84503319 0.03035131 0.01033524 0.00774889
7.85456975 0.01031334 -1.88758015 0.00002036 0.02874119 0.00217847
7.90796673 0.01391711 3.10361408 0.00002103 0.04115467 0.00057563
7.95224825 0.01915342 7.02224375 0.00001059 0.05570214 0.00174753
8.01054300 0.00659157 6.37317911 0.01063885 0.00314392 0.00599194
It’s important to note that also the first excitations change in comparison to the earlier calculation with the 7 eV cutoff energy. As stated earlier, the results must be converged with respect to the cutoff energy, and typically the cutoff energy needs to be much larger than the smallest excitation energy of interest.
Here we plot the photoabsorption and rotatory strength spectra from
the data files (
spectrum_with_07.00eV.dat and
spectrum_with_08.00eV.dat):
We note that the cutoff energy (and the number of unoccupied bands) should be increased more to converge the spectra properly.
Analyzing spectrum¶
Once the response matrix has been calculated, the same input script can be used for calculating the spectrum and analyzing the transitions. But, as all the expensive calculation is done already, it’s sufficient to run the script in serial. Here is an example script for analysis without parallelization settings:
from ase.parallel import paropen
from gpaw import GPAW
from gpaw.lrtddft2 import LrTDDFT2
# Maximum energy difference for Kohn-Sham transitions
# included in the calculation
max_energy_diff = 8.0 # eV
calc = GPAW('unocc.gpw')
lr = LrTDDFT2('lr2', calc,
fxc='LDA',
max_energy_diff=max_energy_diff,
txt='-')
# Get and write spectrum
spec = lr.get_spectrum('spectrum_with_%05.2feV.dat'
% max_energy_diff,
min_energy=0,
max_energy=10,
energy_step=0.01,
width=0.1)
# Get and write transitions
trans = lr.get_transitions('transitions_with_%05.2feV.dat'
% max_energy_diff,
min_energy=0.0,
max_energy=10.0)
# Get and write transition contributions
index = 0
f2 = lr.get_transition_contributions(index_of_transition=index)
with paropen('tc_%03d_with_%05.2feV.txt'
% (index, max_energy_diff), 'w') as f:
f.write('Transition %d at %.2f eV\n' % (index, trans[0][index]))
f.write(' %5s => %5s contribution\n' % ('occ', 'unocc'))
for (ip, val) in enumerate(f2):
if (val > 1e-3):
f.write(' %5d => %5d %8.4f%%\n' %
(lr.ks_singles.kss_list[ip].occ_ind,
lr.ks_singles.kss_list[ip].unocc_ind,
val / 2. * 100))
The script produces the same spectra and transitions as above.
In addition, it demonstrates how to analyze transition contributions.
An example for the first TDDFT excitation
(tc_000_with_08.00eV.txt):
Transition 0 at 5.89 eV
occ => unocc contribution
11 => 12 99.0869%
11 => 13 0.0875%
11 => 14 0.3664%
11 => 15 0.1311%
11 => 19 0.0534%
Note that while this excitation is at about 5.9 eV, it has non-negligible
contributions from Kohn-Sham single-particle transitions above this energy.
This is generally the case, and it is basically the reason why
the max_energy_diff parameter has to be typically much higher
than the highest excitation energies of interest.
Quick reference¶
- class gpaw.lrtddft2.LrTDDFT2(basefilename, gs_calc, fxc=None, min_occ=None, max_occ=None, min_unocc=None, max_unocc=None, max_energy_diff=1000000000.0, recalculate=None, lr_communicators=None, txt='-')[source]¶
Linear response TDDFT (Casida) class with indexed K-matrix storage.
Initialize linear response TDDFT without calculating anything.
Note
Does NOT support spin polarized calculations yet.
Tip
If K_matrix file is too large and you keep running out of memory when trying to calculate spectrum or response wavefunction, you can try
split -l 100000 xxx.K_matrix.ddddddofDDDDDD xxx.K_matrix.ddddddofDDDDDD.- Parameters
basefilename – All files associated with this calculation are stored as <basefilename>.<extension>
gs_calc – Ground state calculator (if you are using eh_communicator, you need to take care that calc has suitable dd_communicator.)
fxc – Name of the exchange-correlation kernel (fxc) used in calculation. (optional)
min_occ – Index of the first occupied state to be included in the calculation. (optional)
max_occ – Index of the last occupied state (inclusive) to be included in the calculation. (optional)
min_unocc – Index of the first unoccupied state to be included in the calculation. (optional)
max_unocc – Index of the last unoccupied state (inclusive) to be included in the calculation. (optional)
max_energy_diff – Noninteracting Kohn-Sham excitations above this value are not included in the calculation. Units: eV (optional)
recalculate –
Force recalculation.’eigen’ : recalculate only eigensystem (useful for on-the-flyspectrum calculations and convergence checking)’matrix’ : recalculate matrix without solving the eigensystem’all’ : recalculate everythingNone : do not recalculate anything if not needed (default)lr_communicators – Communicators for parallelizing over electron-hole pairs (i.e., rows of K-matrix) and domain. Note that ground state calculator must have a matching (domain decomposition) communicator, which can be assured by using lr_communicators to create both communicators.
txt – Filename for text output
- calculate()[source]¶
Calculates linear response matrix and properties of KS electron-hole pairs.
This is called implicitly by get_spectrum, get_transitions, etc. but there is no harm for calling this explicitly.
- calculate_response(excitation_energy, excitation_direction, lorentzian_width, units='eVang')[source]¶
Calculates and returns response using TD-DFPT.
- Parameters
excitation_energy – Energy of the laser in given units
excitation_direction – Vector for direction (will be normalized)
lorentzian_width – Life time or width parameter. Larger width results in wider energy envelope around excitation energy.
- get_spectrum(filename=None, min_energy=0.0, max_energy=30.0, energy_step=0.01, width=0.1, units='eVcgs')[source]¶
Get spectrum for dipole and rotatory strength.
Returns folded spectrum as (w,S,R) where w is an array of frequencies, S is an array of corresponding dipole strengths, and R is an array of corresponding rotatory strengths.
- Parameters
min_energy – Minimum energy
min_energy – Maximum energy
energy_step – Spacing between calculated energies
width – Width of the Gaussian
units – Units for spectrum: ‘au’ or ‘eVcgs’
- get_transition_contributions(index_of_transition)[source]¶
Get contributions of Kohn-Sham singles to a given transition as number of electrons contributing.
Includes population difference.
This method is meant to be used for small number of transitions. It is not suitable for analysing densely packed transitions of large systems. Use transition contribution map (TCM) or similar approach for this.
- Parameters
index_of_transition – index of transition starting from zero
- get_transitions(filename=None, min_energy=0.0, max_energy=30.0, units='eVcgs')[source]¶
Get transitions: energy, dipole strength and rotatory strength.
Returns transitions as (w,S,R, Sx,Sy,Sz) where w is an array of frequencies, S is an array of corresponding dipole strengths, and R is an array of corresponding rotatory strengths.
- Parameters
min_energy – Minimum energy
min_energy – Maximum energy
units – Units for spectrum: ‘au’ or ‘eVcgs’
- class gpaw.lrtddft2.lr_communicators.LrCommunicators(world, dd_size: int, eh_size: Optional[int] = None)[source]¶
Create communicators for LrTDDFT calculation.
- Parameters
world – MPI parent communicator (usually
gpaw.mpi.world)dd_size – Number of domains for domain decomposition
eh_size – Number of groups for parallelization over electron-hole pairs
Note
Sizes must match, i.e., world.size must be equal to dd_size x eh_size, e.g., 1024 = 64*16
Tip
Use enough processes for domain decomposition (dd_size) to fit everything (easily) into memory, and use the remaining processes for electron-hole pairs as K-matrix build is trivially parallel over them.
Pass
lr_comms.dd_commto ground state calc when reading for LrTDDFT.Examples
For 8 MPI processes:
lr_comms = LrCommunicators(gpaw.mpi.world, 4, 2) txt = 'lr_%06d_%06d.txt' % (lr_comms.dd_comm.rank, lr_comms.eh_comm.rank) calc = GPAW('unocc.gpw', communicator=lr_comms.dd_comm) lr = LrTDDFT2(calc, lr_communicators=lr_comms, txt=txt)