请选择 进入手机版 | 继续访问电脑版

计算化学公社

 找回密码
 现在注册!
查看: 2240|回复: 15

[PSI4] 考察SAPT能量分解的能量项随分子二聚体间距变化的简单方法

[复制链接]

2万

帖子

25

威望

3万

eV
积分
56004

管理员

公社社长

发表于 2019-3-4 06:11:58 | 显示全部楼层 |阅读模式
注:本文内容主要是写给一个外国合作者的看的,所以是用英文写的。本文介绍的dimerscan程序也可以结合Gaussian非常容易地实现二聚体间距的刚性扫描,详见《详谈使用Gaussian做势能面扫描》(http://sobereva.com/474)。
对于单一结构做SAPT分析,在此文有详细介绍:《使用PSI4做对称匹配微扰理论(SAPT)能量分解计算》(http://sobereva.com/526)。



考察SAPT能量分解的能量项随分子二聚体间距变化的简单方法
A simple way to investigate the variation of SAPT energy decomposition terms with respect to molecular dimer separation


文/Sobereva@北京科音  2019-Mar-4


1 Foreward

The symmetry adapted perturbation theory (SAPT) is a rigorous and useful way of evaluating intermolecular interaction energy; more importantly, the SAPT interaction energy can be decomposed into various physical terms, which are very useful for understanding the nature of the interaction. PSI4 should be the best code for performing the SAPT calulation, because it is free, fast and relatively easy to use. PSI4 can be obtained from http://psicode.org. I assume readers have already known how to use the PSI4 and meantime PSI4 has been installed on the machine in correct way.

It is often interesting to examine variation of SAPT interaction terms with respect to molecular dimer separation. In very simple cases, such as scanning SAPT terms for Ar-Ar distance of Ar dimer, this can be done by properly making use the rule of PSI4 input file, see this file for example: http://sobereva.com/attach/469/Ar-Ar.inp. However, for most practical cases, e.g. scanning the O...H distance in H2O...HOH dimer, this purpose cannot be easily realized by only using single PSI4 input file. In present article, I will introduce a simple and general way of obtaining all SAPT terms with varying intermolecular distance. Specifically, the "distance" mentioned here denotes spacing of a pair of atoms between the two monomers, and in the variation process the orientation and structure of the monomers keep fixed. The method described in this article consists of three steps:
(1) Generating SAPT input files
(2) Running PSI4 in batch mode
(3) Extracting data and plotting them as curve map

Below I will take water dimer as instance to show how to do. We want to vary H2-O4 distance from 1.5 to 3.5 Angstrom with step size of 0.2, therefore there should be 10 steps.
1.png

All relevant files of the water dimer instance can be downloaded here: http://sobereva.com/attach/469/waterdimer.rar

Note that the method illustrated in this article has been employed in my very systematic research on hydrogen bond, I strongly suggest reading it:
Exploring Nature and Predicting Strength of Hydrogen Bonds: A Correlation Analysis Between Atoms-in-Molecules Descriptors, Binding Energies, and Energy Components of Symmetry-Adapted Perturbation Theory, J. Comput. Chem., 40, 2868–2881 (2019) DOI: 10.1002/jcc.26068


2 Generating SAPT input files

I wrote a code named dimerscan, which is used to generate PSI4 SAPT input files at various intermolecular separations. This code can be downloaded from this page: http://sobereva.com/soft/dimerscan.
PS: Please kindly cite dimerscan like this if dimerscan is utilized in your work: Tian Lu, dimerscan program, http://sobereva.com/soft/dimerscan (accessed month day, year)

In this package, the "dimerscan.exe" is executable file in Windows, the "dimerscan" is executable file in Linux, and "SAPT_template.inp" is template input file of PSI4 SAPT task, it will be loaded by dimerscan. You can manually change memory specification, basis set and keyword of SAPT task in this file, while other parts should keep unchanged without special reasons.
PS: In the template file, the SAPT level is SAPT2+(3)dmp2/aug-cc-pVTZ. According to the benchmark in J. Chem. Phys., 140, 094106 (2014), this level is able to predict binding energy in very high accuracy (usually close to CCSD(T)/aug-cc-pVTZ), therefore it is highly recommended for small dimers.

We first optimize structure of water dimer at reasonable level (e.g. B3LYP-D3(BJ)/def2-TZVP) via your favourite quantum chemistry code, then create a plain text file like below (this is content of the H2Odimer.txt in the waterdimer.rar package):

3 3
0 1 0 1
O                  1.21571500    0.00004000   -0.12596300
H                  0.26191400   -0.00018400    0.04416700
H                  1.62937300   -0.00025200    0.74021400
O                 -1.08643200   -0.00004600    0.11097100
H                 -1.46302100   -0.76572600   -0.33250500
H                 -1.46252900    0.76620800   -0.33193600


The first line is the number of atoms in monomer 1 and monomer 2, the second line is net charge and spin multiplicity of monomer 1 and monomer 2. The content starting from line 3 is dimer coordinate (in Angstrom).

Put the SAPT_template.inp file into the same folder of dimerscan executable file, then boot up dimerscan, input path of the H2Odimer.txt, and input:
2,4  // Index of the atom in monomer 1 (H2) and that in monomer 2 (O4). Their distance will be gradually changed in your specified way
1.5  // Initial distance of H2-O4
10   // Number of steps
0.2  // Stepsize

Now you will find totally 11 .inp files have been generated in current folder. The 0000.inp is PSI4 SAPT input file at initial distance, while 0001.inp, 0002.inp ... 0010.inp correspond to the structure of scan steps 1, 2 ... 10. You can also find a file named scan.xyz in current folder, it is a trajectory file containing all the 11 geometries. You can load it into VMD visualization program to check if the geometry variation meets our expectation. The animation of playing the scan.xyz is shown below, as you can clearly see the 11 geometries have been correctly generated.


2.gif

As an example, the content of the generated 0002.inp is shown below. This file is easy to undertstand, the task use memory up to 2000MB, current H2-O4 distance is 1.9, as denoted by the "dimerdist" variable. After completing the SAPT2+(3)dmp2 energy evaluating, various terms are assigned to user-defined variables, and finally they are printed to the output file together in specific format.

memory 2000 MB

molecule dimer {
  0  1
O      1.215715    0.000040   -0.125963
H      0.261914   -0.000184    0.044167
H      1.629373   -0.000252    0.740214
--
  0  1
O     -1.635758    0.000010    0.138187
H     -2.012347   -0.765670   -0.305289
H     -2.011855    0.766264   -0.304720
}
dimerdist=    1.900000

set {
basis aug-cc-pVTZ
scf_type DF
freeze_core True
}

energy('sapt2+(3)dmp2')
E_disp = get_variable('SAPT DISP ENERGY') * psi_hartree2kcalmol
E_elst = get_variable('SAPT ELST ENERGY') * psi_hartree2kcalmol
E_exch = get_variable('SAPT EXCH ENERGY') * psi_hartree2kcalmol
E_ind = get_variable('SAPT IND ENERGY') * psi_hartree2kcalmol
E_tot = get_variable('SAPT TOTAL ENERGY') * psi_hartree2kcalmol

psi4.print_out("\n")
psi4.print_out(" Summary of SAPT result (kcal/mol)\n")
psi4.print_out(" Distance      E_tot     E_elst     E_exch     E_disp      E_ind\n")
psi4.print_out("%s %6.3f %10.3f %10.3f %10.3f %10.3f %10.3f\n" % ("R=",dimerdist,E_tot,E_elst,E_exch,E_disp,E_ind))


3 Running the SAPT input files in batches

Now, put all the .inp files we just generated to a folder used for running PSI4 calculation, and then put this shell script into the this folder: http://sobereva.com/attach/469/psi4all.sh. Open this file by text editor, change the value after -n argument to actual number of CPU cores that you want to use to conduct PSI4 calculation. Next, enter this folder and run this command: ./psi4all.sh. This script will invoke "psi4" command to execute all .inp files in current folder in turn and produce output files with same name but contain .out suffix.


4 Extracting data and plotting them as curve map

After all SAPT calculations have finished, run this command to extract data from output files to result.txt file: grep "R= " *.out > result.txt

The content of result.txt is shown below. As can be seen from the .inp files, the values should be: distance, total binding energy, electrostatic term, exchange term, dispersion term and induction term.
0000.out:R=  1.500      0.536    -20.989     37.925     -6.180    -10.220
0001.out:R=  1.700     -3.769    -13.378     19.103     -4.034     -5.460
0002.out:R=  1.900     -4.923     -8.860      9.532     -2.660     -2.935
0003.out:R=  2.100     -4.798     -6.144      4.719     -1.771     -1.602
0004.out:R=  2.300     -4.223     -4.460      2.322     -1.189     -0.896
0005.out:R=  2.500     -3.552     -3.370      1.138     -0.805     -0.515
0006.out:R=  2.700     -2.932     -2.630      0.556     -0.551     -0.306
0007.out:R=  2.900     -2.407     -2.108      0.271     -0.381     -0.189
0008.out:R=  3.100     -1.979     -1.723      0.132     -0.267     -0.120
0009.out:R=  3.300     -1.636     -1.431      0.064     -0.190     -0.079
0010.out:R=  3.500     -1.363     -1.203      0.031     -0.137     -0.054

Next we plot the data as curve map via Origin program. Boot up Origin, directly dragging the result.txt into the Origin window to import it, properly specifying column names, and then plotting Line+Symbol map of various SAPT terms with respect to the distance values, you will eventually get below map.
3.png
From the map it can be seen that the minimum point of water dimer should occur around the position of R(H2-O4)=1.9 Angstrom. Electrostatic, dispersion and induction terms play attractive role, while repulsive effect purely comes from exchange term. All terms smoothly converge to zero with respect to increase of H2-O4 separation. The data also clearly shows that the H-bond in water dimer is completely dominated by electrostatic interaction irrespective of the dimer separation.

评分

参与人数 7eV +34 收起 理由
asdf + 5 精品内容
柒月小鱼 + 5 你太可爱!!!!!!!!!
978142355 + 5 GJ!
alonewolfyang + 5 好物!
zsu007 + 5 赞!
lpqfnu + 4 谢谢分享
liyuanhe211 + 5 好物!

查看全部评分

第七届北京科音分子动力学与GROMACS培训班 预计于8月于北京举办,内容见http://www.keinsci.com/workshop/KGMX_content.html,欢迎关注。报名时间将在培训开始前一个月于北京科音官网通知。
北京科音自然科学研究中心http://www.keinsci.com  致力于计算化学的发展和传播,长期开办最高水准的各种量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训,是提升计算化学研究水平的最佳选择。欢迎加入“北京科音”公众号获取培训最新消息和计算化学资讯!培训相关信息见《北京科音办的培训班FAQ》(http://bbs.keinsci.com/thread-5098-1-1.html)。
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群:思想家公社QQ群1号:18616395,2号:466017436。合计约6000人。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一定会被拒绝加入。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

227

帖子

1

威望

3854

eV
积分
4101

Level 6 (一方通行)

发表于 2019-3-4 07:47:31 | 显示全部楼层
谢谢社长的分享!

2470

帖子

10

威望

5285

eV
积分
7955

Level 6 (一方通行)

首席卖萌官

发表于 2019-3-4 08:06:29 | 显示全部楼层
why not upload it onto the Multiwfn forum?
She doesn't love me.
Even so,
my heart has been taken away by her.

2万

帖子

25

威望

3万

eV
积分
56004

管理员

公社社长

 楼主| 发表于 2019-3-4 16:19:54 | 显示全部楼层
我本是个娃娃 发表于 2019-3-4 08:06
why not upload it onto the Multiwfn forum?

Irrelevant to Multiwfn
第七届北京科音分子动力学与GROMACS培训班 预计于8月于北京举办,内容见http://www.keinsci.com/workshop/KGMX_content.html,欢迎关注。报名时间将在培训开始前一个月于北京科音官网通知。
北京科音自然科学研究中心http://www.keinsci.com  致力于计算化学的发展和传播,长期开办最高水准的各种量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训,是提升计算化学研究水平的最佳选择。欢迎加入“北京科音”公众号获取培训最新消息和计算化学资讯!培训相关信息见《北京科音办的培训班FAQ》(http://bbs.keinsci.com/thread-5098-1-1.html)。
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群:思想家公社QQ群1号:18616395,2号:466017436。合计约6000人。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一定会被拒绝加入。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

151

帖子

0

威望

503

eV
积分
654

Level 4 (黑子)

蓝卫兵

发表于 2019-3-22 09:22:27 | 显示全部楼层
Screen Shot 2019-03-21 at 8.21.44 PM.png
B样条插值
个人专栏https://zhuanlan.zhihu.com/p/21936803

527

帖子

8

威望

6423

eV
积分
7110

Level 6 (一方通行)

发表于 2019-8-21 16:06:25 | 显示全部楼层
请问大神
运行程序时 psi4 -i *.dat -o *.dat -n
这个 -n 指的应该是线程数而不是核数吧?
计算化学与分子模拟

2万

帖子

25

威望

3万

eV
积分
56004

管理员

公社社长

 楼主| 发表于 2019-8-22 00:28:19 | 显示全部楼层
captain 发表于 2019-8-21 16:06
请问大神
运行程序时 psi4 -i *.dat -o *.dat -n
这个 -n 指的应该是线程数而不是核数吧?

是。但从实际来讲二者不作区分,设N线程,对于高并行度的阶段,平均利用的CPU核数也相当于N
第七届北京科音分子动力学与GROMACS培训班 预计于8月于北京举办,内容见http://www.keinsci.com/workshop/KGMX_content.html,欢迎关注。报名时间将在培训开始前一个月于北京科音官网通知。
北京科音自然科学研究中心http://www.keinsci.com  致力于计算化学的发展和传播,长期开办最高水准的各种量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训,是提升计算化学研究水平的最佳选择。欢迎加入“北京科音”公众号获取培训最新消息和计算化学资讯!培训相关信息见《北京科音办的培训班FAQ》(http://bbs.keinsci.com/thread-5098-1-1.html)。
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群:思想家公社QQ群1号:18616395,2号:466017436。合计约6000人。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一定会被拒绝加入。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

527

帖子

8

威望

6423

eV
积分
7110

Level 6 (一方通行)

发表于 2019-8-22 07:34:07 | 显示全部楼层
sobereva 发表于 2019-8-22 00:28
是。但从实际来讲二者不作区分,设N线程,对于高并行度的阶段,平均利用的CPU核数也相当于N

明白了 谢大神!
计算化学与分子模拟

527

帖子

8

威望

6423

eV
积分
7110

Level 6 (一方通行)

发表于 2019-12-22 18:40:13 | 显示全部楼层
请问大神
每次SAPT计算后,输出文件和结果看起来都是正常的。
但是终端提示:
Using `psi4.core.get_variable` instead of `psi4.core.variable` (or `psi4.core.scalar_variable` for scalar variables only) is deprecated, and in 1.4 it will stop working
这个提示应该没有什么影响吧?
计算化学与分子模拟

2万

帖子

25

威望

3万

eV
积分
56004

管理员

公社社长

 楼主| 发表于 2019-12-22 19:38:28 | 显示全部楼层
captain 发表于 2019-12-22 18:40
请问大神
每次SAPT计算后,输出文件和结果看起来都是正常的。
但是终端提示:

没事,不用管
第七届北京科音分子动力学与GROMACS培训班 预计于8月于北京举办,内容见http://www.keinsci.com/workshop/KGMX_content.html,欢迎关注。报名时间将在培训开始前一个月于北京科音官网通知。
北京科音自然科学研究中心http://www.keinsci.com  致力于计算化学的发展和传播,长期开办最高水准的各种量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训,是提升计算化学研究水平的最佳选择。欢迎加入“北京科音”公众号获取培训最新消息和计算化学资讯!培训相关信息见《北京科音办的培训班FAQ》(http://bbs.keinsci.com/thread-5098-1-1.html)。
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群:思想家公社QQ群1号:18616395,2号:466017436。合计约6000人。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一定会被拒绝加入。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

527

帖子

8

威望

6423

eV
积分
7110

Level 6 (一方通行)

发表于 2019-12-22 19:56:55 | 显示全部楼层

好的 谢大神!
计算化学与分子模拟

49

帖子

0

威望

307

eV
积分
356

Level 3 能力者

发表于 2020-4-1 06:54:16 | 显示全部楼层
Pleeese, should we follow letterly the tutorial in youtube to install Psi4? Some steps are incomprehensible. We need your assistance Professor Tian about the installation process.

Capture.JPG


2万

帖子

25

威望

3万

eV
积分
56004

管理员

公社社长

 楼主| 发表于 2020-4-4 05:35:06 | 显示全部楼层
zako 发表于 2020-4-1 06:54
Pleeese, should we follow letterly the tutorial in youtube to install Psi4? Some steps are incompreh ...

Steps are described in this article
使用PSI4做对称匹配微扰理论(SAPT)能量分解计算
http://sobereva.com/526http://bbs.keinsci.com/thread-15902-1-1.html
第七届北京科音分子动力学与GROMACS培训班 预计于8月于北京举办,内容见http://www.keinsci.com/workshop/KGMX_content.html,欢迎关注。报名时间将在培训开始前一个月于北京科音官网通知。
北京科音自然科学研究中心http://www.keinsci.com  致力于计算化学的发展和传播,长期开办最高水准的各种量子化学、分子动力学、波函数分析与Multiwfn程序等主题的培训,是提升计算化学研究水平的最佳选择。欢迎加入“北京科音”公众号获取培训最新消息和计算化学资讯!培训相关信息见《北京科音办的培训班FAQ》(http://bbs.keinsci.com/thread-5098-1-1.html)。
欢迎加入人气最高、水准最高的综合性理论与计算化学交流QQ群:思想家公社QQ群1号:18616395,2号:466017436。合计约6000人。两个群讨论范畴相同,可加入任意其一但不可都加入,申请信息必须注明具体研究方向,否则一定会被拒绝加入。
思想家公社的门口Blog:http://sobereva.com(发布大量原创计算化学相关博文)
Multiwfn主页:http://sobereva.com/multiwfn(最强大的量子化学波函数分析程序)
ResearchGate:https://www.researchgate.net/profile/Tian_Lu
Money and papers are rubbish, get a real life!

49

帖子

0

威望

307

eV
积分
356

Level 3 能力者

发表于 2020-4-10 20:23:18 | 显示全部楼层
about the conda package, should i install the miniconda, anaconda, which one ?

446

帖子

1

威望

893

eV
积分
1359

Level 4 (黑子)

发表于 2020-4-10 21:33:59 | 显示全部楼层
zako 发表于 2020-4-10 20:23
about the conda package, should i install the miniconda, anaconda, which one ?

Actually you do not need conda indeed. Just make sure you have numpy (maybe a little others) package of python is properly installed. For Debian/Ubuntu/Linux Mint users you can directly use "sudo apt install psi4" to install psi4. For Fedora/Cent OS/Red Hat users I have not tested, but you can try "sudo yum install psi4".
您需要登录后才可以回帖 登录 | 现在注册!

本版积分规则

手机版|北京科音自然科学研究中心|京公网安备 11010502035419号|计算化学公社 — 北京科音旗下高水平计算化学交流论坛 ( 京ICP备14038949-1号 )

GMT+8, 2020-7-8 13:21 , Processed in 0.172723 second(s), 29 queries .

快速回复 返回顶部 返回列表