樊哲勇《分子动力学模拟》python实例 | 1.一个简单的分子动力学模拟程序

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

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

代码

文本

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

代码

文本

一个简单的分子动力学模拟程序

本节教程目录如下：
1. 原子初始化
2. 边界条件实现（周期性边界条件）
3. 原子间相互作用（**Lennard-Jones势**）
4. 运动方程的数值积分（**Verlet积分算法**）
5. 主控制函数

在本章讨论的分子动力学模拟中，没有外界对系统的干扰，所有粒子的运动完全由粒子间的相互作用力决定。从经典力学的角度看，这样的体系对应哈密顿体系。而从经典统计力学的角度看，这样的体系则属于微正则系综，即粒子数 N、体积 V 和能量 E 保持恒定的 NVE 系综。

一个典型的简单分子动力学模拟流程如下：

初始化：设置系统的初始条件，包括每个粒子的位置和速度。
时间演化：根据粒子间的相互作用规律，确定所有粒子的运动方程（即二阶常微分方程组），并通过数值积分不断更新每个粒子的坐标和速度。最终，我们得到系统在不同时间点上的离散相空间位置，即一条离散的相轨迹。
测量：利用统计力学的方法分析相轨迹中所包含的物理规律。

详细内容请参考樊哲勇老师《分子动力学模拟》

代码

文本

第一部分代码仅供展示，要想运行代码请前往第二部分python总代码实现

代码

文本

1.初始化

初始化是指确定初始的相空间点，包括各个粒子初始的坐标和速度。

坐标初始化是指为系统中每个粒子选定一个初始的位置坐标。分子动力学模拟中如何初始化位置主要取决于所要模拟的体系。例如，如要模拟固态氩，就得让各个氩原子的位置按面心立方结构排列。

最简单的速度初始化方法是产生3N个在某区间均匀分布的随机速度分量，再通过基本条件对其修正：

系统的总动量应该为零。也就是说，我们不希望系统的质心在模拟的过程中跑动。分子间作用力是所谓的内力，不会改变系统的整体动量，即系统的整体动量守恒。只要初始的整体动量为零，在分子动力学模拟的时间演化过程中整体动量将保持为零。
系统的总动能应该与所选定的初始温度对应。我们知道，在经典统计力学中，能量均分定理成立，即粒子的哈密顿量中每一个具有平方形式的能量项的统计平均值都等于kBT /2。
系统的总角动量应该为零，但这是可选条件。这是因为，对于施加周期边界条件（见下面的讲解）的体系，系统的总角动量不守恒，故初始总角动量即使非零也无妨。

代码

文本

[ ]

#Ar原子坐标初始化代码

#生成原子坐标代码（by樊老师书籍）

import numpy as np

r0 = np.array([[0.0, 0.0, 0.5, 0.5],

[0.0, 0.5, 0.0, 0.5],

[0.0, 0.5, 0.5, 0.0]]).T

n0 = r0.shape[0]

nxyz = 6 * np.array([1, 1, 1])

N = nxyz[0] * nxyz[1] * nxyz[2] * n0

a = 5.385 * np.array([1, 1, 1])

box_length = a * nxyz

r = np.zeros((N, 3))

n = 0

for nx in range(nxyz[0]):

for ny in range(nxyz[1]):

for nz in range(nxyz[2]):

for m in range(n0):

n += 1

r[n-1, :] = a * (np.array([nx, ny, nz]) + r0[m, :])

with open('xyz.in', 'w') as fid:

fid.write(f'{N}\n')

fid.write(f'{box_length[0]} {box_length[1]} {box_length[2]}\n')

for n in range(N):

fid.write(f'Ar {r[n, 0]} {r[n, 1]} {r[n, 2]} 40\n')

代码

文本

[ ]

#初始化计算变量

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.mass = self.read_xyz(self.filename)

# Initialize forces, velocities, potential energy, and kinetic energy

self.vx = [0.0] * self.number

self.vy = [0.0] * self.number

self.vz = [0.0] * self.number

self.fx = [0.0] * self.number

self.fy = [0.0] * self.number

self.fz = [0.0] * self.number

self.pe = 0.0

self.ke = self.find_kinetic_energy()

#Ar原子速度初始化代码

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

"""

Initialize the atomic velocity based on a given temperature.

"""

# Random number generator (avoid using time-based seeds in production code)

random.seed()

# Initialize the center of mass velocity to zero

center_of_mass_velocity = [0.0, 0.0, 0.0]

# Initialize total mass

total_mass = 0.0

# Assign random velocities between -1 and 1 and accumulate centers of mass #

for i in range(self.number):

#Find the total mass

total_mass += self.mass[i]

# Initial randomization of velocity for atoms at each position

self.vx[i] = -1.0 + (random.random() * 2.0)

self.vy[i] = -1.0 + (random.random() * 2.0)

self.vz[i] = -1.0 + (random.random() * 2.0)

# Find the momentum in each direction and

center_of_mass_velocity[0] += self.mass[i] * self.vx[i]

center_of_mass_velocity[1] += self.mass[i] * self.vy[i]

center_of_mass_velocity[2] += self.mass[i] * self.vz[i]

# Find the velocity of the overall center of mass

center_of_mass_velocity[0] /= total_mass

center_of_mass_velocity[1] /= total_mass

center_of_mass_velocity[2] /= total_mass

#Delete center of mass motion

for i in range(self.number):

self.vx[i] -= center_of_mass_velocity[0]

self.vy[i] -= center_of_mass_velocity[1]

self.vz[i] -= center_of_mass_velocity[2]

# Scaling speed to reach desired temperature (implementation in scale_velocity)

self.scale_velocity(temperature)

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

"""

Calculate the scaling factor based on the desired and current temperatures

Args.

TEMPERATURE:The desired temperature of the system (in kelvins).

self:Object containing atomic information (mass, velocity).

"""

# Calculation of the current temperature from the kinetic energy of the atoms at the initial randomized velocity (Maxwell distribution)

current_temperature = self.find_kinetic_energy() * 2.0 / (3.0 * Units.kB * self.number)

# Calculate scaling factors based on desired and current temperatures

scale_factor = math.sqrt(temperature / current_temperature)

#Scaling the speed of all atoms

for i in range(self.number):

self.vx[i] *= scale_factor

self.vy[i] *= scale_factor

self.vz[i] *= scale_factor

def find_kinetic_energy(self) -> float:

# This function calculates the total kinetic energy of the system and returns the value

kinetic_energy = 0.0

#Calculate the kinetic energy for each atom and then add the sum

for i in range(self.number):

# Calculate the square of speed

v_squared = self.vx[i]**2+self.vy[i]**2+self.vz[i]**2

kinetic_energy += self.mass[i] * v_squared

return 0.5 * kinetic_energy

代码

文本

2.边界条件实现（周期性边界条件）

边界条件的选取对粒子间作用力的计算也有影响。常用的边界条件有好几种，但我们这里只先讨论其中的一种：周期边界条件。同时在本书的模拟中，总是采用最小镜像约定：在计算两粒子距离时，总是取最小的可能值。

代码

文本

[ ]

#周期性边界条件

def apply_pbc(self) -> None:

"""

Applies a periodic boundary condition (PBC) to all particle coordinates in the system.

This function ensures that all particle coordinates (x, y, z) remain within the simulation box.

If they are outside the boundaries of the simulation box, they are wrapped.

"""

for i in range(self.number):

apply_pbc_one(self.box[0], self.x[i])

apply_pbc_one(self.box[1], self.y[i])

apply_pbc_one(self.box[2], self.z[i])

def apply_pbc_one(box_length, coordinate):

"""

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 coordinate < 0.0:

coordinate += box_length

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

elif coordinate > box_length:

coordinate -= box_length

pass

代码

文本

[ ]

#最小镜像约束条件

def apply_mic(self, xij, yij, zij)

"""

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 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

代码

文本

3.原子间相互作用（Lennard-Jones势）

宏观物质的性质在很大程度上是由微观粒子之间的相互作用力决定的。本章我们只介绍一个称为Lennard-Jones势的简单势函数（简称 LJ 势）。具体公式请参考樊老师2.1.4部分。

代码

文本

[ ]

#Lennard-Jones势

#Lennard-Jones势计算参数

class LJParameters:

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

# LJ potential parameter initialization

self.epsilon: float = epsilon

self.sigma: float = sigma

self.cutoff: float = cutoff

self.cutoff_squared: 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

#使用Lennard-Jones势计算原子受力

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

"""

The function traverses all pairs of atoms and calculates the LJ forces and potentials between them.

"""

# Initialize potential energy and force

self.pe = 0.0

self.fx = [0.0] * self.number

self.fy = [0.0] * self.number

self.fz = [0.0] * self.number

# Iterate over all unique pairs of atoms

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

for j in range(i + 1, self.number):

# Calculate distance vectors and utilize Apply Minimum Image Conventions

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

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

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

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

# Calculate the square distance

r2 = xij**2+ yij**2 + zij**2

# Check for interactions within the cutoff radius (skip long-distance pairs)

if r2 > lj.cutoff_squared:

continue

# Calculate inverse distance and LJ terms (optimized for efficiency)

r2_inv = 1.0 / r2

r4_inv = r2_inv * r2_inv

r6_inv = r2_inv * r4_inv

r8_inv = r4_inv * r4_inv

r12_inv = r4_inv * r8_inv

r14_inv = r6_inv * r8_inv

force_ij = lj.e24s6 * r8_inv - lj.e48s12 * r14_inv

# Update the potential energy and force on the two atoms

self.pe += lj.e4s12 * r12_inv - lj.e4s6 * r6_inv

self.fx[i] += force_ij * xij

self.fx[j] -= force_ij * xij

self.fy[i] += force_ij * yij

self.fy[j] -= force_ij * yij

self.fz[i] += force_ij * zij

self.fz[j] -= force_ij * zij

代码

文本

4.运动方程的数值积分（Verlet积分算法）

给定一个多粒子体系的初始状态（坐标和速度），根据各个粒子之间的相互作用力就可预测该体系的运动状态，即任意时刻各个粒子的坐标和速度。该预测过程本质上就是对运动方程的数值积分。本章使用Verlet积分算法，具体公式请参考樊老师1.1.3部分。

代码

文本

[ ]

#使用Verlet积分算法

