分子吸附材料表面，可能的构象太多，求解！

MrD

研究分子吸附一个材料表面问题，表面是周期性的，分子是不规则的，因此可能的吸附形态太多。

而几何优化的时候，分子不会大幅度翻转或平移，只是在给定初始结构附近，在Z方向上有变化，这样就只能人为做很多种初始构型。

请问有没有更好的方法或算法，能让分子在表面驰豫幅度大一些，尽量找到一个能量最小点。

尝试过第一性原理的动力学，但由于体系大，80多原子，即使降低到最差精度，一天也只能几百步（每步一飞秒），4核8GB的win7机器表示很骄傲。

另外想请教，如果想做第一原理动力学的话，想对这种结构跑平衡，恒温器有没有太大影响，或者使用恒温器应根据什么样的经验来选择呢？时间尺度多少事宜呢？是否还要模拟退火？

我不太懂动力学，感觉在跑平衡的过程中，势能应该是降低然后趋于平衡，但我跑出的势能为什么是波浪线呢？（我用的NVT系宗，300K，simple NH恒温器）


PS：我用的是MS-Dmol, 有人说dmol优化周期性体系不给力，让转VASP，但我始终觉得只要算法完善，不应该有程序不给力的一说。请各位指点。

卡开发发

1、先回答问题。我只说说自己的看法吧：
体系80原子似乎也不算很复杂的体系，不妨把结构发来看看，但我也没有十足的把握。如果构象实在太多，而且彼此的能量也很接近，这种情况优化分子的意义就并不是很大了（各个构型的能量离平均值很近，出现的概率都很高）。由于涉及到成键断键，基于力场化的动力学得不到好的结论，从头算动力学还是挺有必要的。

要是AIMD计算速度比较慢并且非常难于处理的话，不妨试试DFTB+（程序可能要申请才能拿到）来跑PBC分子动力学，速度可能快得多，parameterize可以挑几个可能的或是接近的结构通过DMol3计算来fitting程序的hopping integral。

如果对编程比较熟悉的话，理论上说可以把DFTB+的动力学和DMol3通过Sob前辈编写的MolClus程序（http://www.keinsci.com/research/molclus.html）对接起来（大概相当于DFTB+跑动力学并计算每个点采样并对每一帧的结构进行优化，Molclus对各个构象进行统计排序筛选，之后把挑出的能量比较低的几个结构作为DMol3初始结构来进行优化），但是interface这部分的代码得自己编写，能做出来也算是相当有特色的工作了。不过不知道我理解Molclus的使用方式对不对了@sobereva

恒温器方面的问题我不太懂，好像陈敏伯老师那本《分子模拟》的书上可能有介绍，不过可以问问@ChemiAndy。不过如果势能曲线一直在振荡的话可能应该跑的时长还不够长，体系还没达到平衡。

2、回PS，我对这个问题还是抱有一定的感慨的，非得吐个槽不可：

你说得对，只有算法或方法是否完善，没有程序给不给力的说法。毕竟没有一个程序是万能的，根据体系选择合适的程序还是好的想法。DMol3确实方法不太完善，当然与其文件格式等封闭有关，比如同样是数值轨道方法，SIESTA的基函数便允许自行生成更高角动量的极化函数，而DMol3始终只能采用DNP级别的；比如至今DMol3的杂化泛函、MetaGGA非常慢（难以忍受得慢），而且至今也无法支持周期体系，而且杂化泛函只能做B3LYP，这样就算有TDDFT用处也不大；无法进行DFT+U，处理一些d态为价带的体系甚是棘手，尤其发文章容易被打回来。其他更高级的方法就不提了。。。分析方法由于波函数文件tpdensk等都是加密的，目前来说也只能对grd文件做点简单的分析，MS图形界面下没有的比如Laplacian能够自己写程序凑合做出来，但ELF就没戏了。

