L-J势场的Python代码实现

import matplotlib.pyplot as plt  import numpy as npdef lj_potential(r: np.ndarray, cutoff: float) -> np.ndarray:    """    计算 Lennard-Jones 势能    参数:        r: 距离数组        cutoff: 截断半径    返回:        势能数组    """    # 避免除零错误    r = np.maximum(r, 1e-8)    potential = 4 * (1 / r**12 - 1 / r**6)        potential[r > cutoff] = 0.0    return potentialdef lj_force(r: np.ndarray, cutoff: float) -> np.ndarray:      """    计算 Lennard-Jones 作用力 F = -dU/dr    """    r = np.maximum(r, 1e-8)      force = -4 * (-12 / r**13 + 6 / r**7)    force[r > cutoff] = 0.0      return force# ===================== 可配置参数 =====================CUTOFF_RADIUS = 4.0R_MIN, R_MAX = 1.0, 4.0NUM_POINTS = 500# ======================================================# 生成距离数组r_array = np.linspace(R_MIN, R_MAX, NUM_POINTS)# 计算势能与力u_vals = lj_potential(r_array, CUTOFF_RADIUS)    f_vals = lj_force(r_array, CUTOFF_RADIUS)# 创建画布与轴fig, ax_energy = plt.subplots(figsize=(6, 4))  # --------------------- 绘制势能（左轴）---------------------color_energy = "tab:blue"  ax_energy.plot(r_array, u_vals, color=color_energy, label="Potential")ax_energy.axhline(0, ls=":", c="k", alpha=0.5, zorder=-1)ax_energy.set_xlabel(r"Distance ($\sigma$)")  ax_energy.set_ylabel(r"Energy ($\epsilon$)", color=color_energy)    ax_energy.tick_params(axis="y", labelcolor=color_energy)ax_energy.set_ylim(-1.2, 1.2)ax_energy.grid(True, alpha=0.3)# --------------------- 绘制力（右轴）---------------------ax_force = ax_energy.twinx()    color_force = "tab:orange"ax_force.plot(r_array, f_vals, "--", color=color_force, label="Force")ax_force.set_ylabel(r"Force ($\epsilon/\sigma$)", color=color_force)    ax_force.tick_params(axis="y", labelcolor=color_force)ax_force.set_ylim(-25, 25)# 自动调整布局plt.tight_layout()  # 保存图片plt.savefig("fig-lj.pdf", dpi=300, bbox_inches="tight")      plt.close()

import numpy as npimport matplotlib.pyplot as plt  def U(r, rc):    v = 4 * (1/r**12 - 1/r**6)      v[r>rc] = 0    return v  def F(r, rc):    v = -4 * (-12/r**13 + 6/r**7)          v[r>rc] = 0    return v  rc = 4  r = np.linspace(1,4,500)    fig, ax1 = plt.subplots()  ax1.plot(r, U(r,rc))ax2 = ax1.twinx()  ax2.plot(r, F(r,rc), '--')    plt.savefig("fig.pdf")

 jupyter nbconvert --to html .\6.2.1L-J.ipynb