Lammps

♪ 在清华超算上提交计算

使用360急速浏览器或者ie(管理员模式)进入清华超算https://166.111.143.19:10000/svpn/vpnuser/home.cgi

清华超算采用Xhsell 6版本

进入文件所在的文件夹

运行命令

bsub -a intelmpi -e error.txt -o output.txt -J xx -i xx.in -n nn mpirun.lsf /home/anmeng/WORK1/lammps-3Mar20/src/lmp_mpi
    - xx: your commit id
    - xx.in: your file.in
    - nn: the number of cpu apply to compute
    - ‘/home/anmeng/WORK1/lammps-3Mar20/src/lmp_mpi’: the address of the exe

04.13

更新一个新编译的Lammps路径，同样的版本。之前编译deepmd，结果g++版本不够，但是已经回不去了，无法make clean-all，只好重新编译

bsub -a intelmpi -e error.txt -o output.txt -J xx -i xx.in -n nn mpirun.lsf /work1/anmeng_work/lammps-pure3Mar20/lammps-3Mar20/src/lmp_mpi
    - xx: your commit id
    - xx.in: your file.in
    - nn: the number of cpu apply to compute
    - ‘/home/anmeng/WORK1/lammps-3Mar20/src/lmp_mpi’: the address of the exe

♪ 计算声子谱

感谢泛柏舟、学力学的图图和樊哲勇老师对于计算声子谱所做出的贡献及讲解。

♪ ### VACF(velocity auto-correlation function)

Each atom’s contribution to the VACF is its current velocity vector dotted into its initial velocity vector at the time the compute was specified.

DOS

lammps中也有“compute vacf”的命令，往往也可以得到速度的自关联函数，原则上在这个基础上通过傅里叶变换，可以得到体系中的功率频谱，也就是分子动力学方法得到的声子谱，但是仔细阅读这一命令会发现在这里计算的VACF只有一个时间原点，没有办法如同在GK平衡态方法中计算热导率那般采取许多个时间原点进行平均，一个典型的输出例子如下：

可以发现这个声子谱显得非常凌乱

Manual Note 中强调，如果要重复计算，必须保持所计算的原子ID一样。我理解的也就是他们也明白这个不准，需要取平均，而且取平均的规则很严格。

在樊老师的博客里面给出了一个非常好的具体的代码来实现取多个时间原点平均的声子谱计算方法，只需要通过lammps输出一段时间的位移然后进行计算

在lammps中输入命令：

dump 1 all custom 1 ${simname}_voutput.lammpstrj id type vx vy vz

后处理方式在matlab代码已经讲的十分详细了，根据自己的需要修改即可。

一个处理好的示例文件可以在我的GitHub中找到

♪ 将扁盒子改成方盒子

Redefine lattice

$$ \left[ \begin{matrix} 1 & 0 & 0 \ 1 & 2 & 0 \ 0 & 0 & 1 \ \end{matrix} \right] $$

上面的记着好像不太对，下面的可能是正确的

$$ \left[ \begin{matrix} 1 & 1 & 0 \ -1 & 1 & 0 \ 0 & 0 & 1 \ \end{matrix} \right] $$

♪ 一个课题组干了一件很牛的事（给单体算聚合物的性质）

来自安博去学习的东京的课题组：

简单来说，给一个单体，不用再去淬火，再去驰豫，可以直接得出由单体构成聚合物的

Thermal conductivity 热导率
Thermal diffusivity 热扩散系数
Density 密度
Cp Cp
Cv Cv
Linear expansion coefficient 线性膨胀系数
Volumetric expansion coefficient 体积膨胀系数
Compressibility 压缩系数
Bulk modulus 体积弹性模量
Isentropic compressibility 等熵压缩
Isentropic bulk modulus 等熵体积弹性模量
Static dielectric constant 静态介电常数
Refractive index 折射率
Radius of gyration 回转半径
End-to-end distance 端到端的距离
Nematic order parameter 向列有序参数

没错，你只需要有一个Python 和 Lammps,还有对应的包，就可以把这些打包带走

心动不如行动，快来学习吧↓

文章链接

GitHub链接

♪ MS里建一个球体

首先，你需要一个盒子（尺寸可以小于半径，保证是单体就OK）

Build -> Build nanostructure -> Nanocluster -> Sphere 输入半径就OK

♪ 复杂混合体系中，存在大量单双建的建模问题

存在的问题：

① MS直接使用AC模块，部分会错误成键。且，AC模块无法控制核数，导致过程太慢！

② packmol，只可以识别一级键，无法识别二级建。

③ packmol建模后，使用MS的Calculated bond重新识别键，会存在错误成键，乱成键的问题解决办法：

依然使用packmol建模，只不过，导入MS后，不要全部计算键，只更改键的类型

-> Calculate Bond type

解决了键的类型，这样一来，建模再无大问题

♪ 将LAMMPS跑出的结构，重新导入MS中建模

