Coverage for /builds/ase/ase/ase/calculators/morse.py: 100.00%

1# fmt: off 

2 

3import numpy as np 

4 

5from ase.calculators.calculator import Calculator, all_changes 

6from ase.neighborlist import neighbor_list as ase_neighbor_list 

7from ase.stress import full_3x3_to_voigt_6_stress 

8 

9 

10def fcut(r: np.ndarray, r0: float, r1: float) -> np.ndarray: 

11 """Piecewise quintic C^{2,1} regular polynomial for use as a smooth cutoff. 

12 

13 Ported from JuLIP.jl. 

14 

15 https://github.com/JuliaMolSim/JuLIP.jl/blob/master/src/cutoffs.jl 

16 

17 https://en.wikipedia.org/wiki/Smoothstep 

18 

19 Parameters 

20 ---------- 

21 r : np.ndarray 

22 Distances between atoms. 

23 r0 : float 

24 Inner cutoff radius. 

25 r1 : float 

26 Outder cutoff radius. 

27 

28 Returns 

29 ------- 

30 np.ndarray 

31 Sigmoid-like function smoothly interpolating (r0, 1) and (r1, 0). 

32 

33 """"" 

34 s = 1.0 - (r - r0) / (r1 - r0) 

35 return (s >= 1.0) + ((s > 0.0) & (s < 1.0)) * ( 

36 6.0 * s**5 - 15.0 * s**4 + 10.0 * s**3 

37 ) 

38 

39 

40def fcut_d(r: np.ndarray, r0: float, r1: float) -> np.ndarray: 

41 """Derivative of fcut() function defined above.""" 

42 s = 1.0 - (r - r0) / (r1 - r0) 

43 return -( 

44 ((s > 0.0) & (s < 1.0)) 

45 * (30.0 * s**4 - 60.0 * s**3 + 30.0 * s**2) 

46 / (r1 - r0) 

47 ) 

48 

49 

50class MorsePotential(Calculator): 

51 """Morse potential.""" 

52 

53 implemented_properties = [ 

54 'energy', 'energies', 'free_energy', 'forces', 'stress', 

55 ] 

56 default_parameters = {'epsilon': 1.0, 

57 'rho0': 6.0, 

58 'r0': 1.0, 

59 'rcut1': 1.9, 

60 'rcut2': 2.7} 

61 nolabel = True 

62 

63 def __init__(self, neighbor_list=ase_neighbor_list, **kwargs): 

64 r""" 

65 

66 The pairwise energy between atoms *i* and *j* is given by 

67 

68 .. math:: 

69 

70 V_{ij} = \epsilon \left( 

71 \mathrm{e}^{-2 \rho_0 (r_{ij} / r_0 - 1)} 

72 - 2 \mathrm{e}^{- \rho_0 (r_{ij} / r_0 - 1)} 

73 \right) 

74 

75 Parameters 

76 ---------- 

77 epsilon : float, default 1.0 

78 Absolute minimum depth. 

79 r0 : float, default 1.0 

80 Minimum distance. 

81 rho0 : float, default 6.0 

82 Exponential prefactor. 

83 The force constant in the potential minimum is given by 

84 

85 .. math:: 

86 

87 k = 2 \epsilon \left(\frac{\rho_0}{r_0}\right)^2. 

88 

89 rcut1 : float, default 1.9 

90 Distance starting a smooth cutoff normalized by ``r0``. 

91 rcut2 : float, default 2.7 

92 Distance ending a smooth cutoff normalized by ``r0``. 

93 neighbor_list : callable, optional 

94 neighbor_list function compatible with 

95 ase.neighborlist.neighbor_list 

96 

97 Notes 

98 ----- 

99 The default values are chosen to be similar as Lennard-Jones. 

100 

101 """ 

102 self.neighbor_list = neighbor_list 

103 Calculator.__init__(self, **kwargs) 

104 

105 def calculate(self, atoms=None, properties=['energy'], 

106 system_changes=all_changes): 

107 Calculator.calculate(self, atoms, properties, system_changes) 

108 epsilon = self.parameters['epsilon'] 

109 rho0 = self.parameters['rho0'] 

110 r0 = self.parameters['r0'] 

111 rcut1 = self.parameters['rcut1'] * r0 

112 rcut2 = self.parameters['rcut2'] * r0 

113 

114 number_of_atoms = len(self.atoms) 

115 

116 forces = np.zeros((number_of_atoms, 3)) 

117 

118 i, _j, d, D = self.neighbor_list('ijdD', atoms, rcut2) 

119 dhat = (D / d[:, None]).T 

120 

121 expf = np.exp(rho0 * (1.0 - d / r0)) 

122 

123 cutoff_fn = fcut(d, rcut1, rcut2) 

124 d_cutoff_fn = fcut_d(d, rcut1, rcut2) 

125 

126 pairwise_energies = epsilon * expf * (expf - 2.0) 

127 self.results['energies'] = np.bincount( 

128 i, 

129 weights=0.5 * (pairwise_energies * cutoff_fn), 

130 minlength=number_of_atoms, 

131 ) 

132 self.results['energy'] = self.results['energies'].sum() 

133 self.results['free_energy'] = self.results['energy'] 

134 

135 # derivatives of `pair_energies` with respect to `d` 

136 de = (-2.0 * epsilon * rho0 / r0) * expf * (expf - 1.0) 

137 

138 # smoothened `de` 

139 de = de * cutoff_fn + pairwise_energies * d_cutoff_fn 

140 

141 de_vec = (de * dhat).T 

142 for dim in range(3): 

143 forces[:, dim] = np.bincount( 

144 i, 

145 weights=de_vec[:, dim], 

146 minlength=number_of_atoms, 

147 ) 

148 self.results['forces'] = forces 

149 

150 if self.atoms.cell.rank == 3: 

151 stress = 0.5 * (D.T @ de_vec) / self.atoms.get_volume() 

152 self.results['stress'] = full_3x3_to_voigt_6_stress(stress)