From d823cdcd20a4041b76bdf0d31ee2a0ef0cc76fb8 Mon Sep 17 00:00:00 2001 From: Abdallah Ammar Date: Wed, 28 Aug 2024 18:39:51 +0200 Subject: [PATCH] introduce nBas_MOs in ROHF --- src/AOtoMO/AOtoMO.f90 | 29 +++++++--- src/AOtoMO/MOtoAO.f90 | 39 +++++++++---- src/GF/qsRGF2.f90 | 2 +- src/GF/qsUGF2.f90 | 2 +- src/GT/qsRGTeh.f90 | 2 +- src/GT/qsRGTpp.f90 | 2 +- src/GT/qsUGTpp.f90 | 2 +- src/GW/SRG_qsGW.f90 | 6 +- src/GW/SRG_qsUGW.f90 | 8 +-- src/GW/qsRGW.f90 | 6 +- src/GW/qsUGW.f90 | 8 +-- src/HF/RHF_search.f90 | 2 +- src/HF/ROHF.f90 | 111 +++++++++++++++++++++--------------- src/HF/ROHF_fock_matrix.f90 | 35 +++++++----- src/HF/UHF_search.f90 | 4 +- src/HF/print_ROHF.f90 | 17 +++--- src/QuAcK/RQuAcK.f90 | 9 ++- src/QuAcK/UQuAcK.f90 | 4 +- 18 files changed, 171 insertions(+), 117 deletions(-) diff --git a/src/AOtoMO/AOtoMO.f90 b/src/AOtoMO/AOtoMO.f90 index 2192bc7..8c5ce8f 100644 --- a/src/AOtoMO/AOtoMO.f90 +++ b/src/AOtoMO/AOtoMO.f90 @@ -1,14 +1,15 @@ -subroutine AOtoMO(nBas,C,A,B) +subroutine AOtoMO(nBas_AOs, nBas_MOs, C, M_AOs, M_MOs) -! Perform AO to MO transformation of a matrix A for given coefficients c +! Perform AO to MO transformation of a matrix M_AOs for given coefficients c +! M_MOs = C.T M_AOs C implicit none ! Input variables - integer,intent(in) :: nBas - double precision,intent(in) :: C(nBas,nBas) - double precision,intent(in) :: A(nBas,nBas) + integer,intent(in) :: nBas_AOs, nBas_MOs + double precision,intent(in) :: C(nBas_AOs,nBas_MOs) + double precision,intent(in) :: M_AOs(nBas_AOs,nBas_AOs) ! Local variables @@ -16,11 +17,21 @@ subroutine AOtoMO(nBas,C,A,B) ! Output variables - double precision,intent(out) :: B(nBas,nBas) + double precision,intent(out) :: M_MOs(nBas_MOs,nBas_MOs) - allocate(AC(nBas,nBas)) + allocate(AC(nBas_AOs,nBas_MOs)) - AC = matmul(A,C) - B = matmul(transpose(C),AC) + !AC = matmul(M_AOs, C) + !M_MOs = matmul(transpose(C), AC) + + call dgemm("N", "N", nBas_AOs, nBas_MOs, nBas_AOs, 1.d0, & + M_AOs(1,1), nBas_AOs, C(1,1), nBas_AOs, & + 0.d0, AC(1,1), nBas_AOs) + + call dgemm("T", "N", nBas_MOs, nBas_MOs, nBas_AOs, 1.d0, & + C(1,1), nBas_AOs, AC(1,1), nBas_AOs, & + 0.d0, M_MOs(1,1), nBas_MOs) + + deallocate(AC) end subroutine diff --git a/src/AOtoMO/MOtoAO.f90 b/src/AOtoMO/MOtoAO.f90 index faffa92..06abb38 100644 --- a/src/AOtoMO/MOtoAO.f90 +++ b/src/AOtoMO/MOtoAO.f90 @@ -1,16 +1,19 @@ -subroutine MOtoAO(nBas,S,C,B,A) +subroutine MOtoAO(nBas_AOs, nBas_MOs, S, C, M_MOs, M_AOs) -! Perform MO to AO transformation of a matrix A for a given metric S +! Perform MO to AO transformation of a matrix M_AOs for a given metric S ! and coefficients c +! +! M_AOs = S C M_MOs (S C).T +! implicit none ! Input variables - integer,intent(in) :: nBas - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: C(nBas,nBas) - double precision,intent(in) :: B(nBas,nBas) + integer,intent(in) :: nBas_AOs, nBas_MOs + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: C(nBas_AOs,nBas_MOs) + double precision,intent(in) :: M_MOs(nBas_MOs,nBas_MOs) ! Local variables @@ -18,14 +21,28 @@ subroutine MOtoAO(nBas,S,C,B,A) ! Output variables - double precision,intent(out) :: A(nBas,nBas) + double precision,intent(out) :: M_AOs(nBas_AOs,nBas_AOs) ! Memory allocation - allocate(SC(nBas,nBas),BSC(nBas,nBas)) + allocate(SC(nBas_AOs,nBas_MOs), BSC(nBas_MOs,nBas_AOs)) - SC = matmul(S,C) - BSC = matmul(B,transpose(SC)) - A = matmul(SC,BSC) + !SC = matmul(S, C) + !BSC = matmul(M_MOs, transpose(SC)) + !M_AOs = matmul(SC, BSC) + + call dgemm("N", "N", nBas_AOs, nBas_MOs, nBas_AOs, 0.d0, & + S(1,1), nBas_AOs, C(1,1), nBas_AOs, & + 1.d0, SC(1,1), nBas_AOs) + + call dgemm("N", "T", nBas_MOs, nBas_AOs, nBas_MOs, 0.d0, & + M_MOs(1,1), nBas_MOs, SC(1,1), nBas_AOs, & + 1.d0, BSC(1,1), nBas_MOs) + + call dgemm("N", "N", nBas_AOs, nBas_AOs, nBas_MOs, 0.d0, & + SC(1,1), nBas_AOs, BSC(1,1), nBas_MOs, & + 1.d0, M_AOs(1,1), nBas_AOs) + + deallocate(SC, BSC) end subroutine diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90 index 8222e56..cb744c9 100644 --- a/src/GF/qsRGF2.f90 +++ b/src/GF/qsRGF2.f90 @@ -161,7 +161,7 @@ subroutine qsRGF2(dotest,maxSCF,thresh,max_diis,dophBSE,doppBSE,TDA,dBSE,dTDA,si SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO(nBas,S,c,SigC,SigCp) + call MOtoAO(nBas,nBas,S,c,SigC,SigCp) ! Solve the quasi-particle equation diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index 6fa0bdb..2aedac1 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -197,7 +197,7 @@ subroutine qsUGF2(dotest,maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,spin_conserved end do do is=1,nspin - call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) + call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation diff --git a/src/GT/qsRGTeh.f90 b/src/GT/qsRGTeh.f90 index b814001..60f8d1b 100644 --- a/src/GT/qsRGTeh.f90 +++ b/src/GT/qsRGTeh.f90 @@ -192,7 +192,7 @@ subroutine qsRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas,S,c,Sig,Sigp) + call MOtoAO(nBas,nBas,S,c,Sig,Sigp) ! Solve the quasi-particle equation diff --git a/src/GT/qsRGTpp.f90 b/src/GT/qsRGTpp.f90 index 2ff7165..9cb12da 100644 --- a/src/GT/qsRGTpp.f90 +++ b/src/GT/qsRGTpp.f90 @@ -235,7 +235,7 @@ subroutine qsRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,d Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas,S,c,Sig,Sigp) + call MOtoAO(nBas,nBas,S,c,Sig,Sigp) ! Solve the quasi-particle equation diff --git a/src/GT/qsUGTpp.f90 b/src/GT/qsUGTpp.f90 index 5f539bb..81c1cda 100644 --- a/src/GT/qsUGTpp.f90 +++ b/src/GT/qsUGTpp.f90 @@ -280,7 +280,7 @@ subroutine qsUGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,B end do do ispin=1,nspin - call MOtoAO(nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin)) + call MOtoAO(nBas,nBas,S,c(:,:,ispin),SigT(:,:,ispin),SigTp(:,:,ispin)) end do ! Solve the quasi-particle equation diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 9fbbd8c..84c953f 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -175,7 +175,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call wall_time(tao1) do ixyz=1,ncart - call AOtoMO(nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + call AOtoMO(nBas,nBas,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) @@ -218,7 +218,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, ! Make correlation self-energy Hermitian and transform it back to AO basis call wall_time(tmo1) - call MOtoAO(nBas,S,c,SigC,SigCp) + call MOtoAO(nBas,nBas,S,c,SigC,SigCp) call wall_time(tmo2) tmo = tmo + tmo2 - tmo1 ! Solve the quasi-particle equation @@ -245,7 +245,7 @@ subroutine SRG_qsGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, call diagonalize_matrix(nBas,cp,eGW) c = matmul(X,cp) - call AOtoMO(nBas,c,SigCp,SigC) + call AOtoMO(nBas,nBas,c,SigCp,SigC) ! Compute new density matrix in the AO basis diff --git a/src/GW/SRG_qsUGW.f90 b/src/GW/SRG_qsUGW.f90 index 0705847..6939afb 100644 --- a/src/GW/SRG_qsUGW.f90 +++ b/src/GW/SRG_qsUGW.f90 @@ -184,8 +184,8 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS !-------------------------------------------------- do ixyz=1,ncart - call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block @@ -228,7 +228,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Make correlation self-energy Hermitian and transform it back to AO basis do is=1,nspin - call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) + call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation @@ -279,7 +279,7 @@ subroutine SRG_qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS ! Back-transform self-energy do is=1,nspin - call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) + call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) end do ! Compute density matrix diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index 6e10f75..a0df0e5 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -166,7 +166,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop ! AO to MO transformation of two-electron integrals do ixyz=1,ncart - call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + call AOtoMO(nBas,nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) @@ -191,7 +191,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO(nBas,S,c,SigC,SigCp) + call MOtoAO(nBas,nBas,S,c,SigC,SigCp) ! Solve the quasi-particle equation @@ -238,7 +238,7 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eGW) c = matmul(X,cp) - call AOtoMO(nBas,c,SigCp,SigC) + call AOtoMO(nBas,nBas,c,SigCp,SigC) ! Density matrix diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index ea57d9e..69e3f11 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -184,8 +184,8 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE !-------------------------------------------------- do ixyz=1,ncart - call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block @@ -232,7 +232,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE end do do is=1,nspin - call MOtoAO(nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) + call MOtoAO(nBas,nBas,S,c(:,:,is),SigC(:,:,is),SigCp(:,:,is)) end do ! Solve the quasi-particle equation @@ -283,7 +283,7 @@ subroutine qsUGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE ! Back-transform self-energy do is=1,nspin - call AOtoMO(nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) + call AOtoMO(nBas,nBas,c(:,:,is),SigCp(:,:,is),SigC(:,:,is)) end do ! Compute density matrix diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index b0bdb30..b28469f 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -109,7 +109,7 @@ subroutine RHF_search(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN write(*,*) 'AO to MO transformation... Please be patient' write(*,*) do ixyz=1,ncart - call AOtoMO(nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + call AOtoMO(nBas,nBas,c,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do call AOtoMO_ERI_RHF(nBas,c,ERI_AO,ERI_MO) call wall_time(end_AOtoMO) diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index 0cd1f5f..d67fc0b 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -1,5 +1,8 @@ -subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EROHF,eHF,c,Ptot) + +! --- + +subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI, dipole_int, X, EROHF, eHF, c, Ptot) ! Perform restricted open-shell Hartree-Fock calculation @@ -16,7 +19,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN double precision,intent(in) :: mix double precision,intent(in) :: level_shift double precision,intent(in) :: thresh - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -24,18 +27,18 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN double precision,intent(in) :: ENuc integer,intent(in) :: nO(nspin) - 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(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(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(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) + double precision,intent(in) :: dipole_int(nBas_AOs,nBas_AOs,ncart) ! Local variables integer :: nSCF - integer :: nBasSq + integer :: nBas_AOs_Sq integer :: n_diis double precision :: Conv double precision :: rcond @@ -62,9 +65,9 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Output variables double precision,intent(out) :: EROHF - double precision,intent(out) :: eHF(nBas) - double precision,intent(inout):: c(nBas,nBas) - double precision,intent(out) :: Ptot(nBas,nBas) + double precision,intent(out) :: eHF(nBas_MOs) + double precision,intent(inout):: c(nBas_AOs,nBas_MOs) + double precision,intent(out) :: Ptot(nBas_AOs,nBas_AOs) ! Hello world @@ -76,19 +79,30 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Useful stuff - nBasSq = nBas*nBas + nBas_AOs_Sq = nBas_AOs*nBas_AOs ! Memory allocation - allocate(J(nBas,nBas,nspin),F(nBas,nBas,nspin),Fp(nBas,nBas),Ftot(nBas,nBas), & - P(nBas,nBas,nspin),K(nBas,nBas,nspin),err(nBas,nBas),cp(nBas,nBas), & - err_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) + allocate(J(nBas_AOs,nBas_AOs,nspin)) + allocate(K(nBas_AOs,nBas_AOs,nspin)) + allocate(F(nBas_AOs,nBas_AOs,nspin)) + allocate(Ftot(nBas_AOs,nBas_AOs)) + allocate(P(nBas_AOs,nBas_AOs,nspin)) + allocate(err(nBas_AOs,nBas_AOs)) + + allocate(Fp(nBas_MOs,nBas_MOs)) + allocate(cp(nBas_MOs,nBas_MOs)) + + allocate(err_diis(nBas_AOs_Sq,max_diis)) + allocate(F_diis(nBas_AOs_Sq,max_diis)) ! Guess coefficients and demsity matrices - call mo_guess(nBas,nBas,guess_type,S,Hc,X,c) - do ispin=1,nspin - P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin)))) + call mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) + + do ispin = 1, nspin + !P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin)))) + call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO(ispin), 1.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P(1,1,ispin), nBas_AOs) end do Ptot(:,:) = P(:,:,1) + P(:,:,2) @@ -120,51 +134,51 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN ! Build Hartree repulsion - do ispin=1,nspin - call Hartree_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),J(:,:,ispin)) + do ispin = 1, nspin + call Hartree_matrix_AO_basis(nBas_AOs, P(:,:,ispin), ERI(:,:,:,:), J(:,:,ispin)) end do ! Compute exchange potential - do ispin=1,nspin - call exchange_matrix_AO_basis(nBas,P(:,:,ispin),ERI(:,:,:,:),K(:,:,ispin)) + do ispin = 1, nspin + call exchange_matrix_AO_basis(nBas_AOs, P(:,:,ispin), ERI(:,:,:,:), K(:,:,ispin)) end do ! Build Fock operator - do ispin=1,nspin + do ispin = 1, nspin F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) end do - call ROHF_fock_matrix(nBas,nO(1),nO(2),S,c,F(:,:,1),F(:,:,2),Ftot) + call ROHF_fock_matrix(nBas_AOs, nBas_MOs, nO(1), nO(2), S, c, F(:,:,1), F(:,:,2), Ftot) ! Check convergence - err(:,:) = matmul(Ftot,matmul(Ptot,S)) - matmul(matmul(S,Ptot),Ftot) + err(:,:) = matmul(Ftot, matmul(Ptot, S)) - matmul(matmul(S, Ptot), Ftot) if(nSCF > 1) Conv = maxval(abs(err(:,:))) ! Kinetic energy - do ispin=1,nspin - ET(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),T(:,:))) + do ispin = 1, nspin + ET(ispin) = trace_matrix(nBas_AOs, matmul(P(:,:,ispin), T(:,:))) end do ! Potential energy - do ispin=1,nspin - EV(ispin) = trace_matrix(nBas,matmul(P(:,:,ispin),V(:,:))) + do ispin = 1, nspin + EV(ispin) = trace_matrix(nBas_AOs, matmul(P(:,:,ispin), V(:,:))) end do ! Hartree energy - EJ(1) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,1),J(:,:,1))) - EJ(2) = trace_matrix(nBas,matmul(P(:,:,1),J(:,:,2))) - EJ(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) + EJ(1) = 0.5d0*trace_matrix(nBas_AOs, matmul(P(:,:,1), J(:,:,1))) + EJ(2) = trace_matrix(nBas_AOs, matmul(P(:,:,1), J(:,:,2))) + EJ(3) = 0.5d0*trace_matrix(nBas_AOs, matmul(P(:,:,2), J(:,:,2))) ! Exchange energy - do ispin=1,nspin - EK(ispin) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,ispin),K(:,:,ispin))) + do ispin = 1, nspin + EK(ispin) = 0.5d0*trace_matrix(nBas_AOs, matmul(P(:,:,ispin), K(:,:,ispin))) end do ! Total energy @@ -176,7 +190,7 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,Ftot) + call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,err_diis,F_diis,err,Ftot) end if @@ -185,28 +199,29 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(level_shift,nBas,nBas,maxval(nO),S,c,Ftot) + call level_shifting(level_shift, nBas_AOs, nBas_MOs, maxval(nO), S, c, Ftot) end do end if ! Transform Fock matrix in orthogonal basis - Fp(:,:) = matmul(transpose(X(:,:)),matmul(Ftot(:,:),X(:,:))) + Fp(:,:) = matmul(transpose(X(:,:)), matmul(Ftot(:,:), X(:,:))) ! Diagonalize Fock matrix to get eigenvectors and eigenvalues cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eHF) + call diagonalize_matrix(nBas_MOs, cp, eHF) ! Back-transform eigenvectors in non-orthogonal basis - c(:,:) = matmul(X(:,:),cp(:,:)) + c(:,:) = matmul(X(:,:), cp(:,:)) ! Compute density matrix - do ispin=1,nspin - P(:,:,ispin) = matmul(c(:,1:nO(ispin)),transpose(c(:,1:nO(ispin)))) + do ispin = 1, nspin + !P(:,:,ispin) = matmul(c(:,1:nO(ispin)), transpose(c(:,1:nO(ispin)))) + call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO(ispin), 1.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P(1,1,ispin), nBas_AOs) end do Ptot(:,:) = P(:,:,1) + P(:,:,2) @@ -231,14 +246,16 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' write(*,*) + deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) + stop end if ! Compute final UHF energy - call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) - call print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,EK,EROHF,dipole) + call dipole_moment(nBas_AOs,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_ROHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, EK, EROHF, dipole) ! Print test values @@ -248,4 +265,6 @@ subroutine ROHF(dotest,maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZN end if + deallocate(J, K, F, Ftot, P, err, Fp, cp, err_diis, F_diis) + end subroutine diff --git a/src/HF/ROHF_fock_matrix.f90 b/src/HF/ROHF_fock_matrix.f90 index bb23662..a3c5aad 100644 --- a/src/HF/ROHF_fock_matrix.f90 +++ b/src/HF/ROHF_fock_matrix.f90 @@ -1,4 +1,7 @@ -subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) + +! --- + +subroutine ROHF_fock_matrix(nBas_AOs, nBas_MOs, nOa, nOb, S, c, FaAO, FbAO, FAO) ! Construct the ROHF Fock matrix in the AO basis ! For open shells, the ROHF Fock matrix in the MO basis reads @@ -17,14 +20,14 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) ! Input variables - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nOa integer,intent(in) :: nOb - double precision,intent(in) :: S(nBas,nBas) - double precision,intent(in) :: c(nBas,nBas) - double precision,intent(inout):: FaAO(nBas,nBas) - double precision,intent(inout):: FbAO(nBas,nBas) + double precision,intent(in) :: S(nBas_AOs,nBas_AOs) + double precision,intent(in) :: c(nBas_AOs,nBas_MOs) + double precision,intent(inout):: FaAO(nBas_AOs,nBas_AOs) + double precision,intent(inout):: FbAO(nBas_AOs,nBas_AOs) ! Local variables @@ -42,11 +45,11 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) ! Output variables - double precision,intent(out) :: FAO(nBas,nBas) + double precision,intent(out) :: FAO(nBas_AOs,nBas_AOs) ! Memory allocation - allocate(F(nBas,nBas),Fa(nBas,nBas),Fb(nBas,nBas)) + allocate(F(nBas_MOs,nBas_MOs), Fa(nBas_MOs,nBas_MOs), Fb(nBas_MOs,nBas_MOs)) ! Roothan canonicalization parameters @@ -61,14 +64,14 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) ! Number of closed, open, and virtual orbitals - nC = min(nOa,nOb) + nC = min(nOa, nOb) nO = abs(nOa - nOb) - nV = nBas - nC - nO + nV = nBas_AOs - nC - nO ! Block-by-block Fock matrix - call AOtoMO(nBas,c,FaAO,Fa) - call AOtoMO(nBas,c,FbAO,Fb) + call AOtoMO(nBas_AOs, nBas_MOs, c, FaAO, Fa) + call AOtoMO(nBas_AOs, nBas_MOs, c, FbAO, Fb) F(1:nC, 1:nC ) = aC*Fa(1:nC, 1:nC ) + bC*Fb(1:nC, 1:nC ) F(1:nC, nC+1:nC+nO ) = Fb(1:nC, nC+1:nC+nO ) @@ -82,8 +85,10 @@ subroutine ROHF_fock_matrix(nBas,nOa,nOb,S,c,FaAO,FbAO,FAO) F(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) = Fa(nO+nC+1:nC+nO+nV, nC+1:nC+nO ) F(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) = aV*Fa(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) + bV*Fb(nO+nC+1:nC+nO+nV,nO+nC+1:nC+nO+nV) - call MOtoAO(nBas,S,c,F,FAO) - call MOtoAO(nBas,S,c,Fa,FaAO) - call MOtoAO(nBas,S,c,Fb,FbAO) + call MOtoAO(nBas_AOs, nBas_MOs, S, c, F, FAO) + call MOtoAO(nBas_AOs, nBas_MOs, S, c, Fa, FaAO) + call MOtoAO(nBas_AOs, nBas_MOs, S, c, Fb, FbAO) + + deallocate(F, Fa, Fb) end subroutine diff --git a/src/HF/UHF_search.f90 b/src/HF/UHF_search.f90 index 30ae0e2..dcc1025 100644 --- a/src/HF/UHF_search.f90 +++ b/src/HF/UHF_search.f90 @@ -124,8 +124,8 @@ subroutine UHF_search(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNu ! Transform dipole-related integrals do ixyz=1,ncart - call AOtoMO(nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,c(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90 index 9d29535..09939f6 100644 --- a/src/HF/print_ROHF.f90 +++ b/src/HF/print_ROHF.f90 @@ -1,14 +1,17 @@ -subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) + +! --- + +subroutine print_ROHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, Ex, EROHF, dipole) ! Print one- and two-electron energies and other stuff for RoHF calculation implicit none include 'parameters.h' - integer,intent(in) :: nBas + integer,intent(in) :: nBas_AOs, nBas_MOs integer,intent(in) :: nO(nspin) - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: c(nBas,nBas) + double precision,intent(in) :: eHF(nBas_MOs) + double precision,intent(in) :: c(nBas_AOs,nBas_MOs) double precision,intent(in) :: ENuc double precision,intent(in) :: ET(nspin) double precision,intent(in) :: EV(nspin) @@ -31,7 +34,7 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) do ispin=1,nspin if(nO(ispin) > 0) then HOMO(ispin) = eHF(nO(ispin)) - if(nO(ispin) < nBas) then + if(nO(ispin) < nBas_MOs) then LUMO(ispin) = eHF(nO(ispin)+1) else LUMO(ispin) = 0d0 @@ -102,13 +105,13 @@ subroutine print_ROHF(nBas,nO,eHF,c,ENuc,ET,EV,EJ,Ex,EROHF,dipole) write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') 'ROHF orbital coefficients ' write(*,'(A50)') '-----------------------------------------' - call matout(nBas,nBas,c) + call matout(nBas_AOs, nBas_MOs, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' ROHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas,eHF) + call vecout(nBas_MOs, eHF) write(*,*) end subroutine diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 1c9545b..980a102 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -112,7 +112,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d allocate(cHF(nBas_AOs,nBas_MOs)) allocate(eHF(nBas_MOs)) allocate(PHF(nBas_AOs,nBas_AOs)) - allocate(dipole_int_MO(nBas_AOs,nBas_AOs,ncart)) + allocate(dipole_int_MO(nBas_MOs,nBas_MOs,ncart)) allocate(ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)) !---------------------! @@ -122,7 +122,6 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doRHF) then call wall_time(start_HF) - ! TODO call RHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF) call wall_time(end_HF) @@ -137,8 +136,8 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_HF) ! TODO - call ROHF(dotest,maxSCF_HF,thresh_HF,max_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & - nBas_AOs,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF) + call ROHF(dotest, maxSCF_HF, thresh_HF, max_diis_HF, guess_type, mix, level_shift, nNuc, ZNuc, rNuc, ENuc, & + nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI_AO, dipole_int_AO, X, ERHF, eHF, cHF, PHF) call wall_time(end_HF) t_HF = end_HF - start_HF @@ -161,7 +160,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d do ixyz = 1, ncart ! TODO - call AOtoMO(nBas_AOs,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) + call AOtoMO(nBas_AOs,nBas_MOs,cHF,dipole_int_AO(:,:,ixyz),dipole_int_MO(:,:,ixyz)) end do ! 4-index transform diff --git a/src/QuAcK/UQuAcK.f90 b/src/QuAcK/UQuAcK.f90 index e3453b0..6e45b0d 100644 --- a/src/QuAcK/UQuAcK.f90 +++ b/src/QuAcK/UQuAcK.f90 @@ -163,8 +163,8 @@ subroutine UQuAcK(dotest,doUHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,do ! Read and transform dipole-related integrals do ixyz=1,ncart - call AOtoMO(nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) - call AOtoMO(nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) + call AOtoMO(nBas,nBas,cHF(:,:,1),dipole_int_AO(:,:,ixyz),dipole_int_aa(:,:,ixyz)) + call AOtoMO(nBas,nBas,cHF(:,:,2),dipole_int_AO(:,:,ixyz),dipole_int_bb(:,:,ixyz)) end do ! 4-index transform for (aa|aa) block