KTTK：用于给Gaussian组装片段组合波函数的Python程序

Accelerator

通过片段组合的方法构造复合物的初始波函数是一种比较可控的生成初猜的方法，但Gaussian的guess=fragment功能比较局限，不能进行波函数稳定性检查，也不能读取已经收敛好的片段信息。如果针对某个片段直接收敛出来的波函数并不稳定，最终组装起来的波函数就需要再进行稳定性检查才能用，这往往是十分痛苦的。KTTK是为了解决这一问题而开发的Python程序（名称来源于kataduke片付け）。它分为两个部分：kttk_mixer.py，读取各个片段的fchk文件并组装整体的轨道，然后以fchk形式输出——这当然要求整体的原子坐标和顺序要和各片段完全相同。为了方便地生成满足顺序要求的片段gjf文件，可以使用kttk_separator.py。以下用FeCp2 Cl为例（仅因为它比较简单，不代表有任何化学意义）来展示它的用法。
假如有这样一个dd.gjf，记录整体的结合结构。选取这个例子还有一个原因，是wB97xD直接得到的FeCp2+的波函数是不稳定的，导致没办法用guess=fragment来直接生成正确的片段组合波函数：

1. %chk=dd_mixed.chk
   
2. %mem=24GB
   
3. %nprocshared=24
   
4. # uwb97xd nosymm guess=read g09default scf=xqc
   
5. 
   
6. TC
   
7. 
   
8. 0 2
   
9. H                 -4.49778100    0.71661800    1.28221500
   
10. C                 -5.31465600    0.03561600    1.04781000
    
11. C                 -5.54736900   -0.59233100   -0.21012500
    
12. C                 -6.30355300   -0.42048600    1.95688900
    
13. C                 -6.68257300   -1.44779300   -0.06983100
    
14. H                 -4.95473500   -0.45583000   -1.11330700
    
15. Fe                -5.16530800   -2.06921600    1.17881100
    
16. C                 -7.14298200   -1.33867000    1.27525600
    
17. H                 -6.36571000   -0.16448200    3.01381400
    
18. H                 -7.11211600   -2.08084700   -0.84501900
    
19. H                 -7.97510600   -1.88770500    1.71480600
    
20. C                 -3.39299100   -2.94799800    0.61681100
    
21. C                 -4.45139900   -3.89480800    0.47859800
    
22. C                 -3.31231500   -2.58327100    1.99432000
    
23. H                 -2.76848700   -2.56194300   -0.18729600
    
24. C                 -5.00888300   -4.11958000    1.76821100
    
25. H                 -4.79012500   -4.34732700   -0.45175100
    
26. C                 -4.30458100   -3.31729600    2.69971200
    
27. H                 -2.62319900   -1.85936300    2.42673500
    
28. H                 -5.86499300   -4.75446900    1.99190800
    
29. Cl                -2.04406600    0.01490100   -1.58085900
    
30. H                 -4.52805800   -3.22843700    3.76207200

复制代码

可见其中原子顺序不连续，1到20和22号原子是FeCp2的，21号原子是Cl。接下来执行：

1. $ python3 kttk_separator.py
   
2. Input the gjf file for your whole molecule:
   
3. dd.gjf
   
4. Input your atom indexes for each fragment. E.g. 1-5,10,12-14
   
5. Empty line to end the input, and the unassigned atoms will be automatically assigned into a new fragment.
   
6. 1-20,22
   
7. 
   
8. Input the charge and multiplicity for each fragment, respectively.              Empty line to end the input. E.g. 0 1
   
9. 1 2
   
10. -1 1
    
11. 
    
12. Input the kwds. Do not write guess=xxx.
    
13. wb97xd sto-3g g09default
    
14. Now the gjf files for each fragment is outputted. Modify them to adapt to your task, run,  and use kttk_mixer.py to build the combined fchk.

复制代码

就会生成3个文件（片段的数量可以任意多）：dd_1.gjf, dd_2.gjf和dd_mixed.gjf。其中dd_1和dd_2分别记录了FeCp2+和Cl-片段，而dd_mixed.gjf则是按照片段中的原子顺序重新排列得到的整体gjf文件。为了生成稳定的片段波函数，可以任意修改这些文件，特别是如果想要研究一系列包含FeCp2+的体系，只需要找到一个收敛的稳定波函数，以后每次对于片段1都读取它就行了。
最终得到dd_1.fchk和dd_2.fchk，然后运行：

1. $ python3 kttk_mixer.py
   
2. Input the fragment fchk files, and an empty line to end the input.
   
3. dd_1.fchk
   
4. dd_2.fchk
   
5. 
   
6. Input the flag (1 or -1) for each fragment, respectively. -1 to flip the spin. E.g. 1 -1 1
   
7. 1 1
   
8. Where to save the output fchk file? dd_mixed.fchk

复制代码

从而得到dd_mixed.fchk。unfchk之后，让dd_mixed.gjf读取它，就可以收敛到正确的波函数。在这个例子中，通过调整两个片段的电荷/自旋多重度分别为1 2/-1 1和0 1/0 2，可以收敛到两个态，第一个态的自旋密度完全位于Fe上，第二个态中Cl也分担了一定自旋密度。

为了防止初学者迷路，特此强调：当设法控制收敛到什么态时，应该很清楚自己在做什么。笔者编写KTTK时，是为了用单参考态方法近似地研究一些电子转移反应，需要考虑两个态的交叉问题，并且这两个态的波函数都是稳定的。如果你对自己为何要进行这些控制尚不清楚，应该先想清楚再行动。

程序地址：https://github.com/RimoAccelerator/KTTK

hebrewsnabla

正巧我最近也在写片段组合的程序

冰释之川

hebrewsnabla 发表于 2022-12-5 19:17
正巧我最近也在写片段组合的程序

还不共享一波*.ipynb？