.. _tutorial2-tddftCH4-siesta:
Tutorial 2: calculate the polarizability of CH4 molecule with Siesta
====================================================================
In this second tutorial, we will learn how to calculate the
polarizability of a molecule using TDDFT as in tutorial 1. However in
this case, we will use
`Siesta `__ for the DFT
calculations.
**Requirements:** have PyNAO, PySCF and Siesta installed
Run DFT with Siesta
-------------------
We assume here that Siesta is properly installed and configured on your
system. We will assume also that ``siesta`` executable as been added to
your ``PATH``
Method 1: run Siesta via the command line
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The first way to run Siesta is to call the executable from the command
line. For this, you should make sure to
- to have a siesta ``fdf`` input file on the disk
- to have the necessary pseudopotential in the folder
First, we will define the ``siesta.fdf`` file (we already included the
pseudopotential files in the folder). In the ``siesta.fdf`` file
- take note of the ``SystemLabel`` variable, you will need it for PyNAO
- the option ``COOP.Write`` must be set to ``True``
- the option ``XML.Write`` must be set to ``True``
Let write the ``siesta.fdf`` file to disk
.. code:: ipython3
siestafdf = """SystemName siesta
SystemLabel siesta
COOP.Write True
DM.MixingWeight 0.01
DM.NumberPulay 4
DM.Tolerance 0.0001
MD.MaxForceTol 0.02 eV/Ang
MD.NumCGsteps 0
MaxSCFIterations 10000
PAO.BasisType split
SCFMustConverge False
WriteCoorXmol True
WriteDenchar True
XML.Write True
Spin non-polarized
XC.functional GGA
XC.authors PBE
MeshCutoff 3401.4232530459053 eV
PAO.EnergyShift 0.025 eV
NumberOfSpecies 2
NumberOfAtoms 5
%block ChemicalSpecieslabel
1 6 C
2 1 H
%endblock ChemicalSpecieslabel
%block PAO.BasisSizes
C DZP
H DZP
%endblock PAO.BasisSizes
AtomicCoordinatesFormat Ang
%block AtomicCoordinatesAndAtomicSpecies
0.000000000 0.000000000 0.000000000 1
0.629118000 0.629118000 0.629118000 2
-0.629118000 -0.629118000 0.629118000 2
0.629118000 -0.629118000 -0.629118000 2
-0.629118000 0.629118000 -0.629118000 2
%endblock AtomicCoordinatesAndAtomicSpecies
"""
with open("siesta.fdf", "w") as fl:
fl.write(siestafdf)
We can now launch DFT calculation with Siesta via the command line
.. code:: ipython3
!siesta siesta.out
Method 2: run Siesta via ASE
~~~~~~~~~~~~~~~~~~~~~~~~~~~~
The `Atomic Simulation Environment
(ASE) `__ provides all an environment to
run atomistic calculations. It is possible to run Siesta calculations
via ASE, which is more friendly than running Siesta from the command
line.
The advantage of using ASE is that you do not need to write the siesta
``fdf`` input file by hand, and you don’t need to copy the
pseudopotentials. However, make sure that the `ASE Siesta
calculator `__
is properly setup.
.. code:: ipython3
from ase.calculators.siesta import Siesta
from ase.units import eV, Ry
# First we need to build the molecule
# It can be done directly via ASE
from ase.build import molecule
atoms = molecule('CH4')
# Or you can load the molecule geometry from file
#import ase.io as io
#atoms = io.read(fname)
# set-up the Siesta parameters
siesta = Siesta(
mesh_cutoff=250*Ry,
basis_set='DZP',
pseudo_qualifier='gga',
xc="PBE",
energy_shift=(25 * 10**-3) * eV,
fdf_arguments={
'SCFMustConverge': False,
'COOP.Write': True,
'WriteDenchar': True,
'PAO.BasisType': 'split',
'DM.Tolerance': 1e-4,
'DM.MixingWeight': 0.01,
"MD.NumCGsteps": 0,
"MD.MaxForceTol": (0.02, "eV/Ang"),
'MaxSCFIterations': 10000,
'DM.NumberPulay': 4,
'XML.Write': True,
"WriteCoorXmol": True})
atoms.set_calculator(siesta)
# Run Siesta calculation
e = atoms.get_potential_energy()
print("DFT potential energy", e)
Run TDDFT with PyNAO
--------------------
Now we can run TDDFT calculations similarly to what was done in tutorial
1 with PySCF.
The first step is to initialize the system and calculate the kernel by
calling ``tddft_iter``.
For Siesta inputs, we need to call ``tddft_iter`` with the ``label``
parameter. It must be the same label than the one indicated by
``SystemLabel`` in the Siesta ``fdf`` file.
.. code:: ipython3
from pynao import tddft_iter
td = tddft_iter(label="siesta")
Once the kernel has been calculated, we can calculate the polarizability
of the molecule by calling the method ``comp_polariz_inter_Edir``.
This method takes two main inputs:
- ``freq``: the frequencies at which to calculate the polarizability.
It must be a complex number, the real part is the actual frequencies
points, while the imaginary part is the actual broadening of the
polarizability. In general it will be ``tddft_iter.eps``. The
frequencies must be given in Hartree.
- ``Edir``: the direction of the external electric field
.. code:: ipython3
import numpy as np
from ase.units import Ha
freq = np.arange(0.0, 25.0, 0.05)/Ha + 1j * td.eps
p_mat = td.comp_polariz_inter_Edir(freq, Eext=np.array([1.0, 1.0, 1.0]))
.. parsed-literal::
Total number of iterations: 10853
.. code:: ipython3
import matplotlib.pyplot as plt
h = 7
w = 4*h/3
ft = 20
fig = plt.figure(1, figsize=(w, h))
ax = fig.add_subplot(111)
ax.plot(freq.real*Ha, -p_mat[0,0, :].imag, linewidth=3)
ax.set_xlabel(r"Energy (eV)", fontsize=ft)
ax.set_ylabel(r"$P_{xx}$ (a.u.)", fontsize=ft)
fig.tight_layout()
.. image:: output_10_0.png
Remark
------
If you compare the spectra obtained in tutorial 1 and tutorial 2, you
will see a quite different spectra, while the system is the same. Since
the TDDFT part remained inchanged, it demonstrates how crucial is the
quality of the DFT inputs in order to obtain reliable results.
The purpose of these tutorial is to demonstrate how to use PyNAO with
PySCF and Siesta. It is the responsability of the user to ensure the
quality of the calculations he is running.