def integrate(self, bool_isStepOnen, time_step) -> None:

"""

Performs a velocity Verlet integration step for a system of atoms.

This function updates the velocity and position of the atoms according to the Verlet integration scheme.

Parameters

bool_isStepOnen:Boolean flag indicating whether this is the first step (affects position updates).

time_step:time_step: the time step of the integral.

"""

# Pre-calculate half-time steps for efficiency

time_step_half = time_step * 0.5

for i in range(self.number):

# Updated speed based on force and half-step times

inverse_mass = 1.0 / self.mass[i]

self.vx[i] += self.fx[i] * inverse_mass * time_step_half

self.vy[i] += self.fy[i] * inverse_mass * time_step_half

self.vz[i] += self.fz[i] * inverse_mass * time_step_half

# Update position (only on full step based on flags)

if bool_isStepOnen:

self.x[i] += self.vx[i] * time_step

self.y[i] += self.vy[i] * time_step

self.z[i] += self.vz[i] * time_step

#在main主函数中每一步的迭代更新

for step in range(num_steps):

atom.apply_pbc()#Apply periodic boundary conditions to all atoms

atom.integrate(True, time_step) # Step 1: Update speed and location (half-step)

atom.find_force(lj) # Step 2: Calculate the force acting on each atom

atom.integrate(False, time_step) # Step 3: Update speed and location (full step)

代码

文本

5.主控制函数

主函数控制函数步骤如下：
1. 从文件中读取模拟参数（步数、时间步长、温度）
2. 从单独的 XYZ 文件中读取原子数据（位置、质量）
3. 根据温度初始化原子速度
4. 启动计时器测量模拟时间
5. 打开一个输出文件，用于写入模拟数据 (thermo.out)
6. 按指定步数运行主模拟循环： - 对所有原子应用周期性边界条件 (PBC)
- 使用 Verlet 算法对位置和速度进行积分
- 计算作用在每个原子上的力
- 再次积分更新位置和速度
- 每 Ns 步输出数据（温度、动能、势能）
7.停止计时器并打印模拟耗时。

代码

文本

[ ]

def main():

"""

The main function that performs molecular dynamics simulations.

The function performs the following steps:

1. reads the simulation parameters (number of steps, time step, temperature) from the file.

2. Reads atomic data (position, mass) from a separate XYZ file.

3. Initialize atomic velocity based on temperature.

4. Start a timer to measure the simulation time.

5. Open an output file for writing simulation data (thermo.out).

6. Run the main simulation loop for the specified number of steps:

- Apply periodic boundary conditions (PBC) to all atoms.

- Integrate position and velocity using the Verlet algorithm (step 1).

- Calculate the force acting on each atom.

- Integrate position and velocity again (step 3).

- (Optional) Output data (temperature, kinetic energy, potential energy) every Ns steps.

7. Stop the timer and print the simulated elapsed time.

"""

# Read simulation parameters from file (assumed format)

num_steps, time_step, temperature = read_run()

# Initialize atom and LJ potential parameters

atom = Atom()

lj = LJParameters()

# Initialize atomic velocity based on temperature

atom.initialize_velocity(temperature)

#The clock is ticking

timer_start = time.time()

# Open the output file in write mode

with open("thermo.out", "w") as outfile:

# Main simulation loop

for step in range(num_steps):

atom.apply_pbc()#Apply periodic boundary conditions to all atoms

atom.integrate(True, time_step) # Step 1: Update speed and location (half-step)

atom.find_force(lj) # Step 2: Calculate the force acting on each atom

atom.integrate(False, time_step) # Step 3: Update speed and location (full step)

# Output data (every Ns steps)

if step % Units.Ns == 0:

kinetic_energy = atom.find_kinetic_energy()# Calculating kinetic energy

total_Energy = kinetic_energy + atom.pe# Calculate the total energy

temperature = kinetic_energy / (1.5 * Units.kB * atom.number)# Calculate temperature from kinetic energy

# Write data to output file

outfile.write("{:16.16f} {:16.16f} {:16.16f}\n".format(temperature, kinetic_energy, atom.pe))

print(f'{step:08d} {temperature:.2f} {kinetic_energy:2f} {atom.pe:2f} {total_Energy:2f}')

outfile.close()

timer_end = time.time()

print(f"Total time: {timer_end - timer_start:.1f} s")

代码

文本

双击即可修改

代码

文本

[ ]

代码

文本

双击即可修改

代码

文本

python总代码实现

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

第二版（numpy版本，简洁）

代码

文本

[1]

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

! mkdir -p /personal/01/ && cp -nr /bohr/AI4S101-python-rs2o/v3/chapter-02/ /personal/01/

代码

文本

[2]

import os

# 切换路径

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

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

! pwd && ls && tree -L 2

/personal/01
chapter-02
.
└── chapter-02
    ├── Frist_edition
    └── Second_edition

3 directories, 0 files

代码

文本

[3]

#查看run.in文件

! cd ./chapter-02/Frist_edition/ && cat run.in

velocity  60    # temperature in units of K
time_step 5     # time step in units of fs
run       10000 # run 10000 steps

代码

文本

[4]

#生成原子坐标代码（by樊老师书籍）

import numpy as np

r0 = np.array([[0.0, 0.0, 0.5, 0.5],

[0.0, 0.5, 0.0, 0.5],

[0.0, 0.5, 0.5, 0.0]]).T

n0 = r0.shape[0]

nxyz = 6 * np.array([1, 1, 1])

N = nxyz[0] * nxyz[1] * nxyz[2] * n0

a = 5.385 * np.array([1, 1, 1])

box_length = a * nxyz

r = np.zeros((N, 3))

n = 0

for nx in range(nxyz[0]):

for ny in range(nxyz[1]):

for nz in range(nxyz[2]):

for m in range(n0):

n += 1

r[n-1, :] = a * (np.array([nx, ny, nz]) + r0[m, :])

with open('xyz.in', 'w') as fid:

fid.write(f'{N}\n')

fid.write(f'{box_length[0]} {box_length[1]} {box_length[2]}\n')

for n in range(N):

fid.write(f'Ar {r[n, 0]} {r[n, 1]} {r[n, 2]} 40\n')

代码

文本

[7]

#查看生成的原子坐标

! cd ./chapter-02/Frist_edition/ && cat xyz.in

