Way to the science

VASP中如何消虚频

VASP不同于Gaussian,可以结构优化与频率计算一起进行(opt+freq),VASP中结构优化与频率计算是分开的。需先优化好结构后再进行频率计算。同样的,过渡态也是如此,需要将过渡态结构的CONTCAR复制作为新的POSCAR进行计算,以验证结构的稳定性和过渡态结构的准确性。
然而,在VASP计算中,常常发现优化好的结构会存在虚频,尤其是体系中除固体外,还有气体分子。对此,总结了一些相关的处理办法,具体还请依据个人体系进行调整。

一、虚频较大,应该消除

(1)调整结构
若虚频波数达到几百,建议还是重新调整结构进行计算。除依据体系的特殊性和反应路径外,也可利用脚本消除不想要的虚频以实现消虚频工作。参考网上www.bigbrosci.com 给的脚本,根据虚频的dx、dy、dz微调结构重新计算。在此提供参考,后续可根据自己的需要以完善脚本。

import numpy as np
from ase.io import read, write
import os

# Define the paths
outcar_path = '/home/lmm/ZnGa2O4-highPrecision/species-ads/Ovac1-CO-CH3OH/CHO/freq111/OUTCAR'
poscar_path = '/home/lmm/ZnGa2O4-highPrecision/species-ads/Ovac1-CO-CH3OH/CHO/freq111/POSCAR'

*注:此处频率输出文件依据你的实际情况修改,本人此代码依托Jupyter运行,旨在直接读取OUTCAR中的虚频信息,调整后指定输出新POSCAR的路径。*
#获取虚频开始行
l_position = 0  #虚频振动方向部分在OUTCAR中的起始行数
with open(outcar_path) as f_in:
    lines = f_in.readlines()
    wave_num = 0.0
    for num, line in enumerate(lines):
        if 'f/i' in line:
            wave_tem = float(line.rstrip().split()[6])
            if wave_tem > wave_num: #获取最大的虚频
                wave_num = wave_tem
                l_position = num + 2

# ASE读POSCAR
model = read(poscar_path)
model_positions = model.get_positions()
num_atoms = len(model)

# 获取虚频对应的OUTCAR部分 
vib_lines = lines[l_position:l_position + num_atoms]
vib_dis = []
for line in vib_lines:
    vib_infor = [float(i) for i in line.rstrip().split()[3:]] # dx, dy, dz对应的那三列
    vib_dis.append(vib_infor)
vib_dis = np.array(vib_dis) #将振动部分写到一个array中。

# 微调结构
new_positions = model_positions + vib_dis * 0.4  # 0.4是微调的校正因子,即虚频对应振动位移的0.4,具体多大自己根据经验调。
model.positions = new_positions

#保存结构
write('POSCAR_new', model, vasp5=True)

注:以上脚本只针对所有虚频中的最大虚频进行调整
import numpy as np
from ase.io import read, write
import os

# Define the paths
outcar_path = '/home/lmm/ZnGa2O4-highPrecision/CI-NEB/Ovac1-CH2-CO--CH2CO/test2/Conti-/freq111/OUTCAR'
poscar_path = '/home/lmm/ZnGa2O4-highPrecision/CI-NEB/Ovac1-CH2-CO--CH2CO/test2/Conti-/freq111/POSCAR'

# Read POSCAR
model = read(poscar_path)
model_positions = model.get_positions()
num_atoms = len(model)

# Process OUTCAR to find all frequencies labeled with "f/i" and their displacement vectors
all_freqs = []
with open(outcar_path) as f:
    lines = f.readlines()
    for num, line in enumerate(lines):
        if 'f/i' in line:
            freq_value = float(line.rstrip().split()[6])  # Frequency value
            start_line = num + 2  # Start line for the displacement vectors
            vib_lines = lines[start_line:start_line + num_atoms]
            vib_dis = np.array([[float(i) for i in line.rstrip().split()[3:]] for line in vib_lines])
            all_freqs.append((freq_value, vib_dis))

# Print all frequencies with "f/i"
print("Frequencies with 'f/i' label:")
for idx, (freq, _) in enumerate(all_freqs):
    print(f"Frequency {idx+1}: {freq} cm^-1")
