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