Coverage for gamdpy/interactions/potential_functions/SAAP.py: 100%

1import numba 

2from math import exp 

3 

4def SAAP(dist, params): 

5 """ The SAAP potential 

6 is a pair potential for noble elements, its parameters 

7 are derived from fitting on data obtained from ab initio methods (coupled cluster): 

8 

9 .. math:: 

10 

11 u(x=r/\\sigma)/\\epsilon = (a_0\\exp(a_1 x)/x + a_2\\exp(a_3 x) + a_4) / (1+a_5 x^6) 

12 

13 The six :math:`a_i` parameters are given in units of eps and sigma. 

14 

15 Reference: https://doi.org/10.1063/1.5085420 

16 

17 Parameters 

18 ---------- 

19 

20 dist : float 

21 Distance between particles 

22 

23 params : array-like 

24 a₀, a₁, a₂, a₃, a₄, a₅, σ, ε 

25 

26 """ 

27 

28 # Extract parameters compatibly with numba, in float32 precision 

29 a0 = numba.float32(params[0]) 

30 a1 = numba.float32(params[1]) 

31 a2 = numba.float32(params[2]) 

32 a3 = numba.float32(params[3]) 

33 a4 = numba.float32(params[4]) 

34 a5 = numba.float32(params[5]) 

35 sigma = numba.float32(params[6]) 

36 eps = numba.float32(params[7]) 

37 

38 # Definition of a reduced distance to make all consistent with the fact that 

39 # the various a parameters are given in units of eps and sigma 

40 r = dist / sigma 

41 # Compute helper variables 

42 one = numba.float32(1.0) 

43 inv_r = one/r 

44 

45 exp_a1_r = exp(a1 * r) 

46 a0_exp_r = a0 * exp_a1_r * inv_r # a0 exp(a1 r)/r 

47 a2_exp = a2 * exp(a3 * r) 

48 inverse_one_a5 = one / (one + a5 * r**6) 

49 

50 # Compute pair potential energy, pair force and pair curvature 

51 

52 # SAAP computation for r 

53 u = eps * (a0_exp_r + a2_exp + a4) * inverse_one_a5 

54 

55 # Compute helper variables for s 

56 s1 = -6*a5*r**5*(a0*exp(a1*r)/r + a2*exp(a3*r) + a4)/(a5*r**6 + one)**2 

57 s2 = (a0*a1*exp(a1*r)/r - a0*exp(a1*r)/r**2 + a2*a3*exp(a3*r))/(a5*r**6 + one) 

58 # First derivative of SAAP divided by r with a minus sign 

59 s = - eps * (s1 + s2) / (r * sigma**2) # sigma² because of chain rule (CR) and 1/dist = 1/(r sigma) 

60 

61 # Compute helper variables for u''(r) 

62 d2u_dr2_1 = numba.float32(72.0)*a5**2*r**10*(a0*exp(a1*r)/r 

63 + a2*exp(a3*r) 

64 + a4)/(a5*r**6 + one)**3 

65 d2u_dr2_2 = -numba.float32(12.0)*a5*r**5*(a0*a1*exp(a1*r)/r 

66 - a0*exp(a1*r)/r**2 

67 + a2*a3*exp(a3*r))/(a5*r**6 + one)**2 

68 d2u_dr2_3 = -numba.float32(30.0)*a5*r**4*(a0*exp(a1*r)/r 

69 + a2*exp(a3*r) 

70 + a4)/(a5*r**6 + one)**2 

71 d2u_dr2_4 = (a0*a1**2*exp(a1*r)/r 

72 - numba.float32(2.0)*a0*a1*exp(a1*r)/r**2 

73 + numba.float32(2.0)*a0*exp(a1*r)/r**3 

74 + a2*a3**2*exp(a3*r))/(a5*r**6 + one) 

75 # Second derivative of SAAP 

76 d2u_dr2 = eps * (d2u_dr2_1 + d2u_dr2_2 + d2u_dr2_3 + d2u_dr2_4) / sigma**2 # sigma² because of double CR 

77 

78 return u, s, d2u_dr2 # u(r), -u'(r)/r, u''(r)