新建
使用深度势能分子动力学进行固态电解质研究实战
AnguseZhang
推荐镜像 :DeePMD-kit:2.2.1-cuda11.6-notebook
推荐机型 :c12_m46_1 * NVIDIA GPU B
赞 10
5
29
目录
使用深度势能分子动力学进行固态电解质研究实战
代码
文本
©️ Copyright 2023 @ Authors
作者： 宋哲轩📨 黄剑兴 📨 张与之 📨
日期：2023-05-19
共享协议：本作品采用知识共享署名-非商业性使用-相同方式共享 4.0 国际许可协议进行许可。
代码
文本
🎯 本文档旨在结合J. Chem. Phys. 154, 094703 (2021)`这篇于2021年发表的使用深度势能进行固态电解质研究工作，介绍机器学习与分子动力学的结合是如何发挥作用的。
你可以点击界面上方蓝色按钮 开始连接，选择 deepmd-kit:2.2.1-cuda11.6 镜像及c12_m46_1 * NVIDIA GPU B节点配置，稍等片刻即可运行。
代码
文本
代码
文本
目标
使用 LAMMPS 运行深度势能分子动力学。
在学习本教程后，你将能够：
- 结合分子动力学软件LAMMPS了解分子动力学是如何工作的。
- 计算并绘制体系的径向分布图像(RDF)。
- 计算体系的离子扩散系数D，并绘制均方位移(MSD)。
阅读并动手运行该教程【最多】约需 20 分钟，让我们开始吧！
代码
文本
背景介绍
锂离子在正极-电解质-负极间往复运动。今天的商用锂离子电池主要由正极、电解液、隔膜和负极组成。有机电解液材料的易燃性是造成电池热失控的根本原因。使用不易燃的无机固体电解质的全固态电池体系，可以从根本上解决电池体系的安全问题，同时实现较高的能量密度。
于此同时，要保证“充电一分钟通话两小时”，就需要固体电解质材料具有较高的离子电导率和扩散系数。以LiGePS类无机固体电解质为代表的硫化物材料是目前公认的下一代电池的主要技术路线。要优化并设计新型的固体电解质材料，研究者们需要从根本上理解锂离子在电池材料中的扩散行为，尤其是存在协同扩散行为的固体电解质材料。固体核磁，中子衍射等实验手段可以从原子尺度揭示锂离子在不同位点间扩散的方式，然而，由于成本和研究效率限制，实验表征手段不适宜于高通量的材料设计，例如材料组分精确调控。与之相对，低成本的原子尺度的计算模拟可以提供微观尺度的物理图像，为研究者们理解材料反应，设计优化电池材料提供了理性依据。
目前，对电池材料的扩散行为的模拟主要基于两类方法：基于NEB方法的热力学计算和基于第一性原子分子动力学（AIMD）的扩散过程模拟。
NEB方法通过分析可能的扩散路径上的静态构象，模拟单原子或两原子运动的势能面以推算扩散势垒（Ea）。其计算成本低，适用于分析单原子在低浓度下的扩散行为，但是不适宜于模拟较高离子浓度条件下的协同扩散行为。
基于AIMD的方法可以直接模拟扩散行为，对固体电解质材料中普遍存在的协同扩散行为进行模拟和解释，然而，由于AIMD计算成本高，目前使用AIMD计算模拟的体系大小一般不超过200个原子，模拟时长在100-200 ps左右，并且无法直接模拟室温下的扩散行为，只能通过在高温下模拟扩散过程再线性外推预测室温下的扩散性质。
以固体电解质材料为例，研究已经发现了三类在高温和低温具有不同扩散行为的材料体系。高温外推法只适用于不发生相变或准相变行为的体系。以实验手段相比照而言，目前常用的计算模拟方法是“非原位的”，不能实现对实验条件下电池材料的准确模拟。
简而言之，算的准的算不快（AIMD），算的快的算不准（NEB）。如何才能二者兼具呢？近年来快速发展的机器学习势函数方法为又快又准地模拟扩散行为提供了可能。
机器学习势函数方法通过满足一定几何和物理限制的变换，将每个原子的局域结构（local environment）表示为一个描述符（descriptor），再应用机器学习方法（如高斯过程，神经网络等）拟合其与高精度的第一性原理数据（能量，力等）间的关系，从而获得自适应的势函数（力场）。通过不断的从第一性原理数据中学习迭代提高势函数的精度和泛用性。当势函数迭代收敛后，可以应用于大体系长时间的准确模拟。当前适用于材料体系模拟的机器学习势函数方法主要有Belher-Parrinello Neural Network（BPNN）, Gaussian Approximation Potential（GAP）, aenet, Deep Potential（DP）等。
机器学习势函数方法已被成功应用于材料物理化学性质的模拟，实现千万原子体系纳秒尺度以上快速模拟，极大的拓展了原子尺度高精度模拟的边界。构造高精度的力场以往需要专业的研究人员进行长期开发优化，目前广泛使用的模拟软件如LAMMPS/Gromacs等等自建的力场都来源于专业人士的长期开发和积累。机器学习势函数方法，使得力场构建这一“旧时王谢堂前燕”，现在也能“飞入寻常百姓家”。
尽管如此，在实践中，要快速构造材料的势函数并实现高精度模拟，还需要关注：自动化的计算和数据处理流程。为了实现势函数的自动化构造，计算的所有步骤应当能自动化。由于材料体系往往含有多种元素，基于原子环境的机器学习势函数应当能处理多元素体系且计算成本可控。由于不同机器学习势函数构造方法和代码实现的区别，在计算多元素体系时，计算性能往往有数量级的差异。应用于实际模拟问题时，需要针对所研究问题建立一套高效、稳定的采样手段。保证势函数训练的稳定性。
得益于近年来计算材料学领域的迅猛发展，尤其是材料基因组项目的建设，推动了高效且标准化的开源材料数据处理工具的开发，使得材料计算模拟更加自动化。这些工具为机器学习势函数相关的代码进一步开发提供了坚实基础。目前，在机器学习势函数，基于Deep Potential 的Deep Potential GENerator （DP-GEN）软件包提供了自动训练势函数的开源解决方案。由于其基于End-to-End的神经网络架构，计算复杂度随元素种类数目接近线性增长，快于一般的机器学习势函数方法。因此，我们在研究中选用DP-Gen来构建无机固体电解质体系的势函数并模拟确立可靠的锂离子扩散过程研究protocol。
代码
文本
深度势能分子动力学模拟
在本教程中，我们将以电解质Li10GeP2S12为例，进行 LAMMPS 深度势能分子动力学模拟。
代码
文本
[3]
# 下载教程文件
! cd /data; if [ ! -e colombo-academy-tutorials ];then git clone https://gitee.com/deepmodeling/colombo-academy-tutorials.git;fi;
Cloning into 'colombo-academy-tutorials'... remote: Enumerating objects: 7295, done. remote: Counting objects: 100% (305/305), done. remote: Compressing objects: 100% (265/265), done. remote: Total 7295 (delta 88), reused 175 (delta 34), pack-reused 6990 Receiving objects: 100% (7295/7295), 76.07 MiB | 2.89 MiB/s, done. Resolving deltas: 100% (3388/3388), done. Updating files: 100% (299/299), done.
代码
文本
文件夹包含了400、1000 K下运行LAMMPS所需的输入文件；以400K为例——
代码
文本
[4]
! tree /data/colombo-academy-tutorials/DeePMD-SSE
/data/colombo-academy-tutorials/DeePMD-SSE ├── 1000K │ ├── data.lmp │ ├── dump.lammpstraj │ ├── input.lammps │ ├── log.lammps │ ├── model.pb -> ../400K/model.pb │ ├── target.msd │ └── target.rdf ├── 400K │ ├── data.lmp │ ├── input.lammps │ └── model.pb ├── input_dpa.json └── job_dpa.json 2 directories, 12 files
代码
文本
input.lammps: LAMMPS 输入文件，用于控制 LAMMPS MD 模拟的细节；data.lmp: 用于存放 MD 模拟的初始构型；model.pb: 深度势能模型。
代码
文本
（关于 LAMMPS 分子动力学模拟输入文件的更多信息，可以阅读「超详细「深度势能」材料计算上手指南｜章节 1」
其中只有两行例外：
pair_style deepmd model.pb
pair_coeff * *
其中调用了 DeePMD 的 pair_style，提供了模型文件 model.pb，这意味着原子间相互作用将由存储在文件model.pb中的 DP 模型进行计算。
在具有兼容版本的 LAMMPS 的环境中，可以通过以下命令执行深度势分子动力学模拟：
代码
文本
[5]
! cd /data/colombo-academy-tutorials/DeePMD-SSE/400K && lmp -i input.lammps
Warning: This LAMMPS executable is in a conda environment, but the environment has not been activated. Libraries may fail to load. To activate this environment please see https://conda.io/activation. LAMMPS (23 Jun 2022 - Update 1) OMP_NUM_THREADS environment is not set. Defaulting to 1 thread. (src/comm.cpp:98) using 1 OpenMP thread(s) per MPI task Loaded 1 plugins from /opt/deepmd-kit-2.2.1/lib/deepmd_lmp Reading data file ... orthogonal box = (0 0 0) to (26.67021 26.67021 25.5908) 1 by 1 by 1 MPI processor grid reading atoms ... 900 atoms read_data CPU = 0.008 seconds 360 atoms in group Li 36 atoms in group Ge 72 atoms in group P 432 atoms in group S DeePMD-kit WARNING: Environmental variable TF_INTRA_OP_PARALLELISM_THREADS is not set. Tune TF_INTRA_OP_PARALLELISM_THREADS for the best performance. See https://deepmd.rtfd.io/parallelism/ for more information. DeePMD-kit WARNING: Environmental variable TF_INTER_OP_PARALLELISM_THREADS is not set. Tune TF_INTER_OP_PARALLELISM_THREADS for the best performance. See https://deepmd.rtfd.io/parallelism/ for more information. DeePMD-kit WARNING: Environmental variable OMP_NUM_THREADS is not set. Tune OMP_NUM_THREADS for the best performance. See https://deepmd.rtfd.io/parallelism/ for more information. Summary of lammps deepmd module ... >>> Info of deepmd-kit: installed to: /opt/deepmd-kit-2.2.1 source: v2.2.1 source branch: HEAD source commit: 3ac8c4c7 source commit at: 2023-03-16 12:33:24 +0800 surpport model ver.:1.1 build variant: cuda build with tf inc: /opt/deepmd-kit-2.2.1/include;/opt/deepmd-kit-2.2.1/include build with tf lib: /opt/deepmd-kit-2.2.1/lib/libtensorflow_cc.so set tf intra_op_parallelism_threads: 0 set tf inter_op_parallelism_threads: 0 >>> Info of lammps module: use deepmd-kit at: /opt/deepmd-kit-2.2.1DeePMD-kit WARNING: Environmental variable TF_INTRA_OP_PARALLELISM_THREADS is not set. Tune TF_INTRA_OP_PARALLELISM_THREADS for the best performance. See https://deepmd.rtfd.io/parallelism/ for more information. DeePMD-kit WARNING: Environmental variable TF_INTER_OP_PARALLELISM_THREADS is not set. Tune TF_INTER_OP_PARALLELISM_THREADS for the best performance. See https://deepmd.rtfd.io/parallelism/ for more information. DeePMD-kit WARNING: Environmental variable OMP_NUM_THREADS is not set. Tune OMP_NUM_THREADS for the best performance. See https://deepmd.rtfd.io/parallelism/ for more information. DeePMD-kit: Successfully load libcudart.so 2023-05-20 09:41:01.680872: I tensorflow/core/platform/cpu_feature_guard.cc:193] This TensorFlow binary is optimized with oneAPI Deep Neural Network Library (oneDNN) to use the following CPU instructions in performance-critical operations: SSE4.1 SSE4.2 AVX AVX2 AVX512F FMA To enable them in other operations, rebuild TensorFlow with the appropriate compiler flags. 2023-05-20 09:41:01.692097: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:975] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero 2023-05-20 09:41:01.853637: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:975] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero 2023-05-20 09:41:01.853899: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:975] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero 2023-05-20 09:41:03.497077: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:975] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero 2023-05-20 09:41:03.498086: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:975] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero 2023-05-20 09:41:03.498245: I tensorflow/stream_executor/cuda/cuda_gpu_executor.cc:975] successful NUMA node read from SysFS had negative value (-1), but there must be at least one NUMA node, so returning NUMA node zero 2023-05-20 09:41:03.498394: I tensorflow/core/common_runtime/gpu/gpu_device.cc:1532] Created device /job:localhost/replica:0/task:0/device:GPU:0 with 9910 MB memory: -> device: 0, name: NVIDIA GeForce RTX 2080 Ti, pci bus id: 0000:00:09.0, compute capability: 7.5 2023-05-20 09:41:03.502802: I tensorflow/core/common_runtime/process_util.cc:146] Creating new thread pool with default inter op setting: 2. Tune using inter_op_parallelism_threads for best performance. 2023-05-20 09:41:03.568084: I tensorflow/compiler/mlir/mlir_graph_optimization_pass.cc:354] MLIR V1 optimization pass is not enabled >>> Info of model(s): using 1 model(s): model.pb rcut in model: 6 ntypes in model: 4 using compute id: CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE Your simulation uses code contributions which should be cited: - USER-DEEPMD package: The log file lists these citations in BibTeX format. CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE-CITE Generated 0 of 6 mixed pair_coeff terms from geometric mixing rule Neighbor list info ... update every 10 steps, delay 0 steps, check no max neighbors/atom: 2000, page size: 100000 master list distance cutoff = 7 ghost atom cutoff = 7 binsize = 3.5, bins = 8 8 8 1 neighbor lists, perpetual/occasional/extra = 1 0 0 (1) pair deepmd, perpetual attributes: full, newton on pair build: full/bin/atomonly stencil: full/bin/3d bin: standard Setting up Verlet run ... Unit style : metal Current step : 0 Time step : 0.002 Per MPI rank memory allocation (min/avg/max) = 2.64 | 2.64 | 2.64 Mbytes Step Temp E_pair E_mol TotEng Press 0 400 -3890.7484 0 -3844.2665 -404.48187 2000 414.96465 -3848.873 0 -3800.6521 1643.2348 Loop time of 87.6519 on 1 procs for 2000 steps with 900 atoms Performance: 3.943 ns/day, 6.087 hours/ns, 22.818 timesteps/s 69.3% CPU use with 1 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- Pair | 86.873 | 86.873 | 86.873 | 0.0 | 99.11 Neigh | 0.69058 | 0.69058 | 0.69058 | 0.0 | 0.79 Comm | 0.039376 | 0.039376 | 0.039376 | 0.0 | 0.04 Output | 4.4727e-05 | 4.4727e-05 | 4.4727e-05 | 0.0 | 0.00 Modify | 0.041993 | 0.041993 | 0.041993 | 0.0 | 0.05 Other | | 0.006809 | | | 0.01 Nlocal: 900 ave 900 max 900 min Histogram: 1 0 0 0 0 0 0 0 0 0 Nghost: 2352 ave 2352 max 2352 min Histogram: 1 0 0 0 0 0 0 0 0 0 Neighs: 0 ave 0 max 0 min Histogram: 1 0 0 0 0 0 0 0 0 0 FullNghs: 62510 ave 62510 max 62510 min Histogram: 1 0 0 0 0 0 0 0 0 0 Total # of neighbors = 62510 Ave neighs/atom = 69.455556 Neighbor list builds = 200 Dangerous builds not checked Generated 0 of 6 mixed pair_coeff terms from geometric mixing rule Neighbor list info ... update every 10 steps, delay 0 steps, check no max neighbors/atom: 2000, page size: 100000 master list distance cutoff = 7 ghost atom cutoff = 7 binsize = 3.5, bins = 8 8 8 2 neighbor lists, perpetual/occasional/extra = 1 1 0 (1) pair deepmd, perpetual attributes: full, newton on pair build: full/bin/atomonly stencil: full/bin/3d bin: standard (2) compute rdf, occasional, half/full from (1) attributes: half, newton on pair build: halffull/newton stencil: none bin: none Setting up Verlet run ... Unit style : metal Current step : 0 Time step : 0.002 Per MPI rank memory allocation (min/avg/max) = 6.056 | 6.056 | 6.056 Mbytes Step Time Temp KinEng TotEng Press 0 0 414.96465 48.220913 -3800.6521 1643.2348 100 0.2 396.73006 46.101965 -3802.0645 2453.828 200 0.4 404.00717 46.9476 -3804.4317 4622.9644 300 0.6 390.60969 45.390748 -3806.081 2473.0423 400 0.8 395.53443 45.963027 -3804.9305 1420.8348 500 1 403.84346 46.928577 -3803.8313 3186.2472 600 1.2 403.74749 46.917425 -3803.1558 2818.9726 700 1.4 419.57274 48.756395 -3801.7147 455.40811 800 1.6 409.33217 47.566392 -3801.4872 3850.8294 900 1.8 406.26266 47.209699 -3803.6426 3610.9555 1000 2 389.29442 45.237908 -3805.8068 658.51648 1100 2.2 384.33531 44.661635 -3805.5094 1264.6869 1200 2.4 402.31087 46.750482 -3803.7749 553.25616 1300 2.6 400.70158 46.563476 -3802.8745 1171.1034 1400 2.8 409.6187 47.599688 -3802.8675 2041.7628 1500 3 408.73988 47.497564 -3803.4222 1567.0507 1600 3.2 401.86362 46.698509 -3804.6914 2901.2827 1700 3.4 398.27429 46.281413 -3807.8708 1367.0198 1800 3.6 385.30622 44.77446 -3807.9024 1184.2403 1900 3.8 389.92752 45.311476 -3805.7484 2931.0407 2000 4 390.07137 45.328193 -3804.4805 726.86634 2100 4.2 412.74724 47.96324 -3802.6925 1243.33 2200 4.4 415.74166 48.311206 -3802.0576 3892.4318 2300 4.6 406.21928 47.204659 -3804.391 1833.509 2400 4.8 402.90247 46.819229 -3806.4653 905.23213 2500 5 376.25566 43.72274 -3805.7702 3020.3616 2600 5.2 399.70304 46.44744 -3804.2739 2153.2272 2700 5.4 403.63034 46.903812 -3803.1738 1767.2716 2800 5.6 406.57819 47.246366 -3802.7978 3820.2456 2900 5.8 420.5481 48.869737 -3802.95 1187.1806 3000 6 384.86182 44.722817 -3806.6268 2003.964 3100 6.2 390.30772 45.355658 -3807.4654 2662.2586 3200 6.4 418.98526 48.688128 -3805.2818 862.44607 3300 6.6 415.91672 48.331548 -3803.1299 1554.3029 3400 6.8 400.58742 46.550209 -3802.8476 1986.4111 3500 7 397.41652 46.181735 -3803.4263 2123.9317 3600 7.2 405.49111 47.120042 -3805.2646 2335.4595 3700 7.4 384.21869 44.648083 -3807.321 1155.3073 3800 7.6 392.14723 45.569418 -3805.9379 1637.3393 3900 7.8 412.46583 47.930538 -3803.2359 272.98171 4000 8 397.6087 46.204067 -3803.1318 2276.0607 4100 8.2 403.46435 46.884523 -3802.966 1332.0367 4200 8.4 395.25951 45.93108 -3804.0161 3299.0096 4300 8.6 389.21547 45.228733 -3806.7812 1671.5144 4400 8.8 393.85005 45.767294 -3806.7571 1461.3168 4500 9 404.45485 46.999623 -3803.8941 2199.8231 4600 9.2 407.9495 47.405719 -3802.7351 754.66485 4700 9.4 401.79121 46.690095 -3803.2863 2687.137 4800 9.6 395.11217 45.913959 -3804.8329 1101.1495 4900 9.8 389.69226 45.284139 -3806.0182 3187.5813 5000 10 398.47502 46.304738 -3805.8715 596.19799 5100 10.2 410.79637 47.736539 -3803.2783 2159.5007 5200 10.4 417.32106 48.494739 -3802.1362 3377.6974 5300 10.6 411.6441 47.835049 -3803.4631 735.62014 5400 10.8 387.14209 44.987797 -3806.1587 2683.5597 5500 11 394.53241 45.846587 -3806.3986 2133.0402 5600 11.2 398.71573 46.332709 -3805.5098 2485.7628 5700 11.4 402.61657 46.786006 -3804.3232 -574.32672 5800 11.6 402.84929 46.81305 -3803.9773 2769.5945 5900 11.8 408.02082 47.414006 -3804.5261 2335.2829 6000 12 385.48871 44.795666 -3806.5365 1122.4845 6100 12.2 375.88852 43.680076 -3808.3234 1470.5527 6200 12.4 395.41164 45.948759 -3807.1676 2017.7038 6300 12.6 406.29806 47.213813 -3804.3073 -123.0101 6400 12.8 409.461 47.581363 -3803.8735 3059.5559 6500 13 398.91993 46.356439 -3804.72 2572.3015 6600 13.2 409.94773 47.637923 -3806.3028 1068.0786 6700 13.4 380.12652 44.172553 -3808.5205 2243.954 6800 13.6 392.40697 45.599602 -3808.4786 1823.4312 6900 13.8 400.2379 46.509594 -3804.8889 2743.4092 7000 14 401.02079 46.600569 -3804.0771 3243.1628 7100 14.2 406.67689 47.257835 -3804.4349 3005.8437 7200 14.4 403.01128 46.831873 -3805.6655 989.23654 7300 14.6 401.47637 46.653509 -3807.5742 2056.9532 7400 14.8 386.99829 44.971086 -3808.4332 1929.38 7500 15 395.58113 45.968454 -3805.1327 3345.352 7600 15.2 406.8446 47.277324 -3803.545 2676.5631 7700 15.4 405.84414 47.161066 -3804.0742 1235.6681 7800 15.6 399.91464 46.472029 -3805.3832 3000.7814 7900 15.8 377.35722 43.850747 -3806.9164 3115.0484 8000 16 388.83258 45.184239 -3808.0104 1439.1432 8100 16.2 393.06368 45.675913 -3807.1628 2536.8434 8200 16.4 414.01346 48.11038 -3804.0433 2235.4623 8300 16.6 389.85327 45.302849 -3803.547 3031.1758 8400 16.8 420.26629 48.83699 -3804.5391 1516.8576 8500 17 391.78824 45.527702 -3806.3987 1986.4862 8600 17.2 377.99314 43.924643 -3808.4465 1645.8156 8700 17.4 375.48519 43.633207 -3809.271 3536.4974 8800 17.6 372.1538 43.246084 -3805.9003 2794.7765 8900 17.8 400.85937 46.581812 -3804.1704 2151.4908 9000 18 404.67853 47.025616 -3803.6224 3205.4681 9100 18.2 409.55158 47.591888 -3804.232 1286.6804 9200 18.4 400.24761 46.510722 -3805.7512 1532.6519 9300 18.6 382.97512 44.503574 -3807.8793 945.14596 9400 18.8 385.09617 44.75005 -3807.3921 2210.5693 9500 19 391.36614 45.478652 -3803.8366 1876.0911 9600 19.2 389.49146 45.260805 -3803.2383 3463.6274 9700 19.4 404.73862 47.032599 -3803.4737 1750.2261 9800 19.6 399.45887 46.419066 -3804.6321 -55.941522 9900 19.8 374.73423 43.545943 -3807.0816 1992.0097 10000 20 379.70344 44.123389 -3806.9559 -159.78913 Loop time of 444.987 on 1 procs for 10000 steps with 900 atoms Performance: 3.883 ns/day, 6.180 hours/ns, 22.473 timesteps/s 68.9% CPU use with 1 MPI tasks x 1 OpenMP threads MPI task timing breakdown: Section | min time | avg time | max time |%varavg| %total --------------------------------------------------------------- Pair | 440.68 | 440.68 | 440.68 | 0.0 | 99.03 Neigh | 3.4454 | 3.4454 | 3.4454 | 0.0 | 0.77 Comm | 0.18997 | 0.18997 | 0.18997 | 0.0 | 0.04 Output | 0.16892 | 0.16892 | 0.16892 | 0.0 | 0.04 Modify | 0.47305 | 0.47305 | 0.47305 | 0.0 | 0.11 Other | | 0.03181 | | | 0.01 Nlocal: 900 ave 900 max 900 min Histogram: 1 0 0 0 0 0 0 0 0 0 Nghost: 2356 ave 2356 max 2356 min Histogram: 1 0 0 0 0 0 0 0 0 0 Neighs: 31292 ave 31292 max 31292 min Histogram: 1 0 0 0 0 0 0 0 0 0 FullNghs: 62584 ave 62584 max 62584 min Histogram: 1 0 0 0 0 0 0 0 0 0 Total # of neighbors = 62584 Ave neighs/atom = 69.537778 Neighbor list builds = 1000 Dangerous builds not checked Total wall time: 0:08:58
代码
文本
运行结束后，我们会得到如下输出文件：
dump.lammpstraj: LAMMPS 输出轨迹文件；log.lammps: LAMMPS 输出日志文件；target.msd: 径向分布函数（RDF）文件，记录了每 Nfrequency 步数输出的 RDF；target.rdf: 均方位移（MSD）文件，记录了模拟系统中离子的均方位移随时间的变化。
代码
文本
到这里，我们已经完成了在400 K下的深度势能动力学模拟；
本教程中的案例包中1000K文件夹下已给出了相应的输出文件，计算流程同上，只是修改了目标体系温度。
代码
文本
径向分布函数计算
下面，我们继续使用 Python 脚本进行对分子动力学模拟的结果进行分析，计算我们实际关心的物理性质。
代码
文本
[6]
# 以下命令判断该环境中是否存在所需模块，如果没有，则使用 pip 快速安装。
! if ! pip show matplotlib > /dev/null; then pip install matplotlib; fi;
WARNING: Package(s) not found: matplotlib Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple Collecting matplotlib Using cached https://pypi.tuna.tsinghua.edu.cn/packages/89/f3/84a9a6613ab0d89931d785f13fa2606e03f07252875acc8ebf5b676fa3c5/matplotlib-3.7.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (11.6 MB) Collecting fonttools>=4.22.0 Downloading https://pypi.tuna.tsinghua.edu.cn/packages/ad/5f/20da4f41e33e77723b0100ded6539529bd159319ed49d6459a4647cdc7ee/fonttools-4.39.4-py3-none-any.whl (1.0 MB) ━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━ 1.0/1.0 MB 2.4 MB/s eta 0:00:00a 0:00:01 Collecting cycler>=0.10 Using cached https://pypi.tuna.tsinghua.edu.cn/packages/5c/f9/695d6bedebd747e5eb0fe8fad57b72fdf25411273a39791cde838d5a8f51/cycler-0.11.0-py3-none-any.whl (6.4 kB) Requirement already satisfied: packaging>=20.0 in /opt/deepmd-kit-2.2.1/lib/python3.10/site-packages (from matplotlib) (22.0) Collecting kiwisolver>=1.0.1 Using cached https://pypi.tuna.tsinghua.edu.cn/packages/79/0f/5cc4ca3df66c49817944b9a1c7343ba70aefffc868ddf651d7839cc5dffd/kiwisolver-1.4.4-cp310-cp310-manylinux_2_12_x86_64.manylinux2010_x86_64.whl (1.6 MB) Collecting contourpy>=1.0.1 Using cached https://pypi.tuna.tsinghua.edu.cn/packages/ec/59/5eac40e348a7bf803cea221bcd27f74a49cb81667b400fdfbb680e86e7bb/contourpy-1.0.7-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (300 kB) Collecting pillow>=6.2.0 Using cached https://pypi.tuna.tsinghua.edu.cn/packages/25/6b/d3c35d207c9c0b6c2f855420f62e64ef43d348e8c797ad1c32b9f2106a19/Pillow-9.5.0-cp310-cp310-manylinux_2_28_x86_64.whl (3.4 MB) Collecting pyparsing>=2.3.1 Using cached https://pypi.tuna.tsinghua.edu.cn/packages/6c/10/a7d0fa5baea8fe7b50f448ab742f26f52b80bfca85ac2be9d35cdd9a3246/pyparsing-3.0.9-py3-none-any.whl (98 kB) Requirement already satisfied: python-dateutil>=2.7 in /opt/deepmd-kit-2.2.1/lib/python3.10/site-packages (from matplotlib) (2.8.2) Requirement already satisfied: numpy>=1.20 in /opt/deepmd-kit-2.2.1/lib/python3.10/site-packages (from matplotlib) (1.22.3) Requirement already satisfied: six>=1.5 in /opt/deepmd-kit-2.2.1/lib/python3.10/site-packages (from python-dateutil>=2.7->matplotlib) (1.16.0) Installing collected packages: pyparsing, pillow, kiwisolver, fonttools, cycler, contourpy, matplotlib Successfully installed contourpy-1.0.7 cycler-0.11.0 fonttools-4.39.4 kiwisolver-1.4.4 matplotlib-3.7.1 pillow-9.5.0 pyparsing-3.0.9 WARNING: Running pip as the 'root' user can result in broken permissions and conflicting behaviour with the system package manager. It is recommended to use a virtual environment instead: https://pip.pypa.io/warnings/venv
代码
文本
[7]
import os
os.chdir("/data/colombo-academy-tutorials/DeePMD-SSE")
代码
文本
[9]
import numpy as np
import matplotlib.pyplot as plt
nbins = 100 # define the number of bins in the RDF
T = ['400K', '1000K']
ave_rdf = []
labels = ['Li-Li', 'Li-Ge', 'Li-P', 'Li-S']
fig, axs = plt.subplots(1, 2, figsize=(10, 3))
for j in range(len(T)):
with open(T[j] + '/target.rdf', "r") as f: # read the target.rdf file
lines = f.readlines()
lines = lines[3:]
data = np.zeros((nbins, 9))
count = 0
for line in lines:
nums = line.split()
if len(nums) == 10:
for i in range(1, 10):
data[int(nums[0])-1, i-1] += float(nums[i]) # accumulatie data for each bin
if len(nums) == 2:
count += 1 # count the number of accumulations for each bin
ave_rdf.append(data / count) # calculate the averaged RDF data
np.savetxt( T[j] + '/ave_rdf.txt',ave_rdf[j])
for i, label in zip(range(1, 9, 2), labels):
axs[j].plot(ave_rdf[j][:, 0], ave_rdf[j][:, i], label=label)
axs[j].set_xlabel('r/Å')
axs[j].set_ylabel('g(r)')
axs[j].set_ylim([0,6])
axs[j].set_title('rdf({})'.format(T[j]))
axs[j].legend()
axs[j].plot()
plt.show()
# plt.savefig('rdf.png', dpi=300)
代码
文本
我们可以与文献结果中进行对比
代码
文本
扩散系数计算
代码
文本
[11]
import numpy as np
import matplotlib.pyplot as plt
T = ['400K', '1000K']
fig, axs = plt.subplots(1, 2, figsize=(10, 3))
for j in range(len(T)):
data = np.loadtxt(T[j] + '/target.msd', skiprows=2)
timestep = data[:, 0]
time = timestep*2 # input.lammps：t_step= 2fs
msd1 = data[:, 1]
msd2 = data[:, 2]
msd3 = data[:, 3]
msd4 = data[:, 4]
axs[j].plot(time/1000, msd1, label='Li$^+$') # 1fs= 1/1000ps
axs[j].plot(time/1000, msd2, label='Ge$^{4+}$')
axs[j].plot(time/1000, msd3, label='P$^{5+}$')
axs[j].plot(time/1000, msd4, label='S$^{2-}$')
axs[j].set_xlabel('timestep(ps)')
axs[j].set_ylabel('MSD(Å^2)')
slope1, residuals = np.polyfit(time, msd1, 1)
slope2, residuals = np.polyfit(time, msd2, 1)
slope3, residuals = np.polyfit(time, msd3, 1)
slope4, residuals = np.polyfit(time, msd4, 1)
Diff1 = slope1/6 * 1e-5 # D=1/6*slope; 1 A^2/fs= 1e-5 m^2/s
Diff2 = slope2/6 * 1e-5
Diff3 = slope3/6 * 1e-5
Diff4 = slope4/6 * 1e-5
print(f'temperature: {T[j]}')
print(f"Diffusion Coefficients of Li+: {Diff1} m^2/s")
print(f"Diffusion Coefficients of Ge4+: {Diff2} m^2/s")
print(f"Diffusion Coefficients of P5+: {Diff3} m^2/s")
print(f"Diffusion Coefficients of S2-: {Diff4} m^2/s")
#axs[j].set_ylim([0,300])
axs[j].set_title('msd({})'.format(T[j]))
axs[j].legend()
axs[j].plot()
plt.show()
# plt.savefig('msd.png', dpi=300)
temperature: 400K Diffusion Coefficients of Li+: 2.1051433527470397e-10 m^2/s Diffusion Coefficients of Ge4+: 1.142521034750528e-12 m^2/s Diffusion Coefficients of P5+: 6.579770627062671e-13 m^2/s Diffusion Coefficients of S2-: 1.6257687827606234e-12 m^2/s temperature: 1000K Diffusion Coefficients of Li+: 4.387296831245876e-09 m^2/s Diffusion Coefficients of Ge4+: 1.5165098488282157e-13 m^2/s Diffusion Coefficients of P5+: 1.2384106130369005e-13 m^2/s Diffusion Coefficients of S2-: 6.3535584734912696e-12 m^2/s
代码
文本
以及从均方位移导出的扩散系数的数值：
temperature: 400K
Diffusion Coefficients of Li+: 1.255409630651049e-10 m^2/s
Diffusion Coefficients of Ge4+: 5.933157939101913e-14 m^2/s
Diffusion Coefficients of P5+: 8.503725727827082e-14 m^2/s
Diffusion Coefficients of S2-: 4.7863835675262014e-14 m^2/s
temperature: 1000K
Diffusion Coefficients of Li+: 4.387296831245876e-09 m^2/s
Diffusion Coefficients of Ge4+: 1.5165098488282157e-13 m^2/s
Diffusion Coefficients of P5+: 1.2384106130369005e-13 m^2/s
Diffusion Coefficients of S2-: 6.3535584734912696e-12 m^2/s
对比可以看出：体系的Li+扩散速率最快，占据主导作用；且温度高时，扩散系数明显增大
代码
文本
文献中的计算结果见下，可以看出在在1000K下体系的扩散系数约为 ($\texttt{m^2/s}$), 在400K下体系的扩散系数约为 ($\texttt{m^2/s}$)。与我们的计算结果非常接近。
代码
文本
代码
文本
参考
- Han Wang, Linfeng Zhang, Jiequn Han, and Weinan E. DeePMD-kit: A deep learning package for many-body potential energy representation and molecular dynamics. Comput. Phys. Comm., 228:178–184, 2018. doi:10.1016/j.cpc.2018.03.016.
- Jinzhe Zeng, Duo Zhang, Denghui Lu, Pinghui Mo, Zeyu Li, Yixiao Chen, Marián Rynik, Li'ang Huang, Ziyao Li, Shaochen Shi, Yingze Wang, Haotian Ye, Ping Tuo, Jiabin Yang, Ye Ding, Yifan Li, Davide Tisi, Qiyu Zeng, Han Bao, Yu Xia, Jiameng Huang, Koki Muraoka, Yibo Wang, Junhan Chang, Fengbo Yuan, Sigbjørn Løland Bore, Chun Cai, Yinnian Lin, Bo Wang, Jiayan Xu, Jia-Xin Zhu, Chenxing Luo, Yuzhi Zhang, Rhys E. A. Goodall, Wenshuo Liang, Anurag Kumar Singh, Sikai Yao, Jingchao Zhang, Renata Wentzcovitch, Jiequn Han, Jie Liu, Weile Jia, Darrin M. York, Weinan E, Roberto Car, Linfeng Zhang, and Han Wang. DeePMD-kit v2: A software package for Deep Potential models. 2023. doi:10.48550/arXiv.2304.09409.
- Huang J, Zhang L, Wang H, Zhao J, Cheng J, E W. Deep potential generation scheme and simulation protocol for the Li10GeP2S12-type superionic conductors. J Chem Phys. 2021;154(9):094703. doi:10.1063/5.0041849
- https://docs.deepmodeling.com/projects/deepmd/en/master/index.html
- https://github.com/deepmodeling/deepmd-kit
代码
文本
已赞10
本文被以下合集收录
机器学习与DFT精华帖
gtang
更新于 2024-09-10
38 篇21 人关注
电池计算
微信用户70yQ
更新于 2024-09-06
8 篇1 人关注
推荐阅读
公开
基于图论的电解液溶剂分子生成算法AndyX
发布于 2023-10-27
4 赞2 转存文件
公开
CALYPSO—电池实战之界面结构预测篇 zhanglinshuang
发布于 2023-10-27
1 赞5 转存文件