[PYTHON] Search for stable structures of metal nanoclusters using genetic algorithms

ASE (Atomic Simulation Environment) is equipped with a stable structure search function using a genetic algorithm (GA). Genetic algorithms can be used to efficiently search for stable structures in bulk and nanoclusters.

Let's search the structure of metal nanoclusters by referring to Tutorial of official documents. For simplicity, EMT (Effective Medium Theory) is used for energy calculation.

The following calculation was done on Jupyter Lab.

Create an initial population

Create a random structure by specifying the stoichiometric ratio. Here, 10 clusters of Pt _{15 </ sub> Au 15 </ sub> are created in a 25 × 25 × 25A simulation box.}

import numpy as np
from ase.atoms import Atoms
from ase.ga.startgenerator import StartGenerator
from ase.ga.utilities import closest_distances_generator

L = 25.
population_size = 10
atom_numbers = Atoms('Pt15Au15').get_atomic_numbers()
slab = Atoms(cell=np.array([L, L, L]))

cd = closest_distances_generator(atom_numbers=atom_numbers,
                                 ratio_of_covalent_radii=0.7)
sg = StartGenerator(slab=slab,
                    atom_numbers=atom_numbers,
                    closest_allowed_distances=cd,
                    box_to_place_in=[np.zeros(3), np.eye(3)*L/2.])

initial_population = [sg.get_new_candidate() for _ in range(population_size)]

Use the view function to check the initial structure.

from ase.visualize import view

view(initial_population)

Such an initial population was generated. Since it is generated with random numbers, this result will change from trial to trial.

Write to database

Save the created initial group in an external database file. The ʻase.ga.data` module allows you to efficiently manage large amounts of computational data.

from ase.ga.data import PrepareDB

db_file = 'ga.db'

db = PrepareDB(db_file_name=db_file,
              simulation_cell=slab,
              stoichiometry=atom_numbers,
              population_size=population_size)

for atoms in initial_population:
    db.add_unrelaxed_candidate(atoms)

Structural relaxation of the initial population

Optimize the randomly created structure.

from ase.ga.data import DataConnection
from ase.optimize import BFGS
from ase.calculators.emt import EMT

db = DataConnection(db_file)

while db.get_number_of_unrelaxed_candidates() > 0:
    atoms = db.get_an_unrelaxed_candidate()
    atoms.set_calculator(EMT())
    BFGS(atoms, trajectory=None, logfile=None).run(fmax=0.05, steps=100)
    atoms.info['key_value_pairs']['raw_score'] = -atoms.get_potential_energy()
    db.add_relaxed_step(atoms)

Setting conditions for genetic algorithms

Set the crossover and mutation algorithms.

from random import random
from ase.ga.data import DataConnection
from ase.ga.population import Population
from ase.ga.standard_comparators import InteratomicDistanceComparator
from ase.ga.cutandsplicepairing import CutAndSplicePairing
from ase.ga.offspring_creator import OperationSelector
from ase.ga.standardmutations import (MirrorMutation, RattleMutation)

mutation_probability = 0.3
n_generation = 10

#Atomic structure identity determination algorithm:Used to ensure population diversity
comperator = InteratomicDistanceComparator(n_top=len(atom_numbers),
                                     pair_cor_cum_diff=0.015,
                                     pair_cor_max=0.7,
                                     dE=0.02,
                                     mic=False)

#Crossing
pairing = CutAndSplicePairing(slab, len(atom_numbers), cd)

#mutation:Randomly choose one of two mutation algorithms
mutations = OperationSelector([1., 1.], #Probability ratio to adopt each algorithm
                              [MirrorMutation(cd, len(atom_numbers)),
                               RattleMutation(cd, len(atom_numbers))])
population = Population(data_connection=db,
                        population_size=population_size,
                        comparator=comperator)

Execution of genetic algorithm

Executes a genetic algorithm. This will take a few minutes. If you use DFT instead of EMT, it will take longer.

from random import random

for i in range(population_size * n_generation):
    #Crossing
    atoms1, atoms2 = population.get_two_candidates()
    atoms3, desc = pairing.get_new_individual([atoms1, atoms2])
    if atoms3 is None: continue
    db.add_unrelaxed_candidate(atoms3, description=desc)
    
    #mutation
    if random() < mutation_probability:
        atoms3_mut, desc = mutations.get_new_individual([atoms3])
        if atoms3_mut is not None:
            db.add_unrelaxed_step(atoms3_mut, description=desc)
            atoms3 = atoms3_mut
        
    #Structural relaxation of offspring
    atoms3.set_calculator(EMT())
    BFGS(atoms3, trajectory=None, logfile=None).run(fmax=0.05, steps=100)    
    atoms3.info['key_value_pairs']['raw_score'] = -atoms3.get_potential_energy()
    db.add_relaxed_step(atoms3)
    population.update()
    
    print(f'{i}th structure relaxation completed')

View results

Plot the most stable structure for each generation using matplotlib.

%matplotlib inline

from tempfile import NamedTemporaryFile
from ase.io import write
import matplotlib.pyplot as plt
from matplotlib.offsetbox import (OffsetImage, AnnotationBbox)

#Convert Atoms object to Matplotlib OffsetImage
def atoms_to_img(atoms):
    pngfile = NamedTemporaryFile(suffix='.png').name
    write(pngfile, atoms, scale=10, show_unit_cell=False)
    return OffsetImage(plt.imread(pngfile, format='png'), zoom=1.0)   

plt.rcParams["font.size"] = 24
fig, ax = plt.subplots(figsize=(12,8))    

#Get Atoms object for each generation
X, Y = [], []
for i_gen in range(n_generation):
    atoms = db.get_all_relaxed_candidates_after_generation(i_gen)[0] #Sorted by energy
    e = atoms.get_potential_energy()
    X += [i_gen-0.5, i_gen+0.5]
    Y += [e, e]
    if i_gen % 3 == 0:
        abb = AnnotationBbox(atoms_to_img(atoms), [i_gen+0.5, e+2.], frameon=False)
        ax.add_artist(abb)

ax.plot(X, Y, color='darkblue', linewidth=3.)
ax.set_ylim((min(Y)-1, max(Y)+5))
ax.set_title('Optimal Structure of Pt$_{15}$Au$_{15}$ by Generation')
ax.set_xlabel('GA Generation')
ax.set_ylabel('Total Energy (eV)')
plt.show()

reference

• Optimization with a Genetic Algorithm https://wiki.fysik.dtu.dk/ase/tutorials/ga/ga_optimize.html