In [ ]:
# import libraries and set up the molecule geometry
from ase.units import Ry, eV, Ha
from ase.calculators.siesta import Siesta
from ase import Atoms
import numpy as np
import matplotlib.pyplot as plt
from ase.build import molecule
CH4 = molecule("CH4")
# visualization of the particle
from ase.visualize import view
view(CH4, viewer='x3d')
In [ ]:
# enter siesta input and run siesta
siesta = Siesta(
mesh_cutoff=150 * Ry,
basis_set='DZP',
pseudo_qualifier='lda',
energy_shift=(10 * 10**-3) * eV,
fdf_arguments={
'SCFMustConverge': False,
'COOP.Write': True,
'WriteDenchar': True,
'PAO.BasisType': 'split',
'DM.Tolerance': 1e-4,
'DM.MixingWeight': 0.1,
'MaxSCFIterations': 300,
'DM.NumberPulay': 4,
'XML.Write': True})
CH4.set_calculator(siesta)
e = CH4.get_potential_energy()
In [ ]:
# compute polarizability using pyscf-nao
siesta.pyscf_tddft_eels(label="siesta", jcutoff=7, iter_broadening=0.15/Ha,
xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7, freq = np.arange(0.0, 15.0, 0.05),
velec=np.array([50.0, 0.0, 0.0]), b = np.array([0.0, 0.0, 5.0]))
In [ ]:
# plot polarizability with matplotlib
%matplotlib inline
fig = plt.figure(1)
ax1 = fig.add_subplot(121)
ax2 = fig.add_subplot(122)
ax1.plot(siesta.results["freq range"], siesta.results["eel spectra nonin"].imag)
ax2.plot(siesta.results["freq range"], siesta.results["eel spectra inter"].imag)
ax1.set_xlabel(r"$\omega$ (eV)")
ax2.set_xlabel(r"$\omega$ (eV)")
ax1.set_ylabel(r"$\Gamma_{nonin}$ (au)")
ax2.set_ylabel(r"$\Gamma_{inter}$ (au)")
ax1.set_title(r"Non interacting")
ax2.set_title(r"Interacting")
fig.tight_layout()
In [ ]: