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
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
|
import numpy as np
import pytest
from ase.build import molecule
from ase.utils.ff import Morse, Angle, Dihedral, VdW
from ase.calculators.ff import ForceField
from ase.optimize.precon.neighbors import get_neighbours
from ase.optimize.precon.lbfgs import PreconLBFGS
from ase.optimize.precon import FF
@pytest.fixture(scope='module')
def atoms0():
a = molecule('C60')
a.set_cell(50.0 * np.identity(3))
return a
@pytest.fixture(scope='module')
def forcefield_params(atoms0):
# force field parameters for fulleren, Z. Berkai at al.
# Energy Procedia, 74, 2015, 59-64
a = atoms0
cutoff = 1.5
morse_D = 6.1322
morse_alpha = 1.8502
morse_r0 = 1.4322
angle_k = 10.0
angle_a0 = np.deg2rad(120.0)
dihedral_k = 0.346
vdw_epsilonij = 0.0115
vdw_rminij = 3.4681
neighbor_list = [[] for _ in range(len(a))]
vdw_list = np.ones((len(a), len(a)), dtype=bool)
morses = []
angles = []
dihedrals = []
vdws = []
# create neighbor list
i_list, j_list, d_list, fixed_atoms = get_neighbours(atoms=a, r_cut=cutoff)
for i, j in zip(i_list, j_list):
neighbor_list[i].append(j)
for i in range(len(neighbor_list)):
neighbor_list[i].sort()
# create lists of morse, bending and torsion interactions
for i in range(len(a)):
for jj in range(len(neighbor_list[i])):
j = neighbor_list[i][jj]
if j > i:
morses.append(Morse(atomi=i, atomj=j, D=morse_D,
alpha=morse_alpha, r0=morse_r0))
vdw_list[i, j] = vdw_list[j, i] = False
for kk in range(jj + 1, len(neighbor_list[i])):
k = neighbor_list[i][kk]
angles.append(Angle(atomi=j, atomj=i, atomk=k, k=angle_k,
a0=angle_a0, cos=True))
vdw_list[j, k] = vdw_list[k, j] = False
for ll in range(kk + 1, len(neighbor_list[i])):
l = neighbor_list[i][ll]
dihedrals.append(Dihedral(atomi=j, atomj=i, atomk=k,
atoml=l,
k=dihedral_k))
# create list of van der Waals interactions
for i in range(len(a)):
for j in range(i + 1, len(a)):
if vdw_list[i, j]:
vdws.append(VdW(atomi=i, atomj=j, epsilonij=vdw_epsilonij,
rminij=vdw_rminij))
return dict(morses=morses, angles=angles, dihedrals=dihedrals, vdws=vdws)
@pytest.fixture
def calc(forcefield_params):
return ForceField(**forcefield_params)
@pytest.fixture
def atoms(atoms0, calc):
atoms = atoms0.copy()
atoms.calc = calc
atoms.rattle(0.05)
return atoms
ref_energy = 17.238525
@pytest.mark.slow
def test_opt_no_precon(atoms):
opt = PreconLBFGS(atoms, use_armijo=True, precon='ID')
opt.run(fmax=0.1)
e = atoms.get_potential_energy()
assert abs(e - ref_energy) < 0.01
@pytest.mark.skip('FAILS WITH PYAMG')
@pytest.mark.slow
def test_opt_with_precon(atoms, forcefield_params):
kw = dict(forcefield_params)
kw.pop('vdws')
precon = FF(**kw)
opt = PreconLBFGS(atoms, use_armijo=True, precon=precon)
opt.run(fmax=0.1)
e = atoms.get_potential_energy()
assert abs(e - ref_energy) < 0.01
|
