Commit f1210e0d authored by Malthe Kjær Bisbo's avatar Malthe Kjær Bisbo
Browse files

Initialized sphinx documentation with a few turorials + added hessian to gauss_kernel

parent b9ff7151
import numpy as np
from ase.calculators.dftb import Dftb
from ase.io import read
from candidate_operations.candidate_generation import CandidateGenerator, StartGenerator
from candidate_operations.basic_mutations import RattleMutation, PermutationMutation
from gofee import GOFEE
### Define calculator ###
calc = Dftb(label='TiO2_surface',
Hamiltonian_SCC='No',
Hamiltonian_MaxAngularMomentum_='',
Hamiltonian_MaxAngularMomentum_Ti='"d"',
Hamiltonian_MaxAngularMomentum_O='"p"',
Hamiltonian_Charge='0.000000',
Hamiltonian_Filling ='Fermi {',
Hamiltonian_Filling_empty= 'Temperature [Kelvin] = 0.000000',
kpts=(2,1,1))
### Set up StartGenerator and mutations ###
# read slab
slab = read('TiO2_slab.traj', index='0')
# Stoichiometry of atoms to be placed
stoichiometry = 5*[22]+10*[8]
# Box in which to place atoms
v = slab.get_cell()
v[2,2] = 2.5
p0 = np.array((0.0,0.,8.))
box = [p0, v]
# initialize startgenerator
sg = StartGenerator(slab, stoichiometry, box)
# initialize rattle and permutation mutations
n_to_optimize = len(stoichiometry)
permutation = PermutationMutation(n_to_optimize, Npermute=2)
rattle = RattleMutation(n_to_optimize, Nrattle=3, rattle_range=4)
candidate_generator = CandidateGenerator([0.2, 0.2, 0.6],
[sg, permutation, rattle])
### Initialize and run search ###
search = GOFEE(calc=calc,
startgenerator=sg,
candidate_generator=candidate_generator,
max_steps=100,
population_size=5)
search.run()
\ No newline at end of file
================================================
Searching for the TiO2(001)-(1x4) reconstruction
================================================
For this tutorial we will use the dftb-calculator with
the tiorg parameters.
This tutorial is very similar to the previous one for TiO clusters,
:ref:`searching for TiO clusters <searching-for-TiO-clusters>`. It is
recomended that you do that one before the present one, as it is more
detailed.
The major difference in the present tutorial is that the template will
now not be empty, but contain a number of atoms fixed at bulk positions.
The template is defined in the file :download:`TiO2_slab.traj`. The
following code :download:`TiO2.py` is used to carry out the search:
.. literalinclude:: TiO2.py
If ASE, GPAW and dftb are set up and sourced propperly, you can run
the code as::
mpiexec --mca mpi_warn_on_fork 0 gpaw-python run_search.py
Setting up the system - atoms in template
=========================================
In this case the *template* contains a number of fixed atoms representing the
slap, on top of which we want to optimize a number of atoms given by
*stoichiometry*. The final thing we need to initialize the :class:`StartGenerator`
, used for generation initial structures, is the *box* within which the
:class:`StartGenerator` places atoms randomly.
In this case we choose a box=[p0, v] of height 2.5 starting at p0=(0,0,8), which
is slightly above the slab atoms.
To initialize the startgenerator, we first read in the template::
from ase.io import read
slab = read('TiO2_slab.traj', index='0')
then define the stoichiometry of atoms to be optimized on top of the slab,
in the form of a list of atomic numbers::
stoichiometry = 5*[22]+10*[8]
Then define the *box* within which the :class:`StartGenerator` places atoms randomly::
import numpy as np
v = slab.get_cell()
v[2,2] = 2.5
p0 = np.array((0.0,0.,8.))
box = [p0, v]
Finally the :class:`StartGenerator` can be initialized::
from candidate_operations.candidate_generation import StartGenerator
sg = StartGenerator(slab, stoichiometry, box)
\ No newline at end of file
# Creates: structures.traj
import numpy as np
from ase import Atoms
from ase.calculators.dftb import Dftb
from candidate_operations.candidate_generation import CandidateGenerator, StartGenerator
from candidate_operations.basic_mutations import RattleMutation
from gofee import GOFEE
### Define calculator ###
calc = Dftb(label='TiO2_surface',
Hamiltonian_SCC='No',
Hamiltonian_MaxAngularMomentum_='',
Hamiltonian_MaxAngularMomentum_Ti='"d"',
Hamiltonian_MaxAngularMomentum_O='"p"',
Hamiltonian_Charge='0.000000',
Hamiltonian_Filling ='Fermi {',
Hamiltonian_Filling_empty= 'Temperature [Kelvin] = 0.000000',
kpts=(1,1,1))
### Set up system ###
# make empty cell
template = Atoms('',
cell=[20,20,20],
pbc=[0, 0, 0])
# Stoichiometry of atoms to be placed
stoichiometry = 5*[22]+10*[8]
# Box in which to place atoms randomly
v = 4*np.eye(3)
p0 = np.array((8.0, 8.0, 8.0))
box = [p0, v]
# initialize startgenerator (used to generate initial structures)
sg = StartGenerator(template, stoichiometry, box)
### Set up candidate generation operations ###
# initialize rattle mutation
n_to_optimize = len(stoichiometry)
rattle = RattleMutation(n_to_optimize, Nrattle=3, rattle_range=4)
candidate_generator = CandidateGenerator(probabilities=[0.2, 0.8],
operations=[sg, rattle])
### Initialize and run search ###
search = GOFEE(calc=calc,
startgenerator=sg,
candidate_generator=candidate_generator,
max_steps=100,
population_size=5)
search.run()
\ No newline at end of file
.. _searching-for-TiO-clusters:
==========================
Searching for TiO clusters
==========================
For this tutorial we will use the dftb-calculator with
the tiorg parameters.
In this tutorial we carry out a search for titanium-oxide clusters using
dftb to evaluate energies and forses of the structures.
The following script :download:`Ti5O10.py` is used to carry out the search (the indivitual elements are
explainted further below):
.. literalinclude:: Ti5O10.py
If ASE, GPAW and dftb are set up and sourced propperly, you can run
the code as::
mpiexec --mca mpi_warn_on_fork 0 gpaw-python run_search.py
What follows is a description of the python script above.
Setting up the system
=====================
An important prerequisite for starting a search is to set up the system.
This is done by defining a template and a stoichiometry of the atoms to
optimize.
The *template* is an :class:`Atoms` object, either describing an empty cell or
a cell containing for example a slab of atoms. For most purposes, the atoms
in the template shold be fixed using :class:`ase.constraints.FixAtoms`
constraint, as the template atoms are kept fixed during mutation operation,
but will take part in surrogate-relaxation if not fixed.
In this example the template is taken to be an empty 20Åx20Åx20Å cell, since
we considder isolated TiO-clusters. The code to generate the template is::
from ase import Atoms
template = Atoms('',
cell=[20,20,20],
pbc=[0, 0, 0])
The *stoichiometry* of atoms to optimize is a list of atomic numbers. In this
case 5 titanium (atomic nymber 22) and 10 oxygen (atomic number 8) atoms::
stoichiometry = 5*[22]+10*[8]
Startgenerater - for making initial structures
==============================================
To initialize the search, initial structures need to be generated. This is
carried out using the :class:`StartGenerator`, which in addition to the
*template* and *stoichiometry* defined above, need a *box* in which to randomly
place the atoms defined in the *stoichiometry*.
The *box* is naturally defined by a point *p0* and three spanning vectors going
out from that point. These are defined bu the 3x3 matrix *v* in the example.
In the example a 20Åx20Åx20Å square box in the center of the cell is used::
import numpy as np
v = 4*np.eye(3)
p0 = np.array((8.0, 8.0, 8.0))
box = [p0, v]
The *startgenerator* can then be initialized with the code::
from candidate_operations.candidate_generation import StartGenerator
sg = StartGenerator(template, stoichiometry, box)
CandidateGenerator
==================
In GOFEE, the configurational space is explored by generation new candidate structures.
New candidates can be either completely random structures made using the *startgenerator*
or they can be the result of applying mutation operations to a population of some of the
best structures visited during the search. Examples of mutaion operations are the
:class:'RattleMutation', which randomly shifts some of the atoms and the
:class:`PermutaionMutation` which randomly permutes some atoms of different type.
The rattle mutation in the example, which rattles on average Natoms=3 atom a maximum distance of
rattle_range=4Å, is initialized as::
from candidate_operations.basic_mutations import RattleMutation
n_to_optimize = len(stoichiometry)
rattle = RattleMutation(n_to_optimize, Nrattle=3, rattle_range=4)
Given some of the above described operations. e.g. a :class:`StartGenerator`
and a :class:'RattleMutation', one can initialize a :class:`CandidateGenerator`,
which handles the generation of new candidates by applying the supplied
*operations* with probability specified in the *probabilities* list.
A CandidateGenerator which uses the startgenerator *sg* with 20% probability and
the rattle operation *rattle* with 80% probability, is initialized as follows::
from candidate_operations.candidate_generation import CandidateGenerator
candidate_generator = CandidateGenerator(probabilities=[0.2, 0.8],
operations=[sg, rattle])
Initialize and run GOFEE
========================
With all the above objects defined, we are ready to initialize and run GOFEE.
To run the search for 100 iterations with a population size of 5, use::
from gofee import GOFEE
search = GOFEE(calc=calc,
startgenerator=sg,
candidate_generator=candidate_generator,
max_steps=100,
population_size=5)
search.run()
This tutorial relies on many default settings of GOFEE, which could be changed.
To see how these settings are changed, have a look at the other tutorials.
========
Tutorial
========
This is the turorial for the GOFEE structure search.
.. toctree::
:maxdepth: 1
tio_clusters/tio_clusters
tio2_reconstruction/tio2_reconstruction
modifying_surrogate_model/modifying_surrogate_model
\ No newline at end of file
......@@ -98,12 +98,12 @@ class GOFEE():
kappa=2,
max_steps=200,
Ninit=10,
max_relax_dist=3.5,
max_relax_dist=4,
Ncandidates=30,
population_size=5,
population_size=10,
dualpoint=True,
min_certainty=0.7,
restart=None):
restart='restart.pickl'):
if structures is None:
assert startgenerator is not None
......
......@@ -173,9 +173,43 @@ class gauss_kernel(kernel):
Y = self.apply_eta(Y)
dd_df = 2*(X - Y)
dk_df = np.multiply(dK_dd.T, dd_df)
dk_df = np.multiply(dK_dd, dd_df)
return dk_df
def kernel_hessian(self, X,Y):
""" Jacobian of the kernel with respect to X
"""
X = self.apply_eta(X)
Y = self.apply_eta(Y)
if np.ndim(X) == 1:
X = X.reshape((1,-1))
if np.ndim(Y) == 1:
Y = Y.reshape((1,-1))
d = cdist(X / self.length_scale,
Y / self.length_scale, metric='sqeuclidean')
u = 1/(2*self.length_scale**2)
K = self.amplitude * np.exp(-0.5 * d)
X = self.apply_eta(X)
Y = self.apply_eta(Y)
N_X, Nf = X.shape
N_Y = Y.shape[0]
X_rep = np.tile(X, (N_Y,1,1)).swapaxes(0,1)
Y_rep = np.tile(Y, (N_X,1,1))
dd_dX = -2*(Y_rep-X_rep) # shape: (N_X, N_Y, Nf)
dd_dY = -dd_dX
dd_dX_dd_dY = np.einsum('nmi,nmj->nmij',dd_dX, dd_dY) # shape: (N_X, N_Y, Nf, Nf)
if self.Nsplit_eta is not None:
diag = np.diag(self.Nsplit_eta*[1]+(Nf-self.Nsplit_eta)*[self.eta**2])
else:
diag = np.identity(Nf)
d2d_dXdY = np.tile(-2*diag, (N_X, N_Y,1,1)) # shape: (N_X, N_Y, Nf, Nf)
hess = -u*K.reshape(N_X,N_Y,1,1) * (u*dd_dX_dd_dY - d2d_dXdY)
return hess
@property
def theta(self):
"""Returns the log-transformed hyperparameters of the kernel.
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment