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

1from gpaw.solvation.cavity import get_pbc_positions 

2from ase.build import molecule 

3from ase.units import Bohr 

4import numpy as np 

5 

6 

7def test_solvation_pbc_pos_repeat(): 

8 def get_nrepeats(nrepeats_c): 

9 return np.prod(np.array(nrepeats_c) * 2 + 1) 

10 

11 def check_valid_repeat(pos_v, rpos_v, nrepeats_c, cell_cv): 

12 nrepeats_c = np.array(nrepeats_c) 

13 diff_v = rpos_v - pos_v 

14 actual_nrepeats_c = np.dot(diff_v, np.linalg.inv(cell_cv)) 

15 r_actual_nrepeats_c = np.around(actual_nrepeats_c) 

16 # check for integer multiple 

17 assert np.allclose(r_actual_nrepeats_c, actual_nrepeats_c) 

18 actual_nrepeats_c = np.abs(r_actual_nrepeats_c.astype(int)) 

19 # check for too large repeats 

20 assert (actual_nrepeats_c <= nrepeats_c).all() 

21 

22 def check(pos_av, rpos_aav, nrepeats_c, cell_cv): 

23 assert len(pos_av) == len(rpos_aav) 

24 total_repeats = get_nrepeats(nrepeats_c) 

25 for a, pos_v in enumerate(pos_av): 

26 rpos_av = rpos_aav[a] 

27 assert len(rpos_av) == total_repeats 

28 for rpos_v in rpos_av: 

29 check_valid_repeat(pos_v, rpos_v, nrepeats_c, cell_cv) 

30 

31 cells = ( 

32 ((10., .0, .0), (.0, 12., .0), (.0, .0, 12.)), # orthogonal 

33 ((10., 1., .0), (.0, 12., 1.), (1., .0, 11.5))) # non-orthogonal 

34 

35 for cell in cells: 

36 atoms = molecule('H2O') 

37 atoms.center(vacuum=5.) 

38 atoms.set_cell(cell) 

39 pos_av = atoms.positions / Bohr 

40 cell_cv = atoms.get_cell() / Bohr 

41 Rcell_c = np.sqrt(np.sum(cell_cv ** 2, axis=1)) 

42 

43 # non periodic should never repeat 

44 atoms.pbc = np.zeros((3, ), dtype=bool) 

45 pos_aav = get_pbc_positions(atoms, 1e10 / Bohr) 

46 check(pos_av, pos_aav, (0, 0, 0), cell_cv) 

47 

48 # periodic, zero cutoff should not repeat 

49 atoms.pbc = np.ones((3, ), dtype=bool) 

50 pos_aav = get_pbc_positions(atoms, .0 / Bohr) 

51 check(pos_av, pos_aav, (0, 0, 0), cell_cv) 

52 

53 # periodic, cutoff <= cell size should repeat once 

54 atoms.pbc = np.ones((3, ), dtype=bool) 

55 pos_aav = get_pbc_positions(atoms, 0.99 * Rcell_c.min()) 

56 check(pos_av, pos_aav, (1, 1, 1), cell_cv) 

57 

58 # periodic, cutoff > cell size should repeat twice 

59 atoms.pbc = np.ones((3, ), dtype=bool) 

60 pos_aav = get_pbc_positions(atoms, 1.01 * Rcell_c.max()) 

61 check(pos_av, pos_aav, (2, 2, 2), cell_cv) 

62 

63 # mixed bc, cutoff > 2 * cell size should repeat three times 

64 for i in range(8): 

65 pbc = np.array([int(p) for p in np.binary_repr(i, 3)]) 

66 atoms.pbc = pbc 

67 pos_aav = get_pbc_positions(atoms, 2.01 * Rcell_c.max()) 

68 check(pos_av, pos_aav, pbc * 3, cell_cv)