SCIF20002 – Programming and Data Analysis for Scientists– Final Assessment

AI悦创原创大约 9 分钟...约 2657 字

Deadline: 24th April 2024, 12:00pm

Email any questions to richard.stancliffe@bristol.ac.uk

Scientific problems from the gravitational dynamics of planets to the motion of atoms in a gas require us to be able to compute the behaviour of many interacting bodies. One common approach is to perform a so-called N-body simulation. In an N-body simulation we track the positions and velocities of individual particles, and use the forces between them to calculate their individual accelerations.

For this assessment, you will write an N-body code to simulate the behaviour of argon atoms confined within a box. You will then use this code to investigate how the pressure inside the box changes as you alter the system parameters. You will then write up your work, detailing how you set-up your code, demonstrating that it performs as expected for any relevant test cases, and showing your results. Your report should conclude with a discussion of any deficiencies and possible improvements that could be made.

You may use either C++ and/or Python as you see fit, but your report must justify the choice. You are welcome to use any code that you have produced during this course (including example code we have provided) but it must be properly attributed.

Part 1 – Producing the N-body code

An N-body code must follow the trajectories of a fixed number of particles, N. Each particle, i, has an associated position, $\vec{r}_i$ , and velocity, $\vec{v}_i$ , vector in 3D space. Particles interact with one another via the force between them. If the force on particle i due to particle j is given by $\vec{F}_{ij}$ , we can find the net force acting on particle i by summing over all particles (except for i!), i.e.

\vec{F}_{\text{net}} = \sum_{j \neq i} \vec{F}_{ij}.

The acceleration of particle i, $\vec{a}_i$ , then follows from Newton's second law as

\vec{a}_i = \frac{1}{m_i} \sum_{j \neq i} \vec{F}_{ij},

where $m_i$ is the mass of particle i. From the acceleration, we can then calculate the velocity and position and how they change over a given time-step, $\delta t$ . As will be clear shortly we will consider interactions between particles where the force they exert depends only on their relative positions.

To perform the integration, we use a so-called Verlet 'leap-frog' integration scheme. For each particle i, we update the velocity based on the acceleration at the current time-step using half the step size

\vec{v}_i = \vec{v}_i + \frac{\delta t}{2}a_i.

This must be done for all the particles in our simulation. We then use these new velocity to update the position across the whole time interval

\vec{r}_i = \vec{r}_i + \delta t \vec{v}_i.

Again, this must be done for all the particles. Finally, we use the new positions to recompute the acceleration $\vec{a}_i$ (this is why we had to update the positions of all the particles first!) and then calculate the new velocity after the remaining half time step:

\vec{v}_i = v_i + \frac{\delta t}{2} \vec{a}_i.

To model the interaction between two atoms in a gas we will use the Lennard-Jones potential

V(r) = 4\epsilon \left[ \left( \frac{\sigma}{r} \right)^{12} - \left( \frac{\sigma}{r} \right)^6 \right],

where $r = |\vec{r}_j - \vec{r}_p|$ is the separation between two atoms located at $\vec{r}_j$ and $\vec{r}_p$ , and $\epsilon$ and $\sigma$ are parameters that depend on the atoms. For argon, $\epsilon = 125.7 \times k_B$ ) and $\sigma = 0.3345 \times 10^{-9} m$ . Given this radially symmetric interaction potential, an atom $j$ exerts a force on atom $p$ as

\vec{F}_{ij} = -\nabla V_{LJ}(r) = \frac{24\epsilon}{\sigma} \left[ -2 \left( \frac{\sigma}{r} \right)^{13} + \left( \frac{\sigma}{r} \right)^7 \right] \vec{r}_{ji},

where

\vec{r}_{ji} = \frac{\vec{r}_j - \vec{r}_i}{|\vec{r}_j - \vec{r}_i|},

is a unit vector pointing from atom $i$ to atom $j$ .

Task: Create a piece of code that can perform the Verlet integration described above for a system of N particles. Demonstrate that the code behaves as expected for a system of two particles in close proximity (i.e. with a separation of the order of $\sigma$ ). Place one particle at the origin and the other at a point on the x-axis. Show that the two particles can remain bound to one another, and that they oscillate in place. Your report should include a brief description of your code, and plots to demonstrate the required behaviour.

It is advisable you use a dimensionless formulation of the problem using the fact that the system depends only on the parameters $m_a$ , $\sigma$ and $\epsilon$ . We can then define dimensionless units as

r = \frac{\text{length}}{\sigma}, \quad m = \frac{\text{mass}}{m_a}, \quad t = \frac{\text{time}}{\tau},

where we fix the time scale $\tau$ from $\epsilon$ as

\tau = \sqrt{\frac{m_a \sigma^2}{\epsilon}}.

