Coverage for gpaw/test/solvation/test_forces.py: 68%

1# flake8: noqa 

2import pytest 

3from ase import Atoms 

4from ase.units import Pascal, m 

5from ase.data.vdw import vdw_radii 

6from gpaw.mpi import rank 

7from gpaw import Mixer 

8from gpaw.solvation import ( 

9 SolvationGPAW, 

10 EffectivePotentialCavity, 

11 Power12Potential, 

12 LinearDielectric, 

13 KB51Volume, 

14 GradientSurface, 

15 VolumeInteraction, 

16 SurfaceInteraction, 

17 LeakedDensityInteraction) 

18 

19import numpy as np 

20 

21SKIP_ENERGY_CALCULATION = True 

22F_max_err = 0.005 

23 

24h = 0.2 

25u0 = 0.180 

26epsinf = 80. 

27T = 298.15 

28 

29 

30@pytest.mark.slow 

31def test_solvation_forces(): 

32 atoms = Atoms('NaCl', positions=((5.6, 5.6, 6.8), (5.6, 5.6, 8.8))) 

33 atoms.set_cell((11.2, 11.2, 14.4)) 

34 

35 atoms.calc = SolvationGPAW( 

36 mode='fd', 

37 mixer=Mixer(0.5, 7, 50.0), 

38 eigensolver='dav', 

39 xc='oldPBE', h=h, setups={'Na': '1'}, 

40 cavity=EffectivePotentialCavity( 

41 effective_potential=Power12Potential(u0=u0), 

42 temperature=T, 

43 volume_calculator=KB51Volume(), 

44 surface_calculator=GradientSurface()), 

45 dielectric=LinearDielectric(epsinf=epsinf), 

46 # parameters chosen to give ~ 1eV for each interaction 

47 interactions=[ 

48 VolumeInteraction(pressure=-1e9 * Pascal), 

49 SurfaceInteraction(surface_tension=100. * 1e-3 * Pascal * m), 

50 LeakedDensityInteraction(voltage=10.)]) 

51 

52 def vac(atoms): 

53 return min( 

54 atoms.positions[0][2], 

55 14.4 - atoms.positions[1][2]) 

56 

57 step = 0.05 

58 if not SKIP_ENERGY_CALCULATION: 

59 d = [] 

60 E = [] 

61 F = [] 

62 while vac(atoms) >= 5.6: 

63 d.append(abs(atoms.positions[0][2] - atoms.positions[1][2])) 

64 E.append(atoms.calc.get_potential_energy(atoms, force_consistent=True)) 

65 F.append(atoms.calc.get_forces(atoms)) 

66 atoms.positions[0][2] -= step 

67 

68 d = np.array(d) 

69 E = np.array(E) 

70 F = np.array(F) 

71 

72 if rank == 0: 

73 np.save('d.npy', d) 

74 np.save('E.npy', E) 

75 np.save('F.npy', F) 

76 from pprint import pprint 

77 print('d') 

78 pprint(list(d)) 

79 print() 

80 print('E') 

81 pprint(list(E)) 

82 print() 

83 print('F') 

84 pprint([list([list(l2) for l2 in l1]) for l1 in F]) 

85 else: 

86 # h=0.2, setups: 0.9.11271, analytic gradient for dielectric 

87 d = [ 

88 2.0000000000000009, 

89 2.0500000000000007, 

90 2.1000000000000005, 

91 2.1500000000000004, 

92 2.2000000000000002, 

93 2.25, 

94 2.2999999999999998, 

95 2.3499999999999996, 

96 2.3999999999999995, 

97 2.4499999999999993, 

98 2.4999999999999991, 

99 2.5499999999999989, 

100 2.5999999999999988, 

101 2.6499999999999986, 

102 2.6999999999999984, 

103 2.7499999999999982, 

104 2.799999999999998, 

105 2.8499999999999979, 

106 2.8999999999999977, 

107 2.9499999999999975, 

108 2.9999999999999973, 

109 3.0499999999999972, 

110 3.099999999999997, 

111 3.1499999999999968, 

112 3.1999999999999966 

113 ] 

114 

