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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
|
.. _ExampleScripts:
==================================
Example scripts
==================================
----------------------------------
A basic fitting script
----------------------------------
The below script uses Gaussian descriptors with a neural network backend --- the Behler-Parrinello approach --- to train energies only to a training set made by the script. Note that most of the code is just generating the training data, and the training takes place in a couple of lines.
.. code-block:: python
"""Simple test of the Amp calculator, using Gaussian descriptors and neural
network model. Randomly generates data with the EMT potential in MD
simulations."""
import os
from ase import Atoms, Atom, units
import ase.io
from ase.calculators.emt import EMT
from ase.lattice.surface import fcc110
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md import VelocityVerlet
from ase.constraints import FixAtoms
from amp import Amp
from amp.descriptor.gaussian import Gaussian
from amp.model.neuralnetwork import NeuralNetwork
def generate_data(count, filename='training.traj'):
"""Generates test or training data with a simple MD simulation."""
if os.path.exists(filename):
return
traj = ase.io.Trajectory(filename, 'w')
atoms = fcc110('Pt', (2, 2, 2), vacuum=7.)
atoms.extend(Atoms([Atom('Cu', atoms[7].position + (0., 0., 2.5)),
Atom('Cu', atoms[7].position + (0., 0., 5.))]))
atoms.set_constraint(FixAtoms(indices=[0, 2]))
atoms.set_calculator(EMT())
atoms.get_potential_energy()
traj.write(atoms)
MaxwellBoltzmannDistribution(atoms, 300. * units.kB)
dyn = VelocityVerlet(atoms, dt=1. * units.fs)
for step in range(count - 1):
dyn.run(50)
traj.write(atoms)
generate_data(20)
calc = Amp(descriptor=Gaussian(),
model=NeuralNetwork(hiddenlayers=(10, 10, 10)))
calc.train(images='training.traj')
Note you can monitor the progress of the training by typing `amp-plotconvergence amp-log.txt`, which will create a file called `convergence.pdf`.
----------------------------------
A basic script with forces
----------------------------------
The below script trains both energy and forces to the same training set as above. Note this may take some time to run, which will depend upon the initial guess for the neural network parameters that is randomly generated. Try decreasing the `force_rmse` convergence parameter if you would like faster results.
.. code-block:: python
"""Simple test of the Amp calculator, using Gaussian descriptors and neural
network model. Randomly generates data with the EMT potential in MD
simulations."""
import os
from ase import Atoms, Atom, units
import ase.io
from ase.calculators.emt import EMT
from ase.lattice.surface import fcc110
from ase.md.velocitydistribution import MaxwellBoltzmannDistribution
from ase.md import VelocityVerlet
from ase.constraints import FixAtoms
from amp import Amp
from amp.descriptor.gaussian import Gaussian
from amp.model.neuralnetwork import NeuralNetwork
from amp.model import LossFunction
def generate_data(count, filename='training.traj'):
"""Generates test or training data with a simple MD simulation."""
if os.path.exists(filename):
return
traj = ase.io.Trajectory(filename, 'w')
atoms = fcc110('Pt', (2, 2, 2), vacuum=7.)
atoms.extend(Atoms([Atom('Cu', atoms[7].position + (0., 0., 2.5)),
Atom('Cu', atoms[7].position + (0., 0., 5.))]))
atoms.set_constraint(FixAtoms(indices=[0, 2]))
atoms.set_calculator(EMT())
atoms.get_potential_energy()
traj.write(atoms)
MaxwellBoltzmannDistribution(atoms, 300. * units.kB)
dyn = VelocityVerlet(atoms, dt=1. * units.fs)
for step in range(count - 1):
dyn.run(50)
traj.write(atoms)
generate_data(20)
calc = Amp(descriptor=Gaussian(),
model=NeuralNetwork(hiddenlayers=(10, 10, 10)))
calc.model.lossfunction = LossFunction(convergence={'energy_rmse': 0.02,
'force_rmse': 0.02})
calc.train(images='training.traj')
Note you can monitor the progress of the training by typing `amp-plotconvergence amp-log.txt`, which will create a file called `convergence.pdf`.
----------------------------------
Examining fingerprints
----------------------------------
With the modular nature, it's straightforward to analyze how fingerprints change with changes in images.
The below script makes an animated GIF that shows how a fingerprint about the O atom in water changes as one of the O-H bonds is stretched.
Note that most of the lines of code below are either making the atoms or making the figure; very little effort is needed to produce the fingerprints themselves---this is done in three lines.
.. code-block:: python
# Make a series of images.
import numpy as np
from ase.structure import molecule
from ase import Atoms
atoms = molecule('H2O')
atoms.rotate('y', -np.pi/2.)
atoms.set_pbc(False)
displacements = np.linspace(0.9, 8.0, 20)
vec = atoms[2].position - atoms[0].position
images = []
for displacement in displacements:
atoms = Atoms(atoms)
atoms[2].position = (atoms[0].position + vec * displacement)
images.append(atoms)
# Fingerprint using Amp.
from amp.descriptor.gaussian import Gaussian
descriptor = Gaussian()
from amp.utilities import hash_images
images = hash_images(images, ordered=True)
descriptor.calculate_fingerprints(images)
# Plot the data.
from matplotlib import pyplot
def barplot(hash, name, title):
"""Makes a barplot of the fingerprint about the O atom."""
fp = descriptor.fingerprints[hash][0]
fig, ax = pyplot.subplots()
ax.bar(range(len(fp[1])), fp[1])
ax.set_title(title)
ax.set_ylim(0., 2.)
ax.set_xlabel('fingerprint')
ax.set_ylabel('value')
fig.savefig(name)
for index, hash in enumerate(images.keys()):
barplot(hash, 'bplot-%02i.png' % index,
'%.2f$\\times$ equilibrium O-H bondlength'
% displacements[index])
# For fun, make an animated gif.
import os
filenames = ['bplot-%02i.png' % index for index in range(len(images))]
command = ('convert -delay 100 %s -loop 0 animation.gif' %
' '.join(filenames))
os.system(command)
.. image:: _static/animation.gif
:width: 600 px
:align: center
|
