Phoenix in fire

May 06, 2025 in 沉思录

从上周三开完组会一直到现在都还没怎么缓过神来，不能再这样下去了😟，得主动出击。早上起来忧心忡忡的，高分子物理没心情上，吃完早饭直接去鼎新馆吧。昨晚大概M$入门了，建模后提交了任务，但是可能INCAR里面参数设置得不够合理，今天下午提交的任务才收敛。尝试想了一些方案但是被师兄否定了，提高电解质中离子电导率的机理也许是清楚的，但是没有一篇明晰的文献手把手教我应该挑选那些元素掺杂来提高离子电导率，也许永远不会吧🥺。Keep on going，再试试吧，多读读文献，也许我现在的困境和英语只有60分的同学所面临的困境相似吧？——都是词汇量不够？

Anyway, Just keep on going! 如果走不出迷堆明天就去请教一下魏老师。

合抱之木，生于毫末；九层之台，起于累土；千里之行，始于足下

Comment and share

阳光总在风雨后

May 03, 2025 in 一些杂谈, 歌词赏析

歌词

[00:00.60]许美静 - 阳光总在风雨后
[00:01.60]词：陈佳明
[00:02.60]曲：陈佳明
[00:24.34]人生路上甜苦和喜忧
[00:29.24]愿与你分担所有
[00:34.15]难免曾经跌倒和等候
[00:39.30]要勇敢的抬头
[00:43.95]谁愿藏躲在避风的港口
[00:48.85]宁有波涛汹涌的自由
[00:53.91]愿是你心中灯塔的守候
[00:59.22]在迷雾中让你看透
[01:06.93]阳光总在风雨后
[01:11.89]乌云上有晴空
[01:16.84]珍惜所有的感动
[01:21.89]每一份希望在你手中
[01:26.89]阳光总在风雨后
[01:32.00]请相信有彩虹
[01:37.05]风风雨雨都接受
[01:41.81]我一直会在你的左右
[02:11.74]人生路上甜苦和喜忧
[02:16.98]愿与你分担所有
[02:21.35]难免曾经跌倒和等候
[02:26.71]要勇敢地抬头
[02:31.31]谁愿藏躲在避风的港口
[02:36.26]宁有波涛汹涌的自由
[02:41.36]愿是你心中灯塔的守候
[02:46.67]在迷雾中让你看透
[02:54.22]阳光总在风雨后
[02:59.62]乌云上有晴空
[03:04.38]珍惜所有的感动
[03:09.28]每一份希望在你手中
[03:14.43]阳光总在风雨后
[03:19.39]请相信有彩虹
[03:24.54]风风雨雨都接受
[03:29.24]我一直会在你的左右
[03:37.35]阳光总在风雨后
[03:41.95]乌云上有晴空
[03:46.76]珍惜所有的感动
[03:51.71]每一份希望在你手中
[03:56.86]阳光总在风雨后
[04:01.76]请相信有彩虹
[04:06.82]风风雨雨都接受
[04:11.65]我一直会在你的左右
[04:19.63]风风雨雨都接受
[04:24.13]我一直会在你的左右

歌词深度分析

1. 主题与中心思想

这首歌词的核心主题是逆境中的希望与陪伴，强调人生难免经历风雨（挫折与等待），但要坚持信念（“勇敢抬头”），因为光明终将到来（“阳光总在风雨后”）。中心情感是温暖而坚韧的鼓励，既有对个体勇气的呼唤，也有对陪伴力量的歌颂（“我一直会在你的左右”）。内在冲突体现在“避风港的安全”与“波涛汹涌的自由”之间的选择，最终指向成长必须直面挑战的价值观。

2. 情绪氛围与感受

歌词传递的情绪由低沉渐趋昂扬： - 开篇承认人生的“甜苦喜忧”与“跌倒”，带有淡淡的忧郁； - 副歌通过“阳光”“彩虹”等意象转向明亮，情绪升华成坚定与治愈； - 反复的“风风雨雨都接受”强化了坦然接纳苦难的豁达感。整体是哀而不伤、充满力量的基调。

3. 语言风格与用词特色

口语化与诗意结合：如“跌倒和等候”“避风的港口”通俗易懂，但“灯塔的守候”“乌云上有晴空”又充满诗意。
短句为主：增强节奏感和号召力（如“要勇敢的抬头”）。
关键词重复：“阳光”“风雨”“彩虹”形成记忆点，强化主题。

4. 意象与象征意义

风雨/乌云：象征人生困境；阳光/彩虹：象征希望与转机。
灯塔：象征指引与守护，暗喻陪伴者的角色；避风港口：象征安逸但可能停滞的状态。
文化延伸：呼应中国传统文化中“苦尽甘来”的哲学（如“守得云开见月明”），也符合普世励志价值观。

5. 修辞手法分析

比喻：将人生困境比作“风雨”，希望比作“阳光”“彩虹”。
对比：“避风港口”与“波涛汹涌的自由”形成安全与冒险的张力。
排比：副歌重复“阳光总在风雨后…请相信有彩虹”，增强感染力。
拟人化：乌云“上有晴空”，赋予自然现象情感暗示。

6. 节奏与结构

分段清晰：主歌（叙事铺垫）→ 副歌（情感升华）循环三次，结尾重复强调陪伴。
押韵工整：主歌“忧/有”“候/透”，副歌“后/空/动/中”“虹/受/右”，增强音乐性。
重复设计：副歌重复四次，形成洗脑记忆点，符合流行歌曲传播规律。

7. 文化/社会背景

时代共鸣：1990年代末华语乐坛盛行治愈系歌曲，反映经济快速发展下人们对压力释放的需求。
集体心理：歌词中“分担所有”“在你左右”暗合亚洲文化重视人际支持的特点，区别于西方个人英雄主义。
创作者视角：陈佳明擅长细腻情感表达（如为许美静写的《城里的月光》），本词延续了他“温柔鼓励”的风格。

8. 多重解读角度

个人成长视角：鼓励个体直面挫折，强调自我救赎（“勇敢抬头”）。
人际关系视角：突出陪伴的力量（“与你分担”“在你左右”），适合友情、亲情或爱情解读。
社会隐喻视角：风雨可象征时代变迁（如金融危机），阳光则寄托对未来的集体信心。

总结：这首歌词通过简洁有力的语言、鲜明的自然意象和循环递进的结构，成功构建了一个从困境到希望的叙事闭环，兼具文学美感和大众传播性，成为华语励志经典。

Comment and share

快乐无罪

May 03, 2025 in 一些杂谈, 歌词赏析

歌词

[00:00.43]快乐无罪 - 许美静
[00:01.62]词：许常德/陈佳明
[00:03.13]曲：陈建宁/陈政卿
[00:19.82]你要我快乐
[00:23.76]我非常快乐
[00:27.22]在不想太多的时候
[00:34.32]吻你是烟火 吻后是寂寞
[00:41.32]爱是我犯的差错
[00:49.38]我要你什么 我贪你什么
[00:56.73]在你想消失的时候
[01:04.12]想你是黑洞 它没有尽头
[01:10.72]就算流星也该坠落
[01:17.55]快乐有泪 怎会有泪
[01:21.23]我哭到心碎
[01:25.52]不解幸福会被永恒连累
[01:29.16]不解刹那是你要的美
[01:32.60]快乐有罪 怎会有罪
[01:36.29]我怎么拒绝
[01:40.55]也许到最后你我都反悔
[01:44.21]留下空荡的城市黑夜
[01:47.60]让我撒野
[02:22.36]我要你什么 我贪你什么
[02:29.56]在你想消失的时候
[02:36.93]想你是黑洞 它没有尽头
[02:43.74]就算流星也该坠落
[02:50.47]快乐有泪 怎会有泪
[02:54.10]我哭到心碎
[02:58.48]不解幸福会被永恒连累
[03:02.05]不解刹那是你要的美
[03:05.46]快乐有罪 怎会有罪
[03:09.17]我怎么拒绝
[03:13.44]也许到最后你我都反悔
[03:16.94]留下空荡的城市黑夜 让我撒野
[03:28.09]快乐有泪 怎会有泪
[03:31.54]我哭到心碎
[03:35.86]不解幸福会被永恒连累
[03:39.50]不解刹那是你要的美
[03:42.99]快乐有罪 怎会有罪
[03:46.58]我怎么拒绝
[03:50.86]也许到最后你我都反悔
[03:54.47]留下空荡的城市黑夜
[03:57.90]让我撒野

歌词深度分析

1. 主题与中心思想

这首歌词探讨了爱情中的矛盾与自我欺骗，核心情感是「伪装的快乐」与「真实的痛苦」之间的撕裂感。 - 表面：反复强调“快乐”（“我非常快乐”“快乐有泪/有罪”），实则揭露快乐背后的空虚。 - 内在冲突：渴望爱的永恒与接受爱的短暂（“不解幸福会被永恒连累”），以及自我麻痹（“在不想太多的时候”）与清醒痛苦（“我哭到心碎”）的对抗。

2. 情绪氛围与感受

主调：压抑的忧伤，夹杂自嘲与无力感。
矛盾情绪：
- “吻你是烟火，吻后是寂寞” —— 瞬间炽热与长久孤独的对比。
- “快乐有泪/有罪” —— 快乐被赋予负面色彩，暗示情感中的自我惩罚倾向。

3. 语言风格与用词特色

诗意化隐喻：用“烟火”“黑洞”“流星”等意象抽象化情感体验。
口语化反问：“怎会有泪”“怎会有罪”增强质问感，像自我对话。
重复与矛盾修辞：通过“快乐”与“泪/罪”的并置，制造语义冲突。

4. 意象与象征意义

烟火：象征爱情的短暂绚烂与毁灭性。
黑洞：比喻思念的吞噬性与无解（“没有尽头”）。
流星：暗指注定陨落的爱情，呼应“刹那的美”。
空荡的城市黑夜：投射内心的荒芜感，同时“撒野”暗示崩溃后的释放。

5. 修辞手法分析

矛盾修辞法：如“快乐有泪”“快乐有罪”，颠覆常规认知。
排比与反复：副歌段落的重复强化无力感，像陷入循环的执念。
拟人化：“幸福被永恒连累”将抽象概念人格化，暗示责任转移。

6. 节奏与结构

分段：主歌叙事（具体情境），副歌抒情（情绪爆发）。
押韵：交替押“ei”“ui”韵（如“泪/罪”“碎/累”），音调压抑。
节奏变化：副歌的短句“怎会有泪”“怎会有罪”加快节奏，模拟急促的控诉。

7. 文化/社会背景

90年代华语流行乐：反映都市情感疏离（“空荡的城市”），契合许美静“都市冷调”风格。
女性视角：歌词中被动性（“你怎么拒绝”）与主动宣泄（“让我撒野”）的拉扯，隐含传统性别角色与自我意识的冲突。

8. 多重解读角度

爱情层面：一段不对等的关系，一方逃避（“你想消失的时候”），另一方自我欺骗。
存在主义层面：对“快乐”本质的质疑——是否只是逃避现实的工具？
社会隐喻：“城市黑夜”可延伸为现代人的集体孤独，爱情成为短暂救赎却加剧虚无。