864
32.31 32.31 32.31
Ar 0.0 0.0 0.0 40
Ar 0.0 2.6925 2.6925 40
Ar 2.6925 0.0 2.6925 40
Ar 2.6925 2.6925 0.0 40
Ar 0.0 0.0 5.385 40
Ar 0.0 2.6925 8.0775 40
Ar 2.6925 0.0 8.0775 40
Ar 2.6925 2.6925 5.385 40
Ar 0.0 0.0 10.77 40
Ar 0.0 2.6925 13.462499999999999 40
Ar 2.6925 0.0 13.462499999999999 40
Ar 2.6925 2.6925 10.77 40
Ar 0.0 0.0 16.155 40
Ar 0.0 2.6925 18.8475 40
Ar 2.6925 0.0 18.8475 40
Ar 2.6925 2.6925 16.155 40
Ar 0.0 0.0 21.54 40
Ar 0.0 2.6925 24.232499999999998 40
Ar 2.6925 0.0 24.232499999999998 40
Ar 2.6925 2.6925 21.54 40
Ar 0.0 0.0 26.924999999999997 40
Ar 0.0 2.6925 29.6175 40
Ar 2.6925 0.0 29.6175 40
Ar 2.6925 2.6925 26.924999999999997 40
Ar 0.0 5.385 0.0 40
Ar 0.0 8.0775 2.6925 40
Ar 2.6925 5.385 2.6925 40
Ar 2.6925 8.0775 0.0 40
Ar 0.0 5.385 5.385 40
Ar 0.0 8.0775 8.0775 40
Ar 2.6925 5.385 8.0775 40
Ar 2.6925 8.0775 5.385 40
Ar 0.0 5.385 10.77 40
Ar 0.0 8.0775 13.462499999999999 40
Ar 2.6925 5.385 13.462499999999999 40
Ar 2.6925 8.0775 10.77 40
Ar 0.0 5.385 16.155 40
Ar 0.0 8.0775 18.8475 40
Ar 2.6925 5.385 18.8475 40
Ar 2.6925 8.0775 16.155 40
Ar 0.0 5.385 21.54 40
Ar 0.0 8.0775 24.232499999999998 40
Ar 2.6925 5.385 24.232499999999998 40
Ar 2.6925 8.0775 21.54 40
Ar 0.0 5.385 26.924999999999997 40
Ar 0.0 8.0775 29.6175 40
Ar 2.6925 5.385 29.6175 40
Ar 2.6925 8.0775 26.924999999999997 40
Ar 0.0 10.77 0.0 40
Ar 0.0 13.462499999999999 2.6925 40
Ar 2.6925 10.77 2.6925 40
Ar 2.6925 13.462499999999999 0.0 40
Ar 0.0 10.77 5.385 40
Ar 0.0 13.462499999999999 8.0775 40
Ar 2.6925 10.77 8.0775 40
Ar 2.6925 13.462499999999999 5.385 40
Ar 0.0 10.77 10.77 40
Ar 0.0 13.462499999999999 13.462499999999999 40
Ar 2.6925 10.77 13.462499999999999 40
Ar 2.6925 13.462499999999999 10.77 40
Ar 0.0 10.77 16.155 40
Ar 0.0 13.462499999999999 18.8475 40
Ar 2.6925 10.77 18.8475 40
Ar 2.6925 13.462499999999999 16.155 40
Ar 0.0 10.77 21.54 40
Ar 0.0 13.462499999999999 24.232499999999998 40
Ar 2.6925 10.77 24.232499999999998 40
Ar 2.6925 13.462499999999999 21.54 40
Ar 0.0 10.77 26.924999999999997 40
Ar 0.0 13.462499999999999 29.6175 40
Ar 2.6925 10.77 29.6175 40
Ar 2.6925 13.462499999999999 26.924999999999997 40
Ar 0.0 16.155 0.0 40
Ar 0.0 18.8475 2.6925 40
Ar 2.6925 16.155 2.6925 40
Ar 2.6925 18.8475 0.0 40
Ar 0.0 16.155 5.385 40
Ar 0.0 18.8475 8.0775 40
Ar 2.6925 16.155 8.0775 40
Ar 2.6925 18.8475 5.385 40
Ar 0.0 16.155 10.77 40
Ar 0.0 18.8475 13.462499999999999 40
Ar 2.6925 16.155 13.462499999999999 40
Ar 2.6925 18.8475 10.77 40
Ar 0.0 16.155 16.155 40
Ar 0.0 18.8475 18.8475 40
Ar 2.6925 16.155 18.8475 40
Ar 2.6925 18.8475 16.155 40
Ar 0.0 16.155 21.54 40
Ar 0.0 18.8475 24.232499999999998 40
Ar 2.6925 16.155 24.232499999999998 40
Ar 2.6925 18.8475 21.54 40
Ar 0.0 16.155 26.924999999999997 40
Ar 0.0 18.8475 29.6175 40
Ar 2.6925 16.155 29.6175 40
Ar 2.6925 18.8475 26.924999999999997 40
Ar 0.0 21.54 0.0 40
Ar 0.0 24.232499999999998 2.6925 40
Ar 2.6925 21.54 2.6925 40
Ar 2.6925 24.232499999999998 0.0 40
Ar 0.0 21.54 5.385 40
Ar 0.0 24.232499999999998 8.0775 40
Ar 2.6925 21.54 8.0775 40
Ar 2.6925 24.232499999999998 5.385 40
Ar 0.0 21.54 10.77 40
Ar 0.0 24.232499999999998 13.462499999999999 40
Ar 2.6925 21.54 13.462499999999999 40
Ar 2.6925 24.232499999999998 10.77 40
Ar 0.0 21.54 16.155 40
Ar 0.0 24.232499999999998 18.8475 40
Ar 2.6925 21.54 18.8475 40
Ar 2.6925 24.232499999999998 16.155 40
Ar 0.0 21.54 21.54 40
Ar 0.0 24.232499999999998 24.232499999999998 40
Ar 2.6925 21.54 24.232499999999998 40
Ar 2.6925 24.232499999999998 21.54 40
Ar 0.0 21.54 26.924999999999997 40
Ar 0.0 24.232499999999998 29.6175 40
Ar 2.6925 21.54 29.6175 40
Ar 2.6925 24.232499999999998 26.924999999999997 40
Ar 0.0 26.924999999999997 0.0 40
Ar 0.0 29.6175 2.6925 40
Ar 2.6925 26.924999999999997 2.6925 40
Ar 2.6925 29.6175 0.0 40
Ar 0.0 26.924999999999997 5.385 40
Ar 0.0 29.6175 8.0775 40
Ar 2.6925 26.924999999999997 8.0775 40
Ar 2.6925 29.6175 5.385 40
Ar 0.0 26.924999999999997 10.77 40
Ar 0.0 29.6175 13.462499999999999 40
Ar 2.6925 26.924999999999997 13.462499999999999 40
Ar 2.6925 29.6175 10.77 40
Ar 0.0 26.924999999999997 16.155 40
Ar 0.0 29.6175 18.8475 40
Ar 2.6925 26.924999999999997 18.8475 40
Ar 2.6925 29.6175 16.155 40
Ar 0.0 26.924999999999997 21.54 40
Ar 0.0 29.6175 24.232499999999998 40
Ar 2.6925 26.924999999999997 24.232499999999998 40
Ar 2.6925 29.6175 21.54 40
Ar 0.0 26.924999999999997 26.924999999999997 40
Ar 0.0 29.6175 29.6175 40
Ar 2.6925 26.924999999999997 29.6175 40
Ar 2.6925 29.6175 26.924999999999997 40
Ar 5.385 0.0 0.0 40
Ar 5.385 2.6925 2.6925 40
Ar 8.0775 0.0 2.6925 40
Ar 8.0775 2.6925 0.0 40
Ar 5.385 0.0 5.385 40
Ar 5.385 2.6925 8.0775 40
Ar 8.0775 0.0 8.0775 40
Ar 8.0775 2.6925 5.385 40
Ar 5.385 0.0 10.77 40
Ar 5.385 2.6925 13.462499999999999 40
Ar 8.0775 0.0 13.462499999999999 40
Ar 8.0775 2.6925 10.77 40
Ar 5.385 0.0 16.155 40
Ar 5.385 2.6925 18.8475 40
Ar 8.0775 0.0 18.8475 40
Ar 8.0775 2.6925 16.155 40
Ar 5.385 0.0 21.54 40
Ar 5.385 2.6925 24.232499999999998 40
Ar 8.0775 0.0 24.232499999999998 40
Ar 8.0775 2.6925 21.54 40
Ar 5.385 0.0 26.924999999999997 40
Ar 5.385 2.6925 29.6175 40
Ar 8.0775 0.0 29.6175 40
Ar 8.0775 2.6925 26.924999999999997 40
Ar 5.385 5.385 0.0 40
Ar 5.385 8.0775 2.6925 40
Ar 8.0775 5.385 2.6925 40
Ar 8.0775 8.0775 0.0 40
Ar 5.385 5.385 5.385 40
Ar 5.385 8.0775 8.0775 40
Ar 8.0775 5.385 8.0775 40
Ar 8.0775 8.0775 5.385 40
Ar 5.385 5.385 10.77 40
Ar 5.385 8.0775 13.462499999999999 40
Ar 8.0775 5.385 13.462499999999999 40
Ar 8.0775 8.0775 10.77 40
Ar 5.385 5.385 16.155 40
Ar 5.385 8.0775 18.8475 40
Ar 8.0775 5.385 18.8475 40
Ar 8.0775 8.0775 16.155 40
Ar 5.385 5.385 21.54 40
Ar 5.385 8.0775 24.232499999999998 40
Ar 8.0775 5.385 24.232499999999998 40
Ar 8.0775 8.0775 21.54 40
Ar 5.385 5.385 26.924999999999997 40
Ar 5.385 8.0775 29.6175 40
Ar 8.0775 5.385 29.6175 40
Ar 8.0775 8.0775 26.924999999999997 40
Ar 5.385 10.77 0.0 40
Ar 5.385 13.462499999999999 2.6925 40
Ar 8.0775 10.77 2.6925 40
Ar 8.0775 13.462499999999999 0.0 40
Ar 5.385 10.77 5.385 40
Ar 5.385 13.462499999999999 8.0775 40
Ar 8.0775 10.77 8.0775 40
Ar 8.0775 13.462499999999999 5.385 40
Ar 5.385 10.77 10.77 40
Ar 5.385 13.462499999999999 13.462499999999999 40
Ar 8.0775 10.77 13.462499999999999 40
Ar 8.0775 13.462499999999999 10.77 40
Ar 5.385 10.77 16.155 40
Ar 5.385 13.462499999999999 18.8475 40
Ar 8.0775 10.77 18.8475 40
Ar 8.0775 13.462499999999999 16.155 40
Ar 5.385 10.77 21.54 40
Ar 5.385 13.462499999999999 24.232499999999998 40
Ar 8.0775 10.77 24.232499999999998 40
Ar 8.0775 13.462499999999999 21.54 40
Ar 5.385 10.77 26.924999999999997 40
Ar 5.385 13.462499999999999 29.6175 40
Ar 8.0775 10.77 29.6175 40
Ar 8.0775 13.462499999999999 26.924999999999997 40
Ar 5.385 16.155 0.0 40
Ar 5.385 18.8475 2.6925 40
Ar 8.0775 16.155 2.6925 40
Ar 8.0775 18.8475 0.0 40
Ar 5.385 16.155 5.385 40
Ar 5.385 18.8475 8.0775 40
Ar 8.0775 16.155 8.0775 40
Ar 8.0775 18.8475 5.385 40
Ar 5.385 16.155 10.77 40
Ar 5.385 18.8475 13.462499999999999 40
Ar 8.0775 16.155 13.462499999999999 40
Ar 8.0775 18.8475 10.77 40
Ar 5.385 16.155 16.155 40
Ar 5.385 18.8475 18.8475 40
Ar 8.0775 16.155 18.8475 40
Ar 8.0775 18.8475 16.155 40
Ar 5.385 16.155 21.54 40
Ar 5.385 18.8475 24.232499999999998 40
Ar 8.0775 16.155 24.232499999999998 40
Ar 8.0775 18.8475 21.54 40
Ar 5.385 16.155 26.924999999999997 40
Ar 5.385 18.8475 29.6175 40
Ar 8.0775 16.155 29.6175 40
Ar 8.0775 18.8475 26.924999999999997 40
Ar 5.385 21.54 0.0 40
Ar 5.385 24.232499999999998 2.6925 40
Ar 8.0775 21.54 2.6925 40
Ar 8.0775 24.232499999999998 0.0 40
Ar 5.385 21.54 5.385 40
Ar 5.385 24.232499999999998 8.0775 40
Ar 8.0775 21.54 8.0775 40
Ar 8.0775 24.232499999999998 5.385 40
Ar 5.385 21.54 10.77 40
Ar 5.385 24.232499999999998 13.462499999999999 40
Ar 8.0775 21.54 13.462499999999999 40
Ar 8.0775 24.232499999999998 10.77 40
Ar 5.385 21.54 16.155 40
Ar 5.385 24.232499999999998 18.8475 40
Ar 8.0775 21.54 18.8475 40
Ar 8.0775 24.232499999999998 16.155 40
Ar 5.385 21.54 21.54 40
Ar 5.385 24.232499999999998 24.232499999999998 40
Ar 8.0775 21.54 24.232499999999998 40
Ar 8.0775 24.232499999999998 21.54 40
Ar 5.385 21.54 26.924999999999997 40
Ar 5.385 24.232499999999998 29.6175 40
Ar 8.0775 21.54 29.6175 40
Ar 8.0775 24.232499999999998 26.924999999999997 40
Ar 5.385 26.924999999999997 0.0 40
Ar 5.385 29.6175 2.6925 40
Ar 8.0775 26.924999999999997 2.6925 40
Ar 8.0775 29.6175 0.0 40
Ar 5.385 26.924999999999997 5.385 40
Ar 5.385 29.6175 8.0775 40
Ar 8.0775 26.924999999999997 8.0775 40
Ar 8.0775 29.6175 5.385 40
Ar 5.385 26.924999999999997 10.77 40
Ar 5.385 29.6175 13.462499999999999 40
Ar 8.0775 26.924999999999997 13.462499999999999 40
Ar 8.0775 29.6175 10.77 40
Ar 5.385 26.924999999999997 16.155 40
Ar 5.385 29.6175 18.8475 40
Ar 8.0775 26.924999999999997 18.8475 40
Ar 8.0775 29.6175 16.155 40
Ar 5.385 26.924999999999997 21.54 40
Ar 5.385 29.6175 24.232499999999998 40
Ar 8.0775 26.924999999999997 24.232499999999998 40
Ar 8.0775 29.6175 21.54 40
Ar 5.385 26.924999999999997 26.924999999999997 40
Ar 5.385 29.6175 29.6175 40
Ar 8.0775 26.924999999999997 29.6175 40
Ar 8.0775 29.6175 26.924999999999997 40
Ar 10.77 0.0 0.0 40
Ar 10.77 2.6925 2.6925 40
Ar 13.462499999999999 0.0 2.6925 40
Ar 13.462499999999999 2.6925 0.0 40
Ar 10.77 0.0 5.385 40
Ar 10.77 2.6925 8.0775 40
Ar 13.462499999999999 0.0 8.0775 40
Ar 13.462499999999999 2.6925 5.385 40
Ar 10.77 0.0 10.77 40
Ar 10.77 2.6925 13.462499999999999 40
Ar 13.462499999999999 0.0 13.462499999999999 40
Ar 13.462499999999999 2.6925 10.77 40
Ar 10.77 0.0 16.155 40
Ar 10.77 2.6925 18.8475 40
Ar 13.462499999999999 0.0 18.8475 40
Ar 13.462499999999999 2.6925 16.155 40
Ar 10.77 0.0 21.54 40
Ar 10.77 2.6925 24.232499999999998 40
Ar 13.462499999999999 0.0 24.232499999999998 40
Ar 13.462499999999999 2.6925 21.54 40
Ar 10.77 0.0 26.924999999999997 40
Ar 10.77 2.6925 29.6175 40
Ar 13.462499999999999 0.0 29.6175 40
Ar 13.462499999999999 2.6925 26.924999999999997 40
Ar 10.77 5.385 0.0 40
Ar 10.77 8.0775 2.6925 40
Ar 13.462499999999999 5.385 2.6925 40
Ar 13.462499999999999 8.0775 0.0 40
Ar 10.77 5.385 5.385 40
Ar 10.77 8.0775 8.0775 40
Ar 13.462499999999999 5.385 8.0775 40
Ar 13.462499999999999 8.0775 5.385 40
Ar 10.77 5.385 10.77 40
Ar 10.77 8.0775 13.462499999999999 40
Ar 13.462499999999999 5.385 13.462499999999999 40
Ar 13.462499999999999 8.0775 10.77 40
Ar 10.77 5.385 16.155 40
Ar 10.77 8.0775 18.8475 40
Ar 13.462499999999999 5.385 18.8475 40
Ar 13.462499999999999 8.0775 16.155 40
Ar 10.77 5.385 21.54 40
Ar 10.77 8.0775 24.232499999999998 40
Ar 13.462499999999999 5.385 24.232499999999998 40
Ar 13.462499999999999 8.0775 21.54 40
Ar 10.77 5.385 26.924999999999997 40
Ar 10.77 8.0775 29.6175 40
Ar 13.462499999999999 5.385 29.6175 40
Ar 13.462499999999999 8.0775 26.924999999999997 40
Ar 10.77 10.77 0.0 40
Ar 10.77 13.462499999999999 2.6925 40
Ar 13.462499999999999 10.77 2.6925 40
Ar 13.462499999999999 13.462499999999999 0.0 40
Ar 10.77 10.77 5.385 40
Ar 10.77 13.462499999999999 8.0775 40
Ar 13.462499999999999 10.77 8.0775 40
Ar 13.462499999999999 13.462499999999999 5.385 40
Ar 10.77 10.77 10.77 40
Ar 10.77 13.462499999999999 13.462499999999999 40
Ar 13.462499999999999 10.77 13.462499999999999 40
Ar 13.462499999999999 13.462499999999999 10.77 40
Ar 10.77 10.77 16.155 40
Ar 10.77 13.462499999999999 18.8475 40
Ar 13.462499999999999 10.77 18.8475 40
Ar 13.462499999999999 13.462499999999999 16.155 40
Ar 10.77 10.77 21.54 40
Ar 10.77 13.462499999999999 24.232499999999998 40
Ar 13.462499999999999 10.77 24.232499999999998 40
Ar 13.462499999999999 13.462499999999999 21.54 40
Ar 10.77 10.77 26.924999999999997 40
Ar 10.77 13.462499999999999 29.6175 40
Ar 13.462499999999999 10.77 29.6175 40
Ar 13.462499999999999 13.462499999999999 26.924999999999997 40
Ar 10.77 16.155 0.0 40
Ar 10.77 18.8475 2.6925 40
Ar 13.462499999999999 16.155 2.6925 40
Ar 13.462499999999999 18.8475 0.0 40
Ar 10.77 16.155 5.385 40
Ar 10.77 18.8475 8.0775 40
Ar 13.462499999999999 16.155 8.0775 40
Ar 13.462499999999999 18.8475 5.385 40
Ar 10.77 16.155 10.77 40
Ar 10.77 18.8475 13.462499999999999 40
Ar 13.462499999999999 16.155 13.462499999999999 40
Ar 13.462499999999999 18.8475 10.77 40
Ar 10.77 16.155 16.155 40
Ar 10.77 18.8475 18.8475 40
Ar 13.462499999999999 16.155 18.8475 40
Ar 13.462499999999999 18.8475 16.155 40
Ar 10.77 16.155 21.54 40
Ar 10.77 18.8475 24.232499999999998 40
Ar 13.462499999999999 16.155 24.232499999999998 40
Ar 13.462499999999999 18.8475 21.54 40
Ar 10.77 16.155 26.924999999999997 40
Ar 10.77 18.8475 29.6175 40
Ar 13.462499999999999 16.155 29.6175 40
Ar 13.462499999999999 18.8475 26.924999999999997 40
Ar 10.77 21.54 0.0 40
Ar 10.77 24.232499999999998 2.6925 40
Ar 13.462499999999999 21.54 2.6925 40
Ar 13.462499999999999 24.232499999999998 0.0 40
Ar 10.77 21.54 5.385 40
Ar 10.77 24.232499999999998 8.0775 40
Ar 13.462499999999999 21.54 8.0775 40
Ar 13.462499999999999 24.232499999999998 5.385 40
Ar 10.77 21.54 10.77 40
Ar 10.77 24.232499999999998 13.462499999999999 40
Ar 13.462499999999999 21.54 13.462499999999999 40
Ar 13.462499999999999 24.232499999999998 10.77 40
Ar 10.77 21.54 16.155 40
Ar 10.77 24.232499999999998 18.8475 40
Ar 13.462499999999999 21.54 18.8475 40
Ar 13.462499999999999 24.232499999999998 16.155 40
Ar 10.77 21.54 21.54 40
Ar 10.77 24.232499999999998 24.232499999999998 40
Ar 13.462499999999999 21.54 24.232499999999998 40
Ar 13.462499999999999 24.232499999999998 21.54 40
Ar 10.77 21.54 26.924999999999997 40
Ar 10.77 24.232499999999998 29.6175 40
Ar 13.462499999999999 21.54 29.6175 40
Ar 13.462499999999999 24.232499999999998 26.924999999999997 40
Ar 10.77 26.924999999999997 0.0 40
Ar 10.77 29.6175 2.6925 40
Ar 13.462499999999999 26.924999999999997 2.6925 40
Ar 13.462499999999999 29.6175 0.0 40
Ar 10.77 26.924999999999997 5.385 40
Ar 10.77 29.6175 8.0775 40
Ar 13.462499999999999 26.924999999999997 8.0775 40
Ar 13.462499999999999 29.6175 5.385 40
Ar 10.77 26.924999999999997 10.77 40
Ar 10.77 29.6175 13.462499999999999 40
Ar 13.462499999999999 26.924999999999997 13.462499999999999 40
Ar 13.462499999999999 29.6175 10.77 40
Ar 10.77 26.924999999999997 16.155 40
Ar 10.77 29.6175 18.8475 40
Ar 13.462499999999999 26.924999999999997 18.8475 40
Ar 13.462499999999999 29.6175 16.155 40
Ar 10.77 26.924999999999997 21.54 40
Ar 10.77 29.6175 24.232499999999998 40
Ar 13.462499999999999 26.924999999999997 24.232499999999998 40
Ar 13.462499999999999 29.6175 21.54 40
Ar 10.77 26.924999999999997 26.924999999999997 40
Ar 10.77 29.6175 29.6175 40
Ar 13.462499999999999 26.924999999999997 29.6175 40
Ar 13.462499999999999 29.6175 26.924999999999997 40
Ar 16.155 0.0 0.0 40
Ar 16.155 2.6925 2.6925 40
Ar 18.8475 0.0 2.6925 40
Ar 18.8475 2.6925 0.0 40
Ar 16.155 0.0 5.385 40
Ar 16.155 2.6925 8.0775 40
Ar 18.8475 0.0 8.0775 40
Ar 18.8475 2.6925 5.385 40
Ar 16.155 0.0 10.77 40
Ar 16.155 2.6925 13.462499999999999 40
Ar 18.8475 0.0 13.462499999999999 40
Ar 18.8475 2.6925 10.77 40
Ar 16.155 0.0 16.155 40
Ar 16.155 2.6925 18.8475 40
Ar 18.8475 0.0 18.8475 40
Ar 18.8475 2.6925 16.155 40
Ar 16.155 0.0 21.54 40
Ar 16.155 2.6925 24.232499999999998 40
Ar 18.8475 0.0 24.232499999999998 40
Ar 18.8475 2.6925 21.54 40
Ar 16.155 0.0 26.924999999999997 40
Ar 16.155 2.6925 29.6175 40
Ar 18.8475 0.0 29.6175 40
Ar 18.8475 2.6925 26.924999999999997 40
Ar 16.155 5.385 0.0 40
Ar 16.155 8.0775 2.6925 40
Ar 18.8475 5.385 2.6925 40
Ar 18.8475 8.0775 0.0 40
Ar 16.155 5.385 5.385 40
Ar 16.155 8.0775 8.0775 40
Ar 18.8475 5.385 8.0775 40
Ar 18.8475 8.0775 5.385 40
Ar 16.155 5.385 10.77 40
Ar 16.155 8.0775 13.462499999999999 40
Ar 18.8475 5.385 13.462499999999999 40
Ar 18.8475 8.0775 10.77 40
Ar 16.155 5.385 16.155 40
Ar 16.155 8.0775 18.8475 40
Ar 18.8475 5.385 18.8475 40
Ar 18.8475 8.0775 16.155 40
Ar 16.155 5.385 21.54 40
Ar 16.155 8.0775 24.232499999999998 40
Ar 18.8475 5.385 24.232499999999998 40
Ar 18.8475 8.0775 21.54 40
Ar 16.155 5.385 26.924999999999997 40
Ar 16.155 8.0775 29.6175 40
Ar 18.8475 5.385 29.6175 40
Ar 18.8475 8.0775 26.924999999999997 40
Ar 16.155 10.77 0.0 40
Ar 16.155 13.462499999999999 2.6925 40
Ar 18.8475 10.77 2.6925 40
Ar 18.8475 13.462499999999999 0.0 40
Ar 16.155 10.77 5.385 40
Ar 16.155 13.462499999999999 8.0775 40
Ar 18.8475 10.77 8.0775 40
Ar 18.8475 13.462499999999999 5.385 40
Ar 16.155 10.77 10.77 40
Ar 16.155 13.462499999999999 13.462499999999999 40
Ar 18.8475 10.77 13.462499999999999 40
Ar 18.8475 13.462499999999999 10.77 40
Ar 16.155 10.77 16.155 40
Ar 16.155 13.462499999999999 18.8475 40
Ar 18.8475 10.77 18.8475 40
Ar 18.8475 13.462499999999999 16.155 40
Ar 16.155 10.77 21.54 40
Ar 16.155 13.462499999999999 24.232499999999998 40
Ar 18.8475 10.77 24.232499999999998 40
Ar 18.8475 13.462499999999999 21.54 40
Ar 16.155 10.77 26.924999999999997 40
Ar 16.155 13.462499999999999 29.6175 40
Ar 18.8475 10.77 29.6175 40
Ar 18.8475 13.462499999999999 26.924999999999997 40
Ar 16.155 16.155 0.0 40
Ar 16.155 18.8475 2.6925 40
Ar 18.8475 16.155 2.6925 40
Ar 18.8475 18.8475 0.0 40
Ar 16.155 16.155 5.385 40
Ar 16.155 18.8475 8.0775 40
Ar 18.8475 16.155 8.0775 40
Ar 18.8475 18.8475 5.385 40
Ar 16.155 16.155 10.77 40
Ar 16.155 18.8475 13.462499999999999 40
Ar 18.8475 16.155 13.462499999999999 40
Ar 18.8475 18.8475 10.77 40
Ar 16.155 16.155 16.155 40
Ar 16.155 18.8475 18.8475 40
Ar 18.8475 16.155 18.8475 40
Ar 18.8475 18.8475 16.155 40
Ar 16.155 16.155 21.54 40
Ar 16.155 18.8475 24.232499999999998 40
Ar 18.8475 16.155 24.232499999999998 40
Ar 18.8475 18.8475 21.54 40
Ar 16.155 16.155 26.924999999999997 40
Ar 16.155 18.8475 29.6175 40
Ar 18.8475 16.155 29.6175 40
Ar 18.8475 18.8475 26.924999999999997 40
Ar 16.155 21.54 0.0 40
Ar 16.155 24.232499999999998 2.6925 40
Ar 18.8475 21.54 2.6925 40
Ar 18.8475 24.232499999999998 0.0 40
Ar 16.155 21.54 5.385 40
Ar 16.155 24.232499999999998 8.0775 40
Ar 18.8475 21.54 8.0775 40
Ar 18.8475 24.232499999999998 5.385 40
Ar 16.155 21.54 10.77 40
Ar 16.155 24.232499999999998 13.462499999999999 40
Ar 18.8475 21.54 13.462499999999999 40
Ar 18.8475 24.232499999999998 10.77 40
Ar 16.155 21.54 16.155 40
Ar 16.155 24.232499999999998 18.8475 40
Ar 18.8475 21.54 18.8475 40
Ar 18.8475 24.232499999999998 16.155 40
Ar 16.155 21.54 21.54 40
Ar 16.155 24.232499999999998 24.232499999999998 40
Ar 18.8475 21.54 24.232499999999998 40
Ar 18.8475 24.232499999999998 21.54 40
Ar 16.155 21.54 26.924999999999997 40
Ar 16.155 24.232499999999998 29.6175 40
Ar 18.8475 21.54 29.6175 40
Ar 18.8475 24.232499999999998 26.924999999999997 40
Ar 16.155 26.924999999999997 0.0 40
Ar 16.155 29.6175 2.6925 40
Ar 18.8475 26.924999999999997 2.6925 40
Ar 18.8475 29.6175 0.0 40
Ar 16.155 26.924999999999997 5.385 40
Ar 16.155 29.6175 8.0775 40
Ar 18.8475 26.924999999999997 8.0775 40
Ar 18.8475 29.6175 5.385 40
Ar 16.155 26.924999999999997 10.77 40
Ar 16.155 29.6175 13.462499999999999 40
Ar 18.8475 26.924999999999997 13.462499999999999 40
Ar 18.8475 29.6175 10.77 40
Ar 16.155 26.924999999999997 16.155 40
Ar 16.155 29.6175 18.8475 40
Ar 18.8475 26.924999999999997 18.8475 40
Ar 18.8475 29.6175 16.155 40
Ar 16.155 26.924999999999997 21.54 40
Ar 16.155 29.6175 24.232499999999998 40
Ar 18.8475 26.924999999999997 24.232499999999998 40
Ar 18.8475 29.6175 21.54 40
Ar 16.155 26.924999999999997 26.924999999999997 40
Ar 16.155 29.6175 29.6175 40
Ar 18.8475 26.924999999999997 29.6175 40
Ar 18.8475 29.6175 26.924999999999997 40
Ar 21.54 0.0 0.0 40
Ar 21.54 2.6925 2.6925 40
Ar 24.232499999999998 0.0 2.6925 40
Ar 24.232499999999998 2.6925 0.0 40
Ar 21.54 0.0 5.385 40
Ar 21.54 2.6925 8.0775 40
Ar 24.232499999999998 0.0 8.0775 40
Ar 24.232499999999998 2.6925 5.385 40
Ar 21.54 0.0 10.77 40
Ar 21.54 2.6925 13.462499999999999 40
Ar 24.232499999999998 0.0 13.462499999999999 40
Ar 24.232499999999998 2.6925 10.77 40
Ar 21.54 0.0 16.155 40
Ar 21.54 2.6925 18.8475 40
Ar 24.232499999999998 0.0 18.8475 40
Ar 24.232499999999998 2.6925 16.155 40
Ar 21.54 0.0 21.54 40
Ar 21.54 2.6925 24.232499999999998 40
Ar 24.232499999999998 0.0 24.232499999999998 40
Ar 24.232499999999998 2.6925 21.54 40
Ar 21.54 0.0 26.924999999999997 40
Ar 21.54 2.6925 29.6175 40
Ar 24.232499999999998 0.0 29.6175 40
Ar 24.232499999999998 2.6925 26.924999999999997 40
Ar 21.54 5.385 0.0 40
Ar 21.54 8.0775 2.6925 40
Ar 24.232499999999998 5.385 2.6925 40
Ar 24.232499999999998 8.0775 0.0 40
Ar 21.54 5.385 5.385 40
Ar 21.54 8.0775 8.0775 40
Ar 24.232499999999998 5.385 8.0775 40
Ar 24.232499999999998 8.0775 5.385 40
Ar 21.54 5.385 10.77 40
Ar 21.54 8.0775 13.462499999999999 40
Ar 24.232499999999998 5.385 13.462499999999999 40
Ar 24.232499999999998 8.0775 10.77 40
Ar 21.54 5.385 16.155 40
Ar 21.54 8.0775 18.8475 40
Ar 24.232499999999998 5.385 18.8475 40
Ar 24.232499999999998 8.0775 16.155 40
Ar 21.54 5.385 21.54 40
Ar 21.54 8.0775 24.232499999999998 40
Ar 24.232499999999998 5.385 24.232499999999998 40
Ar 24.232499999999998 8.0775 21.54 40
Ar 21.54 5.385 26.924999999999997 40
Ar 21.54 8.0775 29.6175 40
Ar 24.232499999999998 5.385 29.6175 40
Ar 24.232499999999998 8.0775 26.924999999999997 40
Ar 21.54 10.77 0.0 40
Ar 21.54 13.462499999999999 2.6925 40
Ar 24.232499999999998 10.77 2.6925 40
Ar 24.232499999999998 13.462499999999999 0.0 40
Ar 21.54 10.77 5.385 40
Ar 21.54 13.462499999999999 8.0775 40
Ar 24.232499999999998 10.77 8.0775 40
Ar 24.232499999999998 13.462499999999999 5.385 40
Ar 21.54 10.77 10.77 40
Ar 21.54 13.462499999999999 13.462499999999999 40
Ar 24.232499999999998 10.77 13.462499999999999 40
Ar 24.232499999999998 13.462499999999999 10.77 40
Ar 21.54 10.77 16.155 40
Ar 21.54 13.462499999999999 18.8475 40
Ar 24.232499999999998 10.77 18.8475 40
Ar 24.232499999999998 13.462499999999999 16.155 40
Ar 21.54 10.77 21.54 40
Ar 21.54 13.462499999999999 24.232499999999998 40
Ar 24.232499999999998 10.77 24.232499999999998 40
Ar 24.232499999999998 13.462499999999999 21.54 40
Ar 21.54 10.77 26.924999999999997 40
Ar 21.54 13.462499999999999 29.6175 40
Ar 24.232499999999998 10.77 29.6175 40
Ar 24.232499999999998 13.462499999999999 26.924999999999997 40
Ar 21.54 16.155 0.0 40
Ar 21.54 18.8475 2.6925 40
Ar 24.232499999999998 16.155 2.6925 40
Ar 24.232499999999998 18.8475 0.0 40
Ar 21.54 16.155 5.385 40
Ar 21.54 18.8475 8.0775 40
Ar 24.232499999999998 16.155 8.0775 40
Ar 24.232499999999998 18.8475 5.385 40
Ar 21.54 16.155 10.77 40
Ar 21.54 18.8475 13.462499999999999 40
Ar 24.232499999999998 16.155 13.462499999999999 40
Ar 24.232499999999998 18.8475 10.77 40
Ar 21.54 16.155 16.155 40
Ar 21.54 18.8475 18.8475 40
Ar 24.232499999999998 16.155 18.8475 40
Ar 24.232499999999998 18.8475 16.155 40
Ar 21.54 16.155 21.54 40
Ar 21.54 18.8475 24.232499999999998 40
Ar 24.232499999999998 16.155 24.232499999999998 40
Ar 24.232499999999998 18.8475 21.54 40
Ar 21.54 16.155 26.924999999999997 40
Ar 21.54 18.8475 29.6175 40
Ar 24.232499999999998 16.155 29.6175 40
Ar 24.232499999999998 18.8475 26.924999999999997 40
Ar 21.54 21.54 0.0 40
Ar 21.54 24.232499999999998 2.6925 40
Ar 24.232499999999998 21.54 2.6925 40
Ar 24.232499999999998 24.232499999999998 0.0 40
Ar 21.54 21.54 5.385 40
Ar 21.54 24.232499999999998 8.0775 40
Ar 24.232499999999998 21.54 8.0775 40
Ar 24.232499999999998 24.232499999999998 5.385 40
Ar 21.54 21.54 10.77 40
Ar 21.54 24.232499999999998 13.462499999999999 40
Ar 24.232499999999998 21.54 13.462499999999999 40
Ar 24.232499999999998 24.232499999999998 10.77 40
Ar 21.54 21.54 16.155 40
Ar 21.54 24.232499999999998 18.8475 40
Ar 24.232499999999998 21.54 18.8475 40
Ar 24.232499999999998 24.232499999999998 16.155 40
Ar 21.54 21.54 21.54 40
Ar 21.54 24.232499999999998 24.232499999999998 40
Ar 24.232499999999998 21.54 24.232499999999998 40
Ar 24.232499999999998 24.232499999999998 21.54 40
Ar 21.54 21.54 26.924999999999997 40
Ar 21.54 24.232499999999998 29.6175 40
Ar 24.232499999999998 21.54 29.6175 40
Ar 24.232499999999998 24.232499999999998 26.924999999999997 40
Ar 21.54 26.924999999999997 0.0 40
Ar 21.54 29.6175 2.6925 40
Ar 24.232499999999998 26.924999999999997 2.6925 40
Ar 24.232499999999998 29.6175 0.0 40
Ar 21.54 26.924999999999997 5.385 40
Ar 21.54 29.6175 8.0775 40
Ar 24.232499999999998 26.924999999999997 8.0775 40
Ar 24.232499999999998 29.6175 5.385 40
Ar 21.54 26.924999999999997 10.77 40
Ar 21.54 29.6175 13.462499999999999 40
Ar 24.232499999999998 26.924999999999997 13.462499999999999 40
Ar 24.232499999999998 29.6175 10.77 40
Ar 21.54 26.924999999999997 16.155 40
Ar 21.54 29.6175 18.8475 40
Ar 24.232499999999998 26.924999999999997 18.8475 40
Ar 24.232499999999998 29.6175 16.155 40
Ar 21.54 26.924999999999997 21.54 40
Ar 21.54 29.6175 24.232499999999998 40
Ar 24.232499999999998 26.924999999999997 24.232499999999998 40
Ar 24.232499999999998 29.6175 21.54 40
Ar 21.54 26.924999999999997 26.924999999999997 40
Ar 21.54 29.6175 29.6175 40
Ar 24.232499999999998 26.924999999999997 29.6175 40
Ar 24.232499999999998 29.6175 26.924999999999997 40
Ar 26.924999999999997 0.0 0.0 40
Ar 26.924999999999997 2.6925 2.6925 40
Ar 29.6175 0.0 2.6925 40
Ar 29.6175 2.6925 0.0 40
Ar 26.924999999999997 0.0 5.385 40
Ar 26.924999999999997 2.6925 8.0775 40
Ar 29.6175 0.0 8.0775 40
Ar 29.6175 2.6925 5.385 40
Ar 26.924999999999997 0.0 10.77 40
Ar 26.924999999999997 2.6925 13.462499999999999 40
Ar 29.6175 0.0 13.462499999999999 40
Ar 29.6175 2.6925 10.77 40
Ar 26.924999999999997 0.0 16.155 40
Ar 26.924999999999997 2.6925 18.8475 40
Ar 29.6175 0.0 18.8475 40
Ar 29.6175 2.6925 16.155 40
Ar 26.924999999999997 0.0 21.54 40
Ar 26.924999999999997 2.6925 24.232499999999998 40
Ar 29.6175 0.0 24.232499999999998 40
Ar 29.6175 2.6925 21.54 40
Ar 26.924999999999997 0.0 26.924999999999997 40
Ar 26.924999999999997 2.6925 29.6175 40
Ar 29.6175 0.0 29.6175 40
Ar 29.6175 2.6925 26.924999999999997 40
Ar 26.924999999999997 5.385 0.0 40
Ar 26.924999999999997 8.0775 2.6925 40
Ar 29.6175 5.385 2.6925 40
Ar 29.6175 8.0775 0.0 40
Ar 26.924999999999997 5.385 5.385 40
Ar 26.924999999999997 8.0775 8.0775 40
Ar 29.6175 5.385 8.0775 40
Ar 29.6175 8.0775 5.385 40
Ar 26.924999999999997 5.385 10.77 40
Ar 26.924999999999997 8.0775 13.462499999999999 40
Ar 29.6175 5.385 13.462499999999999 40
Ar 29.6175 8.0775 10.77 40
Ar 26.924999999999997 5.385 16.155 40
Ar 26.924999999999997 8.0775 18.8475 40
Ar 29.6175 5.385 18.8475 40
Ar 29.6175 8.0775 16.155 40
Ar 26.924999999999997 5.385 21.54 40
Ar 26.924999999999997 8.0775 24.232499999999998 40
Ar 29.6175 5.385 24.232499999999998 40
Ar 29.6175 8.0775 21.54 40
Ar 26.924999999999997 5.385 26.924999999999997 40
Ar 26.924999999999997 8.0775 29.6175 40
Ar 29.6175 5.385 29.6175 40
Ar 29.6175 8.0775 26.924999999999997 40
Ar 26.924999999999997 10.77 0.0 40
Ar 26.924999999999997 13.462499999999999 2.6925 40
Ar 29.6175 10.77 2.6925 40
Ar 29.6175 13.462499999999999 0.0 40
Ar 26.924999999999997 10.77 5.385 40
Ar 26.924999999999997 13.462499999999999 8.0775 40
Ar 29.6175 10.77 8.0775 40
Ar 29.6175 13.462499999999999 5.385 40
Ar 26.924999999999997 10.77 10.77 40
Ar 26.924999999999997 13.462499999999999 13.462499999999999 40
Ar 29.6175 10.77 13.462499999999999 40
Ar 29.6175 13.462499999999999 10.77 40
Ar 26.924999999999997 10.77 16.155 40
Ar 26.924999999999997 13.462499999999999 18.8475 40
Ar 29.6175 10.77 18.8475 40
Ar 29.6175 13.462499999999999 16.155 40
Ar 26.924999999999997 10.77 21.54 40
Ar 26.924999999999997 13.462499999999999 24.232499999999998 40
Ar 29.6175 10.77 24.232499999999998 40
Ar 29.6175 13.462499999999999 21.54 40
Ar 26.924999999999997 10.77 26.924999999999997 40
Ar 26.924999999999997 13.462499999999999 29.6175 40
Ar 29.6175 10.77 29.6175 40
Ar 29.6175 13.462499999999999 26.924999999999997 40
Ar 26.924999999999997 16.155 0.0 40
Ar 26.924999999999997 18.8475 2.6925 40
Ar 29.6175 16.155 2.6925 40
Ar 29.6175 18.8475 0.0 40
Ar 26.924999999999997 16.155 5.385 40
Ar 26.924999999999997 18.8475 8.0775 40
Ar 29.6175 16.155 8.0775 40
Ar 29.6175 18.8475 5.385 40
Ar 26.924999999999997 16.155 10.77 40
Ar 26.924999999999997 18.8475 13.462499999999999 40
Ar 29.6175 16.155 13.462499999999999 40
Ar 29.6175 18.8475 10.77 40
Ar 26.924999999999997 16.155 16.155 40
Ar 26.924999999999997 18.8475 18.8475 40
Ar 29.6175 16.155 18.8475 40
Ar 29.6175 18.8475 16.155 40
Ar 26.924999999999997 16.155 21.54 40
Ar 26.924999999999997 18.8475 24.232499999999998 40
Ar 29.6175 16.155 24.232499999999998 40
Ar 29.6175 18.8475 21.54 40
Ar 26.924999999999997 16.155 26.924999999999997 40
Ar 26.924999999999997 18.8475 29.6175 40
Ar 29.6175 16.155 29.6175 40
Ar 29.6175 18.8475 26.924999999999997 40
Ar 26.924999999999997 21.54 0.0 40
Ar 26.924999999999997 24.232499999999998 2.6925 40
Ar 29.6175 21.54 2.6925 40
Ar 29.6175 24.232499999999998 0.0 40
Ar 26.924999999999997 21.54 5.385 40
Ar 26.924999999999997 24.232499999999998 8.0775 40
Ar 29.6175 21.54 8.0775 40
Ar 29.6175 24.232499999999998 5.385 40
Ar 26.924999999999997 21.54 10.77 40
Ar 26.924999999999997 24.232499999999998 13.462499999999999 40
Ar 29.6175 21.54 13.462499999999999 40
Ar 29.6175 24.232499999999998 10.77 40
Ar 26.924999999999997 21.54 16.155 40
Ar 26.924999999999997 24.232499999999998 18.8475 40
Ar 29.6175 21.54 18.8475 40
Ar 29.6175 24.232499999999998 16.155 40
Ar 26.924999999999997 21.54 21.54 40
Ar 26.924999999999997 24.232499999999998 24.232499999999998 40
Ar 29.6175 21.54 24.232499999999998 40
Ar 29.6175 24.232499999999998 21.54 40
Ar 26.924999999999997 21.54 26.924999999999997 40
Ar 26.924999999999997 24.232499999999998 29.6175 40
Ar 29.6175 21.54 29.6175 40
Ar 29.6175 24.232499999999998 26.924999999999997 40
Ar 26.924999999999997 26.924999999999997 0.0 40
Ar 26.924999999999997 29.6175 2.6925 40
Ar 29.6175 26.924999999999997 2.6925 40
Ar 29.6175 29.6175 0.0 40
Ar 26.924999999999997 26.924999999999997 5.385 40
Ar 26.924999999999997 29.6175 8.0775 40
Ar 29.6175 26.924999999999997 8.0775 40
Ar 29.6175 29.6175 5.385 40
Ar 26.924999999999997 26.924999999999997 10.77 40
Ar 26.924999999999997 29.6175 13.462499999999999 40
Ar 29.6175 26.924999999999997 13.462499999999999 40
Ar 29.6175 29.6175 10.77 40
Ar 26.924999999999997 26.924999999999997 16.155 40
Ar 26.924999999999997 29.6175 18.8475 40
Ar 29.6175 26.924999999999997 18.8475 40
Ar 29.6175 29.6175 16.155 40
Ar 26.924999999999997 26.924999999999997 21.54 40
Ar 26.924999999999997 29.6175 24.232499999999998 40
Ar 29.6175 26.924999999999997 24.232499999999998 40
Ar 29.6175 29.6175 21.54 40
Ar 26.924999999999997 26.924999999999997 26.924999999999997 40
Ar 26.924999999999997 29.6175 29.6175 40
Ar 29.6175 26.924999999999997 29.6175 40
Ar 29.6175 29.6175 26.924999999999997 40

