1import pytest
2from numpy.testing import assert_allclose
3from ase.cluster.icosahedron import Icosahedron
4from ase.data import atomic_numbers, atomic_masses
5from ase.optimize import LBFGS
6
7
8@pytest.fixture
9def ar_nc():
10    ar_nc = Icosahedron('Ar', noshells=2)
11    ar_nc.cell = [[300, 0, 0], [0, 300, 0], [0, 0, 300]]
12    ar_nc.pbc = True
13    return ar_nc
14
15
16@pytest.fixture
17def params():
18    params = {}
19    params['pair_style'] = 'lj/cut 8.0'
20    params['pair_coeff'] = ['1 1 0.0108102 3.345']
21    params['masses'] = ['1 {}'.format(atomic_masses[atomic_numbers['Ar']])]
22    return params
23
24
25@pytest.mark.calculator_lite
26@pytest.mark.calculator('lammpsrun')
27def test_Ar_minimize(factory, ar_nc, params):
28    with factory.calc(specorder=['Ar'], **params) as calc:
29        ar_nc.calc = calc
30
31        assert_allclose(ar_nc.get_potential_energy(), -0.468147667942117,
32                        atol=1e-4, rtol=1e-4)
33        assert_allclose(ar_nc.get_forces(),
34                        calc.calculate_numerical_forces(ar_nc),
35                        atol=1e-4, rtol=1e-4)
36
37        with LBFGS(ar_nc, force_consistent=False) as dyn:
38            dyn.run(fmax=1E-6)
39
40        assert_allclose(ar_nc.get_potential_energy(), -0.4791815886953914,
41                        atol=1e-4, rtol=1e-4)
42        assert_allclose(ar_nc.get_forces(),
43                        calc.calculate_numerical_forces(ar_nc),
44                        atol=1e-4, rtol=1e-4)
45
46
47@pytest.mark.calculator_lite
48@pytest.mark.calculator('lammpsrun')
49def test_Ar_minimize_multistep(factory, ar_nc, params):
50    ar_nc = Icosahedron('Ar', noshells=2)
51    ar_nc.cell = [[300, 0, 0], [0, 300, 0], [0, 0, 300]]
52    ar_nc.pbc = True
53
54    with factory.calc(specorder=['Ar'], **params) as calc:
55        ar_nc.calc = calc
56        F1_numer = calc.calculate_numerical_forces(ar_nc)
57
58        assert_allclose(ar_nc.get_potential_energy(), -0.468147667942117,
59                        atol=1e-4, rtol=1e-4)
60        assert_allclose(ar_nc.get_forces(), F1_numer,
61                        atol=1e-4, rtol=1e-4)
62
63        params['minimize'] = '1.0e-15 1.0e-6 2000 4000'   # add minimize
64        calc.parameters = params
65
66        # set_atoms=True to read final coordinates after minimization
67        calc.run(set_atoms=True)
68
69        # get final coordinates after minimization
70        ar_nc.set_positions(calc.atoms.positions)
71
72        assert_allclose(ar_nc.get_potential_energy(), -0.4791815887032201,
73                        atol=1e-4, rtol=1e-4)
74        assert_allclose(ar_nc.get_forces(),
75                        calc.calculate_numerical_forces(ar_nc),
76                        atol=1e-4, rtol=1e-4)
77