Coverage for /builds/kinetik161/ase/ase/build/attach.py: 86.05%

1import numpy as np 

2 

3from ase.parallel import world, broadcast 

4from ase.geometry import get_distances 

5 

6 

7def random_unit_vector(rng): 

8 """Random unit vector equally distributed on the sphere 

9 

10 Parameter 

11 --------- 

12 rng: random number generator object 

13 """ 

14 ct = -1 + 2 * rng.random() 

15 phi = 2 * np.pi * rng.random() 

16 st = np.sqrt(1 - ct**2) 

17 return np.array([st * np.cos(phi), st * np.sin(phi), ct]) 

18 

19 

20def nearest(atoms1, atoms2, cell=None, pbc=None): 

21 """Return indices of nearest atoms""" 

22 p1 = atoms1.get_positions() 

23 p2 = atoms2.get_positions() 

24 vd_aac, d2_aa = get_distances(p1, p2, cell, pbc) 

25 i1, i2 = np.argwhere(d2_aa == d2_aa.min())[0] 

26 return i1, i2, vd_aac[i1, i2] 

27 

28 

29def attach(atoms1, atoms2, distance, direction=(1, 0, 0), 

30 maxiter=50, accuracy=1e-5): 

31 """Attach two structures 

32 

33 Parameters 

34 ---------- 

35 atoms1: Atoms 

36 cell and pbc of this object are used 

37 atoms2: Atoms 

38 distance: float 

39 minimal distance (Angstrom) 

40 direction: unit vector (3 floats) 

41 relative direction between center of masses 

42 maxiter: int 

43 maximal number of iterations to get required distance, default 100 

44 accuracy: float 

45 required accuracy for minimal distance (Angstrom), default 1e-5 

46 

47 Returns 

48 ------- 

49 Joined structure as an atoms object. 

50 """ 

51 atoms = atoms1.copy() 

52 atoms2 = atoms2.copy() 

53 

54 direction = np.array(direction, dtype=float) 

55 direction /= np.linalg.norm(direction) 

56 assert len(direction) == 3 

57 dist2 = distance**2 

58 

59 i1, i2, dv_c = nearest(atoms, atoms2, atoms.cell, atoms.pbc) 

60 

61 for i in range(maxiter): 

62 dv2 = (dv_c**2).sum() 

63 

64 vcost = np.dot(dv_c, direction) 

65 a = np.sqrt(max(0, dist2 - dv2 + vcost**2)) 

66 move = a - vcost 

67 if abs(move) < accuracy: 

68 atoms += atoms2 

69 return atoms 

70 

71 # we need to move 

72 atoms2.translate(direction * move) 

73 i1, i2, dv_c = nearest(atoms, atoms2, atoms.cell, atoms.pbc) 

74 

75 raise RuntimeError('attach did not converge') 

76 

77 

78def attach_randomly(atoms1, atoms2, distance, 

79 rng=np.random): 

80 """Randomly attach two structures with a given minimal distance 

81 

82 Parameters 

83 ---------- 

84 atoms1: Atoms object 

85 atoms2: Atoms object 

86 distance: float 

87 Required distance 

88 rng: random number generator object 

89 defaults to np.random.RandomState() 

90 

91 Returns 

92 ------- 

93 Joined structure as an atoms object. 

94 """ 

95 atoms2 = atoms2.copy() 

96 atoms2.rotate('x', random_unit_vector(rng), 

97 center=atoms2.get_center_of_mass()) 

98 return attach(atoms1, atoms2, distance, 

99 direction=random_unit_vector(rng)) 

100 

101 

102def attach_randomly_and_broadcast(atoms1, atoms2, distance, 

103 rng=np.random, 

104 comm=world): 

105 """Randomly attach two structures with a given minimal distance 

106 and ensure that these are distributed. 

107 

108 Parameters 

109 ---------- 

110 atoms1: Atoms object 

111 atoms2: Atoms object 

112 distance: float 

113 Required distance 

114 rng: random number generator object 

115 defaults to np.random.RandomState() 

116 comm: communicator to distribute 

117 Communicator to distribute the structure, default: world 

118 

119 Returns 

120 ------- 

121 Joined structure as an atoms object. 

122 """ 

123 if comm.rank == 0: 

124 joined = attach_randomly(atoms1, atoms2, distance, rng) 

125 broadcast(joined, 0, comm=comm) 

126 else: 

127 joined = broadcast(None, 0, comm) 

128 return joined