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
|
def test_neb(plt):
from ase import Atoms
from ase.constraints import FixAtoms
import ase.io
from ase.neb import NEB, NEBTools
from ase.calculators.morse import MorsePotential
from ase.optimize import BFGS, QuasiNewton
def calc():
# Common calculator for all images.
return MorsePotential()
# Create and relax initial and final states.
initial = Atoms('H7',
positions=[(0, 0, 0),
(1, 0, 0),
(0, 1, 0),
(1, 1, 0),
(0, 2, 0),
(1, 2, 0),
(0.5, 0.5, 1)],
constraint=[FixAtoms(range(6))],
calculator=calc())
dyn = QuasiNewton(initial)
dyn.run(fmax=0.01)
final = initial.copy()
final.calc = calc()
final.positions[6, 1] = 2 - initial.positions[6, 1]
dyn = QuasiNewton(final)
dyn.run(fmax=0.01)
# Run NEB without climbing image.
fmax = 0.05
nimages = 4
images = [initial]
for index in range(nimages - 2):
images += [initial.copy()]
images[-1].calc = calc()
images += [final]
neb = NEB(images)
neb.interpolate()
with BFGS(neb, trajectory='mep.traj') as dyn:
dyn.run(fmax=fmax)
# Check climbing image.
neb.climb = True
dyn.run(fmax=fmax)
# Check NEB tools.
nt_images = ase.io.read('mep.traj', index='-{:d}:'.format(nimages))
nebtools = NEBTools(nt_images)
nt_fmax = nebtools.get_fmax(climb=True)
Ef, dE = nebtools.get_barrier()
print(Ef, dE, fmax, nt_fmax)
assert nt_fmax < fmax
assert abs(Ef - 1.389) < 0.001
# Plot one band.
nebtools.plot_band()
# Plot many (ok, 2) bands.
nt_images = ase.io.read('mep.traj', index='-{:d}:'.format(2 * nimages))
nebtools = NEBTools(nt_images)
nebtools.plot_bands()
|
