diff --git a/input/methods b/input/methods index 41e9098..f582e3d 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF GHF ROHF - T F F F + F F T F # MP2 MP3 F F # CCD pCCD DCD CCSD CCSD(T) diff --git a/input/options b/input/options index b430666..54c59ac 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS guess mix shift stab search - 1000 0.0000001 5 1 0.0 0.0 F T + 1000 0.0000001 5 1 0.0 0.0 T T # MP: reg F # CC: maxSCF thresh DIIS diff --git a/src/HF/GHF_search.f90 b/src/HF/GHF_search.f90 index b9ce69d..ce2143c 100644 --- a/src/HF/GHF_search.f90 +++ b/src/HF/GHF_search.f90 @@ -52,6 +52,9 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu double precision,allocatable :: Bph(:,:) double precision,allocatable :: AB(:,:) double precision,allocatable :: Om(:) + double precision,allocatable :: R(:,:) + double precision,allocatable :: ExpR(:,:) + integer :: eig double precision :: kick,step @@ -59,9 +62,9 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Output variables double precision,intent(out) :: EHF - double precision,intent(out) :: e(nBas) - double precision,intent(inout):: c(nBas,nBas) - double precision,intent(out) :: P(nBas,nBas) + double precision,intent(out) :: e(nBas2) + double precision,intent(inout):: c(nBas2,nBas2) + double precision,intent(out) :: P(nBas2,nBas2) ! Memory allocation @@ -77,7 +80,8 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu nS = (nO - nC)*(nV - nR) - allocate(ERI_MO(nBas2,nBas2,nBas2,nBas2),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS)) + allocate(ERI_MO(nBas2,nBas2,nBas2,nBas2),Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS), & + R(nBas2,nBas2),ExpR(nBas2,nBas2)) !------------------! ! Search algorithm ! @@ -145,11 +149,11 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu call phLR_A(ispin,.false.,nBas2,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph) call phLR_B(ispin,.false.,nBas2,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) - - AB(:,:) = Aph(:,:) - Bph(:,:) + + AB(:,:) = Aph(:,:) + Bph(:,:) call diagonalize_matrix(nS,AB,Om) - Om(:) = 0.5d0*Om(:) + Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real GHF -> Real GHF |' @@ -182,18 +186,31 @@ subroutine GHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu step = 1d0 - do mu=1,nBas2 - ia = 0 - do i=nC+1,nO - kick = 0d0 - do a=nO+1,nBas2-nR - ia = ia + 1 - kick = kick + AB(ia,eig)*c(mu,a) - end do - c(mu,i) = c(mu,i) + step*kick +! do mu=1,nBas2 +! ia = 0 +! do i=nC+1,nO +! kick = 0d0 +! do a=nO+1,nBas2-nR +! ia = ia + 1 +! kick = kick + AB(ia,eig)*c(mu,a) +! end do +! c(mu,i) = c(mu,i) + step*kick +! end do +! end do + + R(:,:) = 0d0 + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas2-nR + ia = ia + 1 + R(a,i) = +AB(ia,eig) + R(i,a) = -AB(ia,eig) end do end do + call matrix_exponential(nBas2,R,ExpR) + c = matmul(c,ExpR) + else write(*,'(1X,A40,1X)') 'Well done, GHF solution is stable!' diff --git a/src/HF/GHF_stability.f90 b/src/HF/GHF_stability.f90 index 7c076ec..f993147 100644 --- a/src/HF/GHF_stability.f90 +++ b/src/HF/GHF_stability.f90 @@ -44,7 +44,7 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) AB(:,:) = A(:,:) + B(:,:) call diagonalize_matrix(nS,AB,Om) - Om(:) = 0.5d0*Om(:) + Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real GHF -> Real GHF |' @@ -79,7 +79,7 @@ subroutine GHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) AB(:,:) = A(:,:) - B(:,:) call diagonalize_matrix(nS,AB,Om) - Om(:) = 0.5d0*Om(:) + Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real GHF -> Complex GHF |' diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index 2c00bd2..d27efa4 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -123,10 +123,10 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph) call phLR_B(ispin,.false.,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) - AB(:,:) = Aph(:,:) + Bph(:,:) + AB(:,:) = Aph(:,:) - Bph(:,:) call diagonalize_matrix(nS,AB,Om) - Om(:) = 0.5d0*Om(:) + Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real RHF -> Real RHF |' diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90 index a96ca7c..7a72560 100644 --- a/src/HF/RHF_stability.f90 +++ b/src/HF/RHF_stability.f90 @@ -44,7 +44,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) AB(:,:) = A(:,:) + B(:,:) call diagonalize_matrix(nS,AB,Om) - Om(:) = 0.5d0*Om(:) + Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real RHF -> Real RHF |' @@ -81,7 +81,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) AB(:,:) = A(:,:) - B(:,:) call diagonalize_matrix(nS,AB,Om) - Om(:) = 0.5d0*Om(:) + Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real RHF -> Complex RHF |' @@ -121,7 +121,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) AB(:,:) = A(:,:) + B(:,:) call diagonalize_matrix(nS,AB,Om) - Om(:) = 0.5d0*Om(:) + Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real RHF -> Real UHF |' diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90 index a83e19e..ecbf1ce 100644 --- a/src/HF/UHF_search.f90 +++ b/src/HF/UHF_search.f90 @@ -144,7 +144,7 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu AB_sc(:,:) = A_sc(:,:) + B_sc(:,:) call diagonalize_matrix(nS_sc,AB_sc,Om_sc) - Om_sc(:) = 0.5d0*Om_sc(:) + Om_sc(:) = 2d0*Om_sc(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real UHF -> Real UHF |' diff --git a/src/HF/UHF_stability.f90 b/src/HF/UHF_stability.f90 index c1a3767..0d18e14 100644 --- a/src/HF/UHF_stability.f90 +++ b/src/HF/UHF_stability.f90 @@ -57,7 +57,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) AB_sc(:,:) = A_sc(:,:) + B_sc(:,:) call diagonalize_matrix(nS_sc,AB_sc,Om_sc) - Om_sc(:) = 0.5d0*Om_sc(:) + Om_sc(:) = 2d0*Om_sc(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real UHF -> Real UHF |' @@ -92,7 +92,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) AB_sc(:,:) = A_sc(:,:) - B_sc(:,:) call diagonalize_matrix(nS_sc,AB_sc,Om_sc) - Om_sc(:) = 0.5d0*Om_sc(:) + Om_sc(:) = 2d0*Om_sc(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real UHF -> Complex UHF |' @@ -142,7 +142,7 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) AB_sf(:,:) = A_sf(:,:) + B_sf(:,:) call diagonalize_matrix(nS_sf,AB_sf,Om_sf) - Om_sf(:) = 0.5d0*Om_sf(:) + Om_sf(:) = 2d0*Om_sf(:) write(*,*)'-------------------------------------------------------------' write(*,*)'| Stability analysis: Real UHF -> Real GHF |' diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index f1e26e4..345836c 100644 --- a/src/utils/utils.f90 +++ b/src/utils/utils.f90 @@ -72,6 +72,7 @@ subroutine matrix_exponential(N,A,ExpA) implicit none integer,intent(in) :: N + integer :: i double precision,intent(in) :: A(N,N) double precision,allocatable :: W(:,:) double precision,allocatable :: tau(:) @@ -90,7 +91,13 @@ subroutine matrix_exponential(N,A,ExpA) W(:,:) = - matmul(A,A) call diagonalize_matrix(N,W,tau) - tau(:) = sqrt(tau(:)) + + call matout(N,1,tau) +! do i=1,N +! tau(i) = max(abs(tau(i)),1d-14) +! end do + tau(:) = sqrt(abs(tau(:))) + call matout(N,1,tau) ! Construct cos part