24小时热门版块排行榜

返回列表

【奖励】本帖被评价2次，作者beefly增加金币 1.6 个

beefly

专家顾问 (职业作家)

[资源] 用MOLPRO构造原子自然轨道（ANO）基函数

在做高精度计算的时候，经常遇到现有基函数不能满足需要的情况，例如，ECP基组太小，收缩因子未考虑半芯轨道的关联，缺少Rydberg基函数，未考虑相对论效应，等等。这个时候与其怨天尤人，还不如亲自动手，构造符合需要的基函数。

除非有特殊要求，一般量化计算所用的基函数都是收缩的高斯型函数。基函数的收缩方式有分块收缩和一般收缩两类，前者以Pople基组（如STO-3G、6-31G）和Ahlrichs基组（如早期的nZVP和后来的def2-nZVP）为代表，后者以Almolf、Roos等人的ANO基组（ANO-L，ANO-S，ANO-RCC）为代表。至于由Dunning最早开发，后来由Peterson进一步扩充的cc-pVnZ型基组（包括aug-cc-pVnZ，cc-pwCVnZ，cc-pVnZ-PP，cc-pVnZ-DK2等），在收缩方式上介于两者之间，但由于收缩因子的计算方法以ANO为基础，因此通常都作为一般收缩基函数来看待。

本文只涉及一般收缩基函数。关于两种收缩方式的区别，读者可以去查文献，或者翻过去的帖子，这里不再赘述。需要注意的是，一般收缩基函数最好用专门编写的双电子积分代码。由于大多数量子化学程序的积分代码都是为分块收缩基函数编写的，在处理一般收缩基函数的时候往往效率不高，通常要慢一两个数量级，甚至根本算不动。适合一般收缩基函数的程序主要有SEWARD（被MOLPRO和MOLCAS采用）和ABACUS（Dalton程序专用。旧版本也包含在ACES2，CFour，Columbus中）。此外，新版本的Firefly（PC-Gamess）和Gaussian09在处理一般收缩基函数时也较旧版有很大改善。

构造收缩基函数，首先要有未收缩的高斯原函数。原函数有三个来源：文献发表的非收缩函数（如Gaussian程序自带的UGBS），根据公式自动产生原函数（如even-tempered函数），以及对现成的收缩基函数进行修改。建议用第三种方法，因为在相同的精度下，这些经过优化的高斯原函数比其它方法用到的函数少，因此具有较高的计算效率。

有了高斯原函数之后，下一步就是进行收缩，通过从头计算（通常是MRCISD，也可以用CASPT2）获得ANO收缩因子。MOLCAS程序有一个专门的工具GENANO产生ANO收缩因子。而在MOLPRO程序中，利用矩阵操作命令matrop，也可以实现相同的目的。

下面用一个例子来说明如何造ANO基组。文献中，氧原子的ECP基组最高级别是cc-pVTZ-PP。这里我们做一个cc-pVQZ-PP基组。ECP用CRENBL，高斯原函数来自全电子的cc-pVQZ。ECP和基函数可以从这里获得：https://bse.pnl.gov/bse/portal

我们希望这套收缩函数能够描述O原子的基态3Pg（D2h：b1g + b2g + b3g）和成键后带有部分负电荷的O-基态2Pu（D2h：b1u + b2u + b3u），因此将对这些态计算平均的电子密度和ANO。为了描述激发态，还可以再加上O原子的第一激发态1Dg（D2h：2*ag + b1g + b2g + b3g）进行平均，不过作为演示，这里不考虑。

（注：如果原子p、f轨道具有奇数占据，则电子态具有u对称性；否则都是g对称性。O和O-基态的电子组态分别为2s2 2p4和2s2 2p5，因此分别具有g和u对称性。）

为了保持电子态和轨道的球对称性，HF计算闭壳层的O--，接下来做态平均CASSCF计算。输入如下。由于用ECP替代2个1s芯电子，O，O-，O--的电子数分别是6，7，8。在O的全电子cc-pVQZ基函数中，只有s，p函数是收缩的，并且d，f，g轨道都是非占据的，因此这个输入里只有s，p函数，非收缩的极化函数在做完收缩之后补上。

***,O
memory,100,m
gprint,orbital=5

basis={
! O: ECP from
! L.F. Pacios and P.A. Christiansen, J. Chem. Phys. 82, 2664 (1985).
ecp,O,2,1,0;
3; 2, 10.028600, -0.798420;    2, 34.198000, -5.766847;       ! p-term
   1,  100.003900, -1.486456;
4; 2, 2.247900, 11.216304;    2, 2.404900,  -16.344477;       ! s-term
   1, 4.374000, 1.042944;    0, 2.189200, 2.193891;

s, O , 61420.0000000, 9199.0000000, 2091.0000000, 590.9000000, 192.3000000, 69.3200000, 26.9700000, 11.1000000, 4.6820000, 1.4280000, 0.5547000, 0.2067000
p, O , 63.4200000, 14.6600000, 4.4590000, 1.5310000, 0.5302000, 0.1750000

}

