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