Coverage for gpaw/test/solvation/test_solvation_api.py: 17%

1# XXX This test is a use case/acceptance test to help rewrite the api 

2# XXX and not included in the test suite. 

3# XXX TODO: make an example/documentation out of this test 

4# XXX when the api has changed and the test passes 

5import pytest 

6from ase.build import molecule 

7from ase.units import Pascal, m, Bohr 

8from ase.parallel import parprint 

9from gpaw.solvation import ( 

10 # calculator 

11 SolvationGPAW, 

12 # cavities 

13 EffectivePotentialCavity, 

14 ADM12SmoothStepCavity, 

15 FG02SmoothStepCavity, 

16 # custom classes for the cavities 

17 Power12Potential, 

18 ElDensity, 

19 SSS09Density, 

20 # dielectric 

21 LinearDielectric, 

22 # non-el interactions 

23 SurfaceInteraction, 

24 VolumeInteraction, 

25 LeakedDensityInteraction, 

26 # surface and volume calculators 

27 GradientSurface, 

28 KB51Volume, 

29) 

30# poisson solver 

31from gpaw.solvation.poisson import ADM12PoissonSolver 

32 

33 

34# references for custom classes: 

35# KB51 = J. G. Kirkwood and F. P. Buff, 

36# The Journal of Chemical Physics, vol. 19, no. 6, pp. 774--777, 1951 

37# FG02 = J.-L. Fattebert and F. Gygi, 

38# Journal of Computational Chemistry, vol. 23, no. 6, pp. 662--666, 2002 

39# SSS09 = V. M. Sanchez, M. Sued, and D. A. Scherlis, 

40# The Journal of Chemical Physics, vol. 131, no. 17, p. 174108, 2009 

41# ADM12 = O. Andreussi, I. Dabo, and N. Marzari, 

42# The Journal of Chemical Physics, vol. 136, no. 6, p. 064102, 2012 

43 

44 

45@pytest.mark.skip(reason='TODO') 

46def test_solvation_api(): 

47 # define some useful units (all api user units are ASE units!) 

48 dyn_per_cm = 1e-3 * Pascal * m 

49 Giga_Pascal = 1e9 * Pascal 

50 

51 # GPAW params (examples) 

52 # ---------------------- 

53 xc = 'PBE' 

54 h = 0.24 

55 vac = 4.0 

56 

57 # general solvation params (examples) 

58 # ----------------------------------- 

59 # electrostatic 

60 epsinf = 78.36 

61 # other interactions 

62 gamma = 72. * dyn_per_cm # surface tension 

63 p = -0.1 * Giga_Pascal # pressure 

64 V_leak = 1.0 # V (interaction energy E = V_leak * [charge outside cavity]) 

65 # only for volume calculations respecting compressibility 

66 T = 298.15 # K (also used for Boltzmann distribution) 

67 kappa_T = 4.53e-10 / Pascal 

68 

69 # effective potential cavity params (examples) 

70 # -------------------------------------------- 

71 u0 = 0.180 # eV 

72 

73 atomic_radii = {'H': 1.09} 

74 

75 # density cavity params (examples) 

76 # -------------------------------- 

77 # ADM12 

78 rhomin = 0.0001 / Bohr ** 3 

79 rhomax = 0.0050 / Bohr ** 3 

80 # FG02, SSS09 

81 rho0 = 0.0004 / Bohr ** 3 

82 beta = 1.3 

83 rho0_fake = 1.0 / Bohr ** 3 

84 beta_fake = 2.4 

85 

86 atoms = molecule('H2O') 

87 atoms.center(vacuum=vac) 

88 

89 def print_results(atoms): 

90 parprint('E = %.3f eV' % (atoms.get_potential_energy(), )) 

91 parprint('V = %.3f Ang ** 3' % (atoms.calc.get_cavity_volume(), )) 

92 parprint('A = %.3f Ang ** 2' % (atoms.calc.get_cavity_surface(), )) 

93 parprint('Forces:') 

94 parprint(atoms.get_forces()) 

95 parprint('') 

96 

97 def get_base_kwargs() -> dict: 

98 return dict( 

99 mode='fd', xc=xc, h=h, 

100 dielectric=LinearDielectric(epsinf=epsinf), 

101 interactions=[ 

102 SurfaceInteraction(surface_tension=gamma), 

103 VolumeInteraction(pressure=p), 

104 LeakedDensityInteraction(voltage=V_leak)]) 

105 

106 # Cavity from 1 / r ** 12 effective potential 

107 atoms.calc = SolvationGPAW( 

108 **get_base_kwargs(), 

109 cavity=EffectivePotentialCavity( 

110 effective_potential=Power12Potential(atomic_radii=atomic_radii, 

111 u0=u0), 

112 temperature=T, 

113 surface_calculator=GradientSurface(), 

114 volume_calculator=KB51Volume(compressibility=kappa_T, 

115 temperature=T))) 

116 print_results(atoms) 

117 

118 # Cavity from electron density a la ADM12 

119 atoms.calc = SolvationGPAW( 

120 **get_base_kwargs(), 

121 poissonsolver=ADM12PoissonSolver(), 

122 cavity=ADM12SmoothStepCavity( 

123 rhomin, rhomax, epsinf, 

124 density=ElDensity(), 

125 surface_calculator=GradientSurface(), 

126 volume_calculator=KB51Volume(compressibility=kappa_T, 

127 temperature=T))) 

128 print_results(atoms) 

129 

130 # Cavity from electron density a la FG02 

131 atoms.calc = SolvationGPAW( 

132 **get_base_kwargs(), 

133 cavity=FG02SmoothStepCavity( 

134 rho0, beta, 

135 density=ElDensity(), 

136 surface_calculator=GradientSurface(), 

137 volume_calculator=KB51Volume(compressibility=kappa_T, 

138 temperature=T))) 

139 print_results(atoms) 

140 

141 # Cavity from fake electron density a la SSS09 

142 atoms.calc = SolvationGPAW( 

143 **get_base_kwargs(), 

144 cavity=FG02SmoothStepCavity( 

145 rho0_fake, beta_fake, 

146 density=SSS09Density(atomic_radii=atomic_radii), 

147 surface_calculator=GradientSurface(), 

148 volume_calculator=KB51Volume(compressibility=kappa_T, 

149 temperature=T))) 

150 print_results(atoms)