这篇帖子http://bbs.keinsci.com/forum.php ... 87&highlight=Grimme应该能够说明DMol3结果还是可靠吧，VASP算出来0.7eV我用DMol3算出来0.9eV，何况文献中二维结构还没进行优化并且真空层可能不充分（要记得VASP的Grimme校正Cutoff是30A）。DMol3计算表面很大的优势在于速度快并且对晶格常数大小依赖非常小，这个特性足以取非常充分的真空层。VASP也并非那么刀枪不入，有些强关联磁性体系非得LAPW不可，我遇到过这样的体系赝势平面波给不出好结果的时候。

很多人可能认为DMol3只是图形界面下点一点就能给个结果，其实并不然，看了很多人做计算有时候cutoff和multipolar expansion设置得都不太地道。还有work funtion的校正是通过施加外电场进行的，如果你同时自己还额外指定了外电场，抱歉，指定的这部分的外电场完全不起作用。因此，计算给不给力并不在于使用了什么程序，好的结果都是对基于keywords和计算的原理熟悉的基础上，知道哪些东西可靠放心使用而哪些东西不可靠不要去使用，这样就行了。

实际上以前对MS我也有不小的偏见，刚接触的时候DMol3算一些固体表面，SCF各种不收敛，甚至还跟老师吐槽过DMol3，但现在DMol3是我使用频度最高的程序，当然VASP也仍然是一款很出色晶体结构模拟程序。
这点心得（可能也算不上什么有价值的）我也愿意与@jiewei 兄一块分享一下。

sobereva

可以先用molclus或其它分子构象搜索程序（比如MS里的conformer）找出一批能量较低的构象，写个脚本或程序，把分子构象搜索程序产生的那一大批低能构象依次自动摆到界面上方并自动调用dmol3、VASP之类的优化（摆上去的时候还应考虑不同朝向）。

对于通过跑MD产生molclus所用的最初结构，用分子力场就行了，如果你倾向于用MS，用forcite或discovery模块也都能跑，之后按照http://sobereva.com/143把轨迹导出成.xyz文件就能用于molclus了。

jiewei

卡开发发 发表于 2015-9-2 00:46
1、先回答问题。我只说说自己的看法吧：
体系80原子似乎也不算很复杂的体系，不妨把结构发来看看，但我也 ...

先mark一下，有时间来细看。

MrD

卡开发发 发表于 2015-9-2 00:46
1、先回答问题。我只说说自己的看法吧：
体系80原子似乎也不算很复杂的体系，不妨把结构发来看看，但我也 ...

首先，非常感谢您的回复，这么长的回复，啥也不说了，满满的感动。我以前做分子量化计算的，现在转行做第一原理，着实有太多的东西需要学习，感谢您建议的书籍和资料，以及经验分享。

1. 结构就是单层MoS2 （Z轴30埃）吸附Benzene，想把这个做好之后再做别的吸附，可以加您qq么（论坛上好像不能发送附件）。因为S表面实际上体现出一定的负电荷，而Bz的H带有少量正电荷，但C与S表面有色散作用，因此哪个作用多一点是会影响几何构象的。静电吸引强，Bz会立起来，色散强，Bz就躺下，像pai-stacking一样。这样一来，基于立场的动力学，就严重依赖初始电荷的分配了，我尝试过MS中的Adsorption Locator，在分配了电荷的情况下，Bz会立起来，但通过目前的dmol计算来看(考虑了6种平行，12种垂直，没有考虑倾斜)，Bz应该还是平行于底面的。组里有会做gromacs的，那个更需要电荷了。我一会儿尝试下Sob老师说的forcite。

2. DFTB以前用来计算分子过，非常好用，也特别喜欢这个程序。但没有S和Mo的sk文件参数啊。

3. Cutoff的设定着实没有经验。我目前设定是的Fine级别下的4.9，看文献有有设定4.6的，也不知是他打错了还是真调低了。4.9应该可以的吧。多级展开目前我用的是8级，这个多级展开主要对哪些方面有影响呢，还请不吝赐教。


