空间站广场

论文

Notebooks

比赛

课程

Apps

我的主页

我的Notebooks

我的论文库

我的足迹

我的工作空间

任务

节点

文件

数据集

镜像

项目

数据库

公开

樊哲勇《分子动力学模拟》python实例 | 3.Tersoff 势的编程实现

樊哲勇《分子动力学模拟》程序代码python实现

AI4S101

樊哲勇《分子动力学模拟》程序代码python实现AI4S101

yuxiangc22

mosey

更新于 2024-06-26

推荐镜像 :Basic Image:bohrium-notebook:2023-03-26

推荐机型 :c2_m4_cpu

数据集

on实现(v5)

樊哲勇《分子动力学模拟》python实例 | 3.Tersoff 势的编程实现

樊哲勇老师最近出版的书籍《分子动力学模拟》中，使用C++语言进行分子模拟程序的实现。未了更好地理解分子动力学模拟，AI4S1O1-材料小队学习志愿者决定针对樊老师书籍中的内容，进行python语言的程序实现以及注释说明。

本节内容对应樊老师书中的第四章

樊哲勇《分子动力学模拟》python实例 | 概述

代码

文本

©️ Copyright 2023 @ Authors
作者：ChenYuxiang LiDenan(AI4S1O1-材料小队学习志愿者) 📨
日期：2024-6-10
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
快速开始：点击上方的开始连接按钮，选择 bohrium-notebook:2023-03-26 镜像及 c2_m4_cpu 节点配置即可开始运行。

代码

文本

Tersoff 势

本节教程目录如下：
1. Tersoff 势介绍
2. Tersoff 势的编程实现

代码

文本

1.Tersoff 势介绍

Tersoff势是一种广泛应用于材料科学和分子动力学模拟的原子间相互作用势函数，特别适用于描述半导体材料（如碳、硅和锗）中的键合和结构性质。Tersoff势是一种经验势，专门设计用来模拟共价键合系统，能够有效捕捉这些材料的物理和化学性质。

主要特点

局部环境依赖： Tersoff势的一个关键特点是它不仅依赖于原子间的距离，还依赖于原子的局部环境。这使得它能够准确描述不同键合状态（如键角和多体相互作用）的能量变化。

多体相互作用：与简单的两体势不同，Tersoff势考虑了多体相互作用。具体而言，势函数不仅考虑到两个原子之间的直接相互作用，还考虑到周围原子的影响，从而更准确地描述键角和其他几何结构特征。

键次效应： Tersoff势能够模拟键次效应，即一个原子的键合情况会影响其周围键的强度和性质。这对于描述复杂的材料结构和相变行为非常重要。

Tersoff势的参数化依赖于所研究的具体材料系统。参数通常通过拟合实验数据或更高精度的计算（如密度泛函理论计算）得到。这些参数确保势函数能够准确预测材料的物理性质，如晶格常数、弹性模量、表面能等。

代码

文本

Tersoff势函数公式

$V_{ij} = f_{C} (r_{ij}) [f_{R} (r_{ij}) + b_{ij} f_{A} (r_{ij})]$ 其中：

$f_{C} (r_{ij})$ 是截断函数，确保相互作用在一定距离外迅速减小到零。 $f_{R} (r_{ij})$ 是排斥项，描述原子间的排斥力。 $f_{A} (r_{ij})$ 是吸引项，描述原子间的吸引力。 $b_{ij}$ 是键次项，依赖于原子周围的局部环境，描述多体相互作用的修正。

$f_{C} (r_{ij}) = ⎩ ⎨ ⎧ 1 \frac{1}{2} + \frac{1}{2} cos (\frac{π ( r _{ij} - R + D )}{2 D}) 0 r_{ij} < R - D R - D \leq r_{ij} \leq R + D r_{ij} > R + D$

$f_{R} (r_{ij}) = A e^{- λ_{1} r_{ij}}$

$f_{A} (r_{ij}) = - B e^{- λ_{2} r_{ij}}$ 其中，参数 A 和 B 具有能量的量纲，参数 λ1 和 λ2 具有长度倒数的量纲。

Tersoff 势的多体性质包含在因子 bij 中。该因子称为键序（即化学键的某种定量描述），它可表达为如下形式： $b_{ij} = (1 + β^{n} ζ_{ij}^{n})^{- \frac{1}{2 n}}$

$ζ_{ij} = k \neq = i, j \sum f_{C} (r_{ik}) g_{ijk}$ $g_{ijk} = 1 + \frac{c ^{2}}{d ^{2}} - \frac{c ^{2}}{d ^{2} + ( h - cos θ _{ijk} ) ^{2}}$ 其中的参数 β、 n、 c、 d、 h 都是无量纲的。式中有一个角度变量 θijk ，它代表由键ij 和 ik 形成的键角。该键角的余弦为： $cos θ_{ijk} = \frac{r _{ij} \cdot r _{ik}}{r _{ij} r _{ik}}$ Tersoff 势的部分力表达式: $\frac{\partial U _{i}}{\partial r _{ij}} = \frac{1}{2} \frac{\partial f _{C} ( r _{ij} )}{\partial r _{ij}} [f_{R} (r_{ij}) - b_{ij} f_{A} (r_{ij})] \frac{\partial r _{ij}}{\partial r _{ij}} + \frac{1}{2} f_{C} (r_{ij}) [\frac{\partial f _{R} ( r _{ij} )}{\partial r _{ij}} - b_{ij} \frac{\partial f _{A} ( r _{ij} )}{\partial r _{ij}}] \frac{\partial r _{ij}}{\partial r _{ij}} - \frac{1}{2} \frac{\partial f _{C} ( r _{ij} )}{\partial r _{ij}} k \neq = i, j \sum f_{C} (r_{ik}) f_{A} (r_{ik}) \frac{\partial b _{ik}}{\partial ζ _{ik}} g_{ijk} \frac{\partial r _{ij}}{\partial r _{ij}} - \frac{1}{2} f_{C} (r_{ij}) f_{A} (r_{ij}) \frac{\partial b _{ij}}{\partial ζ _{ij}} k \neq = i, j \sum f_{C} (r_{ik}) \frac{\partial g _{ijk}}{\partial cos θ _{ijk}} \frac{\partial cos θ _{ijk}}{\partial r _{ij}} - \frac{1}{2} f_{C} (r_{ij}) k \neq = i, j \sum f_{C} (r_{ik}) f_{A} (r_{ik}) \frac{\partial b _{ik}}{\partial ζ _{ik}} \frac{\partial g _{ijk}}{\partial cos θ _{ijk}} \frac{\partial cos θ _{ijk}}{\partial r _{ij}} .$ 其中 $\frac{\partial cos θ _{ijk}}{\partial r _{ij}} = \frac{1}{r _{ij}} (\frac{r _{ik}}{r _{ik}} - \frac{r _{ij}}{r _{ij}} cos θ_{ijk}) .$ $\frac{\partial b _{ij}}{\partial ζ _{ij}} = \frac{b _{ij} β ^{n} ζ _{ij}^{n}}{2 ( 1 + β ^{n} ζ _{ij}^{n} ) ζ _{ij}}$ 通过对以上参数的求解，可以计算两个原子之间的Tersoff势作用力。

代码

文本

2. Tersoff 势的编程实现

Tersoff势计算主要是各种参数的求解，根据不同原子间的空间信息求解Tersoff势的键序等参数，遍历空间所有原子和其近邻原子得到总体的Tersoff势作用力。

进行Tersoff势作用力计算之后需要进一步的验证计算结果，根据理论体系的自洽性对程序进行检验是一种方法。这里，我们通过能量与力的自洽性进行检验，采用了最低阶的中心差分法。

代码

文本

[ ]

#Tersoff势各种公式参数求解

def find_fr_and_frp(d12):

"""

Calculate the values of fr and frp for a given distance.

Parameters:

d12 (float): Distance value.

"""

a = 1.8308e3

lambda_val = 2.4799

fr = a * math.exp(-lambda_val * d12)

frp = -lambda_val * fr

return fr, frp

def find_fa_and_fap(d12):

"""

Calculate the values of fa and fap for a given distance.

Parameters:

d12 (float): Distance value.

"""

b = 471.18

mu = 1.7322

fa = b * math.exp(-mu * d12)

fap = -mu * fa

return fa, fap

def find_fa(d12):

"""

Calculate the value of fa for a given distance.

Parameters:

d12 (float): Distance value.

"""

b = 471.18

mu = 1.7322

fa = b * math.exp(-mu * d12)

