Skip to article frontmatterSkip to article content
Site not loading correctly?

This may be due to an incorrect BASE_URL configuration. See the MyST Documentation for reference.

Optical activity and dichroism

VeloxChem enables simulations of molecular optical activity by providing access to rotatory strengths, specific rotations, circular dichroism (ECD) spectra, and optical rotatory dispersion (ORD) through density‑functional theory response methods. These chiroptical properties can be obtained either from an eigenvalue‑based TDDFT treatment of excited states or through the complex polarization propagator (CPP) approach, which delivers frequency‑dependent spectra directly comparable to experiment. For systems composed of weakly interacting chromophores, VeloxChem also implements an exciton coupling model to compute ECD spectra from the interactions between local excitations. For further theoretical details, see Norman et al. (2018).

In this section, (S)-methyloxirane is used to showcase the calculations of rotatory strengths via TDDFT, extinction coefficients via CPP, and specific rotations and ORD via CPP. Finally, the exciton coupling model (ECM) is illustrated using a helical stack of ethylene molecules.

Loading...

Rotatory strengths

The strength of an ECD band is given by the anisotropy of the decadic molar extinction coefficient

Δϵ(ω)=16πNAln⁡(10)(4πε0)c2π3ℏ∑n>0f(ω;ωn0,γ) ωn0Rn0\Delta\epsilon(\omega) = \frac{ 16\pi N_\mathrm{A} }{ \ln\left(10\right) \left(4\pi\varepsilon_0\right) c^2 } \frac{\pi}{3 \hbar} \sum_{n>0} f(\omega; \omega_{n0},\gamma)\, \omega_{n0} R_{n0}
(1)

where NAN_\mathrm{A} is Avogadro’s constant, ff is the Cauchy distribution, and Rn0R_{n0} is the rotatory strength defined as

Rn0=∑α=x,y,zℑ⟨0∣μ^α∣n⟩⟨n∣m^α∣0⟩=∑α=x,y,z−emeωn0⟨0∣p^α∣n⟩⟨n∣m^α∣0⟩R_{n0} = \sum_{\alpha = x,y,z} \Im \langle 0 | \hat{\mu}_\alpha | n \rangle \langle n | \hat{m}_\alpha | 0\rangle = \sum_{\alpha = x,y,z} \frac{-e}{m_\mathrm{e} \omega_{n0}} \langle 0 | \hat{p}_\alpha | n \rangle \langle n | \hat{m}_\alpha | 0\rangle
(2)

In VeloxChem, the rotatory strength is evaluated in the velocity gauge as given in the second expression, and it is thereby gauge-origin independent.

Python script

import veloxchem as vlx

molecule = vlx.Molecule.read_xyz_string("""10

O             -1.009866880244         1.299407071912         1.951947754409
C              0.031038818173         2.001294498224         1.244371321259
C              0.339506903623         0.795807201271         2.029917688849
C              0.600316751832        -0.544462572436         1.394186918859
H             -0.026285836651         1.930258644799         0.154459855148
H              0.271804379347         2.990968281608         1.641505800108
H              0.794935917667         0.948100885444         3.014899839267
H              0.113110324011        -0.610670447580         0.412641743927
H              1.681576973773        -0.701096452623         1.264343536564
H              0.213803265264        -1.354101659046         2.029564064183""")


basis = vlx.MolecularBasis.read(molecule, "aug-cc-pvdz")

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.xcfun = "b3lyp"
scf_results = scf_drv.compute(molecule, basis)

rsp_drv = vlx.LinearResponseEigenSolver()
rsp_drv.nstates = 8
rsp_results = rsp_drv.compute(molecule, basis, scf_results)
<Figure size 800x500 with 2 Axes>

Text file

@jobs
task: response
@end

@method settings
basis: def2-svpd
xcfun: b3lyp
@end

@response
property: ecd
nstates: 10
nto: yes
@end

@molecule
charge: 0
multiplicity: 1
xyz:
...
@end

Extinction coefficient

The anisotropy of the decadic molar extinction coefficient can be determined directly from the complex polarization propagator evaluated for mixed electric- and magnetic-dipole operators Jiemchooroj & Norman (2007)

Δϵ(ω)=16πNAω2ln⁡(10)(4πε0)c2 β(ω)\Delta\epsilon(\omega) = \frac{ 16 \pi N_\mathrm{A} \omega^2 }{ \ln(10) \left(4\pi\varepsilon_0\right) c^2 } \, \beta(\omega)
(3)

where the molecular response property, β(ω)\beta(\omega), is defined as

