1import numpy as np
2from ase.calculators.tip3p import TIP3P, rOH, angleHOH
3from ase.optimize import FIRE
4from ase.constraints import FixBondLengths
5from gpaw.utilities.watermodel import FixBondLengthsWaterModel
6from gpaw.utilities.watermodel import TIP3PWaterModel as TIP3Pc
7from ase.data import s22
8
9# check that all combinations of watermodel implementations
10# and constraints give the same results
11
12
13def test_rattle(in_tmp_dir):
14    pairs = [(3 * i + j, 3 * i + (j + 1) % 3)
15             for i in range(2)
16             for j in [0, 1, 2]]
17
18    ct = 0
19    distances = np.zeros(4)
20    for calc in [TIP3P(), TIP3Pc()]:
21        for fix in [FixBondLengths, FixBondLengthsWaterModel]:
22            atoms = s22.create_s22_system('Water_dimer')
23
24            for m in [0, 3]:
25                atoms.set_angle(m + 1, m, m + 2, angleHOH)
26                atoms.set_distance(m, m + 1, rOH, fix=0)
27                atoms.set_distance(m, m + 2, rOH, fix=0)
28                atoms.set_angle(m + 1, m, m + 2, angleHOH)
29
30            atoms.set_constraint(fix(pairs))
31
32            atoms.calc = calc
33
34            opt = FIRE(atoms, logfile='test.log')
35            opt.run(0.05)
36            p = atoms.get_positions()
37            d = np.linalg.norm(p[0] - p[3])
38            distances[ct] = d
39            ct += 1
40
41    assert (np.diff(distances) < 1e-6).all()
42