Coverage for gpaw/test/solvation/test_adm12.py: 91%

35 statements  

« prev     ^ index     » next       coverage.py v7.7.1, created at 2025-07-19 00:19 +0000

1""" 

2tests solvation parameters as in 

3 

4O. Andreussi, I. Dabo, and N. Marzari, 

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

6""" 

7 

8from gpaw import GPAW 

9from gpaw.utilities.adjust_cell import adjust_cell 

10import pytest 

11from ase.build import molecule 

12from ase.units import mol, kcal, Pascal, m, Bohr 

13from gpaw.solvation import ( 

14 SolvationGPAW, 

15 ADM12SmoothStepCavity, 

16 LinearDielectric, 

17 GradientSurface, 

18 SurfaceInteraction, 

19 KB51Volume, 

20 VolumeInteraction, 

21 ElDensity 

22) 

23from gpaw.solvation.poisson import ADM12PoissonSolver 

24import warnings 

25 

26SKIP_VAC_CALC = True 

27 

28h = 0.24 

29vac = 4.0 

30 

31epsinf = 78.36 

32rhomin = 0.0001 / Bohr ** 3 

33rhomax = 0.0050 / Bohr ** 3 

34st = 50. * 1e-3 * Pascal * m 

35p = -0.35 * 1e9 * Pascal 

36convergence = { 

37 'energy': 0.05 / 8., 

38 'density': 10., 

39 'eigenstates': 10., 

40} 

41 

42 

43def test_solvation_adm12(): 

44 atoms = molecule('H2O') 

45 adjust_cell(atoms, vac, h) 

46 

47 if not SKIP_VAC_CALC: 

48 atoms.calc = GPAW(mode='fd', xc='PBE', h=h, convergence=convergence) 

49 Evac = atoms.get_potential_energy() 

50 print(Evac) 

51 else: 

52 # h=0.24, vac=4.0, setups: 0.9.11271, convergence: only energy 0.05 / 8 

53 Evac = -14.857414548 

54 

55 with warnings.catch_warnings(): 

56 # ignore production code warning for ADM12PoissonSolver 

57 warnings.simplefilter("ignore") 

58 psolver = ADM12PoissonSolver() 

59 

60 atoms.calc = SolvationGPAW( 

61 mode='fd', 

62 xc='PBE', h=h, poissonsolver=psolver, convergence=convergence, 

63 cavity=ADM12SmoothStepCavity( 

64 rhomin=rhomin, rhomax=rhomax, epsinf=epsinf, 

65 density=ElDensity(), 

66 surface_calculator=GradientSurface(), 

67 volume_calculator=KB51Volume() 

68 ), 

69 dielectric=LinearDielectric(epsinf=epsinf), 

70 interactions=[ 

71 SurfaceInteraction(surface_tension=st), 

72 VolumeInteraction(pressure=p) 

73 ] 

74 ) 

75 Ewater = atoms.get_potential_energy() 

76 assert atoms.calc.get_number_of_iterations() < 40 

77 atoms.get_forces() 

78 DGSol = (Ewater - Evac) / (kcal / mol) 

79 print('Delta Gsol: %s kcal / mol' % DGSol) 

80 

81 assert DGSol == pytest.approx(-6.3, abs=2.)