Coverage for gpaw/test/solvation/test_vacuum.py: 83%

1import numpy as np 

2import pytest 

3from ase.build import molecule 

4 

5from gpaw import GPAW 

6from gpaw.solvation import (EffectivePotentialCavity, LinearDielectric, 

7 Power12Potential, SolvationGPAW) 

8from gpaw.utilities.adjust_cell import adjust_cell 

9 

10 

11def test_solvation_vacuum(): 

12 SKIP_REF_CALC = True 

13 

14 energy_eps = 0.0005 / 8 

15 forces_eps = 3e-2 

16 

17 h = 0.3 

18 vac = 3.0 

19 u0 = 0.180 

20 T = 298.15 

21 

22 atomic_radii = {'H': 1.09} 

23 

24 atoms = molecule('H2O') 

25 adjust_cell(atoms, vac, h) 

26 

27 convergence = { 

28 'energy': energy_eps, 

29 'forces': forces_eps, 

30 'density': 10.0, 

31 'eigenstates': 10.0} 

32 

33 if not SKIP_REF_CALC: 

34 atoms.calc = GPAW(mode='fd', xc='LDA', h=h, convergence=convergence) 

35 Eref = atoms.get_potential_energy() 

36 print(Eref) 

37 Fref = atoms.get_forces() 

38 print(Fref) 

39 else: 

40 Eref = -11.9929 

41 Fref = np.array([[0.0, 0.0, -6.07500], 

42 [0.0, 1.60924, 0.05999], 

43 [0.0, -1.60924, 0.05999]]) 

44 

45 atoms.calc = SolvationGPAW( 

46 mode='fd', 

47 xc='LDA', 

48 h=h, 

49 convergence=convergence, 

50 cavity=EffectivePotentialCavity( 

51 effective_potential=Power12Potential(atomic_radii=atomic_radii, 

52 u0=u0), 

53 temperature=T), 

54 dielectric=LinearDielectric(epsinf=1.0)) 

55 Etest = atoms.get_potential_energy() 

56 if atoms.calc.old: 

57 Eeltest = atoms.calc.get_electrostatic_energy() 

58 else: 

59 Eeltest = Etest - atoms.calc.environment.interaction_energy() 

60 Ftest = atoms.get_forces() 

61 assert Etest == pytest.approx( 

62 Eref, abs=energy_eps * atoms.calc.get_number_of_electrons()) 

63 assert Ftest == pytest.approx(Fref, abs=forces_eps) 

64 assert Eeltest == pytest.approx(Etest, abs=0.0)