Coverage for gpaw/test/generic/test_proton.py: 100%

1"""Test calculation for bare proton. 

2 

3Also, its interaction with an external potential in the form of a point charge 

4is tested. 

5""" 

6import numpy as np 

7import pytest 

8from ase import Atoms 

9from ase.units import Bohr, Hartree 

10 

11from gpaw import GPAW 

12from gpaw.external import PointChargePotential 

13 

14 

15@pytest.mark.old_gpaw_only 

16def test_generic_proton(in_tmp_dir): 

17 a = 4.5 

18 H = Atoms('H', [(a / 2, a / 2, a / 2)], 

19 pbc=0, 

20 cell=(a, a, a)) 

21 H.calc = GPAW(mode='fd', nbands=1, h=0.2, charge=1, txt='H.txt') 

22 e0 = H.get_potential_energy() 

23 assert abs(e0 + H.calc.get_reference_energy()) < 0.014 

24 

25 # Test the point charge potential with a smooth cutoff: 

26 pcp = PointChargePotential([-1], rc2=5.5, width=1.5) 

27 H.calc = H.calc.new(external=pcp) 

28 E = [] 

29 F = [] 

30 D = np.linspace(2, 6, 30) 

31 for d in D: 

32 pcp.set_positions([[a / 2, a / 2, a / 2 + d]]) 

33 e = H.get_potential_energy() - e0 

34 f1 = H.get_forces() 

35 f2 = pcp.get_forces(H.calc) 

36 eref = -1 / d * Bohr * Hartree 

37 print(d, e, eref, abs(f1 + f2).max()) 

38 if d < 4.0: 

39 error = e + 1 / d * Bohr * Hartree 

40 assert abs(error) < 0.01, error 

41 assert abs(f1 + f2).max() < 0.01 

42 E.append(e) 

43 F.append(f1[0, 2]) 

44 

45 E = np.array(E) 

46 FF = (E[2:] - E[:-2]) / (D[2] - D[0]) # numerical force 

47 print(abs(F[1:-1] - FF).max()) 

48 assert abs(F[1:-1] - FF).max() < 0.1 

49 

50 if not True: 

51 import matplotlib.pyplot as plt 

52 plt.plot(D, E) 

53 plt.plot(D, -1 / D * Bohr * Hartree) 

54 plt.plot(D, F) 

55 plt.plot(D[1:-1], FF) 

56 plt.show()