计算化学公社

标题: 结晶过程晶体生长的分子动力学模拟求助 [打印本页]

作者
Author: syarnold 时间: 2025-5-7 09:45
标题: 结晶过程晶体生长的分子动力学模拟求助
本帖最后由 syarnold 于 2025-5-7 09:45 编辑

各位大佬，我想通过分子动力学来模拟一下结晶中晶体生长的过程，参考的是比较早也比较出名的一篇JACS的文章：
J. Am. Chem. Soc. 2005, 127, 6, 1975–1982

文献中是模拟尿素晶体在水中的溶解和生长过程，模型是两个尿素晶体的超胞中间夹着25埃厚的一层水，首先通过升温并维持的方法实现了晶体分子溶解到水中，然后通过降温并维持的方法实现了水溶液中的尿素分子重新长回到了晶体上。下面是一些文献中模拟的设置和模拟结果：下面第一个图为文献中的控压设置，第二个图为文献中的溶解过程，第三个图文献中的结晶过程。
(, 下载次数 Times of downloads: 11) (, 下载次数 Times of downloads: 17) (, 下载次数 Times of downloads: 16)

控压方法：
文献中控压方法是NPT的anisotropic，描述的是x和y方向是受晶格的影响，几乎不变；z方向受界面作用（溶液和结晶过程）的影响发生改变。但是由于这里我不清楚应该如何设置6个方向的控压参数，我尝试用半各向同性进行控压：
Pcoupl    = berendsen
pcoupltype             = semiisotropic
tau_p                = 1.0
ref_p                = 1.0    0
compressibility       = 0    4.5e-5;

问题：
1. 我想用我要研究的有机小分子的晶胞复现文中的结晶过程。晶胞是从CCDC下载cif后加H后通过cp2k优化晶胞，然后完全按照文献中的建模（两个超胞夹着溶剂层）和模拟方法（除了NPT控压方式之外几乎相同的设置）进行模拟，其中一个面的晶体层溶解过程是可以复现的（如下图左至右），但是尝试了多种方法（尝试不同面、不控压、各向同性和半各向同性的控压、各种温度和退火等设置）都没能实现溶液中分子向晶体分子转变的结晶过程（如下图右至左）。请各位大佬帮忙看下我的模型、力场、拓扑和mdp文件有无问题（附件是这部分的内容）。其他说明：原子电荷是通过Multiwfn计算的RESP2；力场是通过Sobtop软件来指认GAFF原子类型，没法指认的自动用UFF原子类型；晶胞含有水和Na离子，在模型和top文件中也都单独进行了定义。
(, 下载次数 Times of downloads: 16)
2. 考虑到此结构晶体是由有机小分子和Na离子同时按照一定规律排布构成的晶胞（有机小分子见下图），我又尝试了将Na加入到分子中一块构建拓扑文件（RESP2电荷也依此方法进行计算），但是这种方式在MD预平衡过程晶体方式发生明显了歪斜，见下图。请问各位老师，我问题1中难以模拟出结晶过程的原因是否可能是由分子中的Na离子和其他原子分开的缘故呢？如果确实可能是由于这个原因，那我将Na加入到分子中一块进行力场文件生成的方法是否可行？晶体分子发生歪斜的原因又是什么呢？

(, 下载次数 Times of downloads: 14) (, 下载次数 Times of downloads: 17)

作者
Author: sobereva 时间: 2025-5-7 15:30
水要用专门的水模型，不要用通用的拓扑文件产生工具产生。这在http://sobereva.com/soft/Sobtop/#FAQ20明确说了

你的prod.mdp里的控压设置和贴子里写的不同

晶体生长的速度相对于溶解过程慢得多，可能时间不够长。并且要恰当设置温度，特定力场下的物质的熔点和实验可能相差不少，温度设置要留有余量。并且先确保结晶过程的温度设置下，原有的晶体能始终稳定存在。

创建拓扑文件的时候绝对不要把Na+和阴离子部分作为整体。

作者
Author: syarnold 时间: 2025-5-19 14:34
本帖最后由 syarnold 于 2025-5-19 14:56 编辑

sobereva 发表于 2025-5-7 15:30
水要用专门的水模型，不要用通用的拓扑文件产生工具产生。这在http://sobereva.com/soft/Sobtop/#FAQ20明确 ...

首先非常感谢社长回复。

这段时间我根据社长在此帖以及公社群里的回复，针对我之前模拟的一些问题进行了修正，然而依旧无法MD模拟出溶液结晶过程（过饱和溶液中的溶质分子向固体状的晶体分子转变的过程）。进行的修改包括：修正水模型-改用了三点水模型；将Na+离子和主体分子分开进行了拓扑和力场创建；尝试了各种控压方式和设置（包括不控压NVT，isotropic控压以及如一楼所示的semiisotropic控压）；降低了体系温度（实际过程约为283K，模拟使用了263K）；尝试在预平衡过程中加入速度产生（gen_vel = yes; gen_temp = 298.15; gen_seed = -1）；增长了计算时间，最长的跑到了百纳秒还多，但是溶质分子仍未有明显向结晶晶面吸附生长的趋势；精确计算了盒子内溶剂和溶质分子的数目以确保为过饱和溶液；去掉了晶胞中的结晶水（原本为一水合物），重新对晶胞结构进行了优化，然而此模型产生了连之前可以溶解的过程都无法复现的不良结果；我想保持分子刚性，遂将限制设为constraints = all-bonds（后经公社群中社长回复为不可，实际测试了也确实没有任何作用）。

