import ase
from ase import Atoms
from ase.io import write
from ase.vibrations import Vibrations
from .log import reset_logger,add_filehandler
import matlanticgrrm.atomslist as ma
from plotly.subplots import make_subplots
import plotly.graph_objects as go
import pandas as pd
import numpy as np
from scipy.signal import find_peaks
from pathlib import Path
import sys,os,math,pickle
from tqdm.notebook import tqdm
from logging import getLogger,StreamHandler,FileHandler,Formatter,DEBUG,INFO,WARNING,ERROR,CRITICAL
[ドキュメント]class ListVibrations():
"""ase.vibrations.VibrationsをAtomsのリストでも使えるようにしたもの
Parameters:
atoms_list (list(Atoms)) : Atomsのリスト
errorlog (str) : logのファイルパス
indexex (list) : 計算対象をリストで指定 <--ここまでがListVibrations特有の引数
indices (list) : 動かす原子をリストで指定 <-- ここからがVibrationsの引数aseのDocumentを参照
"""
def __init__(self,atoms_list,errorlog,indexes=None,indices=None, name='vib', delta=0.01, nfree=2):
self.atoms_list = atoms_list
self.errorlog = errorlog
self.indexes = indexes if indexes else [i for i in range(len(self.atoms_list))]
self.indices = indices
self.name = Path(name)
self.delta = delta
self.nfree = nfree
self.vib = []
if not Path(name).exists():
os.makedirs(name)
def __iter__(self):
"""Optimizerを使用するためには__iter__を実装する必要がある"""
return iter(self.vib)
def __len__(self):
"""Optimizerを使用するためには__len__を実装する必要がある"""
return len(self.vib)
def __getitem__(self, item):
return self.vib[item]
def mk_vib_func(self,atoms,foldername):
return Vibrations(atoms, indices=self.indices, name=f"{self.name}/{foldername}", delta=self.delta, nfree=self.nfree)
def _start_number(self):
self.logger = reset_logger() #loggerをリセット
if Path("stop.log").exists():
with open("stop.log","r") as f:
n = int(f.read()) #n番目から計算を始める
os.remove("stop.log")
add_filehandler(self.logger,"a",self.errorlog,DEBUG)
else:
n = 0
add_filehandler(self.logger,"w",self.errorlog,DEBUG)
return n
[ドキュメント] def run(self):
"""
| logの内容[logレベル]
| {1} : index未指定 [INFO]
| {2} : Atomsのオブジェクトが存在しない [INFO]
| {5} : エラー [CRITICAL]
Note:
Imaginary_Frequency.htmlとImaginary_Freq_Idx.txtを出力する.
実行中のnotebookと同じディレクトリにstop.logを作成すると計算を中断する.
"""
n = self._start_number() #n番目から計算を始める
self.vib = [None for _ in range(n-1)] if len(self.vib)==n-1 else self.vib
for i, atoms in enumerate(tqdm(self.atoms_list[n:]),start=n):
if Path("stop.log").exists():
if write_steplog is False:
with open("stop.log","w") as f:
f.write(str(i))
write_steplog = True
continue
write_steplog = False
if not i in self.indexes:
self.logger.info(f"i={i} index未指定 return 1")
self.vib.append(None)
elif not atoms:
self.logger.info(f"i={i} objectが存在しない return 2")
self.vib.append(None)
else:
try:
vib = self.mk_vib_func(atoms,str(i))
vib.run()
#runまで到達したら
self.vib.append(vib)
self.logger.info(f"i={i} 完了 return 0")
except Exception as e:
exc_type, exc_obj, exc_tb = sys.exc_info()
fname = os.path.split(exc_tb.tb_frame.f_code.co_filename)[1]
self.logger.critical(f"i={i} {exc_type} {fname} {exc_tb.tb_lineno} return 5")
self.vib.append(None)
def clean(self,empty_files=False, combined=True):
"""self.nameのフォルダ内の全てjsonファイルを削除する"""
for vib_obj in self.vib:
if vib_obj:
vib_obj.clean(empty_files=empty_files,combined=combined)
[ドキュメント] def summary(self,html="vib_table.html",pkl=None):
"""
Parameters:
html: str or path object
htmlパス.htmlに表を保存する. Noneの場合は標準出力だが,数の多い場合は推奨しない
pkl:
pickleパス.指定した場合,各構造の虚数振動のモード番号をリストとしてpickleに保存する.
"""
sample_n = len(self.vib)
cols = 3
rows = math.ceil(sample_n/cols)
specs = [[{"type": "table"} for _ in range(cols)] for _ in range(rows)]
subplot_titles = [f"TS{i}" for i in range(sample_n)]
fig = make_subplots(
rows=rows, cols=cols,
specs=specs,
subplot_titles = subplot_titles
)
div_sample = []
for i in range(0, sample_n, cols):
"""データをcol分割する"""
div_sample.append(self.vib[i: i+cols])
def add_table(fig,df,row,col,fill_color):
fig.add_trace(
go.Table(
header=dict(values=list(df.columns),
fill = {"color":"rgb(68,114,196)"},
font = {"color":"rgb(255,255,255)","family" : "Arial"},
align='center'),
cells=dict(values=[df["<b>#</b>"], df["<b>meV</b>"], df["<b>cm^-1</b>"]],
fill = {"color":"rgb(233,235,245)"},
fill_color = fill_color,
font = {"family":"Arial"},
align='center')
),
row=row, col=col)
self.i_idx = [] #虚数振動のidx
for i, row_vib in enumerate(div_sample,start=1):
for j , vib_obj in enumerate(row_vib,start=1):
if vib_obj:
vib_data = vib_obj.get_vibrations()
text = vib_data.tabulate().split("\n") #summaryを取得
df = pd.DataFrame(
data=np.array([i.split() for i in text[3:-3]]),
columns=['<b>#</b>', '<b>meV</b>', '<b>cm^-1</b>'])
imaginary_vib_list = df.index[df['<b>cm^-1</b>'].str.contains('i')].tolist() #虚振動を含む行番号を抽出
imaginary_vib_list = [i for i in imaginary_vib_list if not i in [1,2,3,4]] #1~4は無視
if len(imaginary_vib_list) == 1: #虚振動が1つだげ存在するとき
self.i_idx.append(imaginary_vib_list[0])
fill_color = [["rgb(233,235,245)" if i!=imaginary_vib_list[0] else "rgb(255,155,155)" for i in range(len(df))]]*3
else:
fill_color = [["rgb(233,235,245)" for _ in range(len(df))]]*3
self.i_idx.append(None)
add_table(fig,df,i,j,fill_color)
else:
self.i_idx.append(None)
fig.update_layout(height=700*rows, width=1000,
title_text="振動")
if html:
fig.write_html(html)
else:
fig.show()
if pkl:
with open(pkl,"wb") as f:
pickle.dump(self.i_idx,f)
def combine(self):
pass
def fold(self):
pass
def get_energies(self):
pass
def get_frequencies(self):
pass
def get_mode(self):
pass
def get_vibrations(self):
pass
def iterdisplace(self):
pass
def iterimages(self):
pass
def split(self):
pass
def write_dos(self):
pass
def write_jmol(self):
pass
[ドキュメント] def write_mode(self,n_list=None,savename="vib_all", kT=0.02585199101165164, nimages=30):
"""trajファイルを作成する
Parameters:
n_list: list of integer
| 各構造のtrajとして書き出したい振動モードをリストで与える.
| Noneの場合, 虚振動を自動判断する.
savename: str
保存ファイル名（拡張子なし）
Note:
n_list=Noneを使う場合,事前にsummary()を実行しておく必要がある.
その他の引数はASEを参照
"""
if n_list is None:
n_list = self.i_idx
status_list = []
with ase.io.Trajectory(f"{savename}.traj", 'w') as traj:
for i,(n,vib_obj) in enumerate(zip(n_list,self.vib)):
if type(n) == int: # Noneの要素は受け付けない
status_list.append(0)
n %= len(vib_obj.get_energies())
for image in (vib_obj.get_vibrations().iter_animated_mode(n,temperature=kT, frames=nimages)):
traj.write(image)
else: # 要素がNoneの場合は空のAtomsを挿入する
status_list.append(2)
for _ in range(nimages):
traj.write(Atoms())
with open(f"{savename}.pickle","wb") as f:
pickle.dump((status_list,nimages),f)
[ドキュメント]def vib_graphs(imagess,html="vib_graph.html",is_traj="Reverse.traj",fs_traj="Forward.traj",ts_traj="TS.traj"):
"""
| Reverse.traj,Forward.traj,TS.trajが作成される.
| Reverse.traj,Forward.trajは構造最適化(IRC)を行なうことで始構造, 終構造となる.
Parameters:
imagess: 2D list of Atoms
Atomsの2次元リスト.ListVibration.write_mode()で作成したtrajファイルから用意する.
"""
energiess = [ma.get_potential_energy(images) if images else None for images in imagess]
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 = data
add_fig(fig,x,y,i,j,'black')
div_data = list(np.array_split(data, 3))[1] #3分割した内，真ん中のデータを抽出
peaks_i, _ = find_peaks(div_data)
if len(peaks_i) == 1: #ピークが1つだけあれば
x = peaks_i[0]+len(div_data)
y = data[x]
add_fig(fig,[x],[y],i,j,'red')
idx.append(x)
else:
idx.append(None)
else:
idx.append(None)
fig.update_layout(height=300*rows, width=1000,
title_text="エネルギーダイアグラム")
if html:
fig.write_html(html)
else:
fig.show()
is_list = []
ts_list = []
fs_list = []
for i, images in zip(idx,imagess):
if i:
is_list.append(images[i-1])
ts_list.append(images[i])
fs_list.append(images[i+1])
else:
is_list.append(Atoms())
ts_list.append(Atoms())
fs_list.append(Atoms())
write(is_traj,is_list)
write(ts_traj,ts_list)
write(fs_traj,fs_list)
return idx