β(ω)=−13ω(Gxx+Gyy+Gzz)\beta(\omega) = -\frac{1}{3 \omega} (G_{xx} + G_{yy} + G_{zz})
(4)

and

Gαβ=−ℜ⟨⟨μ^α;m^β⟩⟩ωγ=−eωmeℑ⟨⟨p^α;m^β⟩⟩ωγG_{\alpha\beta} = - \Re\langle\langle\hat{\mu}_\alpha;\hat{m}_\beta \rangle\rangle_\omega^\gamma = - \frac{e}{\omega m_e} \Im \langle\langle\hat{p}_\alpha; \hat{m}_\beta \rangle\rangle_\omega^\gamma
(5)

The mixed electric–magnetic dipole tensor, GG, is evaluated in the velocity gauge as given in the second expression. Furthermore, it is complex and calculated with a damping term, ℏγ\hbar \gamma, associated with the inverse finite lifetime of the excited states. The default program setting for this parameter is 0.124 eV (or 0.004556 a.u.).

The resulting values for Δϵ(ω)\Delta \epsilon(\omega) are converted from atomic units to units of L mol−1^{-1} cm−1^{-1} by multiplying with a factor of 10 a0210\, a_0^2.

Python script

import numpy as np
import veloxchem as vlx

molecule = vlx.Molecule.read_xyz_string("""10

O             -1.009866880244         1.299407071912         1.951947754409
C              0.031038818173         2.001294498224         1.244371321259
C              0.339506903623         0.795807201271         2.029917688849
C              0.600316751832        -0.544462572436         1.394186918859
H             -0.026285836651         1.930258644799         0.154459855148
H              0.271804379347         2.990968281608         1.641505800108
H              0.794935917667         0.948100885444         3.014899839267
H              0.113110324011        -0.610670447580         0.412641743927
H              1.681576973773        -0.701096452623         1.264343536564
H              0.213803265264        -1.354101659046         2.029564064183""")

basis = vlx.MolecularBasis.read(molecule, "aug-cc-pvdz")

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.xcfun = "b3lyp"
scf_results = scf_drv.compute(molecule, basis)

cpp_drv = vlx.ComplexResponse()
cpp_drv.frequencies = np.arange(0.207, 0.325, 0.0025)
cpp_drv.property = "ecd"

cpp_results = cpp_drv.compute(molecule, basis, scf_results)
<Figure size 800x500 with 1 Axes>

Text file

@jobs
task: response
@end

@method settings
basis: def2-svpd
xcfun: b3lyp
@end

@response
property: ecd (cpp)
# frequency region (and resolution)
frequencies:  0.207-0.325 (0.0025)
@end

@molecule
charge: 0
multiplicity: 1
xyz:
...
@end

Optical Rotatory Dispersion

The specific optical rotation can be determined directly from the complex polarization propagator evaluated for mixed electric- and magnetic-dipole operators Norman et al. (2004)

[α]ω=−18000 NAω2π(4πε0)c2M β′(ω)[\alpha]_\omega = \frac{ -18000\, N_\mathrm{A} \omega^2 }{ \pi \left(4\pi\varepsilon_0\right) c^2 M } \, \beta'(\omega)
(6)

where MM is the molar mass and the molecular response property, β′(ω)\beta'(\omega), is defined as

β′(ω)=13ω(Gxx′+Gyy′+Gzz′)\beta'(\omega) = \frac{1}{3 \omega} (G'_{xx} + G'_{yy} + G'_{zz})
(7)

and

Gαβ′=−ℑ⟨⟨μ^α;m^β⟩⟩ωγ=eωmeℜ⟨⟨p^α;m^β⟩⟩ωγG'_{\alpha\beta} = - \Im\langle\langle\hat{\mu}_\alpha;\hat{m}_\beta \rangle\rangle_\omega^\gamma = \frac{e}{\omega m_e} \Re \langle\langle\hat{p}_\alpha; \hat{m}_\beta \rangle\rangle_\omega^\gamma
(8)

The mixed electric–magnetic dipole tensor, G′G', is evaluated in the velocity gauge as given in the second expression, and is thereby gauge-origin independent. It is complex and calculated with a damping term, ℏγ\hbar \gamma, associated with the inverse finite lifetime of the excited states. The default program setting for this parameter is 0.124 eV (or 0.004556 a.u.).

The resulting values for [α]ω[\alpha]_\omega are reported in units of deg dm−1^{-1} (g cm−3^{-3})−1^{-1}.

Python script