# Optionally, choose a specific frequency to adjust the structure (optional)
# Example: Adjust using the first frequency (index 0)
if all_freqs:
    selected_freq, vib_displacements = all_freqs[0]  # Assuming you still want to use the first one for some operation

*注:0表示根据所有虚频中的第一个虚频进行调整*

    # Adjust the structure based on the chosen frequency
    correction_factor = 0.4  # Adjustment scale factor
    new_positions = model_positions + vib_displacements * correction_factor
    model.positions = new_positions

    # Save the adjusted structure
    write('POSCAR_new', model, vasp5=True)
else:
    print("No frequencies with 'f/i' label found.")
import numpy as np
from ase.io import read, write
import os

# Define the paths
outcar_path = '/home/lmm/CO-O2-oxidation/ZnO2/freq111/OUTCAR'
poscar_path = '/home/lmm/CO-O2-oxidation/ZnO2/freq111/POSCAR'

# Read POSCAR
model = read(poscar_path)
model_positions = model.get_positions()
num_atoms = len(model)

# Process OUTCAR to find all frequencies labeled with "f/i" and their displacement vectors
all_freqs = []
with open(outcar_path) as f:
    lines = f.readlines()
    for num, line in enumerate(lines):
        if 'f/i' in line:
            freq_value = float(line.rstrip().split()[6])  # Frequency value
            start_line = num + 2  # Start line for the displacement vectors
            vib_lines = lines[start_line:start_line + num_atoms]
            vib_dis = np.array([[float(i) for i in line.rstrip().split()[3:]] for line in vib_lines])
            all_freqs.append((freq_value, vib_dis))

# Print all frequencies with "f/i"
print("Frequencies with 'f/i' label:")
for idx, (freq, _) in enumerate(all_freqs):
    print(f"Frequency {idx+1}: {freq} cm^-1")

# Check if there are at least two "f/i" frequencies
if len(all_freqs) >= 2:
    # Get displacement vectors for the first and second frequencies
    _, vib_displacements1 = all_freqs[0]
    _, vib_displacements2 = all_freqs[1]
    _, vib_displacements2 = all_freqs[2]
    _, vib_displacements2 = all_freqs[3]

    # Sum the displacement vectors
    combined_displacements = vib_displacements1 + vib_displacements2

*注:此代码可实现对指定的两个及两个以上的虚频结构同时调整并生成一个新的POSCAR*

    # Apply the combined displacements to the initial positions
    correction_factor = 0.4  # Adjustment scale factor
    new_positions = model_positions + combined_displacements * correction_factor
    model.positions = new_positions

    # Save the adjusted structure
    write('POSCAR_new', model, vasp5=True)
    print("Structure adjusted using the first two 'f/i' frequencies and saved to POSCAR_new.")
else:
    print("Not enough 'f/i' frequencies found to perform the adjustment.")

以上的代码均为参考,请根据自己的体系进行调整。

(2)调整INCAR参数
POTIM更改为0.015;
IBRION =6;
增加K点;
ISYM = -1;
提高EDIFF/EDIFFG精度重新优化

二、小虚频

通过第三方可视化软件(如Jmol)查看结构的频率信息时,发现小虚频往往对应气体分子的受损转动/平动但被近似为振动模式,一般有两种处理方式。
(1)文献中提到如果出现虚频(通常代表系统的不稳定性或过渡态),这些频率会被替换为6.8 meV。这样的处理是基于以下考虑:一个振动模式的熵会随着频率的趋近于零而发散。当振动频率非常低(接近或等于0)时,对应的熵在数学上会变得无限大,这在物理上是不合理的。因此,设定一个最低的能量阈值6.8 meV,这对应于600 K温度下约3kB的熵值,以确保低于这个频率的振动模式的熵被其他效应所限制。
这种处理方法有助于避免在计算熵或其他热力学属性时出现数值问题,同时保证了模型的物理合理性。这种做法是在模拟中常见的一种技术处理,用以处理和解释系统中可能出现的非物理的低频振动模式。(详见文献Thermochemistry and micro-kinetic analysis of methanol synthesis on ZnO (0001))
(2)尽可能消除,实在无法消除可不再进行处理,这可能源于有限差分法普遍存在的误差。www.bigbrosci.com 也曾提及吸附气体的最后5/6个为受损转动/平动。(若为线性分子为最后5个,非线性分子为最后6个)

Leave a comment

Your email address will not be published. Required fields are marked *