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 *
from ase.io import *

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 ase import *
from ase.io import *
from ase.parallel import parprint

from gpaw import *
from gpaw.poisson import *
from gpaw.mpi import world, size, rank

from gpaw.lrtddft2 import *
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 import *
from ase.io import *
from ase.parallel import parprint

from gpaw import *
from gpaw.poisson import *
from gpaw.mpi import world, size, rank

from gpaw.lrtddft2 import *
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