Coverage for gpaw/external.py: 93%
203 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-20 00:19 +0000
1"""This module defines different external potentials."""
2import copy
3import warnings
4from typing import Callable, Dict, Optional
6import numpy as np
7from ase.units import Bohr, Ha
9import gpaw.cgpaw as cgpaw
10from gpaw.typing import Array3D
12__all__ = ['ConstantPotential', 'ConstantElectricField', 'CDFTPotential',
13 'PointChargePotential', 'StepPotentialz',
14 'PotentialCollection']
16known_potentials: Dict[str, Callable] = {}
19def _register_known_potentials():
20 from gpaw.bfield import BField
21 known_potentials['CDFTPotential'] = lambda: None # ???
22 known_potentials['BField'] = BField
23 for name in __all__:
24 known_potentials[name] = globals()[name]
27def create_external_potential(name, **kwargs):
28 """Construct potential from dict."""
29 if not known_potentials:
30 _register_known_potentials()
31 return known_potentials[name](**kwargs)
34class ExternalPotential:
35 vext_g: Optional[Array3D] = None
36 vext_q: Optional[Array3D] = None
38 def get_potential(self, gd):
39 """Get the potential on a regular 3-d grid.
41 Will only call calculate_potential() the first time."""
43 if self.vext_g is None:
44 self.calculate_potential(gd)
45 self.vext_g.flags.writeable = False
46 return self.vext_g
48 def get_potentialq(self, gd, pd3):
49 """Get the potential on a regular 3-d grid in real space.
51 Will only call calculate_potential() the first time."""
53 if self.vext_q is None:
54 vext_g = self.get_potential(gd)
55 self.vext_q = pd3.fft(vext_g)
56 self.vext_q.flags.writeable = False
58 return self.vext_q
60 def calculate_potential(self, gd) -> None:
61 raise NotImplementedError
63 def get_name(self) -> str:
64 return self.__class__.__name__
66 def update_potential_pw(self, ham, dens) -> float:
67 v_q = self.get_potentialq(ham.finegd, ham.pd3).copy()
68 eext = ham.pd3.integrate(v_q, dens.rhot_q, global_integral=False)
69 dens.map23.add_to1(ham.vt_Q, v_q)
70 ham.vt_sG[:] = ham.pd2.ifft(ham.vt_Q)
71 if not ham.collinear:
72 ham.vt_xG[1:] = 0.0
73 return eext
75 def update_atomic_hamiltonians_pw(self, ham, W_aL, dens) -> None:
76 vext_q = self.get_potentialq(ham.finegd, ham.pd3)
77 dens.ghat.integrate(ham.vHt_q + vext_q, W_aL)
79 def paw_correction(self, Delta_p, dH_sp) -> None:
80 pass
82 def derivative_pw(self, ham, ghat_aLv, dens) -> None:
83 vext_q = self.get_potentialq(ham.finegd, ham.pd3)
84 dens.ghat.derivative(ham.vHt_q + vext_q, ghat_aLv)
87class NoExternalPotential(ExternalPotential):
88 vext_g = np.zeros((0, 0, 0))
90 def __str__(self):
91 return 'NoExternalPotential'
93 def update_potential_pw(self, ham, dens) -> float:
94 ham.vt_sG[:] = ham.pd2.ifft(ham.vt_Q)
95 if not ham.collinear:
96 ham.vt_xG[1:] = 0.0
97 return 0.0
99 def update_atomic_hamiltonians_pw(self, ham, W_aL, dens):
100 dens.ghat.integrate(ham.vHt_q, W_aL)
102 def derivative_pw(self, ham, ghat_aLv, dens):
103 dens.ghat.derivative(ham.vHt_q, ghat_aLv)
106class ConstantPotential(ExternalPotential):
107 """Constant potential for tests."""
108 def __init__(self, constant=1.0):
109 self.constant = constant / Ha
110 self.name = 'ConstantPotential'
112 def __str__(self):
113 return f'Constant potential: {(self.constant * Ha):.3f} V'
115 def calculate_potential(self, gd):
116 self.vext_g = gd.zeros() + self.constant
118 def todict(self):
119 return {'name': self.name,
120 'constant': self.constant * Ha}
123class ConstantElectricField(ExternalPotential):
124 def __init__(self, strength, direction=[0, 0, 1], tolerance=1e-7):
125 """External constant electric field.
127 strength: float
128 Field strength in V/Ang.
129 direction: vector
130 Polarization direction.
131 """
132 self.strength = strength * Bohr / Ha
133 self.direction_v = np.array(direction, dtype=float)
134 self.direction_v /= np.linalg.norm(self.direction_v)
135 self.field_v = self.strength * self.direction_v
136 self.tolerance = tolerance
137 self.name = 'ConstantElectricField'
139 def __str__(self):
140 return ('Constant electric field: '
141 '({:.3f}, {:.3f}, {:.3f}) V/Ang'
142 .format(*(self.field_v * Ha / Bohr)))
144 def calculate_potential(self, gd):
145 # Note that PW-mode is periodic in all directions!
146 L_c = abs(gd.cell_cv @ self.direction_v)
147 # eps = self.tolerance
148 # assert (L_c < eps).sum() == 2
150 center_v = 0.5 * gd.cell_cv.sum(0)
151 r_gv = gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
152 f_g = (r_gv - center_v) @ self.direction_v
154 # Set potential to zero at boundary of box (important for PW-mode).
155 # Say we have 8 grid points. Instead of a potential like this:
156 #
157 # -4 -3 -2 -1 0 1 2 3
158 #
159 # we want:
160 #
161 # 0 -3 -2 -1 0 1 2 3
163 L = L_c.sum()
164 f_g[abs(abs(f_g) - L / 2) < 1e-5] = 0.0
166 self.vext_g = f_g * self.strength
168 def todict(self):
169 return {'name': self.name,
170 'strength': self.strength * Ha / Bohr,
171 'direction': self.direction_v}
174class ProductPotential(ExternalPotential):
175 def __init__(self, ext_i):
176 self.ext_i = ext_i
178 def calculate_potential(self, gd):
179 self.vext_g = self.ext_i[0].get_potential(gd).copy()
180 for ext in self.ext_i[1:]:
181 self.vext_g *= ext.get_potential(gd)
183 def __str__(self):
184 return '\n'.join(['Product of potentials:'] +
185 [ext.__str__() for ext in self.ext_i])
187 def todict(self):
188 return {'name': self.__class__.__name__,
189 'ext_i': [ext.todict() for ext in self.ext_i]}
192class PointChargePotential(ExternalPotential):
193 def __init__(self, charges, positions=None,
194 rc=0.2, rc2=np.inf, width=1.0):
195 """Point-charge potential.
197 charges: list of float
198 Charges in units of `|e|`.
199 positions: (N, 3)-shaped array-like of float
200 Positions of charges in Angstrom. Can be set later.
201 rc: float
202 Inner cutoff for Coulomb potential in Angstrom.
203 rc2: float
204 Outer cutoff for Coulomb potential in Angstrom.
205 width: float
206 Width for cutoff function for Coulomb part.
208 For r < rc, 1 / r is replaced by a third order polynomial in r^2 that
209 has matching value, first derivative, second derivative and integral.
211 For rc2 - width < r < rc2, 1 / r is multiplied by a smooth cutoff
212 function (a third order polynomium in r).
214 You can also give rc a negative value. In that case, this formula
215 is used::
217 (r^4 - rc^4) / (r^5 - |rc|^5)
219 for all values of r - no cutoff at rc2!
220 """
221 self._dict = dict(name=self.__class__.__name__,
222 charges=charges, positions=positions,
223 rc=rc, rc2=rc2, width=width)
224 self.q_p = np.ascontiguousarray(charges, float)
225 self.rc = rc / Bohr
226 self.rc2 = rc2 / Bohr
227 self.width = width / Bohr
228 if positions is not None:
229 self.set_positions(positions)
230 else:
231 self.R_pv = None
233 if abs(self.q_p).max() < 1e-14:
234 warnings.warn('No charges!')
235 if self.rc < 0. and self.rc2 < np.inf:
236 warnings.warn('Long range cutoff chosen but will not be applied\
237 for negative inner cutoff values!')
239 def todict(self):
240 return copy.deepcopy(self._dict)
242 def __str__(self):
243 return ('Point-charge potential '
244 '(points: {}, cutoffs: {:.3f}, {:.3f}, {:.3f} Ang)'
245 .format(len(self.q_p),
246 self.rc * Bohr,
247 (self.rc2 - self.width) * Bohr,
248 self.rc2 * Bohr))
250 def set_positions(self, R_pv, com_pv=None):
251 """Update positions."""
252 if com_pv is not None:
253 self.com_pv = np.asarray(com_pv) / Bohr
254 else:
255 self.com_pv = None
257 self.R_pv = np.asarray(R_pv) / Bohr
258 self.vext_g = None
260 def _molecule_distances(self, gd):
261 if self.com_pv is not None:
262 return self.com_pv - gd.cell_cv.sum(0) / 2
264 def calculate_potential(self, gd):
265 assert gd.orthogonal
266 self.vext_g = gd.zeros()
268 dcom_pv = self._molecule_distances(gd)
270 cgpaw.pc_potential(gd.beg_c, gd.h_cv.diagonal().copy(),
271 self.q_p, self.R_pv,
272 self.rc, self.rc2, self.width,
273 self.vext_g, dcom_pv)
275 def get_forces(self, calc):
276 """Calculate forces from QM charge density on point-charges."""
277 dens = calc.density
278 F_pv = np.zeros_like(self.R_pv)
279 gd = dens.finegd
280 dcom_pv = self._molecule_distances(gd)
282 cgpaw.pc_potential(gd.beg_c, gd.h_cv.diagonal().copy(),
283 self.q_p, self.R_pv,
284 self.rc, self.rc2, self.width,
285 self.vext_g, dcom_pv, dens.rhot_g, F_pv)
286 gd.comm.sum(F_pv)
287 return F_pv * Ha / Bohr
290class CDFTPotential(ExternalPotential):
291 # Dummy class to make cDFT compatible with new external
292 # potential class ClassName(object):
293 def __init__(self, regions, constraints, n_charge_regions,
294 difference):
296 self.name = 'CDFTPotential'
297 self.regions = regions
298 self.constraints = constraints
299 self.difference = difference
300 self.n_charge_regions = n_charge_regions
302 def todict(self):
303 return {'name': 'CDFTPotential',
304 # 'regions': self.indices_i,
305 'constraints': self.v_i * Ha,
306 'n_charge_regions': self.n_charge_regions,
307 'difference': self.difference,
308 'regions': self.regions}
311class StepPotentialz(ExternalPotential):
312 def __init__(self, zstep, value_left=0, value_right=0):
313 """Step potential in z-direction
315 zstep: float
316 z-value that splits space into left and right [Angstrom]
317 value_left: float
318 Left side (z < zstep) potentential value [eV]. Default: 0
319 value_right: float
320 Right side (z >= zstep) potentential value [eV]. Default: 0
321 """
322 self.value_left = value_left
323 self.value_right = value_right
324 self.name = 'StepPotentialz'
325 self.zstep = zstep
327 def __str__(self):
328 return f'Step potentialz: {self.value_left:.3f} V to '\
329 f'{self.value_right:.3f} V at z={self.zstep}'
331 def calculate_potential(self, gd):
332 r_vg = gd.get_grid_point_coordinates()
333 self.vext_g = np.where(r_vg[2] < self.zstep / Bohr,
334 gd.zeros() + self.value_left / Ha,
335 gd.zeros() + self.value_right / Ha)
337 def todict(self):
338 return {'name': self.name,
339 'value_left': self.value_left,
340 'value_right': self.value_right,
341 'zstep': self.zstep}
344class PotentialCollection(ExternalPotential):
345 def __init__(self, potentials):
346 """Collection of external potentials to be applied
348 potentials: list
349 List of potentials
350 """
351 self.potentials = []
352 for potential in potentials:
353 if isinstance(potential, dict):
354 potential = create_external_potential(
355 potential.pop('name'), **potential)
356 self.potentials.append(potential)
358 def __str__(self):
359 text = 'PotentialCollection:\n'
360 for pot in self.potentials:
361 text += ' ' + pot.__str__() + '\n'
362 return text
364 def calculate_potential(self, gd):
365 self.potentials[0].calculate_potential(gd)
366 self.vext_g = self.potentials[0].vext_g.copy()
367 for pot in self.potentials[1:]:
368 pot.calculate_potential(gd)
369 self.vext_g += pot.vext_g
371 def todict(self):
372 return {'name': 'PotentialCollection',
373 'potentials': [pot.todict() for pot in self.potentials]}
376def static_polarizability(atoms, strength=0.01):
377 """Calculate polarizability tensor
379 atoms: Atoms object
380 strength: field strength in V/Ang
382 Returns
383 -------
384 polarizability tensor:
385 Unit (e^2 Angstrom^2 / eV).
386 Multiply with Bohr * Ha to get (Angstrom^3)
387 """
388 atoms.get_potential_energy()
389 calc = atoms.calc
390 assert calc.parameters.external is None
391 dipole_gs = calc.get_dipole_moment()
393 alpha = np.zeros((3, 3))
394 for c in range(3):
395 axes = np.zeros(3)
396 axes[c] = 1
397 calc.set(external=ConstantElectricField(strength, axes))
398 calc.get_potential_energy()
399 alpha[c] = (calc.get_dipole_moment() - dipole_gs) / strength
400 calc.set(external=None)
402 return alpha.T