From a4eb01b589d85a7ad0e9a006c7cb154237c8d6b8 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Wed, 28 Aug 2024 19:22:06 +0200 Subject: [PATCH] introduce nBas_MOs in RHF_search --- src/HF/RHF_search.f90 | 67 ++++++++++++++++++++++------------------ src/HF/RHF_stability.f90 | 4 ++- src/QuAcK/RQuAcK.f90 | 8 ++--- src/utils/utils.f90 | 24 +++++++------- 4 files changed, 57 insertions(+), 46 deletions(-) diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index c568fd3..714d11c 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -1,6 +1,9 @@ -subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO, & - X,ERHF,e,c,P) + +! --- + +subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas_AOs, nBas_MOs, nC, nO, nV, nR, S, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, & + X, ERHF, e, c, P) ! Search for RHF solutions @@ -13,7 +16,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,intent(in) :: thresh double precision,intent(in) :: level_shift - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -22,15 +25,15 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: X(nBas,nBas) - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(inout):: ERI_MO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) - double precision,intent(inout):: dipole_int_MO(nBas,nBas,ncart) + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: T(nBas_AOs,nBas_AOs) + double precision,intent(in) :: V(nBas_AOs,nBas_AOs) + double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) + double precision,intent(in) :: X(nBas_AOs,nBas_MOs) + double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) + double precision,intent(inout):: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) + double precision,intent(in) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart) + double precision,intent(inout):: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) ! Local variables @@ -59,9 +62,9 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ! Output variables double precision,intent(out) :: ERHF - 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(nBas_MOs) + double precision,intent(inout):: c(nBas_AOs,nBas_MOs) + double precision,intent(out) :: P(nBas_AOs,nBas_AOs) ! Memory allocation @@ -76,7 +79,8 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN !-------------------! nS = (nO - nC)*(nV - nR) - allocate(Aph(nS,nS),Bph(nS,nS),AB(nS,nS),Om(nS),R(nBas,nBas),ExpR(nBas,nBas)) + allocate(Aph(nS,nS), Bph(nS,nS), AB(nS,nS), Om(nS)) + allocate(R(nBas_MOs,nBas_MOs), ExpR(nBas_MOs,nBas_MOs)) !------------------! ! Search algorithm ! @@ -92,12 +96,12 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN !---------------------! call wall_time(start_HF) - call RHF(.false.,maxSCF,thresh,max_diis,guess,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,e,c,P) + call RHF(.false., maxSCF, thresh, max_diis, guess, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, e, c, P) call wall_time(end_HF) t_HF = end_HF - start_HF - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for RHF = ',t_HF,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for RHF = ',t_HF,' seconds' write(*,*) !----------------------------------! @@ -108,14 +112,14 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN write(*,*) write(*,*) 'AO to MO transformation... Please be patient' write(*,*) - do ixyz=1,ncart - call AOtoMO(nBas,nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + do ixyz = 1, ncart + call AOtoMO(nBas_AOs, nBas_MOs, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas,nBas,c,ERI_AO,ERI_MO) + call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) call wall_time(end_AOtoMO) t_AOtoMO = end_AOtoMO - start_AOtoMO - write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for AO to MO transformation = ',t_AOtoMO,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total wall time for AO to MO transformation = ',t_AOtoMO,' seconds' write(*,*) !-------------------------------------------------------------! @@ -124,12 +128,12 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN ispin = 1 - 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) + call phLR_A(ispin,.false.,nBas_MOs,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph) + call phLR_B(ispin,.false.,nBas_MOs,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) AB(:,:) = Aph(:,:) + Bph(:,:) - call diagonalize_matrix(nS,AB,Om) + call diagonalize_matrix(nS, AB, Om) Om(:) = 2d0*Om(:) write(*,*)'-------------------------------------------------------------' @@ -156,6 +160,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN if(eig < 0 .or. eig > nS) then write(*,'(1X,A40,1X,A10)') 'Invalid option...','Stop...' write(*,*) + deallocate(Aph, Bph, AB, Om, R, ExpR) stop end if @@ -164,15 +169,15 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN R(:,:) = 0d0 ia = 0 do i=nC+1,nO - do a=nO+1,nBas-nR + do a=nO+1,nBas_MOs-nR ia = ia + 1 R(a,i) = +AB(ia,eig) R(i,a) = -AB(ia,eig) end do end do - call matrix_exponential(nBas,R,ExpR) - c = matmul(c,ExpR) + call matrix_exponential(nBas_MOs, R, ExpR) + c = matmul(c, ExpR) else @@ -191,4 +196,6 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN !---------------! end do + deallocate(Aph, Bph, AB, Om, R, ExpR) + end subroutine diff --git a/src/HF/RHF_stability.f90 b/src/HF/RHF_stability.f90 index d7b5385..154a4f9 100644 --- a/src/HF/RHF_stability.f90 +++ b/src/HF/RHF_stability.f90 @@ -30,7 +30,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) ! Memory allocation - allocate(A(nS,nS),B(nS,nS),AB(nS,nS),Om(nS)) + allocate(A(nS,nS), B(nS,nS), AB(nS,nS), Om(nS)) !-------------------------------------------------------------! ! Stability analysis: Real RHF -> Real RHF @@ -148,5 +148,7 @@ subroutine RHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI) end if write(*,*)'-------------------------------------------------------------' write(*,*) + + deallocate(A, B, AB, Om) end subroutine diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 26f817a..18fc00a 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -180,7 +180,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(dostab) then call wall_time(start_stab) - call RHF_stability(nBas_AOs,nC,nO,nV,nR,nS,eHF,ERI_MO) + call RHF_stability(nBas_MOs, nC, nO, nV, nR, nS, eHF, ERI_MO) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -192,9 +192,9 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(dosearch) then call wall_time(start_stab) - ! TODO - call RHF_search(maxSCF_HF,thresh_HF,max_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas_AOs,nC,nO,nV,nR,S,T,V,Hc,ERI_AO,ERI_MO,dipole_int_AO,dipole_int_MO,X,ERHF,eHF,cHF,PHF) + call RHF_search(maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas_AOs, nBas_MOs, nC, nO, nV, nR, S, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, & + dipole_int_MO, X, ERHF, eHF, cHF, PHF) call wall_time(end_stab) t_stab = end_stab - start_stab diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index 80b8322..316a87c 100644 --- a/src/utils/utils.f90 +++ b/src/utils/utils.f90 @@ -65,7 +65,7 @@ subroutine diagonal_matrix(N,D,A) end subroutine !------------------------------------------------------------------------ -subroutine matrix_exponential(N,A,ExpA) +subroutine matrix_exponential(N, A, ExpA) ! Compute Exp(A) @@ -81,7 +81,7 @@ subroutine matrix_exponential(N,A,ExpA) ! Memory allocation - allocate(W(N,N),tau(N),t(N,N)) + allocate(W(N,N), tau(N), t(N,N)) ! Initialize @@ -89,8 +89,8 @@ subroutine matrix_exponential(N,A,ExpA) ! Diagonalize - W(:,:) = - matmul(A,A) - call diagonalize_matrix(N,W,tau) + W(:,:) = - matmul(A, A) + call diagonalize_matrix(N, W, tau) ! do i=1,N ! tau(i) = max(abs(tau(i)),1d-14) @@ -99,16 +99,18 @@ subroutine matrix_exponential(N,A,ExpA) ! Construct cos part - call diagonal_matrix(N,cos(tau),t) - t(:,:) = matmul(t,transpose(W)) - ExpA(:,:) = ExpA(:,:) + matmul(W,t) + call diagonal_matrix(N, cos(tau), t) + t(:,:) = matmul(t, transpose(W)) + ExpA(:,:) = ExpA(:,:) + matmul(W, t) ! Construct sin part - call diagonal_matrix(N,sin(tau)/tau,t) - t(:,:) = matmul(t,transpose(W)) - t(:,:) = matmul(t,A) - ExpA(:,:) = ExpA(:,:) + matmul(W,t) + call diagonal_matrix(N, sin(tau)/tau, t) + t(:,:) = matmul(t, transpose(W)) + t(:,:) = matmul(t, A) + ExpA(:,:) = ExpA(:,:) + matmul(W, t) + + deallocate(W, tau, t) end subroutine