Coverage for gpaw/solvation/sjm.py: 68%
644 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1"""
2The solvated jellium method is contained in this module.
3This enables electronically grand-canonical calculations to be calculated,
4typically for simulating electrochemical interfaces.
6This version the the solvated jellium method has been modified
7and includes the Fluctuation Dissipation Theorem.
8"""
10import copy
11import os
12import textwrap
14import ase.io
15import gpaw.mpi
16import numpy as np
17from ase.calculators.calculator import (InputError, Parameters,
18 PropertyNotPresent, equal)
19from ase.parallel import paropen
20from ase.units import Bohr, Ha, _e, kB
21from gpaw import GPAW_NEW, ConvergenceError
22from gpaw.dipole_correction import DipoleCorrection
23from gpaw.fd_operators import Gradient
24from gpaw.hamiltonian import RealSpaceHamiltonian
25from gpaw.io.logger import indent
26from gpaw.jellium import Jellium, JelliumSlab
27from gpaw.solvation.calculator import OldSolvationGPAW
28from gpaw.solvation.cavity import Power12Potential, get_pbc_positions
29from gpaw.solvation.hamiltonian import SolvationRealSpaceHamiltonian
30from gpaw.solvation.poisson import WeightedFDPoissonSolver
31from scipy.ndimage import uniform_filter1d
32from scipy.signal import find_peaks
33from scipy.stats import linregress
36def SJM(*args, **kwargs):
37 """Backwards compatibility ..."""
38 if GPAW_NEW:
39 from gpaw.new.ase_interface import GPAW
40 from gpaw.new.sjm import SJM
41 environment = SJM(cavity=kwargs.pop('cavity'),
42 dielectric=kwargs.pop('dielectric'),
43 interactions=kwargs.pop('interactions', None),
44 **kwargs.pop('sj'))
45 return GPAW(
46 *args, **kwargs,
47 environment=environment)
48 return OldSJM(*args, **kwargs)
51class OldSJM(OldSolvationGPAW):
52 r"""Solvated Jellium method.
53 (Implemented as a subclass of the SolvationGPAW class.)
55 The method allows the simulation of an electrochemical environment, where
56 the potential can be varied by changing the charging (that is, number of
57 electrons) of the system. For this purpose, it allows the usagelof non-
58 neutral periodic slab systems. Cell neutrality is achieved by adding a
59 background charge in the solvent region above the slab
61 Further details are given in https://doi.org/10.1021/acs.jpcc.8b02465
62 If you use this method, we appreciate it if you cite that work.
64 The method can be run in three modes:
66 - Constant charge: The number of excess electrons in the simulation
67 can be directly specified with the 'excess_electrons' keyword,
68 leaving 'target_potential' set to None.
69 - Constant potential: The target potential (expressed as a work
70 function) can be specified with the 'target_potential' keyword.
71 Optionally, the 'excess_electrons' keyword can be supplied to specify
72 the initial guess of the number of electrons.
73 - Constant inner potential: The target potential (expressed as a inner
74 potential) can be specified with the method='CIP' and
75 'target_potential' keywords. The CIP-DFT method is detailled in
76 https://doi.org/10.1038/s41524-023-01184-4
78 By default, this method writes the grand-potential energy to the output;
79 that is, the energy that has been adjusted with `- \mu N` (in this case,
80 `\mu` is the work function and `N` is the excess electrons). This is the
81 energy that is compatible with the forces in constant-potential mode and
82 thus will behave well with optimizers, NEBs, etc. It is also frequently
83 used in subsequent free-energy calculations.
85 Within this method, the potential is expressed as the top-side work
86 function of the slab or the inner potential of the electrode.
87 In both cases, a potential of 0 V_SHE corresponds to target_potential of
88 roughly 4.4 eV. (That is, the user should specify
89 target_potential as 4.4 in this case.) Because this method is
90 attempting to bring either the work function or the inner potential to a
91 target value, the work function and electrostatics need to be
92 well-converged. For this reason, the 'work function' keyword is
93 automatically added to the SCF convergence dictionary with a value of
94 0.001. This can be overriden by the user.
96 All methods requires a dipole correction, and this is turned on
97 automatically, but can be overridden with the poissonsolver keyword.
99 When using method='CIP', mixed Dirichlet/Neumann boundary conditions
100 for the electrostatic potential are used by default.
102 The SJM class takes a single argument, the sj dictionary. All other
103 arguments are fed to the parent SolvationGPAW (and therefore GPAW)
104 calculator.
106 Parameters:
108 sj: dict
109 Dictionary of parameters for the solvated jellium method, whose
110 possible keys are given below.
112 Parameters in sj dictionary:
114 excess_electrons: float
115 Number of electrons added in the atomic system and (with opposite
116 sign) in the background charge region. If the 'target_potential'
117 keyword is also supplied, 'excess_electrons' is taken as an initial
118 guess for the needed number of electrons.
119 target_potential: float
120 The potential that should be reached or kept in the course of the
121 calculation. If set to 'None' (default) a constant-charge
122 calculation based on the value of 'excess_electrons' is performed.
123 Expressed as a work function, on the top side of the slab; see note
124 above.
125 tol: float
126 Tolerance for the deviation of the target potential.
127 If the potential is outside the defined range 'ne' will be
128 changed in order to get inside again. Default: 0.01 V.
129 max_iters: int
130 In constant-potential mode, the maximum number of iterations to try
131 to equilibrate the potential (by varying ne). Default: 10.
132 jelliumregion: dict
133 Parameters regarding the shape of the counter charge region.
134 Implemented keys:
136 'top': float or None
137 Upper boundary of the counter charge region (z coordinate).
138 A positive float is taken as the upper boundary coordiate.
139 A negative float is taken as the distance below the uppper
140 cell boundary. Default: -1.
141 'bottom': float or 'cavity_like' or None
142 Lower boundary of the counter-charge region (z coordinate).
143 A positive float is taken as the lower boundary coordinate.
144 A negative float is interpreted as the distance above the highest
145 z coordinate of any atom; this can be fixed by setting
146 fix_bottom to True. If 'cavity_like' is given the counter
147 charge will take the form of the cavity up to the 'top'.
148 Default: -3.
149 'thickness': float or None
150 Thickness of the counter charge region in Angstroms.
151 Can only be used if start is not 'cavity_like'. Default: None
152 'fix_bottom': bool
153 Do not allow the lower limit of the jellium region to adapt
154 to a change in the atomic geometry during a relaxation.
155 Omitted in a 'cavity_like' jelliumregion.
156 Default: False
157 grand_output: bool
158 Write the grand-potential energy into output files such as
159 trajectory files. Default: True
160 always_adjust: bool
161 Adjust ne again even when potential is within tolerance.
162 This is useful to set to True along with a loose potential
163 tolerance (tol) to allow the potential and structure to be
164 simultaneously optimized in a geometry optimization, for example.
165 Default: False.
166 slope : float or None
167 Initial guess of the slope, in volts per electron, of the relationship
168 between the electrode potential and excess electrons, used to
169 equilibrate the potential.
170 This only applies to the first step, as the slope is calculated
171 internally at subsequent steps. If None, will be guessed based on
172 apparent capacitance of 10 uF/cm^2.
173 max_step : float
174 When equilibrating the potential, if a step results in the
175 potential being further from the target and changing by more than
176 'max_step' V, then the step size is cut in half and we try again. This
177 is to avoid leaving the linear region. You can set to np.inf to turn
178 this off. Default: 2 V.
179 mixer : float
180 Damping for the slope estimate. Because estimating slopes can
181 sometimes be noisy (particularly for small changes, or when
182 positions have also changed), some fraction of the previous slope
183 estimate can be "mixed" with the current slope estimate before the
184 new slope is established. E.g., new_slope = mixer * old_slope +
185 (1. - mixer) * current_slope_estimate. Set to 0 for no damping.
186 Default: 0.5.
187 fdt : bool or dict
188 Keyboard for switching on/off the computation of the charge using
189 the fluctuation dissipation theorem
190 Default: False.
191 If set to True, use standard parameters for the FDT calculation.
192 If dict, the following keys are implemented:
194 'dt': float
195 Time step for the FDT calculation in fs. Default: 0.5.
196 'po_time': float
197 Relaxation-time constant of potentiostat. Default: 100.
198 'th_temp': float
199 Thermal temperature for the FDT calculation in K. Default: 300.
201 slope_regression_depth : int
202 Number of previous attempts to use for the slope regression.
203 Default: 4.
204 pot_ref: 'wf' or 'CIP'
205 potential reference scale
206 wf: original SJM using workfunction for the absolute potential
207 CIP: use inner potential as the abs. el. pot.
208 cip: dict, parameters when using CIP-DFT
209 inner_region: list of floats: [bottom, top]
210 bottom is the starting point for computing the inner potential,
211 likewise top is the ending point.
212 filter: int
213 number of points for smoothing the el.stat. pot.
214 autoinner: dict with {'nlayers':int, 'threshold':0.001}
215 if inner_region is given, autoinner is automatically disabled
216 nlayers: number of layers
217 threshold: Required threshold of peaks, the innerpotential
218 difference at neighboring grid points
219 mu_pzc: float
220 Fermi level at potential of zero charge
221 Sets the reference scale for the absolute potential level for CIP
222 phi_pzc: float
223 Use only with method='CIP', the value corresponds to the inner
224 potential of the neutral electrode
226 Special SJM methods (in addition to those of GPAW/SolvationGPAW):
228 get_electrode_potential
229 Returns the potential of the simulated electrode, in V, relative
230 to the vacuum.
232 write_sjm_traces
233 Write traces of quantities like electrostatic potential or cavity
234 to disk.
236 """
237 implemented_properties = ['energy', 'forces', 'stress', 'dipole',
238 'magmom', 'magmoms', 'excess_electrons',
239 'electrode_potential']
240 _sj_default_parameters = Parameters(
241 {'excess_electrons': 0.,
242 'jelliumregion': {'top': -1.,
243 'bottom': -3.,
244 'thickness': None,
245 'fix_bottom': False},
246 'target_potential': None,
247 'pot_ref': 'wf',
248 'tol': 0.01,
249 'always_adjust': False,
250 'grand_output': True,
251 'max_iters': 10,
252 'max_step': 2.,
253 'slope': None,
254 'mixer': 0.5,
255 'fdt': False,
256 'previous_electrons': [],
257 'previous_potentials': [],
258 'slope_regression_depth': 4,
259 'cip': {'autoinner': {'nlayers': None,
260 'threshold': 0.0001},
261 'inner_region': None,
262 'mu_pzc': None,
263 'phi_pzc': None,
264 'filter': 10}})
266 _sj_default_parameters.update({'dirichlet': False})
267 default_parameters = copy.deepcopy(OldSolvationGPAW.default_parameters)
268 default_parameters.update({'poissonsolver': {'dipolelayer': 'xy'}})
269 default_parameters['convergence'].update({'work function': 0.001})
270 default_parameters.update({'sj': _sj_default_parameters})
272 def __init__(self, restart=None, **kwargs):
274 deprecated_keys = ['ne', 'potential', 'write_grandcanonical_energy',
275 'potential_equilibration_mode', 'dpot',
276 'max_pot_deviation', 'doublelayer', 'verbose']
277 msg = ('{:s} is no longer a supported keyword argument for the SJM '
278 'class. All SJM arguments should be sent in via the "sj" '
279 'dict.')
280 for key in deprecated_keys:
281 if key in kwargs:
282 raise InputError(textwrap.fill(msg.format(key)))
284 # Note the below line calls self.set().
285 OldSolvationGPAW.__init__(self, restart, **kwargs)
287 def set(self, **kwargs):
288 """Change parameters for calculator.
290 It differs from the standard `set` function in two ways:
291 - Keywords in the `sj` dictionary are handled.
292 - It does not reinitialize and delete `self.wfs` if only the
293 background charge is changed.
295 """
296 p = self.parameters['sj']
298 # We handle 'sj' and 'background_charge' internally; passing to GPAW's
299 # set function will trigger a deletion of the density, etc.
300 sj_changes = kwargs.pop('sj', {})
301 try:
302 sj_changes = {key: value for key, value in sj_changes.items()
303 if not equal(value, p[key])}
304 except KeyError:
305 raise InputError(
306 'Unexpected key(s) provided to sj dict. '
307 'Keys provided were "{}". '
308 'Only keys allowed are "{}".'
309 .format(', '.join(sj_changes),
310 ', '.join(self.default_parameters['sj'])))
311 self.fill_cip_keywords(p.cip)
313 if p.pot_ref == 'CIP':
314 p['dirichlet'] = True
316 if p.cip['mu_pzc'] is None or p.cip['phi_pzc'] is None:
317 p.cip['mu_pzc'] = 0
318 p.cip['phi_pzc'] = 0
319 msg = ('Warning: a CIP calculation has been activated '
320 'but mu_pzc and/or phi_pzc was none. This is fine '
321 'for CIP calibration but meaningful references '
322 'must be provided for production calculations\n')
324 if p.cip['inner_region'] is None and p.cip['autoinner'] is None:
325 raise RuntimeError("The inner region cannot be none" +
326 "when using inner potential as the" +
327 "reference. Please, set up" +
328 "either bottom/top values to define the" +
329 "electrode bulk or set autoinner to True")
331 if p.cip.get('inner_region') and p.cip.get('autoinner'):
332 raise RuntimeError("Only inner_region or autoinner" +
333 "can be set to define the inner potential")
335 if p.cip['inner_region'] is not None:
336 p.cip['inner_region'] = np.array(p.inner_region)
337 p.cip['autoinner'] = None
338 else:
339 assert p.cip['autoinner']['nlayers'] is not None
341 p.update(sj_changes)
343 background_charge = kwargs.pop('background_charge', None)
344 kwargs['_set_ok'] = True
345 OldSolvationGPAW.set(self, **kwargs)
347 # parent_changed checks if GPAW needs to be reinitialized
348 # The following key do not need reinitialization
349 parent_changed = False
350 for key in kwargs:
351 if key not in ['mixer', 'verbose', 'txt', 'hund', 'random',
352 'eigensolver', 'convergence', 'fixdensity',
353 'maxiter', '_set_ok']:
354 parent_changed = True
356 if len(sj_changes):
357 if self.wfs is None:
358 self.log('Non-default Solvated Jellium parameters:')
359 else:
360 self.log('Changed Solvated Jellium parameters:')
361 self.log.print_dict({i: p[i] for i in sj_changes})
362 self.log()
364 if 'target_potential' in sj_changes and p.target_potential is not None:
365 # If target potential is changed by the user and the slope is
366 # known, a step towards the new potential is taken right away.
367 try:
368 true_potential = self.get_electrode_potential()
369 # TypeError is needed for the case of starting from a gpw
370 # file and changing the target potential at the start.
371 except (AttributeError, TypeError):
372 pass
373 else:
374 if self.atoms and p.slope:
376 p.excess_electrons += ((p.target_potential -
377 true_potential) / p.slope)
378 self.log('Number of electrons changed to {:.4f} based '
379 'on slope of {:.4f} V/electron.'
380 .format(p.excess_electrons, p.slope))
382 if (any(key in ['target_potential', 'excess_electrons',
383 'jelliumregion'] for key in sj_changes) and
384 not parent_changed):
385 self.results = {}
386 # SolvationGPAW will not reinitialize anymore if only
387 # 'sj' keywords are set. The lines below will reinitialize and
388 # apply the changed charges.
389 if self.atoms:
390 self.set(background_charge=self._create_jellium())
392 if 'tol' in sj_changes:
393 try:
394 true_potential = self.get_electrode_potential()
395 except (AttributeError, TypeError):
396 pass
397 else:
398 msg = ('Potential tolerance changed to {:1.4f} V: '
399 .format(p.tol))
400 if abs(true_potential - p.target_potential) > p.tol:
401 msg += 'new calculation required.'
402 self.results = {}
403 if self.atoms and p.slope is not None:
404 p.excess_electrons += ((p.target_potential -
405 true_potential) / p.slope)
406 self.set(background_charge=self._create_jellium())
407 msg += ('\n Excess electrons changed to {:.4f} based '
408 'on slope of {:.4f} V/electron.'
409 .format(p.excess_electrons, p.slope))
410 else:
411 msg += 'already within tolerance.'
412 self.log(msg)
414 if background_charge:
415 # background_charge is a GPAW parameter that we handle internally,
416 # as it contains the jellium countercharge. Note if a user tries to
417 # specify an *additional* background charge this will probably
418 # conflict, but we know of no such use cases.
419 if self.wfs is None:
420 kwargs.update({'background_charge': background_charge,
421 '_set_ok': True})
422 OldSolvationGPAW.set(self, **kwargs)
423 else:
424 if parent_changed:
425 self.density = None
426 else:
427 if self.density.background_charge:
428 self.density.background_charge = background_charge
429 self.density.background_charge.set_grid_descriptor(
430 self.density.finegd)
431 self._quick_reinitialization()
433 self.wfs.nvalence = self.setups.nvalence + p.excess_electrons
434 self.log('Number of valence electrons is now {:.5f}'
435 .format(self.wfs.nvalence))
437 if self.parameters['sj']['fdt'] is True:
438 # Default parameters if fdt is True
439 self.parameters['sj']['fdt'] = {
440 'dt': 0.5,
441 'po_time': 100.0,
442 'th_temp': 300.0}
443 elif isinstance(self.parameters['sj']['fdt'], dict):
444 # If fdt is a dict, ensure the dictionary is complete
445 fdt_dict = self.parameters['sj']['fdt']
446 self.parameters['sj']['fdt'] = {
447 'dt': fdt_dict.get('dt', 0.5),
448 'po_time': fdt_dict.get('po_time', 100.0),
449 'th_temp': fdt_dict.get('th_temp', 300.0)}
451 def _quick_reinitialization(self):
452 """Minimal reinitialization of electronic-structure stuff when only
453 background charge changes."""
454 if self.density.nct_G is None:
455 self.initialize_positions()
457 self.density.reset()
458 self._set_atoms(self.atoms)
459 self.density.mixer.reset()
460 self.wfs.initialize(self.density, self.hamiltonian,
461 self.spos_ac)
462 self.wfs.eigensolver.reset()
463 if self.scf:
464 self.scf.reset()
466 def calculate(self, atoms=None, properties=['energy'],
467 system_changes=['cell']):
468 """Perform an electronic structure calculation, with either a
469 constant number of electrons or a target potential, as requested by
470 the user in the 'sj' dict."""
472 if atoms and not self.atoms:
473 # Need to be set before ASE's Calculator.calculate gets to it.
474 self.atoms = atoms.copy()
476 if len(system_changes) == 0 and len(self.results) > 0:
477 # Potential is already equilibrated.
478 OldSolvationGPAW.calculate(self, atoms, properties, system_changes)
479 return
481 self.log('Solvated jellium method (SJM) calculation:')
483 p = self.parameters['sj']
485 if p.target_potential is None:
486 self.log('Constant-charge calculation with {:.5f} excess '
487 'electrons'.format(p.excess_electrons))
488 # Background charge is set here, not earlier, because atoms needed.
489 self.set(background_charge=self._create_jellium())
490 OldSolvationGPAW.calculate(self, atoms, ['energy'], system_changes)
491 self.log('Potential found to be {:.5f} V (with {:+.5f} '
492 'electrons)'.format(self.get_electrode_potential(),
493 p.excess_electrons))
494 else:
495 self.log(' Constant-potential mode.')
496 self.log(' Target potential: {:.5f} +/- {:.5f}'
497 .format(p.target_potential, p.tol))
498 self.log(' Initial guess of excess electrons: {:.5f}'
499 .format(p.excess_electrons))
500 if 'workfunction' in self.parameters.convergence:
501 if self.parameters.convergence['workfunction'] >= p.tol:
502 msg = ('Warning: it appears that your work function '
503 'convergence criterion ({:g}) is higher than your '
504 'desired potential tolerance ({:g}). This may lead '
505 'to issues with potential convergence.'
506 .format(self.parameters.convergence['workfunction'],
507 p.tol))
508 self.log(textwrap.fill(msg))
509 self._equilibrate_potential(atoms, system_changes)
510 if properties != ['energy']:
511 # The equilibration loop only calculated energy, to save
512 # unnecessary computations (mostly of forces) in the loop.
513 OldSolvationGPAW.calculate(self, atoms, properties, [])
515 # Note that grand-potential energies were assembled in summary,
516 # which in turn was called by GPAW.calculate.
518 if p.grand_output:
519 self.results['energy'] = self.omega_extrapolated * Ha
520 self.results['free_energy'] = self.omega_free * Ha
521 self.log('Grand-potential energy was written into results.\n')
522 else:
523 self.log('Canonical energy was written into results.\n')
525 self.results['excess_electrons'] = p.excess_electrons
526 self.results['electrode_potential'] = self.get_electrode_potential()
527 self.log.fd.flush()
529 def _equilibrate_potential(self, atoms, system_changes):
530 """Adjusts the number of electrons until the potential reaches the
531 desired value."""
532 p = self.parameters['sj']
533 iteration = 0
535 rerun = False
536 while iteration <= p.max_iters:
537 self.log('Attempt {:d} to equilibrate potential to {:.3f} +/-'
538 ' {:.3f} V'
539 .format(iteration, p.target_potential, p.tol))
540 self.log('Current guess of excess electrons: {:+.5f}\n'
541 .format(p.excess_electrons))
542 if iteration == 1:
543 self.timer.start('Potential equilibration loop')
544 # We don't want SolvationGPAW to see any more system
545 # changes, like positions, after attempt 0.
546 system_changes = []
548 if any([iteration, rerun, 'positions' in system_changes]):
549 self.set(background_charge=self._create_jellium())
551 # Do the calculation.
552 OldSolvationGPAW.calculate(self, atoms, ['energy'], system_changes)
553 true_potential = self.get_electrode_potential()
554 self.log()
555 msg = (f'Potential found to be {true_potential:.5f} V (with '
556 f'{p.excess_electrons:+.5f} excess electrons, attempt '
557 f'{iteration:d}/{p.max_iters:d}')
558 msg += ' rerun).' if rerun else ').'
559 self.log(msg, flush=True)
561 # Check if we took too big of a step, that moved us further
562 # from the target potential rather than closer. This can
563 # happen when large geometric changes happened in the
564 # last step, the slope is noisy or, at very low
565 # workfunctions, where the electron density starts to spill out
566 # into the vacuum and the slope is unreliable.
567 # The rerun can happen multiple times if needed and the stepsize
568 # will be reduced by factor of 2 every time.
569 # The rerun is disabled if the FDT is used.
571 if len(p.previous_potentials):
573 stepsize = abs(true_potential - p.previous_potentials[-1])
575 if (stepsize > p.max_step and
576 abs(p.previous_potentials[-1] - p.target_potential) <
577 abs(true_potential - p.target_potential)) and not p.fdt:
578 self.log('Step resulted in a potential change of '
579 f'{stepsize:.2f} V, larger than max_step '
580 f'({p.max_step:.2f} V) and\n surpassed the'
581 ' target potential by a dangerous amount.\n'
582 ' The step is rejected and the change in'
583 ' excess_electrons will be halved.')
585 if p.fdt:
586 rerun = False
587 else:
588 pe, ce = p.previous_electrons[-1], p.excess_electrons
589 if abs(pe - ce) < 1e-5:
590 msg = ('Step size is too small to be halved in '
591 'rerun. To avoid this try to change your '
592 'initial guess of excess electrons. '
593 'Potential equilibration failed.')
594 raise PotentialConvergenceError(msg)
596 p.excess_electrons = (pe + ce) / 2.
598 rerun = True
599 continue # back to while
601 # Increase iteration count.
602 iteration += 1
603 rerun = False
605 # Store attempt and calculate slope.
606 p.previous_electrons.append(float(p.excess_electrons))
607 p.previous_potentials.append(float(true_potential))
609 # The following solves a bug, where the code would crash if the
610 # user sets the right number of electrons to reach the target
611 # potential in the first iteration and then changes the target
612 # potential. The code would crash because the slope has not been
613 # calculated yet and so no step is taken towards the new potential.
614 # As two equal charges are added to p.previous_electrons, the
615 # regression of the slope will fail.
616 if len(p.previous_electrons) > 1:
617 if not p.previous_electrons[-2] - p.previous_electrons[-1]:
618 del p.previous_electrons[-2], p.previous_potentials[-2]
620 if len(p.previous_electrons) > 1:
621 slope = _calculate_slope(p.previous_electrons,
622 p.previous_potentials,
623 p.slope_regression_depth)
625 nreg = len(p.previous_electrons[-p.slope_regression_depth:])
626 self.log(f'Slope regressed from last {nreg:d} attempts is '
627 f'{slope:.4f} V/electron,')
628 area = np.linalg.det(atoms.cell[:2, :2])
629 # get capacitance in muF/cm^2
630 capacitance = - _e * 1e22 / (area * slope)
631 self.log(f'or apparent capacitance of {capacitance:.4f} '
632 'muF/cm^2')
634 if p.slope is not None:
635 p.slope = p.mixer * p.slope + (1. - p.mixer) * slope
636 self.log(f'After mixing with {p.mixer:.2f}, new slope is '
637 f'{p.slope:.4f} V/electron.')
638 else:
639 p.slope = slope
641 self.log.flush()
643 # Check if we're equilibrated and exit if always_adjust is False.
644 if abs(true_potential - p.target_potential) < p.tol and not p.fdt:
645 self.log('Potential is within tolerance. Equilibrated.')
646 if iteration >= 2:
647 self.timer.stop('Potential equilibration loop')
648 if not p.always_adjust:
649 return
651 # Guess slope if we don't have enough information yet.
652 if p.slope is None or (p.slope > 0. and p.fdt):
653 area = np.linalg.det(atoms.cell[:2, :2])
654 p.slope = -1.6022e3 / (area * 10.)
655 if p.fdt:
656 self.log('Positive slope! Guessing a slope of '
657 f'{p.slope:.4f} corresponding\nto an apparent '
658 'capacitance of 10 muF/cm^2.')
659 else:
660 self.log('No slope provided, guessing a slope of '
661 f'{p.slope:.4f} corresponding\nto an apparent '
662 'capacitance of 10 muF/cm^2.')
664 if p.fdt:
665 fdt_dict = p['fdt']
666 dt = fdt_dict['dt']
667 po_time = fdt_dict['po_time']
668 th_temp = fdt_dict['th_temp']
670 rn = np.random.standard_normal(1)
672 self.world.broadcast(rn, 0)
673 # set capacitance again
674 area = np.linalg.det(atoms.cell[:2, :2])
675 if abs(p.slope) < 1e-10:
676 raise ValueError(
677 "Slope cannot be zero when calculating capacitance.")
679 capacitance = - _e * 1e22 / (area * p.slope)
681 p.excess_electrons += (
682 capacitance * (true_potential - p.target_potential)
683 * (1 - np.exp(-dt / po_time))
684 + rn[0] * np.sqrt(
685 kB * th_temp * capacitance
686 * (1 - np.exp(-2 * dt / po_time))
687 )
688 )
689 self.log(
690 f'Number of electrons is {p.excess_electrons:.4f} '
691 f'using the FDT, with slope of {p.slope:.4f} V/electron '
692 f'and capacitance of ({capacitance:.4f} muF/cm2).'
693 )
695 else:
696 p.excess_electrons += ((p.target_potential - true_potential)
697 / p.slope)
698 self.log(
699 f'Number of electrons changed to {p.excess_electrons:.4f}'
700 f' based on slope of {p.slope:.4f} V/electron.')
702 # Check if we're equilibrated and exit if always_adjust is True.
703 if (abs(true_potential - p.target_potential) < p.tol
704 and p.always_adjust) or p.fdt:
705 return
707 msg = (f'Potential could not be reached after {iteration - 1:d} '
708 'iterations. This may indicate your workfunction is noisier '
709 'than your potential tol. You may try setting the '
710 'convergence["workfunction"] keyword. The last values of '
711 'excess_electrons and the potential are listed below; '
712 'plotting them could give you insight into the problem.')
713 msg = textwrap.fill(msg) + '\n'
714 for n, p in zip(p.previous_electrons, p.previous_potentials):
715 msg += f'{n:+.6f} {p:.6f}\n'
716 self.log(msg, flush=True)
717 raise PotentialConvergenceError(msg)
719 def write_sjm_traces(self, path='sjm_traces', style='z',
720 props=('potential', 'cavity', 'background_charge')):
721 """Write traces of quantities in `props` to file on disk; traces will
722 be stored within specified path. Default is to save as vertical traces
723 (style 'z'), but can also save as cube (specify `style='cube'`)."""
724 grid = self.density.finegd
725 data = {'cavity': self.hamiltonian.cavity.g_g,
726 'background_charge': self.density.background_charge.mask_g,
727 'potential': (self.hamiltonian.vHt_g * Ha -
728 self.get_fermi_level())}
729 if not os.path.exists(path) and gpaw.mpi.world.rank == 0:
730 os.makedirs(path)
731 for prop in props:
732 if style == 'z':
733 _write_trace_in_z(grid, data[prop], prop + '.txt', path)
734 elif style == 'cube':
735 _write_property_on_grid(grid, data[prop], self.atoms,
736 prop + '.cube', path)
738 def summary(self):
739 """Writes summary information to the log file.
740 This varies from the implementation in gpaw.calculator.GPAW by the
741 inclusion of grand potential quantities."""
742 # Standard GPAW summary.
743 self.hamiltonian.summary(self.wfs, self.log)
744 # Add grand-canonical terms.
745 p = self.parameters['sj']
746 self.log()
747 mu_N = -self.get_electrode_potential()
748 mu_N *= p.excess_electrons / Ha
749 self.omega_free = self.hamiltonian.e_total_free - mu_N
750 self.omega_extrapolated = self.hamiltonian.e_total_extrapolated - mu_N
751 self.log('Legendre-transformed energies (grand potential, '
752 'Omega = E - N mu)')
753 self.log(' N (excess electrons): {:+11.6f}'
754 .format(p.excess_electrons))
755 if p['pot_ref'] == 'wf':
756 self.log(' mu (-workfunction, eV): {:+11.6f}'
757 .format(-self.get_electrode_potential()))
758 elif p['pot_ref'] == 'CIP':
759 self.log('Electrode potential from inner potential: {:+11.6f} [eV]'
760 .format(self.get_electrode_potential()))
761 self.log('The absolute inner potential is: {:+11.6f} [eV]'
762 .format(self.get_inner_potential(self.atoms,
763 p['cip']['inner_region'])))
765 self.log(' (Grand) free energy: {:+11.6f}'
766 .format(Ha * self.omega_free))
767 self.log(' (Grand) extrapolated: {:+11.6f}'
768 .format(Ha * self.omega_extrapolated))
769 self.log()
770 # Back to standard GPAW summary.
771 self.density.summary(self.atoms, self.results.get('magmom', 0.0),
772 self.log)
773 self.wfs.summary(self.log)
774 self._print_gapinfo()
775 self.log.fd.flush()
777 def _create_jellium(self):
778 """Creates the counter charge according to the user's specs."""
779 atoms = self.atoms
781 p = self.parameters['sj']
782 jellium = p['jelliumregion']
783 defaults = self.default_parameters['sj']['jelliumregion']
785 # Populate missing keywords.
786 missing = {'top': None,
787 'bottom': None,
788 'thickness': None,
789 'fix_bottom': False}
790 for key in missing:
791 if key not in jellium:
792 jellium[key] = missing[key]
794 # Catch incompatible specifications.
795 if jellium.get('thickness') and jellium['bottom'] == 'cavity_like':
796 raise InputError("With a cavity-like counter charge only the "
797 "keyword 'top' (not 'thickness') allowed.")
799 # We need 2 of the 3 "specifieds" below.
800 specifieds = [jellium['top'] is not None,
801 jellium['bottom'] is not None,
802 jellium['thickness'] is not None]
804 if sum(specifieds) == 3:
805 raise InputError('The jellium region has been overspecified.')
806 if sum(specifieds) == 0:
807 top = defaults['top']
808 bottom = defaults['bottom']
809 thickness = defaults['thickness']
810 if sum(specifieds) == 2:
811 top = jellium['top']
812 bottom = jellium['bottom']
813 thickness = jellium['thickness']
814 if specifieds == [True, False, False]:
815 top = jellium['top']
816 bottom = defaults['bottom']
817 thickness = None
818 elif specifieds == [False, True, False]:
819 top = defaults['top']
820 bottom = jellium['bottom']
821 thickness = None
822 elif specifieds == [False, False, True]:
823 top = None
824 bottom = defaults['bottom']
825 thickness = jellium['thickness']
827 # Deal with negative specifications of upper and lower limits;
828 # as well as fixed/free lower limit.
829 if top is not None and top < 0.:
830 top = atoms.cell[2][2] + top
831 if bottom not in [None, 'cavity_like'] and bottom < 0.:
832 bottom = (max(atoms.positions[:, 2]) - bottom)
833 if jellium['fix_bottom'] is True:
834 try:
835 bottom = self._fixed_lower_limit
836 except AttributeError:
837 self._fixed_lower_limit = bottom
839 # Use thickness if needed.
840 if top is None:
841 top = bottom + thickness
842 if bottom is None:
843 bottom = top - thickness
845 # Catch unphysical limits.
846 if top > atoms.cell[2][2]:
847 raise InputError('The upper limit of the jellium region lies '
848 'outside of unit cell. If you did not set it '
849 'manually, increase your unit cell size or '
850 'translate the atomic system down along the '
851 'z-axis.')
852 if bottom != 'cavity_like':
853 if bottom > top:
854 raise InputError('Your jellium region has a bottom at {:.3f}'
855 ' AA, which is above the top at {:.3f} AA.'
856 .format(bottom, top))
858 # Finally, make the jellium.
859 if bottom == 'cavity_like':
860 if self.hamiltonian is None:
861 self.initialize(atoms)
862 self.set_positions(atoms)
863 g_g = self.hamiltonian.cavity.g_g.copy()
864 self.wfs = None
865 self.density = None
866 self.hamiltonian = None
867 self.initialized = False
868 else:
869 self.set_positions(atoms)
870 # XXX If you start with a fixed bottom, run a calc,
871 # then hot-switch to cavity_like it crashes here. Not sure
872 # that's common enough to worry about.
873 g_g = self.hamiltonian.cavity.g_g
874 self.log('Jellium counter-charge defined with:\n'
875 ' Bottom: cavity-like\n'
876 ' Top: {:7.3f} AA\n'
877 ' Charge: {:5.4f}\n'
878 .format(top, p.excess_electrons))
879 return CavityShapedJellium(charge=p.excess_electrons,
880 g_g=g_g,
881 z2=top)
883 self.log('Jellium counter-charge defined with:\n'
884 ' Bottom: {:7.3f} AA\n'
885 ' Top: {:7.3f} AA\n'
886 ' Charge: {:5.4f}\n'
887 .format(bottom, top, p.excess_electrons))
888 return JelliumSlab(charge=p.excess_electrons,
889 z1=bottom,
890 z2=top)
892 def get_electrode_potential(self, pot_ref=None,
893 return_referenced=True):
894 """Returns the potential of the simulated electrode, in V, relative
895 to the vacuum. This comes directly from the work function."""
897 if pot_ref is None:
898 pot_ref = self.parameters.sj['pot_ref']
900 if pot_ref == 'wf':
901 try:
902 return Ha * self.hamiltonian.get_workfunctions(self.wfs)[1]
903 except TypeError:
904 # Error happens on freshly-opened *.gpw file.
905 if 'electrode_potential' in self.results:
906 return self.results['electrode_potential']
907 else:
908 msg = ('Electrode potential could not be read. Make sure a'
909 'DFT calculation has been performed before reading '
910 'the potential.')
911 raise PropertyNotPresent(textwrap.fill(msg))
913 elif pot_ref == 'CIP':
914 cip = self.parameters.sj.cip
915 inner_potential = self.get_inner_potential(self.atoms)
916 if return_referenced is True:
917 inner_potential = - (cip['mu_pzc'] - cip['phi_pzc'] +
918 inner_potential)
919 return inner_potential
921 return cip['phi_pzc'] - inner_potential
923 def get_inner_potential(self, atoms, z=None):
925 cip = self.parameters.sj.cip
926 self.fill_cip_keywords(cip)
928 el = self.hamiltonian.vHt_g * Ha
929 gd = self.hamiltonian.finegd
930 elstat = gd.collect(el, broadcast=True)
931 el_stat_z = elstat.mean(0).mean(0)
932 smooth_elstat_z = uniform_filter1d(el_stat_z, size=cip['filter'])
934 if cip['inner_region'] is not None:
935 z_top = (np.abs(gd.coords(2) - (z[1] / Bohr))).argmin()
936 z_bottom = (np.abs(gd.coords(2) - (z[0] / Bohr))).argmin()
937 el_av = np.average(smooth_elstat_z[z_bottom:z_top])
938 else:
939 # find minima of elstat potentials (position of nuclei)
940 peaks, _ = find_peaks(-smooth_elstat_z,
941 threshold=cip['autoinner']['threshold'])
942 # choose maxima of elstatpot using the number of layers
943 if (cip['autoinner']['nlayers'] % 2) == 0:
944 # even number of peaks
945 b = int(cip['autoinner']['nlayers'] / 2)
946 bottom = peaks[b - 1]
947 top = peaks[b]
948 el_av = np.average(smooth_elstat_z[bottom:top])
949 else:
950 # odd number of peaks
951 c = int((cip['autoinner']['nlayers'] - 1) / 2)
952 bottom = peaks[c - 1]
953 top = peaks[c + 1]
954 el_av = np.average(smooth_elstat_z[bottom:top])
956 return el_av
958 def fill_cip_keywords(self, cip):
960 defaults_auto = self.default_parameters['sj']['cip']['autoinner']
961 missing = {'nlayers': None,
962 'threshold': None}
963 for key in missing:
964 if key not in cip['autoinner']:
965 cip['autoinner'][key] = defaults_auto[key]
967 defaults_cip = self.default_parameters['sj']['cip']
968 missing = {'inner_region': None,
969 'mu_pzc': None,
970 'phi_pzc': None,
971 'filter': None}
973 for key in missing:
974 if key not in cip:
975 cip[key] = defaults_cip[key]
977 def initialize(self, atoms=None, reading=False):
978 """Inexpensive initialization.
980 This catches CavityShapedJellium, which GPAW's initialize does not
981 know how to handle on restart. We delete and it will be recreated by
982 the _create_jellium method when needed.
983 """
984 background_charge = self.parameters['background_charge']
985 if isinstance(background_charge, dict):
986 if 'z1' in background_charge:
987 if background_charge['z1'] == 'cavity_like':
988 self.parameters['background_charge'] = None
989 OldSolvationGPAW.initialize(self=self, atoms=atoms, reading=reading)
991 def create_hamiltonian(self, realspace, mode, xc):
992 """This differs from SolvationGPAW's create_hamiltonian method by the
993 ability to use dipole corrections."""
994 if not realspace:
995 raise NotImplementedError(
996 'SJM does not support calculations in reciprocal space yet'
997 ' due to a lack of an implicit solvent module.')
999 dens = self.density
1001 self.hamiltonian = SJM_RealSpaceHamiltonian(
1002 *self.stuff_for_hamiltonian,
1003 gd=dens.gd, finegd=dens.finegd,
1004 nspins=dens.nspins,
1005 collinear=dens.collinear,
1006 setups=dens.setups,
1007 timer=self.timer,
1008 xc=xc,
1009 world=self.world,
1010 redistributor=dens.redistributor,
1011 vext=self.parameters.external,
1012 psolver=self.parameters.poissonsolver,
1013 stencil=mode.interpolation,
1014 dirichlet=self.parameters['sj']['dirichlet'])
1016 xc.set_grid_descriptor(self.hamiltonian.finegd)
1019def _write_trace_in_z(grid, property, name, dir):
1020 """Writes out a property (like electrostatic potential, cavity, or
1021 background charge) as a function of the z coordinate only. `grid` is the
1022 grid descriptor, typically self.density.finegd. `property` is the property
1023 to be output, on the same grid."""
1024 property = grid.collect(property, broadcast=True)
1025 property_z = property.mean(0).mean(0)
1026 with paropen(os.path.join(dir, name), 'w') as f:
1027 for i, val in enumerate(property_z):
1028 f.write(f'{(i + 1) * grid.h_cv[2][2] * Bohr:f} {val:1.8f}\n')
1031def _write_property_on_grid(grid, property, atoms, name, dir):
1032 """Writes out a property (like electrostatic potential, cavity, or
1033 background charge) on the grid, as a cube file. `grid` is the
1034 grid descriptor, typically self.density.finegd. `property` is the property
1035 to be output, on the same grid."""
1036 property = grid.collect(property, broadcast=True)
1037 ase.io.write(os.path.join(dir, name), atoms, data=property)
1040def _calculate_slope(previous_electrons, previous_potentials, n_prev_pot):
1041 """Calculates the slope of potential versus number of electrons;
1042 regresses based on (up to) last four data points to smooth noise."""
1043 # debug
1045 ans = linregress(previous_electrons[-n_prev_pot:],
1046 previous_potentials[-n_prev_pot:])
1047 return ans[0]
1050class SJMPower12Potential(Power12Potential):
1051 r"""Inverse power-law potential.
1052 Inverse power law potential for SJM, inherited from the
1053 Power12Potential of gpaw.solvation. This is a 1/r^{12} repulsive
1054 potential taking the value u0 at the atomic radius. In SJM one also has the
1055 option of removing the solvent from the electrode backside and adding
1056 ghost plane/atoms to remove the solvent from the electrode-water interface.
1058 Parameters:
1060 atomic_radii: function
1061 Callable mapping an ase.Atoms object to an iterable of atomic radii
1062 in Angstroms. If not provided, defaults to van der Waals radii.
1063 u0: float
1064 Strength of the potential at the atomic radius in eV.
1065 Defaults to 0.18 eV, the best-fit value for water from Held &
1066 Walter.
1067 pbc_cutoff: float
1068 Cutoff in eV for including neighbor cells in a calculation with
1069 periodic boundary conditions.
1070 H2O_layer: bool, int or str
1071 True: Exclude the implicit solvent from the interface region
1072 between electrode and water. Ghost atoms will be added below
1073 the water layer.
1074 False: The opposite of True. [default]
1075 int: Explicitly account for the given number of water molecules above
1076 electrode. This is handy if H2O is directly adsorbed and a water layer
1077 is present in the unit cell at the same time.
1078 'plane': Use a plane instead of ghost atoms for freeing the surface.
1079 unsolv_backside: bool
1080 Exclude implicit solvent from the region behind the electrode
1082 """
1084 depends_on_el_density = False
1085 depends_on_atomic_positions = True
1087 def __init__(self, atomic_radii=None, u0=0.180, pbc_cutoff=1e-6,
1088 tiny=1e-10, H2O_layer=False,
1089 unsolv_backside=True, communicator=gpaw.mpi.world):
1090 super().__init__(atomic_radii, u0, pbc_cutoff, tiny)
1091 self.H2O_layer = H2O_layer
1092 self.unsolv_backside = unsolv_backside
1093 self.communicator = communicator
1095 def todict(self):
1096 return {
1097 **super().todict(),
1098 'H2O_layer': self.H2O_layer,
1099 'unsolv_backside': self.unsolv_backside}
1101 def __str__(self):
1102 s = Power12Potential.__str__(self)
1103 s += indent(f' H2O layer: {self.H2O_layer}\n')
1104 s += indent(f' Only solvate front side: {self.unsolv_backside}\n')
1105 return s
1107 def write(self, writer):
1108 writer.write(
1109 name='SJMPower12Potential',
1110 u0=self.u0,
1111 atomic_radii=self.atomic_radii_output,
1112 H2O_layer=self.H2O_layer,
1113 unsolv_backside=self.unsolv_backside)
1115 def update(self, atoms, density):
1116 if atoms is None:
1117 return False
1118 self.r12_a = (self.atomic_radii_output / Bohr) ** 12
1119 r_cutoff = (self.r12_a.max() * self.u0 / self.pbc_cutoff) ** (1. / 12.)
1120 self.pos_aav = get_pbc_positions(atoms, r_cutoff)
1121 self.u_g.fill(.0)
1122 self.grad_u_vg.fill(.0)
1123 na = np.newaxis
1125 if self.unsolv_backside:
1126 # Removing solvent from electrode backside
1127 for z in range(self.u_g.shape[2]):
1128 if (self.r_vg[2, 0, 0, z] - atoms.positions[:, 2].min() /
1129 Bohr < 0):
1130 self.u_g[:, :, z] = np.inf
1131 self.grad_u_vg[:, :, :, z] = 0
1133 if self.H2O_layer:
1134 # Add ghost coordinates and indices to pos_aav dictionary if
1135 # a water layer is present.
1137 all_oxygen_ind = [atom.index for atom in atoms
1138 if atom.symbol == 'O']
1140 # Disregard oxygens that don't belong to the water layer
1141 allwater_oxygen_ind = []
1142 for ox in all_oxygen_ind:
1143 nH = 0
1145 for i, atm in enumerate(atoms):
1146 for period_atm in self.pos_aav[i]:
1147 dist = period_atm * Bohr - atoms[ox].position
1148 if np.linalg.norm(dist) < 1.3 and atm.symbol == 'H':
1149 nH += 1
1151 if nH >= 2:
1152 allwater_oxygen_ind.append(ox)
1154 # If the number of waters in the water layer is given as an input
1155 # (H2O_layer=i) then only the uppermost i water molecules are
1156 # regarded for unsolvating the interface (this is relevant if
1157 # water is adsorbed on the surface)
1158 if not isinstance(self.H2O_layer, (bool, str)):
1159 if self.H2O_layer % 1 < self.tiny:
1160 self.H2O_layer = int(self.H2O_layer)
1161 else:
1162 raise InputError('Only an integer number of water '
1163 'molecules is possible in the water '
1164 'layer')
1166 allwaters = atoms[allwater_oxygen_ind]
1167 indizes_water_ox_ind = np.argsort(allwaters.positions[:, 2],
1168 axis=0)
1170 water_oxygen_ind = []
1171 for i in range(self.H2O_layer):
1172 water_oxygen_ind.append(
1173 allwater_oxygen_ind[indizes_water_ox_ind[-1 - i]])
1175 else:
1176 water_oxygen_ind = allwater_oxygen_ind
1178 oxygen = self.pos_aav[water_oxygen_ind[0]] * Bohr
1179 if len(water_oxygen_ind) > 1:
1180 for windex in water_oxygen_ind[1:]:
1181 oxygen = np.concatenate(
1182 (oxygen, self.pos_aav[windex] * Bohr))
1184 O_layer = []
1185 if isinstance(self.H2O_layer, str):
1186 # Add a virtual plane
1187 if len(self.H2O_layer.split('-')) > 1:
1188 plane_z = float(self.H2O_layer.split('-')[1]) - \
1189 1.0 * self.atomic_radii_output[water_oxygen_ind[0]]
1190 else:
1191 plane_rel_oxygen = -1.5 * self.atomic_radii_output[
1192 water_oxygen_ind[0]]
1193 plane_z = oxygen[:, 2].min() + plane_rel_oxygen
1195 r_diff_zg = self.r_vg[2, :, :, :] - plane_z / Bohr
1196 r_diff_zg[r_diff_zg < self.tiny] = self.tiny
1197 r_diff_zg2 = r_diff_zg ** 2
1198 u_g = self.r12_a[water_oxygen_ind[0]] / r_diff_zg2 ** 6
1199 self.u_g += u_g.copy()
1200 u_g /= r_diff_zg2
1201 r_diff_zg *= u_g.copy()
1202 self.grad_u_vg[2, :, :, :] += r_diff_zg
1204 else:
1205 # Ghost atoms are added below the explicit water layer
1206 cell = atoms.cell.copy() / Bohr
1207 cell[2][2] = 1.
1208 natoms_in_plane = [round(np.linalg.norm(cell[0]) * 1.5),
1209 round(np.linalg.norm(cell[1]) * 1.5)]
1211 plane_z = (oxygen[:, 2].min() - 1.75 *
1212 self.atomic_radii_output[water_oxygen_ind[0]])
1213 nghatoms_z = int(round(oxygen[:, 2].min() -
1214 atoms.positions[:, 2].min()))
1216 for i in range(int(natoms_in_plane[0])):
1217 for j in range(int(natoms_in_plane[1])):
1218 for k in np.linspace(atoms.positions[:, 2].min(),
1219 plane_z, num=nghatoms_z):
1221 O_layer.append(np.dot(np.array(
1222 [(1.5 * i - natoms_in_plane[0] / 4) /
1223 natoms_in_plane[0],
1224 (1.5 * j - natoms_in_plane[1] / 4) /
1225 natoms_in_plane[1],
1226 k / Bohr]), cell))
1228 # Add additional ghost O-atoms below the actual water O atoms
1229 # of water which frees the interface in case of corrugated
1230 # water layers
1231 for ox in oxygen / Bohr:
1232 O_layer.append([ox[0], ox[1], ox[2] - 1.0 *
1233 self.atomic_radii_output[
1234 water_oxygen_ind[0]] / Bohr])
1236 r12_add = []
1237 for i in range(len(O_layer)):
1238 self.pos_aav[len(atoms) + i] = [O_layer[i]]
1239 r12_add.append(self.r12_a[water_oxygen_ind[0]])
1240 r12_add = np.array(r12_add)
1241 # r12_a must have same dimensions as pos_aav items
1242 self.r12_a = np.concatenate((self.r12_a, r12_add))
1244 for index, pos_av in self.pos_aav.items():
1245 pos_av = np.array(pos_av)
1246 r12 = self.r12_a[index]
1247 for pos_v in pos_av:
1248 origin_vg = pos_v[:, na, na, na]
1249 r_diff_vg = self.r_vg - origin_vg
1250 r_diff2_g = (r_diff_vg ** 2).sum(0)
1251 r_diff2_g[r_diff2_g < self.tiny] = self.tiny
1252 u_g = r12 / r_diff2_g ** 6
1253 self.u_g += u_g
1254 u_g /= r_diff2_g
1255 r_diff_vg *= u_g[na, ...]
1256 self.grad_u_vg += r_diff_vg
1258 self.u_g *= self.u0 / Ha
1259 self.grad_u_vg *= -12. * self.u0 / Ha
1260 self.grad_u_vg[self.grad_u_vg < -1e20] = -1e20
1261 self.grad_u_vg[self.grad_u_vg > 1e20] = 1e20
1263 return True
1266class SJM_RealSpaceHamiltonian(SolvationRealSpaceHamiltonian):
1267 """Realspace Hamiltonian with continuum solvent model in the context of
1268 SJM.
1270 See also Section III of
1271 A. Held and M. Walter, J. Chem. Phys. 141, 174108 (2014).
1273 In contrast to the standard implicit solvent model a dipole correction can
1274 also be applied; this is the only difference from its parent.
1275 """
1277 def __init__(self, cavity, dielectric, interactions, gd, finegd, nspins,
1278 setups, timer, xc, world, redistributor, vext=None,
1279 psolver=None, stencil=3, collinear=None, dirichlet=False):
1281 self.cavity = cavity
1282 self.dielectric = dielectric
1283 self.interactions = interactions
1284 cavity.set_grid_descriptor(finegd)
1285 dielectric.set_grid_descriptor(finegd)
1286 for ia in interactions:
1287 ia.set_grid_descriptor(finegd)
1289 if psolver is None:
1290 psolver = WeightedFDPoissonSolver()
1291 self.dipcorr = False
1292 elif isinstance(psolver, dict):
1293 # Sadly the 'dipolelayer' cannot be pop'd because
1294 # CavityShapedJellium calls this twice.
1295 poi_par = {a: b for a, b in psolver.items() if a != 'dipolelayer'}
1296 psolver = SJMDipoleCorrection(WeightedFDPoissonSolver(**poi_par),
1297 psolver['dipolelayer'],
1298 dirichlet=dirichlet)
1299 self.dipcorr = True
1301 if self.dipcorr:
1302 psolver.poissonsolver.set_dielectric(self.dielectric)
1303 else:
1304 psolver.set_dielectric(self.dielectric)
1306 self.gradient = None
1308 RealSpaceHamiltonian.__init__(
1309 self,
1310 gd, finegd, nspins, collinear, setups, timer, xc, world,
1311 vext=vext, psolver=psolver,
1312 stencil=stencil, redistributor=redistributor)
1314 for ia in interactions:
1315 setattr(self, 'e_' + ia.subscript, None)
1316 self.new_atoms = None
1317 self.vt_ia_g = None
1318 self.e_total_free = None
1319 self.e_total_extrapolated = None
1321 def initialize(self):
1322 if self.dipcorr:
1323 self.gradient = [Gradient(self.finegd, i, 1.0,
1324 self.poisson.poissonsolver.nn)
1325 for i in (0, 1, 2)]
1326 else:
1327 self.gradient = [Gradient(self.finegd, i, 1.0,
1328 self.poisson.nn)
1329 for i in (0, 1, 2)]
1331 self.vt_ia_g = self.finegd.zeros()
1332 self.cavity.allocate()
1333 self.dielectric.allocate()
1334 for ia in self.interactions:
1335 ia.allocate()
1336 RealSpaceHamiltonian.initialize(self)
1339class CavityShapedJellium(Jellium):
1340 """The jellium object, where the counter charge takes the form of the
1341 solvent cavity. It puts the jellium background charge where the solvent is
1342 present and z < z2.
1344 Parameters:
1345 ----------
1346 charge: float
1347 The total jellium background charge.
1348 g_g: array
1349 The g function from the implicit solvent model, representing the
1350 percentage of the actual dielectric constant on the grid.
1351 z2: float
1352 Position of upper surface in Angstrom units.
1353 """
1355 def __init__(self, charge, g_g, z2):
1357 Jellium.__init__(self, charge)
1358 self.g_g = g_g
1359 self.z2 = (z2 - 0.0001) / Bohr
1361 def todict(self):
1362 dct = Jellium.todict(self)
1363 dct.update(z2=self.z2 * Bohr + 0.0001,
1364 z1='cavity_like')
1365 return dct
1367 def get_mask(self):
1368 r_gv = self.gd.get_grid_point_coordinates().transpose((1, 2, 3, 0))
1369 mask = np.logical_not(r_gv[:, :, :, 2] > self.z2).astype(float)
1370 mask *= self.g_g
1371 return mask
1374class SJMDipoleCorrection(DipoleCorrection):
1375 """Dipole-correcting wrapper around another PoissonSolver specific for SJM.
1377 Iterative dipole correction class as applied in SJM.
1379 Notes
1380 -----
1382 The modules can easily be incorporated in the trunk version of GPAW
1383 by just adding the `fd_solv_solve` and adapting the `solve` modules
1384 in the `DipoleCorrection` class.
1386 This module is currently calculating the correcting dipole potential
1387 iteratively and we would be very grateful if anybody could
1388 provide an analytical solution.
1390 New Parameters
1391 ---------
1392 corrterm: float
1393 Correction factor for the added countering dipole. This is calculated
1394 iteratively.
1396 last_corrterm: float
1397 Corrterm in the last iteration for getting the change of slope with
1398 change in corrterm
1400 last_slope: float
1401 Same as for `last_corrterm`
1403 """
1404 def __init__(self, poissonsolver, direction, width=1.0, dirichlet=False):
1405 """Construct dipole correction object."""
1407 DipoleCorrection.__init__(self, poissonsolver, direction, width=1.0)
1408 self.corrterm = 1
1409 self.elcorr = None
1410 self.last_corrterm = None
1411 self.dirichlet = dirichlet
1413 def solve(self, pot, dens, **kwargs):
1414 if isinstance(dens, np.ndarray):
1415 # finite-diference Poisson solver:
1416 if hasattr(self.poissonsolver, 'dielectric'):
1417 return self.fd_solv_solve(pot, dens, **kwargs)
1418 else:
1419 return self.fdsolve(pot, dens, **kwargs)
1420 # Plane-wave solver:
1421 self.pwsolve(pot, dens)
1423 def fd_solv_solve(self, vHt_g, rhot_g, **kwargs):
1425 gd = self.poissonsolver.gd
1426 slope_lim = 1e-8
1427 slope = slope_lim * 10
1429 dipmom = gd.calculate_dipole_moment(rhot_g)[2]
1431 if self.elcorr is not None:
1432 vHt_g[:, :] -= self.elcorr
1434 iters2 = self.poissonsolver.solve(vHt_g, rhot_g, **kwargs)
1435 sawtooth_z = self.sjm_sawtooth(dirichlet=self.dirichlet)
1436 L = gd.cell_cv[2, 2]
1438 while abs(slope) > slope_lim:
1439 vHt_g2 = vHt_g.copy()
1440 self.correction = 2 * np.pi * dipmom * L / \
1441 gd.volume * self.corrterm
1442 elcorr = -2 * self.correction
1444 elcorr *= sawtooth_z
1445 elcorr2 = elcorr[gd.beg_c[2]:gd.end_c[2]]
1446 vHt_g2[:, :] += elcorr2
1448 VHt_g = gd.collect(vHt_g2, broadcast=True)
1449 VHt_z = VHt_g.mean(0).mean(0)
1450 slope = VHt_z[2] - VHt_z[10]
1452 if abs(slope) > slope_lim:
1453 if self.last_corrterm is not None:
1454 ds = (slope - self.last_slope) / \
1455 (self.corrterm - self.last_corrterm)
1456 con = slope - (ds * self.corrterm)
1457 self.last_corrterm = self.corrterm
1458 self.corrterm = -con / ds
1459 else:
1460 self.last_corrterm = self.corrterm
1461 self.corrterm -= slope * 10.
1462 self.last_slope = slope
1463 else:
1464 vHt_g[:, :] += elcorr2
1465 self.elcorr = elcorr2
1467 return iters2
1469 def sjm_sawtooth(self, dirichlet):
1470 """Creates a linear function normalized between -0.5 and 0.5 whose
1471 slope is scaled based on the xy-averaged dielectric constant vs z"""
1472 gd = self.poissonsolver.gd
1473 c = self.c
1474 L = gd.cell_cv[c, c]
1475 step = gd.h_cv[c, c] / L
1477 eps_g = gd.collect(self.poissonsolver.dielectric.eps_gradeps[0],
1478 broadcast=True)
1479 eps_z = eps_g.mean(0).mean(0)
1481 saw = np.zeros((int(L / gd.h_cv[c, c])))
1482 for i, eps in enumerate(eps_z):
1483 saw[i + 1] = saw[i] + step / eps
1484 saw /= saw[-1] + step / eps_z[-1] - saw[0]
1486 if dirichlet:
1487 saw -= saw[-1]
1488 else:
1489 saw -= (saw[0] + saw[-1] + step / eps_z[-1]) / 2.
1491 return saw
1494class PotentialConvergenceError(ConvergenceError):
1495 """Raised if potential did not equilibrate."""