Coverage for gpaw/tddft/abc.py: 91%
96 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-19 00:19 +0000
1""" This module implements absorbing boundary conditions to minimize errors
2caused by boundary reflections. The first classes use negative imaginary
3potentials of different shapes which are known to absorb waves. The last
4class PML uses the perfectly matched layer approach.
6 For information about imaginary potentials please see:
8 Takashi Nakatsukasa, Kazuhiro Yabana, J. Chem. Phys. 114, 2550 (2001);
9 doi:10.1063/1.1338527
10 Daniel Neuhasuer and Michael Baer, J. Chem. Phys. 90, 4351 (1989);
11 doi:10.1063/1.456646
13 About PML:
14 Chunxiong Zheng, Journal of Computational Physics,
15 volume 227, Issue 1(Nov 2007) pages 537-556;
16 doi:10.1016/j.jcp.2007.08.004
17"""
19import numpy as np
21from ase.units import Bohr
24class DummyAbsorbingBoundary:
25 """ Virtual and not usable class for abosrbing boundaries."""
26 def __init__(self):
27 self.v_imag = None
28 self.type = None
30 def set_up(self):
31 pass
33 def get_potential_matrix(self):
34 return self.v_imag
37class LinearAbsorbingBoundary(DummyAbsorbingBoundary):
38 """
39 This class uses linear negative imaginary potential for absorption.
40 For distances larger than abc_range, the negative imaginary potential
41 increases linearly. Positive float abc_strength describes the steepness
42 of the slope. High values of abc_strength cause more reflection from
43 the potential and small values do not absorb enough.
45 Parameters:
46 abc_range: Positive float number. Absorbing potential is zero
47 where distance from the middle point or chosen positions
48 in the grid is smaller and starts to increase linearly
49 after this point.
50 abc_strength: Positive float number. A good value for this is usually
51 about 0.01
53 positions: An array of positions in the grid. Each of these points
54 has abc_range of free space around them. If this is left
55 as None, the middle point of the grid is used. You can use
56 the position array of an ASE Atoms object here.
57 """
58 def __init__(self, abc_range, abc_strength, positions=None):
59 DummyAbsorbingBoundary.__init__(self)
60 self.abc_strength = abc_strength
61 self.abc_r = abc_range / Bohr
62 self.positions = positions / Bohr
63 self.type = 'IPOT'
65 def set_up(self, gd):
66 """
67 Creates the potential matrix self.v_imag.
69 Parameters:
71 gd: grid descriptor
72 """
74 assert gd.orthogonal
76 # self.v_imag = np.zeros((gd.n_c[0],gd.n_c[1],gd.n_c[2]),dtype=complex)
77 self.v_imag = gd.zeros(dtype=complex)
79 # If positions array wasn't given, uses the middle point of the grid.
80 if self.positions is None:
81 self.positions = [
82 np.array([
83 gd.N_c[0] * gd.h_cv[0, 0] * 0.5, # middle
84 gd.N_c[1] * gd.h_cv[1, 1] * 0.5,
85 gd.N_c[2] * gd.h_cv[2, 2] * 0.5
86 ])
87 ]
89 for i in range(gd.n_c[0]):
90 x = (i + gd.beg_c[0]) * gd.h_cv[0, 0]
91 for j in range(gd.n_c[1]):
92 y = (j + gd.beg_c[1]) * gd.h_cv[1, 1]
93 for k in range(gd.n_c[2]):
94 z = (k + gd.beg_c[2]) * gd.h_cv[2, 2]
95 position = np.array([x, y, z])
96 # Calculates the distance from the nearest chosen point
97 # in the grid
98 position_vectors = self.positions - position
100 r = np.linalg.norm(position_vectors[0])
101 for vector in position_vectors:
102 if np.linalg.norm(vector) < r:
103 r = np.linalg.norm(vector)
105 if r > self.abc_r:
106 self.v_imag[i][j][k] = \
107 -(0 + 1j) * self.abc_strength * (r - self.abc_r)
110class P4AbsorbingBoundary(DummyAbsorbingBoundary):
111 """
112 The negative imaginary potential used by this class are 4th degree
113 polynomials which are constructed so that the value is zero at the
114 beginning and abc_strength at the end. The derivative is
115 zero at the begininning and zero at the end.
117 Parameters:
118 abc_range: Positive float number. Absorbing potential is zero where
119 distance from the middle point or chosen positions in the
120 grid is smaller and starts to increase after this point.
122 abc_strength: Positive float number.
124 positions: An array of positions in the grid. Each of these points has
125 abc_range of free space around them. If this is left as None,
126 the middle point of the grid is used. You can use the position
127 array of an ASE Atoms object here. You have to define the
128 width parameter if you use this.
130 width: The width of the absorbing layer. If you don't define positions
131 array a value for this is automatically generated so that the
132 width is from abc_range to the end of the grid (works best for
133 cube).
134 If you use the atom-centered model you have to define the
135 width yourself.
136 """
137 def __init__(self, abc_range, abc_strength, positions=None, width=None):
138 DummyAbsorbingBoundary.__init__(self)
139 self.abc_r = abc_range / Bohr
140 self.abc_strength = abc_strength
141 self.positions = positions / Bohr
142 self.width = width / Bohr
143 self.type = 'IPOT'
145 def set_up(self, gd):
147 # self.v_imag = np.zeros((gd.n_c[0],gd.n_c[1],gd.n_c[2]),dtype=complex)
148 self.v_imag = gd.zeros(dtype=complex)
150 # If positions array wasn't given, uses the middle point of the
151 # grid as the center.
152 if self.positions is None:
153 self.positions = [
154 np.array([
155 gd.N_c[0] * gd.h_cv[0, 0] * 0.5, # middle
156 gd.N_c[1] * gd.h_cv[1, 1] * 0.5,
157 gd.N_c[2] * gd.h_cv[2, 2] * 0.5
158 ])
159 ]
161 if self.width is None:
162 self.width = np.linalg.norm(
163 self.positions[0]) / np.sqrt(3) - self.abc_r
165 for i in range(gd.n_c[0]):
166 x = (i + gd.beg_c[0]) * gd.h_cv[0, 0]
167 for j in range(gd.n_c[1]):
168 y = (j + gd.beg_c[1]) * gd.h_cv[1, 1]
169 for k in range(gd.n_c[2]):
170 z = (k + gd.beg_c[2]) * gd.h_cv[2, 2]
171 position = np.array([x, y, z])
173 position_vectors = self.positions - position
175 r = np.linalg.norm(position_vectors[0])
176 for vector in position_vectors:
177 if np.linalg.norm(vector) < r:
178 r = np.linalg.norm(vector)
180 if r > self.abc_r:
181 if r < self.abc_r + self.width:
182 self.v_imag[i][j][k] = (
183 (0 + 1j) *
184 ((np.sqrt(self.abc_strength) -
185 np.sqrt(self.abc_strength) / self.width**2
186 * (r - self.abc_r)**2)**2 -
187 self.abc_strength))
188 else:
189 self.v_imag[i][j][k] = -self.abc_strength * (0 +
190 1j)
193class PML:
194 """
195 Important! You must use the BiCGStab solver or your time propagation
196 will likely crash. Give the parameter solver='BiCGStab' as you create
197 the TDDFT object.
199 Using PML makes the time progation slower so reserve twice more time.
200 As imaginary potential is usually almost equally good, this is mostly
201 for testing (until improved?).
203 And now for something completely different, a Perfectly matched layer.
204 Note that you can't give positions array as parameter for this class.
206 Parameters:
207 abc_range: Starting from the center of the grid, the amount of
208 free space before PML starts to affect.
210 abc_strength: Positive float, good amount should be around 0.1
211 """
212 def __init__(self, abc_range, abc_strength):
214 self.abc_range = abc_range / Bohr
215 self.abc_strength = abc_strength
216 self.type = 'PML'
218 def set_up(self, gd):
219 r"""Set up matrices for PML.
221 Creates the matrixes needed when the PML is applied in tdopers.py
222 when applying time-dependent hamiltonian.
224 Perfectly matched layer is applied as potential Vpml = Tpml-T,
225 Where T = -.5*\nabla^{2}\psi (Use latex for equations)
227 Tpml we get from
229 'A perfectly matched layer approach to the nonlinear Schrodinger wave
230 equations',
231 Journal of Computational Physics,
232 volume 227, Issue 1(Nov 2007) pages 537-556,
233 Author: Chunxiong Zheng
235 T_{pml} = -0.5*(G\nabla G\nabla \psi + G^{2}\nabla^{2}\psi)
237 where G = \frac{1}{1+R\sigma}
239 V_{pml} = -0.5 * (G\nabla G\nabla \psi + (G^{2}-1)\nabla^{2}\psi)
241 This is probably not the most optimal approach and slows the
242 propagation.
244 The matrixes created here are G and the gradients of G in all
245 directions.
247 """
248 R = (0 + 1j) # Complex number, has to be in the first quadrant.
249 self.G = gd.zeros(dtype=complex)
250 self.G[:] = 1.0
252 self.dG = gd.zeros(n=3, dtype=complex)
254 r0 = np.array([
255 gd.N_c[0] * gd.h_cv[0, 0] * 0.5, gd.N_c[1] * gd.h_cv[1, 1] * 0.5,
256 gd.N_c[2] * gd.h_cv[2, 2] * 0.5
257 ]) # middle point
258 for i in range(gd.n_c[0]):
259 x = (i + gd.beg_c[0]) * gd.h_cv[0, 0]
260 for j in range(gd.n_c[1]):
261 y = (j + gd.beg_c[1]) * gd.h_cv[1, 1]
262 for k in range(gd.n_c[2]):
263 z = (k + gd.beg_c[2]) * gd.h_cv[2, 2]
264 position = np.array([x, y, z])
265 r = np.linalg.norm(position - r0)
267 if r > self.abc_range:
268 self.G[i][j][k] = (1.0 + R * self.abc_strength *
269 (r - self.abc_range)**2)**-1.0
270 self.dG[0][i][j][k] = (
271 -(1.0 + R *
272 self.abc_strength *
273 (r - self.abc_range)**2.0)**-2.0 * 2.0 * R *
274 self.abc_strength * (r - self.abc_range) *
275 (x - r0[0]) / r)
277 self.dG[1][i][j][k] = -(
278 1.0 + R * self.abc_strength *
279 (r - self.abc_range)**2.0
280 )**-2.0 * 2.0 * R * self.abc_strength * (
281 r - self.abc_range) * (y - r0[1]) / r
283 self.dG[2][i][j][k] = -(
284 1.0 + R * self.abc_strength *
285 (r - self.abc_range)**2.0
286 )**-2.0 * 2.0 * R * self.abc_strength * (
287 r - self.abc_range) * (z - r0[2]) / r
289 def get_G(self):
290 return self.G
292 def get_dG(self):
293 return self.dG
295 def isPML(self):
296 return True