Case Study: Si interstitial in bulk silicon¶
This tutorial will show an example on how Clinamen can be used to find the grond state configuration of the Si interstitial in bulk Silicon.
Contents
Example 1: using an empirical potential¶
The system energy is calculated using the empirical potential of Pun, Purja and Mishin [PM17]. The computational cost is therefore minimal and the code snippets can be run on a personal computer. In order to use this empirical potential, the OpenKIM API has to be installed first. Once it has been installed, the specific potential can be installed by running the command:
kim-api-collections-management install user Sim_LAMMPS_ModifiedTersoff_PurjaPunMishin_2017_Si__SM_184524061456_000
The potential can then be employed through the ASE KIM calculator, as shown in the snippets below.
Step 0. Import needed modules and classes¶
# general-purpose modules
import numpy as np
import time
# build things and gradient-based optimization with ASE
from ase.calculators.kim import KIM
import ase.build
from ase.optimize.bfgslinesearch import BFGSLineSearch as BFGS
# clinamen classes
from clinamen.evpd.core.individual import Individual
from clinamen.cmaes.population_evolver import RSPopulationEvolverGrad
from clinamen.cmaes.population_evolver import AnalizeRun
from clinamen.cmaes.evolution import TerminationCriteria, StrategyParameters
Step 1. Preparing the population founder¶
The individuals of the population are initialized from a population founder: a structure whose Cartesian coordinates represents the mean of the initial multivariate Gaussian distribution.
In our example, the founder is a 3x3x3 supercell of the cubic Si cell with the addition of an interstitial Si atom. The atom is placed in the cell so that gradient-based optimization algorithms would relax the founder structure to a local minimum of the PES.
## Create the founder
simcell = ase.build.bulk('Si', 'diamond', a=5.434, cubic=True)
# build supercell
supercell = ase.build.supercells.make_supercell(simcell, np.eye(3)*3)
# make interstitial
new_symbols = [x for x in supercell.get_chemical_symbols()]
new_symbols.append('Si')
new_cell = supercell.cell.copy()
atoms = supercell.get_scaled_positions()
# defect position
new_position = np.ones(3)*0.53
new_atoms = np.vstack((atoms, new_position))
# make the founder
founder = Individual(symbols=new_symbols,
scaled_positions=new_atoms,
pbc=True, cell=new_cell)
founder.my_name = 'Founder'
founder.defect_position = new_position # set where the defective site is
# Prepare the parameters for initializing the calculator
calc_name = 'Sim_LAMMPS_ModifiedTersoff_PurjaPunMishin_2017_Si__SM_184524061456_000'
params = {'model_name': calc_name}
# attach the calculator to the founder. Other calculators can be used in the same way
founder.set_calculator_factory(KIM, params)
# calculate founder energy, relax the founder structure and calculate the energy
initial_energy = founder.get_total_energy()
founder_clone = founder.clone()
founder_clone.optimize_structure(optimizer=BFGS) # relax the founder structure
founder_clone.write_poscar() # write POSCAR of the relaxed founder
local_minimum = founder_clone.get_total_energy()
Step 2. Initialize the evolutionary algorithm¶
In this example, the evolutionary process is carried out through the class RSPopulationEvolverGrad.
This class employs the implementation of the CMA-ES that uses a hard cutoff to select which atomic
positions should be optimizes, initializes the initial covariance matrix according with the atomic distances
from the defect positions and updates the distribution mean including information from the forces.
We therefore need to initialize these evolution parameters:
## initialization EA parameters
c_a = 0.6 # coefficient for the gradient
c_cutoff = 4 # hard cut off
c_r = 4 # coefficient for initializing the covariance matrix
step_size = 0.08 # global step-size
random_seed = 10
evolver = RSPopulationEvolverGrad(c_a, c_cutoff, c_r, founder=founder,
step_size=step_size, random_seed=random_seed)
We also need to specify other parameters, as the effective number of degrees of freedom. In general these are 3 times the number of atom in the system, but since we are using a hard cutoff, we have to explicitly inform the evolver to include these degrees of freedom only. By default, the evolver will treat the cutoff as a soft one: the initialization of the covariance matrix is affected only by the atoms within the cutoff, but all atoms participate in the evolutionary process.
# initializing number of degrees of freedom and teermination criteria
strategy_params = StrategyParameters(len(founder.chromosome.flatten()))
evolver.set_strategy_parameters(strategy_params)
terminator = TerminationCriteria(maxiter=200, smallstd=0.05)
evolver.set_termination_criteria(terminator)
# Only atoms within c_cutoff are considered during the evolution
evolver.use_reduced_population_size = True
The employed termination criteria tell the evolver to interrupt the run if the number of generations surpasses 200 or if the population energy standard deviation remains below 0.05 eV for at least 10 consecutive generations.
Step 3. Commence the evolutionary process¶
At this point the evolver has all the necessary information and we can start the evolutionary process. While the process is going on, we will print some data showing its progress.
# start the evolution
epoch = 10 # save the evolution status and the structure
# of the best individual every 10 generations
print('STARTING THE EVOLUTIONARY PROCESS...')
print('ATOMS WITHIN CUTOFF: ', evolver.atoms_within_cutoff)
start = time.time()
analizer = AnalizeRun(evolver)
analizer.initialize()
try:
for j, df in enumerate(analizer.evolve()):
pop = evolver.population
fitness = pop.individuals_fitness
index = np.argsort(fitness)
best_ind = pop[index[-1]]
count = j + 1
print(df.iloc[-1], flush=True)
if count % epoch == 0 or count == 1: # save the CMA-ES status and the
# individual with the best fitness
evolver.cmaes.save_status()
best_ind.write_poscar()
except RuntimeError as e:
print('\n!!! Error: {} !!!\n'.format(e))
end = time.time()
print(f'ELAPSED TIME: {end - start:.5f} s')
Step 4. Relax the energy of the converged solution¶
Once the evolutionary algorithm stops, we can take the best individual of the last generation and relax its structure with gradient-based methods. We can also compare the found minimum with others.
print('RELAXING CONVERGED SOLUTION...')
best_ind.optimize_structure(optimizer=BFGS) # relax converged solution
best_ind.my_name += '_OPT'
best_ind.write_poscar()
print('{:35s} {:.4f}'.format('Initial energy (eV): ',initial_energy))
print('{:35s} {:.4f}'.format('Energy local minimum (eV): ', local_minimum))
print('{:35s} {:.4f}'.format('Energy converged solution (eV): ', best_ind.get_total_energy()))
Example 2: using DFT with VASP¶
All the steps required to run this program are the same as in the previous example. The only difference consists in the use of the VASP calculator (please, consults the relative ASE documentation to learn how to set up a VASP calculator). Additionally, we also change the position of the initial configuration of the interstitial Si atom, as LDA and the empirical potential used above predict different global minima. We also employ a smaller supercell.
To run this example, the use of a computer cluster is required. The algorithm should converge to a solution between a few hours and a few days, depending on the number and capabilities of the cores employed by the VASP executable.
The whole script needed to run the job is pasted below.
import numpy as np
import os
from mpi4py import MPI
import ase
import ase.io
import ase.build
from ase.calculators.vasp import Vasp
from clinamen.evpd.core.individual import Individual
from clinamen.cmaes.population_evolver import RSPopulationEvolverGrad, AnalizeRun
from clinamen.cmaes.evolution import TerminationCriteria, StrategyParameters
from clinamen.cmaes.fitness_calculators import write_train_hdf5
my_rank = MPI.COMM_WORLD.Get_rank()
## Specify the calculator parameters
# CHANGE ACCORDINGLY TO YOUR VASP INSTALLATION
vasp_exec = '/opt/ohpc/pub/vasp/bin/vasp_gam'
vasp_cmd = 'prun ' + vasp_exec + ' 1>vasp.out 2>vasp.err'
incar_params = {'command': vasp_cmd, 'xc': 'lda', 'lwave': True,
'setups': 'minimal', 'addgrid': True, 'encut': 307.1,
'prec': 'accurate', 'ediff': 1e-3,
'ismear': 0, 'sigma': 0.01, 'lscalu': False, 'lplane': True,
'lreal': False, 'npar': 8, 'gamma': True, 'kpts': [1, 1, 1]}
## Create the founder
pos_def = np.array([0.6332521153185382, 0.6185839569765310, 0.3698019833234530])
simcell = ase.build.bulk('Si', 'diamond', a=5.434, cubic=True)
# build supercell
supercell = ase.build.supercells.make_supercell(simcell, np.eye(3)*2)
# make interstitial
new_symbols = [x for x in supercell.get_chemical_symbols()]
new_symbols.append('Si')
new_cell = supercell.cell.copy()
atoms = supercell.get_scaled_positions()
# defect position
new_position = pos_def
new_atoms = np.vstack((atoms, new_position))
# make founder
founder = Individual(symbols=new_symbols,
scaled_positions=new_atoms,
pbc=True, cell=new_cell)
founder.my_name = 'Founder'
founder.defect_position = new_position # set where the defective site is
founder.set_calculator_factory(Vasp, incar_params)
## Specify CMAES initial and termination parameters
strategy_params = StrategyParameters(len(founder.chromosome.flatten()))
terminator = TerminationCriteria(maxiter=300, smallstd=0.05)
# CMAES params
c_a = 0.30
c_cutoff = 4
c_r = 10
sigma = 0.08
runno = 10
# During the EA run we save some data points
# which can eventually be used to train a ML model
dataset = 'Xy_train.hdf5'
with open('cmaes_output.txt', 'w') as out_file:
if my_rank == 0:
out_file.write('Starting CMAES for point defect\n')
out_file.flush()
evolver = RSPopulationEvolverGrad(c_a, c_cutoff, c_r,
founder=founder, step_size=sigma,
random_seed=runno)
evolver.set_strategy_parameters(strategy_params)
evolver.set_termination_criteria(terminator)
evolver.use_reduced_population_size = True
if my_rank == 0:
out_file.write('Initial fitness: {}\n'.format(str(evolver.founder.fitness)))
evolver.founder.write_poscar()
out_file.write('ATOMS WITHIN CUTOFF: {}\n'.format(str(evolver.atoms_within_cutoff)))
out_file.flush()
if not os.path.exists('status'):
os.makedirs('status')
analizer = AnalizeRun(evolver)
analizer.initialize()
try:
for j, df in enumerate(analizer.evolve()):
if my_rank == 0:
out_file.write(str(df.iloc[-1]))
out_file.write('\n')
out_file.flush()
evolver.cmaes.save_status('status')
pop = evolver.population
fitness = pop.individuals_fitness
index = np.argsort(fitness)
pop[index[-1]].write_poscar('status')
# prepare data to be written on the HDF5 file
energies = -np.array(fitness)
forces = [x.get_forces().ravel() for x in pop]
forces = np.vstack(forces)
write_train_hdf5(dataset, pop, energies, forces, 'Si_int')
if my_rank == 0:
df.to_pickle(os.path.join('status', 'summary_df.pkl'))
except RuntimeError as e:
if my_rank == 0:
out_file.write('\n!!! Error: {}\n'.format(e))
out_file.flush()
Example 3: training a metamodel on the fly¶
The process is similar to the one outlined above, but we know need to specify the metamodel
characteristics. Also, the population evolver we are going to use must be one that supports
a metamodel: RSPopulationEvolverMetamodel.
To specify the metamodel, we need:
A metamodel object that implements the Gaussian process regressor. We do so by subclassing the abstract class
GPMetaModeland defining the methodinitialize_model().A
DescriptorsDatabaseobject that will take care to write and read materials descriptors generated on-the-fly through the use of aDescriptorsGeneratorobject.
For example, implementing a Gaussian process regressor employing a kernel made of two squared-exponential kernels we would create the class:
class DoubleRBFMetaModel(GPMetaModel):
def initialize_model(self, X_train, y_train):
mean_function = gpytorch.means.ConstantMean()
x = 0.01
kernel1_kwargs = dict(
lengthscale_constraint=gpytorch.constraints.Interval(1e2, 1e6))
kernel_func1 = gpytorch.kernels.ScaleKernel(
gpytorch.kernels.RBFKernel(**kernel1_kwargs),
outputscale=x, outputscale_constraint=gpytorch.constraints.Interval(x-1e-3, x+1e-3))
kernel2_kwargs = dict(
lengthscale_constraint=gpytorch.constraints.Interval(1e-2, 1e2))
kernel_func2 = gpytorch.kernels.ScaleKernel(
gpytorch.kernels.RBFKernel(**kernel2_kwargs),
outputscale=1-x, outputscale_constraint=gpytorch.constraints.Interval(1-x-1e-3, 1-x+1e-3))
kernel_function = kernel_func1 + kernel_func2
noise_constraint = gpytorch.constraints.Interval(1e-5, 1e-1)
likelihood_kwargs = dict(noise_covar=1e-3, noise_constraint=noise_constraint)
likelihood = gpytorch.likelihoods.GaussianLikelihood(**likelihood_kwargs)
model = ExactGPModel(X_train, y_train, mean_function,
kernel_function, likelihood)
marginal_log_likelihood = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood,
model)
optimizer = torch.optim.Adam(model.parameters(), lr=1, amsgrad=True)
self._mean_function = mean_function
self._kernel_function = kernel_function
self._likelihood = likelihood
self._model = model
self._optimizer = optimizer
self._mll = marginal_log_likelihood
The way the Gaussian process is built follows the implementation in gpytorch, we strongly suggest the users of Clinamen to read gpytorch’s documentation to learn how to build Gaussian process regressor with that library.
ExactGPModel is a conventient wrapper class for initializing and exact Gaussian process.
Note
It is important that the initialize_model() sets the private attributes as shown in the last 6 lines of the code snippet above.
As the name suggests, a DescriptorsGenerator will take care to generate the appropiate descriptors.
For example, if we would like to use the many-body tensor representations, we would write something like this
(please consult the DScribe library for details about descriptors specifications):
descr_kwargs = {'species': ['Si'], 'normalization': 'l2_each',
'k1': {
"geometry": {"function": "atomic_number"},
"grid": {"min": 0, "max": 20, "sigma": 0.2, "n": 100},
},
'k2': {
"geometry": {"function": "inverse_distance"},
"grid": {"min": 0, "max": 1.0, "sigma": 0.02, "n": 200},
"weighting": {"function": "exponential", "scale": 1.0, "cutoff": 1e-3},
},
'k3': {
"geometry": {"function": "cosine"},
"grid": {"min": -1.0, "max": 1.0, "sigma": 0.02, "n": 200},
"weighting": {"function": "exponential", "scale": 1.0, "cutoff": 1e-3},
},
'flatten': True,
'sparse': False,
'periodic': True}
descriptors_generator = DescriptorsGenerator('MBTR', descr_kwargs)
We would then feed the descriptors_generator object to a DescriptorsDatabase instance:
database = DescriptorsDatabase('mbtr_descriptors.hdf5',
descriptors_generator)
We can finally initialize the RSPopulationEvolverMetamodel object. Optionally, we can specify a preprocessing_pipeline which will transform the descriptors before using them to train the metamodel:
scaler_kwargs = dict(with_std=False)
pca_kwargs = dict(n_components=0.95)
scaler = StandardScaler(**scaler_kwargs)
pca = PCA(**pca_kwargs)
steps = [
('pca', pca),
('scaler', scaler)
]
pipeline = Pipeline(steps)
metamodel = DoubleRBFMetaModel(database, preprocessing_pipeline=pipeline,
std_value=5e-3)
In this example, the descriptors dimensionality if reduced by using principal components analysis and then the resulting features are centered. This pipeline is built as a sklearn.pipeline.Pipeline using the from sklearn.decomposition.PCA and
sklearn.preprocessing.StandardScaler objects.
std_value=5e-3 indicates that the metamodel’s prediction will be used only when their standard error will be smalled than 5e-3 eV/atom.
The whole script needed to run the job is pasted below.
import numpy as np
import os
import ase.io
from ase.calculators.vasp import Vasp
import torch
import gpytorch
from sklearn.pipeline import Pipeline
from sklearn.decomposition import PCA
from sklearn.preprocessing import StandardScaler
from clinamen.evpd.core.individual import Individual
from clinamen.cmaes.population_evolver import RSPopulationEvolverMetamodel, AnalizeRun
from clinamen.cmaes.evolution import TerminationCriteria, StrategyParameters
from clinamen.descriptors.descriptors import DescriptorsGenerator
from clinamen.descriptors.descriptors import DescriptorsDatabase
from clinamen.metamodel.metamodel import GPMetaModel
from clinamen.metamodel.gp_models import ExactGPModel
class DoubleRBFMetaModel(GPMetaModel):
def initialize_model(self, X_train, y_train):
mean_function = gpytorch.means.ConstantMean()
x = 0.01
kernel1_kwargs = dict(
lengthscale_constraint=gpytorch.constraints.Interval(1e2, 1e6))
kernel_func1 = gpytorch.kernels.ScaleKernel(
gpytorch.kernels.RBFKernel(**kernel1_kwargs),
outputscale=x, outputscale_constraint=gpytorch.constraints.Interval(x-1e-3, x+1e-3))
kernel2_kwargs = dict(
lengthscale_constraint=gpytorch.constraints.Interval(1e-2, 1e2))
kernel_func2 = gpytorch.kernels.ScaleKernel(
gpytorch.kernels.RBFKernel(**kernel2_kwargs),
outputscale=1-x, outputscale_constraint=gpytorch.constraints.Interval(1-x-1e-3, 1-x+1e-3))
kernel_function = kernel_func1 + kernel_func2
noise_constraint = gpytorch.constraints.Interval(1e-5, 1e-1)
likelihood_kwargs = dict(noise_covar=1e-3, noise_constraint=noise_constraint)
likelihood = gpytorch.likelihoods.GaussianLikelihood(**likelihood_kwargs)
model = ExactGPModel(X_train, y_train, mean_function,
kernel_function, likelihood)
marginal_log_likelihood = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood,
model)
optimizer = torch.optim.Adam(model.parameters(), lr=1, amsgrad=True)
self._mean_function = mean_function
self._kernel_function = kernel_function
self._likelihood = likelihood
self._model = model
self._optimizer = optimizer
self._mll = marginal_log_likelihood
pos_def = np.array([0.6332521153185382, 0.6185839569765310, 0.3698019833234530])
# prepare system
poscar = ase.io.read('POSCAR_tetra')
# prepare calculator
vasp_exec = '/opt/ohpc/pub/vasp/bin/vasp_gam'
vasp_cmd = 'prun ' + vasp_exec + ' 1>vasp.out 2>vasp.err'
incar_params = {'command': vasp_cmd, 'xc': 'lda', 'lwave': True,
'setups': 'minimal', 'addgrid': True, 'encut': 307.1,
'prec': 'accurate', 'ediff': 1e-3,
'ismear': 0, 'sigma': 0.01, 'lscalu': False, 'lplane': True,
'lreal': False, 'npar': 8, 'gamma': True, 'kpts': [1, 1, 1]}
# make founder
founder = Individual.make_individual_from_ase_atoms(poscar)
founder.set_calculator_factory(Vasp, incar_params)
founder.defect_position = pos_def
strategy_params = StrategyParameters(len(founder.chromosome.flatten()))
terminator = TerminationCriteria(maxiter=3000, smallstd=0.01)
# CMAES params
c_a = 0.30
nn_dist = 4
c_r = 10
sigma = 0.08
runno = 10
# prepare descriptors database
symbols = founder.get_chemical_symbols()
descr_kwargs = {'species': symbols, 'normalization': 'l2_each',
'k1': {
"geometry": {"function": "atomic_number"},
"grid": {"min": 0, "max": 20, "sigma": 0.2, "n": 100},
},
'k2': {
"geometry": {"function": "inverse_distance"},
"grid": {"min": 0, "max": 1.0, "sigma": 0.02, "n": 200},
"weighting": {"function": "exponential", "scale": 1.0, "cutoff": 1e-3},
},
'k3': {
"geometry": {"function": "cosine"},
"grid": {"min": -1.0, "max": 1.0, "sigma": 0.02, "n": 200},
"weighting": {"function": "exponential", "scale": 1.0, "cutoff": 1e-3},
},
'flatten': True,
'sparse': False,
'periodic': True}
descriptors_generator = DescriptorsGenerator('MBTR', descr_kwargs)
database = DescriptorsDatabase('mbtr_descriptors.hdf5',
descriptors_generator)
# prepare metamodel
scaler_kwargs = dict(with_std=False)
pca_kwargs = dict(n_components=0.95)
scaler = StandardScaler(**scaler_kwargs)
pca = PCA(**pca_kwargs)
steps = [
('pca', pca),
('scaler', scaler)
]
pipeline = Pipeline(steps)
metamodel = DoubleRBFMetaModel(database, preprocessing_pipeline=pipeline,
std_value=5e-3)
dataset = 'Xy_train.hdf5'
with open('cmaes_output.txt', 'w') as out_file:
out_file.write('Starting CMAES for point defect\n')
out_file.flush()
evolver = RSPopulationEvolverMetamodel(metamodel, dataset, nn_dist, c_r, founder=founder,
step_size=sigma, covariance=1,
random_seed=runno,
train_kwargs={'epochs': 30000})
evolver.set_strategy_parameters(strategy_params)
evolver.set_termination_criteria(terminator)
evolver.use_reduced_population_size = True
out_file.write('Initial fitness: {}\n'.format(str(evolver.founder.fitness)))
evolver.founder.write_poscar()
out_file.write('ATOMS WITHIN CUTOFF: {}\n'.format(str(evolver.atoms_within_cutoff)))
out_file.flush()
if not os.path.exists('status'):
os.makedirs('status')
analizer = AnalizeRun(evolver)
analizer.initialize()
try:
for j, df in enumerate(analizer.evolve()):
out_file.write(str(df.iloc[-1]))
out_file.write('\n')
out_file.flush()
evolver.cmaes.save_status('status')
metamodel.save_model('model_state_' + str(j) + '.pth')
pop = evolver.population
fitness = pop.individuals_fitness
energies = -np.array(fitness)
index = np.argsort(fitness)
pop[index[-1]].write_poscar('status')
df.to_pickle(os.path.join('status', 'summary_df.gzip'))
except RuntimeError as e:
out_file.write('\n!!! Error: {}\n'.format(e))
out_file.flush()
~