return fa

def find_fc_and_fcp(d12):

"""

Calculate the values of fc and fcp for a given distance.

Parameters:

d12 (float): Distance value.

"""

r1 = 2.7

r2 = 3.0

pi = math.pi

# Pre-compute a factor involving pi and the cutoff radii

pi_factor = pi / (r2 - r1)

# Determine the value of fc based on the distance d12

if d12 < r1:

fc = 1.0

fcp = 0.0

elif d12 < r2:

# Between the first and second cutoff radii, fc and fcp are computed using cosine function

fc = math.cos(pi_factor * (d12 - r1)) * 0.5 + 0.5

fcp = -math.sin(pi_factor * (d12 - r1)) * pi_factor * 0.5

else:

fc = 0.0# Beyond the second cutoff radius, fc is 0

fcp = 0.0 # The derivative fcp is 0

return fc, fcp

def find_fc(d12):

"""

Calculate the value of fc for a given distance.

Parameters:

d12 (float): Distance value.

"""

r1 = 2.7

r2 = 3.0

pi = math.pi

# Pre-compute a factor involving pi and the cutoff radii

pi_factor = pi / (r2 - r1)

# Determine the value of fc based on the distance d12

if d12 < r1:

fc = 1.0

elif d12 < r2:

# Between the first and second cutoff radii, fc and fcp are computed using cosine function

fc = math.cos(pi_factor * (d12 - r1)) * 0.5 + 0.5

else:

fc = 0.0 # Beyond the second cutoff radius, fc is 0

return fc

def find_g_and_gp(cos_value):

"""

Calculate the values of g and its derivative gp for a given cosine value.

Parameters:

cos_value (float): Cosine value.

"""

# Parameter for the Tersoff potential

c = 1.0039e5

d = 16.217

h = -0.59825

c2 = c * c

d2 = d * d

# Pre-compute a constant factor

c2overd2 = c2 / d2

# Calculate intermediate term

temp = d2 + (cos_value - h) * (cos_value - h)

# Calculate g and its derivative gp

g = 1.0 + c2overd2 - c2 / temp

gp = 2.0 * c2 * (cos_value - h) / (temp * temp)

return g, gp

def find_g(cos_value):

"""

Calculate the value of g for a given cosine value.

Parameters:

cos_value (float): Cosine value.

"""

# Parameter for the Tersoff potential

c = 1.0039e5

d = 16.217

h = -0.59825

c2 = c * c

d2 = d * d

c2overd2 = c2 / d2

# Calculate intermediate term

temp = d2 + (cos_value - h) * (cos_value - h)

# Calculate g

g = 1.0 + c2overd2 - c2 / temp

return g

代码

文本

[ ]

#计算Tersoff势力

def find_force(self):

"""

计算原子的势能和力。

参数:

atom (Atom): Atom对象，包含原子信息。

"""

self.find_b_and_bp()

self.find_force_tersoff()

def find_b_and_bp(self):

"""

Calculate the bond order term (b) and its derivative (bp).

"""

# Parameter for Tersoff potential

beta = 1.1000e-6

n = 0.78734

minus_half_over_n = -0.5 / n

# Loop over all atoms

for n1 in range(self.number):

# Loop over neighbors of atom n1

for i1 in range(self.NN[n1]):

n2 = self.NL[n1 * self.MN + i1]# Neighbor n2 of atom n1

x12 = self.x[n2] - self.x[n1]

y12 = self.y[n2] - self.y[n1]

z12 = self.z[n2] - self.z[n1]

x12, y12, z12 = self.apply_mic(x12, y12, z12)

d12 = math.sqrt(x12 * x12 + y12 * y12 + z12 * z12)

# Initialize the zeta parameter

zeta = 0.0

# Loop over other neighbors of atom n1

for i2 in range(self.NN[n1]):

n3 = self.NL[n1 * self.MN + i2]

if n3 == n2:

continue

# Calculate the relative position and distance between atom n1 and neighbor n3

x13 = self.x[n3] - self.x[n1]

y13 = self.y[n3] - self.y[n1]

z13 = self.z[n3] - self.z[n1]

x13, y13, z13 = self.apply_mic(x13, y13, z13)

d13 = math.sqrt(x13 * x13 + y13 * y13 + z13 * z13)

# Calculate the cosine of the angle between bonds n1-n2 and n1-n3

cos = (x12 * x13 + y12 * y13 + z12 * z13) / (d12 * d13)

# Calculate the cutoff function for d13

fc13 = find_fc(d13)

# Calculate the angle-dependent function

g123 = find_g(cos)

# Accumulate the zeta parameter

zeta += fc13 * g123

# Calculate the bond order term b12

bzn = math.pow(beta * zeta, n)

b12 = math.pow(1.0 + bzn, minus_half_over_n)

# Store the bond order term

self.b[n1 * self.MN + i1] = b12

# Calculate the derivative of the bond order term bp12

self.bp[n1 * self.MN + i1] = -b12 * bzn * 0.5 / ((1.0 + bzn) * zeta)

def find_force_tersoff(self):

"""

alculate the Tersoff potential and forces.

"""

# Initialize potential energy

self.pe = 0.0

# Initialize force

self.fx = [0.0] * self.number

self.fy = [0.0] * self.number

self.fz = [0.0] * self.number

# Loop over all atoms

for n1 in range(self.number):

neighbors = self.NL[n1 * self.MN : n1 * self.MN + self.NN[n1]]

for i1,n2 in enumerate(neighbors):

# Calculate the relative position and distance between atom n1 and neighbor n2

x12 = self.x[n2] - self.x[n1]

y12 = self.y[n2] - self.y[n1]

z12 = self.z[n2] - self.z[n1]

x12, y12, z12 = self.apply_mic(x12, y12, z12)

d12 = math.sqrt(x12 * x12 + y12 * y12 + z12 * z12)

d12inv = 1.0 / d12

# Calculate force field parameters

fc12, fcp12 = find_fc_and_fcp(d12)

fa12, fap12 = find_fa_and_fap(d12)

fr12, frp12 = find_fr_and_frp(d12)

b12 = self.b[n1 * self.MN + i1]

bp12 = self.bp[n1 * self.MN + i1]

# Initialize force vector

f12 = [0.0, 0.0, 0.0]

factor1 = -b12 * fa12 + fr12

factor2 = -b12 * fap12 + frp12

factor3 = (fcp12 * factor1 + fc12 * factor2) / d12

# Calculate the force components

f12[0] += x12 * factor3 * 0.5

f12[1] += y12 * factor3 * 0.5

f12[2] += z12 * factor3 * 0.5

# Loop over other neighbors of atom n1

for i2 in range(self.NN[n1]):

n3 = self.NL[n1 * self.MN + i2]

if n3 == n2:

continue

# Calculate the relative position and distance between atom n1 and neighbor n3

x13 = self.x[n3] - self.x[n1]

y13 = self.y[n3] - self.y[n1]

z13 = self.z[n3] - self.z[n1]

x13, y13, z13 = self.apply_mic(x13, y13, z13)

d13 = math.sqrt(x13 * x13 + y13 * y13 + z13 * z13)

# Calculate the three-body term

fc13 = find_fc(d13)

fa13 = find_fa(d13)

bp13 = self.bp[n1 * self.MN + i2]

cos123 = (x12 * x13 + y12 * y13 + z12 * z13) / (d12 * d13)

g123, gp123 = find_g_and_gp(cos123)

cos_x = x13 / (d12 * d13) - x12 * cos123 / (d12 * d12)

cos_y = y13 / (d12 * d13) - y12 * cos123 / (d12 * d12)

cos_z = z13 / (d12 * d13) - z12 * cos123 / (d12 * d12)

factor123a = (-bp12 * fc12 * fa12 * fc13 -

bp13 * fc13 * fa13 * fc12) * gp123

factor123b = -bp13 * fc13 * fa13 * fcp12 * g123 * d12inv

# Update force components

f12[0] += (x12 * factor123b + factor123a * cos_x) * 0.5

f12[1] += (y12 * factor123b + factor123a * cos_y) * 0.5

f12[2] += (z12 * factor123b + factor123a * cos_z) * 0.5

# Update potential energy and forces

self.pe += factor1 * fc12 * 0.5

self.fx[n1] += f12[0]

self.fy[n1] += f12[1]

self.fz[n1] += f12[2]

self.fx[n2] -= f12[0]

self.fy[n2] -= f12[1]

self.fz[n2] -= f12[2]

代码

文本

[ ]

