Resonant and non-resonant Raman spectra¶

Note

Raman spectra can also be obtained via siesta calculator.

Raman spectra can be calculated in various approximations [1]. While the examples below are using GPAW explicitly, the modules are intended to work with other calculators also. The strategy is to calculate vibrational properties first and obtain the spectra from these later.

1. Finite difference calculations¶

1a. Forces and excitations¶

The basis for all spectra are finite difference calculations for forces and excited states of our system. These can be performed using the ResonantRamanCalculator.

class ase.vibrations.resonant_raman.ResonantRamanCalculator(atoms, ExcitationsCalculator, *args, exkwargs=None, exext='.ex.gz', overlap=False, **kwargs)[source]¶

Base class for resonant Raman calculators using finite differences.

Parameters:

atoms (Atoms) – The Atoms object
ExcitationsCalculator (object) – Calculator for excited states
exkwargs (dict) – Arguments given to the ExcitationsCalculator object
exext (string) – Extension for filenames of Excitation lists (results of the ExcitationsCalculator).
overlap (function or False) – Function to calculate overlaps between excitation at equilibrium and at a displaced position. Calculators are given as first and second argument, respectively.

Example

>>> from ase.calculators.h2morse import (H2Morse,
...                                      H2MorseExcitedStatesCalculator)
>>> from ase.vibrations.resonant_raman import ResonantRamanCalculator
>>>
>>> atoms = H2Morse()
>>> rmc = ResonantRamanCalculator(atoms, H2MorseExcitedStatesCalculator)
>>> rmc.run()

This produces all necessary data for further analysis.

1b. More accurate forces¶

It is possible to do a vibrational also with a more accurate calculator (or more accurate settings for the forces) using the Vibrations or Infrared modules.

In the example of molecular hydrogen with GPAW this is

from gpaw import GPAW, FermiDirac
from gpaw.cluster import Cluster

from ase import optimize
from ase.build import molecule
from ase.vibrations.infrared import Infrared

h = 0.22

atoms = Cluster(molecule('H2'))
atoms.minimal_box(3.5, h=h)

# relax the molecule
calc = GPAW(h=h, occupations=FermiDirac(width=0.1))
atoms.calc = calc
dyn = optimize.FIRE(atoms)
dyn.run(fmax=0.05)
atoms.write('relaxed.traj')

# finite displacement for vibrations
atoms.calc.set(symmetry={'point_group': False})
ir = Infrared(atoms)
ir.run()

This produces a calculation with rather accurate forces in order to get the Hessian and thus the vibrational frequencies as well as Eigenstates correctly.

In the next step we perform a finite difference optical calculation with less accuracy, where the optical spectra are evaluated using TDDFT

from gpaw import GPAW, FermiDirac
from gpaw.cluster import Cluster
from gpaw.lrtddft import LrTDDFT

from ase.vibrations.resonant_raman import ResonantRamanCalculator

h = 0.25
atoms = Cluster('relaxed.traj')
atoms.minimal_box(3.5, h=h)

# relax the molecule
calc = GPAW(h=h, occupations=FermiDirac(width=0.1),
            eigensolver='cg', symmetry={'point_group': False},
            nbands=10, convergence={'eigenstates': 1.e-5,
                                    'bands': 4})
atoms.calc = calc

# use only the 4 converged states for linear response calculation
rmc = ResonantRamanCalculator(atoms, LrTDDFT,
                              exkwargs={'restrict': {'jend': 3}})
rmc.run()

1c. Overlaps¶

Albrecht B+C terms need wave function overlaps between equilibrium and displaced structures. These are assumed to be calculated in the form

\[o_{ij} = \int d\vec{r} \; \phi_i^{{\rm disp},*}(\vec{r}) \phi_j^{{\rm eq}}(\vec{r})\]

where \(\phi_j^{{\rm eq}}\) is an orbital at equilibrium position and \(\phi_i^{\rm disp}\) is an orbital at displaced position.

The H2MorseExcitedStatesCalculator has a function overlap() for this. We therfore write data including the overlap as

from ase.calculators.h2morse import H2Morse, H2MorseExcitedStatesCalculator
from ase.vibrations.resonant_raman import ResonantRamanCalculator

atoms = H2Morse()
rmc = ResonantRamanCalculator(atoms, H2MorseExcitedStatesCalculator,
                              overlap=lambda x, y: x.overlap(y))
rmc.run()

In GPAW this is implemented in Overlap (approximated by pseudo-wavefunction overlaps) and can be triggered in ResonantRamanCalculator by

from gpaw import GPAW, FermiDirac
from gpaw.analyse.overlap import Overlap
from gpaw.cluster import Cluster
from gpaw.lrtddft import LrTDDFT

from ase.vibrations.resonant_raman import ResonantRamanCalculator

h = 0.25
atoms = Cluster('relaxed.traj')
atoms.minimal_box(3.5, h=h)

# relax the molecule
calc = GPAW(h=h, occupations=FermiDirac(width=0.1),
            symmetry={'point_group': False},
            nbands=10, convergence={'eigenstates': 1.e-5,
                                    'bands': 4})
atoms.calc = calc

# use only the 4 converged states for linear response calculation
rmc = ResonantRamanCalculator(atoms, LrTDDFT, name='rroverlap',
                              exkwargs={'restrict': {'jend': 3}},
                              overlap=lambda x, y: Overlap(x).pseudo(y)[0])
rmc.run()

2. Analysis of the results¶

We assume that the steps above were performed and are able to analyse the results in different approximations.

Placzek¶

The most popular form is the Placzeck approximation that is present in two implementations. The simplest is the direct evaluation from derivatives of the frequency dependent polarizability:

from ase.calculators.h2morse import (H2Morse,
                                     H2MorseExcitedStates)
from ase.vibrations.placzek import Placzek

photonenergy = 7.5  # eV
pz = Placzek(H2Morse(), H2MorseExcitedStates)
pz.summary(photonenergy)

The second implementation evaluates the derivatives differently, allowing for more analysis:

import pylab as plt
from ase.calculators.h2morse import (H2Morse,
                                     H2MorseExcitedStates)
from ase.vibrations.placzek import Profeta

photonenergy = 7.5  # eV
pr = Profeta(H2Morse(), H2MorseExcitedStates, approximation='Placzek')
x, y = pr.get_spectrum(photonenergy, start=4000, end=5000, type='Lorentzian')
plt.plot(x, y)
plt.show()

Both implementations should lead to the same spectrum.

Profeta splits the spectra in two contributions that can be accessed as approximation='P-P' and approximation='Profeta', respectively. Their sum should give approximation='Placzek'. See more details in [1].

Albrecht¶

The more accurate Albrecht approximations partly need overlaps to be present. We therefore have to invoke the Albrecht object as:

from ase.calculators.h2morse import (H2Morse,
                                     H2MorseExcitedStates)
from ase.vibrations.albrecht import Albrecht

photonenergy = 7.5  # eV
al = Albrecht(H2Morse(), H2MorseExcitedStates, approximation='Albrecht', overlap=True)
x, y = al.get_spectrum(photonenergy, start=4000, end=5000, type='Lorentzian')

Albrecht splits the spectra in two contributions that can be accessed as approximation='Albrecht A' and approximation='Albrecht BC', respectively. Their sum should give approximation='Albrecht'. See more details in [1].