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.set_calculator(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
world	Communicator	None	parent communicator (usually gpaw.mpi.world)
dd_size	\(int\)	None	Number of domains for domain decomposition
eh_size	\(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
basefilename	string		Prefix for all files created by LRTDDFT2 (e.g. dir/lr)
gs_calc	GPAW		Ground-state calculator, which has been loaded from a file with communicator=lr_communicators calculation
fxc	string	None	Exchange-correlation kernel
min_occ	int	None	Index of the first occupied state (inclusive)
max_occ	int	None	Index of the last occupied state (inclusive)
min_unocc	int	None	Index of the first unoccupied state (inclusive)
max_unocc	int	None	Index of the last unoccupied state (inclusive)
max_energy_diff	float	None	Maximum Kohn-Sham eigenenergy difference
recalculate	string	None	What should be recalculated. Usually nothing.
lr_communicators	LrCommuncators	None
txt	string	None	Output