CP2K 计算自旋极化单重态%BS输入以及计算结果求教

chands

各位老师，
      小弟初学CP2K，试计算了臭氧（O3）基态的优化。臭氧分子两端的氧原子各有一个未成对电子，基态下这两个电子自旋反向，所以总自旋为0，单重态。
      用了Sob老师的Multiwfn生成CP2K输入文件，参考了CP2K官方手册（顺便吐槽一下，语焉不详，实在难看）还有帖子CP2K如何计算自旋极化单重态？ - 第一性原理 (First Principle) - 计算化学公社 (keinsci.com)，根据自己的理解在输入文件中添加%BS模块。      %BS模块指定Alpha自旋和Beta自旋的电子数，有三个参数，N是主量子数，L是角量子数，NEL官方解释是: "Orbital ccupation change per angular momentum quantum number. In unrestricted calculations applied to spin alpha/beta. This keyword cannot be repeated and it expects a list of integers. Default value: -1”
       ？？？？看了一些例子，NEL似乎是这样的：如果某原子主量子数为N角量子数为L的亚层上有X个电子，那么要指定Alpha（或Beta）自旋的电子数=（X+ NEL)/2，奇怪的是不管角量子数L为何值，一概是（X+NEL）/2，我糊了，”Orbital ccupation change per angular momentum quantum number“是什么意思，大牛指教？
      然后我按照自己的理解，在%BS模块中给三个O原子指定Alpha和Beta自旋的电子数（中间O为O1，两端为O2，O3）
      O1：（2p亚层，N=2, L=1) 三个alpha，一个beta
      O2：（2p亚层，N=2, L=1) 两个alpha，两个beta
      O3：（2p亚层，N=2, L=1) 一个alpha，三个beta
由于2p亚层有四个电子，输入文件就成这个模样（完整输入文件见附件）：
    &KIND O1   
      ELEMENT O
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-BLYP
      &BS
        &ALPHA
           N  2
           L   1
           NEL +2
        &END ALPHA
        &BETA
           N 2
           L  1
           NEL  -2
        &END BETA
      &END BS
    &END KIND
    &KIND O2   
      ELEMENT O
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-BLYP
      &BS
        &ALPHA
           N  2
           L   1
           NEL 0
        &END ALPHA
        &BETA
           N 2
           L  1
           NEL  0
        &END BETA
      &END BS
    &END KIND
    &KIND O3   
      ELEMENT O
      BASIS_SET DZVP-MOLOPT-SR-GTH
      POTENTIAL GTH-BLYP
      &BS
        &ALPHA
           N  2
           L   1
           NEL -2
        &END ALPHA
        &BETA
           N 2
           L  1
           NEL  +2
        &END BETA
      &END BS
    &END KIND

先在我的小破机上跑一跑，out文件中有一列是这样的
  Ideal and single determinant S**2 :                    0.000000       0.554552


至此，有几个问题要请教
（1）我的理解要指定Alpha（或Beta）自旋的电子数=（X+ NEL)/2对不对？
   (2)  官方首次NEL是”Orbital ccupation change per angular momentum quantum number“，是什么意思？
   (3)  网上有些例子连O2-离子都要指定NEL，因为O2-是2p^6，三个alpha, 三个beta, 与O原子2p^4不同，这样看来，是不是所有离子都要指定NEL？
   (4)  O原子基态是2px^2, 2py^1, 2pz^1，是不是也要指定
      &BS
        &ALPHA
           N  2
           L   1
           NEL +2
        &END ALPHA
        &BETA
           N 2
           L  1
           NEL  -2
        &END BETA            &END BS
或者
      &BS
        &ALPHA
           N  2
           L   1
           NEL -2
        &END ALPHA
        &BETA
           N 2
           L  1
           NEL +2
        &END BETA
      &END BS

    （5）输出文件Ideal and single determinant S**2 中的single determinant是什么意思，single determinant不为0算不算自旋污染，CP2K有spin annihilation的方法吗？

sobereva

不要自己在标题上写[CP2K]标签，扰乱论坛秩序
明明第一性原理版有CP2K分类，CP2K明显属于第一性原理程序，相关问题别发到量子化学版。给你移动了并修改了标题，以后注意

sobereva

设置初猜没必要那么讲究。只需要设置MAGNETIZATION让O3两头的两个氧带的单电子自旋相反就完了，电子结构会在SCF过程中自发调节。

另外，算这种体系严重不适合CP2K，效率远低于量化程序，而且还没有波函数稳定性测试的功能。

搞懂不同自旋多重度的S^2算符期望值的理想值是多少（单重态为0），这属于量化基本常识知识。这都不了解的话先学点量化基础知识。
什么叫自旋污染，公社论坛首页全文搜索框搜，有无数回复和讨论。spin annihilation跟你没关系

此文好好看了
谈谈片段组合波函数与自旋极化单重态
http://sobereva.com/82
CASSCF计算双自由基以及双自由基特征的计算
http://sobereva.com/264（http://bbs.keinsci.com/thread-322-1-1.html）

Daniel_Arndt

