import ase
import ase.io
from ase.units import Bohr,Rydberg,kJ,kB,fs,Hartree,mol,kcal
from ase.constraints import FixAtoms
import matlanticgrrm.atomslist as ma
import copy
[ドキュメント]class GrrmData():
"""GRRMのデータをまとめたクラス
Parameters:
eqlog: str or path object
\*_EQ_list.logのパス.
tslog: str or path object
\*_TS_list.logのパス.
ptlog: str or path object
\*_PT_list.logのパス.
comfile: str or path object
comパス.FrozenAtomsがある場合に使用する.
poscar: str or path object
poscarのパス,周期境界条件にする際に使用.
constraint: bool
Trueの場合,FrozenAtomsが固定される.FixAtomsに設定される.
Attributes:
eq: class
:class:`EQ` クラス
ts: class
:class:`TS` クラス
pt: class
:class:`PT` クラス
"""
def __init__(self,eqlog=None,tslog=None,ptlog=None,comfile=None,poscar=None,constraints=False):
self.eq = EQ(eqlog,comfile,poscar,constraints)
self.ts = TS(tslog,comfile,poscar,constraints)
self.pt = PT(ptlog,comfile,poscar,constraints)
def set_eq(self):
pass
def set_ts(self):
pass
def set_pt(self):
pass
def __add__(self,other):
"""足し算できる(+)"""
obj = GrrmData()
if hasattr(self.eq,"atoms") and hasattr(other.eq,"atoms"):
obj.eq = self.eq + other.eq
if hasattr(self.ts,"atoms") and hasattr(other.ts,"atoms"):
obj.ts = self.ts + other.ts
if hasattr(self.pt,"atoms") and hasattr(other.pt,"atoms"):
obj.pt = self.pt + other.pt
return obj
def __iadd__(self,other):
"""足し算できる(+=)"""
return self + other
class _EqTsPt():
"""EQとTSPTクラスの親クラス."""
eqtspt = None #"EQ"or"TS"or"PT"が入る(子クラスでは)
def add(self,other,cls):
if self.eqtspt == other.eqtspt:
obj = cls()
obj.atoms = ma.copy(self.atoms) + ma.copy(other.atoms)
obj.energies = copy.copy(self.energies) + copy.copy(other.energies)
return obj
else:
raise MatranticGrrmError(f"{self.eqtspt}と{other.eqtspt}は結合できません")
[ドキュメント]class EQ(_EqTsPt):
"""EQの情報をまとめたクラス.
Parameters:
logfile: str or pathlib.Path
\*EQ_list.logファイルパス
comfile: str or pathlib.Path
\*comファイルパス
poscar: str or pathlib.Path
POSCARまたはCONTCARファイルパス. 周期境界条件(pbc=True)にする場合に設定
constraints: bool
FrozenAtomsをConstraintsに設定する場合True
"""
eqtspt = "EQ"
def __init__(self,logfile=None,comfile=None,poscar=None,constraints=False):
if logfile:
with open(logfile, "r") as f:
logtext = f.readlines()
if is_EQ_TS_PT(logtext[0]) == self.eqtspt:#EQ,PT,TSファイルである事を確認
(self.atoms,self.energies,
self.move_chemical_symbols,self.move_positions,
self.frozen_chemical_symbols,self.frozen_positions) = read_log(logtext,comfile,poscar,constraints)
else:
raise MatranticGrrmError(f'*{self.eqtspt}_list.logファイルを認識できません.')
def __add__(self,other):
"""足し算できる(+)"""
return super().add(other,EQ)
def __iadd__(self,other):
"""足し算できる(=+)"""
return self + other
class _TSPT(_EqTsPt):
"""TS,PTの親クラス."""
eqtspt = None
def __init__(self,logfile=None,comfile=None,poscar=None,constraints=False):
if logfile:
with open(logfile, "r") as f:
logtext = f.readlines()
if is_EQ_TS_PT(logtext[0]) == self.eqtspt:#EQ,PT,TSファイルである事を確認
(self.atoms,self.energies,
self.move_chemical_symbols,self.move_positions,
self.frozen_chemical_symbols,self.frozen_positions) = read_log(logtext,comfile,poscar,constraints)
self.connections = read_connection(logtext)
else:
raise MatranticGrrmError(f'*{self.eqtspt}_list.logファイルを認識できません.')
[ドキュメント]class TS(_TSPT):
"""TSの情報をまとめたクラス.
Parameters:
constraints: bool
FrozenAtomsをConstraintsに設定する場合True
その他のプロパティは :class:`EQ` クラス を参照
"""
eqtspt = "TS"
def __init__(self,logfile=None,comfile=None,poscar=None,constraints=False):
super().__init__(logfile,comfile,poscar,constraints)
def __add__(self,other):
return super().add(other,TS)
def __iadd__(self,other):
return self + other
[ドキュメント]class PT(_TSPT):
"""PTの情報をまとめたクラス.
Parameters:
constraints: bool
FrozenAtomsをConstraintsに設定する場合True
その他のプロパティは :class:`EQ` クラス を参照
"""
eqtspt = "PT"
def __init__(self,logfile=None,comfile=None,poscar=None,constraints=False):
super().__init__(logfile,comfile,poscar,constraints)
def __add__(self,other):
return super().add(other,PT)
def __iadd__(self,other):
return self + other
[ドキュメント]def log2atoms(logfile,comfile=None,poscar=None,constraints=False):
"""GRRMのlogファイルを読み込んでase.Atomsのリストを返す
Parameters:
logfile: str or path object
\*_list.logファイル
comfile: strr or path object
comファイル.FrozenAtomsがある場合に設定.
poscar: strr or path object
周期境界条件にする場合(pbc=True)の場合に設定.
constraint: bool
Trueの場合,FrozenAtomsがase.FixAtomsに設定される.
Return:
list: ase.Atomsを要素とするリスト.
"""
with open(logfile, "r") as f:
logtext = f.readlines()
if is_EQ_TS_PT(logtext[0]) in ["EQ","TS","PT"]:# logファイルの中身かあっているか確認
atoms, _,_,_,_,_ = read_log(logtext,comfile,poscar,constraints)
return atoms
def is_EQ_TS_PT(ini_text):
"""logファイルの1行目の文字からEQ,TS,PTのどのファイルかを返す"""
if ini_text == "List of Equilibrium Structures\n":
return "EQ"
elif ini_text == "List of Transition Structures\n":
return "TS"
elif ini_text == "List of Path Top (Approximate TS) Structures\n":
return "PT"
else:
raise MatranticGrrmError(f'*_list.logファイルを認識できません.')
def read_log(logtext,comfile=None,poscar=None,constraints=False):
"""\*_list.logファイルのテキストからase.Atomsのリストと(GRRMで計算された)エネルギーを返す"""
def int2float(positions):
return [[float(xyz) for xyz in atom.split()[1:]] for atom in positions]
ini_idx = [i for i,text in enumerate(logtext) if "#" in text]
fin_idx = [i for i,text in enumerate(logtext) if "Energy =" in text]
if len(ini_idx)*2 == len(fin_idx):
"""ReEnergyの場合,`Energy =`が2つ出るため"""
fin_idx = [i for i in fin_idx if i%2 == 0] #偶数個目のEnergy =の値だけ読み取る
move_symbols = [text.split()[0] for text in logtext[3:fin_idx[0]]]
natoms = len(move_symbols)
energies = [float(logtext[i].split()[2])*Hartree for i in fin_idx]
move_positions_list = [logtext[i+1:i+natoms+1] for i in ini_idx]
move_positions_list = [int2float(positions) for positions in move_positions_list] #floatに型変換
if comfile:
with open(comfile) as f:
comtext = f.readlines()
for i,text in enumerate(comtext):
if "Frozen Atoms" in text:
ini_idx = i
elif "Options" in text:
fin_idx = i
break
frozen_position = int2float(comtext[ini_idx+1:fin_idx])
frozen_symbol = [atom.split()[0] for atom in comtext[ini_idx+1:fin_idx]]
positions_list = [positions+frozen_position for positions in move_positions_list]
fix = FixAtoms(mask=[False]*len(move_symbols) + [True]*len(frozen_symbol)) if constraints else None
symbols = move_symbols + frozen_symbol
atoms = [ase.Atoms(symbols,positions,constraint=fix) for positions in positions_list]
else:
atoms = [ase.Atoms(move_symbols,positions) for positions in move_positions_list]
frozen_symbol = []
frozen_position = []
if poscar:
poscar_obj = ase.io.read(poscar)
cell = poscar_obj.get_cell()
for atoms_obj in atoms:
atoms_obj.set_cell(cell)
atoms_obj.set_pbc(True)
return atoms, energies,move_symbols,move_positions_list,frozen_symbol,frozen_position
def read_connection(log_text):
"""\*_list.logファイルからCONNECTIONを読み取る"""
connections = [(int(text.split()[2]) if text.split()[2].isdecimal() else text.split()[2],
int(text.split()[4]) if text.split()[4].isdecimal() else text.split()[4])
for text in log_text
if "CONNECTION" in text]
return connections
[ドキュメント]def get_connections(logfile):
"""GRRMのlogファイルを読み込んでCONNECTIONSを返す
Parameters:
logfile: str or path object
(TS or PT)_list.logファイル
Return: list of tuple of int
CONNECTIONSのリスト
"""
with open(logfile,"r") as f:
l = f.readlines()
return read_connection(l)
[ドキュメント]def atoms2log(name, ini_list,ts_list,fin_list,com=False,calc_func=None):
"""GRRMの\*list.logのようなファイルを作成する
ase.Atomsのリストを読みとり,\*EQ_list.logファイルと\*TS_list.logファイルを作成する
Parameters:
name: str
{name}_EQ_list.logの形でファイルを作成する.
ini_list: list of Atoms
初期構造のase.Atomsのリスト,要素の中にNoneが含まれている場合
ts_list: list of Atoms
TS構造のase.Atomsのリスト,TS計算がうまくいかなかったものは要素をNone
fin_list: list of Atoms
終構造のase.Atomsのリスト
com: bool
Trueにすると,comファイルを作成する. constraintsに設定している部分がFixAtomsになる
calc_func: functions
calculatorを設定していない場合はcalculatorを返す関数を設定する
"""
if calc_func:
ma.set_calculator(ini_list,calc_func)
ma.set_calculator(ts_list,calc_func)
ma.set_calculator(fin_list,calc_func)
eq_file_name = f"{name}_EQ_list.log"
ts_file_name = f"{name}_TS_list.log"
with open(eq_file_name,"w") as f:
f.write("List of Equilibrium Structures\n\n")
with open(ts_file_name,"w") as f:
f.write("List of Transition Structures\n\n")
def write_log(file,text_list):
with open(file,"a") as f:
f.writelines(text_list)
# FrozenAtomaの取得
for eq0 in ini_list:
if eq0:
frozen_idx = []
for const in eq0.constraints:
frozen_idx += const.get_indices().tolist()
frozen_idx = set(frozen_idx)
positions = eq0.get_positions()
symbols = eq0.get_chemical_symbols()
frozen_atoms = [f"{s}\t " + "\t ".join(p.astype(str).tolist()) + "\n"
for i,(s,p) in enumerate(zip(symbols,positions)) if i in frozen_idx]
move_atoms = [f"{s}\t " + "\t ".join(p.astype(str).tolist()) + "\n"
for i,(s,p) in enumerate(zip(symbols,positions)) if not i in frozen_idx]
break
def eq_text(ini, eq_n):
title = [f"# Geometry of EQ {eq_n}, SYMMETRY = -- \n"]
positions = ini.get_positions()
symbols = ini.get_chemical_symbols()
positions_txt = [f"{s}\t " + "\t ".join(p.astype(str).tolist()) + "\n"
for i,(s,p) in enumerate(zip(symbols,positions)) if not i in frozen_idx]
energy = ini.get_potential_energy()
energy_txt = [f"Energy = {energy} ({energy} : 0.000000000000)\n\n"]
text = title + positions_txt + energy_txt
return text
def ts_text(ts, ts_n, ini_idx, fin_idx):
title = [f"# Geometry of TS {ts_n}, SYMMETRY = -- \n"]
positions = ts.get_positions()
symbols = ts.get_chemical_symbols()
positions_txt = [f"{s}\t " + "\t ".join(p.astype(str).tolist()) + "\n"
for i,(s,p) in enumerate(zip(symbols,positions)) if not i in frozen_idx]
energy = ts.get_potential_energy()
energy_txt = [f"Energy = {energy} ({energy} : 0.000000000000)\n"]
connections_txt = [f"CONNECTION : {ini_idx} - {fin_idx}\n\n"]
text = title + positions_txt + energy_txt + connections_txt
return text
eq_n = 0 #EQ番号
ts_n = 0 #TS番号
for ini,ts,fin in zip(ini_list,ts_list,fin_list):
if ini:
write_log(eq_file_name, eq_text(ini,eq_n))
ini_idx = eq_n
eq_n += 1
if fin:
write_log(eq_file_name, eq_text(fin,eq_n))
fin_idx = eq_n
eq_n += 1
if ts:
ini_idx = ini_idx if ini else "??"
fin_idx = fin_idx if fin else "??"
write_log(ts_file_name, ts_text(ts, ts_n, ini_idx, fin_idx))
ts_n += 1
if com:
with open(f"{name}.com","w") as f:
f.write("%link=non-supported\n")
f.write("#SC-AFIR\n\n")
f.write("0 1\n")
f.writelines(move_atoms)
if len(frozen_atoms) > 0:
f.write("Frozen Atoms\n")
f.writelines(frozen_atoms)
f.write("Options\n")
class MatranticGrrmError(Exception):
"""例外クラス"""
pass