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. It is safer to use ‘dav’ or ‘cg’ eigensolver instead of the default ‘rmm-diis’ eigensolver to converge unoccupied states. However, ‘dav’ and ‘cg’ are often too expensive for large systems. In this case, you should use ‘rmm-diis’ with tens or hundreds of extra states in addition to the unoccupied states you wish to converge.
10
http://cccbdb.nist.gov/ Geometry for C3H6O (Propylene oxide), CISD/6-31G*
O 0.8171 -0.7825 -0.2447
C -1.5018 0.1019 -0.1473
H -1.3989 0.3323 -1.2066
H -2.0652 -0.8262 -0.0524
H -2.0715 0.8983 0.3329
C -0.1488 -0.0393 0.4879
H -0.1505 -0.2633 1.5506
C 1.0387 0.6105 -0.0580
H 0.9518 1.2157 -0.9531
H 1.8684 0.8649 0.5908
from gpaw import GPAW
from ase.io import read
L = 14 # ang
atoms = read('r-methyl-oxirane.xyz')
atoms.set_cell((L, L, L))
atoms.center()
calc = GPAW(h=0.25,
nbands=30,
convergence={'bands': 20},
xc='LDA',
maxiter=300,
txt='unocc.txt')
atoms.calc = calc
atoms.get_potential_energy()
calc.write('r-methyl-oxirane.gpw', mode='all')
Calculating response matrix and spectrum¶
The next step is to calculate the response matrix:
from gpaw import GPAW
from gpaw.poisson import PoissonSolver
from gpaw.mpi import world, size
from gpaw.lrtddft2 import LrTDDFT2
from gpaw.lrtddft2.lr_communicators import LrCommunicators
# atoms, gs_calc = restart('r-methyl-oxirane.gpw')
dd_size = 2 * 2 * 2
eh_size = size // dd_size
assert eh_size * dd_size == size
lr_comms = LrCommunicators(world, dd_size, eh_size)
gs_calc = GPAW('r-methyl-oxirane.gpw',
communicator=lr_comms.dd_comm,
txt='lr_gs.txt',
poissonsolver=PoissonSolver(eps=1e-20,
remove_moment=(1 + 3 + 5)))
lr = LrTDDFT2('lrtddft2',
gs_calc,
fxc='LDA',
min_occ=0, # usually zero
max_occ=15, # a few above LUMO
min_unocc=10, # a few before HOMO
max_unocc=20, # number of converged states
max_energy_diff=7.0, # eV
recalculate=None,
lr_communicators=lr_comms,
txt='-')
# This is the expensive part!
lr.calculate()
# transitions
(w, S, R, Sx, Sy, Sz) = lr.get_transitions('trans.dat')
# spectrum
lr.get_spectrum('spectrum.dat', 0, 10.)
Note: Unfortunately, spin is not implemented yet. For now, use ‘lrtddft’.
Restarting, recalculating and analyzing spectrum¶
from ase.parallel import parprint
from gpaw import GPAW
from gpaw.poisson import PoissonSolver
from gpaw.mpi import world, size
from gpaw.lrtddft2 import LrTDDFT2
from gpaw.lrtddft2.lr_communicators import LrCommunicators
# atoms, gs_calc = restart('r-methyl-oxirane.gpw')
dd_size = 2 * 2 * 2
eh_size = size // dd_size
assert eh_size * dd_size == size
lr_comms = LrCommunicators(world, dd_size, eh_size)
gs_calc = GPAW('r-methyl-oxirane.gpw',
communicator=lr_comms.dd_comm,
txt='lr_gs.txt',
poissonsolver=PoissonSolver(eps=1e-20,
remove_moment=(1 + 3 + 5)))
lr = LrTDDFT2('lrtddft2',
gs_calc,
fxc='LDA',
min_occ=0, # usually zero
max_occ=15, # a few above LUMO
min_unocc=10, # a few before HOMO
max_unocc=20, # number of converged states
max_energy_diff=8.0, # eV
recalculate=None,
lr_communicators=lr_comms,
txt='-')
# This is the expensive part!
lr.calculate()
# transitions
(w, S, R, Sx, Sy, Sz) = lr.get_transitions('trans.dat')
# spectrum
lr.get_spectrum('spectrum.dat', 0, 10.)
parprint('')
parprint(' %5s => %5s contribution' % ('occ', 'unocc'))
f2 = lr.lr_transitions.get_transition_contributions(0)
for (ip, val) in enumerate(f2):
if (val > 1e-3):
parprint(' %5d => %5d %lf %%\n ' %
(lr.ks_singles.kss_list[ip].occ_ind,
lr.ks_singles.kss_list[ip].unocc_ind,
val / 2. * 100))
Quick reference¶
Parameters for LrCommunicators:
keyword |
type |
default value |
description |
---|---|---|---|
|
|
None |
parent communicator (usually gpaw.mpi.world) |
|
\(int\) |
None |
Number of domains for domain decomposition |
|
\(int\) |
None |
Number of groups for parallelization over e-h -pairs |
Note: world.size = dd_size x eh_size
Parameters for LrTDDFT2:
keyword |
type |
default value |
description |
---|---|---|---|
|
|
Prefix for all files created by LRTDDFT2
(e.g. |
|
|
|
Ground-state calculator, which has been loaded from a file with communicator=lr_communicators calculation |
|
|
|
None |
Exchange-correlation kernel |
|
|
None |
Index of the first occupied state (inclusive) |
|
|
None |
Index of the last occupied state (inclusive) |
|
|
None |
Index of the first unoccupied state (inclusive) |
|
|
None |
Index of the last unoccupied state (inclusive) |
|
|
None |
Maximum Kohn-Sham eigenenergy difference |
|
|
None |
What should be recalculated. Usually nothing. |
|
|
None |
|
|
|
None |
Output |