新建
ABACUS计算模拟实例 | II. C2H5OH的振动模式与频率计算
FermiNoodles
weiqingzhou@whu.edu.cn
推荐镜像 :abacus:3.6.0-user-guide
推荐机型 :c16_m32_cpu
赞 2
1
4
目录
数据集
ABACUS-cases(v18)
ABACUS计算模拟实例 | II. C2H5OH的振动模式与频率计算
代码
文本
©️ Copyright 2023 @ Authors
作者：Jinghang Wang (AISI电子结构团队实习生) 📨
日期：2024-6-21
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：由于截止目前（2024.6.21）ASE-ABACUS在对ABACUS:3.6.3版本的支持上仍存在问题，因此我们选择ABACUS:3.6.0版本完成本教程中的计算。点击上方的 开始连接 按钮，选择 abacus:3.6.0-user-guide 镜像及 c16_m32_cpu 节点配置，并在连接后点击右上方的切换内核，选择Bash内核稍等片刻即可运行。
代码
文本
🎯 本教程旨在帮助初学者学习使用ABACUS以及ASE-ABACUS接口完成振动模式与频率计算，并且可视化振动模式
代码
文本
本节教程目录如下：
1. 乙醇分子建模及弛豫计算
2. 振动频率的提取和模式分析
代码
文本
在使用ABACUS进行计算时，所有操作都需要在指定的文件夹目录下进行。我们已经准备好了数据集ABACUS-cases，，其中文件夹“2”对应本节教程。
代码
文本
首先我们将数据集中的文件复制到根目录下（否则可能碰到文件没有执行权限的情况）：
代码
文本
[1]
cp -r /bohr/ABACUS-cases-nu6t/v18/ABACUS-cases/2 .
代码
文本
[2]
cd 2 && ls
ase-abacus relax vib
代码
文本
可以看到有一共有两个文件夹，relax对应乙醇分子建模以及弛豫计算，vib文件对应振动频率的提取和模式分析。
代码
文本
1. 计算前的准备
代码
文本
本教程中的计算设计ASE软件以及ASE-ABACUS接口，因此我们先来做一个简单的介绍。
代码
文本
ASE (Atomic Simulation Environment)提供了一组用于设置、运行和分析原子模拟的Python工具。ABACUS作为计算器与ASE工具一起使用，可以通过ASE可视化STRU文件、生成与转换STRU文件、调用ABACUS进行自洽、结构弛豫、振动频率等计算及读取ABACUS计算结果等。它作为ASE的外部项目存在，并由ABACUS开发人员维护。官方网站见：ASE-ABACUS，使用教程见：ASE-ABACUS接口使用方法介绍。
代码
文本
下面我们就先安装最新版本的ASE-ABACUS。
代码
文本
[3]
#!/bin/bash
# 检查是否安装了 ase 软件包
if pip show ase > /dev/null 2>&1; then
echo "ase已安装，正在删除..."
pip uninstall -y ase
echo "ase已成功删除。"
else
echo "ase未安装。"
fi
ase已安装，正在删除... Found existing installation: ase 3.23.1b1 Uninstalling ase-3.23.1b1: Successfully uninstalled ase-3.23.1b1 WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv ase已成功删除。
代码
文本
[4]
#安装最新版的ASE-ABACUS
git clone https://gitlab.com/1041176461/ase-abacus.git
cd ase-abacus
pip install .
fatal: destination path 'ase-abacus' already exists and is not an empty directory. Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple Processing /2/ase-abacus Installing build dependencies ... done Getting requirements to build wheel ... done Preparing metadata (pyproject.toml) ... done Requirement already satisfied: numpy>=1.18.5 in /opt/mamba/lib/python3.10/site-packages (from ase==3.23.1b1) (1.26.4) Requirement already satisfied: scipy>=1.6.0 in /opt/mamba/lib/python3.10/site-packages (from ase==3.23.1b1) (1.13.0) Requirement already satisfied: matplotlib>=3.3.4 in /opt/mamba/lib/python3.10/site-packages/matplotlib-3.9.0rc2-py3.10-linux-x86_64.egg (from ase==3.23.1b1) (3.9.0rc2) Requirement already satisfied: contourpy>=1.0.1 in /opt/mamba/lib/python3.10/site-packages/contourpy-1.2.1-py3.10-linux-x86_64.egg (from matplotlib>=3.3.4->ase==3.23.1b1) (1.2.1) Requirement already satisfied: cycler>=0.10 in /opt/mamba/lib/python3.10/site-packages/cycler-0.12.1-py3.10.egg (from matplotlib>=3.3.4->ase==3.23.1b1) (0.12.1) Requirement already satisfied: fonttools>=4.22.0 in /opt/mamba/lib/python3.10/site-packages/fonttools-4.51.0-py3.10.egg (from matplotlib>=3.3.4->ase==3.23.1b1) (4.51.0) Requirement already satisfied: kiwisolver>=1.3.1 in /opt/mamba/lib/python3.10/site-packages/kiwisolver-1.4.5-py3.10-linux-x86_64.egg (from matplotlib>=3.3.4->ase==3.23.1b1) (1.4.5) Requirement already satisfied: packaging>=20.0 in /opt/mamba/lib/python3.10/site-packages (from matplotlib>=3.3.4->ase==3.23.1b1) (22.0) Requirement already satisfied: pillow>=8 in /opt/mamba/lib/python3.10/site-packages/pillow-10.3.0-py3.10-linux-x86_64.egg (from matplotlib>=3.3.4->ase==3.23.1b1) (10.3.0) Requirement already satisfied: pyparsing>=2.3.1 in /opt/mamba/lib/python3.10/site-packages/pyparsing-3.1.2-py3.10.egg (from matplotlib>=3.3.4->ase==3.23.1b1) (3.1.2) Requirement already satisfied: python-dateutil>=2.7 in /opt/mamba/lib/python3.10/site-packages (from matplotlib>=3.3.4->ase==3.23.1b1) (2.8.2) Requirement already satisfied: six>=1.5 in /opt/mamba/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib>=3.3.4->ase==3.23.1b1) (1.16.0) Building wheels for collected packages: ase Building wheel for ase (pyproject.toml) ... done Created wheel for ase: filename=ase-3.23.1b1-py3-none-any.whl size=2807761 sha256=5e48f67308eab3dbad8314934bfa9684a9d21c54fe070eac135ef87bc8d21db6 Stored in directory: /root/.cache/pip/wheels/46/43/f3/9eebbb6a57e84592c144a7082adc88550e27800901f070c1f5 Successfully built ase Installing collected packages: ase Attempting uninstall: ase Found existing installation: ase 3.23.1b1 Can't uninstall 'ase'. No files were found to uninstall. Successfully installed ase-3.23.1b1 WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv [notice] A new release of pip is available: 24.0 -> 24.1.1 [notice] To update, run: pip install --upgrade pip
代码
文本
2. 乙醇分子建模及弛豫计算
代码
文本
安装完成后就可以开始进行计算。我们首先进入relax文件夹。其中仅有relax.py一个可执行脚本，其余均为输出文件。
代码
文本
[5]
cd ../relax && ls
Ethanol.traj OUT ase_sort.dat relax.py time.json
代码
文本
脚本如下。该脚本首先得到乙醇分子的构型，之后使用该模型进行弛豫计算。
代码
文本
[6]
cat relax.py
import os
from pathlib import Path
import numpy as np
from ase.calculators.abacus import Abacus, AbacusProfile
from ase.collections import g2
from ase.io import Trajectory, read, write
from ase.optimize import QuasiNewton
from ase.visualize import view
# get Ethanol model
atoms = g2['CH3CH2OH']
atoms.center(vacuum=5)
# env setting
out_dir = "OUT"
properties = ["energy", "forces", "stress"]
os.environ['ABACUS_PP_PATH'] = '/abacus-develop/tests/PP_ORB'
os.environ['ABACUS_ORBITAL_PATH'] = '/abacus-develop/tests/PP_ORB'
pseudo_dir = os.environ['ABACUS_PP_PATH']
basis_dir = os.environ['ABACUS_ORBITAL_PATH']
pp = {
'H':'H_ONCV_PBE-1.0.upf',
'C':'C_ONCV_PBE-1.0.upf',
'O':'O_ONCV_PBE-1.0.upf',
}
basis = {
'H': 'H_gga_8au_100Ry_2s1p.orb',
'C': 'C_gga_8au_100Ry_2s2p1d.orb',
'O': 'O_gga_7au_100Ry_2s2p1d.orb',
}
kpts = [1, 1, 1]
# INPUT setting
parameters = {
'suffix': 'Ethanol',
'calculation': 'scf',
'basis_type': 'lcao',
'ks_solver': 'genelpa',
'nspin' : 1,
'xc': 'pbe',
'ecutwfc': 50,
'scf_thr': 1e-7,
'kpts': kpts,
'pp': pp,
'basis': basis,
'pseudo_dir': pseudo_dir,
'basis_dir': basis_dir,
'smearing_method': 'gauss',
'smearing_sigma': 0.0037,
'cal_force': 1,
'cal_stress': 1,
'out_stru': 1,
'out_chg': 0,
}
# attch calculator to the model
profile = AbacusProfile(command='mpirun -n 16 abacus')
atoms.calc = Abacus(profile=profile, directory=out_dir, **parameters)
#relax
dyn = QuasiNewton(atoms, trajectory='Ethanol.traj')
dyn.run(fmax=0.05)
print('Energy:', atoms.get_potential_energy())
代码
文本
Q：为什么我们做的是弛豫计算，但是ABACUS的calculation参数却是scf呢？
A：这是因为，虽然我们进行的是弛豫计算，但是ABACUS是作为calculator计算器被导入并附加在体系上的，作用为计算系统的单点能、力等参数；而执行弛豫（迭代改变原子位置）操作的则是optimizor优化器，在这里为QuasiNewton。前者对应“电子步”，后者对应“离子步”。
代码
文本
运行该脚本即可得到弛豫后的乙醇分子构型。
⚠️这一步使用16核预计耗时5分钟。
代码
文本
[7]
python3 relax.py
Step[ FC] Time Energy fmax BFGSLineSearch: 0[ 0] 16:16:35 -840.685582 0.3340 BFGSLineSearch: 1[ 2] 16:17:24 -840.693133 0.3241 BFGSLineSearch: 2[ 4] 16:18:13 -840.695908 0.1112 BFGSLineSearch: 3[ 6] 16:19:03 -840.696649 0.1408 BFGSLineSearch: 4[ 8] 16:19:51 -840.697181 0.0877 BFGSLineSearch: 5[ 10] 16:20:40 -840.697463 0.0730 BFGSLineSearch: 6[ 12] 16:21:28 -840.697743 0.0535 BFGSLineSearch: 7[ 14] 16:22:17 -840.697816 0.0299 Energy: -840.6978157
代码
文本
弛豫计算结束后会产生Ethanol.traj文件和OUT文件夹，前者为ase专门的轨迹文件，后者为包含弛豫计算过程的文件夹。
代码
文本
[8]
ls
Ethanol.traj OUT ase_sort.dat relax.py time.json
代码
文本
[9]
ls OUT
C_ONCV_PBE-1.0.upf KPT abacus.err C_gga_8au_100Ry_2s2p1d.orb OUT.Ethanol abacus.json H_ONCV_PBE-1.0.upf O_ONCV_PBE-1.0.upf abacus.out H_gga_8au_100Ry_2s1p.orb O_gga_7au_100Ry_2s2p1d.orb time.json INPUT STRU
代码
文本
在结构弛豫计算的输出文件中，STRU_READIN_ADJUST.cif和STRU_SIMPLE.cif文件均为弛豫前构型的.cif格式文件（前者考虑对称性），STRU_NOW.cif文件为弛豫后的cif文件；STRU_ION${i}_D和 STRU_ION_D分别为为第 i 步弛豫后的和最终的STRU文件。
代码
文本
[10]
ls OUT/OUT.Ethanol
INPUT STRU_ION4_D STRU_ION_D running_relax.log NOW_SPIN1_CHG.cube STRU_ION5_D STRU_NOW.cif running_scf.log OLD1_SPIN1_CHG.cube STRU_ION6_D STRU_READIN_ADJUST.cif warning.log STRU_ION1_D STRU_ION7_D STRU_SIMPLE.cif STRU_ION2_D STRU_ION8_D istate.info STRU_ION3_D STRU_ION9_D kpoints
代码
文本
[11]
#将弛豫得到的稳定构型复制到vib文件夹中用于后续振动频率计算
cd ../vib && cp ../relax/OUT/OUT.Ethanol/STRU_ION_D ./STRU
代码
文本
2. 振动频率的提取和模式分析
代码
文本
有了稳定的构型之后，我们就可以开始进行振动频率的计算。振动频率的计算在vib文件夹中进行，其中的输入文件只有STRU, combine_jmol.py和vib.py文件，执行vib.py文件就可得到其余所有输出文件。
代码
文本
[12]
ls
OUT ase_sort.dat combined.xyz mode_21 mode_23 mode_25 vib STRU combine_jmol.py combined_for_jmol.xyz mode_22 mode_24 mode_26 vib.py
代码
文本
vib.py文件的内容和作用如下：
1.进行振动频率计算：这部分的主要内容仍旧是使用ABACUS计算的输入，然后将其作为计算器附加在构型上，调用ase的vibrations方法进行计算。
2.振动频率的可视化：这部分的流程是：对于我们关注的振动模式，首先对每一个模式都创建一个文件夹，然后使用vib.write_mode()方法在其中读取该模式在上一步计算中对应的轨迹文件 (.traj)，通过n_images控制读取的帧数，然后将轨迹文件分解为n_images个.xyz文件，最后再将这些.xyz文件通过combine_jmol.py脚本合并为一个combined_for_jmol.xyz文件，这样就可以在jmol中播放。
注意：这里专门指出是for jmol，这是因为不同的可视化软件对xyz文件的格式要求不同，这部分具体体现在combine_jmol.py脚本中。
代码
文本
[13]
cat vib.py
import os
import subprocess
import shutil
from ase.optimize import QuasiNewton
from ase.vibrations import Vibrations
from ase.thermochemistry import HarmonicThermo
from ase.io import read, write
from ase.io import Trajectory
from ase.visualize import view
from ase.calculators.abacus import Abacus, AbacusProfile
# env setting
os.environ['ABACUS_PP_PATH'] = '/abacus-develop/tests/PP_ORB'
os.environ['ABACUS_ORBITAL_PATH'] = '/abacus-develop/tests/PP_ORB'
## 1. vibration
stru = read('STRU', format='abacus')
pseudo_dir = os.environ['ABACUS_PP_PATH']
basis_dir = os.environ['ABACUS_ORBITAL_PATH']
out_dir = "OUT"
properties = ["energy", "forces", "stress"]
pp = {
'H':'H_ONCV_PBE-1.0.upf',
'C':'C_ONCV_PBE-1.0.upf',
'O':'O_ONCV_PBE-1.0.upf',
}
basis = {
'H': 'H_gga_8au_100Ry_2s1p.orb',
'C': 'C_gga_8au_100Ry_2s2p1d.orb',
'O': 'O_gga_7au_100Ry_2s2p1d.orb',
}
kpts = [1, 1, 1]
# INPUT setting
parameters = {
'suffix': 'Ethanol',
'calculation': 'scf',
'basis_type': 'lcao',
'ks_solver': 'genelpa',
'nspin' : 1,
'xc': 'pbe',
'ecutwfc': 50,
'scf_thr': 1e-5,
'kpts': kpts,
'pp': pp,
'basis': basis,
'pseudo_dir': pseudo_dir,
'basis_dir': basis_dir,
'smearing_method': 'gauss',
'smearing_sigma': 0.0037,
'cal_force': 1,
'cal_stress': 1,
'out_stru': 1,
'out_chg': 0,
}
# attch calculator to the model
profile = AbacusProfile(command='mpirun -n 16 abacus')
stru.calc = Abacus(profile=profile, directory=out_dir, **parameters)
# vibration calculation
vib = Vibrations(stru, nfree=2)
vib.run()
vib.summary()
## 2. trajectory record and visualizaiton
# mkdir for each mode
for i in range(21, 27):
folder_name = f'mode_{i}'
os.makedirs(folder_name, exist_ok=True)
# write and read traj file of a mode, ":" for all configure in traj file
vib.write_mode(i, nimages=20)
traj_file_src = f'vib.{i}.traj'
traj_file_dest = f'mode_{i}/vib.{i}.traj'
shutil.move(traj_file_src, traj_file_dest)
traj = read(traj_file_dest, index=':')
# traj2xyz: split the traj file into several frames in each mode_i file
for j, atoms in enumerate(traj):
frame_file = f'mode_{i}/frame_{j:02d}.xyz'
write(frame_file, atoms)
# combine all frames into a .xyz file for jmol
for i in range(21, 27):
folder_name = f'mode_{i}'
subprocess.run(['python', 'combine_jmol.py', folder_name], check=True)
代码
文本
[14]
cat combine_jmol.py
import glob
import sys
# 检查是否提供了目标文件夹路径作为命令行参数
if len(sys.argv) != 2:
print("Usage: python combine_jmol.py <path_to_xyz_files>")
sys.exit(1)
target_folder = sys.argv[1]
# 更新文件模式以匹配指定文件夹中的XYZ文件
file_pattern = f'{target_folder}/frame_*.xyz'
file_list = sorted(glob.glob(file_pattern))
# 更新输出文件的路径以在指定文件夹内生成
combined_file_path = f'{target_folder}/combined_for_jmol.xyz'
with open(combined_file_path, 'w') as outfile:
for file_name in file_list:
with open(file_name, 'r') as infile:
lines = infile.readlines()
# Write the number of atoms, which is the first line of each XYZ file
outfile.write(lines[0])
# Write a comment line. Here we add a frame number for clarity
frame_number = file_name.split('_')[1].split('.')[0]
outfile.write(f"Frame {frame_number}\n")
# Now write the atom lines from the third line onwards
outfile.writelines(lines[2:])
# Add a newline to separate each frame if it's not the last line already
if not lines[-1].isspace():
outfile.write('\n')
代码
文本
[15]
python3 vib.py
--------------------- # meV cm^-1 --------------------- 0 136.0i 1097.2i 1 43.1i 347.7i 2 40.7i 328.2i 3 34.3i 276.3i 4 17.8i 143.8i 5 4.4i 35.5i 6 1.6 12.9 7 13.0 104.7 8 26.4 213.0 9 78.9 636.6 10 80.8 651.8 11 97.4 785.9 12 120.6 972.8 13 131.2 1058.3 14 132.3 1066.9 15 145.5 1173.8 16 156.4 1261.3 17 158.1 1274.9 18 173.9 1402.7 19 175.6 1416.0 20 175.9 1418.5 21 373.9 3016.0 22 397.0 3202.2 23 403.6 3255.4 24 406.5 3278.3 25 412.7 3329.0 26 578.5 4665.6 --------------------- Zero-point energy: 2.120 eV
代码
文本
计算完成后频率信息会直接打印在命令行界面中。每一行的第一列为振动模式的序号，第二列为能量，第三列为波数。
由于乙醇分子共有九个原子，每个原子在每一个方向上有一个自由度，并且没有对称性，因此一共有二十七个自由度。再减去三个刚性移动的自由度，最终就产生了二十四种振动模式。其中能量为虚数的代表虚频。沿着虚频所对应的振动模式方向，体系的能量会降低，通常表示当前的结构不是最稳定的构型。
代码
文本
输出文件如下：
代码
文本
[16]
ls
OUT ase_sort.dat combined.xyz mode_21 mode_23 mode_25 vib STRU combine_jmol.py combined_for_jmol.xyz mode_22 mode_24 mode_26 vib.py
代码
文本
每个模式对应的文件夹中的内容：1个轨迹文件，20帧对应的xyz文件以及1个合并脚本。
代码
文本
[17]
ls mode_26
combined_for_jmol.xyz frame_04.xyz frame_09.xyz frame_14.xyz frame_19.xyz frame_00.xyz frame_05.xyz frame_10.xyz frame_15.xyz vib.26.traj frame_01.xyz frame_06.xyz frame_11.xyz frame_16.xyz frame_02.xyz frame_07.xyz frame_12.xyz frame_17.xyz frame_03.xyz frame_08.xyz frame_13.xyz frame_18.xyz
代码
文本
运行完上述脚本后，我们就可以进入每一个模式对应的文件夹中使用jmol可视化振动模式。在数据集中，我们已经保存了乙醇分子的六种振动频率最高的模式，对应的图片如下：
代码
文本
(a)
(b)
(c)
(d)
(e)
(f)
图 1 六种频率最高的振动模式。
代码
文本
振动模式分析：
第一种模式的振动频率为3741.6 cm-1，是乙醇分子中振动频率最大的模式。在该模式中，振动主要在羟基中的O原子和H原子的连线方向上进行。
第二至第四种振动模式主要体现在乙醇分子的第二个C原子(C2)和其相连的H原子上,振动频率分别为3067.8 cm-1，3060.4 cm-1和2982.3 cm-1。在第二种振动模式中，振动主要在甲基中的碳氢键方向上进行。当其中的一个C-H键压缩时，另外两个C-H键拉长；在第三种振动模式中，C2上的一个C-H键长度不变，另外两个C-H键交替伸缩，同时C2原子在垂直于C-C键的方向上进行小幅度的旋转；在第四种振动模式中，三个C-H键共同进行伸缩和拉长，同时C2原子在C-C键方向上进行小幅度的振动。
第五种和第六种振动模式主要体现在乙醇分子的第一个C原子和其相连的H原子上，振动频率分别为2938.2 cm-1和2906.5 cm-1。在第五种振动模式中，C1上的两个C-H键交替伸缩，同时C1原子在垂直于C-C键的方向上进行小幅度的旋转；在第六种振动模式中，C1上的C-H键保持同频率的伸缩，同时C1原子在C-C键方向上进行小幅度的振动。
代码
文本
恭喜你完成了第二个ABACUS计算模拟实例🎉
代码
文本
如果你还想尝试更多的计算实例，可以点击下方了解ABACUS计算模拟合集链接： ABACUS计算模拟实例 | 概述
代码
文本
已赞2
本文被以下合集收录
ABACUS@计算材料学 | 计算模拟实例
weiqingzhou@whu.edu.cn
更新于 2024-08-09
13 篇8 人关注
ABACUS计算模拟实例
FermiNoodles
更新于 2024-07-15
13 篇1 人关注
推荐阅读
公开
ABACUS计算模拟实例 | VIII. 基于HSE06的态密度与能带计算FermiNoodles
发布于 2024-05-23
1 转存文件
公开
ABACUS计算模拟实例 | IX. 表面能的计算FermiNoodles
发布于 2024-05-21
2 转存文件
评论
yuxiangc22
FermiNoodles