Coverage for gpaw/test/ext_potential/test_point_charge.py: 100%

1import numpy as np 

2import pytest 

3from ase import Atoms 

4from gpaw.test import calculate_numerical_forces 

5 

6import gpaw.cgpaw as cgpaw 

7from gpaw import GPAW, Mixer, PoissonSolver 

8from gpaw.external import PointChargePotential 

9 

10 

11@pytest.mark.old_gpaw_only 

12def test_ext_potential_point_charge(in_tmp_dir): 

13 # Find coefs for polynomial: 

14 c = np.linalg.solve([[1, 1, 1, 1], 

15 [0, 2, 4, 6], 

16 [0, 2, 12, 30], 

17 [1 / 3, 1 / 5, 1 / 7, 1 / 9]], 

18 [1, -1, 2, 0.5]) 

19 print(c) 

20 print(c * 32) 

21 x = np.linspace(0, 1, 101) 

22 v = np.polyval(c[::-1], x**2) 

23 assert (v * x**2 - x).sum() / 100 == pytest.approx(0, abs=1e-5) 

24 

25 if 0: 

26 import matplotlib.pyplot as plt 

27 plt.plot(x, v) 

28 x = np.linspace(0.2, 1.5, 101) 

29 plt.plot(x, 1 / x) 

30 plt.show() 

31 

32 # Low-level test: 

33 h = 0.1 

34 q = 2.2 

35 beg_v = np.zeros(3, int) 

36 h_v = np.ones(3) * h 

37 q_p = np.ones(1) * q 

38 R_pv = np.array([[0.1, 0.2, -0.3]]) 

39 vext_G = np.zeros((1, 1, 1)) 

40 rhot_G = np.ones((1, 1, 1)) 

41 

42 def f(rc): 

43 vext_G[:] = 0.0 

44 cgpaw.pc_potential(beg_v, h_v, q_p, R_pv, rc, 

45 np.inf, 1.0, vext_G, None) 

46 return -vext_G[0, 0, 0] 

47 

48 d = (R_pv[0]**2).sum()**0.5 

49 

50 for rc in [0.9 * d, 1.1 * d]: 

51 if d > rc: 

52 v0 = q / d 

53 else: 

54 v0 = np.polyval(c[::-1], (d / rc)**2) * q / rc 

55 

56 assert f(rc) == pytest.approx(v0, abs=1e-12) 

57 

58 F_pv = np.zeros((1, 3)) 

59 cgpaw.pc_potential(beg_v, h_v, q_p, R_pv, rc, np.inf, 1.0, 

60 vext_G, None, rhot_G, F_pv) 

61 eps = 0.0001 

62 for v in range(3): 

63 R_pv[0, v] += eps 

64 ep = f(rc) 

65 R_pv[0, v] -= 2 * eps 

66 em = f(rc) 

67 R_pv[0, v] += eps 

68 F = -(ep - em) / (2 * eps) * h**3 

69 assert F == pytest.approx(-F_pv[0, v], abs=1e-9) 

70 

71 # High-level test: 

72 lih = Atoms('LiH') 

73 lih[1].y += 1.64 

74 lih.center(vacuum=3) 

75 

76 pos = lih.cell.sum(axis=0) 

77 print(pos) 

78 pc = PointChargePotential([-1.0], [pos]) 

79 lih.calc = GPAW(mode='fd', 

80 external=pc, 

81 mixer=Mixer(0.8, 5, 20.0), 

82 xc='oldLDA', 

83 poissonsolver=PoissonSolver(), 

84 convergence=dict(density=3e-6), 

85 txt='lih-pc.txt') 

86 f1 = lih.get_forces() 

87 fpc1 = pc.get_forces(lih.calc) 

88 print(fpc1) 

89 print(fpc1 + f1.sum(0)) 

90 

91 f2 = calculate_numerical_forces(lih, 0.001) 

92 print(f1) 

93 print(f1 - f2) 

94 err = abs(f1 - f2).max() 

95 print('err', err) 

96 assert err < 2e-3, err 

97 

98 x = 0.0001 

99 for v in range(3): 

100 pos[v] += x 

101 pc.set_positions([pos]) 

102 ep = lih.get_potential_energy() 

103 pos[v] -= 2 * x 

104 pc.set_positions([pos]) 

105 em = lih.get_potential_energy() 

106 pos[v] += x 

107 pc.set_positions([pos]) 

108 error = (em - ep) / (2 * x) - fpc1[0, v] 

109 print(v, error) 

110 assert abs(error) < 0.006