01-Homework 01

# 导入需要的库
import numpy as np
import matplotlib.pyplot as plt

# 定义晶格类
class Lattice:
    # 初始化方法
    def __init__(self):
        # 用一个列表来存储所有的原子位置
        self.atoms = []

    # 方法用于生成并存储原子的位置
    def populate_atoms(self):
        # 通过嵌套循环确定原子的x和y位置
        for x in np.arange(-2, 3):
            for y in np.arange(-2, 3):
                # 计算并添加原子的位置到atoms列表
                self.atoms.extend([
                    [4 * x, 2 * y * np.sqrt(3)],
                    [2 + 4 * x, 2 * y * np.sqrt(3)],
                    [1 + 4 * x, np.sqrt(3) + 2 * y * np.sqrt(3)],
                    [3 + 4 * x, np.sqrt(3) + 2 * y * np.sqrt(3)]
                ])

    # 用于在图上绘制原子的方法
    def plot_atoms(self, ax):
        # 遍历每个原子的位置并在图上绘制
        for atom in self.atoms:
            ax.plot(atom[0], atom[1], "yo")

    # 用于在图上绘制圆形的方法
    def plot_circle(self, ax, radius, color="r", fill=False):
        # 创建一个圆形并添加到图上
        circle = plt.Circle(xy=(0, 0), radius=radius, color=color, fill=fill)
        ax.add_patch(circle)

    # 可视化方法
    def visualize(self):
        # 从原子位置中提取x和y的值
        x_vals, y_vals = zip(*self.atoms)
        # 计算x和y的范围
        x_dim = max(x_vals) - min(x_vals)
        y_dim = max(y_vals) - min(y_vals)
        
        # 创建一个图形
        fig, ax = plt.subplots(figsize=(9 * x_dim / y_dim, 9))
        
        # 使用上面定义的方法在图上绘制原子和圆形
        self.plot_atoms(ax)
        self.plot_circle(ax, 2, "r")
        self.plot_circle(ax, 2 * np.sqrt(3), "g")
        self.plot_circle(ax, 4, "b")
        
        # 关闭坐标轴并显示图形
        plt.axis("off")
        plt.show()

# 使用上面定义的Lattice类
lattice = Lattice()
# 生成原子位置
lattice.populate_atoms()
# 显示可视化图形
lattice.visualize()

# 导入matplotlib的pyplot模块，用于绘图
import matplotlib.pyplot as plt

# 定义Lennard-Jones势函数
def L_J_potential(r):
    # 如果距离小于等于7.0，使用以下公式计算势能
    if r <= 7.0:
        return 4 * 0.010323 * ((3.405/r)**12 - (3.405/r)**6)
    # 如果距离小于等于7.5，使用以下公式计算势能
    elif r <= 7.5:
        return -6.8102128 * 0.001 * (r-7.5)**3 - 5.5640876 * 0.001 * (r-7.5)**2
    # 其他情况下势能为0
    else:
        return 0

# 定义一个函数来计算指定距离范围内的能量并绘制图形
def compute_energy_and_plot(r_start=3, r_end=10, delta_r=0.001):
    # 使用frange函数生成一个浮点数范围的列表，并存储在r_values中
    r_values = [r for r in frange(r_start, r_end, delta_r)]
    # 根据r_values计算势能，并乘以6
    E_total = [6 * L_J_potential(r) for r in r_values]
    # 获取最小势能值
    E_min = min(E_total)
    # 获取与最小势能值对应的距离
    r_min = r_values[E_total.index(E_min)]
    
    # 绘制势能随距离变化的曲线
    plt.plot(r_values, E_total)
    plt.xlabel('distance/Angstrom')  # x轴标签
    plt.ylabel('L-J potential energy/eV')  # y轴标签
    plt.title('the first shell')  # 图形标题
    plt.show()  # 显示图形
    
    return E_min, r_min

# 定义一个生成浮点数范围的函数，作用类似于内置的range函数，但可以处理浮点数
def frange(start, stop, step):
    r = start
    while r < stop:
        yield r
        r += step

# 调用上面定义的函数来计算并绘图
E_min_1, r_min_1 = compute_energy_and_plot()
# 打印最小势能值及其对应的距离
print(E_min_1, r_min_1)

import numpy as np
import matplotlib.pyplot as plt

# Constants
A = 0.010323
SIGMA = 3.405
R_CUTOFF_1 = 7.0
R_CUTOFF_2 = 7.5
DELTA_R = 0.001
START_R = 3
END_R = 10

