| 查看: 2637 | 回复: 3 | |||
| 【奖励】 本帖被评价2次，作者beefly增加金币 1.6 个 | |||
| 当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖 | |||
[资源]
用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 ] |
回复此楼» 猜你喜欢
304求调剂
已经有5人回复
-
1
已经有21人回复
-
物理化学论文润色/翻译怎么收费?
已经有55人回复
-
0703 调剂
已经有5人回复
-
上海华东理工大学机械与动力工程学院先进材料与制造研究组招聘公告
已经有0人回复
-
304求调剂
已经有6人回复
-
西南某211招27届博士
已经有30人回复
-
大连工业杰青/长江团队-生物质材料-储能电池方向招收2026级博士生
已经有0人回复
-
304求调剂
已经有1人回复
-
中药学306，求调剂
已经有0人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 在Gaussian09计算中考虑标量相对论和旋轨耦合
已经有9人回复
★★★★★ 五星级,优秀推荐
|
请问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
4楼2014-01-30 12:26:47
简单回复
chrinide2楼
2014-01-27 19:28
回复
五星好评 