115 E = [ 

116 -3.563293139400153, 

117 -3.8483728941723685, 

118 -4.0711353674406672, 

119 -4.2422702451403147, 

120 -4.370843143863639, 

121 -4.4644180141621437, 

122 -4.5291689227132208, 

123 -4.5703343350581118, 

124 -4.5923649740023134, 

125 -4.5988279473461287, 

126 -4.5927583650683008, 

127 -4.5766206693403326, 

128 -4.552512679729916, 

129 -4.5222435353897978, 

130 -4.4871564101028012, 

131 -4.4484418158342498, 

132 -4.4070815207527305, 

133 -4.3639583908891826, 

134 -4.3196547998536206, 

135 -4.2746896647162842, 

136 -4.2295530542966002, 

137 -4.1846071284157169, 

138 -4.1401332574837273, 

139 -4.0962878505444289, 

140 -4.0533187524468151 

141 ] 

142 

143 F = [ 

144 [[2.6695821634315027e-14, -1.9815222433664118e-14, -6.4048490688641513], 

145 [-3.4082241002543622e-12, -1.6240575950714694e-12, 6.4041543958163452]], 

146 [[5.0237752987734333e-14, -1.8243293058344225e-14, -5.0393351704101388], 

147 [-1.5274896963914418e-12, -2.6398128880388972e-12, 5.0384209277642569]], 

148 [[-3.3648148693779676e-14, -5.9189391352451007e-15, -3.9060170441670885], 

149 [-8.4063321767924475e-13, -1.8821508970704685e-12, 3.905543132953595]], 

150 [[8.033887557888106e-14, 2.2292480376612319e-14, -2.9689893614346223], 

151 [-1.223736631764757e-12, -3.0969447122690189e-12, 2.9698074624482889]], 

152 [[-7.1961182128217342e-14, -9.2369061233814621e-15, -2.1978946263288037], 

153 [-8.8389572960386298e-13, -1.7994648800211978e-12, 2.1973890330684451]], 

154 [[6.839819262929247e-14, -5.8813682508237254e-14, -1.5641714869085164], 

155 [-1.7622743497617748e-12, -3.0046379594525875e-12, 1.5635650982049627]], 

156 [[-9.2964056597828368e-15, -7.5214119452955559e-15, -1.0432778856648208], 

157 [-1.5739852469080523e-12, -2.4988123414408888e-12, 1.0427858515240367]], 

158 [[4.9365637419280049e-14, -1.2806637110762107e-14, -0.61777588826245566], 

159 [-4.4633495248199538e-12, -8.4857800635478483e-13, 0.6191354280714596]], 

160 [[-2.7465001491431585e-15, 2.2194992830262815e-15, -0.27335940186763846], 

161 [2.2159953452415821e-12, -1.6412130688619773e-12, 0.27444763375180131]], 

162 [[-1.1307851878181378e-15, 1.1152866935660596e-14, 0.0049642333143336626], 

163 [-4.8799705828538529e-13, -2.6066991522648012e-12, -0.0056198136590138066]], 

164 [[-3.347924840460358e-14, -3.2847274823028079e-14, 0.22997686507849791], 

165 [-1.3848503589740797e-12, -3.1093262854030672e-12, -0.23005715760038931]], 

166 [[-1.8200707826858258e-14, -4.3512573943996744e-14, 0.40941568333028211], 

167 [1.1291695481871819e-12, -1.7474358989696015e-12, -0.40852810117309907]], 

168 [[-1.8441970643900321e-14, 5.61993561755351e-15, 0.54932435973885352], 

169 [-2.9861492319230874e-12, -3.0227468023234651e-12, -0.54910391544690951]], 

170 [[5.7679940133610544e-14, -2.8084539456163582e-14, 0.65784317706542283], 

171 [-1.0234224221710884e-12, -2.5375102240498517e-12, -0.65870527795968759]], 

172 [[-4.036664664589713e-15, -4.8205738644069875e-15, 0.7418625322893313], 

173 [-1.2904684108015449e-12, -2.1348233844072512e-12, -0.74210985960596609]], 

174 [[1.5318647653685363e-14, -3.421407576047577e-14, 0.80443218412113771], 

175 [8.3926234818054389e-13, -2.6560910148144288e-12, -0.80377843042267316]], 

176 [[6.9157565324685627e-15, -2.0265394116808891e-14, 0.84757769602691035], 

177 [2.2564924910397676e-12, -2.2656092749597723e-12, -0.84788230013436461]], 

178 [[3.6816279559139879e-14, -5.6217408874846806e-14, 0.87630143837759911], 

179 [-2.8919485514440065e-12, -1.5785323374781975e-12, -0.8773089495263976]], 

180 [[-2.7265578221877411e-15, 1.8112097101058187e-14, 0.89457926961794421], 

181 [-4.3811750564077208e-12, -2.461714836405928e-12, -0.89468520749309732]], 

182 [[3.7482219205057949e-14, -4.1131307258530919e-14, 0.90318154899617831], 

183 [-1.0179375774172992e-12, -1.3810014419535958e-12, -0.90288363369205138]], 

184 [[6.8901452284609408e-14, -1.8368129697299518e-14, 0.90212377386324472], 

185 [-2.5725988597478254e-12, -2.4508096630572398e-12, -0.90278569263072728]], 

186 [[-2.9987853255786194e-14, -8.7933798499792121e-15, 0.89490384223117858], 

187 [-1.9656418164895607e-12, -2.497775366538778e-12, -0.89593617789571922]], 

188 [[3.5493492557854517e-14, -9.8132290416331934e-15, 0.88412548228462806], 

189 [-2.1672889098350696e-13, -1.9186968970762135e-12, -0.88434959210747588]], 

190 [[-4.2258218665175417e-14, -7.2197458334904237e-15, 0.86940669722572905], 

191 [-6.3029668473160416e-13, -2.8462694786350656e-12, -0.86883783222557032]], 

192 [[-1.3554138996860724e-14, 7.8037843761066284e-16, 0.84970515380081257], 

193 [-4.5784249812786949e-12, -2.0773781762455308e-12, -0.8501771195912331]] 

194 ] 

