如何计算电子的Wannier轨道中心？

gingaorb

各位小哥哥，各位大佬们，请教一个问题呀~
是这样的，最近我在重复一篇关于水溶液中溶解有两个电子的体系的文章，文章的题目是“

Molecular Dynamics Characterization of Dielectron Hydration in Liquid Water with Unique Double Proton Transfers

”，发在JCTC上，之前在站长大大的CP2K培训讲义中也提到了这篇文章，详细教了我们如何用CP2K设置参数去模拟水中溶解电子的体系。
现在我希望能像文章图1B中那样，去计算水合电子的径向分布函数（RDF），文章中提到“通过过剩电子（S3中的图）的Wannier轨道中心（WOC）与溶剂H2O的O或H原子之间的径向分布函数（RDF）”。意思是说通过计算水合电子的Wannier轨道中心的位置，作为径向分布函数计算的原点位置。请问这个Wannier轨道中心有没有什么比较好的办法去计算和确定呀？比如用CP2K联合站长大大开发的Multiwfn软件能不能做到呢？


谢谢大家~

以下是我所参考的文章：

sobereva

Multiwfn做轨道定域化的时候直接就可以要求把定域化轨道的质心位置给出，下文说了。基于CP2K的molden文件可以对周期性体系做
Multiwfn的轨道定域化功能的使用以及与NBO、AdNDP分析的对比
http://sobereva.com/380（http://bbs.keinsci.com/thread-6053-1-1.html）

highlight

https://manual.cp2k.org/trunk/CP ... ANNIER_CENTERS.html

gingaorb

highlight 发表于 2023-10-19 08:13
https://manual.cp2k.org/trunk/CP2K_INPUT/FORCE_EVAL/DFT/LOCALIZE/PRINT/WANNIER_CENTERS.html

你好呀你好呀！我按照您的提示设置了WANNIER_CENTERS，跑出来体系全部的WANNIER_CENTERS，以下是文件的头几行内容：
   273
1_0
X     -8.60466311     -8.50704784      9.75358398      1.13542678
X    -10.03166259     -8.66810850      7.81756337      1.46142051
X     -6.98420241    -11.35492078    -10.78278863      0.86987786
X     -9.11487975     11.10467697     -9.18618679      1.12202325
X    -10.80799432     11.99922239     -6.76559405      1.06902915
X     -9.86097161     -8.45320591    -11.17988242      0.89247883
X     -9.05482060     11.34349378     12.45059165      0.90257299
X    -12.37699676    -10.58238531     -9.27741512      0.92725393
我有三个问题想请教以下，
首先是，我的体系在一个边长25埃的三个方向周期格子中，所有坐标都在0到25的范围内，为什么文件中的坐标会有负值？是因为体系的周期性吗？
另外，打印出这么多轨道中心，它的排列顺序有什么规则吗？我如何找到我要的水合电子的轨道中心呢？
最后，文件中的第四列内容又是什么含义呢？
非常感谢！期待您的回复！

highlight

1. 是
2. 应该就是 Wannier 轨道顺序吧。我的体系里水合电子在最高占据轨道，使用https://manual.cp2k.org/trunk/CP ... L.DFT.LOCALIZE.LIST 就可以设置只计算那些轨道了。
3. 不知道。我开了 IONS_CENTERS，X 也只显示 3 坐标，也方便直接用 vmd 计算 X-O 的 RDF。
PS：其实使用 WANNIER_CENTERS 定义过剩电子中心，是考虑周期性导致的，簇的话直接用 http://sobereva.com/658 中的 r_c 就行。
不知道 Multiwfn 有没有这方面考量，有的话直接用也行。

gingaorb

highlight 发表于 2024-1-4 16:36
1. 是
2. 应该就是 Wannier 轨道顺序吧。我的体系里水合电子在最高占据轨道，使用https://manual.cp2k.org ...

真的非常感谢您的解答！关于您的解答我说一下我的理解，如果理解有误还请您指正！
关于您的“2.”的回答，按照我的理解，请问输入文件中的这段内容是否应该设置成以下内容：
    &LOCALIZE
      METHOD CRAZY
      EPS_LOCALIZATION 1.0E-8
      LIST 1 2
      &PRINT
        ! Output the centers of all Wannier functions as xyz-file
        &WANNIER_CENTERS
           IONS+CENTERS .TRUE.
        &END
      &END
    &END      

其中，METHOD CRAZY这个关键词我是参考链接：https://www.cp2k.org/exercises:2014_ethz_mmm:wannier?s[]=wannier&s[]=centers设置的，我的理解是，既然要算WANNIER_CENTERS，那么首先应该要使用某种方法计算WANNIER函数才可以吧？所以我设置了METHOD CRAZY以及EPS_LOCALIZATION 1.0E-8。

另外，按照您的提示，意思是不是加上LIST 1 2就只打印最高占据轨道了呢？我不是很懂其中LIST 1 2的含义，希望您能指导一下。

然后根据您“3.”的指导，我查阅官网手册添加了IONS+CENTERS .TRUE.关键词，应该是您提到的开了“IONS_CENTERS”吧？这样一来，计算得到的输出将包括离子（原子核）位置和 Wannier 中心。想请问一下，最后您提到的用 vmd 计算 X-O 的 RDF，具体应该如何操作呢？

最后关于您“PS”提到的内容，我也同样做了尝试，根据站长大大的链接http://sobereva.com/658，我确实计算得到了r_c，因为我的体系中有两个自旋配对的水合电子，所以积分的时候需要取绝对值。但是接下来我计算回旋半径得到的结果却很奇怪，Multiwfn计算得到的是11 a.u.多的值，如果按照Bohr和埃之间的换算关系1Bohr约0.529埃，和文献中的二电子回旋半径约2.68埃，相差甚大，所以我想通过设置WANNIER_CENTERS计算的结果进行比较，排查问题出在哪里。

再次向您表达我诚挚的感谢，期待您的回复！

highlight

“METHOD CRAZY以及EPS_LOCALIZATION 1.0E-8” 可以
“LIST 1 2” 不对，你计算的是两个最低占据轨道，不是 HOMO。
LIST 写多少你得根据体系的电子数算。比如 64 个水有 64*8 =512 个电子，中性体系就有 256 个占据轨道。
加入一个过剩电子并占据 HOMO 就是 257；两个的话，单重态是 257，三重是 257 258，以此类推。
偷懒的话可以看输出轨道 cube 时的文件名。
是 IONS+CENTERS。vmd 计算 RDF 的教程论坛里以前就有不少。
“两个自旋配对的水合电子，所以积分的时候需要取绝对值”不太能理解，是 BS 态还是 CS ? 谁的绝对值？自旋密度？我用的轨道空间分布，并没用自旋密度。
至于 Multiwfn 计算的结果是不是有问题，需要具体分析。
一方面你是不是真得到了定域化的水合电子结构，这是要优先怀疑的；
另一方面就像我之前说的，不清楚Multiwfn有没有对周期性体系做特殊处理，如果没有，当 r_c 处于某周期性边界附近是，可能会有问题。