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