总结：这首歌词以“快乐”为伪装，揭露爱情中的痛苦本质。通过意象碰撞与矛盾修辞，完成了一场从自我欺骗到崩溃宣泄的情感解剖，最终指向现代人面对爱与孤独的永恒困境。

Comment and share

离子电导率的计算

May 01, 2025 in 从零开始的科研之旅, VASP

使用VASP计算离子电导率，似乎比较麻烦

一个能量较低、结构合理的起始三维构型是进行MD模拟的第一步

分子动力学模拟的基本步骤

确定起始构型
选用适当的力场和模拟软件
构建体系和能量最小化
平衡过程
数据采集过程
结果分析

根据系统的传播方式，分子模拟方法可以分为分子动力学MD和蒙特卡洛MC两类方法

分子动力学通过数值积分求解牛顿运动方程，生成系统轨迹，可研究结构、动力学和热力学性质。
蒙特卡洛方法基于概率规则生成新构型，侧重采样平衡态结构和热力学性质，不涉及时间动力学。

力场是一组描述分子中原子间相互作用的势能函数和参数的集合，用于计算系统的总能量和原子受力，从而驱动分子动力学（MD）模拟或蒙特卡罗（MC）采样：

将量子力学描述的分子相互作用简化为经典力学模型，降低计算成本。
通过参数化的势能函数，预测分子的结构、动力学和热力学性质（如键长、自由能、扩散系数等）。

固定电荷力场：固定原子电荷，不考虑环境引起的电子云变形（极化），计算效率高

生物分子力场
- AMBER
  - 广泛用于蛋白质、核酸模拟，参数优化针对溶液环境，如 ff14SB、ff19SB
- CHARMM
  - 广泛用于蛋白质、核酸模拟，参数优化针对溶液环境，如 ff14SB、ff19SB
- GROMOS
  - 用于膜蛋白和聚合物，参数简洁，计算速度快
小分子力场
- GAFF
  - 可参数化多种有机小分子，支持自定义分子
- UFF
  - 覆盖周期表多数元素，适合无机材料和有机分子

极化力场：显示考虑原子极化率（电子云随电场变化），通过附加自由度（如 Drude 振子、诱导偶极）模拟极化效应，精度更高但计算成本显著增加

AMOEBA：基于原子极化率和多极矩，适用于极性体系（如水、蛋白质 - 配体相互作用）。
Drude 模型：每个原子附加虚拟振子模拟极化，如 CHARMM-Drude 力场。
Polarizable AMBER (PAMBS)：在 AMBER 框架中引入极化参数。

粗粒化力场（Coarse-Grained Force Fields）：将多个原子简化为单个 “珠子”（如蛋白质的每个残基、脂质的疏水尾），降低计算复杂度，适用于大尺度系统（如细胞膜、聚合物）

Martini：广泛用于膜蛋白、病毒颗粒和软物质，4:1 粗粒化（4 个原子→1 个珠子）。
MARTINI 3.0：支持多尺度模拟，兼容全原子模型。

使用pymatgen提取数据

'''
分析AIMD结果，计算MSD 和 conductivity
'''
import os
from pymatgen.core.trajectory import Trajectory
from pymatgen.io.vasp.outputs import Xdatcar
from pymatgen.core import Structure
from pymatgen.analysis.diffusion.analyzer import DiffusionAnalyzer
import numpy as np
import pickle
import matplotlib.pyplot as plt

# 这一步是读取 XDATCAR，得到一系列结构信息
traj = Trajectory.from_file('XDATCAR')

# 这一步是实例化 DiffusionAnalyzer 的类
# 并用 from_structures 方法初始化这个类； 900 是温度，2 是POTIM 的值，1是间隔步数
# 间隔步数（step_skip）不太容易理解，但是根据官方教程:
# dt = timesteps * self.time_step * self.step_skip

diff = DiffusionAnalyzer.from_structures(traj,'Li',900,2,1)

# 可以用内置的 plot_msd 方法画出 MSD 图像
# 有些终端不能显示图像，这时候可以调用 export_msdt() 方法，得到数据后再自己作图
ax = diff.export_msdt("msd.dat")
#plt.show()

# 接下来直接得到 离子迁移率， 单位是 mS/cm
C = diff.conductivity

with open('result.dat','w') as f:
    f.write('# AIMD result for Li-ion\n')
    f.write('temp\tconductivity\n')
    f.write('%d\t%.2f\n' %(900,C))

画图

import numpy as np
import matplotlib.pyplot as plt

# 读取数据
data = np.loadtxt('msd.txt', comments='#')  # 假设你的文件名为 msd.txt

# 提取时间步长和 MSD 数据
time = data[:, 0]  # 第一列是时间步长
msd = data[:, 1]   # 第二列是总的 MSD

# 绘制 MSD 曲线
plt.figure(figsize=(10, 6))
plt.plot(time, msd, label='MSD', color='blue', linewidth=2)
plt.xlabel('Time (fs)', fontsize=14)
plt.ylabel('Mean Square Displacement (Å²)', fontsize=14)
plt.title('MSD vs Time', fontsize=16)
plt.legend(fontsize=12)
plt.grid(True)
plt.show()

System preparation
Minimization/Relaxation
Equilibration
Production

系统准备可以说是模拟中最关键的阶段，但在许多情况下却得到了最少的关注；如果你的系统准备有缺陷，这些缺陷可能会证明是致命的。最坏的可能结果是，准备好的系统并非你想要的（例如，它包含不正确的分子或质子化状态），但在化学上却是有效的，并且可以被你的力场很好地描述，因此没有错误地通过了剩余的步骤——实际上，这是系统准备问题中常见的后果。如果系统在随后的平衡步骤中行为良好，不应假设系统已正确准备；这里应该格外小心。

生产流程如下：

为什么需要有恒温器？

分子动力学模拟用于观察和获取某个研究系统的有趣性质。在很多情况下，为了模拟实验室条件下进行的实验，希望从正则（恒温）配分函数中进行抽样。一般来说，如果系统在模拟过程中必须保持温度不变，将采用某种温控算法。

分子动力学模拟的温度通常使用根据等分配定理定义的动能来测量： $$ \frac{3}{2}Nk_BT = <\sum_{i=1}^{N}\frac{1}{2}m_iv_i^2> $$ 带角标的括号表示温度是一个时间平均值。如果我们使用等分配定理来计算分子动力学模拟单时间点的温度快照，而不是时间平均，这个量称为瞬时温度。瞬时温度不一定等于目标温度；实际上，在正则配分函数中，瞬时温度应该围绕目标温度波动。

温控算法通过改变本质上微观可逆的牛顿运动方程来工作。因此，如果希望计算扩散系数等动力学性质，则最好不使用温控器；相反，在系统达到所需温度后，应关闭温控器。然而，虽然所有温控器都会给出非物理动力学，但有些研究发现它们对特定动力学性质的计算影响很小，并且它们在生产模拟中也常被使用。

一些恒温器

Gaussian
- 使用了高斯最小约束原理来确定维持瞬时温度所需的最小扰动力。这种温控器不采样正则分布；相反，它采样的是等速（恒动能）配分函数。但是，等速配分函数采样与正则配分函数相同的配置相空间，因此，使用任一配分函数都可以等效地获得位置依赖（结构）平衡性质。然而，速度依赖（动力学）性质在两个配分函数之间不会等效。这种温控器仅适用于某些高级应用。
Simple Velocity Rescaling
- 实施起来最容易的温控器之一,也是最不物理的温控器之一。这种温控器依赖于重新标粒子的动量，使模拟的瞬时温度恰好与目标温度相匹配，不建议使用
Berendsen
- 既不采样正则分布也不采样等速分布，它的使用可能导致模拟伪象，不建议使用
Bussi-Donadio-Parrinello (Canonical Sampling Velocity Rescaling)
- 正确采样了正则配分函数，也被称为“V-rescale”热浴方法，非常常用的热浴方法
Andersen
- 只能用于采样结构性质，因为动力学性质可能会受到突然碰撞的很大影响，不建议使用
Langevin
- 不了解
Noés-Hoover
- 不了解

压浴类型	原理	特点	适用性	局限性
Simple volume rescaling	每次修改系统体积使瞬时压强等于目标压强	操作简单直接	仅可用于模拟初始阶段，帮助系统初步接近目标压强	不能采样正确配分函数，不能用于生产采样；不平滑地接近目标压强，可能导致系统出现不真实的物理问题
Berendsen	将系统耦合到弱相互作用的压力浴，通过缩放因子定期缩放体积	比简单体积缩放更真实地接近目标压强，产生更真实的压力波动	对于平衡开始阶段可能有用	采样配分函数定义不佳，不能保证是 NPT 或 NPH，不应用于生产采样
Andersen	在运动方程中增加额外自由度，将系统耦合到假想压力浴（类似各向同性活塞作用）	采样正确的配分函数	适用于对配分函数采样要求高的场景	是各向同性的，不能对系统的一部分施加各向异性压力
Parrinello - Rahman	对 Andersen 压缩计的扩展，支持模拟盒大小和形状的各向异性缩放	具有 Andersen 压缩计的性质，增加了各向异性支持，适用于固体模拟（可模拟晶体格子中相变形状变化）等	在固体模拟等领域应用优势明显	并非所有模拟器都已实现该压缩计，限制了用户选择
Martyna - Tuckerman - Tobias - Klein (MTTK)	引入替代运动方程，正确采样较小系统配分函数	通常被视为对 Parrinello - Rahman 的改进，适用于小系统	为小系统的模拟提供了更准确的配分函数采样方法	未提及明显独特的局限性（与其他扩展自由度压缩计类似，受模拟器实现情况等限制）

Comment and share

截断能ENCUT的选取

Apr 29, 2025 in 从零开始的科研之旅, VASP

ENCUT是用来指定平面波基组能量截止值的参数

以上是wiki对ENCUT的描述，POTCAR文件中已经包含了元素的ENMAX和ENMIN，这也是互联网上对ENCUT的选择方案对通常是1.3倍的ENMAX

侯柱峰老师提供了另外一种方法并指出通过该方法选择的ENCUT通常能满足1.3倍ENMAX，供参考使用

截断能指定了用于波函数展开的平面波基组的截断能量，此能量越大则用来描述波函数的平面波基组越多，精度越高，但计算也越耗时。

可以通过计算测试来选择合适的ENCUT。

以下是用于测试的脚本：

#!/bin/sh
rm WAVECAR
for i in 150 200 250 300 350 400
do
cat > INCAR <<!
SYSTEM = Si-Diamond
ENCUT = $i
ISTART = 0 ; ICHARG = 2
ISMEAR = -5
PREC = Accurate
!
echo "ENCUT = $i eV" ; time vasp
E=`grep "TOTEN" OUTCAR | tail -1 | awk '{printf "%12.6f \n", $5 }'`
echo $i $E >>comment
done

计算后会得到comment文件，它列出了在不同ENCUT下计算得到的总能量

150 -11.900655
200 -11.938864
250 -11.944599
300 -11.945248
350 -11.945503
400 -11.945622

