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

1import numpy as np 

2 

3from ase.calculators.calculator import Calculator 

4from ase.neighborlist import neighbor_list 

5 

6 

7def fcut(r, r0, r1): 

8 """ 

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

10 Ported from JuLIP.jl, https://github.com/JuliaMolSim/JuLIP.jl 

11 

12 Parameters 

13 ---------- 

14 r0 - inner cutoff radius 

15 r1 - outder cutoff radius 

16 """"" 

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

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

19 (6.0 * s**5 - 15.0 * s**4 + 10.0 * s**3)) 

20 

21 

22def fcut_d(r, r0, r1): 

23 """ 

24 Derivative of fcut() function defined above 

25 """ 

26 s = 1 - (r - r0) / (r1 - r0) 

27 return -(((0.0 < s) & (s < 1.0)) * 

28 ((30 * s**4 - 60 * s**3 + 30 * s**2) / (r1 - r0))) 

29 

30 

31class MorsePotential(Calculator): 

32 """Morse potential. 

33 

34 Default values chosen to be similar as Lennard-Jones. 

35 """ 

36 

37 implemented_properties = ['energy', 'forces'] 

38 default_parameters = {'epsilon': 1.0, 

39 'rho0': 6.0, 

40 'r0': 1.0, 

41 'rcut1': 1.9, 

42 'rcut2': 2.7} 

43 nolabel = True 

44 

45 def __init__(self, **kwargs): 

46 """ 

47 Parameters 

48 ---------- 

49 epsilon: float 

50 Absolute minimum depth, default 1.0 

51 r0: float 

52 Minimum distance, default 1.0 

53 rho0: float 

54 Exponential prefactor. The force constant in the potential minimum 

55 is k = 2 * epsilon * (rho0 / r0)**2, default 6.0 

56 """ 

57 Calculator.__init__(self, **kwargs) 

58 

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

60 system_changes=['positions', 'numbers', 'cell', 

61 'pbc', 'charges', 'magmoms']): 

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

63 epsilon = self.parameters.epsilon 

64 rho0 = self.parameters.rho0 

65 r0 = self.parameters.r0 

66 rcut1 = self.parameters.rcut1 * r0 

67 rcut2 = self.parameters.rcut2 * r0 

68 

69 forces = np.zeros((len(self.atoms), 3)) 

70 preF = - 2 * epsilon * rho0 / r0 

71 

72 i, j, d, D = neighbor_list('ijdD', atoms, rcut2) 

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

74 

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

76 fc = fcut(d, rcut1, rcut2) 

77 

78 E = epsilon * expf * (expf - 2) 

79 dE = preF * expf * (expf - 1) * dhat 

80 energy = 0.5 * (E * fc).sum() 

81 

82 F = (dE * fc + E * fcut_d(d, rcut1, rcut2) * dhat).T 

83 for dim in range(3): 

84 forces[:, dim] = np.bincount(i, weights=F[:, dim], 

85 minlength=len(atoms)) 

86 

87 self.results['energy'] = energy 

88 self.results['forces'] = forces