molpro中SAPT能量分解求助

jobs2016 · 发表于 Post on 2016-12-5 11:10:15

初次做能量分解，用DFT-SAPT方法。体系是由两个同种金属有机配合物及其反荷离子组成的二聚体。输入文件是照着molpro手册中SAPT能量分解给出的例子写的，例子见http://www.molpro.net/info/current/examples/near_sapt_acdft.com
因为例子中是两个惰性气体间的能量分解，每一个单体的原子就只有一个，而我的体系中每一个单体就有56个原子，所以我参照例子写的输入文件如下
gthresh,energy=1.d-8,orbital=1.d-8,grid=1.d-8
symmetry,nosym
orient,noorient
geometry={

28,Pt                0.07863700 0.44148700 -1.57899800
29,C                1.54895400 1.62862800 -1.37370800
30,N                2.40894600 2.41398200 -1.24533900
31,N                -1.45142600 -0.80653600 -1.87320100
32,N                -1.57370600 1.77195400 -1.47222200
33,C                3.46777700 3.26477400 -1.00224100
34,C                3.18847400 4.58211700 -0.59707800
35,C                4.76859400 2.74908000 -1.13110700
36,C                4.27605000 5.38533600 -0.26057800
37,C                5.82143400 3.59251400 -0.77016200
38,C                5.57890900 4.89125500 -0.33277800
39,H                4.09658200 6.41299200 0.06154600
40,H                6.84026500 3.20663400 -0.83506200
41,H                6.41647000 5.53391900 -0.05529800
42,C                5.00159400 1.35682200 -1.60667900
43,H                4.74421700 0.62402000 -0.82806000
44,H                4.39409600 1.13148600 -2.49567000
45,H                6.05640700 1.17896400 -1.83551000
46,C                1.77938200 5.08412700 -0.53339800
47,H                1.32140900 5.10025700 -1.53596400
48,H                1.14906700 4.44177700 0.10202200
49,H                1.74421800 6.10484400 -0.13133800
50,P                -7.24025400 1.29594500 -0.00374700
51,F                -6.80621500 2.39245800 -1.16923900
52,F                -7.92020800 2.47510800 0.88155600
53,F                -5.77374200 1.57341500 0.74861200
54,F                -6.47263800 0.09836200 -0.90477800
55,F                -7.60491500 0.17247300 1.13585200
56,F                -8.62668700 0.99040500 -0.79139000

84,Pt             -0.07870200 -0.44153500 1.57881100
85,C                -1.54901400 -1.62870800 1.37344100
86,N                -2.40895700 -2.41409900 1.24499100
87,N                1.45136400 0.80651100 1.87299500
88,N                1.57372100 -1.77194700 1.47204600
89,C                -3.46768400 -3.26499100 1.00191700
90,C                -3.18827000 -4.58230700 0.59676200
91,C                -4.76854800 -2.74941800 1.13079600
92,C                -4.27577300 -5.38563200 0.26029200
93,C                -5.82131900 -3.59295900 0.76988800
94,C                -5.57868200 -4.89168500 0.33250800
95,H                -4.09621800 -6.41326800 -0.06184500
96,H                -6.84019200 -3.20718200 0.83482000
97,H                -6.41618400 -5.53443200 0.05503400
98,C                -5.00168200 -1.35716600 1.60632600
99,H                -4.39438900 -1.13183000 2.49545800
100,H                -4.74414300 -0.62435000 0.82778200
101,H                -6.05654000 -1.17933500 1.83493800
102,C                -1.77915100 -5.08422000 0.53306400
103,H                -1.14878800 -4.44159800 -0.10202900
104,H                -1.32129700 -5.10072200 1.53567700
105,H                -1.74388100 -6.10477500 0.13060300
106,P                7.24034400 -1.29556800 0.00422400
107,F                7.60484200 -0.17181500 -1.13516200
108,F                5.77393500 -1.57319800 -0.74822400
109,F                6.47248700 -0.09822600 0.90539200
110,F                6.80640400 -2.39229200 1.16954600
111,F                7.92056000 -2.47446400 -0.88122100
112,F                8.62668500 -0.98990300 0.79197600
}

basis={
spd,pt,def2-tzvp,c;
spd,p,vdz,c;
spd,f,vdz,c;
spd,c,vdz,c;
spd,n,vdz,c;
sp,h,vdz,c;
}

!=========delta(HF) contribution for higher order interaction terms====
!sapt files
ca=2101.2
cb=2102.2

!dimer
hf
edm=energy

!monomer A
dummy,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112
{hf; save,$ca}
ema=energy
{sapt;monomerA}

!monomer Bdummy,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56
{hf; start,atdens; save,$cb}
emb=energy
{sapt;monomerB}

!interaction contributions
{sapt,SAPT_LEVEL=3;intermol,ca=$ca,cb=$cb,icpks=1}