不讲这一问题的重要意义了，直接上干货：

dump xx all custom Nfreq dump.lammpstrj id mol type element x y z ix iy iz
dump_modify element A B C D

已知需要输出的步数，输出这一步数的结构信息（①修改type 的123为具体元素，或②dump_modify 加入元素识别信息）

删除Timestep 0步的结构，之后导入VMD

这样，就可以保存好pdb结构，导入MS后，重新添加原子建模

♪ 切111切面并改盒子

选择build > surface > cleave surface

Cleave plane 1 1 1 Then Cleave

切后的结构是没有"封顶"的，需要手动设置一个真空层

build > crystal > build vacuum slab

vacuum thickness选择0 Å, 就可以得到最小的顶部盖子

可以看到当前坐标系与晶格并不平行

Lattice Parameters - Advanced - Re-oriented to standard

转为正交晶格，然后再参考之前 Redefine lattice,变为方盒子，然后重复上方转正操作即可。

♪ TIP4P水模型

这是TIP4P水模型的始祖文章，点这里

从上图可以看到，体系中水分子即使存在键角及相关能量参数，但对于体系的能量都是没有贡献的。图中蓝色为加入了除了水分子之外的其他的键

这是由于，TIP4P/TIP3P水模型中，键角能量是不统计进总能量的。从统计水分子自身的能量式可以看到

其中涉及了电荷以及带电原子(氧原子)间的距离。并不涉及键角参数。

也可以看到，并没有设计到原子的能量参数，那么能量参数用在哪呢？当然是用在静电作用了。所以也可以得出另一个结果，如果只统计水分子间作用的话，哪怕不给氧原子的能量参数也是可以的。

Wikipedia also has a nice article on water models

♪ 为什么需要Voronoi体积？

应力，是固体力学乃至连续介质力学无法绕开的议题，整个固体力学在讨论本构关系的时候，其中一极就是应力。在lammps中是不存在这样的应力概念的，因为lammps中的界面是离散的。但是在计算材料科学中，我们往往会用lammps计算一种新的材料的模量或者强度，那么如何将lammps中的“应力”转化为宏观固体力学意义上的应力呢。

阅读lammps mannual中关于compute stress/atom command中提到，关于lammps中应力的量纲，提到

意思是说，lammps的应力单位是“体积×压强”，这是一个能量项，但绝不是eV

以Metal单位举例：距离的单位为Å，那么体积单位为$Å^3$，压强的单位为bars，也就是100KPa，

把这些转化为国际单位就是

$$bars(A^3)=10^5(N/m^2)(10^{-10})^3(m^3) = 10^{-25}(N*m)$$，

这就是lammps中应力的单位，而并非是

$$ 1ev=1.610^{-19}J=1.610^{-19}N*m $$

既然lammps应力的单位表现出能量的含义，那么转化为宏观意义上的应力则需要除以一个体积，那么这个体积该如何得到呢?

这是没有定论的！

lammps手册也提到（红色标记）：这很难定义出来。

我认为，有以下三种计算体积的办法：

1 采用compute voronoi/atom命令（下文会将如何安装Voronoi包）

2 自己计算每个原子的体积，乘以总数。然后统计组内原子的应力，除以原子个数和单个体积。

compute         perCstress cnt stress/atom NULL
compute         szzx cntstresszz reduce sum c_perCstress[1]
variable        szzxbar equal c_szzx/(8.784*count(cntstresszz))*(10^-4)

其中8.784是单个碳原子的体积 (1)，单位为$Å^3$,这个脚本计算了一个group cntstresszz在x方向的平均应力

3 还有一种方法文献(2)，就是采用连续假设，统计出所有应力和，然后除以整个结构的应力张量的截面的面积。比如把纳米线抽象成一个圆柱，求圆柱的体积。这样的连续假设也常见于计算石墨烯或碳纳米管的应力，往往假设石墨烯或碳纳米管厚度为0.35nm，正好是石墨烯的堆叠平衡距离。

(1)Lee, J. G. (2016). Computational materials science: an introduction, Crc Press.

(2)Roman, R. E., et al. (2015). "Mechanical properties and defect sensitivity of diamond nanothreads." Nano Lett 15(3): 1585-1590.

♪ 安装Voronoi包

1.首先下载voro++包

2.在lammps/lib/voronoi的文件夹下，解压voro++-0.4.6

命令：tar -xvf voro++-0.4.6

3.解压结束后，打开voro++-0.4.6，输入命令:make

4.make结束后，直接输入命令：sudo make install

5.创造链接：输入命令：ln -s voro++-0.4.6/src includelink

                  `ln -s voro++-0.4.6/src liblink`

6.编译lib/voronoi文件下的Makefile.lammps,将Makefile.lammps改成如下格式

voronoi_SYSINC = -I/usr/local/include/voro++
voronoi_SYSLIB = -lvoro++
voronoi_SYSPATH = -L/usr/local/lib

