Coverage for /builds/hweiske/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) + ((s > 0.0) & (s < 1.0)) * ( 

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

20 ) 

21 

22 

23def fcut_d(r, r0, r1): 

24 """ 

25 Derivative of fcut() function defined above 

26 """ 

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

28 return -( 

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

30 * (30 * s**4 - 60 * s**3 + 30 * s**2) 

31 / (r1 - r0) 

32 ) 

33 

34 

35class MorsePotential(Calculator): 

36 """Morse potential. 

37 

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

39 """ 

40 

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

42 default_parameters = {'epsilon': 1.0, 

43 'rho0': 6.0, 

44 'r0': 1.0, 

45 'rcut1': 1.9, 

46 'rcut2': 2.7} 

47 nolabel = True 

48 

49 def __init__(self, **kwargs): 

50 """ 

51 Parameters 

52 ---------- 

53 epsilon: float 

54 Absolute minimum depth, default 1.0 

55 r0: float 

56 Minimum distance, default 1.0 

57 rho0: float 

58 Exponential prefactor. The force constant in the potential minimum 

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

60 """ 

61 Calculator.__init__(self, **kwargs) 

62 

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

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

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

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

67 epsilon = self.parameters.epsilon 

68 rho0 = self.parameters.rho0 

69 r0 = self.parameters.r0 

70 rcut1 = self.parameters.rcut1 * r0 

71 rcut2 = self.parameters.rcut2 * r0 

72 

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

74 preF = - 2 * epsilon * rho0 / r0 

75 

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

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

78 

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

80 fc = fcut(d, rcut1, rcut2) 

81 

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

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

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

85 

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

87 for dim in range(3): 

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

89 minlength=len(atoms)) 

90 

91 self.results['energy'] = energy 

92 self.results['forces'] = forces