#总控制函数，采用了最低阶的中心差分法，进行自洽性检验，

def main():

atom = Atom()

atom.find_neighbor_on2()

print("neighbor list is built.")

# Calculate analytical forces and write to file

with open("force.out", "w") as ofile:

atom.find_force()

for n in range(atom.number):

ofile.write(f"{atom.fx[n]:.16f} {atom.fy[n]:.16f} {atom.fz[n]:.16f}\n")

print("analytical force is calculated.")

# Calculate finite-difference forces and append to file

delta = 2.0e-5

with open("force.out", "a") as ofile:

for n in range(atom.number):

# Perturb x coordinate and calculate forces

atom.x[n] = atom.x0[n] + delta

atom.find_force()

pe_positive = atom.pe

atom.x[n] = atom.x0[n] - delta

atom.find_force()

pe_negative = atom.pe

atom.x[n] = atom.x0[n]

# Calculate finite-difference force derivative and write to file

ofile.write(f"{(pe_negative - pe_positive) / (delta * 2.0):.16f} ")

# Perturb y coordinate and calculate forces

atom.y[n] = atom.y0[n] + delta

atom.find_force()

pe_positive = atom.pe

atom.y[n] = atom.y0[n] - delta

atom.find_force()

pe_negative = atom.pe

atom.y[n] = atom.y0[n]

ofile.write(f"{(pe_negative - pe_positive) / (delta * 2.0):.16f} ")

# Perturb z coordinate and calculate forces

atom.z[n] = atom.z0[n] + delta

atom.find_force()

pe_positive = atom.pe

atom.z[n] = atom.z0[n] - delta

atom.find_force()

pe_negative = atom.pe

atom.z[n] = atom.z0[n]

ofile.write(f"{(pe_negative - pe_positive) / (delta * 2.0):.16f}\n")

print("finite-difference force is calculated.")

代码

文本

python总代码实现

第一版（python-list版本，简单）

第二版（numpy版本，简洁）

代码

文本

[1]

# 出于安全考虑，我们没有数据集所在文件夹的写入权限，因此我们将其复制到 `/personal/` 目录下,预计耗时 15s

! mkdir -p /personal/04/ && cp -nr /bohr/AI4S101-python-rs2o/v5/chapter-04/ /personal/04/

代码

文本

[2]

import os

# 切换路径

os.chdir("/personal/04/")

# 展示路径，文件，树状图

! pwd && ls && tree -L 2

/personal/04
chapter-04
.
└── chapter-04
    ├── frist_edition
    ├── second_edition
    └── tersoff.cpp

3 directories, 1 file

代码

文本

[4]

#查看run.in文件

! cd ./chapter-04/frist_edition/ && cat xyz.in

64
10.9 0.0 0.0 0.0 10.9 0.0 0.0 0.0 10.9
Si 0.10868992837464586 -0.032134702424903905 0.2187017222224601
Si -0.16377199226355477 2.8344268955061516 2.7686329237480223
Si 2.924303081565661 -0.06653451020698886 2.8233959400451725
Si 2.7289325506164364 2.9182767349309238 0.23076689830057823
Si 1.3738843150062376 1.1181285271430546 1.31887686428363
Si 1.2868552257796053 4.124342119490165 4.280519161225453
Si 4.213409782739129 1.425180316855076 4.27522285707026
Si 4.2429721295422205 4.008002267523735 1.5187917541087543
Si 0.07685711575733634 -0.012800324450596368 5.41637196544254
Si -0.21608477123736125 2.5090385975786944 7.999899164294128
Si 2.551903290764507 -0.014486795034510891 8.34004463810667
Si 2.9062705303823044 2.7079582434385676 5.326609309891285
Si 1.3305932675174954 1.2951164177647623 6.8371949632025615
Si 1.1388471183834088 4.176075735777663 9.297888026569774
Si 4.127246795498884 1.5738614329837757 9.32033674724465
Si 4.258179509720246 3.889443413733163 6.857007438325161
Si 0.00028480261417884734 5.39989154147427 0.07467423236715748
Si -0.23011497149448473 8.33893486821361 2.756758441814191
Si 2.8089428202296274 5.424982552756125 2.643884130462187
Si 2.693796840096173 7.974953816394539 -0.1007762329143288
Si 1.1542958572194892 6.915264128243375 1.47775719008434
Si 1.1316044291509757 9.542647964418103 3.8551757126103574
Si 3.985149774107869 6.696459621950303 3.970409507562836
Si 4.245554867574212 9.640483822821084 1.4426732028435962
Si -0.16323136850147296 5.52850990542994 5.503905435605338
Si 0.2488726990509471 8.41829039048415 8.340552382035298
Si 2.519421330465362 5.299886660374831 8.395733895294441
Si 2.73958912905055 8.304906283426151 5.5614974355016225
Si 1.26081135206559 6.927036262825779 6.612339496825972
Si 1.3385326583097048 9.744855681305083 9.354309467000677
Si 3.845964976019473 6.650675393471734 9.481840088005475
Si 3.9585356562836926 9.334458196398666 6.749167969527589
Si 5.471134122230594 -0.10546039350252456 0.21597026172130723
Si 5.241538513883275 2.7301414932748087 2.675063174189857
Si 8.185952893716127 0.12941557704810952 2.8874602250512824
Si 8.200760660351413 2.71776806398955 -0.007730313413334022
Si 6.783740185592028 1.4160278428546986 1.128340273710156
Si 6.663181007819496 4.198835728376649 4.145400222972699
Si 9.39231395361778 1.216901898033485 3.838629296843775
Si 9.397876436704726 4.314557597802333 1.3975670606176072
Si 5.384860228109241 0.07772491670947224 5.299844982238505
Si 5.41583578304503 2.8549363995740706 8.192050742410647
Si 7.9494183299810866 0.15658945491139714 8.309055273585361
Si 7.964226679620323 2.9571477546762517 5.447500146461646
Si 6.762292592994396 1.1343980301669259 7.020073142488122
Si 6.643399930613076 3.9800729554690957 9.300463279354044
Si 9.472566872684194 1.136509405415022 9.76912625686094
Si 9.449543851600112 3.9727335365420204 6.92676811445636
Si 5.6345352994554325 5.206204033662746 0.04300494869048027
Si 5.496991752033784 7.988102715150229 2.618668895313702
Si 8.195776944486434 5.295575843923728 2.7010305029835644
Si 8.232677701928264 8.138219661159287 0.03791326418649349
Si 6.737507448216503 6.565278813947499 1.372282451897921
Si 6.594637673494488 9.72158152613276 4.107727774258703
Si 9.351145232682786 6.650201452895643 4.081851392932531
Si 9.46167183715384 9.465721942938783 1.229205999136364
Si 5.56864177791515 5.37176308800035 5.2614019237384095
Si 5.663568711548033 8.213523828818152 7.996988113940924
Si 8.187708744236973 5.688781240691262 7.955952234812813
Si 8.23428798826512 8.298305192328604 5.397228005675858
Si 6.952267493476326 7.04956351935475 6.966104001181325
Si 6.6738433971675315 9.613399924237317 9.603138974494346
Si 9.775877284385647 6.983931621664527 9.340341760319745
Si 9.512267526432813 9.490936916245106 6.919173061208715

代码

文本

[5]

# 尝试在文件夹下运行数据集中的主程序代码

! cd ./chapter-04/frist_edition/ && python Tersoff.py

neighbor list is built.
analytical force is calculated.
finite-difference force is calculated.

代码

文本

[9]

! cd ./chapter-04/frist_edition/ && ls

Tersoff.py     fig-chapter-8-tersoff-validation.pdf  plot_force.py
create_xyz.py  force.out			     xyz.in

代码

文本

[11]

import matplotlib.pyplot as plt

import numpy as np

force = np.loadtxt('./chapter-04/frist_edition/force.out')

plt.figure()

plt.plot(force[64:,:]-force[0:64,:], 'o')

plt.xlabel('atom index')

plt.ylabel('Force difference (eV/Å)')

plt.savefig('fig-chapter-8-tersoff-validation.pdf')

plt.show()

代码

文本

第一版（python-list版本，简单）

代码

文本

[ ]

999

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

›

⌄

#第一版（python-list版本，简单）

import math

class Atom:

"""

This class represents an atom in a molecular simulation.

"""

def __init__(self,filename: str='xyz.in') -> None:

self.filename = filename

# Read data from the XYZ file

self.number, self.box, self.x, self.y, self.z, self.x0, self.y0, self.z0= self.read_xyz(self.filename)