7.进入到lammps/src文件夹下，

make yes-voronoi

make mpi

**超算由于没有权限，无法sudo。 what a pity

♪ 计算透射系数

声子态密度表征了声子的本征特征，表明了可以在哪些频率进行传输，而传输能力则取决于界面左右两侧声子的透射系数。透射系数越强，传输能力越强。

感谢梁挺博士的贡献，根据他的工作，我在我的Github上给出了一个算例。

引用本代码的话，以下是一些相关工作：

(1) K. Sääskilahti, J. Oksanen, J. Tulkki, and S. Volz, Phys. Rev. B 90, 134312 (2014)

(2) K. Sääskilahti, J. Oksanen, S. Volz, and J. Tulkki, Phys. Rev. B 91, 115426 (2015)

(3) Xu K, Deng S, Liang T, et al, Efficient mechanical modulation of the phonon thermal conductivity of Mo6S6 nanowires. Nanoscale, 2022

(4) An M, Chen D, Ma W, et al. Directly visualizing the crossover from incoherent to coherent phonons in two-dimensional periodic MoS2/MoSe2 arrayed heterostructure. International Journal of Heat and Mass Transfer, 2021

1 首先，编译Scriptscompactify_vels.cpp文件，并将生成的compactify_vels添加到命令路径中。

g++  compactify_vels.cpp  -o  compactify_vels

2 跑lammps文件l2%_relax_thermal.in，以期得到界面左右两侧原子的速度文件

本例研究了压缩的BC3，首先进行了压缩

3 转换速度文件vel.dat为透射系数计算所需的文件

compactify_vels vels.dat

4 由于NPT和压缩过程，界面左右两边原子会有微小变化，单纯采用距离统计界面左右两端原子，会造成原子数目少量增减情况 因此，本例对forces.in做了改进。

这样一来，直接输出界面左右两边原子ID，并替换距离定义的group。在接下来的计算中就不会存在报错情况。

5 运行SHC_generate.py

可以看到，生成了相关图片。为了方便做图，本例将生成的图片分别保存为.txt文件，

光谱热导率(frequencies_and_ITC.txt)

透射系数(frequencies_T_W.txt)

透射系数积分(frequencies_accumulated_ITC.txt)

♪ 计算周期性条件下的氢键

为什么不使用VMD进行方便快捷的计算？是因为他无法考虑周期性。

以下是各种程序计算氢键的准则，详见这里。

在此，我们采用MDAnalysis方法。

首先，安装他

pip install --upgrade MDAnalysis

撰写python文件。Anyway，你需要一个开头。

import MDAnalysis as mda
import numpy as np
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis
from matplotlib import pyplot as plt


u = mda.Universe('finish.data','NPT_D.trr')

print (u)
print (u.trajectory)

这里包含了画图以及后处理的import。需要注意的是u的定义，finish.data是初始帧的文件，NPT_D.trr是轨迹文件。

这里修改了元素的lammpstrij轨迹文件无法被读入。因此需要将其导入vmd进行转换。

经过测试，可以识别的文件类型包括：mol2, xyz, dcd, trr and gro。

But！mol2和xyz的文件计算结果是2.815个，gro的计算结果是3.12（但是这里只包含了一帧，难以计算其寿命），dcd和trr计算结果是3.13个。

因此，采用trr或dcd文件均可。

这里给出一个完整的例子

import MDAnalysis as mda
import numpy as np
from MDAnalysis.analysis.hydrogenbonds import HydrogenBondAnalysis
from matplotlib import pyplot as plt


u = mda.Universe('finish.data','NPT_D.trr')

print (u)
print (u.trajectory)

O_atom = u.select_atoms('type 1')


hbonds = HydrogenBondAnalysis(
    universe=u,
    donors_sel='type 1 2',
    hydrogens_sel='type 2',
    acceptors_sel='type 1',
    d_a_cutoff=3.5,
    d_h_a_angle_cutoff=138.1,
    update_selections=False
)
hbonds.run()

plt.plot(hbonds.times, hbonds.count_by_time(), lw=2)

plt.title("Number of hydrogon bonds over time", weight="bold")
plt.xlabel("Time (ns)")
plt.ylabel(r"$N_{HB}$")

plt.show()

print(hbonds.count_by_time())

with open('output.txt', 'w') as f:

    print(hbonds.count_by_time(),file=f)

其中：type 1或name O均可，这取决于定型的finish.data中如何定义。

donors_sel 表示供体是谁；hydrogens_sel 表示受体氢的选择；acceptors_sel 表示受体氧的选择。

这是最终的结果，表现了平均氢键数目随时间的变化关系。可以看到NPT阶段末期氢键数目已经趋于稳定。

注意：通过核对print后u的值可以检验总原子个数以及多少帧（包括1+N，1为定型文件，N为轨迹数目）。这里总原子数是19200包括了水和其中的其他原子，表示只要选对了类型，无需只dump水的轨迹。