Coverage for gpaw/test/solvation/test_modes.py: 82%

1import numpy as np 

2import pytest 

3from ase import Atoms 

4from gpaw import FermiDirac 

5from gpaw.new.ase_interface import GPAW 

6from gpaw.new.sjm import SJM 

7from gpaw.solvation.sjm import SJM as OldSJM 

8 

9 

10@pytest.mark.parametrize('mode', ['pw', 'fd']) 

11def test_h(gpaw_new, mode, in_tmp_dir): 

12 if mode == 'pw' and not gpaw_new: 

13 pytest.skip('PW-mode not implemented for old GPAW') 

14 

15 a = 1.4 

16 atoms = Atoms('H', cell=[a, a, 11.0], pbc=(1, 1, 0)) 

17 atoms.positions[0, 2] = 4.0 

18 

19 k = 2 

20 params = dict( 

21 mode=mode, 

22 kpts=(k, k, 1), 

23 occupations=FermiDirac(0.2)) 

24 sjm = {'target_potential': 7.5, 

25 'excess_electrons': -0.0045, 

26 'jelliumregion': {'top': 9.0, 'bottom': 7.0}, 

27 'tol': 0.01} 

28 solvation = dict(cavity=NoCavity(), dielectric=Vacuum(), interactions=[]) 

29 

30 if gpaw_new: 

31 atoms.calc = GPAW(environment=SJM(**sjm, **solvation), **params) 

32 atoms.get_potential_energy() 

33 pot = -atoms.calc.get_fermi_level() 

34 else: 

35 atoms.calc = OldSJM(**solvation, sj=sjm, **params) 

36 atoms.get_potential_energy() 

37 pot = atoms.calc.get_electrode_potential() 

38 

39 assert pot == pytest.approx(7.5, abs=0.01) 

40 

41 atoms.write('h.traj') 

42 if gpaw_new: 

43 atoms.calc.write('h.gpw') 

44 

45 def hook(dct): 

46 return SJM(**sjm, **solvation) 

47 

48 GPAW('h.gpw', object_hooks={'environment': hook}) 

49 

50 if 0: 

51 v = atoms.calc.get_electrostatic_potential() 

52 import matplotlib.pyplot as plt 

53 plt.plot(np.linspace(0, 11, v.shape[2], 0), v[0, 0]) 

54 plt.show() 

55 

56 

57class NoCavity: 

58 depends_on_el_density = False 

59 

60 def set_grid_descriptor(self, gd): 

61 pass 

62 

63 def allocate(self): 

64 pass 

65 

66 def update_atoms(self, atoms, log): 

67 pass 

68 

69 def update(self, atoms, density): 

70 pass 

71 

72 def communicate_vol_surf(self, comm): 

73 pass 

74 

75 def summary(self, log): 

76 pass 

77 

78 def todict(self): 

79 return {} 

80 

81 

82class Vacuum: 

83 def set_grid_descriptor(self, gd): 

84 self.eps_gradeps = [gd.zeros() for _ in range(4)] 

85 self.eps_gradeps[0][:] = 1.0 

86 

87 def allocate(self): 

88 pass 

89 

90 def todict(self): 

91 return {} 

92 

93 

94if __name__ == '__main__': 

95 import sys 

96 test_h(int(sys.argv[1]), sys.argv[2], 1)