# Initialize neighbor list parameters

self.cutoffNeighbor = 3.1

self.MN = 4

self.NN = [0] * self.number

self.NL = [0] * (self.number*self.MN)

# Initialize position, velocity, and force arrays

self.fx = [0.0] * self.number

self.fy = [0.0] * self.number

self.fz = [0.0] * self.number

self.pe = 0.0

self.b = [0.0] * (self.number*self.MN)

self.bp = [0.0] * (self.number*self.MN)

def read_xyz(self,filename: str):

number = 0

# Extract the number of atoms

for tokens in get_tokens(filename):

if len(tokens) == 1:

number = int(tokens[0])

break # Exit the loop after processing the first line

# Allocate memory for arrays storing atomic data

box = [0.0] * 6

x = [0.0] * number

y = [0.0] * number

z = [0.0] * number

x0 = [0.0] * number

y0 = [0.0] * number

z0 = [0.0] * number

for tokens in get_tokens(filename):

if len(tokens) == 9:

box[0] = float(tokens[0])

box[1] = float(tokens[4])

box[2] = float(tokens[8])

for n in range(3):

box[3 + n] = box[n] * 0.5

# Temporary lists are used to store rows containing 4 tokens (assuming atomic information)

length_4_lists = []

for tokens in get_tokens(filename):

if len(tokens) == 4:

length_4_lists.append(tokens)

# Process rows containing information about atoms (assuming 4 tokens per atom)

for n in range(number):

x[n] =x0[n] = float(length_4_lists[n][1])

y[n] =y0[n] = float(length_4_lists[n][2])

z[n] = z0[n] = float(length_4_lists[n][3])

return number, box, x, y, z, x0, y0, z0

def find_neighbor_on2(self):

"""

This function finds and stores the neighbors of each atom within a specified cutoff distance

using an O(N^2) algorithm.

Returns:

None: Modifies the neighbor lists and neighbor counts in-place.

"""

# Square of the neighbor cutoff distance

cutoff_square = self.cutoffNeighbor ** 2

# Initialize neighbor lists with zeros

self.NN = [0] * self.number

# Loop over all pairs of atoms (i, j)

for i in range(self.number):

for j in range(self.number):

if j == i:

continue

xij = self.x[j] - self.x[i]

yij = self.y[j] - self.y[i]

zij = self.z[j] - self.z[i]

# Apply Minimum Image Convention (MIC)

xij, yij, zij = self.apply_mic(xij, yij, zij)

distance_square = xij * xij + yij * yij + zij * zij

# Check if the distance is within the cutoff radius

if distance_square < cutoff_square:

# Add atom j to the neighbor list of atom i

self.NL[i * self.MN + self.NN[i]] = j

self.NN[i] += 1

# Check if the number of neighbors exceeds the maximum allowed number of neighbors (MN)

if self.NN[i] > self.MN:

print(f"Error: number of neighbors for atom {i} exceeds {self.MN}")

def apply_mic(self, xij, yij, zij) -> None:

"""

Apply Minimum Image Conventions (MIC) to individual coordinates.

"""

xij = apply_mic_one(self.box[0], self.box[3], xij)

yij = apply_mic_one(self.box[1], self.box[4], yij)

zij = apply_mic_one(self.box[2], self.box[5], zij)

return xij, yij, zij

def find_force_tersoff(self):

"""

alculate the Tersoff potential and forces.

"""

# Initialize potential energy

self.pe = 0.0

# Initialize force

self.fx = [0.0] * self.number

self.fy = [0.0] * self.number

self.fz = [0.0] * self.number

# Loop over all atoms

for n1 in range(self.number):

neighbors = self.NL[n1 * self.MN : n1 * self.MN + self.NN[n1]]

for i1,n2 in enumerate(neighbors):

# Calculate the relative position and distance between atom n1 and neighbor n2

x12 = self.x[n2] - self.x[n1]

y12 = self.y[n2] - self.y[n1]

z12 = self.z[n2] - self.z[n1]

x12, y12, z12 = self.apply_mic(x12, y12, z12)

d12 = math.sqrt(x12 * x12 + y12 * y12 + z12 * z12)

d12inv = 1.0 / d12

# Calculate force field parameters

fc12, fcp12 = find_fc_and_fcp(d12)

fa12, fap12 = find_fa_and_fap(d12)

fr12, frp12 = find_fr_and_frp(d12)

b12 = self.b[n1 * self.MN + i1]

bp12 = self.bp[n1 * self.MN + i1]

# Initialize force vector

f12 = [0.0, 0.0, 0.0]

factor1 = -b12 * fa12 + fr12

factor2 = -b12 * fap12 + frp12

factor3 = (fcp12 * factor1 + fc12 * factor2) / d12

# Calculate the force components

f12[0] += x12 * factor3 * 0.5

f12[1] += y12 * factor3 * 0.5

f12[2] += z12 * factor3 * 0.5

# Loop over other neighbors of atom n1

for i2 in range(self.NN[n1]):

n3 = self.NL[n1 * self.MN + i2]

if n3 == n2:

continue

# Calculate the relative position and distance between atom n1 and neighbor n3

x13 = self.x[n3] - self.x[n1]

y13 = self.y[n3] - self.y[n1]

z13 = self.z[n3] - self.z[n1]

x13, y13, z13 = self.apply_mic(x13, y13, z13)

d13 = math.sqrt(x13 * x13 + y13 * y13 + z13 * z13)

# Calculate the three-body term

fc13 = find_fc(d13)

fa13 = find_fa(d13)

bp13 = self.bp[n1 * self.MN + i2]

cos123 = (x12 * x13 + y12 * y13 + z12 * z13) / (d12 * d13)

g123, gp123 = find_g_and_gp(cos123)

cos_x = x13 / (d12 * d13) - x12 * cos123 / (d12 * d12)

cos_y = y13 / (d12 * d13) - y12 * cos123 / (d12 * d12)

cos_z = z13 / (d12 * d13) - z12 * cos123 / (d12 * d12)

factor123a = (-bp12 * fc12 * fa12 * fc13 -

bp13 * fc13 * fa13 * fc12) * gp123

factor123b = -bp13 * fc13 * fa13 * fcp12 * g123 * d12inv

# Update force components

f12[0] += (x12 * factor123b + factor123a * cos_x) * 0.5

f12[1] += (y12 * factor123b + factor123a * cos_y) * 0.5

f12[2] += (z12 * factor123b + factor123a * cos_z) * 0.5

# Update potential energy and forces

self.pe += factor1 * fc12 * 0.5

self.fx[n1] += f12[0]

self.fy[n1] += f12[1]

self.fz[n1] += f12[2]

self.fx[n2] -= f12[0]

self.fy[n2] -= f12[1]

self.fz[n2] -= f12[2]

def find_b_and_bp(self):

"""

Calculate the bond order term (b) and its derivative (bp).

"""

# Parameter for Tersoff potential

beta = 1.1000e-6

n = 0.78734

minus_half_over_n = -0.5 / n

# Loop over all atoms

for n1 in range(self.number):

# Loop over neighbors of atom n1

for i1 in range(self.NN[n1]):

n2 = self.NL[n1 * self.MN + i1]# Neighbor n2 of atom n1

x12 = self.x[n2] - self.x[n1]

y12 = self.y[n2] - self.y[n1]

z12 = self.z[n2] - self.z[n1]

x12, y12, z12 = self.apply_mic(x12, y12, z12)

d12 = math.sqrt(x12 * x12 + y12 * y12 + z12 * z12)

# Initialize the zeta parameter

zeta = 0.0

# Loop over other neighbors of atom n1

for i2 in range(self.NN[n1]):

n3 = self.NL[n1 * self.MN + i2]

if n3 == n2:

continue

# Calculate the relative position and distance between atom n1 and neighbor n3

x13 = self.x[n3] - self.x[n1]

y13 = self.y[n3] - self.y[n1]

z13 = self.z[n3] - self.z[n1]

x13, y13, z13 = self.apply_mic(x13, y13, z13)

d13 = math.sqrt(x13 * x13 + y13 * y13 + z13 * z13)

# Calculate the cosine of the angle between bonds n1-n2 and n1-n3

cos = (x12 * x13 + y12 * y13 + z12 * z13) / (d12 * d13)

# Calculate the cutoff function for d13

fc13 = find_fc(d13)

# Calculate the angle-dependent function

g123 = find_g(cos)

# Accumulate the zeta parameter

zeta += fc13 * g123