代码

文本

双击即可修改

代码

文本

[25]

# 尝试在文件夹下运行数据集中的主程序代码（比较慢，建议下载文件后本地电脑跑）

! cd ./chapter-02/Frist_edition/ && python main.py

00000000 59.92 6.692420 -71.255508 -64.563087
00000100 31.08 3.470918 -68.035226 -64.564308
00000200 29.34 3.276910 -67.846829 -64.569920
00000300 30.75 3.433715 -68.000568 -64.566854
00000400 31.01 3.463713 -68.036443 -64.572730
00000500 29.18 3.258472 -67.826802 -64.568329
00000600 29.56 3.300967 -67.870418 -64.569450
00000700 29.41 3.284413 -67.855155 -64.570742
00000800 30.50 3.405912 -67.976061 -64.570149
00000900 30.92 3.453024 -68.021811 -64.568786
00001000 32.39 3.617094 -68.186024 -64.568931
00001100 31.09 3.471843 -68.039995 -64.568151
00001200 30.59 3.416782 -67.986246 -64.569464
^C
Traceback (most recent call last):
  File "main.py", line 371, in <module>
    main()
  File "main.py", line 358, in main
    atom.find_force(lj)              # Step 2: Calculate the force acting on each atom
  File "main.py", line 192, in find_force
    xij,yij,zij = self.apply_mic(xij, yij, zij)  
  File "main.py", line 172, in apply_mic
    zij = apply_mic_one(self.box[2], self.box[5], zij)
