Coverage for gpaw/test/solvation/test_sjm.py: 80%

1import pytest 

2from ase.build import fcc111 

3from gpaw import FermiDirac 

4from gpaw.mpi import size 

5from gpaw.new.ase_interface import GPAW 

6from gpaw.new.sjm import SJM 

7from gpaw.solvation import (EffectivePotentialCavity, GradientSurface, 

8 LinearDielectric, SurfaceInteraction) 

9from gpaw.solvation.sjm import SJM as OldSJM 

10from gpaw.solvation.sjm import SJMPower12Potential 

11 

12 

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

14def test_sjm(gpaw_new, in_tmp_dir, mode): 

15 if mode == 'pw': 

16 pytest.skip('Not working at the moment!') 

17 if not gpaw_new and size > 1: 

18 pytest.skip('https://gitlab.com/gpaw/gpaw/-/issues/1381') 

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

20 pytest.skip('Not implemented') 

21 # Solvent parameters 

22 u0 = 0.180 # eV 

23 epsinf = 78.36 # dielectric constant of water at 298 K 

24 gamma = 0.00114843767916 # 18.4*1e-3 * Pascal* m 

25 T = 298.15 # K 

26 

27 # Structure is created 

28 atoms = fcc111('Au', size=(1, 1, 3)) 

29 atoms.cell[2][2] = 15 

30 atoms.translate([0, 0, 6 - min(atoms.positions[:, 2])]) 

31 

32 # SJM parameters 

33 potential = 4.5 

34 tol = 0.02 

35 sj = {'target_potential': potential, 

36 'excess_electrons': -0.045, 

37 'jelliumregion': {'top': 14.5}, 

38 'tol': tol} 

39 

40 convergence = { 

41 'energy': 0.05 / 8., 

42 'density': 1e-4, 

43 'eigenstates': 1e-4} 

44 

45 params = dict( 

46 mode=mode, 

47 kpts=(2, 2, 1), 

48 xc='PBE', 

49 convergence=convergence, 

50 occupations=FermiDirac(0.1), 

51 txt=f'{gpaw_new}-{mode}.txt') 

52 

53 solvation = dict( 

54 cavity=EffectivePotentialCavity( 

55 effective_potential=SJMPower12Potential(u0=u0, 

56 unsolv_backside=False), 

57 temperature=T, 

58 surface_calculator=GradientSurface()), 

59 dielectric=LinearDielectric(epsinf=epsinf), 

60 interactions=[SurfaceInteraction(surface_tension=gamma)]) 

61 

62 if not gpaw_new: 

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

64 atoms.get_potential_energy() 

65 pot = atoms.calc.get_electrode_potential() 

66 else: 

67 atoms.calc = GPAW( 

68 **params, 

69 environment=SJM(**sj, **solvation)) 

70 atoms.get_potential_energy() 

71 pot = -atoms.calc.get_fermi_level() 

72 

73 assert abs(pot - potential) < tol 

74 

75 atoms.write('Au.traj') 

76 

77 atoms.calc.write(f'Au-{gpaw_new}-{mode}.gpw') 

78 if gpaw_new: 

79 calc = GPAW(f'Au-{gpaw_new}-{mode}.gpw') 

80 print(atoms.calc.environment) 

81 print(calc.environment) 

82 

83 if 0: # gpaw_new: 

84 import matplotlib.pyplot as plt 

85 import numpy as np 

86 x, y = np.array(atoms.calc.environment.jellium.history).T 

87 plt.plot(x, y) 

88 plt.show() 

89 if 0: 

90 v = atoms.calc.get_electrostatic_potential() 

91 import matplotlib.pyplot as plt 

92 import numpy as np 

93 plt.plot(np.linspace(0, atoms.cell[2, 2], v.shape[2], 0), v[0, 0]) 

94 plt.show() 

95 

96 

97if __name__ == '__main__': 

98 import sys 

99 test_sjm(sys.argv[1] == 'new', None, sys.argv[2])