Coverage for /builds/ase/ase/ase/calculators/siesta/siesta_lrtddft.py: 19.64%
56 statements
« prev ^ index » next coverage.py v7.5.3, created at 2025-08-02 00:12 +0000
1# fmt: off
3import numpy as np
5import ase.units as un
6from ase.calculators.polarizability import StaticPolarizabilityCalculator
9class SiestaLRTDDFT:
10 """Interface for linear response TDDFT for Siesta via `PyNAO`_
12 When using PyNAO please cite the papers indicated in the `documentation \
13<https://mbarbrywebsite.ddns.net/pynao/doc/html/references.html>`_
14 """
16 def __init__(self, initialize=False, **kw):
17 """
18 Parameters
19 ----------
20 initialize: bool
21 To initialize the tddft calculations before
22 calculating the polarizability
23 Can be useful to calculate multiple frequency range
24 without the need to recalculate the kernel
25 kw: dictionary
26 keywords for the tddft_iter function from PyNAO
27 """
29 try:
30 from pynao import tddft_iter
31 except ModuleNotFoundError as err:
32 msg = ("running lrtddft with Siesta calculator "
33 "requires pynao package")
34 raise ModuleNotFoundError(msg) from err
36 self.initialize = initialize
37 self.lrtddft_params = kw
38 self.tddft = None
40 # convert iter_broadening to Ha
41 if "iter_broadening" in self.lrtddft_params:
42 self.lrtddft_params["iter_broadening"] /= un.Ha
44 if self.initialize:
45 self.tddft = tddft_iter(**self.lrtddft_params)
47 def get_ground_state(self, atoms, **kw):
48 """
49 Run siesta calculations in order to get ground state properties.
50 Makes sure that the proper parameters are unsed in order to be able
51 to run pynao afterward, i.e.,
53 COOP.Write = True
54 WriteDenchar = True
55 XML.Write = True
56 """
57 from ase.calculators.siesta import Siesta
59 if "fdf_arguments" not in kw.keys():
60 kw["fdf_arguments"] = {"COOP.Write": True,
61 "WriteDenchar": True,
62 "XML.Write": True}
63 else:
64 for param in ["COOP.Write", "WriteDenchar", "XML.Write"]:
65 kw["fdf_arguments"][param] = True
67 siesta = Siesta(**kw)
68 atoms.calc = siesta
69 atoms.get_potential_energy()
71 def get_polarizability(self, omega, Eext=np.array(
72 [1.0, 1.0, 1.0]), inter=True):
73 """
74 Calculate the polarizability of a molecule via linear response TDDFT
75 calculation.
77 Parameters
78 ----------
79 omega: float or array like
80 frequency range for which the polarizability should be
81 computed, in eV
83 Returns
84 -------
85 polarizability: array like (complex)
86 array of dimension (3, 3, nff) with nff the number of frequency,
87 the first and second dimension are the matrix elements of the
88 polarizability in atomic units::
90 P_xx, P_xy, P_xz, Pyx, .......
92 Example
93 -------
95 from ase.calculators.siesta.siesta_lrtddft import siestaLRTDDFT
96 from ase.build import molecule
97 import numpy as np
98 import matplotlib.pyplot as plt
100 # Define the systems
101 CH4 = molecule('CH4')
103 lr = siestaLRTDDFT(label="siesta", jcutoff=7, iter_broadening=0.15,
104 xc_code='LDA,PZ', tol_loc=1e-6, tol_biloc=1e-7)
106 # run DFT calculation with Siesta
107 lr.get_ground_state(CH4)
109 # run TDDFT calculation with PyNAO
110 freq=np.arange(0.0, 25.0, 0.05)
111 pmat = lr.get_polarizability(freq)
112 """
113 from pynao import tddft_iter
115 if not self.initialize:
116 self.tddft = tddft_iter(**self.lrtddft_params)
118 if isinstance(omega, float):
119 freq = np.array([omega])
120 elif isinstance(omega, list):
121 freq = np.array([omega])
122 elif isinstance(omega, np.ndarray):
123 freq = omega
124 else:
125 raise ValueError("omega soulf")
127 freq_cmplx = freq / un.Ha + 1j * self.tddft.eps
128 if inter:
129 pmat = -self.tddft.comp_polariz_inter_Edir(freq_cmplx, Eext=Eext)
130 self.dn = self.tddft.dn
131 else:
132 pmat = -self.tddft.comp_polariz_nonin_Edir(freq_cmplx, Eext=Eext)
133 self.dn = self.tddft.dn0
135 return pmat
138class RamanCalculatorInterface(SiestaLRTDDFT, StaticPolarizabilityCalculator):
139 """Raman interface for Siesta calculator.
140 When using the Raman calculator, please cite
142 M. Walter and M. Moseler, Ab Initio Wavelength-Dependent Raman
143 Spectra: Placzek Approximation and Beyond, J. Chem. Theory
144 Comput. 2020, 16, 1, 576–586"""
146 def __init__(self, omega=0.0, **kw):
147 """
148 Parameters
149 ----------
150 omega: float
151 frequency at which the Raman intensity should be computed, in eV
153 kw: dictionary
154 The parameter for the siesta_lrtddft object
155 """
157 self.omega = omega
158 super().__init__(**kw)
160 def calculate(self, atoms):
161 """
162 Calculate the polarizability for frequency omega
164 Parameters
165 ----------
166 atoms: atoms class
167 The atoms definition of the system. Not used but required by Raman
168 calculator
169 """
170 pmat = self.get_polarizability(
171 self.omega, Eext=np.array([1.0, 1.0, 1.0]))
173 # Specific for raman calls, it expects just the tensor for a single
174 # frequency and need only the real part
176 # For static raman, imaginary part is zero??
177 # Answer from Michael Walter: Yes, in the case of finite systems you may
178 # choose the wavefunctions to be real valued. Then also the density
179 # response function and hence the polarizability are real.
181 # Convert from atomic units to e**2 Ang**2/eV
182 return pmat[:, :, 0].real * (un.Bohr**2) / un.Ha
185def pol2cross_sec(p, omg):
186 """
187 Convert the polarizability in au to cross section in nm**2
189 Input parameters:
190 -----------------
191 p (np array): polarizability from mbpt_lcao calc
192 omg (np.array): frequency range in eV
194 Output parameters:
195 ------------------
196 sigma (np array): cross section in nm**2
197 """
199 c = 1 / un.alpha # speed of the light in au
200 omg /= un.Ha # to convert from eV to Hartree
201 sigma = 4 * np.pi * omg * p / (c) # bohr**2
202 return sigma * (0.1 * un.Bohr)**2 # nm**2