geometry={O}

{hf;occ,1,1,1,0,1,0,0,0;wf,8,1,0;}

! CASSCF: 2s2p are active
{multi;NATORB,ci;start,2100.2;
frozen,0;closed,0;occ,1,1,1,0,1,0,0,0;
! O:3Pg
wf,6,4,2;state,1;
wf,6,6,2;state,1;
wf,6,7,2;state,1;
! O-: 2Pu
wf,7,2,1;state,1;
wf,7,3,1;state,1;
wf,7,5,1;state,1;
}
---

计算结束后检查轨道因子，发现由于使用了ECP，导致一些s函数的占据轨道因子非常小。这些函数必须删除。一般情况下，在占据轨道中因子小于0.0001的（甚至可以扩大到0.001）都可以删除。这里删除前4个s函数。CASSCF之后，对6个根分别执行MRCI计算，并计算平均电子密度及其ANO因子。输入文件如下：

***,O
memory,100,m
gprint,orbital=5

basis={
! O: ECP from
! L.F. Pacios and P.A. Christiansen, J. Chem. Phys. 82, 2664 (1985).
ecp,O,2,1,0;
3; 2, 10.028600, -0.798420;    2, 34.198000, -5.766847;       ! p-term
   1,  100.003900, -1.486456;
4; 2, 2.247900, 11.216304;    2, 2.404900,  -16.344477;       ! s-term
   1, 4.374000, 1.042944;    0, 2.189200, 2.193891;

s, O , 192.3000000, 69.3200000, 26.9700000, 11.1000000, 4.6820000, 1.4280000, 0.5547000, 0.2067000
p, O , 63.4200000, 14.6600000, 4.4590000, 1.5310000, 0.5302000, 0.1750000

}

geometry={O}

{hf;occ,1,1,1,0,1,0,0,0;wf,8,1,0;}

! CASSCF: 2s2p are active
{multi;NATORB,ci;start,2100.2;
frozen,0;closed,0;occ,1,1,1,0,1,0,0,0;
! O:3Pg
wf,6,4,2;state,1;
wf,6,6,2;state,1;
wf,6,7,2;state,1;
! O-: 2Pu
wf,7,2,1;state,1;
wf,7,3,1;state,1;
wf,7,5,1;state,1;
}

{ci;wf,6,4,2;natorb,2351.2,state=1.4;dm,2351.2;}
{ci;wf,6,6,2;natorb,2352.2,state=1.6;dm,2352.2;}
{ci;wf,6,7,2;natorb,2353.2,state=1.7;dm,2353.2;}
{ci;wf,7,2,1;natorb,2354.2,state=1.2;dm,2354.2;}
{ci;wf,7,3,1;natorb,2355.2,state=1.3;dm,2355.2;}
{ci;wf,7,5,1;natorb,2356.2,state=1.5;dm,2356.2;}

cf=1/6;

{matrop
load,den01,den,2351.2;
load,den02,den,2352.2;
load,den03,den,2353.2;
load,den04,den,2354.2;
load,den05,den,2355.2;
load,den06,den,2356.2;
add, densum, cf,den01, cf,den02, cf,den03, cf,den04, cf,den05, cf,den06;
natorb,norb1,densum,1.d-10;
}
---

MOLPRO最后打印的ANO占据数和因子如下：

Orb    Occ          Coefficients

                     1 1s    1 1s    1 1s    1 1s    1 1s    1 1s    1 1s    1 1s

1.1 1.9862178 -0.000502 -0.001844 -0.005858 -0.045350 -0.130504  0.302721  0.524298  0.334656

2.1 0.0119301 0.000435  0.004099  0.001586  0.090359  0.147742 -0.969071 -0.457808  1.284988

（略）

                     1 2px    1 2px    1 2px    1 2px    1 2px    1 2px

1.2 1.4804766 0.005798  0.040341  0.153455  0.343784  0.392650  0.341316

2.2 0.0188067 -0.007237 -0.046058 -0.213095 -0.468113 -0.184618  0.963171

（略）

                     1 2py    1 2py    1 2py    1 2py    1 2py    1 2py

1.3 1.4804766 0.005798  0.040341  0.153455  0.343784  0.392650  0.341316

2.3 0.0188067 -0.007237 -0.046058 -0.213094 -0.468114 -0.184618  0.963171

（略）

                     1 2pz    1 2pz    1 2pz    1 2pz    1 2pz    1 2pz

1.5 1.4804766 0.005798  0.040341  0.153455  0.343784  0.392650  0.341316

2.5 0.0188067 -0.007237 -0.046058 -0.213094 -0.468114 -0.184619  0.963172

（略）

