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

1import numpy as np 

2import pytest 

3from ase import Atoms 

4from ase.units import Bohr, Ha 

5 

6from gpaw import GPAW 

7from gpaw.external import ExternalPotential, known_potentials 

8from gpaw.poisson import NoInteractionPoissonSolver 

9 

10 

11class HarmonicPotential(ExternalPotential): 

12 def calculate_potential(self, gd): 

13 a = gd.cell_cv[0, 0] 

14 r_vg = gd.get_grid_point_coordinates() 

15 self.vext_g = 0.5 * ((r_vg - a / 2)**2).sum(0) 

16 

17 def todict(self): 

18 return {'name': 'HarmonicPotential'} 

19 

20 

21@pytest.mark.old_gpaw_only 

22def test_ext_potential_harmonic(in_tmp_dir): 

23 """Test against the analytic result (no xc, no Coulomb).""" 

24 a = 4.0 

25 x = Atoms(cell=(a, a, a)) # no atoms 

26 

27 calc = GPAW(mode='fd', 

28 charge=-8, 

29 nbands=4, 

30 h=0.2, 

31 xc={'name': 'null'}, 

32 external=HarmonicPotential(), 

33 poissonsolver=NoInteractionPoissonSolver(), 

34 eigensolver='cg') 

35 

36 x.calc = calc 

37 x.get_potential_energy() 

38 

39 eigs = calc.get_eigenvalues() 

40 assert eigs[0] == pytest.approx(1.5 * Ha, abs=0.002) 

41 assert abs(eigs[1:] - 2.5 * Ha).max() == pytest.approx(0, abs=0.003) 

42 

43 # Check write + read: 

44 calc.write('harmonic.gpw') 

45 known_potentials['HarmonicPotential'] = HarmonicPotential 

46 GPAW('harmonic.gpw') 

47 

48 

49class PöschlTellerPotential(ExternalPotential): 

50 """Slab with Pöschel-Teller well along z-direction. 

51 

52 See: 

53 

54 https://en.wikipedia.org/wiki/P%C3%B6schl%E2%80%93Teller_potential 

55 """ 

56 def calculate_potential(self, gd): 

57 a = gd.cell_cv[2, 2] 

58 r_vg = gd.get_grid_point_coordinates() 

59 lam = 2 

60 self.vext_g = -lam * (lam + 1) / 2 * np.cosh(r_vg[2] - a / 2)**-2 

61 

62 def todict(self): 

63 return {'name': 'HarmonicPotential'} 

64 

65 

66@pytest.mark.old_gpaw_only 

67def test_pt_potential(): 

68 """Test againts analytic result (no xc, no Coulomb).""" 

69 d = 6.0 

70 a = 2 

71 x = Atoms(cell=(a, a, d), pbc=[1, 1, 0]) # no atoms 

72 

73 calc = GPAW(charge=-12, 

74 nbands=8, 

75 mode='pw', 

76 xc={'name': 'null'}, 

77 external=PöschlTellerPotential(), 

78 poissonsolver=NoInteractionPoissonSolver()) 

79 

80 x.calc = calc 

81 x.get_potential_energy() 

82 

83 eigs = calc.get_eigenvalues() / Ha 

84 print(eigs) 

85 k = 2 * np.pi / (a / Bohr) 

86 e0 = -2 

87 e1234 = -2 + 0.5 * k**2 

88 e5 = -0.5 

89 assert eigs[0] == pytest.approx(e0, abs=0.0002) 

90 assert eigs[1] == pytest.approx(e1234, abs=0.001) 

91 assert np.ptp(eigs[1:5]) == pytest.approx(0) 

92 assert eigs[5] == pytest.approx(e5, abs=0.001)