Having done this all energies are measured relative to $\epsilon$ , all speeds are measured relative to $\frac{\sigma}{\tau} = \sqrt{\frac{\epsilon}{m_a}}$ , and accelerations are relative to $\frac{\sigma}{\tau^2} = \frac{\epsilon}{m_\alpha \sigma}$ . Consequently, Eq. (1) reduces to

V_{LJ}(r) = 4 \left[ \left(\frac{1}{r}\right)^{12} - \left(\frac{1}{r}\right)^{6} \right],

while Eq. (2) reduces to

\vec{F}_{pj} = 24 \left[ -2 \left( \frac{1}{r} \right)^{13} + \left( \frac{1}{r} \right)^7 \right] \vec{r}_{jp},

so the equation of motion become

\ddot{\vec{r}}_p(t) = \sum_{\substack{j=1 , j \neq p}}^{N} \vec{F}_{pj}.

Thus, in our chosen system of units there is no explicit reference to any of the parameters $m_a$ , $\sigma$ , and $\epsilon$ . What are the characteristic timescale $(\tau)$ and speed $(\sigma/\tau)$ associated with a gas of argon atoms?

You are strongly recommended to consider the efficiency of the code! Use array operations rather than for loops where possible – they are much, much faster!

Part 2 – Particles in a box

We now require some boundary conditions for the simulation. Assume our particles are going to be confined to a cubic box whose sides have a length of $L$ . We will assume that the centre of the box is centred on the origin. Consider a particle travelling only in the $x$ -direction. When it reaches the wall of the box at $x = L/2$ , it rebounds elastically – i.e. its velocity in the $x$ -direction is reversed, but its speed (and thus its kinetic energy is unchanged).

Task: Add a cubic box to your simulation so that all particles are confined within $-L/2 \leq xL/2$ , $-L/2 \leq yL/2$ and $-L/2 \leq zL/2$ . Describe how you have implemented the particle-wall collision in your report. Show that a particle hitting the walls of the box behaves as expected by including an appropriate plot.

Part 3 – Investigate!

You now have all the components you need to investigate the behaviour of the gas. The temperature of the gas is defined by the total kinetic energy of the particles as $E = \sum_{i}^{N} \frac{1}{2} m_i v_i^2 = \frac{3}{2} N k_B T,$

where $v_i$ is the speed of the $i$ th particle. The pressure $P$ is linked to the momentum flowing through a unit area in a unit time period (note you only need to consider particles passing in one direction through the area). Suppose we have an area, $A$ , lying in the $y-z$ plane. The momentum crossing this plane in time $\Delta t$ is the $x$ -component of the momentum of any particle crossing this plane, $p_x = mv_x$ where $v_x$ is the speed of the particle in the $x$ -direction. The pressure in the $x$ -direction is therefore $P_x = \frac{1}{A \Delta t} \sum_{i \text{ crossing}} m_i v_{x,i},$

where the sum is performed only over those particles $i$ crossing the plane in the time-interval $\Delta t$ . If we assume the pressure is isotropic, $P_x = P_y = P_z = \frac{1}{3}P$ where $P$ is the total pressure (you might like to test if this is true in your simulations!). Note you also only need to consider particles travelling in one direction, as pressure should be the same in both positive and negative directions when the system is in equilibrium.

Task: How does the pressure of the gas change with changes in the volume of the box and the temperature of the gas? You will need to calculate the pressure of your gas and its temperature. It is up to you to decide on how many particles $N$ to use, what time-step $\delta t$ to use, how long to run the simulations for, and what initial conditions you use. Your choices should be justified in your report based on valid scientific/computational reasons.

Hint: It is strongly recommended that you initialize the positions of your particles on a grid whose spacing is no smaller than the value of $\sigma$ . The steep dependence of the potential for small $r$ can lead to extremely large forces if you initialize the particles too close to one another. Also, the temperature and pressure estimated from your simulations will be very stochastic on the microscopic time-scales you will examine. Consider performing time-averaging to extract the expected behaviour.

Submisssion

Your submission must consist of both your (documented) code and a report into your investigation of how the gas behaves. You must include working examples for the first two tasks which we must be able to run. It should be clear how these test cases are linked to the code you use in the final investigation of the gas. Your report should include details of how you tested your code to ensure that it is working as intended (e.g. does the code conserve energy sufficiently well?), how many particles you chose to use and the length of time you chose to run your simulations for.

Please submit your report as a PDF document via the submission point on Blackboard, along with a file containing your code and everything needed to run it.

Solution

Part1

import numpy as np
import matplotlib.pyplot as plt


