import ase
from ase import Atoms
from ase.io import write
from ase.neb import NEB,interpolate
from ase.units import Bohr,Rydberg,kJ,kB,fs,Hartree,mol,kcal
from ase.geometry.analysis import Analysis
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import math,pickle
import matlanticgrrm.atomslist as ma
import pickle
from scipy.signal import find_peaks
import itertools
import numpy as np
from .geometry import _calc_max_distance, Geometries
[ドキュメント]class ListNEB():
def __init__(self,imagess,k=0.1,climb=False,parallel=False,remove_rotation_and_translation=False,
world=None,method='aseneb',allow_shared_calculator=False,precon=None):
self.imagess = imagess
self.neb = [NEB(images, k=k, climb=climb, parallel=parallel,
remove_rotation_and_translation=remove_rotation_and_translation,
world=world, method=method,
allow_shared_calculator=allow_shared_calculator,
precon=precon) if images else None for images in self.imagess]
for images in imagess:
if images:
self.n_images = len(images)
def __iter__(self):
"""Optimizerを使用するためには__iter__を実装する必要がある"""
return iter(self.neb)
def __len__(self):
"""Optimizerを使用するためには__len__を実装する必要がある"""
return len(self.imagess)
def __getitem__(self, item):
return self.neb[item]
[ドキュメント]def calc_image_numbers(ini, fin, indices=None, divide_per=0.1, min_images=8):
"""始状態と終状態からNEB計算に最適なimageの数を算出する
Paramaters:
ini: Atoms
始状態のAtomsオブジェクト
fin: Atoms
終状態のAtomsオブジェクト
indices: List of integer
着目する原子のインデックス番号
divide_per: float
最も大きく動いた原子の移動距離をこの距離で割ってイメージ数を定める（Å）
min_images: int
最小のイメージ数を定める
"""
number_of_images = min_images
maxd = _calc_max_distance(ini, fin, indices)
inumber = int((maxd / divide_per) // 2) # 奇数，偶数の場合がある
inumber = inumber if inumber%2==0 else inumber+1 # 奇数の場合は1を足して偶数にする
if inumber > min_images:
number_of_images = inumber
return number_of_images
[ドキュメント]def make_images(eq_list,connections,n_images=None, name="ReactionsSummary",
indices=None, divide_per=0.1, min_images=8,mic=True):
"""NEB計算のImagesを作成する.
| 不要な反応のImageは作成しない. 例えば,同じini,finの組み合わせを持つ反応は1つしかImagesを作らない.
| 作成されなかったImageの理由をHTMLファイル,またはlogファイルに出力する.
| eq_listをGeometryで与えた場合HTMLファイル,Atomsリストで与えた場合logファイルが作成される.
Parameters:
eq_list: list of Atoms or matlanticgrrm.geometry.Geometry
| EQ構造
| Atomsのリストの場合: 同じCONNECTIONの反応のimagesは作らない.
| matlanticgrrm.geometry.Geometryの場合: 同じ構造と判断された反応のimagesは作らない(構造のグループ化)
connections: list of tuple of interger
| CONNECTIONのリスト
| ex) [(0,1),(2,4),(3,5)]
n_images: interger or list or None
| Imageの数.
| intの場合: 全てintで指定したImage数で作成
| listの場合:
| Noneの場合: iniとfinの構造の差から最適なImage数を算出し作成.
indices:
| indicesに指定した原子の移動距離から適切なimages数を決定する.
| Noneの場合，全ての原子の移動距離を見る
divide_per: float
最も大きく動いた原子の移動距離をこの距離で割ってイメージ数を定める（Å）
min_images: int
最小のイメージ数を定める
"""
if type(eq_list) == list:
#計算するTS番号と移動距離の取得&Log出力
indexes = make_nebmakelog(eq_list,connections)
if type(n_images) == int:
n_images = [n_images if i in indexes else None for i in range(len(eq_list))]
elif type(n_images) == list:
pass # なにもしない(n_images=n_images)
elif n_images is None:
n_images = [calc_image_numbers(eq_list[ini],eq_list[fin],indices,divide_per,min_images)
if type(ini)!=str and type(fin)!=str else None
for ini,fin in connections]
else:
raise TypeError("n_imagesはint,list,Noneのいずれか")
elif type(eq_list) == Geometries:
geometry = eq_list
eq_list = geometry.atoms_list
#計算するTS番号と移動距離の取得&HTML出力
indexes_and_distance = geometry.get_valid_reaction_and_summary(connections,name+".html")
indexes = indexes_and_distance.keys()
distance = indexes_and_distance.values()
if type(n_images) == int:
n_images = [n_images if i in indexes else None for i in range(len(geometry))]
elif type(n_images) == list:
pass # なにもしない(n_images=n_images)
elif n_images is None:
indices = geometry.target0+geometry.target1
n_images = [calc_image_numbers(eq_list[ini],eq_list[fin],indices,divide_per,min_images)
if type(ini)!=str and type(fin)!=str else None
for ini,fin in connections]
else:
raise TypeError("n_imagesはint,list,Noneのいずれか")
else:
raise TypeError("eq_listはAtomsのリスト,matlanticgrrm.geometry.Geometryのいずれか")
return _make_images(eq_list,connections,indexes,n_images,mic)
def _make_images(eq_list:list,connections,indexes,n_images:list,mic=True):
imagess = []
for i,((ini_idx,fin_idx),n_image) in enumerate(zip(connections,n_images)):
if i in indexes:
images=[eq_list[ini_idx].copy()]+[eq_list[ini_idx].copy() for _ in range(n_image)]+[eq_list[fin_idx].copy()]
interpolate(images,mic=mic)
else:
images=None
imagess.append(images)
return imagess
def make_nebmakelog(eq,connections,ts=None,same_react=False,same_connect=False,log="nebmake.log"):
connect_copy = [] #重複したini,finの組み合わせがないか探すため
indexes = []
text = ""
for i,(ini,fin) in enumerate(connections):
if type(ini)!= int or type(fin)!= int:
text += f"i={i} ({ini},{fin})\n"
elif ini == fin:
if same_connect:
text += f"i={i} imageを作成({ini},{fin})\n"
indexes.append(i)
else:
text += f"i={i} iniとfinが同じ({ini},{fin})\n"
elif (ini,fin) in connect_copy:
if same_react or ts:
text += f"i={i} imageを作成({ini},{fin})\n"
indexes.append(i)
else:
index = connect_copy.index((ini,fin))
text += f"i={i} i={index}と同じ組み合わせ({ini},{fin})\n"
elif (fin,ini) in connect_copy:
if same_react or ts:
text += f"i={i} imageを作成({ini},{fin})\n"
indexes.append(i)
else:
index = connect_copy.index((fin,ini))
text += f"i={i} i={index}と同じ組み合わせ({ini},{fin})\n"
elif not eq[ini]:
text += f"i={i} EQ{ini}の構造がない({ini},{fin})\n"
elif not eq[fin]:
text += f"i={i} EQ{fin}の構造がない({ini},{fin})\n"
else:
text += f"i={i} imageを作成({ini},{fin})\n"
indexes.append(i)
connect_copy.append((ini,fin))
if log == "-":
print(text)
elif log is not None:
with open(log,"w") as f:
f.write(text)
return indexes
[ドキュメント]def neb_graphs(imagess,html=None,traj=None,pkl=None):
"""NEBのエネルギーダイアグラムを作成
Parameters:
imagess: list of images (2D list of Atoms)
imagesのリスト(Atomsの2次元リスト)
html: str or path object
| htmlのパス
| htmlにグラフを作成する.Noneの場合は標準出力(非推奨).
traj: str or path object
| trajのパス
| TS構造のtrajファイルを作成する. TS構造の存在しない要素は空のAtomsオブジェクト
pkl: str or path object
TS構造のindex番号, TS構造が存在しない要素はNone
"""
energiess = [ma.get_potential_energy(images) if images else None for images in imagess]
for images_e in energiess:
if images_e:
for image_e in images_e:
image_e#*mol/kJ
sample_n = len(energiess)
cols = 3 #3列
rows = math.ceil(sample_n/cols)
subplot_titles = [f"TS{i}" for i in range(sample_n)]
"""データをcol個ずつに分割する"""
div_sample = []
for i in range(0, sample_n, cols):
div_sample.append(energiess[i: i+cols])
def add_fig(fig,x,y,row,col,color):
fig.add_trace(go.Scatter(x=x, y=y,
line=dict(width=2,color=color)),
row=row, col=col)
fig = make_subplots(
rows=rows, cols=cols,
subplot_titles=subplot_titles)
idx = []
for i, row_data in enumerate(div_sample,start=1):
for j, data in enumerate(row_data,start=1):
if data:
x = [i for i in range(len(data))]
y = [(i-data[0])*mol/kJ for i in data]
max_y = max(y)
min_y = min(y)
add_fig(fig,x,y,i,j,'black')
if max_y-min_y > 0.1: # エネルギーの幅が0.1以下の時はグラフがガタついている場合はほとんどなので無視
peaks_i, _ = find_peaks(y)
peaks_eneegies = [(i,y[i]) for i in peaks_i]
ts_energy = sorted(peaks_eneegies, key=lambda x:x[1])
add_fig(fig,[ts_energy[-1][0]],[ts_energy[-1][1]],i,j,'red')
idx.append(ts_energy[0][0])
else:
idx.append(None)
else:
idx.append(None)
fig.update_layout(height=300*rows, width=1000,
title_text="エネルギーダイアグラム")
if pkl:
with open(pkl,"wb") as f:
pickle.dump(idx,f)
if traj:
atoms_list = [images[n] if n is not None else Atoms() for n,images in zip(idx,imagess)]
write(traj,atoms_list)
if html:
fig.write_html(html)
else:
fig.show()
return idx
# def neb_graph(image,calculator=None):
# if calculator:
# ma.set_calculator(calculator)
def make_connections(eq_list,target:list,n=2):
"""EQ構造のAtomsリストから反応しそうなCONNECTIONを返す
Parameter:
eq_list: list of Atoms
EQ構造のAtomsリスト
n:
nを大きくするとより多くのCONNECTIONを返す(その分反応しづらいCONNECTIONが含まれる)
"""
eq_list= _extract_target_atoms(eq_list,target) # targetに指定した原子のみを抽出する
similes = ma.get_smiles(eq_list)
group, bond_dict = _atoms_matrix_dict(eq_list)
conv = []
for (a_key,a_mat),(b_key,b_mat) in itertools.combinations(bond_dict.items(), 2):
dif = (a_mat-b_mat).data # 構造Aと構造Bの結合の差を計算
if np.count_nonzero(dif < 0) <= n and np.count_nonzero(dif > 0) <= 2:
"""反応しそうな組み合わせ"""
if similes[group[a_key][0]] != similes[group[b_key][0]]:
"""SMILES的に違う構造の時"""
min_distance_conbination = None
distance = float('inf')
for ini,fin in itertools.product(group[a_key], group[b_key]):
"""最も動く距離の短いものを算出"""
ini_structure = eq_list[ini]
fin_structure = eq_list[fin]
d = ase.geometry.distance(ini_structure, fin_structure)
if d < distance:
distance = d
min_distance_conbination = (ini,fin)
conv.append(min_distance_conbination)
return conv
def _atoms_matrix_dict(atoms_list):
"""ユニークなbondの行列を作成する
bond_dict = {hashableな行列:adjacency_matrix}
group = {hashableな行列:[idx]}
"""
group = {}
bond_dict = {}
for i, atoms in enumerate(atoms_list):
ana = Analysis(atoms)
bond_mat = ana.adjacency_matrix[0]
hashable_mat = tuple(bond_mat.items())
if not hashable_mat in group.keys():
bond_dict[hashable_mat] = bond_mat
group[hashable_mat] = [i]
else:
group[hashable_mat].append(i)
return group, bond_dict
def _extract_target_atoms(atoms_list,target):
"""targetに指定した原子のみを抽出する"""
atoms_list = ma.copy(atoms_list)
atoms_n = len(atoms_list[0])
del_idx = [i for i in range(atoms_n) if i not in target]
for atoms in atoms_list:
del atoms [del_idx]
return atoms_list