多场耦合优化-主题088-纳米技术与分子工程
·
第八十八篇：纳米技术与分子工程
摘要
纳米技术与分子工程是现代科学与工程的前沿交叉领域，涉及从原子尺度到宏观尺度的多尺度建模与仿真。本教程系统介绍纳米技术的仿真方法，包括分子动力学模拟、密度泛函理论、连续介质模型以及多尺度耦合方法。通过六个典型案例——碳纳米管力学性能仿真、纳米颗粒自组装模拟、分子机器运动仿真、石墨烯电子性质计算、纳米流体传热仿真和DNA分子动力学模拟，展示纳米尺度下的物理化学现象和工程应用。本教程将帮助读者掌握纳米尺度仿真的核心方法，理解从分子到宏观的跨尺度传递机制，为纳米材料设计和分子工程应用奠定基础。
关键词
纳米技术, 分子动力学, 密度泛函理论, 自组装, 分子机器, 多尺度建模, 纳米材料, 分子工程
1. 纳米技术基础理论
1.1 纳米尺度特征
纳米技术研究的尺度范围通常在1-100纳米之间，这一尺度具有独特的物理化学性质：
表面效应：纳米材料的比表面积显著增大，表面原子比例急剧上升。对于球形纳米颗粒，表面原子比例可近似表示为：
fsurface=4πr2⋅d43πr3⋅ρ=3drρf_{surface} = \frac{4\pi r^2 \cdot d}{\frac{4}{3}\pi r^3 \cdot \rho} = \frac{3d}{r\rho}fsurface=34πr3⋅ρ4πr2⋅d=rρ3d
其中ddd为原子直径，rrr为颗粒半径，ρ\rhoρ为原子数密度。当颗粒尺寸减小到10纳米以下时，表面原子可占总原子数的50%以上。
量子尺寸效应：当材料尺寸减小到与电子德布罗意波长或激子玻尔半径相当时，电子能级由连续变为离散，导致光学、电学和磁学性质发生显著变化。量子限域效应可用粒子在势阱中的能级公式描述：
En=n2h28m∗L2E_n = \frac{n^2 h^2}{8m^* L^2}En=8m∗L2n2h2
其中nnn为量子数，hhh为普朗克常数，m∗m^*m∗为有效质量，LLL为限域尺寸。
宏观量子隧道效应：微观粒子具有贯穿势垒的能力，在纳米尺度下这一效应更加显著，影响电子输运和磁化反转等过程。
1.2 纳米材料分类
纳米材料按维度可分为四类：
零维纳米材料：量子点、纳米颗粒、富勒烯。在所有三个维度上都处于纳米尺度，具有离散的电子态。
一维纳米材料：纳米线、纳米管、纳米棒。在两个维度上受限，电子在一个方向上自由运动，表现出独特的量子输运性质。
二维纳米材料：石墨烯、纳米薄膜、二维材料。在一个维度上受限，电子在平面内自由运动，具有优异的电子和力学性能。
三维纳米材料：纳米多孔材料、纳米复合材料。由纳米结构单元组装而成的宏观材料，兼具纳米效应和宏观可加工性。
1.3 分子工程原理
分子工程是通过设计和操控分子结构来实现特定功能的科学。核心原理包括：
自下而上组装：从原子、分子出发，通过化学合成或自组装构建复杂结构。这种方法可以精确控制材料的组成和结构。
分子识别与自组装：利用分子间的非共价相互作用（氢键、范德华力、静电作用、疏水作用）实现分子的选择性识别和有序组装。
分子机器：能够在外界刺激下执行机械运动的分子器件，包括分子马达、分子开关、分子电梯等。2016年诺贝尔化学奖授予了分子机器的设计与合成研究。
2. 分子动力学模拟方法
2.1 基本原理
分子动力学（Molecular Dynamics, MD）模拟通过数值求解牛顿运动方程来追踪系统中所有原子的运动轨迹。对于第iii个原子，运动方程为：
mid2ridt2=Fi=−∇iU(r1,r2,...,rN)m_i \frac{d^2 \mathbf{r}_i}{dt^2} = \mathbf{F}_i = -\nabla_i U(\mathbf{r}_1, \mathbf{r}_2, ..., \mathbf{r}_N)midt2d2ri=Fi=−∇iU(r1,r2,...,rN)
其中mim_imi为原子质量，ri\mathbf{r}_iri为位置矢量，Fi\mathbf{F}_iFi为作用力，UUU为系统势能函数。
2.2 势能函数
势能函数（力场）是分子动力学模拟的核心，决定了计算的精度和效率。常用的势能函数包括：
Lennard-Jones势：描述中性原子或分子间的范德华相互作用：
ULJ(r)=4ε[(σr)12−(σr)6]U_{LJ}(r) = 4\varepsilon \left[\left(\frac{\sigma}{r}\right)^{12} - \left(\frac{\sigma}{r}\right)^6\right]ULJ(r)=4ε[(rσ)12−(rσ)6]
其中ε\varepsilonε为势阱深度，σ\sigmaσ为原子直径。r−12r^{-12}r−12项描述短程排斥，r−6r^{-6}r−6项描述长程吸引。
键合势：描述化学键的伸缩、弯曲和扭转：
Ubond=12kb(r−r0)2U_{bond} = \frac{1}{2}k_b(r - r_0)^2Ubond=21kb(r−r0)2
Uangle=12kθ(θ−θ0)2U_{angle} = \frac{1}{2}k_\theta(\theta - \theta_0)^2Uangle=21kθ(θ−θ0)2
Utorsion=∑n=13Vn[1+cos(nϕ−ϕn)]U_{torsion} = \sum_{n=1}^{3} V_n[1 + \cos(n\phi - \phi_n)]Utorsion=n=1∑3Vn[1+cos(nϕ−ϕn)]
反应力场（ReaxFF）：能够描述化学键的断裂和形成，适用于化学反应模拟。
2.3 积分算法
常用的数值积分算法包括：
Verlet算法：
r(t+Δt)=2r(t)−r(t−Δt)+F(t)mΔt2\mathbf{r}(t + \Delta t) = 2\mathbf{r}(t) - \mathbf{r}(t - \Delta t) + \frac{\mathbf{F}(t)}{m}\Delta t^2r(t+Δt)=2r(t)−r(t−Δt)+mF(t)Δt2
速度Verlet算法：
r(t+Δt)=r(t)+v(t)Δt+F(t)2mΔt2\mathbf{r}(t + \Delta t) = \mathbf{r}(t) + \mathbf{v}(t)\Delta t + \frac{\mathbf{F}(t)}{2m}\Delta t^2r(t+Δt)=r(t)+v(t)Δt+2mF(t)Δt2
v(t+Δt)=v(t)+F(t)+F(t+Δt)2mΔt\mathbf{v}(t + \Delta t) = \mathbf{v}(t) + \frac{\mathbf{F}(t) + \mathbf{F}(t + \Delta t)}{2m}\Delta tv(t+Δt)=v(t)+2mF(t)+F(t+Δt)Δt
速度Verlet算法具有时间可逆性和较好的能量守恒特性，是目前最常用的算法。
2.4 系综控制
分子动力学模拟可以在不同统计系综下进行：
NVE系综：微正则系综，粒子数NNN、体积VVV、能量EEE恒定。适用于研究系统的本征动力学行为。
NVT系综：正则系综，粒子数NNN、体积VVV、温度TTT恒定。通过Nose-Hoover热浴或Berendsen热浴控制温度。
NPT系综：等温等压系综，粒子数NNN、压强PPP、温度TTT恒定。通过Parrinello-Rahman方法或Berendsen压强耦合控制压强。
3. 密度泛函理论
3.1 Hohenberg-Kohn定理
密度泛函理论（Density Functional Theory, DFT）基于两个基本定理：
第一定理：外势v(r)v(\mathbf{r})v(r)（从而哈密顿量）由基态电子密度ρ(r)\rho(\mathbf{r})ρ(r)唯一确定（至多相差一个常数）。
第二定理：对于给定的外势，基态能量可以通过能量泛函对密度的变分极小化得到：
E[ρ]=FHK[ρ]+∫v(r)ρ(r)drE[\rho] = F_{HK}[\rho] + \int v(\mathbf{r})\rho(\mathbf{r})d\mathbf{r}E[ρ]=FHK[ρ]+∫v(r)ρ(r)dr
其中FHK[ρ]F_{HK}[\rho]FHK[ρ]是普适的密度泛函，包含动能和电子相互作用能。
3.2 Kohn-Sham方程
Kohn和Sham将多电子问题转化为等效的单电子问题。Kohn-Sham方程为：
[−ℏ22m∇2+veff(r)]ϕi(r)=εiϕi(r)\left[-\frac{\hbar^2}{2m}\nabla^2 + v_{eff}(\mathbf{r})\right]\phi_i(\mathbf{r}) = \varepsilon_i\phi_i(\mathbf{r})[−2mℏ2∇2+veff(r)]ϕi(r)=εiϕi(r)
其中有效势包含外势、Hartree势和交换关联势：
veff(r)=vext(r)+∫ρ(r′)∣r−r′∣dr′+vxc(r)v_{eff}(\mathbf{r}) = v_{ext}(\mathbf{r}) + \int \frac{\rho(\mathbf{r}')}{|\mathbf{r} - \mathbf{r}'|}d\mathbf{r}' + v_{xc}(\mathbf{r})veff(r)=vext(r)+∫∣r−r′∣ρ(r′)dr′+vxc(r)
电子密度由Kohn-Sham轨道构造：
ρ(r)=∑i=1N∣ϕi(r)∣2\rho(\mathbf{r}) = \sum_{i=1}^{N} |\phi_i(\mathbf{r})|^2ρ(r)=i=1∑N∣ϕi(r)∣2
3.3 交换关联泛函
交换关联泛函的近似是DFT计算精度的关键：
LDA（局域密度近似）：假设电子密度在空间缓慢变化，使用均匀电子气的交换关联能。
GGA（广义梯度近似）：考虑电子密度的梯度修正，如PBE、BLYP等泛函。
杂化泛函：混合Hartree-Fock交换和DFT交换关联，如B3LYP、HSE06等，通常能给出更精确的结果。
4. 纳米尺度传热与输运
4.1 声子输运
在纳米尺度，声子（晶格振动的量子）的输运行为与宏观尺度显著不同。声子输运可以用玻尔兹曼输运方程描述：
∂f∂t+vg⋅∇f=(∂f∂t)coll\frac{\partial f}{\partial t} + \mathbf{v}_g \cdot \nabla f = \left(\frac{\partial f}{\partial t}\right)_{coll}∂t∂f+vg⋅∇f=(∂t∂f)coll
其中fff为声子分布函数，vg\mathbf{v}_gvg为群速度。
声子平均自由程：声子在两次散射之间传播的平均距离。当结构尺寸小于声子平均自由程时，出现弹道输运，热导率显著降低。
界面热阻（Kapitza热阻）：由于声子在界面处的反射和透射不完全，导致界面处存在温度跃变：
RK=ΔTQ/AR_K = \frac{\Delta T}{Q/A}RK=Q/AΔT
其中ΔT\Delta TΔT为温度跃变，QQQ为热流，AAA为界面面积。
4.2 电子输运
纳米结构的电子输运表现出量子特性：
弹道输运：当通道长度小于电子平均自由程时，电子不发生散射，电导量子化为：
G=2e2h∑nTnG = \frac{2e^2}{h} \sum_{n} T_nG=h2e2n∑Tn
其中TnT_nTn为第nnn个模式的透射系数。
库仑阻塞：在量子点等小型导体中，单个电子的充电能EC=e2/2CE_C = e^2/2CEC=e2/2C可能远大于热能kBTk_BTkBT，导致电子输运呈现离散台阶（库仑台阶）。
5. 典型案例分析
案例1：碳纳米管力学性能仿真
碳纳米管是由石墨烯片层卷曲而成的管状结构，具有优异的力学性能。本案例通过分子动力学方法模拟单壁碳纳米管的拉伸行为。
物理模型：
- 采用(10,10)扶手椅型单壁碳纳米管
- 使用AIREBO势函数描述碳-碳相互作用
- 应用周期性边界条件
- 在恒温恒容（NVT）系综下进行模拟
仿真步骤：
- 构建碳纳米管初始结构
- 能量最小化弛豫
- 一端固定，另一端施加拉伸位移
- 记录应力和应变，计算杨氏模量和断裂强度
Python实现：
# -*- coding: utf-8 -*-
"""
案例1：碳纳米管力学性能仿真
使用简化分子动力学模型模拟碳纳米管拉伸
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyArrowPatch
import warnings
warnings.filterwarnings('ignore')
import os
# 创建输出目录
output_dir = r'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题088\\output'
os.makedirs(output_dir, exist_ok=True)
# 设置字体
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("案例1：碳纳米管力学性能仿真")
print("=" * 60)
class CarbonNanotubeMD:
"""碳纳米管分子动力学模拟器"""
def __init__(self, n=10, m=10, length=10):
"""
初始化碳纳米管
n, m: 手性指数 (n,n)为扶手椅型
length: 管长 (nm)
"""
self.n = n
self.m = m
self.length = length
# 碳纳米管几何参数
self.a_cc = 0.142 # C-C键长 (nm)
self.diameter = self.a_cc * np.sqrt(n**2 + m**2 + n*m) / np.pi * np.sqrt(3)
self.radius = self.diameter / 2
# 力场参数 (简化Lennard-Jones + 谐振键)
self.epsilon = 0.00284 # eV
self.sigma = 0.34 # nm
self.k_bond = 469.0 # eV/nm^2 (键伸缩刚度)
self.k_angle = 7.0 # eV/rad^2 (键角弯曲刚度)
self.r0 = self.a_cc # 平衡键长
# 原子质量
self.m_C = 12.011 # amu
# 初始化原子位置
self.build_structure()
def build_structure(self):
"""构建扶手椅型碳纳米管原子结构"""
# 每个单胞有4个碳原子
num_unit_cells = int(self.length / (self.a_cc * np.sqrt(3)))
atoms_per_cell = 4
self.N_atoms = num_unit_cells * atoms_per_cell
# 分配位置数组
self.positions = np.zeros((self.N_atoms, 3))
self.velocities = np.zeros((self.N_atoms, 3))
self.forces = np.zeros((self.N_atoms, 3))
self.masses = np.ones(self.N_atoms) * self.m_C
# 构建原子位置
idx = 0
for i in range(num_unit_cells):
z = i * self.a_cc * np.sqrt(3)
# 单胞内4个原子的位置
angles = np.array([0, np.pi/3, np.pi, 4*np.pi/3])
z_offsets = np.array([0, self.a_cc, self.a_cc*np.sqrt(3)/2, self.a_cc*(1+np.sqrt(3)/2)])
for j in range(4):
if idx < self.N_atoms:
theta = angles[j % 4]
self.positions[idx, 0] = self.radius * np.cos(theta)
self.positions[idx, 1] = self.radius * np.sin(theta)
self.positions[idx, 2] = z + z_offsets[j % 4]
idx += 1
# 构建近邻列表（简化：距离小于1.5倍键长）
self.build_neighbor_list()
def build_neighbor_list(self):
"""构建原子近邻列表"""
self.neighbors = []
cutoff = 1.5 * self.a_cc
for i in range(self.N_atoms):
neighbors_i = []
for j in range(self.N_atoms):
if i != j:
r_ij = np.linalg.norm(self.positions[i] - self.positions[j])
if r_ij < cutoff:
neighbors_i.append(j)
self.neighbors.append(neighbors_i)
def compute_forces(self):
"""计算原子受力"""
self.forces.fill(0)
potential_energy = 0
# 键伸缩力 (谐振子模型)
for i in range(self.N_atoms):
for j in self.neighbors[i]:
if j > i: # 避免重复计算
r_vec = self.positions[j] - self.positions[i]
r = np.linalg.norm(r_vec)
if r > 0:
# 谐振键势能
dr = r - self.r0
f_mag = -self.k_bond * dr
f_vec = f_mag * r_vec / r
self.forces[i] -= f_vec
self.forces[j] += f_vec
potential_energy += 0.5 * self.k_bond * dr**2
# 非键相互作用 (Lennard-Jones)
for i in range(self.N_atoms):
for j in range(i+1, self.N_atoms):
r_vec = self.positions[j] - self.positions[i]
r = np.linalg.norm(r_vec)
if r < 2.5 * self.sigma and r > 0.15: # 排除键合原子
# LJ势
sr6 = (self.sigma / r)**6
sr12 = sr6**2
f_mag = 24 * self.epsilon * (2*sr12 - sr6) / r
f_vec = f_mag * r_vec / r
self.forces[i] -= f_vec
self.forces[j] += f_vec
potential_energy += 4 * self.epsilon * (sr12 - sr6)
return potential_energy
def velocity_verlet(self, dt):
"""速度Verlet积分"""
# 更新位置
self.positions += self.velocities * dt + 0.5 * self.forces / self.masses[:, np.newaxis] * dt**2
# 保存旧力
forces_old = self.forces.copy()
# 计算新力
pe = self.compute_forces()
# 更新速度
self.velocities += 0.5 * (forces_old + self.forces) / self.masses[:, np.newaxis] * dt
# 计算动能
ke = 0.5 * np.sum(self.masses * np.sum(self.velocities**2, axis=1))
return pe, ke
def berendsen_thermostat(self, target_temp, tau, dt):
"""Berendsen热浴控温"""
kinetic_energy = 0.5 * np.sum(self.masses * np.sum(self.velocities**2, axis=1))
current_temp = 2 * kinetic_energy / (3 * self.N_atoms * 8.617e-5) # eV to K
if current_temp > 0:
scale = np.sqrt(1 + (dt/tau) * (target_temp/current_temp - 1))
self.velocities *= scale
def apply_strain(self, strain_rate, dt):
"""施加拉伸应变"""
# 固定底部原子，拉伸顶部原子
z_max = np.max(self.positions[:, 2])
z_min = np.min(self.positions[:, 2])
for i in range(self.N_atoms):
z = self.positions[i, 2]
# 顶部10%原子受拉伸
if z > z_max - 0.5:
self.positions[i, 2] += strain_rate * dt * self.length
self.velocities[i, 2] = strain_rate * self.length
# 底部10%原子固定
elif z < z_min + 0.5:
self.velocities[i, :] = 0
def compute_stress(self):
"""计算应力"""
# 简化应力计算：轴向力除以截面积
z_max = np.max(self.positions[:, 2])
z_min = np.min(self.positions[:, 2])
# 计算顶部原子受到的总力
top_force = 0
for i in range(self.N_atoms):
if self.positions[i, 2] > z_max - 0.5:
top_force += self.forces[i, 2]
# 截面积 (环形)
thickness = 0.34 # nm (石墨烯层间距)
area = np.pi * ((self.radius + thickness/2)**2 - (self.radius - thickness/2)**2)
# 应力 (GPa)
stress = abs(top_force) / area * 160.2 # eV/nm^3 to GPa
return stress
# 运行拉伸模拟
print("\n初始化碳纳米管结构...")
nanotube = CarbonNanotubeMD(n=10, m=10, length=10)
print(f"原子数: {nanotube.N_atoms}")
print(f"直径: {nanotube.diameter:.3f} nm")
print(f"长度: {nanotube.length:.1f} nm")
# 能量最小化
print("\n能量最小化...")
dt = 0.001 # ps
temp = 300 # K
for step in range(1000):
pe, ke = nanotube.velocity_verlet(dt)
nanotube.berendsen_thermostat(temp, 0.1, dt)
if step % 200 == 0:
print(f" Step {step}: PE = {pe:.2f} eV, KE = {ke:.2f} eV")
# 拉伸模拟
print("\n开始拉伸模拟...")
strain_rate = 0.001 # 1/ps
max_strain = 0.3
num_steps = int(max_strain / strain_rate / dt)
# 记录数据
strains = []
stresses = []
positions_history = []
sample_interval = max(1, num_steps // 100)
for step in range(num_steps):
pe, ke = nanotube.velocity_verlet(dt)
nanotube.apply_strain(strain_rate, dt)
nanotube.berendsen_thermostat(temp, 0.1, dt)
if step % sample_interval == 0:
strain = step * strain_rate * dt
stress = nanotube.compute_stress()
strains.append(strain)
stresses.append(stress)
if step % (10 * sample_interval) == 0:
positions_history.append(nanotube.positions.copy())
if step % 5000 == 0:
print(f" Step {step}/{num_steps}, Strain = {step*strain_rate*dt:.3f}, Stress = {stress:.2f} GPa")
strains = np.array(strains)
stresses = np.array(stresses)
print(f"\n最大应力: {np.max(stresses):.2f} GPa")
print(f"杨氏模量: {stresses[10]/strains[10]:.2f} GPa")
# 可视化结果
fig = plt.figure(figsize=(16, 12))
# 1. 应力-应变曲线
ax1 = fig.add_subplot(2, 3, 1)
ax1.plot(strains * 100, stresses, 'b-', linewidth=2)
ax1.set_xlabel('Strain (%)', fontsize=11)
ax1.set_ylabel('Stress (GPa)', fontsize=11)
ax1.set_title('Stress-Strain Curve', fontsize=12, fontweight='bold')
ax1.grid(True, alpha=0.3)
ax1.axhline(y=0, color='k', linewidth=0.5)
# 2. 碳纳米管结构示意图
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
# 绘制初始结构
pos_init = positions_history[0] if positions_history else nanotube.positions
ax2.scatter(pos_init[:, 0], pos_init[:, 1], pos_init[:, 2],
c='gray', s=50, alpha=0.8, label='Initial')
ax2.set_xlabel('X (nm)')
ax2.set_ylabel('Y (nm)')
ax2.set_zlabel('Z (nm)')
ax2.set_title('Initial Structure', fontsize=12, fontweight='bold')
# 3. 变形后的结构
ax3 = fig.add_subplot(2, 3, 3, projection='3d')
if len(positions_history) > 5:
pos_final = positions_history[-1]
ax3.scatter(pos_final[:, 0], pos_final[:, 1], pos_final[:, 2],
c='red', s=50, alpha=0.8, label='Deformed')
ax3.set_xlabel('X (nm)')
ax3.set_ylabel('Y (nm)')
ax3.set_zlabel('Z (nm)')
ax3.set_title('Deformed Structure', fontsize=12, fontweight='bold')
# 4. 能量演化
ax4 = fig.add_subplot(2, 3, 4)
# 重新计算能量
pe_list = []
ke_list = []
total_energy = []
nanotube_temp = CarbonNanotubeMD(n=10, m=10, length=10)
for step in range(1000):
pe, ke = nanotube_temp.velocity_verlet(dt)
nanotube_temp.berendsen_thermostat(temp, 0.1, dt)
pe_list.append(pe)
ke_list.append(ke)
total_energy.append(pe + ke)
steps = np.arange(len(pe_list))
ax4.plot(steps, pe_list, 'b-', label='Potential', linewidth=1.5)
ax4.plot(steps, ke_list, 'r-', label='Kinetic', linewidth=1.5)
ax4.plot(steps, total_energy, 'g--', label='Total', linewidth=1.5)
ax4.set_xlabel('Step', fontsize=11)
ax4.set_ylabel('Energy (eV)', fontsize=11)
ax4.set_title('Energy Minimization', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)
# 5. 杨氏模量计算
ax5 = fig.add_subplot(2, 3, 5)
# 线性拟合计算杨氏模量
linear_region = strains < 0.05
if np.sum(linear_region) > 5:
E_fit = np.polyfit(strains[linear_region], stresses[linear_region], 1)
E_young = E_fit[0]
ax5.plot(strains[linear_region] * 100, stresses[linear_region], 'bo', markersize=4, label='Data')
ax5.plot(strains[linear_region] * 100,
np.polyval(E_fit, strains[linear_region]),
'r-', linewidth=2, label=f'E = {E_young:.1f} GPa')
ax5.set_xlabel('Strain (%)', fontsize=11)
ax5.set_ylabel('Stress (GPa)', fontsize=11)
ax5.set_title('Young\'s Modulus', fontsize=12, fontweight='bold')
ax5.legend()
ax5.grid(True, alpha=0.3)
# 6. 径向分布函数
ax6 = fig.add_subplot(2, 3, 6)
# 计算RDF
pos = nanotube.positions
rdf_r = np.linspace(0.1, 1.0, 100)
rdf_g = np.zeros_like(rdf_r)
dr = rdf_r[1] - rdf_r[0]
for i in range(nanotube.N_atoms):
for j in range(i+1, nanotube.N_atoms):
r_ij = np.linalg.norm(pos[i] - pos[j])
idx = int((r_ij - rdf_r[0]) / dr)
if 0 <= idx < len(rdf_r):
rdf_g[idx] += 2 # 每对原子贡献2
# 归一化
shell_volume = 4 * np.pi * rdf_r**2 * dr
rho = nanotube.N_atoms / (np.pi * nanotube.radius**2 * nanotube.length)
rdf_g = rdf_g / (nanotube.N_atoms * rho * shell_volume)
ax6.plot(rdf_r, rdf_g, 'b-', linewidth=2)
ax6.set_xlabel('r (nm)', fontsize=11)
ax6.set_ylabel('g(r)', fontsize=11)
ax6.set_title('Radial Distribution Function', fontsize=12, fontweight='bold')
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/case1_nanotube_mechanics.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例1完成：碳纳米管力学性能仿真")
print(f" 结果保存在: {output_dir}/case1_nanotube_mechanics.png")
案例2：纳米颗粒自组装模拟
纳米颗粒的自组装是构建有序纳米结构的重要方法。本案例模拟疏水性纳米颗粒在溶液中的自组装过程。
物理模型：
- 纳米颗粒用粗粒化珠子表示
- 颗粒间通过Lennard-Jones势相互作用
- 疏水相互作用通过调节LJ参数实现
- 溶剂效应隐式处理
仿真目标：
- 观察纳米颗粒从无序到有序的相变
- 分析团簇尺寸分布
- 研究自组装动力学
Python实现：
# -*- coding: utf-8 -*-
"""
案例2：纳米颗粒自组装模拟
模拟疏水性纳米颗粒的聚集和自组装过程
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
from matplotlib.collections import PatchCollection
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题088\\output'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("\n案例2：纳米颗粒自组装模拟")
print("=" * 60)
class NanoparticleSelfAssembly:
"""纳米颗粒自组装模拟器"""
def __init__(self, N_particles=200, box_size=50.0, temperature=300):
"""
初始化系统
N_particles: 颗粒数量
box_size: 模拟盒子尺寸 (nm)
temperature: 温度 (K)
"""
self.N = N_particles
self.L = box_size
self.T = temperature
# 颗粒参数
self.sigma = 2.0 # 颗粒直径 (nm)
self.epsilon = 0.5 # 相互作用强度 (eV)
self.mass = 100.0 # 颗粒质量 (amu)
# 初始化位置 (随机分布)
self.positions = np.random.rand(self.N, 3) * self.L
self.velocities = np.random.randn(self.N, 3) * np.sqrt(self.T / self.mass)
self.forces = np.zeros((self.N, 3))
# 记录团簇信息
self.clusters_history = []
self.potential_energy_history = []
def minimum_image(self, r):
"""最小镜像约定处理周期性边界"""
return r - self.L * np.round(r / self.L)
def compute_forces(self):
"""计算Lennard-Jones力"""
self.forces.fill(0)
potential = 0
for i in range(self.N):
for j in range(i+1, self.N):
r_vec = self.positions[j] - self.positions[i]
r_vec = self.minimum_image(r_vec)
r = np.linalg.norm(r_vec)
if r < 2.5 * self.sigma and r > 0:
# LJ势和力
sr = self.sigma / r
sr6 = sr**6
sr12 = sr6**2
f_mag = 24 * self.epsilon * (2*sr12 - sr6) / r
f_vec = f_mag * r_vec / r
self.forces[i] -= f_vec
self.forces[j] += f_vec
potential += 4 * self.epsilon * (sr12 - sr6)
return potential
def velocity_verlet(self, dt):
"""速度Verlet积分"""
self.positions += self.velocities * dt + 0.5 * self.forces / self.mass * dt**2
# 周期性边界条件
self.positions = self.positions % self.L
forces_old = self.forces.copy()
pe = self.compute_forces()
self.velocities += 0.5 * (forces_old + self.forces) / self.mass * dt
ke = 0.5 * self.mass * np.sum(self.velocities**2)
return pe, ke
def andersen_thermostat(self, dt, nu=0.1):
"""Andersen热浴"""
for i in range(self.N):
if np.random.rand() < nu * dt:
self.velocities[i] = np.random.randn(3) * np.sqrt(self.T / self.mass)
def find_clusters(self, cutoff=None):
"""识别团簇 (基于距离)"""
if cutoff is None:
cutoff = 1.5 * self.sigma
clusters = []
visited = np.zeros(self.N, dtype=bool)
for i in range(self.N):
if not visited[i]:
cluster = []
stack = [i]
while stack:
j = stack.pop()
if not visited[j]:
visited[j] = True
cluster.append(j)
# 寻找近邻
for k in range(self.N):
if not visited[k]:
r_vec = self.positions[k] - self.positions[j]
r_vec = self.minimum_image(r_vec)
if np.linalg.norm(r_vec) < cutoff:
stack.append(k)
clusters.append(cluster)
return clusters
def run_simulation(self, n_steps=10000, dt=0.001):
"""运行模拟"""
print(f"\n运行自组装模拟...")
print(f"颗粒数: {self.N}")
print(f"盒子尺寸: {self.L} nm")
print(f"温度: {self.T} K")
# 平衡阶段
print("\n平衡阶段...")
for step in range(2000):
pe, ke = self.velocity_verlet(dt)
self.andersen_thermostat(dt)
# 生产阶段
print("生产阶段...")
sample_interval = 100
for step in range(n_steps):
pe, ke = self.velocity_verlet(dt)
self.andersen_thermostat(dt)
if step % sample_interval == 0:
clusters = self.find_clusters()
self.clusters_history.append([len(c) for c in clusters])
self.potential_energy_history.append(pe)
if step % 2000 == 0:
print(f" Step {step}/{n_steps}, PE = {pe:.2f} eV, Clusters = {len(clusters)}")
print(f"\n模拟完成！")
print(f"最终团簇数: {len(clusters)}")
return clusters
# 运行模拟
assembly = NanoparticleSelfAssembly(N_particles=200, box_size=50.0, temperature=300)
final_clusters = assembly.run_simulation(n_steps=10000, dt=0.001)
# 可视化结果
fig = plt.figure(figsize=(16, 12))
# 1. 初始构型
ax1 = fig.add_subplot(2, 3, 1)
initial_pos = np.random.rand(200, 3) * 50.0 # 重新生成初始位置用于展示
ax1.scatter(initial_pos[:, 0], initial_pos[:, 1], s=50, c='blue', alpha=0.6)
ax1.set_xlim(0, 50)
ax1.set_ylim(0, 50)
ax1.set_xlabel('X (nm)', fontsize=11)
ax1.set_ylabel('Y (nm)', fontsize=11)
ax1.set_title('Initial Configuration', fontsize=12, fontweight='bold')
ax1.set_aspect('equal')
# 2. 最终构型
ax2 = fig.add_subplot(2, 3, 2)
# 根据团簇着色
colors = plt.cm.tab20(np.linspace(0, 1, len(final_clusters)))
for idx, cluster in enumerate(final_clusters):
cluster_pos = assembly.positions[cluster]
ax2.scatter(cluster_pos[:, 0], cluster_pos[:, 1],
s=80, c=[colors[idx]], alpha=0.8, edgecolors='black', linewidth=0.5)
ax2.set_xlim(0, 50)
ax2.set_ylim(0, 50)
ax2.set_xlabel('X (nm)', fontsize=11)
ax2.set_ylabel('Y (nm)', fontsize=11)
ax2.set_title('Final Configuration (Clusters)', fontsize=12, fontweight='bold')
ax2.set_aspect('equal')
# 3. 势能演化
ax3 = fig.add_subplot(2, 3, 3)
pe_history = np.array(assembly.potential_energy_history)
steps = np.arange(len(pe_history)) * 100
colors_pe = plt.cm.coolwarm((pe_history - pe_history.min()) / (pe_history.max() - pe_history.min()))
ax3.scatter(steps, pe_history, c=colors_pe, s=20, alpha=0.6)
ax3.plot(steps, pe_history, 'b-', alpha=0.3, linewidth=1)
ax3.set_xlabel('Step', fontsize=11)
ax3.set_ylabel('Potential Energy (eV)', fontsize=11)
ax3.set_title('Potential Energy Evolution', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
# 4. 团簇尺寸分布
ax4 = fig.add_subplot(2, 3, 4)
final_sizes = [len(c) for c in final_clusters]
unique_sizes, counts = np.unique(final_sizes, return_counts=True)
ax4.bar(unique_sizes, counts, color='steelblue', edgecolor='black', alpha=0.7)
ax4.set_xlabel('Cluster Size', fontsize=11)
ax4.set_ylabel('Count', fontsize=11)
ax4.set_title('Cluster Size Distribution', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3, axis='y')
# 5. 团簇数量演化
ax5 = fig.add_subplot(2, 3, 5)
n_clusters_history = [len(sizes) for sizes in assembly.clusters_history]
steps_clusters = np.arange(len(n_clusters_history)) * 100
ax5.plot(steps_clusters, n_clusters_history, 'g-', linewidth=2)
ax5.set_xlabel('Step', fontsize=11)
ax5.set_ylabel('Number of Clusters', fontsize=11)
ax5.set_title('Cluster Count Evolution', fontsize=12, fontweight='bold')
ax5.grid(True, alpha=0.3)
# 6. 平均团簇尺寸演化
ax6 = fig.add_subplot(2, 3, 6)
avg_sizes = [np.mean(sizes) for sizes in assembly.clusters_history]
ax6.plot(steps_clusters, avg_sizes, 'r-', linewidth=2)
ax6.set_xlabel('Step', fontsize=11)
ax6.set_ylabel('Average Cluster Size', fontsize=11)
ax6.set_title('Average Cluster Size Evolution', fontsize=12, fontweight='bold')
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/case2_self_assembly.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例2完成：纳米颗粒自组装模拟")
print(f" 结果保存在: {output_dir}/case2_self_assembly.png")
案例3：分子机器运动仿真
分子机器是能够在分子尺度执行机械功的人工分子器件。本案例模拟一个简化的分子马达模型。
物理模型：
- 转子：由3个质量点构成的三角形结构
- 定子：固定的环形轨道
- 驱动：通过周期性改变电场方向驱动转子转动
- 相互作用：Lennard-Jones势 + 库仑势
仿真目标：
- 观察转子的定向旋转
- 分析转速与驱动频率的关系
- 研究能量转换效率
Python实现：
# -*- coding: utf-8 -*-
"""
案例3：分子机器运动仿真
模拟光驱动分子马达的旋转运动
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle, FancyArrowPatch, Wedge
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题088\\output'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("\n案例3：分子机器运动仿真")
print("=" * 60)
class MolecularMotor:
"""分子马达模拟器"""
def __init__(self):
"""初始化分子马达系统"""
# 定子参数 (环形轨道)
self.stator_radius = 2.0 # nm
self.n_stator_sites = 12 # 定子位点数
# 转子参数 (三角形结构)
self.rotor_arm_length = 1.5 # nm
self.rotor_mass = 100.0 # amu
self.rotor_inertia = self.rotor_mass * self.rotor_arm_length**2 # 转动惯量
# 初始状态
self.theta = 0.0 # 转子角度
self.omega = 0.0 # 角速度
# 驱动参数
self.drive_frequency = 1.0 # GHz
self.drive_amplitude = 0.5 # eV/nm
# 相互作用参数
self.epsilon_lj = 0.1 # LJ势阱深度
self.sigma_lj = 0.3 # LJ特征长度
self.k_spring = 10.0 # 弹簧常数
# 记录数据
self.theta_history = []
self.omega_history = []
self.energy_history = []
self.torque_history = []
def compute_stator_potential(self, theta):
"""计算定子产生的周期性势场"""
# 余弦势场，n_stator_sites个势阱
V = -self.epsilon_lj * np.cos(self.n_stator_sites * theta)
return V
def compute_stator_torque(self, theta):
"""计算定子产生的力矩"""
# 力矩 = -dV/dtheta
tau = -self.epsilon_lj * self.n_stator_sites * np.sin(self.n_stator_sites * theta)
return tau
def compute_drive_torque(self, t):
"""计算驱动力矩 (周期性驱动)"""
# 通过改变电场方向驱动
drive_phase = 2 * np.pi * self.drive_frequency * t
tau_drive = self.drive_amplitude * np.sin(drive_phase)
return tau_drive
def compute_friction_torque(self):
"""计算摩擦力矩"""
gamma = 0.1 # 摩擦系数
tau_friction = -gamma * self.omega
return tau_friction
def step(self, dt, t):
"""时间步进"""
# 计算总力矩
tau_stator = self.compute_stator_torque(self.theta)
tau_drive = self.compute_drive_torque(t)
tau_friction = self.compute_friction_torque()
tau_total = tau_stator + tau_drive + tau_friction
# 更新角速度和角度
alpha = tau_total / self.rotor_inertia # 角加速度
self.omega += alpha * dt
self.theta += self.omega * dt
# 保持角度在[0, 2π)
self.theta = self.theta % (2 * np.pi)
# 计算能量
V_stator = self.compute_stator_potential(self.theta)
K_rot = 0.5 * self.rotor_inertia * self.omega**2
E_total = V_stator + K_rot
# 记录
self.theta_history.append(self.theta)
self.omega_history.append(self.omega)
self.energy_history.append(E_total)
self.torque_history.append(tau_total)
return E_total, tau_total
def get_rotor_positions(self):
"""获取转子三个臂的位置"""
positions = []
for i in range(3):
angle = self.theta + i * 2 * np.pi / 3
x = self.rotor_arm_length * np.cos(angle)
y = self.rotor_arm_length * np.sin(angle)
positions.append([x, y])
return np.array(positions)
# 运行模拟
print("\n初始化分子马达...")
motor = MolecularMotor()
print(f"定子半径: {motor.stator_radius} nm")
print(f"转子臂长: {motor.rotor_arm_length} nm")
print(f"驱动频率: {motor.drive_frequency} GHz")
# 模拟参数
total_time = 10.0 # ns
dt = 0.001 # ns
n_steps = int(total_time / dt)
print(f"\n运行模拟: {total_time} ns, {n_steps} 步")
for step in range(n_steps):
t = step * dt
motor.step(dt, t)
if step % 10000 == 0:
print(f" t = {t:.2f} ns, theta = {motor.theta:.3f} rad, omega = {motor.omega:.3f} rad/ns")
# 转换为numpy数组
motor.theta_history = np.array(motor.theta_history)
motor.omega_history = np.array(motor.omega_history)
motor.energy_history = np.array(motor.energy_history)
motor.torque_history = np.array(motor.torque_history)
time_array = np.arange(len(motor.theta_history)) * dt
print(f"\n模拟完成!")
print(f"平均转速: {np.mean(np.abs(motor.omega_history)):.3f} rad/ns")
print(f"总转数: {motor.theta_history[-1] / (2*np.pi):.1f} 圈")
# 可视化结果
fig = plt.figure(figsize=(16, 12))
# 1. 分子马达结构示意图
ax1 = fig.add_subplot(2, 3, 1)
# 绘制定子
for i in range(motor.n_stator_sites):
angle = 2 * np.pi * i / motor.n_stator_sites
x = motor.stator_radius * np.cos(angle)
y = motor.stator_radius * np.sin(angle)
circle = Circle((x, y), 0.15, color='blue', alpha=0.6)
ax1.add_patch(circle)
# 绘制转子
rotor_pos = motor.get_rotor_positions()
for pos in rotor_pos:
circle = Circle(pos, 0.2, color='red', alpha=0.8)
ax1.add_patch(circle)
# 连接线
ax1.plot(np.append(rotor_pos[:, 0], rotor_pos[0, 0]),
np.append(rotor_pos[:, 1], rotor_pos[0, 1]),
'r-', linewidth=2)
ax1.set_xlim(-3, 3)
ax1.set_ylim(-3, 3)
ax1.set_aspect('equal')
ax1.set_xlabel('X (nm)', fontsize=11)
ax1.set_ylabel('Y (nm)', fontsize=11)
ax1.set_title('Molecular Motor Structure', fontsize=12, fontweight='bold')
# 2. 角度演化
ax2 = fig.add_subplot(2, 3, 2)
# 展开角度（去除2π跳跃）
theta_unwrapped = np.unwrap(motor.theta_history)
ax2.plot(time_array, theta_unwrapped, 'b-', linewidth=1.5)
ax2.set_xlabel('Time (ns)', fontsize=11)
ax2.set_ylabel('Angle (rad)', fontsize=11)
ax2.set_title('Rotor Angle Evolution', fontsize=12, fontweight='bold')
ax2.grid(True, alpha=0.3)
# 3. 角速度演化
ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(time_array, motor.omega_history, 'r-', linewidth=1.5, alpha=0.7)
# 平滑曲线
window = 1000
if len(motor.omega_history) > window:
omega_smooth = np.convolve(motor.omega_history, np.ones(window)/window, mode='valid')
ax3.plot(time_array[window-1:], omega_smooth, 'g-', linewidth=2, label='Smoothed')
ax3.set_xlabel('Time (ns)', fontsize=11)
ax3.set_ylabel('Angular Velocity (rad/ns)', fontsize=11)
ax3.set_title('Angular Velocity', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
ax3.legend()
# 4. 能量演化
ax4 = fig.add_subplot(2, 3, 4)
ax4.plot(time_array, motor.energy_history, 'purple', linewidth=1.5, alpha=0.7)
ax4.set_xlabel('Time (ns)', fontsize=11)
ax4.set_ylabel('Energy (eV)', fontsize=11)
ax4.set_title('Total Energy', fontsize=12, fontweight='bold')
ax4.grid(True, alpha=0.3)
# 5. 力矩分析
ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(time_array, motor.torque_history, 'orange', linewidth=1, alpha=0.5)
# 计算平均力矩
window_torque = 5000
if len(motor.torque_history) > window_torque:
torque_smooth = np.convolve(motor.torque_history, np.ones(window_torque)/window_torque, mode='valid')
ax5.plot(time_array[window_torque-1:], torque_smooth, 'darkorange', linewidth=2, label='Average')
ax5.axhline(y=0, color='k', linestyle='--', linewidth=0.5)
ax5.set_xlabel('Time (ns)', fontsize=11)
ax5.set_ylabel('Torque (eV)', fontsize=11)
ax5.set_title('Applied Torque', fontsize=12, fontweight='bold')
ax5.grid(True, alpha=0.3)
ax5.legend()
# 6. 相空间轨迹
ax6 = fig.add_subplot(2, 3, 6)
ax6.plot(motor.theta_history, motor.omega_history, 'g.', markersize=1, alpha=0.3)
ax6.set_xlabel('Angle (rad)', fontsize=11)
ax6.set_ylabel('Angular Velocity (rad/ns)', fontsize=11)
ax6.set_title('Phase Space Trajectory', fontsize=12, fontweight='bold')
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/case3_molecular_motor.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例3完成：分子机器运动仿真")
print(f" 结果保存在: {output_dir}/case3_molecular_motor.png")
案例4：石墨烯电子性质计算
石墨烯是由碳原子以sp²杂化形成的二维蜂窝状晶格，具有独特的电子性质。本案例使用紧束缚模型计算石墨烯的能带结构。
物理模型：
- 紧束缚近似：只考虑最近邻相互作用
- 每个碳原子贡献一个p_z电子
- 哈密顿量矩阵在k空间构建
紧束缚模型：
H=−t∑<i,j>(ci†cj+h.c.)H = -t\sum_{<i,j>} (c_i^\dagger c_j + h.c.)H=−t<i,j>∑(ci†cj+h.c.)
其中t≈2.7t \approx 2.7t≈2.7 eV为跃迁积分，ci†c_i^\daggerci†和cjc_jcj为产生和湮灭算符。
Python实现：
# -*- coding: utf-8 -*-
"""
案例4：石墨烯电子性质计算
使用紧束缚模型计算能带结构和态密度
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题088\\output'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("\n案例4：石墨烯电子性质计算")
print("=" * 60)
class GrapheneTightBinding:
"""石墨烯紧束缚模型"""
def __init__(self, t=2.7, a=0.246):
"""
初始化参数
t: 最近邻跃迁积分 (eV)
a: 晶格常数 (nm)
"""
self.t = t
self.a = a
# 倒格子基矢
self.b1 = 2 * np.pi / a * np.array([1, 1/np.sqrt(3)])
self.b2 = 2 * np.pi / a * np.array([0, 2/np.sqrt(3)])
# 高对称点
self.Gamma = np.array([0, 0])
self.K = 2 * np.pi / a * np.array([1/3, 1/np.sqrt(3)])
self.M = 2 * np.pi / a * np.array([1/2, 1/(2*np.sqrt(3))])
def hamiltonian(self, kx, ky):
"""
构建k点哈密顿量
对于石墨烯，每个原胞有2个原子（A和B子晶格）
"""
# 最近邻矢量在实空间
delta1 = self.a * np.array([0, 1/np.sqrt(3)])
delta2 = self.a * np.array([-1/2, -1/(2*np.sqrt(3))])
delta3 = self.a * np.array([1/2, -1/(2*np.sqrt(3))])
# 跃迁项 f(k) = sum_j exp(i*k·delta_j)
k_vec = np.array([kx, ky])
f = np.exp(1j * np.dot(k_vec, delta1)) + \
np.exp(1j * np.dot(k_vec, delta2)) + \
np.exp(1j * np.dot(k_vec, delta3))
# 2x2 哈密顿量矩阵
H = np.array([[0, -self.t * f],
[-self.t * np.conj(f), 0]])
return H
def compute_bands(self, kx_array, ky_array):
"""计算能带"""
E_plus = np.zeros((len(kx_array), len(ky_array)))
E_minus = np.zeros((len(kx_array), len(ky_array)))
for i, kx in enumerate(kx_array):
for j, ky in enumerate(ky_array):
H = self.hamiltonian(kx, ky)
eigenvalues = np.linalg.eigvalsh(H)
E_plus[i, j] = eigenvalues[1]
E_minus[i, j] = eigenvalues[0]
return E_plus, E_minus
def band_structure_path(self, n_points=100):
"""计算高对称路径上的能带"""
# 路径: Gamma -> K -> M -> Gamma
path_segments = [
(self.Gamma, self.K),
(self.K, self.M),
(self.M, self.Gamma)
]
k_points = []
k_distances = []
energies_plus = []
energies_minus = []
total_dist = 0
for start, end in path_segments:
for i in range(n_points):
alpha = i / (n_points - 1)
k = (1 - alpha) * start + alpha * end
k_points.append(k)
H = self.hamiltonian(k[0], k[1])
eigenvalues = np.linalg.eigvalsh(H)
energies_plus.append(eigenvalues[1])
energies_minus.append(eigenvalues[0])
if i == 0 and len(k_distances) == 0:
k_distances.append(0)
elif i > 0:
dk = np.linalg.norm(k - k_points[-2])
total_dist += dk
k_distances.append(total_dist)
return np.array(k_distances), np.array(energies_plus), np.array(energies_minus)
def density_of_states(self, n_k=100, n_E=200):
"""计算态密度"""
# 在第一布里渊区采样
kx_vals = np.linspace(-np.pi/self.a, np.pi/self.a, n_k)
ky_vals = np.linspace(-np.pi/self.a, np.pi/self.a, n_k)
energies = []
for kx in kx_vals:
for ky in ky_vals:
H = self.hamiltonian(kx, ky)
eigenvalues = np.linalg.eigvalsh(H)
energies.extend(eigenvalues)
energies = np.array(energies)
# 直方图统计
E_min, E_max = energies.min(), energies.max()
E_bins = np.linspace(E_min, E_max, n_E)
dos, _ = np.histogram(energies, bins=E_bins, density=True)
E_centers = (E_bins[:-1] + E_bins[1:]) / 2
return E_centers, dos
# 创建模型
print("\n初始化石墨烯紧束缚模型...")
graphene = GrapheneTightBinding(t=2.7, a=0.246)
print(f"跃迁积分 t = {graphene.t} eV")
print(f"晶格常数 a = {graphene.a} nm")
# 计算布里渊区能带
print("\n计算能带结构...")
kx_array = np.linspace(-np.pi/graphene.a, np.pi/graphene.a, 100)
ky_array = np.linspace(-np.pi/graphene.a, np.pi/graphene.a, 100)
KX, KY = np.meshgrid(kx_array, ky_array)
E_plus, E_minus = graphene.compute_bands(kx_array, ky_array)
# 高对称路径能带
print("计算高对称路径能带...")
k_path, E_plus_path, E_minus_path = graphene.band_structure_path(n_points=150)
# 态密度
print("计算态密度...")
E_dos, dos = graphene.density_of_states(n_k=80, n_E=150)
print("\n能带计算完成!")
print(f"导带底能量: {E_plus.min():.3f} eV")
print(f"价带顶能量: {E_minus.max():.3f} eV")
print(f"能隙: {E_plus.min() - E_minus.max():.3f} eV (应为0，Dirac点)")
# 可视化结果
fig = plt.figure(figsize=(16, 12))
# 1. 3D能带图 - 导带
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
surf1 = ax1.plot_surface(KX, KY, E_plus, cmap='viridis', alpha=0.8)
ax1.set_xlabel('$k_x$ (nm$^{-1}$)', fontsize=10)
ax1.set_ylabel('$k_y$ (nm$^{-1}$)', fontsize=10)
ax1.set_zlabel('E (eV)', fontsize=10)
ax1.set_title('Conduction Band', fontsize=12, fontweight='bold')
fig.colorbar(surf1, ax=ax1, shrink=0.5)
# 2. 3D能带图 - 价带
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
surf2 = ax2.plot_surface(KX, KY, E_minus, cmap='plasma', alpha=0.8)
ax2.set_xlabel('$k_x$ (nm$^{-1}$)', fontsize=10)
ax2.set_ylabel('$k_y$ (nm$^{-1}$)', fontsize=10)
ax2.set_zlabel('E (eV)', fontsize=10)
ax2.set_title('Valence Band', fontsize=12, fontweight='bold')
fig.colorbar(surf2, ax=ax2, shrink=0.5)
# 3. 高对称路径能带
ax3 = fig.add_subplot(2, 3, 3)
ax3.plot(k_path, E_plus_path, 'b-', linewidth=2, label='Conduction')
ax3.plot(k_path, E_minus_path, 'r-', linewidth=2, label='Valence')
ax3.axhline(y=0, color='k', linestyle='--', linewidth=0.5)
# 标记高对称点
n_per_segment = 100
ax3.axvline(x=k_path[n_per_segment-1], color='gray', linestyle=':', alpha=0.5)
ax3.axvline(x=k_path[2*n_per_segment-1], color='gray', linestyle=':', alpha=0.5)
ax3.text(k_path[0], E_minus_path.min()-0.3, 'Γ', fontsize=12, ha='center')
ax3.text(k_path[n_per_segment-1], E_minus_path.min()-0.3, 'K', fontsize=12, ha='center')
ax3.text(k_path[2*n_per_segment-1], E_minus_path.min()-0.3, 'M', fontsize=12, ha='center')
ax3.text(k_path[-1], E_minus_path.min()-0.3, 'Γ', fontsize=12, ha='center')
ax3.set_xlabel('k-path', fontsize=11)
ax3.set_ylabel('Energy (eV)', fontsize=11)
ax3.set_title('Band Structure', fontsize=12, fontweight='bold')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 4. 态密度
ax4 = fig.add_subplot(2, 3, 4)
ax4.fill_between(E_dos, dos, alpha=0.5, color='green')
ax4.plot(E_dos, dos, 'g-', linewidth=2)
ax4.axvline(x=0, color='k', linestyle='--', linewidth=0.5, label='Fermi level')
ax4.set_xlabel('Energy (eV)', fontsize=11)
ax4.set_ylabel('DOS (a.u.)', fontsize=11)
ax4.set_title('Density of States', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)
# 5. Dirac锥放大图
ax5 = fig.add_subplot(2, 3, 5)
# 在K点附近采样
K_point = graphene.K
kx_dirac = np.linspace(K_point[0]-0.5, K_point[0]+0.5, 80)
ky_dirac = np.linspace(K_point[1]-0.5, K_point[1]+0.5, 80)
KX_D, KY_D = np.meshgrid(kx_dirac, ky_dirac)
E_plus_D, E_minus_D = graphene.compute_bands(kx_dirac, ky_dirac)
# 绘制等高线
levels = np.linspace(-1, 1, 21)
contour = ax5.contour(KX_D - K_point[0], KY_D - K_point[1], E_plus_D,
levels=levels, cmap='RdBu_r')
ax5.contour(KX_D - K_point[0], KY_D - K_point[1], E_minus_D,
levels=levels, cmap='RdBu_r')
ax5.set_xlabel('$k_x - K_x$ (nm$^{-1}$)', fontsize=11)
ax5.set_ylabel('$k_y - K_y$ (nm$^{-1}$)', fontsize=11)
ax5.set_title('Dirac Cone (Zoom near K)', fontsize=12, fontweight='bold')
ax5.set_aspect('equal')
plt.colorbar(contour, ax=ax5)
# 6. 费米速度计算
ax6 = fig.add_subplot(2, 3, 6)
# 沿K-Gamma方向计算色散
k_dir = np.linspace(0, 1.0, 100)
E_dirac = []
for dk in k_dir:
k_point = graphene.K + dk * (graphene.Gamma - graphene.K)
H = graphene.hamiltonian(k_point[0], k_point[1])
eigenvalues = np.linalg.eigvalsh(H)
E_dirac.append(eigenvalues)
E_dirac = np.array(E_dirac)
ax6.plot(k_dir, E_dirac[:, 0], 'ro-', markersize=3, linewidth=1, label='Valence')
ax6.plot(k_dir, E_dirac[:, 1], 'bo-', markersize=3, linewidth=1, label='Conduction')
# 线性拟合
linear_fit = np.polyfit(k_dir[:20], E_dirac[:20, 1], 1)
ax6.plot(k_dir[:30], np.polyval(linear_fit, k_dir[:30]), 'g--', linewidth=2, label=f'Linear fit, slope={linear_fit[0]:.2f}')
ax6.set_xlabel('|k - K| (a.u.)', fontsize=11)
ax6.set_ylabel('Energy (eV)', fontsize=11)
ax6.set_title('Linear Dispersion near Dirac Point', fontsize=12, fontweight='bold')
ax6.legend()
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/case4_graphene_bands.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例4完成：石墨烯电子性质计算")
print(f" 结果保存在: {output_dir}/case4_graphene_bands.png")
案例5：纳米流体传热仿真
纳米流体是将纳米颗粒分散在基液中形成的悬浮液，具有优异的传热性能。本案例模拟纳米流体的有效热导率。
物理模型：
- 采用Maxwell-Garnett有效介质理论
- 考虑颗粒-基液界面热阻
- 布朗运动对传热的贡献
有效热导率模型：
keffkf=kp+2kf+2ϕ(kp−kf)kp+2kf−ϕ(kp−kf)\frac{k_{eff}}{k_f} = \frac{k_p + 2k_f + 2\phi(k_p - k_f)}{k_p + 2k_f - \phi(k_p - k_f)}kfkeff=kp+2kf−ϕ(kp−kf)kp+2kf+2ϕ(kp−kf)
其中keffk_{eff}keff为有效热导率，kpk_pkp和kfk_fkf分别为颗粒和基液的热导率，ϕ\phiϕ为体积分数。
Python实现：
# -*- coding: utf-8 -*-
"""
案例5：纳米流体传热仿真
计算纳米流体的有效热导率和温度分布
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题088\\output'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("\n案例5：纳米流体传热仿真")
print("=" * 60)
class NanofluidHeatTransfer:
"""纳米流体传热仿真器"""
def __init__(self, k_f=0.6, k_p=40, rho_f=1000, rho_p=3900,
cp_f=4180, cp_p=765, mu_f=0.001):
"""
初始化纳米流体参数
k_f: 基液热导率 (W/m·K)
k_p: 颗粒热导率 (W/m·K)
rho_f, rho_p: 密度 (kg/m³)
cp_f, cp_p: 比热容 (J/kg·K)
mu_f: 基液粘度 (Pa·s)
"""
self.k_f = k_f
self.k_p = k_p
self.rho_f = rho_f
self.rho_p = rho_p
self.cp_f = cp_f
self.cp_p = cp_p
self.mu_f = mu_f
def maxwell_garnett(self, phi):
"""Maxwell-Garnett有效热导率模型"""
numerator = self.k_p + 2*self.k_f + 2*phi*(self.k_p - self.k_f)
denominator = self.k_p + 2*self.k_f - phi*(self.k_p - self.k_f)
k_eff = self.k_f * numerator / denominator
return k_eff
def hamilton_crosser(self, phi, n=3):
"""Hamilton-Crosser模型 (非球形颗粒)"""
psi = self.k_p / self.k_f
k_eff = self.k_f * (psi + (n-1) + (n-1)*phi*(psi-1)) / \
(psi + (n-1) - phi*(psi-1))
return k_eff
def brownian_motion_contribution(self, phi, T=300, d_p=50e-9):
"""布朗运动对热导率的贡献"""
k_B = 1.38e-23 # 玻尔兹曼常数
# 简化模型
k_brownian = 5 * k_B * T * phi * self.rho_p * self.cp_p / \
(2 * np.pi * self.mu_f * d_p)
return k_brownian
def effective_properties(self, phi):
"""计算纳米流体的有效物性"""
# 密度
rho_eff = (1 - phi) * self.rho_f + phi * self.rho_p
# 比热容 (质量加权)
cp_eff = ((1-phi)*self.rho_f*self.cp_f + phi*self.rho_p*self.cp_p) / rho_eff
# 热导率
k_eff = self.maxwell_garnett(phi)
# 热扩散系数
alpha_eff = k_eff / (rho_eff * cp_eff)
return rho_eff, cp_eff, k_eff, alpha_eff
def heat_conduction_1d(self, phi, L=0.01, T_hot=350, T_cold=300, nx=100):
"""
一维稳态热传导
L: 通道长度 (m)
T_hot, T_cold: 热端和冷端温度 (K)
"""
_, _, k_eff, _ = self.effective_properties(phi)
x = np.linspace(0, L, nx)
# 线性温度分布
T = T_hot - (T_hot - T_cold) * x / L
# 热流密度
q = k_eff * (T_hot - T_cold) / L
return x, T, q, k_eff
def transient_cooling(self, phi, L=0.01, T_init=350, T_ambient=300,
h=100, t_max=100, nt=1000):
"""
瞬态冷却过程
h: 对流换热系数 (W/m²·K)
"""
rho, cp, k, alpha = self.effective_properties(phi)
dx = L / 50
dt = 0.1
nx = int(L / dx) + 1
nt = int(t_max / dt) + 1
x = np.linspace(0, L, nx)
T = np.ones(nx) * T_init
T[0] = T_ambient # 边界条件
T[-1] = T_ambient
T_history = [T.copy()]
for n in range(1, nt):
T_new = T.copy()
for i in range(1, nx-1):
# 有限差分
T_new[i] = T[i] + alpha * dt / dx**2 * (T[i+1] - 2*T[i] + T[i-1])
T = T_new
if n % 100 == 0:
T_history.append(T.copy())
return x, T_history
# 创建纳米流体模型
print("\n初始化纳米流体模型...")
nanofluid = NanofluidHeatTransfer()
print(f"基液热导率: {nanofluid.k_f} W/m·K")
print(f"颗粒热导率: {nanofluid.k_p} W/m·K")
# 计算不同体积分数下的热导率
print("\n计算有效热导率...")
phi_range = np.linspace(0, 0.1, 50) # 0-10%体积分数
k_maxwell = []
k_hamilton = []
for phi in phi_range:
k_maxwell.append(nanofluid.maxwell_garnett(phi))
k_hamilton.append(nanofluid.hamilton_crosser(phi, n=3))
k_maxwell = np.array(k_maxwell)
k_hamilton = np.array(k_hamilton)
# 一维稳态热传导
print("计算稳态温度分布...")
x_0, T_0, q_0, k_0 = nanofluid.heat_conduction_1d(phi=0)
x_5, T_5, q_5, k_5 = nanofluid.heat_conduction_1d(phi=0.05)
x_10, T_10, q_10, k_10 = nanofluid.heat_conduction_1d(phi=0.10)
print(f"纯基液热导率: {k_0:.3f} W/m·K")
print(f"5%纳米流体热导率: {k_5:.3f} W/m·K (提升 {(k_5/k_0-1)*100:.1f}%)")
print(f"10%纳米流体热导率: {k_10:.3f} W/m·K (提升 {(k_10/k_0-1)*100:.1f}%)")
# 瞬态冷却
print("\n计算瞬态冷却过程...")
x_trans, T_history = nanofluid.transient_cooling(phi=0.05, t_max=50)
# 可视化结果
fig = plt.figure(figsize=(16, 12))
# 1. 有效热导率 vs 体积分数
ax1 = fig.add_subplot(2, 3, 1)
ax1.plot(phi_range * 100, k_maxwell / nanofluid.k_f, 'b-', linewidth=2, label='Maxwell-Garnett')
ax1.plot(phi_range * 100, k_hamilton / nanofluid.k_f, 'r--', linewidth=2, label='Hamilton-Crosser')
ax1.set_xlabel('Volume Fraction (%)', fontsize=11)
ax1.set_ylabel('$k_{eff}/k_f$', fontsize=11)
ax1.set_title('Effective Thermal Conductivity', fontsize=12, fontweight='bold')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 2. 稳态温度分布
ax2 = fig.add_subplot(2, 3, 2)
ax2.plot(x_0 * 1000, T_0, 'b-', linewidth=2, label=f'Pure fluid (k={k_0:.2f})')
ax2.plot(x_5 * 1000, T_5, 'r-', linewidth=2, label=f'5% nanofluid (k={k_5:.2f})')
ax2.plot(x_10 * 1000, T_10, 'g-', linewidth=2, label=f'10% nanofluid (k={k_10:.2f})')
ax2.set_xlabel('Position (mm)', fontsize=11)
ax2.set_ylabel('Temperature (K)', fontsize=11)
ax2.set_title('Steady-State Temperature', fontsize=12, fontweight='bold')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 3. 热流密度对比
ax3 = fig.add_subplot(2, 3, 3)
phi_test = np.array([0, 0.02, 0.04, 0.06, 0.08, 0.10])
q_values = []
for phi in phi_test:
_, _, q, _ = nanofluid.heat_conduction_1d(phi)
q_values.append(q)
q_values = np.array(q_values)
ax3.bar(phi_test * 100, q_values / q_values[0], color='steelblue', edgecolor='black', alpha=0.7)
ax3.set_xlabel('Volume Fraction (%)', fontsize=11)
ax3.set_ylabel('Heat Flux / Base', fontsize=11)
ax3.set_title('Heat Transfer Enhancement', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3, axis='y')
# 4. 瞬态温度分布
ax4 = fig.add_subplot(2, 3, 4)
times = [0, 2, 5, 10, 20, len(T_history)-1]
colors = plt.cm.coolwarm(np.linspace(0, 1, len(times)))
for i, t_idx in enumerate(times):
if t_idx < len(T_history):
ax4.plot(x_trans * 1000, T_history[t_idx], color=colors[i],
linewidth=2, label=f't = {t_idx*0.1*100:.0f} s')
ax4.set_xlabel('Position (mm)', fontsize=11)
ax4.set_ylabel('Temperature (K)', fontsize=11)
ax4.set_title('Transient Cooling (5% nanofluid)', fontsize=12, fontweight='bold')
ax4.legend(fontsize=9)
ax4.grid(True, alpha=0.3)
# 5. 有效物性汇总
ax5 = fig.add_subplot(2, 3, 5)
phi_props = np.linspace(0, 0.1, 20)
rho_eff_list = []
cp_eff_list = []
alpha_eff_list = []
for phi in phi_props:
rho, cp, k, alpha = nanofluid.effective_properties(phi)
rho_eff_list.append(rho)
cp_eff_list.append(cp)
alpha_eff_list.append(alpha * 1e6) # 放大显示
ax5_twin = ax5.twinx()
line1 = ax5.plot(phi_props * 100, rho_eff_list, 'b-', linewidth=2, label='Density')
line2 = ax5.plot(phi_props * 100, cp_eff_list, 'r-', linewidth=2, label='Specific Heat')
line3 = ax5_twin.plot(phi_props * 100, alpha_eff_list, 'g--', linewidth=2, label='Thermal Diffusivity (×10⁶)')
ax5.set_xlabel('Volume Fraction (%)', fontsize=11)
ax5.set_ylabel('ρ (kg/m³), cp (J/kg·K)', fontsize=11)
ax5_twin.set_ylabel('α (×10⁻⁶ m²/s)', fontsize=11, color='g')
ax5.set_title('Effective Properties', fontsize=12, fontweight='bold')
ax5.grid(True, alpha=0.3)
# 合并图例
lines = line1 + line2 + line3
labels = [l.get_label() for l in lines]
ax5.legend(lines, labels, loc='center right')
# 6. 纳米颗粒分布示意图
ax6 = fig.add_subplot(2, 3, 6)
# 绘制简化的纳米流体微观结构
np.random.seed(42)
n_particles = 30
box_size = 10
particle_positions = np.random.rand(n_particles, 2) * box_size
particle_radii = np.random.uniform(0.2, 0.4, n_particles)
for pos, radius in zip(particle_positions, particle_radii):
circle = Circle(pos, radius, color='red', alpha=0.6, edgecolor='darkred')
ax6.add_patch(circle)
ax6.set_xlim(0, box_size)
ax6.set_ylim(0, box_size)
ax6.set_aspect('equal')
ax6.set_xlabel('X (nm)', fontsize=11)
ax6.set_ylabel('Y (nm)', fontsize=11)
ax6.set_title('Nanoparticle Distribution', fontsize=12, fontweight='bold')
plt.tight_layout()
plt.savefig(f'{output_dir}/case5_nanofluid_heat.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例5完成：纳米流体传热仿真")
print(f" 结果保存在: {output_dir}/case5_nanofluid_heat.png")
案例6：DNA分子动力学模拟
DNA是存储遗传信息的生物大分子，其结构和动力学对理解生命过程至关重要。本案例模拟DNA双螺旋的简化模型。
物理模型：
- 粗粒化模型：每个碱基对用一个珠子表示
- 碱基对间通过谐振势连接
- 考虑氢键相互作用
- 溶剂效应通过随机力模拟
Python实现：
# -*- coding: utf-8 -*-
"""
案例6：DNA分子动力学模拟
模拟DNA双螺旋的简化粗粒化模型
"""
import matplotlib
matplotlib.use('Agg')
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.patches import FancyArrowPatch
from mpl_toolkits.mplot3d import proj3d
import warnings
warnings.filterwarnings('ignore')
import os
output_dir = r'd:\\文档\\500仿真领域\\工程仿真\\多场耦合优化\\主题088\\output'
os.makedirs(output_dir, exist_ok=True)
plt.rcParams['font.sans-serif'] = ['SimHei', 'DejaVu Sans']
plt.rcParams['axes.unicode_minus'] = False
print("\n案例6：DNA分子动力学模拟")
print("=" * 60)
class DNAMolecularDynamics:
"""DNA分子动力学模拟器（粗粒化模型）"""
def __init__(self, n_base_pairs=50):
"""
初始化DNA
n_base_pairs: 碱基对数量
"""
self.N = n_base_pairs
# DNA几何参数
self.rise = 0.34 # 碱基对上升距离 (nm)
self.twist = 36 * np.pi / 180 # 每步扭转角 (rad)
self.radius = 1.0 # 螺旋半径 (nm)
# 力场参数
self.k_backbone = 100.0 # 骨架弹簧常数
self.k_hbond = 50.0 # 氢键弹簧常数
self.k_bending = 30.0 # 弯曲刚度
self.k_torsion = 20.0 # 扭转刚度
# 质量
self.mass = 300.0 # amu (每个碱基对)
# 初始化位置
self.build_dna_structure()
# 速度初始化
self.velocities = np.zeros((self.N, 3))
self.forces = np.zeros((self.N, 3))
# 记录
self.trajectory = []
self.energy_history = []
def build_dna_structure(self):
"""构建B-DNA双螺旋结构"""
self.positions = np.zeros((self.N, 3))
for i in range(self.N):
theta = i * self.twist
z = i * self.rise
# 双螺旋位置
self.positions[i, 0] = self.radius * np.cos(theta)
self.positions[i, 1] = self.radius * np.sin(theta)
self.positions[i, 2] = z
# 保存初始结构
self.initial_positions = self.positions.copy()
def compute_forces(self):
"""计算力"""
self.forces.fill(0)
potential = 0
# 1. 骨架连接力 (相邻碱基对)
for i in range(self.N - 1):
r_vec = self.positions[i+1] - self.positions[i]
r = np.linalg.norm(r_vec)
if r > 0:
# 谐振势
dr = r - self.rise
f_mag = -self.k_backbone * dr
f_vec = f_mag * r_vec / r
self.forces[i] -= f_vec
self.forces[i+1] += f_vec
potential += 0.5 * self.k_backbone * dr**2
# 2. 弯曲力 (三体相互作用)
for i in range(1, self.N - 1):
r1 = self.positions[i] - self.positions[i-1]
r2 = self.positions[i+1] - self.positions[i]
# 计算角度
cos_angle = np.dot(r1, r2) / (np.linalg.norm(r1) * np.linalg.norm(r2))
cos_angle = np.clip(cos_angle, -1, 1)
angle = np.arccos(cos_angle)
# 弯曲势能
potential += 0.5 * self.k_bending * angle**2
# 3. 扭转力 (恢复原始扭转)
for i in range(self.N - 1):
# 简化的扭转恢复力
current_twist = np.arctan2(self.positions[i+1, 1], self.positions[i+1, 0]) - \
np.arctan2(self.positions[i, 1], self.positions[i, 0])
twist_diff = current_twist - self.twist
# 切向力
tangent = np.array([-self.positions[i, 1], self.positions[i, 0], 0])
tangent = tangent / (np.linalg.norm(tangent) + 1e-10)
f_torsion = -self.k_torsion * twist_diff * tangent
self.forces[i] += f_torsion
self.forces[i+1] -= f_torsion
# 4. 氢键力 (保持双链配对 - 简化为径向恢复力)
for i in range(self.N):
r_xy = np.array([self.positions[i, 0], self.positions[i, 1]])
current_radius = np.linalg.norm(r_xy)
if current_radius > 0:
dr = current_radius - self.radius
f_radial = -self.k_hbond * dr * r_xy / current_radius
self.forces[i, 0] += f_radial[0]
self.forces[i, 1] += f_radial[1]
potential += 0.5 * self.k_hbond * dr**2
return potential
def velocity_verlet(self, dt):
"""速度Verlet积分"""
self.positions += self.velocities * dt + 0.5 * self.forces / self.mass * dt**2
forces_old = self.forces.copy()
pe = self.compute_forces()
self.velocities += 0.5 * (forces_old + self.forces) / self.mass * dt
ke = 0.5 * self.mass * np.sum(self.velocities**2)
return pe, ke
def langevin_dynamics(self, dt, gamma=0.1, T=300):
"""朗之万动力学（包含溶剂效应）"""
k_B = 8.617e-5 # eV/K
# 确定性力
pe = self.compute_forces()
# 随机力 (白噪声)
random_force = np.random.randn(self.N, 3) * np.sqrt(2 * gamma * k_B * T / dt)
# 总力
total_force = self.forces + random_force - gamma * self.velocities * self.mass
# 更新
self.velocities += total_force / self.mass * dt
self.positions += self.velocities * dt
ke = 0.5 * self.mass * np.sum(self.velocities**2)
return pe, ke
def compute_persistence_length(self):
"""计算持久长度"""
# 计算切向量相关函数
tangents = np.zeros((self.N - 1, 3))
for i in range(self.N - 1):
t = self.positions[i+1] - self.positions[i]
tangents[i] = t / np.linalg.norm(t)
# 相关函数 <t(s)·t(0)> = exp(-s/Lp)
correlations = []
distances = []
for s in range(1, min(20, self.N//2)):
corr = np.mean([np.dot(tangents[i], tangents[i+s]) for i in range(self.N-1-s)])
correlations.append(corr)
distances.append(s * self.rise)
# 拟合指数衰减
if len(correlations) > 2:
log_corr = np.log(np.maximum(correlations, 0.01))
fit = np.polyfit(distances, log_corr, 1)
Lp = -1 / fit[0] # 持久长度
else:
Lp = 50 # 默认值
return Lp, distances, correlations
def run_simulation(self, n_steps=5000, dt=0.001):
"""运行模拟"""
print(f"\n运行DNA动力学模拟...")
print(f"碱基对数: {self.N}")
print(f"模拟步数: {n_steps}")
sample_interval = 50
for step in range(n_steps):
pe, ke = self.langevin_dynamics(dt, gamma=0.5, T=300)
if step % sample_interval == 0:
self.trajectory.append(self.positions.copy())
self.energy_history.append(pe + ke)
if step % 1000 == 0:
print(f" Step {step}/{n_steps}, Energy = {pe+ke:.2f} eV")
print(f"\n模拟完成!")
return self.trajectory
# 运行模拟
print("\n初始化DNA结构...")
dna = DNAMolecularDynamics(n_base_pairs=40)
print(f"碱基对数量: {dna.N}")
print(f"螺旋半径: {dna.radius} nm")
print(f"每步上升: {dna.rise} nm")
print(f"每步扭转: {dna.twist * 180/np.pi:.1f}°")
trajectory = dna.run_simulation(n_steps=3000, dt=0.001)
# 计算持久长度
print("\n计算持久长度...")
Lp, distances, correlations = dna.compute_persistence_length()
print(f"持久长度: {Lp:.1f} nm")
# 可视化结果
fig = plt.figure(figsize=(16, 12))
# 1. DNA 3D结构 - 初始
ax1 = fig.add_subplot(2, 3, 1, projection='3d')
pos_init = dna.initial_positions
ax1.plot(pos_init[:, 0], pos_init[:, 1], pos_init[:, 2],
'b-', linewidth=3, label='Backbone')
ax1.scatter(pos_init[:, 0], pos_init[:, 1], pos_init[:, 2],
c=np.arange(dna.N), cmap='viridis', s=50)
ax1.set_xlabel('X (nm)', fontsize=10)
ax1.set_ylabel('Y (nm)', fontsize=10)
ax1.set_zlabel('Z (nm)', fontsize=10)
ax1.set_title('Initial DNA Structure', fontsize=12, fontweight='bold')
# 2. DNA 3D结构 - 最终
ax2 = fig.add_subplot(2, 3, 2, projection='3d')
pos_final = dna.positions
ax2.plot(pos_final[:, 0], pos_final[:, 1], pos_final[:, 2],
'r-', linewidth=3, label='Backbone')
ax2.scatter(pos_final[:, 0], pos_final[:, 1], pos_final[:, 2],
c=np.arange(dna.N), cmap='plasma', s=50)
ax2.set_xlabel('X (nm)', fontsize=10)
ax2.set_ylabel('Y (nm)', fontsize=10)
ax2.set_zlabel('Z (nm)', fontsize=10)
ax2.set_title('Final DNA Structure', fontsize=12, fontweight='bold')
# 3. 能量演化
ax3 = fig.add_subplot(2, 3, 3)
energy_array = np.array(dna.energy_history)
steps = np.arange(len(energy_array)) * 50
ax3.plot(steps, energy_array, 'g-', linewidth=1.5)
ax3.set_xlabel('Step', fontsize=11)
ax3.set_ylabel('Total Energy (eV)', fontsize=11)
ax3.set_title('Energy Evolution', fontsize=12, fontweight='bold')
ax3.grid(True, alpha=0.3)
# 4. 持久长度计算
ax4 = fig.add_subplot(2, 3, 4)
if len(distances) > 0:
ax4.plot(distances, correlations, 'bo', markersize=6, label='Simulation')
# 指数拟合
fit_exp = np.exp(-np.array(distances) / Lp)
ax4.plot(distances, fit_exp, 'r-', linewidth=2, label=f'Fit: Lp={Lp:.1f} nm')
ax4.set_xlabel('Distance (nm)', fontsize=11)
ax4.set_ylabel('Correlation', fontsize=11)
ax4.set_title('Tangent Correlation', fontsize=12, fontweight='bold')
ax4.legend()
ax4.grid(True, alpha=0.3)
# 5. 末端距分布
ax5 = fig.add_subplot(2, 3, 5)
end_to_end_distances = []
for pos in trajectory:
r_end = pos[-1] - pos[0]
end_to_end_distances.append(np.linalg.norm(r_end))
end_to_end_distances = np.array(end_to_end_distances)
ax5.hist(end_to_end_distances, bins=20, color='steelblue', edgecolor='black', alpha=0.7)
ax5.axvline(x=np.mean(end_to_end_distances), color='r', linestyle='--',
linewidth=2, label=f'Mean: {np.mean(end_to_end_distances):.1f} nm')
ax5.set_xlabel('End-to-End Distance (nm)', fontsize=11)
ax5.set_ylabel('Frequency', fontsize=11)
ax5.set_title('End-to-End Distance Distribution', fontsize=12, fontweight='bold')
ax5.legend()
ax5.grid(True, alpha=0.3, axis='y')
# 6. 螺旋参数演化
ax6 = fig.add_subplot(2, 3, 6)
# 计算每步的半径和扭转
radii = []
twists = []
for pos in trajectory[::10]: # 每10帧采样
# 平均半径
r_xy = np.sqrt(pos[:, 0]**2 + pos[:, 1]**2)
radii.append(np.mean(r_xy))
# 平均扭转角
angles = np.arctan2(pos[:, 1], pos[:, 0])
twist_per_step = np.mean(np.diff(angles))
twists.append(twist_per_step * 180 / np.pi)
frames = np.arange(len(radii)) * 10 * 50
ax6.plot(frames, radii, 'b-', linewidth=2, label='Radius')
ax6.axhline(y=dna.radius, color='b', linestyle='--', alpha=0.5)
ax6_twin = ax6.twinx()
ax6_twin.plot(frames, twists, 'r-', linewidth=2, label='Twist angle')
ax6_twin.axhline(y=dna.twist*180/np.pi, color='r', linestyle='--', alpha=0.5)
ax6.set_xlabel('Step', fontsize=11)
ax6.set_ylabel('Radius (nm)', fontsize=11, color='b')
ax6_twin.set_ylabel('Twist angle (°)', fontsize=11, color='r')
ax6.set_title('Helical Parameters', fontsize=12, fontweight='bold')
ax6.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(f'{output_dir}/case6_dna_dynamics.png', dpi=150, bbox_inches='tight')
plt.close()
print("\n✓ 案例6完成：DNA分子动力学模拟")
print(f" 结果保存在: {output_dir}/case6_dna_dynamics.png")
AtomGit 是由开放原子开源基金会联合 CSDN 等生态伙伴共同推出的新一代开源与人工智能协作平台。平台坚持“开放、中立、公益”的理念，把代码托管、模型共享、数据集托管、智能体开发体验和算力服务整合在一起，为开发者提供从开发、训练到部署的一站式体验。
更多推荐
0
0- 0
已为社区贡献83条内容
所有评论(0)