# Calculate the bond order term b12

bzn = math.pow(beta * zeta, n)

b12 = math.pow(1.0 + bzn, minus_half_over_n)

# Store the bond order term

self.b[n1 * self.MN + i1] = b12

# Calculate the derivative of the bond order term bp12

self.bp[n1 * self.MN + i1] = -b12 * bzn * 0.5 / ((1.0 + bzn) * zeta)

def find_force(self):

"""

计算原子的势能和力。

参数:

atom (Atom): Atom对象，包含原子信息。

"""

self.find_b_and_bp()

self.find_force_tersoff()

def find_neighbor_on2(self):

"""

This function finds and stores the neighbors of each atom within a specified cutoff distance

using an O(N^2) algorithm.

Returns:

None: Modifies the neighbor lists and neighbor counts in-place.

"""

# Square of the neighbor cutoff distance

cutoff_square = self.cutoffNeighbor ** 2

# Initialize neighbor lists with zeros

self.NN = [0] * self.number

# Loop over all pairs of atoms (i, j)

for i in range(self.number):

x1, y1, z1 = self.x[i], self.y[i], self.z[i]

for j in range(self.number):

if j == i:

continue

xij = self.x[j] - x1

yij = self.y[j] - y1

zij = self.z[j] - z1

# Apply Minimum Image Convention (MIC)

xij, yij, zij = self.apply_mic(xij, yij, zij)

distance_square = xij * xij + yij * yij + zij * zij

# Check if the distance is within the cutoff radius

if distance_square < cutoff_square:

if self.NN[i] < self.MN:

# Add atom j to the neighbor list of atom i

self.NL[i * self.MN + self.NN[i]] = j

self.NN[i] += 1

else:

print(f"Error: number of neighbors for atom {i} exceeds {self.MN}")

exit(1)

def get_tokens(filename):

"""

Reads a text file and generates markers (words) line by line.

This function simulates the behavior of reading markers from a file by taking the filename as input and splitting the contents of the file into words.

Takes the filename as input and splits the contents into words. It removes leading/trailing whitespace from each line and skips empty lines.

Parameters:

filename:path of the text file to be read.

result:a list of words (phrases) in each line of the file.

"""

try:

with open(filename, 'r') as input_file:

# Iterate over each line in the file

for line in input_file:

line = line.strip() # Remove leading/trailing spaces from a line

if line:

yield line.split()# Skip blank lines (only blank lines)

# Handle file not found cases

except FileNotFoundError:

print(f"Error: Could not open file {filename}")

def apply_mic_one(box_length, half_box_length, distance):

"""

Apply Minimum Image Conventions (MIC) to individual coordinates.

Parameters

box_length:length of the simulation box in this dimension.

half_box_length:half the length of the simulation box (pre-calculated for efficiency).

distance:the distance value to be wrapped

"""

if distance < -half_box_length:

distance += box_length

elif distance > half_box_length:

distance -= box_length

return distance

def apply_pbc_one(sx):

"""

Apply periodic boundary conditions (PBC) to individual coordinates.

parameter:

box_length: the length of the simulation box in this dimension.

coordinate: the coordinate value to be encapsulated

"""

# If the coordinate is less than zero, the particle has moved beyond the lower boundary of the box.

if sx < 0.0:

sx += 1.0

# If the coordinate is larger than box_length, the particle has moved beyond the upper boundary of the box.

elif sx > 1.0:

sx -= 1.0

pass

def find_fr_and_frp(d12):

"""

Calculate the values of fr and frp for a given distance.

Parameters:

d12 (float): Distance value.

"""

a = 1.8308e3

lambda_val = 2.4799

fr = a * math.exp(-lambda_val * d12)

frp = -lambda_val * fr

return fr, frp

def find_fa_and_fap(d12):

"""

Calculate the values of fa and fap for a given distance.

Parameters:

d12 (float): Distance value.

"""

b = 471.18

mu = 1.7322

fa = b * math.exp(-mu * d12)

fap = -mu * fa

return fa, fap

def find_fa(d12):

"""

Calculate the value of fa for a given distance.

Parameters:

d12 (float): Distance value.

"""

b = 471.18

mu = 1.7322

fa = b * math.exp(-mu * d12)

return fa

def find_fc_and_fcp(d12):

"""

Calculate the values of fc and fcp for a given distance.

Parameters:

d12 (float): Distance value.

"""

r1 = 2.7

r2 = 3.0

pi = math.pi

# Pre-compute a factor involving pi and the cutoff radii

pi_factor = pi / (r2 - r1)

# Determine the value of fc based on the distance d12

if d12 < r1:

fc = 1.0

fcp = 0.0

elif d12 < r2:

# Between the first and second cutoff radii, fc and fcp are computed using cosine function

fc = math.cos(pi_factor * (d12 - r1)) * 0.5 + 0.5

fcp = -math.sin(pi_factor * (d12 - r1)) * pi_factor * 0.5

else:

fc = 0.0# Beyond the second cutoff radius, fc is 0

fcp = 0.0 # The derivative fcp is 0

return fc, fcp

def find_fc(d12):

"""

Calculate the value of fc for a given distance.

Parameters:

d12 (float): Distance value.

"""

r1 = 2.7

r2 = 3.0

pi = math.pi

# Pre-compute a factor involving pi and the cutoff radii

pi_factor = pi / (r2 - r1)

# Determine the value of fc based on the distance d12

if d12 < r1:

fc = 1.0

elif d12 < r2:

# Between the first and second cutoff radii, fc and fcp are computed using cosine function

fc = math.cos(pi_factor * (d12 - r1)) * 0.5 + 0.5

else:

fc = 0.0 # Beyond the second cutoff radius, fc is 0

return fc

def find_g_and_gp(cos_value):

"""

Calculate the values of g and its derivative gp for a given cosine value.

Parameters:

cos_value (float): Cosine value.

"""

# Parameter for the Tersoff potential

c = 1.0039e5

d = 16.217

h = -0.59825

c2 = c * c

d2 = d * d

# Pre-compute a constant factor

c2overd2 = c2 / d2

# Calculate intermediate term

temp = d2 + (cos_value - h) * (cos_value - h)

# Calculate g and its derivative gp

g = 1.0 + c2overd2 - c2 / temp

gp = 2.0 * c2 * (cos_value - h) / (temp * temp)

return g, gp

def find_g(cos_value):

"""

Calculate the value of g for a given cosine value.

Parameters:

cos_value (float): Cosine value.

"""

# Parameter for the Tersoff potential

c = 1.0039e5

d = 16.217

h = -0.59825

c2 = c * c

d2 = d * d

c2overd2 = c2 / d2

# Calculate intermediate term

temp = d2 + (cos_value - h) * (cos_value - h)

# Calculate g

g = 1.0 + c2overd2 - c2 / temp

return g

def main():

atom = Atom()

atom.find_neighbor_on2()

print("neighbor list is built.")

# Calculate analytical forces and write to file

with open("force.out", "w") as ofile:

atom.find_force()

for n in range(atom.number):

ofile.write(f"{atom.fx[n]:.16f} {atom.fy[n]:.16f} {atom.fz[n]:.16f}\n")

print("analytical force is calculated.")

# Calculate finite-difference forces and append to file

delta = 2.0e-5

with open("force.out", "a") as ofile:

for n in range(atom.number):

# Perturb x coordinate and calculate forces

atom.x[n] = atom.x0[n] + delta

atom.find_force()

pe_positive = atom.pe

atom.x[n] = atom.x0[n] - delta

atom.find_force()

pe_negative = atom.pe

atom.x[n] = atom.x0[n]

# Calculate finite-difference force derivative and write to file

ofile.write(f"{(pe_negative - pe_positive) / (delta * 2.0):.16f} ")

# Perturb y coordinate and calculate forces

atom.y[n] = atom.y0[n] + delta

atom.find_force()

pe_positive = atom.pe

atom.y[n] = atom.y0[n] - delta

atom.find_force()

pe_negative = atom.pe

atom.y[n] = atom.y0[n]

ofile.write(f"{(pe_negative - pe_positive) / (delta * 2.0):.16f} ")

# Perturb z coordinate and calculate forces

atom.z[n] = atom.z0[n] + delta

atom.find_force()

pe_positive = atom.pe

atom.z[n] = atom.z0[n] - delta

atom.find_force()

pe_negative = atom.pe

atom.z[n] = atom.z0[n]

ofile.write(f"{(pe_negative - pe_positive) / (delta * 2.0):.16f}\n")