import numpy as np
import veloxchem as vlx

molecule = vlx.Molecule.read_xyz_string("""10

O             -1.009866880244         1.299407071912         1.951947754409
C              0.031038818173         2.001294498224         1.244371321259
C              0.339506903623         0.795807201271         2.029917688849
C              0.600316751832        -0.544462572436         1.394186918859
H             -0.026285836651         1.930258644799         0.154459855148
H              0.271804379347         2.990968281608         1.641505800108
H              0.794935917667         0.948100885444         3.014899839267
H              0.113110324011        -0.610670447580         0.412641743927
H              1.681576973773        -0.701096452623         1.264343536564
H              0.213803265264        -1.354101659046         2.029564064183""")

basis = vlx.MolecularBasis.read(molecule, "aug-cc-pvdz")

scf_drv = vlx.ScfRestrictedDriver()
scf_drv.xcfun = "b3lyp"
scf_results = scf_drv.compute(molecule, basis)

cpp_drv = vlx.ComplexResponse()
cpp_drv.frequencies = np.arange(0.207, 0.325, 0.0025)
cpp_drv.property = "ord"

cpp_results = cpp_drv.compute(molecule, basis, scf_results)
<Figure size 800x500 with 1 Axes>

Text file

@jobs
task: response
@end

@method settings
basis: def2-svpd
xcfun: b3lyp
@end

@response
property: ord (cpp)
# frequency region (and resolution)
frequencies: 0.207-0.325 (0.0025)
@end

@molecule
charge: 0
multiplicity: 1
xyz:
...
@end

Exciton coupling model

VeloxChem implements the exciton coupling model (ECM) by Li et al. (2017) to determine circular dichroism spectra.

The ECM describes optical activity in molecular aggregates by representing the excited states of the full system in a reduced basis of localized monomer excitations and intermolecular charge-transfer (CT) configurations. In this framework, the electronic Hamiltonian is projected onto a subspace spanned by a user-defined number of excited states on each monomer, augmented by CT states constructed from selected occupied and virtual orbitals on neighboring units. The resulting excitonic Hamiltonian captures both Coulombic coupling between local excitations and explicit electron-hole transfer between monomers. Diagonalization of this Hamiltonian yields delocalized excited states from which rotatory strengths and circular dichroism spectra can be evaluated efficiently while retaining a clear connection to the underlying monomer electronic structure.

Loading...

Figure: A helical stack of 10 ethylene molecules (or fragments).

Python script

import veloxchem as vlx

# the molecule object can be built from a molecular string listing all fragments 
# in the helical stack and all atoms in each fragment
molecule = vlx.Molecule.read_str(mol_str)
basis = vlx.MolecularBasis.read(molecule, "def2-svp")

exciton_drv = vlx.ExcitonModelDriver()

exciton_drv.update_settings({"fragments": "10", "atoms_per_fragment": "6"})

exciton_drv.xcfun_label = "b3lyp"
exciton_drv.nstates = 1
exciton_drv.ct_nocc = 0
exciton_drv.ct_nvir = 0

exciton_results = exciton_drv.compute(molecule, basis)
<Figure size 700x400 with 1 Axes>

Text file

@jobs
task: exciton
@end

@method settings
xcfun: b3lyp
basis: cc-pvdz
@end

@exciton
fragments: 40
atoms_per_fragment: 55
nstates: 5
ct_nocc: 0
ct_nvir: 0
@end

@molecule
charge: 0
multiplicity: 1
xyz:
...
# XYZ coordinates for 40 x 55 atoms
...
@end
References
  1. Norman, P., Ruud, K., & Saue, T. (2018). Principles and practices of molecular properties. John Wiley & Sons, Ltd.
  2. Jiemchooroj, A., & Norman, P. (2007). Electronic circular dichroism spectra from the complex polarization propagator. J. Chem. Phys., 126(13), 134102. 10.1063/1.2716660
  3. Norman, P., Ruud, K., & Helgaker, T. (2004). Density-functional theory calculations of optical rotatory dispersion in the nonresonant and resonant frequency regions. J. Chem. Phys., 120(11), 5027–5035. 10.1063/1.1647515
  4. Li, X., Parrish, R. M., Liu, F., Kokkila Schumacher, S. I. L., & Martinez, T. J. (2017). An Ab Initio Exciton Model Including Charge-Transfer Excited States. J. Chem. Theory Comput., 13(8), 3493–3504. 10.1021/acs.jctc.7b00171
Properties
Polarizability
Properties
Vibrational spectroscopies