PS1. 根据以前做分子量化时候的经验，dmol的速度确实堪称世界上最快的了，也只有dmol配得上我现在的机器了，先坚持用dmol，多掌握些理论。

PS2. 小眼神很销魂@卡开发发 ，一看就是大神。

MrD

sobereva 发表于 2015-9-2 06:20
可以先用molclus或其它分子构象搜索程序（比如MS里的conformer）找出一批能量较低的构象，写个脚本或程序， ...

非常感谢Sob老师，做计算差不多有5年了，一直都跟着Sob老师的博文学习，受益良多。

结构就是单层MoS2 （Z轴30埃）吸附Benzene。因为S表面实际上体现出一定的负电荷，而Bz的H带有少量正电荷，但C与S表面有色散作用，因此哪个作用多一点是会影响几何构象的。静电吸引强，Bz会立起来，色散强，Bz就躺下，像pai-stacking一样。这样一来，基于立场的动力学，就严重依赖初始电荷的分配了，我尝试过MS中的Adsorption Locator，在分配了电荷的情况下，Bz会立起来，但通过目前的dmol计算来看(考虑了6种平行，12种垂直，没有考虑倾斜)，Bz应该还是平行于底面的。组里有会做gromacs的，那个更需要电荷了。

Forcite，discover在设定的时候好像不涉及电荷，但是在结构上是可以分配电荷的呀，不知道分配电荷与否是不是有大的影响，我正在尝试。

zhangj158

sobereva 发表于 2015-9-2 06:20
可以先用molclus或其它分子构象搜索程序（比如MS里的conformer）找出一批能量较低的构象，写个脚本或程序， ...

老师，那存不存在分子构象搜索的时候1和2两个构象的能量差很大，1的能量更低，但是在吸附以后发现构象2这种情况吸附的总能量更低呢

卡开发发

MrD 发表于 2015-9-2 10:46
首先，非常感谢您的回复，这么长的回复，啥也不说了，满满的感动。我以前做分子量化计算的，现在转行做第 ...

1、DFTB+如果fitting了DMol3的结果作为SK文件的话可能靠谱点，个人觉得力场有时候不太可靠；
2、严格来说，你的体系应该造一个S、Mo、C、H比较小的体系，然后这些元素都得进行fitting才能够给出对于你的体系比较可靠的结果；
3、cutoff对应的精度是fine为0.1eV/atom，medium位0.3，coarse为1，一般fine就够了，但是注意，A+B=C的体系，A、B、C必须全局截断的数值相同而不是都取fine，你可以取体系中原子的cutoff最大的作为参考标准；

回PS1：DMol3做分子量化反而没啥优势，GGA对于分子体系普遍认可度不高，现在都杂化泛函；
PS2不知道怎么回，绝不是大神，只是老菜鸟而已，更期望做有价值的讨论。

BTW，QQ还是私下给吧。

sobereva

MrD 发表于 2015-9-2 10:46
首先，非常感谢您的回复，这么长的回复，啥也不说了，满满的感动。我以前做分子量化计算的，现在转行做第 ...

最快的绝对不是dmol3，而是ORCA开RI，见
大体系弱相互作用计算的解决之道
http://sobereva.com/214

sobereva

zhangj158 发表于 2015-9-2 16:18
老师，那存不存在分子构象搜索的时候1和2两个构象的能量差很大，1的能量更低，但是在吸附以后发现构象2这 ...

很可能这样，所以不能光取构象搜索时能量最低的那个，而是最低的一批都应当考虑，看谁吸附后能量最低。

sobereva

MrD 发表于 2015-9-2 11:08
非常感谢Sob老师，做计算差不多有5年了，一直都跟着Sob老师的博文学习，受益良多。

结构就是单层MoS2  ...

我说的基于力场的动力学是指对于被吸附的分子自身做动力学，然后再去摆到表面上通过量化来优化，所以不会有这个问题。

即便是力场下的表面+吸附分子的动力学，力场也只是提供一批初猜而已，最终优化时的描述还是量化级别的，因此也没这个问题。