print("finite-difference force is calculated.")

main()

代码

文本

第二版（numpy版本，简洁）

代码

文本

[ ]

999

100

101

102

103

104

105

106

107

108

109

110

111

112

113

114

115

116

117

118

119

120

121

122

123

124

125

126

127

128

129

130

131

132

133

134

135

136

137

138

139

140

141

142

143

144

145

146

147

148

149

150

151

152

153

154

155

156

157

158

159

160

161

162

163

164

165

166

167

168

169

170

171

172

173

174

175

176

177

178

179

180

181

182

183

184

185

186

187

188

189

190

191

192

193

194

195

196

197

198

199

200

201

202

203

204

205

206

207

208

209

210

211

212

213

214

215

216

217

218

219

220

221

222

223

224

225

226

227

228

229

230

231

232

233

234

235

236

237

238

239

240

241

242

243

244

245

246

247

248

249

250

251

252

253

254

255

256

257

258

259

260

261

262

263

264

265

266

267

268

269

270

271

272

273

274

275

276

277

278

279

280

281

282

283

284

285

286

287

288

289

290

291

292

293

294

295

296

297

298

299

300

301

302

303

304

305

306

307

308

309

310

311

312

313

314

315

316

317

318

319

320

321

322

323

324

325

326

327

328

329

330

331

332

333

334

335

336

337

338

339

340

341

342

343

344

345

346

347

348

349

350

351

352

353

354

355

356

357

358

359

360

361

362

363

364

365

366

367

368

369

370

371

372

373

374

375

376

377

378

379

380

381

382

383

384

385

386

387

388

389

390

391

392

393

394

395

396

397

398

399

400

401

402

403

404

405

406

407

408

409

410

411

412

413

414

415

416

417

418

419

420

421

422

423

424

425

426

427

428

429

430

431

432

433

434

435

436

437

438

439

440

441

442

443

444

445

446

447

448

449

450

451

452

453

454

455

456

457

458

459

460

461

462

463

464

465

466

467

468

469

470

471

472

473

474

475

476

477

478

479

480

481

482

483

484

485

486

487

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

511

512

513

514

515

516

517

518

519

520

521

522

523

524

525

526

527

528

529

530

531

532

533

534

535

536

537

538

539

540

541

542

543

544

545

546

547

548

549

550

551

552

553

554

555

556

557

558

559

560

561

562

563

564

565

566

567

568

569

570

571

572

573

574

575

576

577

578

579

580

581

582

583

584

585

586

587

588

589

590

591

592

593

594

595

596

597

598

599

600

601

602

603

604

605

606

607

608

609

610

›

⌄

#第二版（numpy版本，简洁）

import numpy as np

from time import time

def timer(func):

def func_wrapper(*args, **kwargs):

time_start = time()

result = func(*args, **kwargs)

time_end = time()

time_spend = time_end - time_start

print('%s cost time: %.3f s' % (func.__name__, time_spend))

return result

return func_wrapper

class Units:

k_B: float = 8.617343e-5

TIME_UNIT_CONVERSION: float = 1.018051e+1

class LJParameters:

def __init__(self,epsilon: float = 1.032e-2,sigma: float = 3.405,cutoff: float = 9.0):

# LJ势参数初始化

self.epsilon: float = epsilon

self.sigma: float = sigma

self.cutoff: float = cutoff

self.cutoffSquare: float = cutoff**2

self.sigma3: float = sigma**3

self.sigma6: float = sigma**6

self.sigma12: float = sigma**12

self.e24s6: float = 24.0 * epsilon * self.sigma6

self.e48s12: float = 48.0 * epsilon * self.sigma12

self.e4s6: float = 4.0 * epsilon * self.sigma6

self.e4s12: float = 4.0 * epsilon * self.sigma12

class TersoffParameters:

def __init__(self, A: float=1830.8, B: float=471.18, lamb: float=2.4799, mu: float=1.7322, beta: float=1.1e-6,

n: float=0.78734, c: float=1.0039e5, d: float=16.217, h: float=-0.59825,

R: float=2.7, S: float=3.0) -> None:

self.A = A

self.B = B

self.lamb = lamb

self.mu = mu

self.beta = beta

self.n = n

self.c = c

self.d = d

self.h = h

self.R = R

self.S = S

self.minus_half_over_n = -0.5 / n

self.pi = np.pi

self.pi_factor = np.pi / (S-R)

self.c2 = c**2

self.d2 = d**2

self.c2overd2 = self.c2 / self.d2

pass

class Atom:

def __init__(self, filename: str='xyz.in',

cutoffNeighbor: float=10.0, MaxNeighbor: int=1000,

neighborFlag: int=0, tersoff: TersoffParameters=None) -> None:

'''

读取xyz文件，初始化原子的坐标，质量，速度，势能，动能

:param filename: xyz文件名

:returns: None

'''

self.filename = filename

# 读取xyz文件

self.number, self.box, self.coords, self.boxInv = self.parseXyzFile(self.filename)

self.mass = np.full(self.number, 28)

# 初始化原子的力、速度、势能、动能

self.forces = np.zeros((self.number, 3))

self.velocities = np.zeros((self.number, 3))

self.pe = 0.0

self.ke = self.getKineticEnergy()

# Neighbor list parameters

self.MaxNeighbor: int = MaxNeighbor

self.cutoffNeighbor: float = cutoffNeighbor

self.NeighborNumber: np.ndarray = np.zeros(self.number, dtype=int)

self.NeighborList: np.ndarray = np.zeros((self.number, self.MaxNeighbor), dtype=int)

self.NeighborFlag: int = neighborFlag

self.numUpdates: int = 0

self.coord_ref: np.ndarray = np.zeros((self.number, 3))

# Tersoff potential

self.b = np.zeros((self.number, self.MaxNeighbor))

self.bp = np.zeros((self.number, self.MaxNeighbor))

if tersoff is not None:

self.tersoff = tersoff

def parseXyzFile(self, filename: str) -> tuple[int, np.ndarray, np.ndarray, np.ndarray]:

'''

读取xyz文件，返回原子数，盒子大小，原子坐标，原子质量

:param filename: xyz文件名

:returns: 原子数，盒子大小，原子坐标，原子质量, 盒子逆矩阵

'''

with open(filename, 'r') as f:

# 读取第一行的原子数

number = int(f.readline().strip())

# 初始化坐标和质量

coords = np.zeros((number, 3))

# 读取box的大小

box = [float(x) for x in f.readline().split()]

box = np.array(box).reshape((3,3))

# 遍历读取原子坐标和质量

for i in range(number):

line = f.readline().split()

coords[i] = [float(line[1]), float(line[2]), float(line[3])]

# 求box的逆矩阵

boxInv = np.linalg.inv(box)

return number, box, coords, boxInv

def getKineticEnergy(self) -> float:

'''

返回体系的动能：K = 0.5 * sum(m_i * v_i^2)

:returns: 动能

'''

return 0.5 * np.sum(np.sum(self.velocities**2, axis=1) * self.mass)

def initializeVelocities(self, temperature: float) -> None:

'''

初始化给定温度下的原子速度

:param temperature: 目标温度

:returns: None

'''

# 计算总质量

totalMass = np.sum(self.mass)

# 初始化速度，速度大小在-1到1之间均匀分布

self.velocities = np.random.uniform(-1, 1, (self.number, 3))

# 计算质心速度

centerOfMassVelocity = np.sum(self.velocities * self.mass[:, np.newaxis], axis=0) / totalMass

# 去除质心速度，防止整体运动

self.velocities -= centerOfMassVelocity

# 调整速度，使得温度符合目标温度

self.scaleVelocities(temperature)

def scaleVelocities(self, temperature: float) -> None:

'''

调整速度，使得温度符合目标温度

:param temperature: 目标温度

:returns: None

'''

# 计算当前动能，并计算当前温度

self.ke = self.getKineticEnergy()

currentTemperature = 2.0 * self.ke / (3.0 * Units.k_B * self.number)

# 计算调整因子

scalingFactor = np.sqrt(temperature / currentTemperature)

# 调整速度

self.velocities *= scalingFactor

def applyMic(self, rij: np.ndarray) -> np.ndarray:

'''

对于给定的两个原子间的距离，应用最小镜像约定

:param rij: 两个原子间的距离

:returns: None

'''

# rij转换为分数坐标

rijFractional = np.dot(rij, self.boxInv)

# 对于每一个维度，如果分数坐标小于-0.5，就加1，如果分数坐标大于0.5，就减1

for i in range(3):

