Coverage for gpaw/test/solvation/test_cip.py: 100%

1import pytest 

2from gpaw.solvation.sjm import SJM, SJMPower12Potential 

3import numpy as np 

4from ase.build import fcc111 

5from gpaw import FermiDirac 

6 

7# Import solvation modules 

8from gpaw.solvation import ( 

9 EffectivePotentialCavity, 

10 LinearDielectric, 

11 GradientSurface, 

12 SurfaceInteraction) 

13 

14 

15@pytest.mark.slow 

16@pytest.mark.old_gpaw_only 

17def test_cip(in_tmp_dir): 

18 # Solvent parameters 

19 u0 = 0.180 # eV 

20 epsinf = 78.36 # Dielectric constant of water at 298 K 

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

22 T = 298.15 # K 

23 

24 # Structure is created 

25 atoms = fcc111('Au', size=(1, 1, 4)) 

26 atoms.cell[2][2] = 18 

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

28 

29 # SJM parameters 

30 potential = 4.5 

31 tol = 0.02 

32 sj_calib = {'excess_electrons': 0, 

33 'pot_ref': 'CIP', 

34 'cip': {'autoinner': {'nlayers': 4, 

35 'threshold': 0.01}}} 

36 

37 convergence = { 

38 'energy': 0.05 / 8., 

39 'density': 1e-4, 

40 'eigenstates': 1e-4} 

41 

42 # Calculator 

43 calc = SJM( 

44 mode='lcao', 

45 sj=sj_calib, 

46 txt='sjmwf.test', 

47 gpts=(8, 8, 48), 

48 kpts=(2, 2, 1), 

49 xc='PBE', 

50 convergence=convergence, 

51 occupations=FermiDirac(0.1), 

52 cavity=EffectivePotentialCavity( 

53 effective_potential=SJMPower12Potential(u0=u0), 

54 temperature=T, 

55 surface_calculator=GradientSurface()), 

56 dielectric=LinearDielectric(epsinf=epsinf), 

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

58 # Run calibration calculations with zero charge conditions 

59 atoms.calc = calc 

60 atoms.get_potential_energy() 

61 phi_pzc = calc.get_inner_potential(atoms) 

62 mu_pzc = calc.get_fermi_level() 

63 

64 assert np.isclose(mu_pzc, -4.39, 1e-1, 1e-1) 

65 # Reset for CIP using computed potential of zero charge and fermi level 

66 potential = 4.5 

67 tol = 0.02 

68 sj_cip = {'target_potential': potential, 

69 'excess_electrons': 0., 

70 'jelliumregion': {'top': 17.9}, 

71 'tol': tol, 

72 'pot_ref': 'CIP', 

73 'cip': {'autoinner': {'nlayers': 4, 

74 'threshold': 0.01}, 

75 'inner_region': None, 

76 'mu_pzc': mu_pzc, 

77 'phi_pzc': phi_pzc, 

78 'filter': 10}} 

79 

80 calc = SJM(mode='lcao', 

81 sj=sj_cip, 

82 txt='sjmcip.test', 

83 gpts=(8, 8, 48), 

84 kpts=(2, 2, 1), 

85 xc='PBE', 

86 convergence=convergence, 

87 occupations=FermiDirac(0.1), 

88 cavity=EffectivePotentialCavity( 

89 effective_potential=SJMPower12Potential(u0=u0), 

90 temperature=T, 

91 surface_calculator=GradientSurface()), 

92 dielectric=LinearDielectric(epsinf=epsinf), 

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

94 

95 atoms.calc = calc 

96 atoms.get_potential_energy() 

97 assert (abs(calc.get_electrode_potential(sj_cip['pot_ref'], 

98 return_referenced=True)) 

99 - potential) < tol