Coverage for gpaw/solvation/interactions.py: 83%

1import numpy as np 

2from ase.units import Bohr, Hartree 

3 

4from gpaw.solvation.gridmem import NeedsGD 

5 

6 

7class Interaction(NeedsGD): 

8 """Base class for non-electrostatic solvent solute interactions. 

9 

10 See also Section III of 

11 A. Held and M. Walter, J. Chem. Phys. 141, 174108 (2014). 

12 

13 Attributes: 

14 subscript -- Short label used to reference the interaction. 

15 E -- Local interaction energy in Hartree 

16 delta_E_delta_n_g -- Functional derivative of the interaction energy 

17 with respect to the electron density. 

18 delta_E_delta_n_g -- Functional derivative of the interaction energy 

19 with respect to the cavity. 

20 """ 

21 

22 subscript = 'unnamed' 

23 

24 def __init__(self): 

25 NeedsGD.__init__(self) 

26 self.E = None 

27 self.delta_E_delta_n_g = None 

28 self.delta_E_delta_g_g = None 

29 

30 @classmethod 

31 def from_dict(self, dct): 

32 if not isinstance(dct, dict): 

33 return dct 

34 dct = dct.copy() 

35 name = dct.pop('name') 

36 if name == 'SurfaceInteraction': 

37 return SurfaceInteraction(**dct) 

38 if name == 'VolumeInteraction': 

39 return VolumeInteraction(**dct) 

40 if name == 'LeakedDensityInteraction': 

41 return LeakedDensityInteraction(**dct) 

42 raise ValueError(name) 

43 

44 def write(self, writer): 

45 pass 

46 

47 def update(self, atoms, density, cavity): 

48 """Update the Kohn-Sham potential and the energy. 

49 

50 atoms and/or cavity are None iff they have not changed 

51 since the last call 

52 

53 Return whether the interaction has changed. 

54 """ 

55 raise NotImplementedError 

56 

57 def estimate_memory(self, mem): 

58 ngrids = 1 + self.depends_on_el_density 

59 mem.subnode('Functional Derivatives', ngrids * self.gd.bytecount()) 

60 

61 def allocate(self): 

62 NeedsGD.allocate(self) 

63 self.delta_E_delta_g_g = self.gd.empty() 

64 if self.depends_on_el_density: 

65 self.delta_E_delta_n_g = self.gd.empty() 

66 

67 @property 

68 def depends_on_atomic_positions(self): 

69 """Return whether the ia depends explicitly on atomic positions.""" 

70 raise NotImplementedError 

71 

72 @property 

73 def depends_on_el_density(self): 

74 """Return whether the ia depends explicitly on the electron density.""" 

75 raise NotImplementedError 

76 

77 def get_del_r_vg(self, atomic_index, density): 

78 """Return spatial derivatives with respect to atomic position.""" 

79 raise NotImplementedError 

80 

81 def __str__(self): 

82 s = f"Interaction: {self.__class__.__name__}\n" 

83 s += f" subscript: {self.subscript}\n" 

84 return s 

85 

86 def update_atoms(self, atoms, log): 

87 """Inexpensive initialization when atoms change.""" 

88 pass 

89 

90 

91class SurfaceInteraction(Interaction): 

92 """An interaction with energy proportional to the cavity surface area.""" 

93 

94 subscript = 'surf' 

95 depends_on_el_density = False 

96 depends_on_atomic_positions = False 

97 

98 def __init__(self, surface_tension): 

99 """Constructor for SurfaceInteraction class. 

100 

101 Arguments: 

102 surface_tension -- Proportionality factor to calculate 

103 energy from surface area in eV / Angstrom ** 2. 

104 """ 

105 Interaction.__init__(self) 

106 self.surface_tension = float(surface_tension) 

107 

108 def todict(self): 

109 return {'surface_tension': self.surface_tension} 

110 

111 def write(self, writer): 

112 writer.write(name='SurfaceInteraction', 

113 surface_tension=self.surface_tension) 

114 

115 def update(self, atoms, density, cavity): 

116 if cavity is None: 

117 return False 

118 acalc = cavity.surface_calculator 

119 st = self.surface_tension * Bohr ** 2 / Hartree 

120 self.E = st * acalc.A 

121 np.multiply(st, acalc.delta_A_delta_g_g, self.delta_E_delta_g_g) 

122 return True 

123 

124 def __str__(self): 

125 s = Interaction.__str__(self) 

126 s += f' Surface_tension: {self.surface_tension} eV/(Angstrom^2)' 

127 return s 

128 

129 

130class VolumeInteraction(Interaction): 

131 """An interaction with energy proportional to the cavity volume""" 

132 

133 subscript = 'vol' 

134 depends_on_el_density = False 

135 depends_on_atomic_positions = False 

136 

137 def __init__(self, pressure): 

138 """Constructor for VolumeInteraction class. 

139 

140 Arguments: 

141 pressure -- Proportionality factor to calculate 

142 energy from volume in eV / Angstrom ** 3. 

143 """ 

144 Interaction.__init__(self) 

145 self.pressure = float(pressure) 

146 

147 def todict(self): 

148 return {'pressure': self.pressure} 

149 

150 def update(self, atoms, density, cavity): 

151 if cavity is None: 

152 return False 

153 vcalc = cavity.volume_calculator 

154 pressure = self.pressure * Bohr ** 3 / Hartree 

155 self.E = pressure * vcalc.V 

156 np.multiply(pressure, vcalc.delta_V_delta_g_g, self.delta_E_delta_g_g) 

157 return True 

158 

159 def __str__(self): 

160 s = Interaction.__str__(self) 

161 s += f' Pressure: {self.pressure}' 

162 return s 

163 

164 

165class LeakedDensityInteraction(Interaction): 

166 """Interaction proportional to charge leaking outside cavity. 

167 

168 The charge outside the cavity is calculated as 

169 """ 

170 

171 subscript = 'leak' 

172 depends_on_el_density = True 

173 depends_on_atomic_positions = False 

174 

175 def __init__(self, voltage): 

176 """Constructor for LeakedDensityInteraction class. 

177 

178 Arguments: 

179 voltage -- Proportionality factor to calculate 

180 energy from integrated electron density in Volts. 

181 A positive value of the voltage leads to a 

182 positive interaction energy. 

183 """ 

184 Interaction.__init__(self) 

185 self.voltage = float(voltage) 

186 

187 def todict(self): 

188 return {'voltage': self.voltage} 

189 

190 def update(self, atoms, density, cavity): 

191 E0 = self.voltage / Hartree 

192 if cavity is not None: 

193 np.multiply(E0, cavity.g_g, self.delta_E_delta_n_g) 

194 np.multiply(E0, density.nt_g, self.delta_E_delta_g_g) 

195 self.E = self.gd.integrate( 

196 density.nt_g * self.delta_E_delta_n_g, 

197 global_integral=False 

198 ) 

199 return True 

200 

201 def __str__(self): 

202 s = Interaction.__str__(self) 

203 s += f' Voltage: {self.voltage:.4f}' 

204 return s