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
|
"""Test TIP3P forces."""
from math import cos, sin, pi
from ase import Atoms
from ase.calculators.tip3p import TIP3P, rOH, angleHOH
from ase.calculators.tip4p import TIP4P
r = rOH
a = angleHOH * pi / 180
dimer = Atoms('H2OH2O',
[(r * cos(a), 0, r * sin(a)),
(r, 0, 0),
(0, 0, 0),
(r * cos(a / 2), r * sin(a / 2), 0),
(r * cos(a / 2), -r * sin(a / 2), 0),
(0, 0, 0)])
dimer = dimer[[2, 0, 1, 5, 3, 4]]
dimer.positions[3:, 0] += 2.8
for TIPnP in [TIP3P, TIP4P]:
# put O-O distance in the cutoff range
dimer.calc = TIPnP(rc=4.0, width=2.0)
F = dimer.get_forces()
print(F)
dF = dimer.calc.calculate_numerical_forces(dimer) - F
print(dF)
assert abs(dF).max() < 2e-6
|
