Coverage for /builds/debichem-team/python-ase/ase/calculators/morse.py: 100.00%

1import numpy as np 

2 

3from ase.calculators.calculator import Calculator, all_changes 

4from ase.neighborlist import neighbor_list 

5from ase.stress import full_3x3_to_voigt_6_stress 

6 

7 

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

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

10 

11 Ported from JuLIP.jl. 

12 

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

14 

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

16 

17 Parameters 

18 ---------- 

19 r : np.ndarray 

20 Distances between atoms. 

21 r0 : float 

22 Inner cutoff radius. 

23 r1 : float 

24 Outder cutoff radius. 

25 

26 Returns 

27 ------- 

28 np.ndarray 

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

30 

31 """"" 

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

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

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

35 ) 

36 

37 

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

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

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

41 return -( 

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

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

44 / (r1 - r0) 

45 ) 

46 

47 

48class MorsePotential(Calculator): 

49 """Morse potential.""" 

50 

51 implemented_properties = [ 

52 'energy', 'energies', 'free_energy', 'forces', 'stress', 

53 ] 

54 default_parameters = {'epsilon': 1.0, 

55 'rho0': 6.0, 

56 'r0': 1.0, 

57 'rcut1': 1.9, 

58 'rcut2': 2.7} 

59 nolabel = True 

60 

61 def __init__(self, **kwargs): 

62 r""" 

63 

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

65 

66 .. math:: 

67 

68 V_{ij} = \epsilon \left( 

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

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

71 \right) 

72 

73 Parameters 

74 ---------- 

75 epsilon : float, default 1.0 

76 Absolute minimum depth. 

77 r0 : float, default 1.0 

78 Minimum distance. 

79 rho0 : float, default 6.0 

80 Exponential prefactor. 

81 The force constant in the potential minimum is given by 

82 

83 .. math:: 

84 

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

86 

87 rcut1 : float, default 1.9 

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

89 rcut2 : float, default 2.7 

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

91 

92 Notes 

93 ----- 

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

95 

96 """ 

97 Calculator.__init__(self, **kwargs) 

98 

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

100 system_changes=all_changes): 

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

102 epsilon = self.parameters['epsilon'] 

103 rho0 = self.parameters['rho0'] 

104 r0 = self.parameters['r0'] 

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

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

107 

108 number_of_atoms = len(self.atoms) 

109 

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

111 

112 i, _j, d, D = neighbor_list('ijdD', atoms, rcut2) 

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

114 

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

116 

117 cutoff_fn = fcut(d, rcut1, rcut2) 

118 d_cutoff_fn = fcut_d(d, rcut1, rcut2) 

119 

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

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

122 i, 

123 weights=0.5 * (pairwise_energies * cutoff_fn), 

124 minlength=number_of_atoms, 

125 ) 

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

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

128 

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

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

131 

132 # smoothened `de` 

133 de = de * cutoff_fn + pairwise_energies * d_cutoff_fn 

134 

135 de_vec = (de * dhat).T 

136 for dim in range(3): 

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

138 i, 

139 weights=de_vec[:, dim], 

140 minlength=number_of_atoms, 

141 ) 

142 self.results['forces'] = forces 

143 

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

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

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