195 d = np.array(d) 

196 E = np.array(E) 

197 F = np.array(F) 

198 

199 

200 # test for orthogonal forces equal zero: 

201 assert F[..., :2] == pytest.approx(.0, abs=1e-7) 

202 

203 stencil = 2 # 1 is too rough, 3 does not change compared to 2 

204 FNa, FCl = F[..., 2].T 

205 FNa *= -1. 

206 # test symmetry 

207 assert FNa == pytest.approx(FCl, abs=F_max_err) 

208 dd = np.diff(d)[0] 

209 kernel = { 

210 1: np.array((0.5, 0, -0.5)), 

211 2: np.array((-1. / 12., 2. / 3., 0, -2. / 3., 1. / 12.)), 

212 3: np.array((1. / 60., -0.15, 0.75, 0, -0.75, 0.15, -1. / 60.))} 

213 

214 dEdz = np.convolve(E, kernel[stencil] / step, 'valid') 

215 

216 err = np.maximum( 

217 np.abs(-dEdz - FNa[stencil:-stencil]), 

218 np.abs(-dEdz - FCl[stencil:-stencil])) 

219 

220 # test forces against -dE / dd finite difference 

221 print(err) 

222 assert err == pytest.approx(.0, abs=F_max_err) 

223 

224 if SKIP_ENERGY_CALCULATION: 

225 # check only selected points: 

226 

227 def check(index): 

228 atoms.positions[0][2] = 6.8 - index * step 

229 F_check = atoms.get_forces() 

230 assert F_check[..., :2] == pytest.approx(.0, abs=1e-7) 

231 FNa_check, FCl_check = F_check[..., 2].T 

232 FNa_check *= -1. 

233 assert FNa_check == pytest.approx(FCl_check, abs=F_max_err) 

234 err = np.maximum( 

235 np.abs(-dEdz[index - stencil] - FNa_check), 

236 np.abs(-dEdz[index - stencil] - FCl_check)) 

237 print(err) 

238 assert err == pytest.approx(.0, abs=F_max_err) 

239 

240 l = len(FNa) 

241 # check(stencil) 

242 check(l // 2) 

243 # check(l - 1 - stencil)