KeyboardInterrupt

代码

文本

[8]

#结果可视化代码（可以运行，查看结果）

import numpy as np

import matplotlib.pyplot as plt

thermo = np.loadtxt('./chapter-02/Frist_edition/thermo.out')

timeStep = 5 / 1000 # ps

sampleInterval = 100

timeInterval = timeStep * sampleInterval

numData = thermo.shape[0]

time = np.arange(1, numData + 1) * timeInterval

totalEnergy = thermo[:, 1] + thermo[:, 2]

relativeEnergy = (totalEnergy - np.mean(totalEnergy)) / np.abs(np.mean(totalEnergy))

plt.figure(figsize=(12, 8))

plt.subplot(2, 2, 1)

plt.plot(time, thermo[:, 1], '-', linewidth=2)

plt.xlabel('Time (ps)', fontsize=15)

plt.ylabel('Kinetic Energy (eV)', fontsize=15)

plt.title('(a)', fontsize=15)

plt.subplot(2, 2, 2)

plt.plot(time, thermo[:, 2], '-', linewidth=2)

plt.xlabel('Time (ps)', fontsize=15)

plt.ylabel('Potential Energy (eV)', fontsize=15)

plt.title('(b)', fontsize=15)

plt.subplot(2, 2, 3)

plt.plot(time, totalEnergy, '-', linewidth=2)