!HF supermolecular interaction energy and delta(HF) contribution
eint_hf=(edm-ema-emb)*1000 mH
delta_hf=eint_hf-e1pol-e1ex-e2ind-e2exind

!=========DFT-SAPT at second order intermol. perturbation theory====
!sapt files
ca=2103.2
cb=2104.2

!shifts for asymptotic correction to xc potentialeps_homo_pbe0_57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112=-9.16    !HOMO(Ar)/PBE0 functionaleps_homo_pbe0_1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56=-9.16    !HOMO(Ne)/PBE0ip_57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112=-9.16=10.34                   !exp. ionisation potentialip_1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56=10.34                   !exp. ionisation potentialshift_57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112=ip_57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112+eps_homo_pbe0_57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112 !shift for bulk xc potential (Ar)shift_1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56=ip_1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56+eps_homo_pbe0_1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56 !shift for bulk xc potential (Ne)

!monomer Adummy,57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112{ks,pbe0; asymp,shift_1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56; save,$ca}
{sapt;monomerA}

!monomer Bdummy,1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19,20,21,22,23,24,25,26,27,28,29,30,31,32,33,34,35,36,37,38,39,40,41,42,43,44,45,46,47,48,49,50,51,52,53,54,55,56{ks,pbe0; start,atdens; asymp,shift_57,58,59,60,61,62,63,64,65,66,67,68,69,70,71,72,73,74,75,76,77,78,79,80,81,82,83,84,85,86,87,88,89,90,91,92,93,94,95,96,97,98,99,100,101,102,103,104,105,106,107,108,109,110,111,112; save,$cb}{sapt;monomerB}
!interaction contributions
{sapt;intermol,ca=$ca,cb=$cb,icpks=0}
!add high-order approximation to obtain the total interaction energy
eint_dftsapt=e12tot+delta_hf

几何坐标部分在这里我没有写全，请大神帮我看看怎样写才是正确！！！万分感谢！

sobereva · 发表于 Post on 2016-12-6 00:40:06

这么大体系想做DFT-SAPT基本没戏，算不动
建议你用GAMESS-US做LMOEDA能量分解，这种规模还算能算得动
Gaussian+Multiwfn也可以做简单的能量分解，看multiwfn手册4.100.8节的例子。

jobs2016 · 发表于 Post on 2016-12-6 10:27:01

sobereva 发表于 2016-12-6 00:40
这么大体系想做DFT-SAPT基本没戏，算不动
建议你用GAMESS-US做LMOEDA能量分解，这种规模还算能算得动
Gau ...

谢谢sob老师，已经在研究用Gaussian+Multiwfn能量分解了。既然molpro中DFT-SAPT对于该体系没戏，那我就放弃这种方法。

daogege · 发表于 Post on 2017-6-21 20:36:26

请问，这里面的交换相关渐进式的矫正是必须的吗？没有实验值的话，可以用计算值吗

sobereva · 发表于 Post on 2017-6-22 01:57:19

daogege 发表于 2017-6-21 20:36
请问，这里面的交换相关渐进式的矫正是必须的吗？没有实验值的话，可以用计算值吗

是。
可以用高精度计算值，比如CCSD(T)/cc-pVQZ

才雪 · 发表于 Post on 2022-1-14 20:39:25

请问这里的HF能换成PBE么

lcdamoy · 发表于 Post on 2022-1-15 21:34:29

才雪发表于 2022-1-14 20:39
请问这里的HF能换成PBE么

不可以

桂圆柿子 · 发表于 Post on 2023-5-26 11:28:47

sobereva 发表于 2016-12-6 00:40
这么大体系想做DFT-SAPT基本没戏，算不动
建议你用GAMESS-US做LMOEDA能量分解，这种规模还算能算得动
Gau ...

咨询一下老师，50个原子左右的体系做SAPT0可以算动吗

sobereva · 发表于 Post on 2023-5-27 07:25:21

桂圆柿子发表于 2023-5-26 11:28
咨询一下老师，50个原子左右的体系做SAPT0可以算动吗

可以

桂圆柿子 · 发表于 Post on 2023-6-14 14:18:14

sobereva 发表于 2023-5-27 07:25
可以

好的，谢谢老师，我这里实践是可以的。

秋心 · 发表于 Post on 2024-7-30 21:59:01

老师您好，想问一下对于一个中性分子一个带一个负电荷的离子（硫酸氢根）构成的二聚体计算SAPT时，电荷应该怎么写？我看手册上说WF,[NELEC=nelec],[SYM[METRY]=irrep],[spin= spin],[CHARGE=charge]，这一部分可以设置分子的电荷以及自选多重度，我自己尝试写了一下，拿不准，麻烦老师能帮我看一下吗？附件中分子B就是负离子。

		自动登录 Automatic login	找回密码 Forget password
密码 Password			注册 Register

[Molpro] molpro中SAPT能量分解求助