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