plt.xlabel('Time (ps)', fontsize=15)

plt.ylabel('Total Energy (eV)', fontsize=15)

plt.title('(c)', fontsize=15)

plt.subplot(2, 2, 4)

plt.plot(time[1:], relativeEnergy[1:], '-', linewidth=2)

plt.xlabel('Time (ps)', fontsize=15)

plt.ylabel('Relative fluctuations', fontsize=15)

plt.ylim((-1e-3,1e-3))

plt.title('(d)', fontsize=15)

plt.tight_layout()

plt.show

<function matplotlib.pyplot.show(close=None, block=None)>

代码

文本

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

代码

文本

[ ]

999

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

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

›

⌄

#第一版本主程序代码（建议下载文件夹后在本地运行）

import math

import time

import random

def timer(func):

"""

This function is a decorator that measures the execution time of another function.

Args:

func: The function to be decorated.

Returns:

A wrapper function that measures the execution time of the decorated function.

"""

def func_wrapper(*args, **kwargs):

from time import time

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:

# Output frequency (per Ns step)

Ns = 100

# Boltzmann constant (assuming natural units)

kB = 8.617343e-5

# Time unit conversion factor (from natural units to fs)

TIME_UNIT_CONVERSION = 1.018051e+1

class LJParameters:

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

# LJ potential parameter initialization

