计算化学公社

标题: 考察SAPT能量分解的能量项随分子二聚体间距变化的简单方法 [打印本页]

作者
Author: sobereva 时间: 2019-3-4 06:11
标题: 考察SAPT能量分解的能量项随分子二聚体间距变化的简单方法
注：如今非常推荐用sobEDA方法做能量项随二聚体间距离变化的扫描，sobEDA比SAPT普适得多，而且速度快得多，内存消耗少得多，结合我给的脚本做扫描极为方便。具体做法在《使用sobEDA和sobEDAw方法做非常准确、快速、方便、普适的能量分解分析》（http://sobereva.com/685）里提到的sobEDA原文附带的教程里给出了。

注：本文内容主要是写给一个外国合作者的看的，所以是用英文写的。本文介绍的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.

(, 下载次数 Times of downloads: 148)

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.

(, 下载次数 Times of downloads: 153)

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.

(, 下载次数 Times of downloads: 151)

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.

作者
Author: zsu007 时间: 2019-3-4 07:47
谢谢社长的分享！

作者
Author: 我本是个娃娃 时间: 2019-3-4 08:06
why not upload it onto the Multiwfn forum?

作者
Author: sobereva 时间: 2019-3-4 16:19

我本是个娃娃发表于 2019-3-4 08:06
why not upload it onto the Multiwfn forum?

Irrelevant to Multiwfn

作者
Author: pyscf 时间: 2019-3-22 09:22
(, 下载次数 Times of downloads: 152)

作者
Author: captain 时间: 2019-8-21 16:06
请问大神
运行程序时 psi4 -i *.dat -o *.dat -n
这个 -n 指的应该是线程数而不是核数吧？

作者
Author: sobereva 时间: 2019-8-22 00:28

captain 发表于 2019-8-21 16:06
请问大神
运行程序时 psi4 -i *.dat -o *.dat -n
这个 -n 指的应该是线程数而不是核数吧？

是。但从实际来讲二者不作区分，设N线程，对于高并行度的阶段，平均利用的CPU核数也相当于N

作者
Author: captain 时间: 2019-8-22 07:34

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

明白了谢大神！

作者
Author: captain 时间: 2019-12-22 18:40
请问大神
每次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
这个提示应该没有什么影响吧？

作者
Author: sobereva 时间: 2019-12-22 19:38

captain 发表于 2019-12-22 18:40
请问大神
每次SAPT计算后，输出文件和结果看起来都是正常的。
但是终端提示：

没事，不用管

作者
Author: captain 时间: 2019-12-22 19:56

sobereva 发表于 2019-12-22 19:38
没事，不用管

好的谢大神！

作者
Author: zako 时间: 2020-4-1 06:54
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.

(, 下载次数 Times of downloads: 106)

作者
Author: sobereva 时间: 2020-4-4 05:35

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/526（http://bbs.keinsci.com/thread-15902-1-1.html）

作者
Author: zako 时间: 2020-4-10 20:23
about the conda package, should i install the miniconda, anaconda, which one ?

作者
Author: snljty 时间: 2020-4-10 21:33

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".

作者
Author: zako 时间: 2020-4-19 04:48
Dear snljty,
Thank you very very much.

作者
Author: mekelleche 时间: 2021-1-10 16:15
本帖最后由 mekelleche 于 2021-1-10 16:44 编辑

Is it possible to analyse a transition state of an organic reaction (example : cycloaddition) by sapt using Psi4 ?

作者
Author: sobereva 时间: 2021-1-11 02:37

mekelleche 发表于 2021-1-10 16:15
Is it possible to analyse a transition state of an organic reaction (example : cycloaddition) by sap ...

SAPT assumes that interfragment interaction is relatively weak (since it is based on perturbation theory), the two fragments should not interact with each other via chemical bonding (very strong interaction), hence SAPT is not well-suited for analyzing TS.

作者
Author: aitouna 时间: 2021-1-11 05:36
本帖最后由 aitouna 于 2021-1-11 06:12 编辑

sobereva 发表于 2021-1-11 02:37
SAPT assumes that interfragment interaction is relatively weak (since it is based on perturbation ...

Dear mekelleche, thank you for asking that question, and thank you very much Professor Tian LU for your convincing answer, could you please nominate some suited tools to analyse the transition state by energy decomposition.

作者
Author: sobereva 时间: 2021-1-12 07:30

aitouna 发表于 2021-1-11 05:36
Dear mekelleche, thank you for asking that question, and thank you very much Professor Tian LU for ...

Try LMOEDA in GAMESS-US

作者
Author: lqhenwunai 时间: 2021-2-5 11:55
请问sapt0里面的
Disp20 (SS)       = -0.001622876335 [Eh]
Disp20 (OS)       = -0.001622876335 [Eh]
Exch-Disp20       =    0.000540619581 [Eh]
Exch-Disp20 (SS) =    0.000295146561 [Eh]
Exch-Disp20 (OS) =    0.000245473020 [Eh]

SS, OS分别代表的是什么？same spin 和opposite spin么？
谢谢。

作者
Author: sobereva 时间: 2021-2-6 03:45

lqhenwunai 发表于 2021-2-5 11:55
请问sapt0里面的
Disp20 (SS) = -0.001622876335 [Eh]
Disp20 (OS) = -0. ...

是

作者
Author: Xiaol 时间: 2024-10-23 14:15
请问dimerscan能否生成考察A+B二聚体SAPT能量分解的能量项随A分子甲基转动的变化的输入文件

作者
Author: sobereva 时间: 2024-10-24 12:25

Xiaol 发表于 2024-10-23 14:15
请问dimerscan能否生成考察A+B二聚体SAPT能量分解的能量项随A分子甲基转动的变化的输入文件

不能

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