==总能量变化在0.001 eV即可==，因此在这个例子中ENCUT可以选择250

Comment and share

表面能的计算

Apr 29, 2025 in 从零开始的科研之旅, VASP

关于表面能

根据VASP官方手册Hands on Session III，表面能的定义如下，大师兄的文章也提到了这个定义( Nanoscale, 2017, 9, 13089–13094) $$ \sigma=\sigma^{unrelax}+E^{relax}=\dfrac{1}{2}(E_{surf}-N_{atoms}\cdot E_{bulk})+E^{relax} $$

$$ \gamma_s=\dfrac{1}{2A}(E_{surf}^{unrelax}-N\cdot E_{bulk})+\dfrac{1}{A}(E_{surf}^{relax}-E_{surf}^{unrelax}) $$

其中，弛豫能量被定义为： E^relax = E_surf^relax − E_surf^unrelax

手册的公式和大师兄文章中公式的区别在于表面积A，表面能应当除以表面积，所以后者更合理一些
E_bulk为体相中单个原子的能量，所以N ⋅ E_bulk为这个体相的能量，即为OUTCAR直接读出的能量

关于表面积A的计算

这部分有些难理解，涉及到一点高等数学和线性代数🤕

表面积可通过晶格基矢量的叉积确定，在POSCAR或CONTCAR文件中，晶格基矢量通常以三行的形式定义

对于面积的计算，可设一组不共线的基矢量为a₁ = (a₁₁, a₁₂, 0), a₂ = (a₂₁, a₂₂, 0)

下面这张图应该挺形象

a⃗ × b⃗ = (0, 0, a₁₁ × a₂₂ − a₁₂ × a₂₁)

面积为外积模乘以缩放银子实际面积 = s² ⋅ |a₁₁a₂₂ − a₁₂a₂₁| PS：对于三维的体积体积 = s³ ⋅ |a₁ ⋅ (a₂ × a₃)|

p(1x1)Cu(111)的表面能

前面已经算得差不多了

slab模型的能量信息

