遺伝的アルゴリズムを用いて金属ナノクラスターの安定構造を探索する

ASE (Atomic Simulation Environment) には遺伝的アルゴリズム (GA) を用いた安定構造探索の機能が実装されています．
遺伝的アルゴリズムを用いると，バルクやナノクラスターの安定構造を効率的に探索することができます．

公式ドキュメントのチュートリアルを参考に，金属ナノクラスターの構造探索を行ってみます．簡単のためエネルギー計算にはEMT (Effective Medium Theory) を用います．

以下の計算はJupyterLab上で行いました．

初期集団を作成する

化学量論比を指定してランダムな構造を作成します．ここでは25×25×25Aのシミュレーションボックスの中にPt₁₅Au₁₅のクラスターを10個作成しています．

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)]

view 関数を用いて，初期構造の確認します．

from ase.visualize import view

view(initial_population)

このような初期集団が生成されました．乱数で生成しているので，この結果は試行ごとに変化します．

データベースへの書き込み

作成した初期集団を外部のデータベースファイルに保存します．
ase.ga.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)

初期集団の構造緩和

ランダムに作成した構造を最適化します．

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)

遺伝的アルゴリズムの条件設定

交叉や突然変異のアルゴリズムを設定します．

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

# 原子構造の同一性判定アルゴリズム: 集団の多様性を確保するために用いる
comperator = InteratomicDistanceComparator(n_top=len(atom_numbers),
                                     pair_cor_cum_diff=0.015,
                                     pair_cor_max=0.7,
                                     dE=0.02,
                                     mic=False)

# 交叉
pairing = CutAndSplicePairing(slab, len(atom_numbers), cd)

# 突然変異: 2種類の突然変異アルゴリズムのうちからランダムに1つを選ぶ
mutations = OperationSelector([1., 1.], # 各アルゴリズムを採択する確率比
                              [MirrorMutation(cd, len(atom_numbers)),
                               RattleMutation(cd, len(atom_numbers))])
population = Population(data_connection=db,
                        population_size=population_size,
                        comparator=comperator)

遺伝的アルゴリズムの実行

遺伝的アルゴリズムを実行します．これには数分を要します．もしEMTではなくDFTを用いた場合はさらに長い時間がかかるでしょう．

from random import random

for i in range(population_size * n_generation):
    # 交叉
    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)
    
    # 突然変異
    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
        
    # 子個体の構造緩和
    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')

結果の表示

matplotlibを用いて世代ごとの最安定構造をプロットします．

%matplotlib inline

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

# Atomsオブジェクトを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))    

# 世代ごとのAtomsオブジェクトを取得 
X, Y = [], []
for i_gen in range(n_generation):
    atoms = db.get_all_relaxed_candidates_after_generation(i_gen)[0] # エネルギーごとにソートされている
    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()

参考

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