当前位置:首页>java>[VibeCoding,26,快速总结代码和复现结果]聚合物粗粒化模拟: Optuna+BO+Lammps

[VibeCoding,26,快速总结代码和复现结果]聚合物粗粒化模拟: Optuna+BO+Lammps

  • 2026-09-08 13:03:30
[VibeCoding,26,快速总结代码和复现结果]聚合物粗粒化模拟: Optuna+BO+Lammps
本文旨在通过VibeCoding的方式快速了解和复现开源项目或者文献的结果Martins, C. A., Jr. et al. Bayesian Optimization for High-Dimensional Coarse-Grained Model Parameterization: A Case Study on Pebax Polymer. Journal of Chemical Theory and Computation (2026). https://doi.org:10.1021/acs.jctc.5c01500为例,原代码地址为https://github.com/camjjr/bo_cgff

字数 1945,阅读大约需 10 分钟

1. 项目核心目标

一句话总结:使用贝叶斯优化(Bayesian Optimization)自动优化粗粒度(Coarse-Grained, CG)分子动力学模型的力场参数,使其在密度(ρ)、回转半径(Rg)和玻璃化转变温度(Tg)三个物理性质上与原子级参考数据匹配。


2. 采用的主要算法/方法/技术栈

类别
技术/方法
用途
优化算法
Optuna (TPE采样器)
高维参数空间的贝叶斯优化
分子动力学引擎
LAMMPS (mpirun并行)
执行CG-MD模拟计算物理性质
力场类型
Mie势 + 谐波键/角势
描述非键相互作用和成键相互作用
数据分析
scikit-learn LinearRegression
从密度-温度数据提取Tg(斜率分析)
参数注入
Unix sed命令
动态替换LAMMPS输入文件中的占位符
目标函数设计
多目标加权误差平方和
同时优化密度、Rg、Tg三个性质

3. 整体代码目录结构

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

bo_cgff/├── bo_cg_search.py           # 主优化脚本(211行)- 核心入口├── utils.py                  # 工具函数(17行)- Rg和均值计算├── requirements.txt          # Python依赖├── parameters_search.dat     # LAMMPS参数模板(含占位符)├── input.dat                 # LAMMPS主模拟输入脚本├── input_vaccum.dat          # LAMMPS真空参考模拟输入├── structure.dat             # 分子结构文件├── parameters.dat            # 初始参数值├── README.md                 # 项目文档│├── Atomistic_PEBAX_50chains/ # 原子级参考数据│   ├── input.dat│   ├── structure.dat│   ├── parameters.dat│   └── README.md│└── optimal_parameters_simulation/  # 优化结果验证    ├── get_mean.py           # 计算平均值工具    ├── input.dat    ├── structure.dat    ├── parameters.dat    │    └── Rg_vaccum/            # 真空Rg计算        ├── get_rg.py         # 提取回转半径        └── get_mean.py


4. 最核心的执行流程 / Pipeline

否

是

启动Optuna优化循环n_trials=3000
生成候选参数new_parameters
使用sed注入参数parameters_search_aux.dat
并行执行LAMMPS模拟mpirun -np 6/1
主模拟: input.dat真空模拟: input_vaccum.dat
解析输出文件*_averages.txt, *_gyration_0_vaccum.txt
计算目标函数objective
密度误差线性回归提取斜率
Rg误差与原子级参考比较
Tg惩罚斜率差异
加权求和objective = error_density + error_slope + error_rg
返回目标值Optuna更新代理模型
是否达到n_trials?
输出最优参数sqlite数据库

详细步骤:

  1. 1. 初始化:创建Optuna Study,配置TPE采样器,SQLite持久化存储
  2. 2. 参数生成(new_parameters):定义37维搜索空间(6键+14角+5σ+5ε+5γ+6键长),使用suggest_float采样
  3. 3. 参数注入:复制模板文件→使用sed批量替换37个占位符→生成parameters_search_aux.dat
  4. 4. 并行模拟:
    • • 主模拟(6核):计算不同温度下的密度
    • • 真空模拟(1核):计算孤立链的Rg
  5. 5. 数据提取:
    • • 密度:读取2.28-2.31和2.39-2.42共8个温度点的平均值文件
    • • Rg:读取2.28-2.42共15个链的回转半径
  6. 6. 目标计算:
    • • 密度误差:低温区(150-225K)和高温区(425-500K)线性插值后比较
    • • 斜率惩罚:(a_sim - a_target)² / a_target²,反映Tg准确性
    • • Rg误差:15个链的相对误差平方和,权重因子2×len(density_range)/len(rg_chains)
  7. 7. 迭代优化:Optuna根据历史结果更新代理模型,指导下一次采样
  8. 8. 结果保存:最优参数保存至SQLite数据库,支持中断恢复

5. 关键类的职责或主要函数的作用