if rijFractional[i] < -0.5:

rijFractional[i] += 1.0

elif rijFractional[i] > 0.5:

rijFractional[i] -= 1.0

# 转换为笛卡尔坐标

rij = np.dot(rijFractional, self.box)

return rij

def getForce_lj(self, lj: LJParameters) -> None:

self.pe =0.0

self.forces = np.zeros((self.number, 3))

if self.NeighborFlag == 0:

index = np.triu_indices(self.number, 1)

rij = self.coords[index[1]] - self.coords[index[0]]

rij_frac = np.dot(rij, self.boxInv)

rij_frac = np.where(rij_frac < -0.5, rij_frac + 1.0, rij_frac)

rij_frac = np.where(rij_frac > 0.5, rij_frac - 1.0, rij_frac)

rij = np.dot(rij_frac, self.box)

r2 = np.sum(rij**2, axis=1)

mask = r2 < lj.cutoffSquare

rij = rij[mask]

r2 = r2[mask]

r2i = 1.0 / r2

r6i = r2i**3

r8i = r6i * r2i

r12i = r6i**2

r14i = r12i * r2i

f_ij = lj.e24s6 * r8i - lj.e48s12*r14i

self.pe += np.sum(lj.e4s12*r12i - lj.e4s6*r6i)

force = f_ij[:, np.newaxis] * rij

np.add.at(self.forces, index[0][mask], force)

np.subtract.at(self.forces, index[1][mask], force)

else:

for i in range(self.number-1):

neighbors = self.NeighborList[i, :self.NeighborNumber[i]]

rij = self.coords[neighbors] - self.coords[i]

rij_frac = np.dot(rij, self.boxInv)

rij_frac = np.where(rij_frac < -0.5, rij_frac + 1.0, rij_frac)

rij_frac = np.where(rij_frac > 0.5, rij_frac - 1.0, rij_frac)

rij = np.dot(rij_frac, self.box)

r2 = np.sum(rij**2, axis=1)

mask = r2 < lj.cutoffSquare

rij = rij[mask]

r2 = r2[mask]

r2i = 1.0 / r2

r6i = r2i**3

r8i = r6i * r2i

r12i = r6i**2

r14i = r12i * r2i

f_ij = lj.e24s6 * r8i - lj.e48s12*r14i

self.pe += np.sum(lj.e4s12*r12i - lj.e4s6*r6i)

force = f_ij[:, np.newaxis] * rij

np.add.at(self.forces, [i], np.sum(force, axis=0))

np.subtract.at(self.forces, neighbors[mask], force)

def applyPbc(self) -> None:

'''

对原子坐标应用周期性边界条件，支持三斜盒子

:returns: None

'''

# 转换为分数坐标

coordsFractional = np.dot(self.coords, self.boxInv)

# 对于每一个维度，如果分数坐标小于0，就加1，如果分数坐标大于1，就减1

for i in range(3):

coordsFractional[:, i] -= np.floor(coordsFractional[:, i])

# 转换为笛卡尔坐标

self.coords = np.dot(coordsFractional, self.box)

def findNeighborON2(self):

cutoffSquare = self.cutoffNeighbor**2

index = np.triu_indices(self.number, 1)

rij = self.coords[index[1]] - self.coords[index[0]]

rij_frac = np.dot(rij, self.boxInv)

rij_frac = np.where(rij_frac < -0.5, rij_frac + 1.0, rij_frac)

rij_frac = np.where(rij_frac > 0.5, rij_frac - 1.0, rij_frac)

rij = np.dot(rij_frac, self.box)

r2 = np.sum(rij**2, axis=1)

mask = r2 < cutoffSquare

pairs = np.column_stack((index[0][mask], index[1][mask]))

# build neighbor list

for i, j in pairs:

self.NeighborList[i, self.NeighborNumber[i]] = j

self.NeighborNumber[i] += 1

self.NeighborList[j, self.NeighborNumber[j]] = i

self.NeighborNumber[j] += 1

if np.any(self.NeighborNumber > self.MaxNeighbor):

raise ValueError(f'Error: number of neighbors exceeds the maximum value {self.MaxNeighbor}')

def getThickness(self) -> np.ndarray:

'''

计算盒子的厚度

:returns: 盒子的厚度

'''

volume = np.abs(np.linalg.det(self.box))

area1 = np.linalg.norm(np.cross(self.box[0], self.box[1]))

area2 = np.linalg.norm(np.cross(self.box[1], self.box[2]))

area3 = np.linalg.norm(np.cross(self.box[2], self.box[0]))

return volume / np.array([area1, area2, area3])

def getCell(self, thickness: np.ndarray, cutoffInv: float, numCells: np.ndarray) -> np.ndarray:

'''

得到每个盒子中原子数目

:param thickness: 盒子的厚度

:param cutoffInv: 截断距离的倒数

:param numCells: 每个方向盒子的个数

:returns: 每个原子所在的盒子索引

'''

coordFractional = np.dot(self.coords, self.boxInv)

cellIndex = np.floor(coordFractional * thickness * cutoffInv).astype(int)

cellIndex = np.mod(cellIndex, numCells)

return cellIndex

def findNeighborON1(self):

cutoffInv = 1.0 / self.cutoffNeighbor

cutoffSquare = self.cutoffNeighbor**2

# 计算盒子的厚度

thickness = self.getThickness()

# 计算每个方向盒子的个数

numCells = np.floor(thickness * cutoffInv).astype(int)

totalNumCells = numCells[0] * numCells[1] * numCells[2]

# 获得每个原子所在的盒子索引

cellIndex = self.getCell(thickness, cutoffInv, numCells)

# 计算每个盒子中有哪些原子

cellAtoms = {tuple(idx): [] for idx in np.ndindex(*numCells)}

for atom_idx, cell_idx in enumerate(cellIndex):

cellAtoms[tuple(cell_idx)].append(atom_idx)

# 计算近邻盒子

offsets = np.array([[-1,0,1]]*3)

neighbor_offsets = np.array(np.meshgrid(*offsets)).T.reshape(-1, 3)

# walk through each atom

for n in range(self.number):

currentCell = tuple(cellIndex[n])

neighbor_cells = (currentCell + neighbor_offsets) % numCells

neighbor_cells = [tuple(cell) for cell in neighbor_cells]

all_neighbors = np.array([m for cell in neighbor_cells if cell in cellAtoms for m in cellAtoms[cell] if n < m])

if all_neighbors.size == 0:

continue

rij = self.coords[all_neighbors] - self.coords[n]

rij_frac = np.dot(rij, self.boxInv)

rij_frac = np.where(rij_frac < -0.5, rij_frac + 1.0, rij_frac)

rij_frac = np.where(rij_frac > 0.5, rij_frac - 1.0, rij_frac)

rij = np.dot(rij_frac, self.box)

r2 = np.sum(rij**2, axis=1)

valid_neighbors = all_neighbors[r2 < cutoffSquare]

self.NeighborList[n, :len(valid_neighbors)] = valid_neighbors

self.NeighborNumber[n] = len(valid_neighbors)

if len(valid_neighbors) > self.MaxNeighbor:

raise ValueError(f'Error: number of neighbors exceeds the maximum value {self.MaxNeighbor}')

def checkIfNeedUpdate(self) -> bool:

'''

检查是否需要更新NeighborList

:returns: 是否需要更新

'''

diff = self.coords - self.coord_ref

displacement_squared = np.sum(diff**2, axis=1)

if np.any(displacement_squared > 0.25):

return True

return False

def findNeighbor(self):

'''

更新NeighborList

'''

if self.checkIfNeedUpdate():

self.numUpdates += 1

self.applyPbc()

# 重置NeighborList

self.NeighborNumber: np.ndarray = np.zeros(self.number, dtype=int)

self.NeighborList: np.ndarray = np.zeros((self.number, self.MaxNeighbor), dtype=int)

if self.NeighborFlag == 1:

self.findNeighborON1()

elif self.NeighborFlag == 2:

self.findNeighborON2()

self.coord_ref = np.copy(self.coords)

def update(self, isStepOne: bool , dt: float) -> None:

'''

基于Verlet算法更新原子的坐标和速度

:param isStepOne: 是否是第一步

:param dt: 时间步长

:returns: None

'''

# 如果是第一步，就只更新速度的一半

halfdt = 0.5 * dt

self.velocities += halfdt * self.forces / self.mass[:, np.newaxis]

# 完全更新坐标

if isStepOne:

self.coords += dt * self.velocities

def getForce_tersoff(self) -> None:

