Coverage for /builds/debichem-team/python-ase/ase/calculators/morse.py : 100.00%
Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1import numpy as np
3from ase.calculators.calculator import Calculator
4from ase.neighborlist import neighbor_list
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
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))
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)))
31class MorsePotential(Calculator):
32 """Morse potential.
34 Default values chosen to be similar as Lennard-Jones.
35 """
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
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)
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
69 forces = np.zeros((len(self.atoms), 3))
70 preF = - 2 * epsilon * rho0 / r0
72 i, j, d, D = neighbor_list('ijdD', atoms, rcut2)
73 dhat = (D / d[:, None]).T
75 expf = np.exp(rho0 * (1.0 - d / r0))
76 fc = fcut(d, rcut1, rcut2)
78 E = epsilon * expf * (expf - 2)
79 dE = preF * expf * (expf - 1) * dhat
80 energy = 0.5 * (E * fc).sum()
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))
87 self.results['energy'] = energy
88 self.results['forces'] = forces