Coverage for /builds/kinetik161/ase/ase/optimize/precon/neighbors.py: 97.30%

1import numpy as np 

2 

3from ase.constraints import FixAtoms 

4from ase.filters import Filter 

5from ase.geometry.cell import cell_to_cellpar 

6from ase.neighborlist import neighbor_list 

7 

8 

9def get_neighbours(atoms, r_cut, self_interaction=False, 

10 neighbor_list=neighbor_list): 

11 """Return a list of pairs of atoms within a given distance of each other. 

12 

13 Uses ase.neighborlist.neighbour_list to compute neighbors. 

14 

15 Args: 

16 atoms: ase.atoms object to calculate neighbours for 

17 r_cut: cutoff radius (float). Pairs of atoms are considered neighbours 

18 if they are within a distance r_cut of each other (note that this 

19 is double the parameter used in the ASE's neighborlist module) 

20 neighbor_list: function (optional). Optionally replace the built-in 

21 ASE neighbour list with an alternative with the same call 

22 signature, e.g. `matscipy.neighbours.neighbour_list`. 

23 

24 Returns: a tuple (i_list, j_list, d_list, fixed_atoms): 

25 i_list, j_list: i and j indices of each neighbour pair 

26 d_list: absolute distance between the corresponding pair 

27 fixed_atoms: indices of any fixed atoms 

28 """ 

29 

30 if isinstance(atoms, Filter): 

31 atoms = atoms.atoms 

32 

33 i_list, j_list, d_list = neighbor_list('ijd', atoms, r_cut) 

34 

35 # filter out self-interactions (across PBC) 

36 if not self_interaction: 

37 mask = i_list != j_list 

38 i_list = i_list[mask] 

39 j_list = j_list[mask] 

40 d_list = d_list[mask] 

41 

42 # filter out bonds where 1st atom (i) in pair is fixed 

43 fixed_atoms = [] 

44 for constraint in atoms.constraints: 

45 if isinstance(constraint, FixAtoms): 

46 fixed_atoms.extend(list(constraint.index)) 

47 

48 return i_list, j_list, d_list, fixed_atoms 

49 

50 

51def estimate_nearest_neighbour_distance(atoms, 

52 neighbor_list=neighbor_list): 

53 """ 

54 Estimate nearest neighbour distance r_NN 

55 

56 Args: 

57 atoms: Atoms object 

58 neighbor_list: function (optional). Optionally replace the built-in 

59 ASE neighbour list with an alternative with the same call 

60 signature, e.g. `matscipy.neighbours.neighbour_list`. 

61 

62 Returns: 

63 rNN: float 

64 Nearest neighbour distance 

65 """ 

66 

67 if isinstance(atoms, Filter): 

68 atoms = atoms.atoms 

69 

70 # start_time = time.time() 

71 # compute number of neighbours of each atom. If any atom doesn't 

72 # have a neighbour we increase the cutoff and try again, until our 

73 # cutoff exceeds the size of the system 

74 r_cut = 1.0 

75 phi = (1.0 + np.sqrt(5.0)) / 2.0 # Golden ratio 

76 

77 # cell lengths and angles 

78 a, b, c, alpha, beta, gamma = cell_to_cellpar(atoms.cell) 

79 extent = [a, b, c] 

80 # print('estimate_nearest_neighbour_distance(): extent=%r' % extent) 

81 

82 while r_cut < 2.0 * max(extent): 

83 # print('estimate_nearest_neighbour_distance(): ' 

84 # 'calling neighbour_list with r_cut=%.2f A' % r_cut) 

85 i, j, rij, fixed_atoms = get_neighbours( 

86 atoms, r_cut, self_interaction=True, 

87 neighbor_list=neighbor_list) 

88 if len(i) != 0: 

89 nn_i = np.bincount(i, minlength=len(atoms)) 

90 if (nn_i != 0).all(): 

91 break 

92 r_cut *= phi 

93 else: 

94 raise RuntimeError('increased r_cut to twice system extent without ' 

95 'finding neighbours for all atoms. This can ' 

96 'happen if your system is too small; try ' 

97 'setting r_cut manually') 

98 

99 # maximum of nearest neighbour distances 

100 nn_distances = [np.min(rij[i == I]) for I in range(len(atoms))] 

101 r_NN = np.max(nn_distances) 

102 

103 # print('estimate_nearest_neighbour_distance(): got r_NN=%.3f in %s s' % 

104 # (r_NN, time.time() - start_time)) 

105 return r_NN