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