def Lennard_Jones_Force(r, epsilon=1.0, sigma=1.0):
    """
        计算两个粒子之间的 Lennard-Jones 力。
        参数:
        - r: 粒子之间的距离
        - epsilon: 势能深度，描述粒子间的吸引力强度
        - sigma: 粒子间的有效距离，当粒子间距离为 sigma 时，势能为 0

        返回:
        - LJ_force: Lennard-Jones 力
        """
    # 使用  Lennard-Jones 公式计算力
    LJ_force = 24 * epsilon * ((2 * (sigma / r) ** 13) - ((sigma / r) ** 7))
    return LJ_force


def Verlet_Integration(N, positions, velocities, mass, delta_t, epsilon=1.0, sigma=1.0):
    """
        执行 Verlet 积分更新粒子的位置和速度。
        参数:
        - N: 粒子数量
        - positions: 所有粒子的位置数组
        - velocities: 所有粒子的速度数组
        - mass: 粒子质量
        - delta_t: 时间步长
        - epsilon: Lennard-Jones 势的 epsilon 参数
        - sigma: Lennard-Jones 势的 sigma 参数

        返回:
        - new_positions: 更新后的位置数组
        - new_velocities: 更新后的速度数组
        """
    new_positions = np.copy(positions)
    new_velocities = np.copy(velocities)
    forces = np.zeros_like(positions)

    # 计算所有粒子之间的力
    for i in range(N):
        for j in range(i + 1, N):
            r_ij = np.linalg.norm(positions[i] - positions[j])
            force_direction = (positions[i] - positions[j]) / r_ij
            force_magnitude = Lennard_Jones_Force(r_ij, epsilon, sigma)
            force_vector = force_magnitude * force_direction
            forces[i] -= force_vector
            forces[j] += force_vector

    accelerations = forces / mass

    # 使用 Verlet 积分更新位置和速度
    for i in range(N):
        new_positions[i] = positions[i] + velocities[i] * delta_t + 0.5 * accelerations[i] * delta_t ** 2
        new_velocities[i] = velocities[i] + 0.5 * accelerations[i] * delta_t
    return new_positions, new_velocities


# demo: 模拟两个粒子
N = 2  # 粒子数
mass = 1.0  # 粒子质量
delta_t = 0.01  # 时间步长
steps = 10000  # 模拟步数
epsilon = 1.0  # Lennard-Jones 参数 epsilon
sigma = 1.0  # Lennard-Jones 参数 sigma
# 初始位置和速度
positions = np.array([[0.0, 0.0, 0.0], [1.5 * sigma, 0.0, 0.0]])
velocities = np.array([[0.0, 0.0, 0.0], [0.0, 0.0, 0.0]])

# 记录粒子位置的历史以用于绘图
positions_history = [positions.copy()]

for step in range(steps):
    positions, velocities = Verlet_Integration(N, positions, velocities, mass, delta_t, epsilon, sigma)
    positions_history.append(positions.copy())

# 绘图展示两个粒子的运动轨迹
positions_history = np.array(positions_history)
plt.plot(positions_history[:, 0, 0], positions_history[:, 0, 1], label='Particle 1')
plt.plot(positions_history[:, 1, 0], positions_history[:, 1, 1], label='Particle 2')
plt.xlabel('X Position')
plt.ylabel('Y Position')
plt.legend()
plt.show()

摘要：

首先，我定义了 Lennard- Jones 势函数来描述粒子之间的相互作用力。然后，我们实现了 Verlet 积分方法来更新粒子的位置和速度。这种方法考虑了时间步长内的加速度变化。

通过上面的代码，演示了两个粒子在空间中的运动，其中粒子之间的相互作用力由 Lennard-Jones 势给出。可以观察到，当两个粒子靠近到一定距离时，它们会因相互作用力而产生振荡运动，这符合预期的物理行为。

选择 Python 的原因是其强大的科学计算库（如 NumPy）可以简化数组和矩阵运算，而且 Python 的代码更易于阅读和维护。此外，利用 Matplotlib 库来可视化粒子运动的轨迹，进一步验证了模拟的正确性。

Part2

公众号：AI悦创【二维码】

AI悦创·编程一对一

AI悦创·推出辅导班啦，包括「Python 语言辅导班、C++ 辅导班、java 辅导班、算法/数据结构辅导班、少儿编程、pygame 游戏开发、Web、Linux」，全部都是一对一教学：一对一辅导 + 一对一答疑 + 布置作业 + 项目实践等。当然，还有线下线上摄影课程、Photoshop、Premiere 一对一教学、QQ、微信在线，随时响应！微信：Jiabcdefh

C++ 信息奥赛题解，长期更新！长期招收一对一中小学信息奥赛集训，莆田、厦门地区有机会线下上门，其他地区线上。微信：Jiabcdefh

方法一：QQ

方法二：微信：Jiabcdefh