函数/模块
位置
职责
关键实现细节
new_parameters(trial)
bo_cg_search.py:11-112
参数空间定义与注入
37个suggest_float调用,6类参数(bond_k, angle_k/theta, sigma, epsilon, gamma_r, bond_l),使用sed进行37次占位符替换
objective(trial)
bo_cg_search.py:115-190
目标函数计算
调用LAMMPS→解析输出→计算三项目标→返回标量损失
get_mean(filename)
utils.py:15-17
密度数据读取np.loadtxt
读取LAMMPS输出,取第二列均值
get_rg()
utils.py:5-13
回转半径提取
遍历15个文件(2.28-2.42),跳过前5行取均值
LinearRegression
sklearn
Tg提取
对密度-温度数据拟合直线,斜率反映Tg
optuna.create_study()
main块
优化器管理
TPE采样、SQLite持久化、load_if_exists支持断点续跑
subprocess.run(sed...)
bo_cg_search.py:65-111
模板渲染
37个sed命令顺序执行,替换parameters_search_aux.dat

6. 参数/配置/搜索空间的组织

搜索空间维度:37维连续参数

参数类别
数量
变量名模式
搜索范围
步长
键力常数
6
bond_tX_tY_k
[1.0, 12.0]
0.001
角力常数
7
angle_tX_tY_tZ_k
[1.0, 12.0]
0.001
平衡角
7
angle_tX_tY_tZ_theta
[70.0, 180.0]
0.001
Mie σ
5
sigma_tN
[2.5, 5.0]
0.001
Mie ε
5
epsilon_tN
[0.3, 1.5]
0.001
Mie γ
5
gamma_r_tN
[8, 16]
0.001
键平衡长度
6
l_tXtY
[2.5, 5.0]
0.001

配置组织方式:

  • • 模板文件:parameters_search.dat使用占位符(如k_bond_t1t2, sigma_t1)
  • • 运行时生成:每次trial复制模板→sed替换→生成parameters_search_aux.dat
  • • 硬编码目标值:密度参考值直接写在代码中(y_1=[1.2344, 1.2288, 1.2187, 1.2115])
  • • 原子级Rg参考:15个链的Rg值硬编码在objective()函数中

7. 设计亮点与取舍

亮点

  1. 1. 模板-注入架构的简洁性
    • • 使用sed进行参数替换而非Python字符串操作,保持LAMMPS文件格式完整性
    • • 37个参数通过命名占位符管理,可读性强于索引数组
  2. 2. 多物理量联合优化的权重设计
    • • 密度误差和Rg误差通过2×len(density_range)/len(rg_chains)因子平衡量级
    • • 斜率惩罚项直接关联Tg物理意义,避免单独优化Tg点
  3. 3. 断点续跑机制
    • • load_if_exists=True + SQLite持久化,支持3000次trial的长期运行中断恢复

取舍

  1. 1. 硬编码 vs 配置文件
    • • 目标密度值、Rg参考值、温度点全部硬编码,修改需改源码
    • • 取舍原因:减少文件IO,但牺牲灵活性
  2. 2. 同步执行 vs 异步队列
    • • 每次trial阻塞等待LAMMPS完成(~3分钟×3000次≈150小时)
    • • 取舍原因:实现简单,但未利用Optuna的异步优化能力
  3. 3. 全局搜索 vs 分阶段优化
    • • 37维同时优化,未采用先优化键/角参数再优化非键参数的分层策略
    • • 取舍原因:依赖TPE采样器处理高维空间,但收敛可能较慢

8. 可能的优化方向

  1. 1. 参数外置化
    • • 将目标值(密度、Rg参考数据)、温度点、权重因子移至YAML/JSON配置文件
    • • 支持不同材料体系无需修改源码
  2. 2. 并行Trial执行
    • • 使用Optuna的n_jobs参数或集成Joblib/Dask,同时运行多个trial
    • • 潜在加速比:6-10×(取决于CPU核心数)
  3. 3. 早停机制(Early Stopping)
    • • 集成Optuna的MedianPruner或HyperbandPruner,在trial中途判断 hopeless 并终止
    • • 可节省30-50%计算资源(基于MD模拟中途物理量已明显偏离目标)
  4. 4. 搜索空间降维
    • • 分析参数相关性(如相同原子类型的σ/ε可能相关),使用PCA或手动分组减少独立维度
    • • 或采用两阶段优化:先粗粒度网格搜索确定大致范围,再BO精细优化

总结:这是一个针对聚合物材料(Pebax)CG模型参数化的专业工具,采用Optuna+LAMMPS的经典组合,代码结构清晰但存在硬编码和串行执行的可优化空间。核心创新在于通过线性回归斜率间接优化Tg,避免了玻璃化转变点的精确定义难题。


最新文章

随机文章