1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
|
def test_atoms_angle():
from ase.geometry import (get_angles, get_dihedrals,
get_distances_derivatives,
get_angles_derivatives,
get_dihedrals_derivatives)
import numpy as np
import pytest
# example: positions for H2O2 molecule
pos = np.asarray([[0.840, 0.881, 0.237], # H
[0., 0.734, -0.237], # O
[0., -0.734, -0.237], # O
[-0.840, -0.881, 0.237]]) # H
def get_numerical_derivatives(positions, mode, epsilon):
if mode == 'distance':
mode_n = 0 # get derivative for bond 0-1
elif mode == 'angle':
mode_n = 1 # get derivative for angle 0-1-2
elif mode == 'dihedral':
mode_n = 2 # get derivative for dihedral 0-1-2-3
derivs = np.zeros((2 + mode_n, 3))
for i in range(2 + mode_n):
for j in range(3):
pos = positions.copy()
pos[i, j] -= epsilon
if mode == 'distance':
minus = np.linalg.norm(pos[1] - pos[0])
elif mode == 'angle':
minus = get_angles([pos[0] - pos[1]], [pos[2] - pos[1]])
elif mode == 'dihedral':
minus = get_dihedrals([pos[1] - pos[0]], [pos[2] - pos[1]],
[pos[3] - pos[2]])
pos[i, j] += 2 * epsilon
if mode == 'distance':
plus = np.linalg.norm(pos[1] - pos[0])
elif mode == 'angle':
plus = get_angles([pos[0] - pos[1]], [pos[2] - pos[1]])
elif mode == 'dihedral':
plus = get_dihedrals([pos[1] - pos[0]], [pos[2] - pos[1]],
[pos[3] - pos[2]])
derivs[i, j] = (plus - minus) / (2 * epsilon)
return derivs
# analytical derivatives in Angstrom/Angstrom, i.e. degrees/Angstrom
distances_derivs = get_distances_derivatives([pos[1] - pos[0]])[0]
angles_derivs = get_angles_derivatives([pos[0] - pos[1]],
[pos[2] - pos[1]])[0]
dihedrals_derivs = get_dihedrals_derivatives([pos[1] - pos[0]],
[pos[2] - pos[1]],
[pos[3] - pos[2]])[0]
# numerical approximations to derivatives using finite differences
epsilon = 1e-5
num_distances_derivs = get_numerical_derivatives(pos, mode='distance',
epsilon=epsilon)
num_angles_derivs = get_numerical_derivatives(pos, mode='angle',
epsilon=epsilon)
num_dihedrals_derivs = get_numerical_derivatives(pos, mode='dihedral',
epsilon=epsilon)
# print(distances_derivs - num_distances_derivs)
# print(angles_derivs - num_angles_derivs)
# print(dihedrals_derivs - num_dihedrals_derivs)
# finite differences versus analytical results
assert num_distances_derivs == pytest.approx(distances_derivs, abs=1e-8)
assert num_angles_derivs == pytest.approx(angles_derivs, abs=1e-8)
assert num_dihedrals_derivs == pytest.approx(dihedrals_derivs, abs=1e-8)
# derivatives of multiple internal coordinates at once
assert (distances_derivs ==
get_distances_derivatives([pos[1] - pos[0],
pos[1] - pos[0]])[0]).all()
assert (angles_derivs ==
get_angles_derivatives([pos[0] - pos[1], pos[0] - pos[1]],
[pos[2] - pos[1],
pos[2] - pos[1]])[0]).all()
assert (dihedrals_derivs ==
get_dihedrals_derivatives([pos[1] - pos[0], pos[1] - pos[0]],
[pos[2] - pos[1], pos[2] - pos[1]],
[pos[3] - pos[2],
pos[3] - pos[2]])[0]).all()
|
