Coverage for gpaw/utilities/gpts.py: 96%

1import warnings 

2 

3import numpy as np 

4 

5from ase.units import Bohr, Hartree 

6 

7from gpaw.fftw import get_efficient_fft_size 

8from gpaw.utilities import h2gpts 

9from gpaw.wavefunctions.fd import FD 

10 

11 

12def get_number_of_grid_points(cell_cv, 

13 h=None, 

14 mode=None, 

15 realspace=None, 

16 symmetry=None): 

17 if mode is None: 

18 mode = FD() 

19 

20 if realspace is None: 

21 realspace = mode.name != 'pw' 

22 

23 if h is None: 

24 if mode.name == 'pw': 

25 h = np.pi / (4 * mode.ecut)**0.5 

26 elif mode.name == 'lcao' and not realspace: 

27 h = np.pi / (4 * 340 / Hartree)**0.5 

28 else: 

29 h = 0.2 / Bohr 

30 

31 if realspace or mode.name == 'fd': 

32 N_c = h2gpts(h, cell_cv, 4) 

33 else: 

34 N_c = np.ceil((cell_cv**2).sum(1)**0.5 / h).astype(int) 

35 if symmetry is None: 

36 N_c = np.array([get_efficient_fft_size(N) for N in N_c]) 

37 else: 

38 N_c = np.array([get_efficient_fft_size(N, n) 

39 for N, n in zip(N_c, symmetry.gcd_c)]) 

40 

41 if symmetry is not None: 

42 ok = symmetry.check_grid(N_c) 

43 if not ok: 

44 warnings.warn( 

45 'Initial realspace grid ' 

46 '({},{},{}) inconsistent with symmetries.'.format(*N_c)) 

47 # The grid is not symmetric enough. The essential problem 

48 # is that we should start at some other Nmin_c and possibly with 

49 # other gcd_c when getting the most efficient fft grids 

50 gcd_c = symmetry.gcd_c.copy() 

51 Nmin_c = N_c.copy() 

52 for i in range(3): 

53 for op_cc in symmetry.op_scc: 

54 for i, o in enumerate((op_cc.T).flat): 

55 i1, i2 = np.unravel_index(i, (3, 3)) 

56 if i1 == i2: 

57 continue 

58 if o: 

59 # The axes are related and therefore they share 

60 # lowest common multiple of gcd 

61 gcd = np.lcm.reduce(gcd_c[[i1, i2]]) 

62 gcd_c[[i1, i2]] = gcd 

63 # We just take the maximum of the two axes to make 

64 # sure that they are divisible always 

65 Nmin = np.max([Nmin_c[i1], Nmin_c[i2]]) 

66 Nmin_c[i1] = Nmin 

67 Nmin_c[i2] = Nmin 

68 

69 N_c = np.array([get_efficient_fft_size(N, n) 

70 for N, n in zip(Nmin_c, gcd_c)]) 

71 warnings.warn('Using symmetrized grid: ({},{},{}).\n'.format(*N_c)) 

72 ok = symmetry.check_grid(N_c) 

73 assert ok, ('Grid still not constistent with symmetries ' 

74 '({},{},{})'.format(*N_c)) 

75 

76 return N_c