Merge pull request #223 from Ydrnan/qp2_add

Change coupling parameters ROHF
2024-12-22 12:23:43 +01:00 · 2022-11-29 15:28:33 +01:00 · 2022-11-29 15:28:33 +01:00 · 0a0b266eb7
commit 0a0b266eb7
parent 21bbc6911c 0b27f42c96
1 changed files with 138 additions and 54 deletions
--- a/src/scf_utils/fock_matrix.irp.f
+++ b/src/scf_utils/fock_matrix.irp.f
@ -5,6 +5,90 @@
   ! Fock matrix on the MO basis.
   ! For open shells, the ROHF Fock Matrix is ::
   !
+   !       |  Rcc  |  F^b  |  Fcv  |
+   !       |-----------------------|
+   !       |  F^b  |  Roo  |  F^a  |
+   !       |-----------------------|
+   !       |  Fcv  |  F^a  |  Rvv  |
+   !
+   ! C: Core, O: Open, V: Virtual 
+   ! 
+   ! Rcc = Acc Fcc^a + Bcc Fcc^b
+   ! Roo = Aoo Foo^a + Boo Foo^b
+   ! Rcc = Avv Fvv^a + Bvv Fvv^b
+   ! Fcv = (F^a + F^b)/2
+   ! 
+   ! F^a: Fock matrix alpha (MO), F^b: Fock matrix beta (MO)
+   ! A,B: Coupling parameters
+   !
+   ! J. Chem. Phys. 133, 141102 (2010), https://doi.org/10.1063/1.3503173
+   ! Coupling parameters from J. Chem. Phys. 125, 204110 (2006); https://doi.org/10.1063/1.2393223.
+   !       cc   oo   vv
+   !  A  -0.5  0.5  1.5
+   !  B   1.5  0.5 -0.5
+   ! 
+   END_DOC
+   integer                        :: i,j,n
+   if (elec_alpha_num == elec_beta_num) then
+     Fock_matrix_mo = Fock_matrix_mo_alpha
+   else
+     ! Core
+     do j = 1, elec_beta_num
+       ! Core
+       do i = 1, elec_beta_num
+         fock_matrix_mo(i,j) = - 0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 1.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_beta(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         fock_matrix_mo(i,j) =   0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+     enddo
+     ! Open
+     do j = elec_beta_num+1, elec_alpha_num
+       ! Core
+       do i = 1, elec_beta_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_beta(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         fock_matrix_mo(i,j) =   0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_alpha(i,j)
+       enddo
+     enddo
+     ! Virtual
+     do j = elec_alpha_num+1, mo_num
+       ! Core
+       do i = 1, elec_beta_num
+         fock_matrix_mo(i,j) =   0.5d0 * fock_matrix_mo_alpha(i,j) &
+                               + 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+       ! Open
+       do i = elec_beta_num+1, elec_alpha_num
+         fock_matrix_mo(i,j) = fock_matrix_mo_alpha(i,j)
+       enddo
+       ! Virtual
+       do i = elec_alpha_num+1, mo_num
+         fock_matrix_mo(i,j) =   1.5d0 * fock_matrix_mo_alpha(i,j) &
+                               - 0.5d0 * fock_matrix_mo_beta(i,j)
+       enddo
+     enddo
+   endif
+
+   ! Old
+   ! BEGIN_DOC
+   ! Fock matrix on the MO basis.
+   ! For open shells, the ROHF Fock Matrix is ::
+   !
   !       |   F-K    |  F + K/2  |    F     |
   !       |---------------------------------|
   !       | F + K/2  |     F     |  F - K/2 |
@ -16,64 +100,64 @@
   !
   ! K = Fb - Fa
   !
-   END_DOC
-   integer                        :: i,j,n
-   if (elec_alpha_num == elec_beta_num) then
-     Fock_matrix_mo = Fock_matrix_mo_alpha
-   else
+   ! END_DOC
+   !integer                        :: i,j,n
+   !if (elec_alpha_num == elec_beta_num) then
+   !  Fock_matrix_mo = Fock_matrix_mo_alpha
+   !else

-     do j=1,elec_beta_num
-       ! F-K
-       do i=1,elec_beta_num !CC
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
-       ! F+K/2
-       do i=elec_beta_num+1,elec_alpha_num  !CA
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
-       ! F
-       do i=elec_alpha_num+1, mo_num !CV
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
-       enddo
-     enddo
+   !  do j=1,elec_beta_num
+   !    ! F-K
+   !    do i=1,elec_beta_num !CC
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !          - (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !    enddo
+   !    ! F+K/2
+   !    do i=elec_beta_num+1,elec_alpha_num  !CA
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !          + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !    enddo
+   !    ! F
+   !    do i=elec_alpha_num+1, mo_num !CV
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
+   !    enddo
+   !  enddo

-     do j=elec_beta_num+1,elec_alpha_num
-       ! F+K/2
-       do i=1,elec_beta_num !AC
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
-       ! F
-       do i=elec_beta_num+1,elec_alpha_num !AA
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
-       enddo
-       ! F-K/2
-       do i=elec_alpha_num+1, mo_num !AV
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
-     enddo
+   !  do j=elec_beta_num+1,elec_alpha_num
+   !    ! F+K/2
+   !    do i=1,elec_beta_num !AC
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !          + 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !    enddo
+   !    ! F
+   !    do i=elec_beta_num+1,elec_alpha_num !AA
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
+   !    enddo
+   !    ! F-K/2
+   !    do i=elec_alpha_num+1, mo_num !AV
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !          - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !    enddo
+   !  enddo

-     do j=elec_alpha_num+1, mo_num
-       ! F
-       do i=1,elec_beta_num !VC
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
-       enddo
-       ! F-K/2
-       do i=elec_beta_num+1,elec_alpha_num !VA
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
-             - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
-       ! F+K
-       do i=elec_alpha_num+1,mo_num !VV
-         Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) &
-             + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
-       enddo
-     enddo
+   !  do j=elec_alpha_num+1, mo_num
+   !    ! F
+   !    do i=1,elec_beta_num !VC
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))
+   !    enddo
+   !    ! F-K/2
+   !    do i=elec_beta_num+1,elec_alpha_num !VA
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j))&
+   !          - 0.5d0*(Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !    enddo
+   !    ! F+K
+   !    do i=elec_alpha_num+1,mo_num !VV
+   !      Fock_matrix_mo(i,j) = 0.5d0*(Fock_matrix_mo_alpha(i,j)+Fock_matrix_mo_beta(i,j)) &
+   !          + (Fock_matrix_mo_beta(i,j) - Fock_matrix_mo_alpha(i,j))
+   !    enddo
+   !  enddo

-   endif
+   !endif

   do i = 1, mo_num
     Fock_matrix_diag_mo(i) = Fock_matrix_mo(i,i)