可以参考 https://www.cp2k.org/_media/even ... al:ling_hybrids.pdf 和 https://github.com/cp2k/cp2k/tree/master/tests/QS/regtest-bs 。
第三个问题，你可能是看到了第一个链接中的三氧化二铁的资料。其实很简单，默认情况下初猜是不带电荷的氧原子，为了方便设置对称破缺，把不带电荷的氧原子改成带两个单位的负电荷的离子（对应于第34页的内容）。相应地，不带电荷的铁原子就得变成带三个单位的正电荷的离子。有一层铁原子，是3d轨道上的5个电子都是alpha（对应于第35页的内容）；紧邻的一层铁原子，就是3d轨道上的5个电子都是beta（对应于第36页的内容）。

chands

Daniel_Arndt 发表于 2022-1-9 05:50
可以参考 https://www.cp2k.org/_media/even ... al:ling_hybrids.pdf 和 https://github.com/cp2k/cp2k/tr ...

谢谢老师这么早起来指教。我还有两个问题
（1）如果按照ling_hybrids.pdf上的例子，设置O2-是为了更好设置3d轨道半空的Fe3+的对称破缺，那么如果要计算MgO，Mg2+和O2-都是轨道全满，没有对称破缺，是不是根本就不用初猜指定自旋？
（2）您的第二个链接是计算[Cu2Cl6]2-的，其中cu2cl6_m1_std.inp文件中指定自旋初猜的部分是
    &KIND  Cu1
      ELEMENT Cu
      BASIS_SET  DZVP-MOLOPT-SR-GTH
      POTENTIAL  GTH-BLYP-q11
      &BS
        &ALPHA
          NEL -1   0
          L    0   2
          N    4   3
        &END
        &BETA
          NEL -1  -2
          L    0   2
          N    4   3
        &END
      &END
    &END
    &KIND  Cu2
      ELEMENT Cu
      BASIS_SET  DZVP-MOLOPT-SR-GTH
      POTENTIAL  GTH-BLYP-q11
      &BS
        &ALPHA
          NEL -1  -2
          L    0   2
          N    4   3
        &END
        &BETA
          NEL -1   0
          L    0   2
          N    4   3
        &END
      &END
    &END
    &KIND  Cl
      BASIS_SET  DZVP-MOLOPT-GTH
      POTENTIAL  GTH-BLYP-q7
      &BS
        &ALPHA
          NEL 2
          L   1
          N   3
        &END
        &BETA
          NEL 2
          L   1
          N   3
        &END
      &END
    &END
    &KIND  Cl1
      ELEMENT Cl
      BASIS_SET  DZVP-MOLOPT-GTH
      POTENTIAL  GTH-BLYP-q7
      &BS
        &ALPHA
          NEL 2
          L   1
          N   3
        &END
        &BETA
          NEL 0
          L   1
          N   3
        &END
      &END
    &END
    &KIND  Cl2
      ELEMENT Cl
      BASIS_SET  DZVP-MOLOPT-GTH
      POTENTIAL  GTH-BLYP-q7
      &BS
        &ALPHA
          NEL 0
          L   1
          N   3
        &END
        &BETA
          NEL 2
          L   1
          N   3
        &END
      &END
    &END

Cu2+离子的设定我明白，但是Cl-离子我有点糊，Cl原子是3p^5, 如果设定Cl-, 3p^6, 三个alpha三个beta,是不是应该为
  &KIND  Cl1
      ELEMENT Cl
      BASIS_SET  DZVP-MOLOPT-GTH
      POTENTIAL  GTH-BLYP-q7
      &BS
        &ALPHA
          NEL 1
          L   1
          N   3
        &END
        &BETA
          NEL 1
          L   1
          N   3
        &END
      &END
    &END

chands

sobereva 发表于 2022-1-9 05:48
设置初猜没必要那么讲究。只需要设置MAGNETIZATION让O3两头的两个氧带的单电子自旋相反就完了，电子结构会 ...

谢谢社长指教。我用Gaussian优化O3几秒钟就好了，然后将优化好的结构用CP2K算，程序蹦啊跳啊好久才算好。我是想跑一个水溶液中的臭氧反应的MD，试过ORCA，发现用AIMD来算ORCA根本跑不动，于是想改为QM/MM，但是ORCA用MM需要水的力场文件，这个文件不好弄，于是想用CP2K来计算QM/MM-MD，不知这种情况社长有何建议。

关于S^2，谢谢社长解惑，以前真的不知道这种情况用UDFT或UHF算出来S^2会有很大偏差，我看S^2怎么跟0差这莫大都吓坏了。

Daniel_Arndt

chands 发表于 2022-1-9 11:23
谢谢老师这么早起来指教。我还有两个问题
（1）如果按照ling_hybrids.pdf上的例子，设置O2-是为了更好设 ...

如果是MgO的话，你用默认的初猜，也就是不带电荷的镁原子和不带电荷的氧原子，应该是可以收敛的。你要是愿意折腾的话，可以手动用BS设置成带两个单位正电荷的镁离子和带两个单位负电荷的氧离子，对结果的正确性应该不会有影响；第二个链接其实就是CP2K自带的tests/QS/regtest-bs。这应该是一个分子体系（而不是晶体），但我自己算过的对称破缺体系中，没有这么复杂的体系。我现在也来不及检查，你可以先用可视化软件看看，是不是有两个氯离子其实是桥联的；如果是这种情况的话，可以看看 https://gaussian.com/afc/ ，然后你再思考一下，应该能想通的。

chands

Daniel_Arndt 发表于 2022-1-10 12:00
如果是MgO的话，你用默认的初猜，也就是不带电荷的镁原子和不带电荷的氧原子，应该是可以收敛的。你要是 ...

谢谢老师！有问题我再求教。