我猜测了一些这个问题可能的原因，还请社长和各位大佬看一下是否问题出在这里面：

1. 刚刚发现单胞中的分子是有两种构型（很细小的差别，见下图），我仔细查看了晶体库中的cif文件的模型，确实是不太一样。我打算将两种分子分别进行拓扑和力场的生成，不知道这样是否会有帮助。
(, 下载次数 Times of downloads: 17)
2. 此分子在晶体中的构型长宽比例较大，有点类似于条状。MD模拟中发现溶质分子的构型会发生变形，是否是因为这个原因导致很难按照晶体中的构型排列规律地吸附到晶体上去？

3. 由于Na+离子和溶质分子的力场是分开的，溶液中它们的分布也并不是十分统一（并不是1：1相邻），是否会因此导致溶液中的溶质分子和Na+离子想按照晶体中的排布方式“生长”上去会很困难（见下图）？
(, 下载次数 Times of downloads: 13)
4. 如果通过力场将分子设置为刚性分子是否对此问题有所帮助？

5. 根据社长使用Sobtop超级方便地创建二茂铁的GROMACS的拓扑文件帖子，将这个体系，按照全用谐振势的分子的拓扑来产生全刚性的分子拓扑文件是否可行？

6. 晶体中水分子的存在是否会影响MD模拟中溶质分子往晶体上的吸附和生长过程？我的思考是：这里若使得一个溶质分子按照晶胞排布规律吸附到晶面上，也需要一个Na+离子和一个水分子在特定位置被吸附，但整个溶液结晶生长的模拟过程变得概率很低也很困难。因为就之前模拟的观察发现，Na+离子和水分子并不会明显向晶体中的Na离子层或水分子层位置聚集，因此我在这里有疑问的是：是我力场拓扑等设置不当导致上述原因，还是这个过程本身（Na+离子或水分子被吸附到晶面特定区域）在MD模拟中就很难产生？

作者
Author: sobereva 时间: 2025-5-20 01:31

syarnold 发表于 2025-5-19 14:34
首先非常感谢社长回复。

这段时间我根据社长在此帖以及公社群里的回复，针对我之前模拟的一些问题进行 ...

1 不需要分别产生拓扑信息

2 溶剂中的构象和堆叠进入晶体后的构象差异大很正常，分子间相互作用情况明显不同。溶液中构象正常的话就不需要考虑这个

3 有可能导致晶体长得较慢，因为MD过程中分子和离子自发产生那种排列的概率较低

4、5 如果你是指设成和晶体中的构象一样的刚性结构，可以一试，但物理意义不清楚，就算跑出来了也可能被审稿人comment

6 可能在微观尺度来看就是很难产生。虽然你也可以尝试不同原子电荷、力场看看情况

作者
Author: syarnold 时间: 2025-5-20 15:45
本帖最后由 syarnold 于 2025-5-20 16:17 编辑

sobereva 发表于 2025-5-20 01:31
1 不需要分别产生拓扑信息

2 溶剂中的构象和堆叠进入晶体后的构象差异大很正常，分子间相互作用情况明 ...

感谢社长回复！

我想尝试一下sobtop用谐振势方法生成拓扑（二茂铁例子：使用Sobtop超级方便地创建二茂铁的GROMACS的拓扑文件），但是力场结果中的能量k全部为0：
(, 下载次数 Times of downloads: 12)

模型是从cif晶胞文件中扣出来的单个分子，加H后进行了限制性几何优化（冻结了除H之外的所有原子）和振动计算，方法和基组见下图（一开始用nosymm是想将结构优化后H的坐标导回到晶胞文件中，后面尝试了去掉nosymm，也并未解决问题），然后将优化后的结构去掉Na+离子后进行了RESP2电荷计算，将准备好的.mol2、.chg和.fchk按照示例中的方法（见下图）通过sobtop生成了拓扑和力场文件（结果见1-nona.itp和FFderiv.txt）。
(, 下载次数 Times of downloads: 16)

(, 下载次数 Times of downloads: 17)
(, 下载次数 Times of downloads: 16) (, 下载次数 Times of downloads: 13) (, 下载次数 Times of downloads: 14)

结果中原子类型和原子index部分的有几个乱码，应该是.fchk文件中有Na而.mol2模型文件中没有Na导致的（加上Na就没有乱码了）。主要问题是生成的力场中所有的k都是0，请问社长，这个问题可能是什么原因造成的？检查了振动也没有虚频：
(, 下载次数 Times of downloads: 16)

作者
Author: sobereva 时间: 2025-5-21 03:54

syarnold 发表于 2025-5-20 15:45
感谢社长回复！

我想尝试一下sobtop用谐振势方法生成拓扑（二茂铁例子：使用Sobtop超级方便地创建二茂 ...

我没时间仔细检查。对于维持刚性的目的，本质上就是用谐振势进行限制，都不需要基于Hessian来算严格的k，直接让sobtop直接猜力常数就够了

欢迎光临计算化学公社 (http://bbs.keinsci.com/)