Coverage for gpaw/test/solvation/test_swap_atoms.py: 98%

1from gpaw.utilities.adjust_cell import adjust_cell 

2from ase.build import molecule 

3from ase.units import Pascal, m 

4from gpaw.solvation import ( 

5 SolvationGPAW, 

6 EffectivePotentialCavity, 

7 Power12Potential, 

8 LinearDielectric, 

9 GradientSurface, 

10 SurfaceInteraction) 

11 

12 

13def test_solvation_swap_atoms(): 

14 h = 0.3 

15 vac = 3.0 

16 u0 = 0.180 

17 epsinf = 80. 

18 st = 18.4 * 1e-3 * Pascal * m 

19 T = 298.15 

20 

21 atomic_radii = {'H': 1.09} 

22 

23 convergence = { 

24 'energy': 0.1 / 8., 

25 'density': 10., 

26 'eigenstates': 10., 

27 } 

28 

29 atoms = molecule('H2O') 

30 adjust_cell(atoms, vac, h) 

31 

32 calc = SolvationGPAW( 

33 mode='fd', xc='LDA', h=h, convergence=convergence, 

34 cavity=EffectivePotentialCavity( 

35 effective_potential=Power12Potential(atomic_radii, u0), 

36 temperature=T, 

37 surface_calculator=GradientSurface() 

38 ), 

39 dielectric=LinearDielectric(epsinf=epsinf), 

40 interactions=[SurfaceInteraction(surface_tension=st)] 

41 ) 

42 atoms.calc = calc 

43 atoms.get_potential_energy() 

44 atoms.get_forces() 

45 

46 def env(calc): 

47 if calc.old: 

48 return calc.hamiltonian 

49 return calc.environment 

50 

51 eps_gradeps = env(calc).dielectric.eps_gradeps 

52 

53 # same molecules, different cell, reallocate 

54 atoms = molecule('H2O') 

55 atoms.positions[0][0] = atoms.positions[0][0] - 1. 

56 adjust_cell(atoms, vac, h) 

57 atoms.calc = calc 

58 atoms.get_potential_energy() 

59 atoms.get_forces() 

60 

61 assert env(calc).dielectric.eps_gradeps is not eps_gradeps 

62 eps_gradeps = env(calc).dielectric.eps_gradeps 

63 

64 # small position change, no reallocate 

65 atoms.positions[0][0] = atoms.positions[0][0] + 1e-2 

66 atoms.get_potential_energy() 

67 atoms.get_forces() 

68 assert env(calc).dielectric.eps_gradeps is eps_gradeps 

69 eps_gradeps = env(calc).dielectric.eps_gradeps 

70 radii = env(calc).cavity.effective_potential.r12_a 

71 

72 # completely different atoms object, reallocate, read new radii 

73 atoms = molecule('NH3') 

74 adjust_cell(atoms, vac, h) 

75 atoms.calc = calc 

76 atoms.get_potential_energy() 

77 atoms.get_forces() 

78 assert env(calc).dielectric.eps_gradeps is not eps_gradeps 

79 assert env(calc).cavity.effective_potential.r12_a is not radii