可以看到，p轨道的简并性保持得很好。已知cc-pVQZ的收缩形式为(12s,6p,3d,2f,1g)/[5s,4p,3d,2f,1g]，并且只有前2个s和第一个p是做了ANO收缩的。那么去掉1s轨道及其对应的4个s函数后，cc-pVQZ-PP的收缩形式应当是(8s,6p,...)/[4s,4p,...]，并且只有第一个s和p是收缩的，收缩因子来自第一个s型ANO和第一个p型ANO。非收缩的函数仍然保持。基组部分经过整理，就得到cc-pVQZ-PP。如果再加上来自aug-cc-pVQZ的弥散函数，就得到aug-cc-pVQZ-PP。如下：

basis={
! O: ECP from
! L.F. Pacios and P.A. Christiansen, J. Chem. Phys. 82, 2664 (1985).
ecp,O,2,1,0;
3; 2, 10.028600, -0.798420;    2, 34.198000, -5.766847;       ! p-term
   1,  100.003900, -1.486456;
4; 2, 2.247900, 11.216304;    2, 2.404900,  -16.344477;       ! s-term
   1, 4.374000, 1.042944;    0, 2.189200, 2.193891;

! (8s,6p,3d,2f,1g) -> [4s,4p,3d,2f,1g] + [1s,1p,1d,1f,1g]
! from the exponents of aug-cc-pVQZ by deleting first 4 s-functions, and recontracted by MR-CISD.
s, O , 192.3000000, 69.3200000, 26.9700000, 11.1000000, 4.6820000, 1.4280000, 0.5547000, 0.2067000;
c, 1.8, -0.000502,-0.001844,-0.005858,-0.045350,-0.130504, 0.302721, 0.524298, 0.334656;
c, 6.6, 1.0;
c, 7.7, 1.0;
c, 8.8, 1.0;
p, O , 63.4200000, 14.6600000, 4.4590000, 1.5310000, 0.5302000, 0.1750000;
c, 1.6, 0.005798, 0.040341, 0.153455, 0.343784, 0.392650, 0.341316;
c, 4.4, 1.0;
c, 5.5, 1.0;
c, 6.6, 1.0;
d, O , 3.7750000, 1.3000000, 0.4440000;
f, O , 2.6660000, 0.8590000;
g, O , 1.8460000;
! Diffuse
s, O , 0.0695900;
p, O , 0.0534800;
d, O , 0.1540000;
f, O , 0.3240000;
g, O , 0.7140000;
}

[ Last edited by beefly on 2014-1-27 at 20:59 ]

回复此楼

» 猜你喜欢

» 本主题相关价值贴推荐，对您同样有帮助:

在Gaussian09计算中考虑标量相对论和旋轨耦合已经有9人回复

1楼 2014-01-27 15:29:26

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

luqing6879

木虫 (著名写手)

应助: 8 (幼儿园)
金币: 3563.5
帖子: 1907
在线: 792.8小时
虫号: 614833

★★★★★ 五星级,优秀推荐

请问ECP这部分的输入中，这些整数是什么意思？
3; 2, 10.028600, -0.798420;    2, 34.198000, -5.766847;       ! p-term
   1,  100.003900, -1.486456;
4; 2, 2.247900, 11.216304;    2, 2.404900,  -16.344477;       ! s-term
   1, 4.374000, 1.042944;    0, 2.189200, 2.193891;

另外请问(8s,6p,...)/[4s,4p,...]，为什么是收缩成4s,而不是5s？不好意思没看懂。

还有最后的recontraction。怎么得到的？

[ Last edited by luqing6879 on 2014-1-28 at 08:29 ]

赞一下

回复此楼

3楼2014-01-28 04:03:28

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

beefly

专家顾问 (职业作家)

引用回帖:

3楼: Originally posted by luqing6879 at 2014-01-28 04:03:28
请问ECP这部分的输入中，这些整数是什么意思？
3; 2, 10.028600, -0.798420; 2, 34.198000, -5.766847; ! p-term
1, 100.003900, -1.486456;
4; 2, 2.247900, 11.216304; ...

第一部分是L=P的赝势函数，作为局域项。之后是0,1,...,L-1的函数在局域项上投影。在P赝势数据区中，第一个3表示有3项P赝势函数的线性组合。每一个赝势函数由三个数定义：第一个整数定义r的幂指数（实际用到的幂指数是r^(n-2)）,然后是e的幂指数因子gamma，和线性组合因子。

因为全电子基组收缩成5个s，现在1s轨道不存在了，故删除一个收缩的s，变成4个s。

从自然轨道中找到第一个s和第一个p的因子，替换原来的收缩因子。非收缩函数不变。由于s型原函数数量少了，因此需要修改编号。

赞一下(2人)

回复此楼

4楼2014-01-30 12:26:47

已阅回复此楼关注TA 给TA发消息送TA红花 TA的回帖

简单回复

chrinide2楼

2014-01-27 19:28 回复

五星好评

相关版块跳转我要订阅楼主 beefly 的主题更新

返回列表

☆ 无星级 ★ 一星级 ★★★ 三星级 ★★★★★ 五星级

普通表情龙兔虎猫高级回复 (可上传附件)