self.epsilon: float = epsilon

self.sigma: float = sigma

self.cutoff: float = cutoff

self.cutoff_squared: 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 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.mass = self.read_xyz(self.filename)

# Initialize forces, velocities, potential energy, and kinetic energy

self.vx = [0.0] * self.number

self.vy = [0.0] * self.number

self.vz = [0.0] * self.number

self.fx = [0.0] * self.number

self.fy = [0.0] * self.number

self.fz = [0.0] * self.number

self.pe = 0.0

self.ke = self.find_kinetic_energy()

def read_xyz(self,filename: str):

"""

Parses an XYZ file containing atomic information in a simulated system.

It creates an Atom object and populates its attributes with data extracted from the file.

Args:

filename (str, optional): The filename of the input file (default: "xyz.in").

Returns:

Atom: An Atom object containing atomic information (positions, masses, box size).

"""

number = 0

# Extract atomic numbe

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

mass = [0.0] * number

box = [0.0,0.0,0.0,0.0,0.0,0.0]

x = [0.0] * number

y = [0.0] * number

z = [0.0] * number

# Extract box size

for tokens in get_tokens(filename): # Use get_tokens here

if len(tokens) == 3:

for d in range(3):

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

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

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

length_5_lists = []

for tokens in get_tokens(filename):

if len(tokens) == 5:

length_5_lists.append(tokens)

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

for n in range(number):

x[n] = float(length_5_lists[n][1])

y[n] = float(length_5_lists[n][2])

z[n] = float(length_5_lists[n][3])

mass[n] = float(length_5_lists[n][4])

return number, box, x, y, z, mass

def find_kinetic_energy(self) -> float:

# This function calculates the total kinetic energy of the system and returns the value

kinetic_energy = 0.0

#Calculate the kinetic energy for each atom and then add the sum

for i in range(self.number):

# Calculate the square of speed

v_squared = self.vx[i]**2+self.vy[i]**2+self.vz[i]**2

kinetic_energy += self.mass[i] * v_squared

return 0.5 * kinetic_energy

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

"""

Calculate the scaling factor based on the desired and current temperatures

Args.

TEMPERATURE:The desired temperature of the system (in kelvins).

self:Object containing atomic information (mass, velocity).

"""

# Calculation of the current temperature from the kinetic energy of the atoms at the initial randomized velocity (Maxwell distribution)

current_temperature = self.find_kinetic_energy() * 2.0 / (3.0 * Units.kB * self.number)

# Calculate scaling factors based on desired and current temperatures

scale_factor = math.sqrt(temperature / current_temperature)

#Scaling the speed of all atoms

for i in range(self.number):

self.vx[i] *= scale_factor

self.vy[i] *= scale_factor

self.vz[i] *= scale_factor

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

"""

Initialize the atomic velocity based on a given temperature.

"""

# Random number generator (avoid using time-based seeds in production code)

random.seed()

# Initialize the center of mass velocity to zero

center_of_mass_velocity = [0.0, 0.0, 0.0]

# Initialize total mass

total_mass = 0.0

# Assign random velocities between -1 and 1 and accumulate centers of mass #

for i in range(self.number):

#Find the total mass

total_mass += self.mass[i]

# Initial randomization of velocity for atoms at each position

self.vx[i] = -1.0 + (random.random() * 2.0)

self.vy[i] = -1.0 + (random.random() * 2.0)

self.vz[i] = -1.0 + (random.random() * 2.0)

# Find the momentum in each direction and

center_of_mass_velocity[0] += self.mass[i] * self.vx[i]

center_of_mass_velocity[1] += self.mass[i] * self.vy[i]

center_of_mass_velocity[2] += self.mass[i] * self.vz[i]

# Find the velocity of the overall center of mass

center_of_mass_velocity[0] /= total_mass

center_of_mass_velocity[1] /= total_mass

center_of_mass_velocity[2] /= total_mass

#Delete center of mass motion

for i in range(self.number):

self.vx[i] -= center_of_mass_velocity[0]

self.vy[i] -= center_of_mass_velocity[1]

self.vz[i] -= center_of_mass_velocity[2]

# Scaling speed to reach desired temperature (implementation in scale_velocity)

self.scale_velocity(temperature)

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(self, lj: LJParameters) -> None:

"""

The function traverses all pairs of atoms and calculates the LJ forces and potentials between them.

"""

# Initialize potential energy and force

self.pe = 0.0

self.fx = [0.0] * self.number

self.fy = [0.0] * self.number