jiewei

首先，你要明白一点， MoS2也算是一种衬底，分子在衬底上的吸附本身就有很多种可能， 包括平躺（一般计算都选的是这种形式，为什么？ 作用原子多作用力大啊），斜靠，和直立。 没中作用方式在实验上（也就是很多分子吸附在衬底上）都是有实验依据的。 所以构型多不是问题，问题是你要研究什么作用方式？ 因为我上面说的三种作用方式都是每一种都会带来不同的性质，有些利于电荷传输有些利于发光，有些利于。。。。
所以，能量最大原理这些都是在这种建模过程中 理论化学家 不懂实验的结果，只要明白了 你想做哪个放向，就建立那一种类型的模型。 这样构象是否就少了？ 这样是否就没有上面那些问题了？

另外再说一些，做理论别脱离了实验，尤其是计算化学，当然如果你是纯粹的理论化学可能要好一点， 不过理论最终目的是解释实验现象或者知道实验进行，就算是开辟新的理论也是为了后续的实验发现做准备。
最后，请记住上面一句话：理论是为实验服务的，模型是根据实验来的。

卡开发发

jiewei 发表于 2015-9-4 09:52
首先，你要明白一点， MoS2也算是一种衬底，分子在衬底上的吸附本身就有很多种可能， 包括平躺（一般计算都 ...

能量最大原理这些都是在这种建模过程中 理论化学家 不懂实验的结果

这个并不赞同，不知道算不算是以讹传讹，之所以做结构优化要求能量最低还是有原因的，并非理论家不了解。按照统计理论，低温条件下能量越低出现的概率越高，此时能量较低的构型对体系的贡献比较大。由于我们大多情况关注的是平均值，以及考虑到能量较低的情形贡献较大，所以对能量较低的结构的研究就更为关注。从这点出发，对于高温情形以及构象多且能量接近，优化也就完全无意义了http://sobereva.com/165。比较好的方式是，把几个能量较低的结构以及关注的（具有特殊性质的）结构挑出来参考。

后面的

理论是为实验服务的，模型是根据实验来的

我赞同，模型不对模拟就是乱做了，要依据实验依据原理作为指导来构建模型选取方法，这样的工作才有意义。

jiewei

卡开发发 发表于 2015-9-4 12:09
这个并不赞同，不知道算不算是以讹传讹，之所以做结构优化要求能量最低还是有原因的，并非理论家不了解 ...

可能是我表达有误， 能量最低肯定是有原因， 不过在某些条件下的分子构象未必是以能量最低的形态存在。 比如有些能长单晶的分子， 在旋涂的时候就很难会成为纯粹的单晶。而是会以amorphous的状态存在。
此时，如果我们要建立其中每种相与衬底的作用，则要建立分子和衬底的不同接触方式。根据能量最低原理（量化优化的时候）在某些接触方式时 一 优化结构确实就去了能量低的接触构型了，然而这个构型 却不能反应真实的实验条件，即某些分子和衬底接触不是以能量最低的方式接触的

卡开发发

jiewei 发表于 2015-9-4 13:48
可能是我表达有误， 能量最低肯定是有原因， 不过在某些条件下的分子构象未必是以能量最低的形态存在。  ...

一般来说应该很多情况都是polycrystal，如果是amorphous计算就真的没办法按照表面去做了。polycrystal就可以按照上次的想法，实际上每个晶粒及各个暴露面还是很大的，足以使用Slab模型进行近似。原则来说确实像兄台所说的，甚至不同的面的不同的吸附方式也得考虑。amorphous有时候用一些disorder处理（比如VCA）还凑合，但也很难说可不可靠。

当然，确实有不少相当复杂的情况其实理论非常难于研究，理论不是万能的。顺带吐槽一个：有时候这样的体系老师不了解计算还想让学生做，学生自己也没判断能力来判断这样的体系能不能通过ab initio的方式实现，这种情况非常让人头疼。