# reset forces and potential energy

self.pe = 0.0

self.forces = np.zeros((self.number, 3))

for n1 in range(self.number):

neighbors = self.NeighborList[n1, :self.NeighborNumber[n1]]

for idx2, n2 in enumerate(neighbors):

r12 = self.coords[n2] - self.coords[n1]

r12_frac = np.dot(r12, self.boxInv)

r12_frac = np.where(r12_frac < -0.5, r12_frac + 1.0, r12_frac)

r12_frac = np.where(r12_frac > 0.5, r12_frac - 1.0, r12_frac)

r12 = np.dot(r12_frac, self.box)

d12 = np.linalg.norm(r12)

d12_inv = 1.0 / d12

fc12, fcp12, fa12, fap12, fr12, frp12 = self.find_all(d12)

b12 = self.b[n1, idx2]

bp12 = self.bp[n1, idx2]

f12 = np.zeros(3)

factor1 = -b12 * fa12 + fr12

factor2 = -b12 * fap12 + frp12

factor3 = (fcp12 * factor1 + fc12 * factor2) / d12

f12 += 0.5 * factor3 * r12

for idx3, n3 in enumerate(neighbors):

if n3 == n2:

continue

r13 = self.coords[n3] - self.coords[n1]

r13_frac = np.dot(r13, self.boxInv)

r13_frac = np.where(r13_frac < -0.5, r13_frac + 1.0, r13_frac)

r13_frac = np.where(r13_frac > 0.5, r13_frac - 1.0, r13_frac)

r13 = np.dot(r13_frac, self.box)

d13 = np.linalg.norm(r13)

fc13 = self.find_fc(d13)

fa13 = self.find_fa(d13)

bp13 = self.bp[n1, idx3]

cos123 = np.dot(r12, r13) / (d12 * d13)

g123, gp123 = self.find_g_and_gp(cos123)

cos = r13/(d12 * d13) - r12*cos123/(d12**2)

factor123a = (-bp12 * fc12 * fa12 * fc13 - bp13 * fc13 * fa13 * fc12) * gp123

factor123b = -bp13 * fc13 * fa13 * fcp12 * g123 * d12_inv

f12 += 0.5 * (r12 * factor123b + factor123a * cos)

self.pe += factor1 * fc12 * 0.5

self.forces[n1] += f12

self.forces[n2] -= f12

def find_b_and_bp(self) -> None:

for n1 in range(self.number):

neighbors = self.NeighborList[n1, :self.NeighborNumber[n1]]

for idx2, n2 in enumerate(neighbors):

coord_12 = self.coords[n2] - self.coords[n1]

coord_12_frac = np.dot(coord_12, self.boxInv)

coord_12_frac = np.where(coord_12_frac < -0.5, coord_12_frac + 1.0, coord_12_frac)

coord_12_frac = np.where(coord_12_frac > 0.5, coord_12_frac - 1.0, coord_12_frac)

coord_12 = np.dot(coord_12_frac, self.box)

d_12 = np.linalg.norm(coord_12)

zeta = 0.0

for n3 in neighbors:

if n3 == n2:

continue

coord_13 = self.coords[n3] - self.coords[n1]

coord_13_frac = np.dot(coord_13, self.boxInv)

coord_13_frac = np.where(coord_13_frac < -0.5, coord_13_frac + 1.0, coord_13_frac)

coord_13_frac = np.where(coord_13_frac > 0.5, coord_13_frac - 1.0, coord_13_frac)

coord_13 = np.dot(coord_13_frac, self.box)

d_13 = np.linalg.norm(coord_13)

cos = np.dot(coord_12, coord_13) / (d_12 * d_13)

zeta += self.find_g(cos) * self.find_fc(d_13)

bzn = np.power(self.tersoff.beta * zeta, self.tersoff.n)

b12 = np.power(1.0 + bzn, self.tersoff.minus_half_over_n)

self.b[n1, idx2] = b12

self.bp[n1, idx2] = -b12 * bzn * 0.5 / ((1.0 + bzn) * zeta)

def find_fc(self, d_12: float) -> float:

if d_12 < self.tersoff.R:

return 1.0

elif d_12 < self.tersoff.S:

return 0.5 * np.cos((self.tersoff.pi_factor * (d_12 - self.tersoff.R))) + 0.5

else:

return 0.0

def find_fa(self, d12: float) -> float:

return self.tersoff.B * np.exp(-self.tersoff.mu * d12)

def find_g_and_gp(self, cos: float) -> tuple[float, float]:

temp = self.tersoff.d2 + (cos - self.tersoff.h)**2

g = 1.0 + self.tersoff.c2overd2 - self.tersoff.c2 / temp

gp = 2.0 * self.tersoff.c2 * (cos - self.tersoff.h) / temp**2

return g, gp

def find_g(self, cos: float) -> float:

temp = self.tersoff.d2 + (cos - self.tersoff.h)**2

g = 1.0 + self.tersoff.c2overd2 - self.tersoff.c2 / temp

return g

def find_all(self, d12: float) -> tuple[float, float, float, float, float, float]:

fr = self.tersoff.A * np.exp(-self.tersoff.lamb * d12)

frp = -self.tersoff.lamb * fr

fa = self.tersoff.B * np.exp(-self.tersoff.mu * d12)

fap = -self.tersoff.mu * fa

if d12 < self.tersoff.R:

fc = 1.0

fcp = 0.0

elif d12 < self.tersoff.S:

fc = 0.5 * np.cos(self.tersoff.pi_factor * (d12 - self.tersoff.R)) + 0.5

fcp = -0.5 * self.tersoff.pi_factor * np.sin(self.tersoff.pi_factor * (d12 - self.tersoff.R))

else:

fc = 0.0

fcp = 0.0

return fc, fcp, fa, fap, fr, frp

def getForce_warp(self) -> None:

self.find_b_and_bp()

self.getForce_tersoff()

def readRun(filename: str='run.in') -> tuple[float, float, int, int]:

'''

读取run文件，返回速度，时间步长，总步数

:param filename: run文件名

:returns: 速度（即温度），时间步长，总步数, neighbor_flag

'''

with open(filename, 'r') as f:

for line in f:

if line.startswith('velocity'):

velocity = float(line.split()[1])

if line.startswith('time_step'):

time_step = float(line.split()[1]) / Units.TIME_UNIT_CONVERSION

if line.startswith('run'):

run = int(line.split()[1])

if line.startswith('neighbor_flag'):

neighbor_flag = int(line.split()[1])

return velocity, time_step, run, neighbor_flag

def main():

atom =Atom(filename='model.xyz', cutoffNeighbor=3.1, MaxNeighbor=4, neighborFlag=2, tersoff=TersoffParameters())

atom.findNeighbor()

atom.getForce_warp()

analytical = np.copy(atom.forces)

delta = 2.0e-5

f = np.zeros((atom.number, 3))

for n in range(atom.number):

atom.coords[n, 0] = atom.coord_ref[n, 0] + delta

atom.getForce_warp()

pePostive = atom.pe

atom.coords[n, 0] = atom.coord_ref[n, 0] - delta

atom.getForce_warp()

peNegative = atom.pe

atom.coords[n, 0] = atom.coord_ref[n, 0]

fx = (peNegative - pePostive) / (2.0 * delta)

atom.coords[n, 1] = atom.coord_ref[n, 1] + delta

atom.getForce_warp()

pePostive = atom.pe

atom.coords[n, 1] = atom.coord_ref[n, 1] - delta

atom.getForce_warp()

peNegative = atom.pe

atom.coords[n, 1] = atom.coord_ref[n, 1]

fy = (peNegative - pePostive) / (2.0 * delta)

atom.coords[n, 2] = atom.coord_ref[n, 2] + delta

atom.getForce_warp()

pePostive = atom.pe

atom.coords[n, 2] = atom.coord_ref[n, 2] - delta

atom.getForce_warp()

peNegative = atom.pe

atom.coords[n, 2] = atom.coord_ref[n, 2]

fz = (peNegative - pePostive) / (2.0 * delta)

f[n] = np.array([fx, fy, fz])

data = np.concatenate((analytical, f), axis=0)

np.save('force.npy', data)

if __name__ == '__main__':

main()

代码

文本

代码

文本

樊哲勇《分子动力学模拟》程序代码python实现

AI4S101

樊哲勇《分子动力学模拟》程序代码python实现AI4S101

点个赞吧

本文被以下合集收录

chuwei

更新于 2024-06-18

2 篇0 人关注