Coverage for gpaw/solvation/calculator.py: 86%
88 statements
« prev ^ index » next coverage.py v7.7.1, created at 2025-07-14 00:18 +0000
1from ase.data import chemical_symbols
2from ase.units import Bohr, Hartree
3from gpaw import GPAW_NEW
4from gpaw.calculator import GPAW as OldGPAW
5from gpaw.io import Reader
6from gpaw.solvation.hamiltonian import SolvationRealSpaceHamiltonian
9def SolvationGPAW(*args, **kwargs):
10 if GPAW_NEW:
11 from gpaw.new.ase_interface import GPAW
12 solvation = dict(name='solvation',
13 cavity=kwargs.pop('cavity'),
14 dielectric=kwargs.pop('dielectric'),
15 interactions=kwargs.pop('interactions', None))
16 return GPAW(*args, **kwargs, environment=solvation)
17 return OldSolvationGPAW(*args, **kwargs)
20class OldSolvationGPAW(OldGPAW):
21 """Subclass of gpaw.GPAW calculator with continuum solvent model.
23 See also Section III of
24 A. Held and M. Walter, J. Chem. Phys. 141, 174108 (2014).
25 """
27 def __init__(self, restart=None, cavity=None, dielectric=None,
28 interactions=None, **gpaw_kwargs):
29 """Constructor for SolvationGPAW class.
31 Additional arguments not present in GPAW class:
32 cavity -- A Cavity instance.
33 dielectric -- A Dielectric instance.
34 interactions -- A list of Interaction instances.
35 """
36 if interactions is None:
37 interactions = []
39 # if not all([cavity, dielectric]):
40 # raise IOError('Cavity and dielectric modules need to be '
41 # 'defined in the calculator')
43 self.stuff_for_hamiltonian = (cavity, dielectric, interactions)
45 OldGPAW.__init__(self, restart, **gpaw_kwargs)
47 self.log('Implicit solvation parameters:')
48 for stuff in self.stuff_for_hamiltonian:
49 if isinstance(stuff, list):
50 for instuff in stuff:
51 self.log(instuff)
52 else:
53 self.log(stuff)
54 self.log()
56 def read(self, filename):
57 """Read yourself from a file"""
58 self.reader = reader = Reader(filename)
59 if 'implicit_solvent' in reader:
60 impl_in = reader.implicit_solvent
61 if 'name' in impl_in.cavity.effective_potential:
62 efpot = impl_in.cavity.effective_potential
64 atomic_radii = {}
65 for Z, r in zip(reader.atoms.numbers, efpot.atomic_radii):
66 symbol = chemical_symbols[Z]
67 if symbol in atomic_radii:
68 assert atomic_radii[symbol] == r
69 else:
70 atomic_radii[symbol] = r
72 if efpot.name == 'SJMPower12Potential':
73 from gpaw.solvation.sjm import SJMPower12Potential
75 effective_potential = SJMPower12Potential(
76 atomic_radii=atomic_radii,
77 u0=efpot.u0,
78 H2O_layer=efpot.H2O_layer,
79 unsolv_backside=efpot.unsolv_backside)
80 elif efpot.name == 'Power12Potential':
81 from gpaw.solvation.sjm import Power12Potential
82 effective_potential = Power12Potential(
83 atomic_radii=atomic_radii,
84 u0=efpot.u0)
85 else:
86 raise OSError('Reading the given effective potential '
87 'is not implemented yet')
89 if 'name' in impl_in.cavity.surface_calculator:
90 suca = impl_in.cavity.surface_calculator
91 if suca.name == 'GradientSurface':
92 from gpaw.solvation.cavity import GradientSurface
93 surface_calculator = GradientSurface(suca.nn)
94 else:
95 raise OSError('Reading in the given used surface '
96 'calculator is not implemented')
98 T = impl_in.cavity.temperature
100 from gpaw.solvation.cavity import EffectivePotentialCavity
101 cavity = EffectivePotentialCavity(
102 effective_potential=effective_potential,
103 temperature=T,
104 surface_calculator=surface_calculator)
106 if impl_in.dielectric.name == 'LinearDielectric':
107 from gpaw.solvation.dielectric import LinearDielectric
108 dielectric = LinearDielectric(epsinf=impl_in.dielectric.epsinf)
110 if impl_in.interactions.name == 'SurfaceInteraction':
111 suin = impl_in.interactions
112 from gpaw.solvation.interactions import SurfaceInteraction
113 interactions = [SurfaceInteraction(suin.surface_tension)]
115 self.stuff_for_hamiltonian = (cavity, dielectric, interactions)
117 reader = OldGPAW.read(self, filename)
118 return reader
120 def _write(self, writer, mode):
121 OldGPAW._write(self, writer, mode)
122 stuff = self.stuff_for_hamiltonian
123 writer.child('implicit_solvent').write(cavity=stuff[0],
124 dielectric=stuff[1],
125 interactions=stuff[2][0])
127 def create_hamiltonian(self, realspace, mode, xc):
128 if not realspace:
129 raise NotImplementedError(
130 'SolvationGPAW does not support '
131 'calculations in reciprocal space yet.')
133 dens = self.density
134 self.hamiltonian = SolvationRealSpaceHamiltonian(
135 *self.stuff_for_hamiltonian,
136 gd=dens.gd, finegd=dens.finegd,
137 nspins=dens.nspins,
138 collinear=dens.collinear,
139 setups=dens.setups,
140 timer=self.timer,
141 xc=xc,
142 world=self.world,
143 redistributor=dens.redistributor,
144 vext=self.parameters.external,
145 psolver=self.parameters.poissonsolver,
146 stencil=mode.interpolation)
148 self.log(self.hamiltonian)
150 xc.set_grid_descriptor(self.hamiltonian.finegd)
152 def initialize_positions(self, atoms=None):
153 spos_ac = OldGPAW.initialize_positions(self, atoms)
154 self.hamiltonian.update_atoms(self.atoms, self.log)
155 return spos_ac
157 def get_electrostatic_energy(self):
158 """Return electrostatic part of the total energy.
160 The electrostatic part consists of everything except
161 the short-range interactions defined in the interactions list.
162 """
163 # Energy extrapolated to zero width:
164 return Hartree * self.hamiltonian.e_el_extrapolated
166 def get_solvation_interaction_energy(self, subscript, atoms=None):
167 """Return a specific part of the solvation interaction energy.
169 The subscript parameter defines which part is to be returned.
170 It has to match the value of a subscript attribute of one of
171 the interactions in the interactions list.
172 """
173 return Hartree * getattr(self.hamiltonian, 'e_' + subscript)
175 def get_cavity_volume(self, atoms=None):
176 """Return the cavity volume in Angstrom ** 3.
178 In case no volume calculator has been set for the cavity, None
179 is returned.
180 """
181 V = self.hamiltonian.cavity.V
182 return V and V * Bohr ** 3
184 def get_cavity_surface(self, atoms=None):
185 """Return the cavity surface area in Angstrom ** 2.
187 In case no surface calculator has been set for the cavity,
188 None is returned.
189 """
190 A = self.hamiltonian.cavity.A
191 return A and A * Bohr ** 2