RASSI模块EJOB关键词在新版本OpenMolcas-v22.02的改动

叫个啥名字发布于 2022-4-20 21:06 阅读：348 linux

发现过程：

整理数据发现，用新版本一直重复不了旧版本的光谱图，一度以为自己参数改了，找了半天发现毛都没变。所有结果检查以后发现只有版本和编译时用的的编译器不同，换成同版本编译器结果依旧不同，最后老版本算了一遍，结果和原来一样。

RASSI模块的变化：

新版本的变化主要是添加了新功能，比如NTO，SNTO等等，其中还有一个关键词描述是不同的。

旧版本v19.11

新版本v22.02

EJOB

The spin-free effective Hamiltonian is assumed to be diagonal, with energies being read from a JOBIPH or JOBMIX file. If this keyword is used together with HEFF, or if the input file is an HDF5 file for which the effective Hamiltonian is automatically read, only the diagonal elements will be read and off-diagonal elements will be set to zero. This can be useful to use the SS-CASPT2 energies from a MS-CASTP2 calculation.

EJOB

The spin-free effective Hamiltonian’s diagonal is filled with energies read from a JOBIPH or JOBMIX file. If an effective Hamiltonian is read (using HEFF or reading from an HDF5 file), the diagonal elements are taken from the stored Hamiltonian; this can be useful for using the SS-CASPT2 energies from a MS-CASTP2 calculation. The off-diagonal elements are approximated as H_ij~ 0.5 * S_ij(H_ii+H_ij) , where ??? is the overlap between two states; so if the input states are orthogonal, the effective Hamiltonian will be diagonal.

按照描述的意思是旧版本读取哈密顿量只读取对角元，非对角元全部设置0。新版本则是将非对角元近似的给了一个数值。

这些改变在软件的源代码中的体现在两个文件中的不同，gtdmctl.f以及eigctl.f

gtdmctl.f中修改了对角元的表达形式，eigctl.f修改了循环的位置。

gtdmctl.f 修改
新比旧删除了：
          IF(IFHAM.AND..NOT.(IFHEXT.or.IFHEFF.or.IFEJOB))THEN
            HZERO              = ECORE*SIJ
            HIJ                = HZERO+HONE+HTWO
            HAM(ISTATE,JSTATE) = HIJ
            HAM(JSTATE,ISTATE) = HIJ
添加了：
        IF (IFEJOB.and.(ISTATE.ne.JSTATE)) THEN 
          HAM(ISTATE,JSTATE) = SIJ  
          HAM(JSTATE,ISTATE) = SIJ
        END IF                 

  IF (IFEJOB) THEN
      DO JST=1,NSTAT(JOB2)
         JSTATE=ISTAT(JOB2)-1+JST
         DO IST=1,NSTAT(JOB1)
           ISTATE=ISTAT(JOB1)-1+IST
           IF(ISTATE.LE.JSTATE) CYCLE
          SIJ=HAM(ISTATE,JSTATE)
           HII=HAM(ISTATE,ISTATE)
           HJJ=HAM(JSTATE,JSTATE)
           HAM(ISTATE,JSTATE) = SIJ*(HII+HJJ)*0.5D0
          HAM(JSTATE,ISTATE) = SIJ*(HII+HJJ)*0.5D0
         END DO
       END DO
     END IF

eigctl.f修改
新比旧删除了：
      IJ=0
      DO I=1,MSTATE
        DO J=1,I
          IJ=IJ+1
          If (I.NE.J .AND.
     &    ABS(WORK(LHSQ-1+I+MSTATE*(J-1))).gt.1.0D-9) Then
            DIAGONAL=.FALSE.
            go to 11
          EndIf
        End Do
      End Do
 11   Continue
添加了：
         If (I.NE.J .AND.
      &    (ABS(ovlp(i,j)).gt.1.0D-9.or.ABS(ham(i,j)).gt.1.0D-9)) Then
             DIAGONAL=.FALSE.
           EndIf
并且循环的位置改了。详细变化可以去查看两个版本的源代码位于/src/rassi

开发者回复：

追踪看了下这个EJOB code，在19年三月份的时候更新过一次，当时投文章遇到审稿人提出相关问题，为了该文章做过一次改动（没有根本性解决问题，因为文章中的计算是一个特例）。后来的改动就是你现在看到的是在2020年五月之后的更新。之前的EJOB code assumes Hamitonian to be diagonal butnon-diagonal overlap,所以ejob code 会diagonalize the overlap，从而有时候描述Hamiltonian会很差。最准确的方式解决orthonormality 是CASSCF off-diagonal Hamitonia + (0.5(E_PT2(I) + E_PT2(J))* SIJ)），这个就是算法非常费时，后面的这次改动是对off-diagonal energy 采用校正 (0.5(E_total(I) + E_total(J))* SIJ)。至于对property计算的影响，主要体现在Hard x-ray需要计算multiple expansion，比如second order expansion对non-orthonarmality比较敏感。其他的例子目前就不大清楚，也没有别的人report bug。

This makes it feel like the idea is that the Hamiltonian (CASSCF level) is computed completely and the PT2 correction is added to diagonal and off-diagonal terms Such an implementation would thus tolerate non-orthonormality without much issue. If you think this is correct, it would take about 2 min fixing it in the code to try it, since it is just a flag to activate or not.

结论：

虽然将Ejob关键词的相关改动部分改回来，新旧结果相同。但是新版本的是更准确且合理的。

OpenMolcas