self.fz = [0.0] * self.number

# Iterate over all unique pairs of atoms

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

for j in range(i + 1, self.number):

# Calculate distance vectors and utilize Apply Minimum Image Conventions

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

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

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

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

# Calculate the square distance

r2 = xij**2+ yij**2 + zij**2

# Check for interactions within the cutoff radius (skip long-distance pairs)

if r2 > lj.cutoff_squared:

continue

# Calculate inverse distance and LJ terms (optimized for efficiency)

r2_inv = 1.0 / r2

r4_inv = r2_inv * r2_inv

r6_inv = r2_inv * r4_inv

r8_inv = r4_inv * r4_inv

r12_inv = r4_inv * r8_inv

r14_inv = r6_inv * r8_inv

force_ij = lj.e24s6 * r8_inv - lj.e48s12 * r14_inv

# Update the potential energy and force on the two atoms

self.pe += lj.e4s12 * r12_inv - lj.e4s6 * r6_inv

self.fx[i] += force_ij * xij

self.fx[j] -= force_ij * xij

self.fy[i] += force_ij * yij

self.fy[j] -= force_ij * yij

self.fz[i] += force_ij * zij

self.fz[j] -= force_ij * zij

def apply_pbc(self) -> None:

"""

Applies a periodic boundary condition (PBC) to all particle coordinates in the system.

This function ensures that all particle coordinates (x, y, z) remain within the simulation box.

If they are outside the boundaries of the simulation box, they are wrapped.

"""

for i in range(self.number):

apply_pbc_one(self.box[0], self.x[i])

apply_pbc_one(self.box[1], self.y[i])

apply_pbc_one(self.box[2], self.z[i])

def integrate(self, bool_isStepOnen, time_step) -> None:

"""

Performs a velocity Verlet integration step for a system of atoms.

This function updates the velocity and position of the atoms according to the Verlet integration scheme.

Parameters

bool_isStepOnen:Boolean flag indicating whether this is the first step (affects position updates).

time_step:time_step: the time step of the integral.

"""

# Pre-calculate half-time steps for efficiency

time_step_half = time_step * 0.5

for i in range(self.number):

# Updated speed based on force and half-step times

inverse_mass = 1.0 / self.mass[i]

self.vx[i] += self.fx[i] * inverse_mass * time_step_half

self.vy[i] += self.fy[i] * inverse_mass * time_step_half

self.vz[i] += self.fz[i] * inverse_mass * time_step_half

# Update position (only on full step based on flags)

if bool_isStepOnen:

self.x[i] += self.vx[i] * time_step

self.y[i] += self.vy[i] * time_step

self.z[i] += self.vz[i] * time_step

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(box_length, coordinate):

"""

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 coordinate < 0.0:

coordinate += box_length

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

elif coordinate > box_length:

coordinate -= box_length

pass

def read_run(filename="run.in"):

"""

Parses a text file for simulation parameters.

This function reads a text file (in the assumed format) and extracts three simulation parameters:

- Number of steps

- time step

- Temperature

"""

# Initialize parameters

num_steps = None

time_step = None

temperature = None

# Loop through the rows (indicated by the tokens in the get_tokens function)

for tokens in get_tokens(filename):

# Check the first element of each line (keyword)

if tokens[0] == 'velocity':

temperature = float(tokens[1])

elif tokens[0] == 'run':

num_steps = int(tokens[1])

elif tokens[0] == 'time_step':

time_step = int(tokens[1])

time_step /= Units.TIME_UNIT_CONVERSION

# Return the extracted parameters

return num_steps, time_step, temperature

def main():

"""

The main function that performs molecular dynamics simulations.

The function performs the following steps:

1. reads the simulation parameters (number of steps, time step, temperature) from the file.

2. Reads atomic data (position, mass) from a separate XYZ file.

3. Initialize atomic velocity based on temperature.

4. Start a timer to measure the simulation time.

5. Open an output file for writing simulation data (thermo.out).

6. Run the main simulation loop for the specified number of steps:

- Apply periodic boundary conditions (PBC) to all atoms.

- Integrate position and velocity using the Verlet algorithm (step 1).

- Calculate the force acting on each atom.

- Integrate position and velocity again (step 3).

- (Optional) Output data (temperature, kinetic energy, potential energy) every Ns steps.

7. Stop the timer and print the simulated elapsed time.

"""

# Read simulation parameters from file (assumed format)

num_steps, time_step, temperature = read_run()

# Initialize atom and LJ potential parameters

atom = Atom()

lj = LJParameters()

# Initialize atomic velocity based on temperature

atom.initialize_velocity(temperature)

#The clock is ticking

timer_start = time.time()

# Open the output file in write mode

with open("thermo.out", "w") as outfile:

# Main simulation loop

for step in range(num_steps):

atom.apply_pbc()#Apply periodic boundary conditions to all atoms

atom.integrate(True, time_step) # Step 1: Update speed and location (half-step)

atom.find_force(lj) # Step 2: Calculate the force acting on each atom

atom.integrate(False, time_step) # Step 3: Update speed and location (full step)

# Output data (every Ns steps)

if step % Units.Ns == 0:

kinetic_energy = atom.find_kinetic_energy()# Calculating kinetic energy

total_Energy = kinetic_energy + atom.pe# Calculate the total energy

temperature = kinetic_energy / (1.5 * Units.kB * atom.number)# Calculate temperature from kinetic energy

# Write data to output file

outfile.write("{:16.16f} {:16.16f} {:16.16f}\n".format(temperature, kinetic_energy, atom.pe))

print(f'{step:08d} {temperature:.2f} {kinetic_energy:2f} {atom.pe:2f} {total_Energy:2f}')

outfile.close()

timer_end = time.time()

print(f"Total time: {timer_end - timer_start:.1f} s")

main()

代码

文本

第二版（numpy版本，简洁）

代码

文本

[ ]

999

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

38

39

40

41

42

43

44

45

46

47

48

49

50

51

52

53

54

55

56

57

58

59

60

61

62

63

64

65

66

67

68

69

70

71

72

73

74

75

76

77

78

79

80

81

82

83

84

85

86

87

88

89

90

91

92

93

94

95

96

97

98

99

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

›

⌄

#第二版本主程序代码（建议下载文件夹后在本地运行）

import numpy as np

import time

from line_profiler import LineProfiler

profile = LineProfiler()

def timer(func):

def func_wrapper(*args, **kwargs):

from time import time

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 Atom:

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

'''

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

:param filename: xyz文件名

:returns: None

'''

self.filename = filename

# 读取xyz文件

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

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

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

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

self.pe = 0.0

self.ke = self.getKineticEnergy()

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))

mass = np.zeros(number)

# 读取box的大小

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

box = [abc[0], abc[1], abc[2], abc[0]/2, abc[1]/2, abc[2]/2]

box = np.array(box)

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

for i in range(number):

line = f.readline().split()

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

mass[i] = float(line[4])

return number, box, coords, mass

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) -> None:

'''

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

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

:returns: None

'''

# 对于每一个维度，如果距离大于盒子的一半，就减去盒子的大小，如果距离小于盒子的一半，就加上盒子的大小

for i in range(3):

if rij[i] > self.box[i+3]:

rij[i] -= self.box[i]

elif rij[i] < -self.box[i+3]:

rij[i] += self.box[i]

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

'''

计算原子间的力和势能

:param lj: LJ势参数

:returns: None

'''

# 初始化势能和力

self.pe = 0.0

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

# 遍历所有原子对

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

for j in range(i+1, self.number):

# 计算两个原子间的距离,并应用最小镜像约定

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

self.applyMic(rij)

r2 = rij[0]**2 + rij[1]**2 + rij[2]**2

# 如果距离大于截断距离，就跳过

if r2 > lj.cutoffSquare:

continue

# 计算一些常量

r2_inv = 1.0 / r2

r4_inv = r2_inv * r2_inv

r6_inv = r2_inv * r4_inv

r8_inv = r4_inv * r4_inv

r12_inv = r4_inv * r8_inv

r14_inv = r6_inv * r8_inv

force_ij = lj.e24s6 * r8_inv - lj.e48s12 * r14_inv

# 更新势能和力

self.pe += lj.e4s12 * r12_inv - lj.e4s6 * r6_inv

self.forces[i] += force_ij * rij

self.forces[j] -= force_ij * rij

def applyPbc(self) -> None:

'''

对原子坐标应用周期性边界条件

:returns: None

'''

# 对于每一个维度，如果坐标小于0，就加上盒子大小，如果坐标大于盒子大小，就减去盒子大小

for i in range(3):

self.coords[:,i] = np.where(self.coords[:,i] < 0, self.coords[:,i] + self.box[i], self.coords[:,i])

self.coords[:,i] = np.where(self.coords[:,i] > self.box[i], self.coords[:,i] - self.box[i], self.coords[:,i])

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 readRun(filename: str='run.in') -> tuple[float, float, int]:

'''

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

:param filename: run文件名

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

'''

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])

if line.startswith('run'):

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

return velocity, time_step, run

def main():

timer_start = time.time()

# 读取run文件

velocity, time_step, run = readRun()

# 初始化原子和LJ势参数

atom = Atom()

lj = LJParameters()

# 输出热力学量的频率

thermo_file = 'thermo.out'

f = open(thermo_file, 'w')

thermo_freq = 100

# 初始化速度

atom.initializeVelocities(velocity)

# 开始模拟

for i in range(run):

atom.applyPbc()

atom.update(True, time_step)

atom.getForce(lj)

atom.update(False, time_step)

# 输出热力学量

if i % thermo_freq == 0:

ke = atom.getKineticEnergy()

temp = 2.0 * ke / (3.0 * Units.k_B * atom.number)

#print(f'{i:04d} {temp:16.16f} {atom.pe:16.16f} {ke:16.16f} {atom.pe + ke:16.16f}')

#f.write(f'{i:04d} {temp:16.16f} {atom.pe:16.16f} {ke:16.16f} {atom.pe + ke:16.16f}\n')

print(f"{temp:16.16f} {ke:16.16f} {atom.pe:16.16f}")

f.write(f"{temp:16.16f} {ke:16.16f} {atom.pe:16.16f}\n")

f.close()

timer_end = time.time()

print(f"Total time: {timer_end - timer_start:.1f} s")

if __name__ == '__main__':

main()

代码

文本