def L_J_potential(r):
    """Calculate L-J potential for a given distance r."""
    if r <= R_CUTOFF_1:
        return 4 * A * ( (SIGMA/r)**12 - (SIGMA/r)**6 )
    elif r <= R_CUTOFF_2:
        return -6.8102128 * 0.001 * (r-R_CUTOFF_2)**3 - 5.5640876 * 0.001 * (r-R_CUTOFF_2)**2
    else:
        return 0

def compute_energy(r):
    """Compute energy considering the nearest two shells."""
    r_2 = np.sqrt(3) * r
    return 6 * L_J_potential(r) + 6 * L_J_potential(r_2)

def find_min_energy_and_r():
    """Iterate over r values to find the minimum energy and its corresponding r value."""
    r_values = np.arange(START_R, END_R, DELTA_R)
    energies = [compute_energy(r) for r in r_values]

    r_min = r_values[np.argmin(energies)]
    E_min = min(energies)
    return r_min, E_min, r_values, energies

def plot_energy_vs_distance(r_values, energies):
    """Plot L-J potential energy against distance."""
    plt.plot(r_values, energies)
    plt.xlabel('distance/Angstrom')
    plt.ylabel('L-J potential energy/eV')
    plt.title('the nearest two shells')
    plt.show()

def main():
    r_min, E_min, r_values, energies = find_min_energy_and_r()
    plot_energy_vs_distance(r_values, energies)
    print(E_min, r_min)

if __name__ == "__main__":
    main()

import numpy as np                 # 引入NumPy库，用于数值计算
import matplotlib.pyplot as plt   # 引入matplotlib库，用于绘图

# Constants
# 定义常量
A = 0.010323                      # Lennard-Jones势能的A值
SIGMA = 3.405                     # Lennard-Jones势能的σ值
R_CUTOFF_1 = 7.0                  # 第一个距离截断值
R_CUTOFF_2 = 7.5                  # 第二个距离截断值
DELTA_R = 0.001                  # r值的增量
START_R = 3                      # 开始的r值
END_R = 10                       # 结束的r值

def L_J_potential(r):
    """Calculate L-J potential for a given distance r."""
    # 根据给定的距离r计算L-J势能
    if r <= R_CUTOFF_1:
        # 如果r小于等于第一个截断值
        return 4 * A * ( (SIGMA/r)**12 - (SIGMA/r)**6 )
    elif r <= R_CUTOFF_2:
        # 如果r在两个截断值之间
        return -6.8102128 * 0.001 * (r-R_CUTOFF_2)**3 - 5.5640876 * 0.001 * (r-R_CUTOFF_2)**2
    else:
        # 如果r大于第二个截断值
        return 0

def compute_energy(r):
    """Compute energy considering the nearest two shells."""
    # 计算考虑最近的两个shells的能量
    r_2 = np.sqrt(3) * r            # 计算第二个shell的r值
    return 6 * L_J_potential(r) + 6 * L_J_potential(r_2)

def find_min_energy_and_r():
    """Iterate over r values to find the minimum energy and its corresponding r value."""
    # 遍历r值以找到最小的能量及其对应的r值
    r_values = np.arange(START_R, END_R, DELTA_R)     # 创建一个r值的数组
    energies = [compute_energy(r) for r in r_values]  # 计算每个r值的能量

    r_min = r_values[np.argmin(energies)]             # 找到能量最小值对应的r值
    E_min = min(energies)                             # 找到最小的能量值
    return r_min, E_min, r_values, energies

def plot_energy_vs_distance(r_values, energies):
    """Plot L-J potential energy against distance."""
    # 绘制L-J势能与距离的关系图
    plt.plot(r_values, energies)                      # 绘制能量与距离的图
    plt.xlabel('distance/Angstrom')                   # 设置x轴标签
    plt.ylabel('L-J potential energy/eV')             # 设置y轴标签
    plt.title('the nearest two shells')               # 设置图表标题
    plt.show()                                       # 显示图表

def main():
    # 主函数
    r_min, E_min, r_values, energies = find_min_energy_and_r()  # 获取最小能量、对应的r值和所有r值与能量的列表
    plot_energy_vs_distance(r_values, energies)                 # 绘制L-J势能与距离的图表
    print(E_min, r_min)                                         # 打印最小的能量值及其对应的r值

if __name__ == "__main__":
    # 如果此脚本作为主程序运行，则执行main()函数
    main()