1from math import radians, sin, cos
2
3import pytest
4
5from ase import Atoms
6from ase.neb import NEB
7from ase.constraints import FixAtoms
8from ase.optimize import QuasiNewton, BFGS
9from ase.visualize import view
10
11
12@pytest.mark.calculator('nwchem')
13def test_h3o2m(factory):
14    # http://jcp.aip.org/resource/1/jcpsa6/v97/i10/p7507_s1
15    doo = 2.74
16    doht = 0.957
17    doh = 0.977
18    angle = radians(104.5)
19    initial = Atoms('HOHOH',
20                    positions=[(-sin(angle) * doht, 0, cos(angle) * doht),
21                               (0., 0., 0.),
22                               (0., 0., doh),
23                               (0., 0., doo),
24                               (sin(angle) * doht, 0., doo - cos(angle) * doht)])
25    if 0:
26        view(initial)
27
28    final = Atoms('HOHOH',
29                  positions=[(- sin(angle) * doht, 0., cos(angle) * doht),
30                             (0., 0., 0.),
31                             (0., 0., doo - doh),
32                             (0., 0., doo),
33                             (sin(angle) * doht, 0., doo - cos(angle) * doht)])
34    if 0:
35        view(final)
36
37    # Make band:
38    images = [initial.copy()]
39    for i in range(3):
40        images.append(initial.copy())
41    images.append(final.copy())
42    neb = NEB(images, climb=True)
43
44    def calculator():
45        return factory.calc(
46            task='gradient',
47            theory='scf',
48            charge=-1
49        )
50
51    # Set constraints and calculator:
52    constraint = FixAtoms(indices=[1, 3])  # fix OO
53    for image in images:
54        image.calc = calculator()
55        image.set_constraint(constraint)
56
57    # Relax initial and final states:
58    with QuasiNewton(images[0]) as dyn1:
59        dyn1.run(fmax=0.10)
60    with QuasiNewton(images[-1]) as dyn2:
61        dyn2.run(fmax=0.10)
62
63    # Interpolate positions between initial and final states:
64    neb.interpolate()
65
66    for image in images:
67        print(image.get_distance(1, 2), image.get_potential_energy())
68
69    with BFGS(neb) as dyn:
70        # use better basis (e.g. aug-cc-pvdz) for NEB to converge
71        dyn.run(fmax=0.10)
72
73    for image in images:
74        print(image.get_distance(1, 2), image.get_potential_energy())
75