¶ 氯化钠晶体的模拟
¶ 引入
本文主要复刻物理学报上的一篇文献的核心模拟部分
崔守鑫, 蔡灵仓, 胡海泉, 郭永新, 向士凯, 经福谦. 氯化钠晶体在高温高压下热物理参数的分子动力学计算. 物理学报, 2005, 54(6): 2826-2831. doi: 10.7498/aps.54.2826
¶ 建模
原文使用力场较为复杂，本文只讨论基础MD，并非为了完整重现结果，因此力场使用Xponge中TIP3P水模型中附带的钠离子和氯离子力场
In Suk Joung and Thomas E. Cheatham, Determination of Alkali and Halide Monovalent Ion Parameters for Use in Explicitly Solvated Biomolecular Simulations, The Journal of Physical Chemistry B 2008 112 (30), 9020-9041, DOI: 10.1021/jp8001614
参照原文献
使用Xponge建模，编写python脚本build.py如下
#导入Xponge模块
import Xponge
#引入力场
import Xponge.forcefield.amber.tip3p
#定义结构基元，(0,0,0)的NA和(4,0,0)的CL
#此处4埃取值较为任意，不要太近即可，MD模拟能将其平衡至正确位置
mol0 = NA + CL
mol0.residues[-1].CL.x = 4
#定义从(0,0,0)到(64,64,64)的正方体填充区域
region = Xponge.BlockRegion(0, 0, 0, 64, 64, 64)
#定义从(0,0,0)到(70,70,70)的周期性盒子
box = Xponge.BlockRegion(0, 0, 0, 65, 65, 65)
#定义点阵，单位点阵间距8埃的fcc晶体
lattice = Xponge.Lattice("fcc", mol0, 8)
#创建元胞
mol = lattice.Create(box, region)
#保存
Save_PDB(mol, "NaCl.pdb")
Save_SPONGE_Input(mol, "NaCl")
使用命令行执行
python build.py
vmd NaCl.pdb
调整部分可视化设置，获得可视化结果：
¶ 动力学模拟
¶ 最小化
在命令行中输入以下命令，进行5000步的梯度下降最小化
SPONGE -default_in_file_prefix NaCl -mode minimization -step_limit 5000 -rst min
各部分的含义：
SPONGE
CUDA SPONGE的程序-default_in_file_prefix NaCl
输入文件的默认前缀为NaCl，对应于建模过程中Save_SPONGE_Input(mol, “NaCl”)-mode minimization
运行最小化，默认采用变步长的梯度下降算法-step_limit 5000
总步数5000步-rst min
生成的坐标和速度文件的前缀是min，也即生成的坐标文件为min_coordinate.txt，生成的速度文件为min_velocity.txt（对最小化而言速度文件无明显意义）
¶ 300 K, 40 GPa下的NPT模拟
在命令行中输入以下命令，进行50000步的NPT系综模拟
SPONGE -default_in_file_prefix NaCl -coordinate_in_file min_coordinate.txt -mode NPT -dt 1e-3 -step_limit 55000 -thermostat middle_langevin -middle_langevin_gamma 10 -barostat berendsen_barostat -target_temperature 300 -target_pressure 400000 -mdout 40GPa.out -rst 40GPa -crd 40GPa.dat -box 40GPa.box
各部分的含义：
SPONGE
CUDA SPONGE的程序-default_in_file_prefix NaCl
输入文件的默认前缀为NaCl，对应于建模过程中Save_SPONGE_Input(mol, “NaCl”)-coordinate_in_file min_coordinate.txt
使用最小化的结果的坐标作为初始坐标进行模拟-mode NPT
运行NPT模拟-dt 1e-3
模拟的步长1e-3皮秒，也即1飞秒-step_limit 5000
总步数50000步-thermostat middle_langevin
使用middle langevin方法控温-middle_langevin_gamma 10
控温的频率为10 ps^-1。注意此处默认值为1，而本模拟中压强过大，使用默认值频率太低可能无法控制住温度-barostat berendsen_barostat
使用berendsen方法控压-target_temperature 300
控温至300开-target_pressure 400000
控压至400000 bar（40 GPa）-mdout 40GPa.out
记录文件打印至40GPa.out
¶ 300 K, 0 GPa下的NPT模拟
在命令行中输入以下命令，进行50000步的NPT系综模拟。此处以40 GPa下获得的结果作为输入条件
SPONGE -default_in_file_prefix NaCl -coordinate_in_file 40GPa_coordinate.txt -velocity_in_file 40GPa_velocity.txt -mode NPT -dt 1e-3 -step_limit 55000 -thermostat middle_langevin -middle_langevin_gamma 10 -barostat berendsen_barostat -target_temperature 300 -target_pressure 1 -mdout 0GPa.out -crd 0GPa.dat -box 0GPa.box
各部分的含义同上，只额外添加了命令-velocity_in_file 40GPa_velocity.txt，添加了初始速度。
¶ 分析与讨论
¶ 轨迹可视化
使用VMD进行可视化，以40GPa的轨迹为例
vmd -sponge_mass NaCl_mass.txt 40GPa.dat
¶ 结果分析
编写分析脚本analysis.py
import numpy as np
import matplotlib.pyplot as plt
import Xponge
from Xponge.analysis import MdoutReader
mdout1 = MdoutReader("0GPa.out")
mdout2 = MdoutReader("40GPa.out")
plt.plot(mdout1.time, mdout1.density, label="0 GPa")
plt.plot(mdout1.time, mdout2.density, label="40 GPa")
plt.legend()
plt.show()
density1 = np.mean(mdout1.density)
density2 = np.mean(mdout2.density)
print(density1, density2, density1/density2)
获得结果分别为：
密度平均值分别为与，密度相对值为。
¶ 讨论
¶ 密度绝对值
氯化钠是常见晶体，我们熟悉其密度常温常压下约为，而本次模拟中获得的结果为
¶ 密度相对值
原文献中主要包含体积变化曲线，对比40 GPa和0 GPa下的数据，约为的倍数，本次模拟数值为偏高
¶ 分析
本次模拟结果均存在一定误差，因为本教程只是粗略模拟，而非科研文章。误差主要来源应来自力场，本文使用的力场主要对常温常压进行拟合，在高温高压场景下误差较大。