[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' OUTCAR
  energy  without entropy=      -13.96941183  energy(sigma->0) =      -13.97673168
  energy  without entropy=      -13.96945528  energy(sigma->0) =      -13.97676958
  energy  without entropy=      -13.97006413  energy(sigma->0) =      -13.97726241
  energy  without entropy=      -13.97020643  energy(sigma->0) =      -13.97735554
  energy  without entropy=      -13.97017969  energy(sigma->0) =      -13.97735013

slab模型优化后的结构信息

CONTCAR\(1\1\1)
   1.00000000000000
     2.5701000690000000    0.0000000000000000    0.0000000000000000
    -1.2850500345000000    2.2257719501000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
               4
Selective dynamics
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000   F   F   F
  0.6666700000000034  0.3333299999999966  0.0985399999999998   F   F   F
  0.3333311322281505  0.6666688677718494  0.1972481653092927   T   T   T
 -0.0000004145140272  0.0000004145140272  0.2951613337135311   T   T   T

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

bulk模型的能量信息

[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' .bulk-calc/OUTCAR
  energy  without entropy=      -14.90258332  energy(sigma->0) =      -14.90537746
  energy  without entropy=      -14.90568039  energy(sigma->0) =      -14.90844732
  energy  without entropy=      -14.90859467  energy(sigma->0) =      -14.91129762

A = 2.570100069 * 2.2257719501 − 0.0000000000000000 * (−1.2850500345) = 5.720456642530275 Å²
E_surf^unrelax = −13.9767 eV
E_bulk ⋅ N = −14.9112 eV
E_surface^relax = −13.9773 eV

所以表面能γ_s = 0.08158 eV/Å² ,又1 eV/Å²=16.0218 J/m²，所以表面能γ_s= 1.3070 J/m²

c(1x1)Cu(111)的表面能

之前大师兄提到过，对于面心立方最密堆积的形式，slab模型的建立可以使用原胞也可以使用单胞，那么接下来看看用单胞计算表面能吧

c(1x1)Cu(111)优化计算的能量信息：

[ctan@baifq-hpc141 conventional-slab-opt]$ grep '  without' OUTCAR
  energy  without entropy=      -55.94972046  energy(sigma->0) =      -55.94434499
  energy  without entropy=      -55.94997517  energy(sigma->0) =      -55.94458201
  energy  without entropy=      -55.95396000  energy(sigma->0) =      -55.94820485
  energy  without entropy=      -55.95519341  energy(sigma->0) =      -55.94928561
  energy  without entropy=      -55.95554732  energy(sigma->0) =      -55.94959634
  energy  without entropy=      -55.95713778  energy(sigma->0) =      -55.95095507
  energy  without entropy=      -55.95788514  energy(sigma->0) =      -55.95158140
  energy  without entropy=      -55.95893404  energy(sigma->0) =      -55.95242747
  energy  without entropy=      -55.95963672  energy(sigma->0) =      -55.95296199
  energy  without entropy=      -55.96031694  energy(sigma->0) =      -55.95340564
  energy  without entropy=      -55.96050719  energy(sigma->0) =      -55.95349958
  energy  without entropy=      -55.96053048  energy(sigma->0) =      -55.95316488
  energy  without entropy=      -55.96067995  energy(sigma->0) =      -55.95346716
  energy  without entropy=      -55.96069855  energy(sigma->0) =      -55.95356576
  energy  without entropy=      -55.96051601  energy(sigma->0) =      -55.95360396
  energy  without entropy=      -55.96050802  energy(sigma->0) =      -55.95361124
  energy  without entropy=      -55.96041586  energy(sigma->0) =      -55.95360050
  energy  without entropy=      -55.96018279  energy(sigma->0) =      -55.95347398

结构信息：

[ctan@baifq-hpc141 conventional-slab-opt]$ cat CONTCAR
CONTCAR\(1\1\1)
   1.00000000000000
     5.1402001381000000    0.0000000000000000    0.0000000000000000
    -2.5701000690000000    4.4515439000999999    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
              16
Direct
  0.0000006270807907 -0.0000006270807907  0.0011656342599333
  0.6666672048402060  0.3333327951597896  0.0988195745854847
  0.3333327951597896  0.6666672048402060  0.1968004254145135
 -0.0000006270807907  0.0000006270807907  0.2944543657400676
  0.0000006270807907  0.4999993729192098  0.0011656342599333
  0.6666672048402060  0.8333327951597940  0.0988195745854847
  0.3333327951597896  0.1666672048402087  0.1968004254145135
 -0.0000006270807907  0.5000006270807924  0.2944543657400676
  0.5000006270807924  0.4999993729192098  0.0011656342599333
  0.1666672048402087  0.8333327951597940  0.0988195745854847
  0.8333327951597940  0.1666672048402087  0.1968004254145135
  0.4999993729192098  0.5000006270807924  0.2944543657400676
  0.5000006270807924 -0.0000006270807907  0.0011656342599333
  0.1666672048402087  0.3333327951597896  0.0988195745854847
  0.8333327951597940  0.6666672048402060  0.1968004254145135
  0.4999993729192098  0.0000006270807907  0.2944543657400676

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

Cu(111)的体相能量：

[ctan@baifq-hpc141 conventional-slab-opt]$ grep '  without' .bulk-calc/OUTCAR
  energy  without entropy=      -14.90258332  energy(sigma->0) =      -14.90537746
  energy  without entropy=      -14.90568039  energy(sigma->0) =      -14.90844732
  energy  without entropy=      -14.90859467  energy(sigma->0) =      -14.91129762

表面积A=5.1402001381000000 * 4.4515439000999999 = 22.8818 Å²
E_surf^unrelax = -55.9443 eV
N ⋅ E_bulk = -59.6452 eV
- N指的是slab模型中的原子个数而非bulk中的原子个数
E_surf^relax = -55.9535 eV

计算得表面能γ_s=0.08047 eV/Å² = 1.2892 J/m²

这与p(1x1)Cu(111)计算所得的表面能有一定的差距，目前现阶段我还无法解决这个问题🤔😵‍💫

p(1x1)Ni(100)的表面能

这个例子来自于vasp官方的手册Hands on Session III

general:
SYSTEM = clean Ni(100) surface
ISTART = 0 ; ICHARG=2
ENCUT = 270
ISMEAR = 2 ; SIGMA = 0.2
spin:
ISPIN=2
MAGMOM = 5*1
dynamic:
IBRION = 1
NSW = 100
POTIM = 0.2

有个疑惑的点是手册中的INCAR没有加上偶极矫正，我在自己算的时候加上了偶极矫正

System = Nickel-slab-opt
ISMEAR = 1
ISTART = 0
ICHARG = 2
SIGMA = 0.05
ENCUT = 360
ALGO = Normal
ISPIN = 2
MAGMOM = 5*1
NSW = 200
NELM = 200
EDIFF = 1E-6
EDIFFG = -0.01
IDIPOL = 3
LDIPOL = .TRUE.
LWAVE = .FALSE.
LCHARG = .FALSE.
IBRION = 1
POTIM = 0.05

体相能量信息：

[ctan@baifq-hpc141 primitive-opt]$ grep '  without' .bulk/OUTCAR
  energy  without entropy=      -21.88865316  energy(sigma->0) =      -21.88759144
  energy  without entropy=      -21.88930785  energy(sigma->0) =      -21.88825056
  energy  without entropy=      -21.89056178  energy(sigma->0) =      -21.88952707

弛豫与未弛豫的slab模型的能量信息：

[ctan@baifq-hpc141 primitive-opt]$ grep '  without' OUTCAR
  energy  without entropy=      -25.97155386  energy(sigma->0) =      -25.97291371
  energy  without entropy=      -25.97233802  energy(sigma->0) =      -25.97366730
  energy  without entropy=      -25.97412653  energy(sigma->0) =      -25.97534340
  energy  without entropy=      -25.97511519  energy(sigma->0) =      -25.97602399
  energy  without entropy=      -25.97516432  energy(sigma->0) =      -25.97619350
  energy  without entropy=      -25.97516747  energy(sigma->0) =      -25.97619463
  energy  without entropy=      -25.97517140  energy(sigma->0) =      -25.97620156

优化后的结构：

CONTCAR\(1\0\0)
   1.00000000000000
     2.4846000671000001    0.0000000000000000    0.0000000000000000
    -1.2423000336000001    2.1517267763999999    0.0000000000000000
     0.0000000000000000    0.0000000000000000   18.7047996521000002
  Ni_pv/368cd815
               5
Direct
  0.6666685456044492  0.3333314543955508  0.0005454007373249
  0.0000003831810720 -0.0000003831810720  0.1078816295897802
  0.3333305745814666  0.6666694254185335  0.2169200598414882
  0.6666681017676376  0.3333318982323624  0.3259932620996842
  0.0000023948653782 -0.0000023948653782  0.4332496477317282

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

表面积A=2.4846000671000001*2.1517267763999999=5.346180493024307 eV
E_surf^unrelax = -25.97291371 eV
E_surf^relax = -25.97620156 eV
N ⋅ E_bulk = -21.88952707 eV

表面能 γ_s = 0.1293 eV/Å² = 2.0715 J/m²

小结

本节主要是联系表面能的计算
Ni(100)表面能的计算选用的赝势为Ni_sv赝势

Comment and share

POTCAR赝势的选取

Apr 29, 2025 in 从零开始的科研之旅, VASP

本文为一个查阅性的记录文章，也可以说是搬运，内容取自侯柱锋老师的使用VASP的个人经验手册

元素	普通赝势 (eV)	高硬赝势 (eV)	软赝势 (eV)	d 电子作为半芯态来处理赝势 (eV)
B	318	700 (B_h)	250 (B_s)	-
C	400	700 (C_h)	273 (C_s)	-
N	400	700 (N_h)	250 (N_s)	-
O	400	700 (O_h)	250 (O_s)	-
F	400	700 (F_h)	250 (F_s)	-
Al	240	295 (Al_h)	-	-
Si	245	380 (Si_h)	-	-
P	270	390 (P_h)	-	-
S	280	402 (S_h)	-	-
Cl	280	409 (Cl_h)	-	-
Ga	134	404 (Ga_h)	-	282 (Ga_d)
Ge	173	410 (Ge_h)	-	287 (Ge_d)
As	208	-	-	-
Se	211	-	-	-
Br	216	-	-	-
In	95	-	-	239 (In_d)
Sn	103	-	-	241 (Sn_d)
Sb	172	-	-	-
Te	174	-	-	-
I	175	-	-	-
Tl	95	-	-	239 (Tl_d)
Pb	98	-	-	237 (Pb_d)
Bi	105	-	-	242 (Bi_d)

元素	X——普通赝势 (eV)	X_h——硬赝势 (eV)	其他赝势 (eV)
H	250	700 (H_h)	-
Li	140	-	271 (Li_sv)
Be	300	-	308 (Be_sv)
Na	81	-	300 (Na_pv), 700 (Na_sv)
Mg	210	-	265 (Mg_pv)
K	-	-	150 (K_pv), 259 (K_sv)
Ca	-	-	150 (Ca_pv), 290 (Ca_sv)
Rb	-	-	121 (Rb_pv), 220 (Rb_sv)
Sr	-	-	226 (Sr_sv)
Cs	-	-	220 (Cs_sv)
Ba	-	-	187 (Ba_sv)

元素	普通赝势 (eV)	X_pv赝势 (eV)	X_sv赝势 (eV)
Sc	-	-	222 (Sc_sv)
Ti	178	222 (Ti_pv)	-
V	192	263 (V_pv)	-
Cr	227	265 (Cr_pv)	-
Mn	269	269 (Mn_pv)	-
Fe	267	293 (Fe_pv)	-
Co	267	-	-
Ni	269	367 (Ni_pv)	-
Cu	273	368 (Cu_pv)	-
Zn	276	-	-
Y	-	-	211 (Y_sv)
Zr	-	-	229 (Zr_sv)
Nb	-	207 (Nb_pv)	-
Mo	224	224 (Mo_pv)	-
Tc	228	228 (Tc_pv)	-
Ru	213	230 (Ru_pv)	-
Rh	228	271 (Rh_pv)	-
Pd	250	350 (Pd_pv)	-
Ag	249	-	-
Cd	274	-	-
Hf	220	220 (Hf_pv)	-
Ta	223	223 (Ta_pv)	-
W	223	223 (W_pv)	-
Re	226	226 (Re_pv)	-
Os	228	228 (Os_pv)	-
Ir	210	-	-
Pt	230	-	-
Au	229	-	-
Hg	233	-	-

Comment and share

以p(1x1)Cu(111)为例分析表面弛豫的结果

Apr 27, 2025 in 从零开始的科研之旅, VASP

本文其实没啥新的内容，只是对上两节的一个总结分析

在表面弛豫前（slab模型优化计算前）

SYSTEM = p(1x1)Cu(111)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Selective Dynamics
Direct
     0.000000000         0.000000000         0.000000000 F F F
     0.666670000         0.333330000         0.098540000 F F F
     0.333330000         0.666670000         0.197080000 T T T
     0.000000000         0.000000000         0.295620000 T T T

优化计算后

SYSTEM = p(1x1)Cu(111)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Selective Dynamics
Direct
     0.000000000         0.000000000         0.000000000 F F F
     0.666670000         0.333330000         0.098540000 F F F
     0.333330000         0.666670000         0.197080000 T T T
     0.000000000         0.000000000         0.295620000 T T T
[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat CONTCAR
CONTCAR\(1\1\1)
   1.00000000000000
     2.5701000690000000    0.0000000000000000    0.0000000000000000
    -1.2850500345000000    2.2257719501000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
               4
Selective dynamics
Direct
  0.0000000000000000  0.0000000000000000  0.0000000000000000   F   F   F
  0.6666700000000034  0.3333299999999966  0.0985399999999998   F   F   F
  0.3333311322281505  0.6666688677718494  0.1972481653092927   T   T   T
 -0.0000004145140272  0.0000004145140272  0.2951613337135311   T   T   T

  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00
  0.00000000E+00  0.00000000E+00  0.00000000E+00

工作站离线状态下numpy on python2地安装

这两个文件中的坐标都是分数坐标，不太利于分析位置变化，需要将分数坐标转化为笛卡尔坐标，使用大师兄的脚本来进行转换，

用的是python2

#!/usr/bin/env python
# -*- coding: utf-8 -*-
#Convert direc coordiation to cartesian Writen By Qiang and Xu Nan
"""
The right input way to run this script is : 
      
python dire2cart.py CONTCAR 3 

CONTCAR will be converted from Direct to Cartesian and 3 layers will be fixed.  

If you do not want to fix layers and keep the same as before, run command:
    
python dire2cat.py CONTCAR    

"""

import sys
import os
import numpy as np 

print 'Please read the head part of this script and get more information!'
print """
###################################
#                                 #
#for VASP 5.2 or higher versions  #
#                                 #
###################################
"""

if len(sys.argv) <= 1:
    print('\n' + ' Warning! ' * 3 + '\n')
    print('You did not select the inputfile to be converted. \n By defalut, we are going to convert your CONTCAR.\n')
    
    if not os.path.isfile("POSCAR") and not os.path.isfile("CONTCAR"):
        print("Error:" * 3 + "\n Can not find neither POSCAR nor CONTCAR!!\n")
        exit()
    else:
        if os.path.isfile("CONTCAR"):
            script = sys.argv[0]
            file_to_be_converted = "CONTCAR"
            print "\n Convertion starts........\n"
    
        else:
            script = sys.argv[0]
            file_to_be_converted = "POSCAR"
            print "\n There is no CONTCAR in your directory. \n \n POSCAR Convertion starts........\n"
            
if len(sys.argv) == 2:
    print("\n%s Convertion starts......" %sys.argv[1]) 
    script, file_to_be_converted = sys.argv[:2]
else:
    print("\n%s Convertion starts......\n" %sys.argv[1]) 
    script, file_to_be_converted, fixedlayer = sys.argv[:3]
    fixedlayer = int(fixedlayer)

def get_infor():
    f = open(file_to_be_converted, 'r')
    lines = f.readlines()
    f.close()
    num_atoms = sum([int(x) for x in lines[6].split()])
    if lines[7][0]  == 'S' or lines[7][0]  == 's':  # # With Selected T T T, coordination starts from line 9
        start_num = 9
        if lines[8][0]  == 'D' or lines[8][0]  == 'd':
            is_direct = True 
        else:
            is_direct = False
    else: 
        start_num = 8
        print '----------------------------------------------------'
        print 'Pay Attetion! There is no TTT in  %s   ' %(file_to_be_converted)
        print '----------------------------------------------------'
        if  lines[7][0]  == 'D' or lines[7][0]  == 'd':
            is_direct = True 
        else:
            is_direct = False
    a = []
    b = []
    c = []
    if is_direct:         
        for i in np.arange(2,5):
            line = [float(i) for i in lines[i].strip().split()]
            a.append(line[0])
            b.append(line[1])
            c.append(line[2])
        vector = np.array([a,b,c])
    else: 
        vector = np.array([[1, 0 , 0], [0, 1, 0], [0, 0, 1]])
    return vector, lines, start_num, num_atoms, is_direct

def determinelayers(z_cartesian):
    threshold = 0.5  
    seq = sorted(z_cartesian)
    min = seq[0]
    layerscount = {}
    sets = [min]
    for j in range(len(seq)):
        if abs(seq[j]-min) >= threshold:
            min = seq[j]
            sets.append(min)

    for i in range(1,len(sets)+1):
        layerscount[i] = []            
        for k in range(len(z_cartesian)):   
            if abs(z_cartesian[k]-sets[i-1]) <= threshold:
                layerscount[i].append(k)

    return layerscount

def convert():
    x_cartesian = []
    y_cartesian = []
    z_cartesian = []
    tf = []
    for i in range(start_num, num_atoms + start_num):
        line_data =  [float(ele) for ele in lines[i].split()[0:3]]
        line_data = np.array([line_data])
        x, y, z = [sum(k) for k in line_data * vector ]     
        x_cartesian.append(x)
        y_cartesian.append(y)
        z_cartesian.append(z)
         
        if start_num == 9 : # if  T T T exist, the start_num will be 9        
            tf.append((lines[i].split()[3:]))
        else:
            tf.append(' ')   # if there is no T T T, use space instead. 

    layerscount =determinelayers(z_cartesian)
    
    file_out = open(file_to_be_converted+'_C', 'w')
    for i in range(0,7):
        file_out.write(lines[i].rstrip() + '\n')  # first 7 lines are kept the same 
    
    print '\n Find %s layers!' * 3 %(len(layerscount), len(layerscount), len(layerscount))
    
    if len(sys.argv) >= 3:  # This means that the nuber for fixing layers is given by the user.
        file_out.write('Selective\n')  
        file_out.write('Cartesian' + '\n') 
        for i in range(1,len(layerscount)+1):
            if i <= fixedlayer: 
                for j in layerscount[i]:
                    tf[j] = ['F','F','F']
            else:
                for k in layerscount[i]:
                    tf[k] = ['T','T','T']                    
    else:
        if start_num == 9:  # 9 means there are T T T or F F F in the file
            file_out.write('Selective\n')
            file_out.write('Cartesian' + '\n')        
        else:
            file_out.write('Cartesian' + '\n') 
    
    for i in range(0,len(x_cartesian)):
        file_out.write("\t%+-3.10f   %+-3.10f   %+-3.10f  %s\n" %(x_cartesian[i], y_cartesian[i], z_cartesian[i], ' '.join(tf[i])))
        
    file_out.close()
    
vector, lines, start_num, num_atoms, is_direct = get_infor()

if is_direct : 
    print "\n%s has Direct Coordinates, Contersion starts.... "  %(file_to_be_converted)
    convert()
else:
    print "\n%s has Cartesian Coordinates Already! We are going to fix layers only."  %(file_to_be_converted)
    convert()
    
print '-----------------------------------------------------\n'
print '\n %s with Cartesian Coordiates is named as  %s_C\n' %(file_to_be_converted, file_to_be_converted)
print '-----------------------------------------------------\n'

脚本会调用到numpy库，然而集群上并没有连接互联网，所以没有办法通过pip2 install numpy来安装，只能手动上传.whl文件或者源码编译等。

进入python2环境，输入import pip;print(pip.pep425tags.get_supported())来查看支持的.whl文件版本

[ctan@baifq-hpc141 primitive-slab-opt-1]$ python2
Python 2.7.18 (default, Oct 23 2023, 20:01:03)
[GCC 8.5.0 20210514 (Red Hat 8.5.0-18)] on linux2
Type "help", "copyright", "credits" or "license" for more information.
>>> import pip;print(pip.pep425tags.get_supported())
[('cp27', 'cp27mu', 'manylinux1_x86_64'), ('cp27', 'cp27mu', 'linux_x86_64'), ('cp27', 'none', 'manylinux1_x86_64'), ('cp27', 'none', 'linux_x86_64'), ('py2', 'none', 'manylinux1_x86_64'), ('py2', 'none', 'linux_x86_64'), ('cp27', 'none', 'any'), ('cp2', 'none', 'any'), ('py27', 'none', 'any'), ('py2', 'none', 'any'), ('py26', 'none', 'any'), ('py25', 'none', 'any'), ('py24', 'none', 'any'), ('py23', 'none', 'any'), ('py22', 'none', 'any'), ('py21', 'none', 'any'), ('py20', 'none', 'any')]

这里下载的是numpy-1.16.5-cp27-cp27mu-manylinux1_x86_64.whl（不可下载过高版本的numpy，它支持python2的最后一个版本似乎是1.16.6）

然后安装，由于我不在sudoer中，没有root权限，所以需要手动指定安装目录

1	pip2 install --prefix=~/python2/numpy_installation numpy-1.16.5-cp27-cp27mu-manylinux1_x86_64.whl

然后添加环境变量，注意numpy库在/lib64/python2.7/site-packages目录中的

1	export PYTHONPATH=/home/ctan/python2/numpy_installation/lib64/python2.7/site-packages:$PYTHONPATH

ok啦，以上是题外话，安装完numpy之后就可以运行脚本了，将POSCAR和CONTCAR的坐标从分数坐标转换为笛卡尔坐标

POSCAR

[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat POSCAR_C
CONTCAR\(1\1\1)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Selective
Cartesian
        +0.0000000000   +0.0000000000   +0.0000000000  F F F
        +1.2850628850   +0.7419165641   +2.0984487770  F F F
        -0.0000128505   +1.4838553860   +4.1968975541  T T T
        +0.0000000000   +0.0000000000   +6.2953463311  T T T

CONTCAR_C

[ctan@baifq-hpc141 primitive-slab-opt-1]$ cat CONTCAR_C
CONTCAR\(1\1\1)
   1.00000000000000
     2.5701000690000000    0.0000000000000000    0.0000000000000000
    -1.2850500345000000    2.2257719501000000    0.0000000000000000
     0.0000000000000000    0.0000000000000000   21.2954006195000005
  Cu/0e71e558f37
               4
Selective
Cartesian
        +0.0000000000   +0.0000000000   +0.0000000000  F F F
        +1.2850628850   +0.7419165641   +2.0984487770  F F F
        -0.0000084856   +1.4838528659   +4.2004787017  T T T
        -0.0000015980   +0.0000009226   +6.2855788488  T T T

Alright, I’m back. 突然想起昨天球磨的样品今天结束了，连忙赶回实验室取出来装进样品管里等待封管。

By the way，我真是和鱼一样，记忆只有七秒、、、😅。师兄教的球磨罐和球珠的洗涤方法忘记下来了，今天不出意外地给忘了😂。汲取教训，每次出实验室后都更加努力地总结和复盘，仔细回想每一个实验细节

优化前后对比

Alright，扯远了，分析一下顶层和亚顶层优化前后地结构变化吧：

POSCAR: (+6.2953463311) - (+4.1968975541) = 2.0984487769999998
CONTCAR: (+6.2855788488) - (+4.2004787017) = 2.0851001471000004

前后变化了：(2.0851001471000004-2.0984487769999998)/2.0984487769999998 = -0.006361189296735172

可以看到，表层原子向体相收缩

能量变化

[ctan@baifq-hpc141 primitive-slab-opt-1]$ grep '  without' OUTCAR
  energy  without entropy=      -13.96941183  energy(sigma->0) =      -13.97673168
  energy  without entropy=      -13.96945528  energy(sigma->0) =      -13.97676958
  energy  without entropy=      -13.97006413  energy(sigma->0) =      -13.97726241
  energy  without entropy=      -13.97020643  energy(sigma->0) =      -13.97735554
  energy  without entropy=      -13.97017969  energy(sigma->0) =      -13.97735013

能量从-13.97673168变化为-13.97735013，降低了0.000618449999999271 eV，即弛豫过程是放热，刚切开的表面不稳定，表层原子向体相收缩后，体系能量降低，变得更稳定

小Tips：可以使用python环境便捷地进行基本运算

Comment and share

p(1x1)Cu(111)的优化计算

Apr 27, 2025 in 从零开始的科研之旅, VASP

`POSCAR`的修改

与单点计算不同，优化计算可能涉及到离子步迭代、晶格常数的改变。但Cu(111)晶面的优化计算并==不只是==在INCAR里添加参数IBRION，优化计算确实允许原子弛豫，但在slab模型中，只允许表面的两层原子允许弛豫，而下面的那些原子则直接固定住

这样做的物理意义在于：在真实的环境下，催化剂体相看做是不变的，只有表面的原子参与催化反应

前面在乙醇分子的振动频率分析里面已经学过了，添加以下选择性动力学标记Slective Dynamics即可

PS：看分数坐标即可知道哪些是表层哪些是体相

CONTCAR\(1\1\1)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Selective Dynamics
Direct
     0.000000000         0.000000000         0.000000000 F F F
     0.666670000         0.333330000         0.098540000 F F F
     0.333330000         0.666670000         0.197080000 T T T
     0.000000000         0.000000000         0.295620000 T T T

`INCAR`的修改

现在是优化计算，就得加上算法了，大师兄用的IBRION = 2，那我就对比一下，用IBRION = 1试一下吧

`IBRION = 1`

官方wiki说准牛顿迭代法对POTIM比较敏感，那我就调小一点吧

SYSTEM = Cu-concentional-slab-statistic
ISMEAR = 1
SIGMA = 0.5
ALGO = Normal
ENCUT = 700
EDIFF = 1E-6
LDIPOL = .TRUE.
IDIPOL = 3
NCORE = 20
NSW = 500
NELM = 200
IBRION = 1
POTIM = 0.05

成功收敛，迭代情况

[ctan@baifq-hpc141 primitive-slab-opt]$ cat OSZICAR
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1     0.380094340720E+03    0.38009E+03   -0.23566E+04  2548   0.269E+03
DAV:   2     0.163549324393E+02   -0.36374E+03   -0.35588E+03  2548   0.438E+02
DAV:   3    -0.154087935417E+02   -0.31764E+02   -0.30239E+02  2816   0.186E+02
DAV:   4    -0.171245615883E+02   -0.17158E+01   -0.16227E+01  2856   0.463E+01
DAV:   5    -0.171794690835E+02   -0.54907E-01   -0.54864E-01  2954   0.690E+00    0.146E+01
DAV:   6    -0.156158003013E+02    0.15637E+01   -0.42034E+01  2712   0.133E+02    0.829E+00
DAV:   7    -0.140078760986E+02    0.16079E+01   -0.79564E+00  2992   0.596E+01    0.790E-01
DAV:   8    -0.140167702165E+02   -0.88941E-02   -0.53048E-01  3968   0.226E+01    0.122E+00
DAV:   9    -0.140083290175E+02    0.84412E-02   -0.45257E-01  3848   0.291E+01    0.104E+00
DAV:  10    -0.139869448329E+02    0.21384E-01   -0.20408E-01  4800   0.187E+01    0.615E-01
DAV:  11    -0.139810722283E+02    0.58726E-02   -0.60166E-02  4016   0.780E+00    0.113E-01
DAV:  12    -0.139806486192E+02    0.42361E-03   -0.10641E-02  3786   0.272E+00    0.980E-02
DAV:  13    -0.139810968626E+02   -0.44824E-03   -0.87951E-03  4868   0.367E+00    0.844E-02
DAV:  14    -0.139803135502E+02    0.78331E-03   -0.42901E-03  4920   0.238E+00    0.308E-02
DAV:  15    -0.139802543403E+02    0.59210E-04   -0.32281E-04  3540   0.324E-01    0.182E-02
DAV:  16    -0.139803209356E+02   -0.66595E-04   -0.83429E-04  4558   0.937E-01    0.301E-02
DAV:  17    -0.139803228923E+02   -0.19567E-05   -0.15361E-04  3074   0.264E-01    0.168E-02
DAV:  18    -0.139803417113E+02   -0.18819E-04   -0.36368E-04  4750   0.744E-01    0.954E-03
DAV:  19    -0.139803412500E+02    0.46126E-06   -0.10659E-04  2728   0.404E-01    0.517E-03
DAV:  20    -0.139803676358E+02   -0.26386E-04   -0.80023E-05  3638   0.357E-01    0.818E-03
DAV:  21    -0.139803787022E+02   -0.11066E-04   -0.13797E-05  3336   0.920E-02    0.101E-02
DAV:  22    -0.139803765945E+02    0.21077E-05   -0.39937E-05  2614   0.265E-01    0.222E-03
DAV:  23    -0.139803888096E+02   -0.12215E-04   -0.60877E-05  2528   0.297E-01    0.725E-03
DAV:  24    -0.139803907613E+02   -0.19517E-05   -0.71066E-06  3156   0.553E-02    0.611E-03
DAV:  25    -0.139803916036E+02   -0.84239E-06   -0.13171E-06  1634   0.198E-02
   1 F= -.13980392E+02 E0= -.13976732E+02  d E =-.139804E+02
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139807173282E+02   -0.32657E-03    0.14112E+02  2602   0.612E+00    0.467E-03
DAV:   2    -0.139804390721E+02    0.27826E-03   -0.16468E-04  4156   0.597E-01    0.165E-02
DAV:   3    -0.139804261532E+02    0.12919E-04   -0.10604E-04  3756   0.481E-01    0.167E-03
DAV:   4    -0.139804267349E+02   -0.58173E-06   -0.46237E-06  2500   0.225E-02
   2 F= -.13980427E+02 E0= -.13976770E+02  d E =-.351312E-04
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139804330907E+02   -0.69375E-05    0.16389E+02  2632   0.107E+01    0.452E-02
DAV:   2    -0.139809280760E+02   -0.49499E-03   -0.57833E-03  3892   0.164E+00    0.586E-02
DAV:   3    -0.139808976764E+02    0.30400E-04   -0.10914E-03  4308   0.154E+00    0.348E-02
DAV:   4    -0.139808745784E+02    0.23098E-04   -0.21714E-04  4572   0.607E-01    0.278E-02
DAV:   5    -0.139808671461E+02    0.74323E-05   -0.56092E-05  4164   0.315E-01    0.187E-02
DAV:   6    -0.139808630630E+02    0.40832E-05   -0.44161E-05  4470   0.289E-01    0.698E-03
DAV:   7    -0.139808615417E+02    0.15213E-05   -0.76178E-06  1946   0.121E-01    0.290E-03
DAV:   8    -0.139808615453E+02   -0.35725E-08   -0.25204E-06  1800   0.234E-02
   3 F= -.13980862E+02 E0= -.13977262E+02  d E =-.434810E-03
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139808810490E+02   -0.19507E-04    0.54240E+01  2632   0.376E+00    0.136E-02
DAV:   2    -0.139809351461E+02   -0.54097E-04   -0.54142E-04  3890   0.509E-01    0.173E-02
DAV:   3    -0.139809327068E+02    0.24393E-05   -0.10930E-04  4292   0.484E-01    0.104E-02
DAV:   4    -0.139809303964E+02    0.23104E-05   -0.18479E-05  2078   0.186E-01    0.828E-03
DAV:   5    -0.139809300877E+02    0.30864E-06   -0.87511E-06  2780   0.122E-01
   4 F= -.13980930E+02 E0= -.13977356E+02  d E =-.685425E-04
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139808203045E+02    0.11009E-03    0.13819E+02  3808   0.576E+00    0.330E-03
DAV:   2    -0.139808743367E+02   -0.54032E-04    0.77025E+00  2180   0.218E+00    0.818E-03
DAV:   3    -0.139809353212E+02   -0.60985E-04   -0.25088E-05  2014   0.239E-01    0.103E-03
DAV:   4    -0.139809353447E+02   -0.23522E-07   -0.62565E-07  1516   0.240E-02
   5 F= -.13980935E+02 E0= -.13977350E+02  d E =-.525697E-05

能量：

[ctan@baifq-hpc141 primitive-slab-opt]$ grep '  without' OUTCAR
  energy  without entropy=      -13.96941183  energy(sigma->0) =      -13.97673168
  energy  without entropy=      -13.96945528  energy(sigma->0) =      -13.97676958
  energy  without entropy=      -13.97006413  energy(sigma->0) =      -13.97726241
  energy  without entropy=      -13.97020643  energy(sigma->0) =      -13.97735554
  energy  without entropy=      -13.97017969  energy(sigma->0) =      -13.97735013

计算时间：

General timing and accounting informations for this job:
========================================================

                 Total CPU time used (sec):       53.453
                           User time (sec):       49.556
                         System time (sec):        3.896
                        Elapsed time (sec):       55.466

                  Maximum memory used (kb):      134456.
                  Average memory used (kb):          N/A

                         Minor page faults:        12550
                         Major page faults:          686
                Voluntary context switches:         1587

`IBRION =2`

SYSTEM = Cu-concentional-slab-statistic
ISMEAR = 1
SIGMA = 0.5
ALGO = Normal
ENCUT = 700
EDIFF = 1E-6
LDIPOL = .TRUE.
IDIPOL = 3
NCORE = 20
NSW = 500
NELM = 200
IBRION = 2
POTIM = 0.1

成功收敛，迭代情况：

       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1     0.380094340720E+03    0.38009E+03   -0.23566E+04  2548   0.269E+03
DAV:   2     0.163549324393E+02   -0.36374E+03   -0.35588E+03  2548   0.438E+02
DAV:   3    -0.154087935417E+02   -0.31764E+02   -0.30239E+02  2816   0.186E+02
DAV:   4    -0.171245615883E+02   -0.17158E+01   -0.16227E+01  2856   0.463E+01
DAV:   5    -0.171794690835E+02   -0.54907E-01   -0.54864E-01  2954   0.690E+00    0.146E+01
DAV:   6    -0.156158003013E+02    0.15637E+01   -0.42034E+01  2712   0.133E+02    0.829E+00
DAV:   7    -0.140078760986E+02    0.16079E+01   -0.79564E+00  2992   0.596E+01    0.790E-01
DAV:   8    -0.140167702165E+02   -0.88941E-02   -0.53048E-01  3968   0.226E+01    0.122E+00
DAV:   9    -0.140083290175E+02    0.84412E-02   -0.45257E-01  3848   0.291E+01    0.104E+00
DAV:  10    -0.139869448329E+02    0.21384E-01   -0.20408E-01  4800   0.187E+01    0.615E-01
DAV:  11    -0.139810722283E+02    0.58726E-02   -0.60166E-02  4016   0.780E+00    0.113E-01
DAV:  12    -0.139806486192E+02    0.42361E-03   -0.10641E-02  3786   0.272E+00    0.980E-02
DAV:  13    -0.139810968626E+02   -0.44824E-03   -0.87951E-03  4868   0.367E+00    0.844E-02
DAV:  14    -0.139803135502E+02    0.78331E-03   -0.42901E-03  4920   0.238E+00    0.308E-02
DAV:  15    -0.139802543403E+02    0.59210E-04   -0.32281E-04  3540   0.324E-01    0.182E-02
DAV:  16    -0.139803209356E+02   -0.66595E-04   -0.83429E-04  4558   0.937E-01    0.301E-02
DAV:  17    -0.139803228923E+02   -0.19567E-05   -0.15361E-04  3074   0.264E-01    0.168E-02
DAV:  18    -0.139803417113E+02   -0.18819E-04   -0.36368E-04  4750   0.744E-01    0.954E-03
DAV:  19    -0.139803412500E+02    0.46126E-06   -0.10659E-04  2728   0.404E-01    0.517E-03
DAV:  20    -0.139803676358E+02   -0.26386E-04   -0.80023E-05  3638   0.357E-01    0.818E-03
DAV:  21    -0.139803787022E+02   -0.11066E-04   -0.13797E-05  3336   0.920E-02    0.101E-02
DAV:  22    -0.139803765945E+02    0.21077E-05   -0.39937E-05  2614   0.265E-01    0.222E-03
DAV:  23    -0.139803888096E+02   -0.12215E-04   -0.60877E-05  2528   0.297E-01    0.725E-03
DAV:  24    -0.139803907613E+02   -0.19517E-05   -0.71066E-06  3156   0.553E-02    0.611E-03
DAV:  25    -0.139803916036E+02   -0.84239E-06   -0.13171E-06  1634   0.198E-02
   1 F= -.13980392E+02 E0= -.13976732E+02  d E =-.139804E+02
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139807447433E+02   -0.35398E-03    0.14111E+02  2590   0.632E+00    0.528E-03
DAV:   2    -0.139804678790E+02    0.27686E-03   -0.17889E-04  4964   0.563E-01    0.155E-02
DAV:   3    -0.139804572701E+02    0.10609E-04   -0.99595E-05  4260   0.455E-01    0.347E-03
DAV:   4    -0.139804576704E+02   -0.40025E-06   -0.38794E-06  2444   0.252E-02
   2 F= -.13980458E+02 E0= -.13976804E+02  d E =-.660667E-04
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139806303531E+02   -0.17308E-03    0.10590E+01  2646   0.165E+00    0.137E-02
DAV:   2    -0.139806282211E+02    0.21320E-05   -0.51586E-04  3892   0.441E-01    0.161E-02
DAV:   3    -0.139806267938E+02    0.14273E-05   -0.86921E-05  4092   0.433E-01    0.112E-02
DAV:   4    -0.139806245851E+02    0.22087E-05   -0.18757E-05  2156   0.193E-01    0.871E-03
DAV:   5    -0.139806244737E+02    0.11138E-06   -0.11052E-05  3344   0.135E-01    0.639E-03
DAV:   6    -0.139806241801E+02    0.29361E-06   -0.44497E-06  1850   0.106E-01
   3 F= -.13980624E+02 E0= -.13976988E+02  d E =-.232576E-03
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139806399034E+02   -0.15430E-04    0.10673E+02  2632   0.655E+00    0.273E-02
DAV:   2    -0.139808756284E+02   -0.23573E-03   -0.22093E-03  3886   0.102E+00    0.353E-02
DAV:   3    -0.139808651090E+02    0.10519E-04   -0.43844E-04  4364   0.978E-01    0.208E-02
DAV:   4    -0.139808558492E+02    0.92598E-05   -0.82890E-05  3866   0.380E-01    0.163E-02
DAV:   5    -0.139808537696E+02    0.20796E-05   -0.22746E-05  3796   0.195E-01    0.113E-02
DAV:   6    -0.139808525623E+02    0.12073E-05   -0.15366E-05  3678   0.171E-01    0.403E-03
DAV:   7    -0.139808519218E+02    0.64052E-06   -0.31561E-06  1528   0.772E-02
   4 F= -.13980852E+02 E0= -.13977253E+02  d E =-.460318E-03
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139800417341E+02    0.81083E-03    0.23560E+01  2630   0.588E+00    0.573E-02
DAV:   2    -0.139809885034E+02   -0.94677E-03   -0.91190E-03  3890   0.208E+00    0.745E-02
DAV:   3    -0.139809389586E+02    0.49545E-04   -0.17575E-03  4372   0.196E+00    0.428E-02
DAV:   4    -0.139809007116E+02    0.38247E-04   -0.35205E-04  4684   0.779E-01    0.339E-02
DAV:   5    -0.139808904102E+02    0.10301E-04   -0.81050E-05  4124   0.377E-01    0.239E-02
DAV:   6    -0.139808824251E+02    0.79850E-05   -0.67930E-05  4464   0.342E-01    0.861E-03
DAV:   7    -0.139808804209E+02    0.20042E-05   -0.13063E-05  2516   0.156E-01    0.364E-03
DAV:   8    -0.139808803022E+02    0.11871E-06   -0.29028E-06  1874   0.341E-02
   5 F= -.13980880E+02 E0= -.13977356E+02  d E =-.488699E-03
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139807672910E+02    0.11313E-03    0.15468E+02  2630   0.858E+00    0.205E-02
DAV:   2    -0.139809519135E+02   -0.18462E-03   -0.16290E-03  3860   0.606E-01    0.224E-02
DAV:   3    -0.139809495045E+02    0.24089E-05   -0.16101E-04  3760   0.572E-01    0.195E-02
DAV:   4    -0.139809446950E+02    0.48095E-05   -0.34715E-05  3338   0.266E-01    0.149E-02
DAV:   5    -0.139809410051E+02    0.36899E-05   -0.17275E-05  3888   0.130E-01    0.678E-03
DAV:   6    -0.139809401472E+02    0.85792E-06   -0.69054E-06  2716   0.164E-01
   6 F= -.13980940E+02 E0= -.13977383E+02  d E =-.548544E-03
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139808959432E+02    0.45062E-04    0.81868E+01  2640   0.491E+00    0.921E-03
DAV:   2    -0.139810234583E+02   -0.12752E-03   -0.40347E-04  3868   0.335E-01    0.124E-02
DAV:   3    -0.139810222522E+02    0.12061E-05   -0.37434E-05  3890   0.257E-01    0.934E-03
DAV:   4    -0.139810216716E+02    0.58055E-06   -0.56497E-06  1910   0.105E-01
   7 F= -.13981022E+02 E0= -.13977454E+02  d E =-.815244E-04
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139809168510E+02    0.10540E-03    0.16544E+02  2640   0.100E+01    0.295E-02
DAV:   2    -0.139812498696E+02   -0.33302E-03   -0.35344E-03  3874   0.942E-01    0.353E-02
DAV:   3    -0.139812400259E+02    0.98437E-05   -0.31428E-04  3728   0.768E-01    0.285E-02
DAV:   4    -0.139812339348E+02    0.60911E-05   -0.66313E-05  3954   0.364E-01    0.240E-02
DAV:   5    -0.139812229797E+02    0.10955E-04   -0.53085E-05  4306   0.249E-01    0.134E-02
DAV:   6    -0.139812204144E+02    0.25653E-05   -0.20199E-05  3228   0.260E-01    0.363E-03
DAV:   7    -0.139812202373E+02    0.17704E-06   -0.19687E-06  1482   0.408E-02
   8 F= -.13981220E+02 E0= -.13977623E+02  d E =-.280090E-03
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139799374972E+02    0.12829E-02    0.13429E+02  2638   0.119E+01    0.609E-02
DAV:   2    -0.139815536956E+02   -0.16162E-02   -0.15321E-02  3874   0.237E+00    0.834E-02
DAV:   3    -0.139814840776E+02    0.69618E-04   -0.19373E-03  4340   0.194E+00    0.553E-02
DAV:   4    -0.139814504163E+02    0.33661E-04   -0.33985E-04  4524   0.775E-01    0.477E-02
DAV:   5    -0.139814168287E+02    0.33588E-04   -0.15051E-04  4446   0.480E-01    0.325E-02
DAV:   6    -0.139813975448E+02    0.19284E-04   -0.12765E-04  4566   0.480E-01    0.853E-03
DAV:   7    -0.139813962132E+02    0.13316E-05   -0.15668E-05  3464   0.152E-01    0.507E-03
DAV:   8    -0.139813960169E+02    0.19635E-06   -0.22934E-06  1728   0.311E-02
   9 F= -.13981396E+02 E0= -.13977735E+02  d E =-.455870E-03
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139815084243E+02   -0.11221E-03    0.83969E+01  2670   0.421E+00    0.434E-03
DAV:   2    -0.139813984003E+02    0.11002E-03   -0.73291E-05  3840   0.104E-01    0.437E-03
DAV:   3    -0.139813981891E+02    0.21124E-06   -0.11427E-06  1510   0.490E-02
  10 F= -.13981398E+02 E0= -.13977729E+02  d E =-.217219E-05
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.139815087603E+02   -0.11036E-03    0.15888E+02  2670   0.681E+00    0.323E-03
DAV:   2    -0.139813984696E+02    0.11029E-03   -0.26471E-05  2428   0.210E-01    0.588E-03
DAV:   3    -0.139813976417E+02    0.82788E-06   -0.27049E-05  3780   0.204E-01    0.141E-03
DAV:   4    -0.139813978006E+02   -0.15885E-06   -0.17146E-06  1590   0.262E-02
  11 F= -.13981398E+02 E0= -.13977732E+02  d E =-.178370E-05

计算时间：

General timing and accounting informations for this job:
========================================================

                 Total CPU time used (sec):       93.113
                           User time (sec):       87.076
                         System time (sec):        6.036
                        Elapsed time (sec):       94.073

                  Maximum memory used (kb):      126136.
                  Average memory used (kb):          N/A

                         Minor page faults:        11497
                         Major page faults:          687
                Voluntary context switches:         1456

能量信息：

[ctan@baifq-hpc141 primitive-slab-opt-2]$ grep '  without' OUTCAR
  energy  without entropy=      -13.96941183  energy(sigma->0) =      -13.97673168
  energy  without entropy=      -13.96949582  energy(sigma->0) =      -13.97680372
  energy  without entropy=      -13.96971572  energy(sigma->0) =      -13.97698803
  energy  without entropy=      -13.97005607  energy(sigma->0) =      -13.97725330
  energy  without entropy=      -13.97030617  energy(sigma->0) =      -13.97735559
  energy  without entropy=      -13.97026927  energy(sigma->0) =      -13.97738319
  energy  without entropy=      -13.97031853  energy(sigma->0) =      -13.97745396
  energy  without entropy=      -13.97042910  energy(sigma->0) =      -13.97762319
  energy  without entropy=      -13.97041285  energy(sigma->0) =      -13.97773496
  energy  without entropy=      -13.97039032  energy(sigma->0) =      -13.97772890
  energy  without entropy=      -13.97040171  energy(sigma->0) =      -13.97773244

对比一下发现，正如wiki所说，IBRION = 1适合初始结构较好的情况，IBRION - 2的适用性更广，更可靠，尽管IBRION = 2 的POTIM更大，但所用时间更长。二者最终的能量相差不大

Comment and share

Cu(111) slab模型的单点计算

Apr 23, 2025 in 从零开始的科研之旅, VASP

首先对Cu单胞进行结构优化（bulk计算），然后使用优化的结构进行切割，再进行计算（由于晶胞已经优化过，所以本次计算其实是单点计算）。大师兄提到：

Cu(111) slab 模型（POSCAR）的制作流程（复习一下）。注意：前面我们用的是Conventional cell，下面用的是Primitive Cell。一般来说用Conventional cell， FCC的金属可以用Primitive Cell 但对于其他体系，通过Primitive Cell切出来表面模型有问题

我准备两种都验证一下。

使用Conventional Cell建立slab模型

切割结构

将bulk计算完毕的CONTCAR文件用Vesta打开，再转换为.cif文件使用Materials Studio打开

Step1：Build-Surface-Cleave Surface

Angstrom是单位埃，需要切割的晶面是(111)，相应的平面方程为x + y + z − 1 = 0，使用点到直线的距离公式即可求出晶面间距离，，即厚度thickness

PS：之前111晶面是啥都给整忘了，算个点到直线的距离也能算错，真该抽自己两巴掌回去恶补高等数学😭

以下是来自维基百科的关于米勒指数以及111晶面的介绍，很形象了

平面ABC为：x + y + z − 1 = 0，点D为(1,1,0)

算一下厚度： $$ d=\dfrac{|Ax_0+By_0+Cz_0|}{\sqrt{A^2+B^2+C^2}}=\dfrac{1}{\sqrt{3}}a\approx \dfrac{1}{\sqrt{3}}\times 3.63468 \approx 2.0984 $$ Cleave Surface in Materials Studio

在切割晶面时，thickness设置为4.0 (4.0 = 1.0 × 4)，即四倍的晶胞参数、四层原子，按一下TAB，软件软件会自动计算四层原子的真是厚度（单位：Angstrom埃），大约是8.394 ≈ 4 × 2.0984，这与我们手动计算的晶面间距离2.0984的四倍相符

Step2：添加真空层，Build-Crystal-Build vacuum slab-Vacuum orientation：C

Step3：右键-lattice parameters

导出为.cif，vesta打开转换为POSCAR,

INCAR如下

System = Cu-conventional-slab
ISMEAR = 1
SIGMA = 0.1
ALGO = Fast
ENCUT = 450 # the system is too huge to adopt huge ENCUT
EDIFF = 1E-5
LDIPOL = .TRUE.
IDIPOL = 3
NWRITE = 0
LWAVE = .FALSE.
LCHARG = .FALSE.
NCORE = 20

Slab模型虽然是代表表面，但是实际上在z方向是固体-真空-固体-真空-…的交替。如果我们建立的Slab模型在z方向是非对称的，模型就会产生一个沿z方向的偶极。偶极会产生静电势，静电势接着会影响模型的镜像（周期性边界条件）。最后算出来的模型的总能量和力与真实情况是不相符的。因此我们需要方法去矫正这种虚假的静电影响。

一种常用的方法就是偶极矫正，在真空部分加入一个超窄的但是方向相反的偶极。这样一来，固体模型产生的偶极和真空中的偶极就会相互抵消。模型和其镜像之间的静电势影响就会抵消

微调KPOINTS

K-POINTS
 0
Gamma
  13 13 1
  0 0 0

计算结果

1 2	[ctan@baifq-hpc141 Cu-conventional-slab]$ grep ' without' OUTCAR energy without entropy= -55.91176854 energy(sigma->0) = -55.91397643

PS：最近学迷糊了，关于收敛脑海中有一个极其不正且的理解

[ctan@baifq-hpc141 bulk-calc]$ grep accuracy OUTCAR
 reached required accuracy - stopping structural energy minimisation
[ctan@baifq-hpc141 bulk-calc]$ cat OSZICAR
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1     0.103476848532E+03    0.10348E+03   -0.27847E+04  6720   0.246E+03
DAV:   2    -0.180387697818E+02   -0.12152E+03   -0.12091E+03  7280   0.155E+02
DAV:   3    -0.185227503579E+02   -0.48398E+00   -0.48396E+00 10160   0.132E+01
DAV:   4    -0.185240376489E+02   -0.12873E-02   -0.12873E-02  8200   0.464E-01
DAV:   5    -0.185240414620E+02   -0.38131E-05   -0.38134E-05 10480   0.143E-02    0.154E+01
DAV:   6    -0.164394189736E+02    0.20846E+01   -0.72731E+01  6800   0.881E+01    0.783E+00
DAV:   7    -0.149577606472E+02    0.14817E+01   -0.11165E+01  6720   0.344E+01    0.111E+00
DAV:   8    -0.149148685722E+02    0.42892E-01   -0.18431E-02  7280   0.139E+00    0.620E-01
DAV:   9    -0.149081819010E+02    0.66867E-02   -0.51836E-02  8440   0.220E+00    0.636E-02
DAV:  10    -0.149081753161E+02    0.65849E-05   -0.22202E-04  9000   0.149E-01    0.115E-02
DAV:  11    -0.149081715196E+02    0.37964E-05   -0.16862E-06  8520   0.968E-03    0.844E-04
DAV:  12    -0.149081716052E+02   -0.85577E-07   -0.88627E-08  4640   0.248E-03
   1 F= -.14908172E+02 E0= -.14905377E+02  d E =-.149082E+02
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.149363274123E+02   -0.28156E-01   -0.51946E-01  6720   0.830E+00    0.131E+00
DAV:   2    -0.149199449776E+02    0.16382E-01   -0.41291E-01  6720   0.646E+00    0.538E-01
DAV:   3    -0.149119014518E+02    0.80435E-02   -0.68998E-02  6720   0.271E+00    0.192E-01
DAV:   4    -0.149112107699E+02    0.69068E-03   -0.11896E-03  9800   0.216E-01    0.117E-02
DAV:   5    -0.149112123403E+02   -0.15704E-05   -0.80430E-06  7320   0.267E-02    0.594E-03
DAV:   6    -0.149112142891E+02   -0.19488E-05   -0.13229E-06  4760   0.997E-03    0.128E-04
DAV:   7    -0.149112142406E+02    0.48522E-07   -0.22610E-08  3560   0.101E-03
   2 F= -.14911214E+02 E0= -.14908447E+02  d E =-.304264E-02
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1    -0.150389301673E+02   -0.12772E+00   -0.26022E+00  6720   0.187E+01    0.290E+00
DAV:   2    -0.149566578541E+02    0.82272E-01   -0.20265E+00  6720   0.143E+01    0.119E+00
DAV:   3    -0.149175257623E+02    0.39132E-01   -0.34317E-01  6720   0.604E+00    0.437E-01
DAV:   4    -0.149140056079E+02    0.35202E-02   -0.72028E-03 10080   0.565E-01    0.378E-02
DAV:   5    -0.149139977907E+02    0.78172E-05   -0.11528E-04  7840   0.108E-01    0.134E-02
DAV:   6    -0.149140005388E+02   -0.27480E-05   -0.66656E-06  7680   0.220E-02    0.371E-04
DAV:   7    -0.149140005655E+02   -0.26792E-07   -0.64306E-08  3960   0.196E-03
   3 F= -.14914001E+02 E0= -.14911298E+02  d E =-.582896E-02

这是上一节里的关于bulk计算（结构优化）的部分内容，可以看到，reached required accuracy - stopping structural energy minimisation——离子步已经收敛，在之后的OSZICAR里，电子步迭代最多7步，并未超过设置的最大电子步数目。因此对于bulk计算，==电子步和离子步都收敛==

而对于slab模型的单点计算，我也傻傻的这么做😢

[ctan@baifq-hpc141 conventional-slab-static]$ grep accuracy OUTCAR
[ctan@baifq-hpc141 conventional-slab-static]$ cat OSZICAR
       N       E                     dE             d eps       ncg     rms          rms(c)
DAV:   1     0.178654188479E+04    0.17865E+04   -0.88940E+04 10388   0.264E+03
DAV:   2     0.200572973408E+03   -0.15860E+04   -0.15355E+04 10388   0.519E+02
DAV:   3    -0.535557174826E+02   -0.25413E+03   -0.23428E+03 10438   0.249E+02
DAV:   4    -0.683234131101E+02   -0.14768E+02   -0.14274E+02 11572   0.815E+01
DAV:   5    -0.687193157556E+02   -0.39590E+00   -0.39405E+00 12400   0.117E+01    0.292E+01
RMM:   6    -0.611494529970E+02    0.75699E+01   -0.15131E+02 11054   0.104E+02    0.161E+01
RMM:   7    -0.558685188476E+02    0.52809E+01   -0.32335E+01 10772   0.530E+01    0.226E+00
RMM:   8    -0.559296519539E+02   -0.61133E-01   -0.22789E+00 11552   0.614E+00    0.111E+00
RMM:   9    -0.559488617332E+02   -0.19210E-01   -0.37622E-01 11539   0.361E+00    0.509E-01
RMM:  10    -0.559454147061E+02    0.34470E-02   -0.13294E-01 12844   0.193E+00    0.214E-01
RMM:  11    -0.559437301097E+02    0.16846E-02   -0.22655E-02 12875   0.691E-01    0.208E-01
RMM:  12    -0.559413428432E+02    0.23873E-02   -0.55145E-03 11960   0.354E-01    0.586E-02
RMM:  13    -0.559414638149E+02   -0.12097E-03   -0.17465E-03 13023   0.187E-01    0.323E-02
RMM:  14    -0.559411451635E+02    0.31865E-03   -0.83089E-04 11756   0.134E-01    0.373E-02
RMM:  15    -0.559412390394E+02   -0.93876E-04   -0.47990E-04 11372   0.168E-01    0.284E-02
RMM:  16    -0.559414153722E+02   -0.17633E-03   -0.21175E-04 11475   0.107E-01    0.129E-02
RMM:  17    -0.559414836424E+02   -0.68270E-04   -0.10022E-04 11751   0.526E-02    0.889E-03
RMM:  18    -0.559415856898E+02   -0.10205E-03   -0.51724E-05 11898   0.366E-02    0.682E-03
RMM:  19    -0.559416259149E+02   -0.40225E-04   -0.26297E-05 11504   0.350E-02    0.861E-03
RMM:  20    -0.559416406670E+02   -0.14752E-04   -0.10748E-05 11263   0.208E-02    0.367E-03
RMM:  21    -0.559416498741E+02   -0.92071E-05   -0.47420E-06  8911   0.102E-02    0.119E-03
RMM:  22    -0.559416543655E+02   -0.44914E-05   -0.22030E-06  7155   0.977E-03    0.193E-03
RMM:  23    -0.559416573370E+02   -0.29715E-05   -0.14726E-06  6702   0.894E-03    0.135E-03
RMM:  24    -0.559416575122E+02   -0.17515E-06   -0.44209E-07  6318   0.507E-03
   1 F= -.55941658E+02 E0= -.55944345E+02  d E =0.806296E-02

不管怎么修改ALGO、ENCUT等参数，grep accuracy OUTCAR始终不返回reached required accuracy - stopping structural energy minimisation，后来我才意识到，这是==单点计算==，NSW = 0，并不会进行离子步迭代呀，所以我们只需关注电子步收敛即可，很明显电子步收敛，OUTCAR里的信息也可说明

1 2	[ctan@baifq-hpc141 conventional-slab-static]$ grep abort OUTCAR ------------------------ aborting loop because EDIFF is reached ----------------------------------------

This is definitely a stupid, stupid and stupid mistake😭

就写到这，吃饭去了，待会还要回实验室

使用Primitive Cell

使用bulk计算且收敛的CONTCAR进行slab模型的建立，操作还是类似的

值得注意的是，需要先==find symmetry==，然后才能转换为primitive cell

然后impose symmetry，接着再转换为primitive cell，然后建立slab模型

POSCAR会有些不同

CONTCAR\(1\1\1)
1.0
        2.5701000690         0.0000000000         0.0000000000
       -1.2850500345         2.2257719501         0.0000000000
        0.0000000000         0.0000000000        21.2954006195
   Cu
    4
Direct
     0.000000000         0.000000000         0.000000000
     0.666670000         0.333330000         0.098540000
     0.333330000         0.666670000         0.197080000
     0.000000000         0.000000000         0.295620000

计算后的能量也有所不同：

1 2	(base) storm@DESKTOP-HE4FQ8Q:~/my-learn/ex42/Cu-primitive-slab$ grep ' without' OUTCAR energy without entropy= -13.97087892 energy(sigma->0) = -13.97120814

值得一提的是，在白老师的计算集群上运行时报错了(vasp.6.4.2)

如果截断能ENCUT过小可能报错Error EDDDAV: Call to ZHEGV failed. Returncode = 7 1 8

[ctan@baifq-hpc141 Cu-primmitive-slab]$ ls
INCAR  KPOINTS  POSCAR  POTCAR
[ctan@baifq-hpc141 Cu-primmitive-slab]$ clear
[ctan@baifq-hpc141 Cu-primmitive-slab]$ sbatch vasp
Submitted batch job 4321
[ctan@baifq-hpc141 Cu-primmitive-slab]$ tail -f OUTCAR
--------------------------------------- Ionic step        1  -------------------------------------------




--------------------------------------- Iteration      1(   1)  ---------------------------------------


    POTLOK:  cpu time      0.0215: real time      0.0258
    SETDIJ:  cpu time      0.0109: real time      0.0118
     EDDAV:  cpu time      0.7451: real time      0.7613
       DOS:  cpu time      0.0052: real time      0.0060
    --------------------------------------------
      LOOP:  cpu time      0.7826: real time      0.8050

 eigenvalue-minimisations  :  2548
 total energy-change (2. order) : 0.3445715E+03  (-0.1837296E+04)
 number of electron      44.0000000 magnetization
 augmentation part       44.0000000 magnetization

 DIPCOR: dipole corrections for dipol
 direction  3 min pos   155,
 dipolmoment           0.000000      0.000000      0.000026 electrons x Angstroem
 Tr[quadrupol]       -31.061780

 energy correction for charged system         0.000000 eV
 dipol+quadrupol energy correction           -0.000000 eV
 added-field ion interaction          0.001451 eV  (added to PSCEN)


 Free energy of the ion-electron system (eV)
  ---------------------------------------------------
  alpha Z        PSCENC =        96.14159556
  Ewald energy   TEWEN  =     16384.30326894
  -Hartree energ DENC   =    -21746.47963692
  -exchange      EXHF   =         0.00000000
  -V(xc)+E(xc)   XCENC  =       182.01436818
  PAW double counting   =      4912.82663951    -5336.19386850
  entropy T*S    EENTRO =         0.00034874
  eigenvalues    EBANDS =       288.07434016
  atomic energy  EATOM  =      5563.88447798
  Solvation  Ediel_sol  =         0.00000000
  ---------------------------------------------------
  free energy    TOTEN  =       344.57153365 eV

  energy without entropy =      344.57118491  energy(sigma->0) =      344.57141740


--------------------------------------------------------------------------------------------------------




--------------------------------------- Iteration      1(   2)  ---------------------------------------


 -----------------------------------------------------------------------------
|                                                                             |
|     EEEEEEE  RRRRRR   RRRRRR   OOOOOOO  RRRRRR      ###     ###     ###     |
|     E        R     R  R     R  O     O  R     R     ###     ###     ###     |
|     E        R     R  R     R  O     O  R     R     ###     ###     ###     |
|     EEEEE    RRRRRR   RRRRRR   O     O  RRRRRR       #       #       #      |
|     E        R   R    R   R    O     O  R   R                               |
|     E        R    R   R    R   O     O  R    R      ###     ###     ###     |
|     EEEEEEE  R     R  R     R  OOOOOOO  R     R     ###     ###     ###     |
|                                                                             |
|     Error EDDDAV: Call to ZHEGV failed. Returncode = 7 1 8                  |
|                                                                             |
|       ---->  I REFUSE TO CONTINUE WITH THIS SICK JOB ... BYE!!! <----       |
|                                                                             |
 -----------------------------------------------------------------------------

小结

此计算是在bulk计算的基础上(最优结构)进行的单点计算，计算时不能进行离子弛豫，不能进行离子步迭代，默认NSW=0
不进行离子弛豫，那么自然不能设置离子步收敛条件，不能设置EDIFFG
如果不报错运行但结果不收敛可试着增大电子步数NSW以及截断能ENCUT

Comment and share

NEWER POSTS
OLDER POSTS
page 3 of 8

歌词

歌词深度分析

1. 主题与中心思想

2. 情绪氛围与感受

3. 语言风格与用词特色

4. 意象与象征意义

5. 修辞手法分析

6. 节奏与结构

7. 文化/社会背景

8. 多重解读角度

歌词

歌词深度分析

1. 主题与中心思想

2. 情绪氛围与感受

3. 语言风格与用词特色

4. 意象与象征意义

5. 修辞手法分析

6. 节奏与结构

7. 文化/社会背景

8. 多重解读角度

一些恒温器

关于表面能

关于表面积A的计算

p(1x1)Cu(111)的表面能

c(1x1)Cu(111)的表面能

p(1x1)Ni(100)的表面能

小结

工作站离线状态下numpy on python2地安装

优化前后对比

POSCAR的修改

INCAR的修改

IBRION = 1

IBRION =2

使用Conventional Cell建立slab模型

切割结构

计算结果

使用Primitive Cell

小结

John Doe

`POSCAR`的修改

`INCAR`的修改

`IBRION = 1`

`IBRION =2`