diff --git a/src/AOtoMO/AOtoMO.f90 b/src/AOtoMO/AOtoMO.f90 index 8c5ce8f..8383273 100644 --- a/src/AOtoMO/AOtoMO.f90 +++ b/src/AOtoMO/AOtoMO.f90 @@ -1,36 +1,30 @@ -subroutine AOtoMO(nBas_AOs, nBas_MOs, C, M_AOs, M_MOs) +subroutine AOtoMO(nBas, nOrb, C, M_AOs, M_MOs) -! Perform AO to MO transformation of a matrix M_AOs for given coefficients c -! M_MOs = C.T M_AOs 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, nOrb + double precision, intent(in) :: C(nBas,nOrb) + double precision, intent(in) :: M_AOs(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) + double precision, intent(out) :: M_MOs(nOrb,nOrb) -! Local variables + double precision, allocatable :: AC(:,:) - double precision,allocatable :: AC(:,:) - -! Output variables - - double precision,intent(out) :: M_MOs(nBas_MOs,nBas_MOs) - - allocate(AC(nBas_AOs,nBas_MOs)) + allocate(AC(nBas,nOrb)) !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("N", "N", nBas, nOrb, nBas, 1.d0, & + M_AOs(1,1), nBas, C(1,1), nBas, & + 0.d0, AC(1,1), nBas) - 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) + call dgemm("T", "N", nOrb, nOrb, nBas, 1.d0, & + C(1,1), nBas, AC(1,1), nBas, & + 0.d0, M_MOs(1,1), nOrb) deallocate(AC) diff --git a/src/AOtoMO/AOtoMO_ERI_RHF.f90 b/src/AOtoMO/AOtoMO_ERI_RHF.f90 index 6b1b95b..f9f64a2 100644 --- a/src/AOtoMO/AOtoMO_ERI_RHF.f90 +++ b/src/AOtoMO/AOtoMO_ERI_RHF.f90 @@ -1,7 +1,7 @@ ! --- -subroutine AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) +subroutine AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm @@ -10,9 +10,9 @@ subroutine AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs - double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) - double precision,intent(in) :: c(nBas_AOs,nBas_MOs) + integer,intent(in) :: nBas, nOrb + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: c(nBas,nOrb) ! Local variables @@ -21,35 +21,35 @@ subroutine AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) ! Output variables - double precision,intent(out) :: ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) + double precision,intent(out) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) ! Memory allocation - allocate(a2(nBas_AOs,nBas_AOs,nBas_AOs,nBas_MOs)) - allocate(a1(nBas_AOs,nBas_AOs,nBas_MOs,nBas_MOs)) + allocate(a2(nBas,nBas,nBas,nOrb)) + allocate(a1(nBas,nBas,nOrb,nOrb)) ! Four-index transform via semi-direct O(N^5) algorithm - call dgemm( 'T', 'N', nBas_AOs*nBas_AOs*nBas_AOs, nBas_MOs, nBas_AOs, 1.d0 & - , ERI_AO(1,1,1,1), nBas_AOs, c(1,1), nBas_AOs & - , 0.d0, a2(1,1,1,1), nBas_AOs*nBas_AOs*nBas_AOs) + call dgemm( 'T', 'N', nBas*nBas*nBas, nOrb, nBas, 1.d0 & + , ERI_AO(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, a2(1,1,1,1), nBas*nBas*nBas) - call dgemm( 'T', 'N', nBas_AOs*nBas_AOs*nBas_MOs, nBas_MOs, nBas_AOs, 1.d0 & - , a2(1,1,1,1), nBas_AOs, c(1,1), nBas_AOs & - , 0.d0, a1(1,1,1,1), nBas_AOs*nBas_AOs*nBas_MOs) + call dgemm( 'T', 'N', nBas*nBas*nOrb, nOrb, nBas, 1.d0 & + , a2(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, a1(1,1,1,1), nBas*nBas*nOrb) deallocate(a2) - allocate(a2(nBas_AOs,nBas_MOs,nBas_MOs,nBas_MOs)) + allocate(a2(nBas,nOrb,nOrb,nOrb)) - call dgemm( 'T', 'N', nBas_AOs*nBas_MOs*nBas_MOs, nBas_MOs, nBas_AOs, 1.d0 & - , a1(1,1,1,1), nBas_AOs, c(1,1), nBas_AOs & - , 0.d0, a2(1,1,1,1), nBas_AOs*nBas_MOs*nBas_MOs) + call dgemm( 'T', 'N', nBas*nOrb*nOrb, nOrb, nBas, 1.d0 & + , a1(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, a2(1,1,1,1), nBas*nOrb*nOrb) deallocate(a1) - call dgemm( 'T', 'N', nBas_MOs*nBas_MOs*nBas_MOs, nBas_MOs, nBas_AOs, 1.d0 & - , a2(1,1,1,1), nBas_AOs, c(1,1), nBas_AOs & - , 0.d0, ERI_MO(1,1,1,1), nBas_MOs*nBas_MOs*nBas_MOs) + call dgemm( 'T', 'N', nOrb*nOrb*nOrb, nOrb, nBas, 1.d0 & + , a2(1,1,1,1), nBas, c(1,1), nBas & + , 0.d0, ERI_MO(1,1,1,1), nOrb*nOrb*nOrb) deallocate(a2) diff --git a/src/AOtoMO/MOtoAO.f90 b/src/AOtoMO/MOtoAO.f90 index 06abb38..a5ffaed 100644 --- a/src/AOtoMO/MOtoAO.f90 +++ b/src/AOtoMO/MOtoAO.f90 @@ -1,47 +1,38 @@ -subroutine MOtoAO(nBas_AOs, nBas_MOs, S, C, M_MOs, M_AOs) +subroutine MOtoAO(nBas, nOrb, S, C, M_MOs, M_AOs) -! 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 -! + ! 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, nOrb + double precision, intent(in) :: S(nBas,nBas) + double precision, intent(in) :: C(nBas,nOrb) + double precision, intent(in) :: M_MOs(nOrb,nOrb) + double precision, intent(out) :: M_AOs(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) + double precision, allocatable :: SC(:,:),BSC(:,:) -! Local variables - double precision,allocatable :: SC(:,:),BSC(:,:) - -! Output variables - - double precision,intent(out) :: M_AOs(nBas_AOs,nBas_AOs) - -! Memory allocation - - allocate(SC(nBas_AOs,nBas_MOs), BSC(nBas_MOs,nBas_AOs)) + allocate(SC(nBas,nOrb), BSC(nOrb,nBas)) !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", "N", nBas, nOrb, nBas, 1.d0, & + S(1,1), nBas, C(1,1), nBas, & + 0.d0, SC(1,1), nBas) - 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", "T", nOrb, nBas, nOrb, 1.d0, & + M_MOs(1,1), nOrb, SC(1,1), nBas, & + 0.d0, BSC(1,1), nOrb) - 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) + call dgemm("N", "N", nBas, nBas, nOrb, 1.d0, & + SC(1,1), nBas, BSC(1,1), nOrb, & + 0.d0, M_AOs(1,1), nBas) deallocate(SC, BSC) diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90 index a95121b..57538ea 100644 --- a/src/CC/RCC.f90 +++ b/src/CC/RCC.f90 @@ -2,7 +2,7 @@ ! --- subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, & - maxSCF, thresh, max_diis, nBas_AOs, nBas_MOs, nC, nO, nV, nR, Hc, ERI, ENuc, ERHF, eHF, cHF) + maxSCF, thresh, max_diis, nBas, nOrb, nC, nO, nV, nR, Hc, ERI, ENuc, ERHF, eHF, cHF) ! Coupled-cluster module @@ -27,17 +27,17 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) - double precision,intent(in) :: ERI(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) ! Local variables @@ -50,7 +50,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(doCCD) then call wall_time(start_CC) - call CCD(dotest,maxSCF,thresh,max_diis,nBas_MOs,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call CCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -66,7 +66,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(doDCD) then call wall_time(start_CC) - call DCD(dotest,maxSCF,thresh,max_diis,nBas_MOs,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call DCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -84,7 +84,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(doCCSD) then call wall_time(start_CC) - call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas_MOs,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -100,7 +100,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(dodrCCD) then call wall_time(start_CC) - call drCCD(dotest,maxSCF,thresh,max_diis,nBas_MOs,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call drCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -116,7 +116,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(dorCCD) then call wall_time(start_CC) - call rCCD(dotest,maxSCF,thresh,max_diis,nBas_MOs,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call rCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -132,7 +132,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(docrCCD) then call wall_time(start_CC) - call crCCD(dotest,maxSCF,thresh,max_diis,nBas_MOs,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call crCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -148,7 +148,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(dolCCD) then call wall_time(start_CC) - call lCCD(dotest,maxSCF,thresh,max_diis,nBas_MOs,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call lCCD(dotest,maxSCF,thresh,max_diis,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -164,7 +164,7 @@ subroutine RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, d if(dopCCD) then call wall_time(start_CC) - call pCCD(dotest, maxSCF, thresh, max_diis, nBas_AOs, nBas_MOs, & + call pCCD(dotest, maxSCF, thresh, max_diis, nBas, nOrb, & nC, nO, nV, nR, Hc, ERI, ENuc, ERHF, eHF, cHF) call wall_time(end_CC) diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index 2f5afb2..7dfd3cd 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -1,7 +1,7 @@ ! --- -subroutine pCCD(dotest, maxSCF, thresh, max_diis, nBas_AOs, nBas_MOs, & +subroutine pCCD(dotest, maxSCF, thresh, max_diis, nBas, nOrb, & nC, nO, nV, nR, Hc, ERI, ENuc, ERHF, eHF, cHF) ! pair CCD module @@ -16,12 +16,12 @@ subroutine pCCD(dotest, maxSCF, thresh, max_diis, nBas_AOs, nBas_MOs, & integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas_AOs, nBas_MOs, nC, nO, nV, nR + integer,intent(in) :: nBas, nOrb, nC, nO, nV, nR double precision,intent(in) :: ENuc,ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) - double precision,intent(in) :: ERI(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) ! Local variables @@ -94,7 +94,7 @@ subroutine pCCD(dotest, maxSCF, thresh, max_diis, nBas_AOs, nBas_MOs, & allocate(eO(O),eV(V),delta_OV(O,V)) eO(:) = eHF(nC+1:nO) - eV(:) = eHF(nO+1:nBas_MOs-nR) + eV(:) = eHF(nO+1:nOrb-nR) call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV) diff --git a/src/CI/RCI.f90 b/src/CI/RCI.f90 index 1d068f3..3a9905e 100644 --- a/src/CI/RCI.f90 +++ b/src/CI/RCI.f90 @@ -1,7 +1,7 @@ ! --- -subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nBas_MOs, & +subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nOrb, & nC, nO, nV, nR, nS, ERI, dipole_int, epsHF, EHF) ! Configuration interaction module @@ -21,16 +21,16 @@ subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, n logical,intent(in) :: singlet logical,intent(in) :: triplet - integer,intent(in) :: nBas_MOs + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(nBas_MOs) - double precision,intent(in) :: ERI(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs) - double precision,intent(in) :: dipole_int(nBas_MOs,nBas_MOs,ncart) + double precision,intent(in) :: epsHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables @@ -43,7 +43,7 @@ subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, n if(doCIS) then call wall_time(start_CI) - call RCIS(dotest,singlet,triplet,doCIS_D,nBas_MOs,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF) + call RCIS(dotest,singlet,triplet,doCIS_D,nOrb,nC,nO,nV,nR,nS,ERI,dipole_int,epsHF) call wall_time(end_CI) t_CI = end_CI - start_CI @@ -59,7 +59,7 @@ subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, n if(doCID) then call wall_time(start_CI) - call CID(dotest,singlet,triplet,nBas_MOs,nC,nO,nV,nR,ERI,epsHF,EHF) + call CID(dotest,singlet,triplet,nOrb,nC,nO,nV,nR,ERI,epsHF,EHF) call wall_time(end_CI) t_CI = end_CI - start_CI @@ -75,7 +75,7 @@ subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, n if(doCISD) then call wall_time(start_CI) - call CISD(dotest,singlet,triplet,nBas_MOs,nC,nO,nV,nR,ERI,epsHF,EHF) + call CISD(dotest,singlet,triplet,nOrb,nC,nO,nV,nR,ERI,epsHF,EHF) call wall_time(end_CI) t_CI = end_CI - start_CI @@ -92,7 +92,7 @@ subroutine RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, n call wall_time(start_CI) write(*,*) ' FCI is not yet implemented! Sorry.' -! call FCI(nBas_MOs,nC,nO,nV,nR,ERI,epsHF) +! call FCI(nOrb,nC,nO,nV,nR,ERI,epsHF) call wall_time(end_CI) t_CI = end_CI - start_CI diff --git a/src/GF/RGF.f90 b/src/GF/RGF.f90 index 72911e4..0f1deaf 100644 --- a/src/GF/RGF.f90 +++ b/src/GF/RGF.f90 @@ -3,7 +3,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm, maxSCF, & thresh, max_diis, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, linearize, & - eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, EHF, & + eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, & S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF) ! Green's function module @@ -42,7 +42,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -50,18 +50,18 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren integer,intent(in) :: nS double precision,intent(in) :: EHF - double precision,intent(in) :: epsHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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(in) :: 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(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) + double precision,intent(in) :: epsHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -75,7 +75,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren call wall_time(start_GF) call RG0F2(dotest, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, & - linearize, eta, regularize, nBas_MOs, nC, nO, nV, nR, nS, & + linearize, eta, regularize, nOrb, nC, nO, nV, nR, nS, & ENuc, EHF, ERI_MO, dipole_int_MO, epsHF) call wall_time(end_GF) @@ -93,7 +93,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren call wall_time(start_GF) call evRGF2(dotest,dophBSE,doppBSE,TDA,dBSE,dTDA,maxSCF,thresh,max_diis, & - singlet,triplet,linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,EHF, & + singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,EHF, & ERI_MO,dipole_int_MO,epsHF) call wall_time(end_GF) @@ -112,7 +112,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren call wall_time(start_GF) call qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, & - rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, EHF, S, & + rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, EHF, S, & X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, epsHF) call wall_time(end_GF) @@ -129,7 +129,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren if(doufG0F02) then call wall_time(start_GF) - call ufRG0F02(dotest, nBas_MOs, nC, nO, nV, nR, nS, ENuc, EHF, ERI_MO, epsHF) + call ufRG0F02(dotest, nOrb, nC, nO, nV, nR, nS, ENuc, EHF, ERI_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF @@ -145,7 +145,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren if(doG0F3) then call wall_time(start_GF) - call RG0F3(dotest, renorm, nBas_MOs, nC, nO, nV, nR, ERI_MO, epsHF) + call RG0F3(dotest, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF @@ -161,7 +161,7 @@ subroutine RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, ren if(doevGF3) then call wall_time(start_GF) - call evRGF3(dotest, maxSCF, thresh, max_diis, renorm, nBas_MOs, nC, nO, nV, nR, ERI_MO, epsHF) + call evRGF3(dotest, maxSCF, thresh, max_diis, renorm, nOrb, nC, nO, nV, nR, ERI_MO, epsHF) call wall_time(end_GF) t_GF = end_GF - start_GF diff --git a/src/GF/print_qsRGF2.f90 b/src/GF/print_qsRGF2.f90 index 42132c0..feb590f 100644 --- a/src/GF/print_qsRGF2.f90 +++ b/src/GF/print_qsRGF2.f90 @@ -1,7 +1,7 @@ ! --- -subroutine print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, c, & +subroutine print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, c, & SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF, dipole) ! Print one-electron energies and other stuff for qsGF2 @@ -11,17 +11,17 @@ subroutine print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, c, ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: eGF(nBas_MOs) - double precision,intent(in) :: c(nBas_AOs,nBas_MOs) - double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs) - double precision,intent(in) :: Z(nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGF(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: ET double precision,intent(in) :: EV double precision,intent(in) :: EJ @@ -57,7 +57,7 @@ subroutine print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, c, '|','#','|','e_HF (eV)','|','Sig_c (eV)','|','Z','|','e_QP (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do q = 1, nBas_MOs + do q = 1, nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',q,'|',eHF(q)*HaToeV,'|',SigC(q,q)*HaToeV,'|',Z(q),'|',eGF(q)*HaToeV,'|' end do @@ -106,12 +106,12 @@ subroutine print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, c, write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGF2 MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas_AOs, nBas_MOs, c) + call matout(nBas, nOrb, c) write(*,*) write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGF2 MO energies' write(*,'(A50)') '---------------------------------------' - call matout(nBas_MOs, 1, eGF) + call matout(nOrb, 1, eGF) write(*,*) end if diff --git a/src/GF/qsRGF2.f90 b/src/GF/qsRGF2.f90 index 4287676..651c876 100644 --- a/src/GF/qsRGF2.f90 +++ b/src/GF/qsRGF2.f90 @@ -3,7 +3,7 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, & - rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, & + rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, & S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GF2 calculation @@ -33,25 +33,25 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs,nBas_MOs,nC,nO,nV,nR,nS + integer,intent(in) :: nBas,nOrb,nC,nO,nV,nR,nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBas_AOs_Sq + integer :: nBas_Sq integer :: ispin integer :: n_diis double precision :: EqsGF2 @@ -98,7 +98,7 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & ! Stuff - nBas_AOs_Sq = nBas_AOs*nBas_AOs + nBas_Sq = nBas*nBas ! TDA @@ -109,27 +109,27 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & ! Memory allocation - allocate(eGF(nBas_MOs)) - allocate(eOld(nBas_MOs)) + allocate(eGF(nOrb)) + allocate(eOld(nOrb)) - allocate(c(nBas_AOs,nBas_MOs)) + allocate(c(nBas,nOrb)) - allocate(cp(nBas_MOs,nBas_MOs)) - allocate(Fp(nBas_MOs,nBas_MOs)) + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) - allocate(P(nBas_AOs,nBas_AOs)) - allocate(F(nBas_AOs,nBas_AOs)) - allocate(J(nBas_AOs,nBas_AOs)) - allocate(K(nBas_AOs,nBas_AOs)) - allocate(error(nBas_AOs,nBas_AOs)) + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(error(nBas,nBas)) - allocate(Z(nBas_MOs)) - allocate(SigC(nBas_MOs,nBas_MOs)) + allocate(Z(nOrb)) + allocate(SigC(nOrb,nOrb)) - allocate(SigCp(nBas_AOs,nBas_AOs)) + allocate(SigCp(nBas,nBas)) - allocate(error_diis(nBas_AOs_Sq,max_diis)) - allocate(F_diis(nBas_AOs_Sq,max_diis)) + allocate(error_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Initialization @@ -157,25 +157,25 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & ! Buid Hartree matrix - call Hartree_matrix_AO_basis(nBas_AOs, P, ERI_AO, J) + call Hartree_matrix_AO_basis(nBas, P, ERI_AO, J) ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas_AOs, P, ERI_AO, K) + call exchange_matrix_AO_basis(nBas, P, ERI_AO, K) ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute self-energy and renormalization factor if(regularize) then - call GF2_reg_self_energy(eta, nBas_MOs, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) + call GF2_reg_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) else - call GF2_self_energy(eta, nBas_MOs, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) + call GF2_self_energy(eta, nOrb, nC, nO, nV, nR, eGF, ERI_MO, SigC, Z) end if @@ -183,7 +183,7 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO(nBas_AOs, nBas_MOs, S, c, SigC, SigCp) + call MOtoAO(nBas, nOrb, S, c, SigC, SigCp) ! Solve the quasi-particle equation @@ -197,7 +197,7 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & n_diis = min(n_diis+1, max_diis) if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F) else n_diis = 0 end if @@ -206,7 +206,7 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas_MOs, cp, eGF) + call diagonalize_matrix(nOrb, cp, eGF) c = matmul(X, cp) ! Compute new density matrix in the AO basis @@ -224,23 +224,23 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & ! Kinetic energy - ET = trace_matrix(nBas_AOs, matmul(P, T)) + ET = trace_matrix(nBas, matmul(P, T)) ! Potential energy - EV = trace_matrix(nBas_AOs, matmul(P, V)) + EV = trace_matrix(nBas, matmul(P, V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas_AOs, matmul(P, J)) + EJ = 0.5d0*trace_matrix(nBas, matmul(P, J)) ! Exchange energy - Ex = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K)) + Ex = 0.25d0*trace_matrix(nBas, matmul(P, K)) ! Correlation energy - call RMP2(.false., regularize, nBas_MOs, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec) + call RMP2(.false., regularize, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, EqsGF2, eGF, Ec) ! Total energy @@ -251,8 +251,8 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & ! Print results !------------------------------------------------------------------------ - call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) - call print_qsRGF2(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGF, & + call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) + call print_qsRGF2(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGF, & c, SigC, Z, ENuc, ET, EV, EJ, Ex, Ec, EqsGF2, dipole) end do @@ -283,7 +283,7 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & if(dophBSE) then - call GF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, nC, nO, & + call GF2_phBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, nC, nO, & nV, nR, nS, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) @@ -302,7 +302,7 @@ subroutine qsRGF2(dotest, maxSCF, thresh, max_diis, dophBSE, doppBSE, TDA, & if(doppBSE) then - call GF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, & + call GF2_ppBSE2(TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & nC, nO, nV, nR, ERI_MO, dipole_int_MO, eGF, EcBSE) write(*,*) diff --git a/src/GT/RGT.f90 b/src/GT/RGT.f90 index a4f46b6..74d1d12 100644 --- a/src/GT/RGT.f90 +++ b/src/GT/RGT.f90 @@ -4,7 +4,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & doppBSE, TDA_T, TDA, dBSE, dTDA, singlet, triplet, linearize, eta, regularize, & - nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, & + nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, & V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! T-matrix module @@ -48,7 +48,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -56,18 +56,18 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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(in) :: 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(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -82,7 +82,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG call wall_time(start_GT) call RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -99,7 +99,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG call wall_time(start_GT) call evRGTpp(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,singlet,triplet, & - linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -116,7 +116,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG call wall_time(start_GT) call qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, dBSE, & - dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, & + dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, & nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) @@ -133,7 +133,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG if(doufG0T0pp) then call wall_time(start_GT) - call ufG0T0pp(dotest,TDA_T,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufG0T0pp(dotest,TDA_T,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -150,7 +150,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG call wall_time(start_GT) call RG0T0eh(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -167,7 +167,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG call wall_time(start_GT) call evRGTeh(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_T,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + singlet,triplet,linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GT) t_GT = end_GT - start_GT @@ -185,7 +185,7 @@ subroutine RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevG call wall_time(start_GT) call qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, & dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, & - ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & + ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) diff --git a/src/GT/print_qsRGTeh.f90 b/src/GT/print_qsRGTeh.f90 index 72ffe99..e1e4f95 100644 --- a/src/GT/print_qsRGTeh.f90 +++ b/src/GT/print_qsRGTeh.f90 @@ -1,7 +1,7 @@ ! --- -subroutine print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, & +subroutine print_qsRGTeh(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, & Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) ! Print one-electron energies and other stuff for qsGTeh @@ -11,7 +11,7 @@ subroutine print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -23,11 +23,11 @@ subroutine print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: eGT(nBas_MOs) - double precision,intent(in) :: c(nBas_AOs,nBas_MOs) - double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs) - double precision,intent(in) :: Z(nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGT(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: EqsGT double precision,intent(in) :: dipole(ncart) @@ -62,7 +62,7 @@ subroutine print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c '|','#','|','e_HF (eV)','|','Sig_GTeh (eV)','|','Z','|','e_GTeh (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas_MOs + do p=1,nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' end do @@ -113,13 +113,13 @@ subroutine print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTeh MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas_AOs, nBas_MOs, c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTeh MO energies' write(*,'(A50)') '---------------------------------------' - call vecout(nBas_MOs, eGT) + call vecout(nOrb, eGT) write(*,*) end if diff --git a/src/GT/print_qsRGTpp.f90 b/src/GT/print_qsRGTpp.f90 index 9d1479d..8eab4c5 100644 --- a/src/GT/print_qsRGTpp.f90 +++ b/src/GT/print_qsRGTpp.f90 @@ -1,7 +1,7 @@ ! --- -subroutine print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, Z, & +subroutine print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, SigC, Z, & ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) ! Print one-electron energies and other stuff for qsGT @@ -11,7 +11,7 @@ subroutine print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -23,11 +23,11 @@ subroutine print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: eGT(nBas_MOs) - double precision,intent(in) :: c(nBas_AOs,nBas_MOs) - double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs) - double precision,intent(in) :: Z(nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGT(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: EqsGT double precision,intent(in) :: dipole(ncart) @@ -62,7 +62,7 @@ subroutine print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c '|','#','|','e_HF (eV)','|','Sig_GTpp (eV)','|','Z','|','e_GTpp (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas_MOs + do p=1,nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGT(p)*HaToeV,'|' end do @@ -113,13 +113,13 @@ subroutine print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTpp MO coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas_AOs, nBas_MOs, c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A32)') ' qsGTpp MO energies' write(*,'(A50)') '---------------------------------------' - call vecout(nBas_MOs, eGT) + call vecout(nOrb, eGT) write(*,*) end if diff --git a/src/GT/qsRGTeh.f90 b/src/GT/qsRGTeh.f90 index e608d9a..5f6acab 100644 --- a/src/GT/qsRGTeh.f90 +++ b/src/GT/qsRGTeh.f90 @@ -3,7 +3,7 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, & dophBSE2, TDA_T, TDA, dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, & - ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & + ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, & Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GTeh calculation @@ -37,31 +37,31 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables logical :: dRPA = .false. integer :: nSCF - integer :: nBas_AOs_Sq + integer :: nBas_Sq integer :: ispin integer :: n_diis double precision :: ET @@ -117,7 +117,7 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Stuff - nBas_AOs_Sq = nBas_AOs*nBas_AOs + nBas_Sq = nBas*nBas ! TDA for T @@ -137,26 +137,26 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d allocate(Aph(nS,nS), Bph(nS,nS), Om(nS), XpY(nS,nS), XmY(nS,nS)) - allocate(eGT(nBas_MOs)) - allocate(eOld(nBas_MOs)) - allocate(Z(nBas_MOs)) + allocate(eGT(nOrb)) + allocate(eOld(nOrb)) + allocate(Z(nOrb)) - allocate(c(nBas_AOs,nBas_MOs)) + allocate(c(nBas,nOrb)) - allocate(cp(nBas_MOs,nBas_MOs)) - allocate(Fp(nBas_MOs,nBas_MOs)) - allocate(Sig(nBas_MOs,nBas_MOs)) + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + allocate(Sig(nOrb,nOrb)) - allocate(P(nBas_AOs,nBas_AOs)) - allocate(F(nBas_AOs,nBas_AOs)) - allocate(J(nBas_AOs,nBas_AOs)) - allocate(K(nBas_AOs,nBas_AOs)) - allocate(Sigp(nBas_AOs,nBas_AOs)) - allocate(err(nBas_AOs,nBas_AOs)) + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(Sigp(nBas,nBas)) + allocate(err(nBas,nBas)) - allocate(err_diis(nBas_AOs_Sq,max_diis), F_diis(nBas_AOs_Sq,max_diis)) + allocate(err_diis(nBas_Sq,max_diis), F_diis(nBas_Sq,max_diis)) - allocate(rhoL(nBas_MOs,nBas_MOs,nS), rhoR(nBas_MOs,nBas_MOs,nS)) + allocate(rhoL(nOrb,nOrb,nS), rhoR(nOrb,nOrb,nS)) ! Initialization @@ -185,20 +185,20 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Buid Hartree matrix - call Hartree_matrix_AO_basis(nBas_AOs,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas_AOs,P,ERI_AO,K) + call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute linear response - call phLR_A(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) - if(.not.TDA_T) call phLR_B(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGT,ERI_MO,Aph) + if(.not.TDA_T) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_T,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -206,17 +206,17 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Compute correlation part of the self-energy - call GTeh_excitation_density(nBas_MOs,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) + call GTeh_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,XmY,rhoL,rhoR) - if(regularize) call GTeh_regularization(nBas_MOs,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) + if(regularize) call GTeh_regularization(nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR) - call GTeh_self_energy(eta,nBas_MOs,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z) + call GTeh_self_energy(eta,nOrb,nC,nO,nV,nR,nS,eGT,Om,rhoL,rhoR,EcGM,Sig,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas_AOs, nBas_MOs, S, c, Sig, Sigp) + call MOtoAO(nBas, nOrb, S, c, Sig, Sigp) ! Solve the quasi-particle equation @@ -231,7 +231,7 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,err_diis,F_diis,err,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) end if @@ -239,7 +239,7 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas_MOs, cp, eGT) + call diagonalize_matrix(nOrb, cp, eGT) c = matmul(X,cp) ! Compute new density matrix in the AO basis @@ -257,19 +257,19 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Kinetic energy - ET = trace_matrix(nBas_AOs,matmul(P,T)) + ET = trace_matrix(nBas,matmul(P,T)) ! Potential energy - EV = trace_matrix(nBas_AOs,matmul(P,V)) + EV = trace_matrix(nBas,matmul(P,V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas_AOs,matmul(P,J)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) ! Exchange energy - Ex = 0.25d0*trace_matrix(nBas_AOs,matmul(P,K)) + Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) ! Total energy @@ -277,8 +277,8 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Print results - call dipole_moment(nBas_AOs,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGTeh(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGT, c, Sig, & + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsRGTeh(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGT, c, Sig, & Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGT, dipole) end do @@ -310,7 +310,7 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! if(BSE) then -! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & +! call Bethe_Salpeter(BSE2,TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO, & ! eGT,eGT,EcBSE) ! if(exchange_kernel) then @@ -345,7 +345,7 @@ subroutine qsRGTeh(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! end if -! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas_AOs,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) +! call ACFDT(exchange_kernel,doXBS,.true.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) ! write(*,*) ! write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GT/qsRGTpp.f90 b/src/GT/qsRGTpp.f90 index 1f8e17e..15a29dc 100644 --- a/src/GT/qsRGTpp.f90 +++ b/src/GT/qsRGTpp.f90 @@ -2,7 +2,7 @@ ! --- subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, TDA_T, TDA, & - dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, & + dBSE, dTDA, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, & nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GT calculation @@ -34,26 +34,26 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC,nO,nV,nR,nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBas_AOs_Sq + integer :: nBas_Sq integer :: ispin integer :: iblock integer :: n_diis @@ -123,7 +123,7 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Stuff - nBas_AOs_Sq = nBas_AOs*nBas_AOs + nBas_Sq = nBas*nBas ! TDA for T @@ -141,30 +141,30 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Memory allocation - allocate(eGT(nBas_MOs)) - allocate(eOld(nBas_MOs)) - allocate(Z(nBas_MOs)) + allocate(eGT(nOrb)) + allocate(eOld(nOrb)) + allocate(Z(nOrb)) - allocate(c(nBas_AOs,nBas_MOs)) + allocate(c(nBas,nOrb)) - allocate(Fp(nBas_MOs,nBas_MOs)) - allocate(cp(nBas_MOs,nBas_MOs)) - allocate(Sig(nBas_MOs,nBas_MOs)) + allocate(Fp(nOrb,nOrb)) + allocate(cp(nOrb,nOrb)) + allocate(Sig(nOrb,nOrb)) - allocate(P(nBas_AOs,nBas_AOs)) - allocate(F(nBas_AOs,nBas_AOs)) - allocate(J(nBas_AOs,nBas_AOs)) - allocate(K(nBas_AOs,nBas_AOs)) - allocate(error(nBas_AOs,nBas_AOs)) - allocate(Sigp(nBas_AOs,nBas_AOs)) + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(error(nBas,nBas)) + allocate(Sigp(nBas,nBas)) - allocate(error_diis(nBas_AOs_Sq,max_diis)) - allocate(F_diis(nBas_AOs_Sq,max_diis)) + allocate(error_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) - allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nBas_MOs,nBas_MOs,nVVs)) - allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nBas_MOs,nBas_MOs,nOOs)) - allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nBas_MOs,nBas_MOs,nVVt)) - allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nBas_MOs,nBas_MOs,nOOt)) + allocate(Om1s(nVVs), X1s(nVVs,nVVs), Y1s(nOOs,nVVs), rho1s(nOrb,nOrb,nVVs)) + allocate(Om2s(nOOs), X2s(nVVs,nOOs), Y2s(nOOs,nOOs), rho2s(nOrb,nOrb,nOOs)) + allocate(Om1t(nVVt), X1t(nVVt,nVVt), Y1t(nOOt,nVVt), rho1t(nOrb,nOrb,nVVt)) + allocate(Om2t(nOOt), X2t(nVVt,nOOt), Y2t(nOOt,nOOt), rho2t(nOrb,nOrb,nOOt)) ! Initialization @@ -192,15 +192,15 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Buid Hartree matrix - call Hartree_matrix_AO_basis(nBas_AOs,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas_AOs,P,ERI_AO,K) + call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) ! AO to MO transformation of two-electron integrals - call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute linear response @@ -209,9 +209,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - call ppLR_C(iblock,nBas_MOs,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI_MO,Bpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -222,9 +222,9 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR_C(iblock,nBas_MOs,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) - call ppLR_D(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI_MO,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI_MO,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI_MO,Bpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) @@ -236,24 +236,24 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Compute correlation part of the self-energy iblock = 3 - call GTpp_excitation_density(iblock,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + call GTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) iblock = 4 - call GTpp_excitation_density(iblock,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + call GTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) if(regularize) then - call GTpp_regularization(eta,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) - call GTpp_regularization(eta,nBas_MOs,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t) + call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,eGT,Om1s,rho1s,Om2s,rho2s) + call GTpp_regularization(eta,nOrb,nC,nO,nV,nR,nOOt,nVVt,eGT,Om1t,rho1t,Om2t,rho2t) end if - call GTpp_self_energy(eta,nBas_MOs,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & + call GTpp_self_energy(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eGT,Om1s,rho1s,Om2s,rho2s, & Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) ! Make correlation self-energy Hermitian and transform it back to AO basis Sig = 0.5d0*(Sig + transpose(Sig)) - call MOtoAO(nBas_AOs, nBas_MOs, S, c, Sig, Sigp) + call MOtoAO(nBas, nOrb, S, c, Sig, Sigp) ! Solve the quasi-particle equation @@ -267,7 +267,7 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d n_diis = min(n_diis+1,max_diis) if(abs(rcond) > 1d-7) then - call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F) else n_diis = 0 end if @@ -276,7 +276,7 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas_MOs, cp, eGT) + call diagonalize_matrix(nOrb, cp, eGT) c = matmul(X, cp) ! Compute new density matrix in the AO basis @@ -294,19 +294,19 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Kinetic energy - ET = trace_matrix(nBas_AOs,matmul(P,T)) + ET = trace_matrix(nBas,matmul(P,T)) ! Potential energy - EV = trace_matrix(nBas_AOs,matmul(P,V)) + EV = trace_matrix(nBas,matmul(P,V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas_AOs,matmul(P,J)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) ! Exchange energy - Ex = 0.25d0*trace_matrix(nBas_AOs,matmul(P,K)) + Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) ! Total energy @@ -314,8 +314,8 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d ! Print results - call dipole_moment(nBas_AOs,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGTpp(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, & + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsRGTpp(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, & eGT, c, Sig, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, & EqsGT, dipole) @@ -351,7 +351,7 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d if(dophBSE) then - call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas_MOs,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + call GTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI_MO,dipole_int_MO,eGT,eGT,EcBSE) @@ -387,7 +387,7 @@ subroutine qsRGTpp(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, d end if - call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas_MOs,nC,nO,nV,nR,nS, & + call GTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & Om2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcBSE) diff --git a/src/GW/RGW.f90 b/src/GW/RGW.f90 index 8498485..48408f7 100644 --- a/src/GW/RGW.f90 +++ b/src/GW/RGW.f90 @@ -3,7 +3,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF, thresh, max_diis, doACFDT, & exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, & - linearize, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, & + linearize, eta, regularize, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, & S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Restricted GW module @@ -46,7 +46,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -54,18 +54,18 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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(in) :: 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(in) :: dipole_int_MO(nBas_MOs,nBas_MOs,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -79,7 +79,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS call wall_time(start_GW) call RG0W0(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas_MOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -96,7 +96,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS call wall_time(start_GW) call evRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,doppBSE, & - singlet,triplet,linearize,eta,regularize,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,dipole_int_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -114,7 +114,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS call wall_time(start_GW) call qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, & - ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, & + ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, & dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GW) @@ -133,7 +133,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS call wall_time(start_GW) call SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, & dophBSE, dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, & - nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, & + nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, & ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, & PHF, cHF, eHF) call wall_time(end_GW) @@ -152,7 +152,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS call wall_time(start_GW) ! TODO - call ufG0W0(dotest,TDA_W,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW @@ -169,7 +169,7 @@ subroutine RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxS call wall_time(start_GW) ! TODO - call ufGW(dotest,TDA_W,nBas_AOs,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) + call ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF) call wall_time(end_GW) t_GW = end_GW - start_GW diff --git a/src/GW/SRG_qsGW.f90 b/src/GW/SRG_qsGW.f90 index 32c2a19..807991b 100644 --- a/src/GW/SRG_qsGW.f90 +++ b/src/GW/SRG_qsGW.f90 @@ -3,7 +3,7 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, & BSE, BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nNuc, & - ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, & + ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, & X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -36,30 +36,30 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBas_AOs_Sq + integer :: nBas_Sq integer :: ispin integer :: ixyz integer :: n_diis @@ -118,7 +118,7 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, ! Stuff - nBas_AOs_Sq = nBas_AOs*nBas_AOs + nBas_Sq = nBas*nBas ! TDA for W @@ -136,32 +136,32 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, ! Memory allocation - allocate(eGW(nBas_MOs)) - allocate(eOld(nBas_MOs)) - allocate(Z(nBas_MOs)) + allocate(eGW(nOrb)) + allocate(eOld(nOrb)) + allocate(Z(nOrb)) - allocate(c(nBas_AOs,nBas_MOs)) + allocate(c(nBas,nOrb)) - allocate(cp(nBas_MOs,nBas_MOs)) - allocate(Fp(nBas_MOs,nBas_MOs)) - allocate(SigC(nBas_MOs,nBas_MOs)) + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + allocate(SigC(nOrb,nOrb)) - allocate(P(nBas_AOs,nBas_AOs)) - allocate(F(nBas_AOs,nBas_AOs)) - allocate(J(nBas_AOs,nBas_AOs)) - allocate(K(nBas_AOs,nBas_AOs)) - allocate(error(nBas_AOs,nBas_AOs)) - allocate(SigCp(nBas_AOs,nBas_AOs)) + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(error(nBas,nBas)) + allocate(SigCp(nBas,nBas)) allocate(Aph(nS,nS)) allocate(Bph(nS,nS)) allocate(Om(nS)) allocate(XpY(nS,nS)) allocate(XmY(nS,nS)) - allocate(rho(nBas_MOs,nBas_MOs,nS)) + allocate(rho(nOrb,nOrb,nS)) - allocate(error_diis(nBas_AOs_Sq,max_diis)) - allocate(F_diis(nBas_AOs_Sq,max_diis)) + allocate(error_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Initialization @@ -189,11 +189,11 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, ! Buid Hartree matrix call wall_time(t1) - call Hartree_matrix_AO_basis(nBas_AOs,P,ERI_AO,J) + call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas_AOs,P,ERI_AO,K) + call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) call wall_time(t2) tt=tt+t2-t1 @@ -202,10 +202,10 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, call wall_time(tao1) do ixyz = 1, ncart - call AOtoMO(nBas_AOs, nBas_MOs, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) + call AOtoMO(nBas, nOrb, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) call wall_time(tao2) @@ -215,8 +215,8 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, call wall_time(tlr1) - call phLR_A(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) - if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas_MOs,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + call phLR_A(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) @@ -230,13 +230,13 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, call wall_time(tex1) - call GW_excitation_density(nBas_MOs,nC,nO,nR,nS,ERI_MO,XpY,rho) + call GW_excitation_density(nOrb,nC,nO,nR,nS,ERI_MO,XpY,rho) call wall_time(tex2) tex=tex+tex2-tex1 call wall_time(tsrg1) - call SRG_self_energy(flow,nBas_MOs,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) + call SRG_self_energy(flow,nOrb,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) call wall_time(tsrg2) @@ -245,7 +245,7 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, ! Make correlation self-energy Hermitian and transform it back to AO basis call wall_time(tmo1) - call MOtoAO(nBas_AOs, nBas_MOs, S, c, SigC, SigCp) + call MOtoAO(nBas, nOrb, S, c, SigC, SigCp) call wall_time(tmo2) tmo = tmo + tmo2 - tmo1 ! Solve the quasi-particle equation @@ -261,7 +261,7 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,error_diis,F_diis,error,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,error_diis,F_diis,error,F) end if @@ -269,10 +269,10 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas_MOs, cp, eGW) + call diagonalize_matrix(nOrb, cp, eGW) c = matmul(X, cp) - call AOtoMO(nBas_AOs, nBas_MOs, c, SigCp, SigC) + call AOtoMO(nBas, nOrb, c, SigCp, SigC) ! Compute new density matrix in the AO basis @@ -289,19 +289,19 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, ! Kinetic energy - ET = trace_matrix(nBas_AOs,matmul(P,T)) + ET = trace_matrix(nBas,matmul(P,T)) ! Potential energy - EV = trace_matrix(nBas_AOs,matmul(P,V)) + EV = trace_matrix(nBas,matmul(P,V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas_AOs,matmul(P,J)) + EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) ! Exchange energy - Ex = 0.25d0*trace_matrix(nBas_AOs,matmul(P,K)) + Ex = 0.25d0*trace_matrix(nBas,matmul(P,K)) ! Total energy @@ -309,8 +309,8 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, ! Print results - call dipole_moment(nBas_AOs,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, & + call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) + call print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, & SigC, Z, ENuc, ET, EV, EJ, Ex, EcGM, EcRPA, EqsGW, dipole) end do @@ -343,7 +343,7 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, ! Cumulant expansion - call RGWC(dotest,eta,nBas_MOs,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) + call RGWC(dotest,eta,nOrb,nC,nO,nV,nR,nS,Om,rho,eHF,eGW,eGW,Z) ! Deallocate memory @@ -353,7 +353,7 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, if(BSE) then - call GW_phBSE(BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, & + call GW_phBSE(BSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) if(exchange_kernel) then @@ -389,7 +389,7 @@ subroutine SRG_qsGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, end if call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, BSE, singlet, triplet, & - eta, nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) + eta, nOrb, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GW/print_qsRGW.f90 b/src/GW/print_qsRGW.f90 index 7e90ce7..9fd695c 100644 --- a/src/GW/print_qsRGW.f90 +++ b/src/GW/print_qsRGW.f90 @@ -1,7 +1,7 @@ ! --- -subroutine print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, & +subroutine print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, & Z, ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole) ! Print useful information about qsRGW calculation @@ -11,7 +11,7 @@ subroutine print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nSCF double precision,intent(in) :: ENuc @@ -23,11 +23,11 @@ subroutine print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, double precision,intent(in) :: EcRPA double precision,intent(in) :: Conv double precision,intent(in) :: thresh - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: eGW(nBas_MOs) - double precision,intent(in) :: c(nBas_AOs,nBas_MOs) - double precision,intent(in) :: SigC(nBas_MOs,nBas_MOs) - double precision,intent(in) :: Z(nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: eGW(nOrb) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(in) :: SigC(nOrb,nOrb) + double precision,intent(in) :: Z(nOrb) double precision,intent(in) :: EqsGW double precision,intent(in) :: dipole(ncart) @@ -63,7 +63,7 @@ subroutine print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, '|','#','|','e_HF (eV)','|','Sig_GW (eV)','|','Z','|','e_GW (eV)','|' write(*,*)'-------------------------------------------------------------------------------' - do p=1,nBas_MOs + do p=1,nOrb write(*,'(1X,A1,1X,I3,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X,F15.6,1X,A1,1X)') & '|',p,'|',eHF(p)*HaToeV,'|',SigC(p,p)*HaToeV,'|',Z(p),'|',eGW(p)*HaToeV,'|' end do @@ -114,13 +114,13 @@ subroutine print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' Restricted qsGW orbital coefficients' write(*,'(A50)') '---------------------------------------' - call matout(nBas_AOs, nBas_MOs, c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' Restricted qsGW orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas_MOs, eGW) + call vecout(nOrb, eGW) write(*,*) end if diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index 395c30e..1756a0b 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -3,7 +3,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, & TDA_W, TDA, dBSE, dTDA, doppBSE, singlet, triplet, eta, regularize, nNuc, ZNuc, rNuc, & - ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, & + ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, & ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) ! Perform a quasiparticle self-consistent GW calculation @@ -38,30 +38,30 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nS double precision,intent(in) :: ERHF - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) - double precision,intent(in) :: PHF(nBas_AOs,nBas_AOs) - 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_AOs) - 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) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) + double precision,intent(in) :: PHF(nBas,nBas) + 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(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables integer :: nSCF - integer :: nBas_AOs_Sq + integer :: nBas_Sq integer :: ispin integer :: ixyz integer :: n_diis @@ -116,7 +116,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX ! Stuff - nBas_AOs_Sq = nBas_AOs*nBas_AOs + nBas_Sq = nBas*nBas ! TDA for W @@ -134,31 +134,31 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX ! Memory allocation - allocate(eGW(nBas_MOs)) - allocate(Z(nBas_MOs)) + allocate(eGW(nOrb)) + allocate(Z(nOrb)) - allocate(c(nBas_AOs,nBas_MOs)) + allocate(c(nBas,nOrb)) - allocate(cp(nBas_MOs,nBas_MOs)) - allocate(Fp(nBas_MOs,nBas_MOs)) - allocate(SigC(nBas_MOs,nBas_MOs)) + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) + allocate(SigC(nOrb,nOrb)) - allocate(P(nBas_AOs,nBas_AOs)) - allocate(F(nBas_AOs,nBas_AOs)) - allocate(J(nBas_AOs,nBas_AOs)) - allocate(K(nBas_AOs,nBas_AOs)) - allocate(err(nBas_AOs,nBas_AOs)) - allocate(SigCp(nBas_AOs,nBas_AOs)) + allocate(P(nBas,nBas)) + allocate(F(nBas,nBas)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) + allocate(err(nBas,nBas)) + allocate(SigCp(nBas,nBas)) allocate(Aph(nS,nS)) allocate(Bph(nS,nS)) allocate(Om(nS)) allocate(XpY(nS,nS)) allocate(XmY(nS,nS)) - allocate(rho(nBas_MOs,nBas_MOs,nS)) + allocate(rho(nOrb,nOrb,nS)) - allocate(err_diis(nBas_AOs_Sq,max_diis)) - allocate(F_diis(nBas_AOs_Sq,max_diis)) + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Initialization @@ -185,38 +185,38 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX ! Build Hartree-exchange matrix - call Hartree_matrix_AO_basis(nBas_AOs, P, ERI_AO, J) - call exchange_matrix_AO_basis(nBas_AOs, P, ERI_AO, K) + call Hartree_matrix_AO_basis(nBas, P, ERI_AO, J) + call exchange_matrix_AO_basis(nBas, P, ERI_AO, K) ! AO to MO transformation of two-electron integrals do ixyz = 1, ncart - call AOtoMO(nBas_AOs, nBas_MOs, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) + call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) ! Compute linear response - call phLR_A(ispin, dRPA, nBas_MOs, nC, nO, nV, nR, nS, 1d0, eGW, ERI_MO, Aph) - if(.not.TDA_W) call phLR_B(ispin, dRPA, nBas_MOs, nC, nO, nV, nR, nS, 1d0, ERI_MO, Bph) + call phLR_A(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, eGW, ERI_MO, Aph) + if(.not.TDA_W) call phLR_B(ispin, dRPA, nOrb, nC, nO, nV, nR, nS, 1d0, ERI_MO, Bph) call phLR(TDA_W, nS, Aph, Bph, EcRPA, Om, XpY, XmY) if(print_W) call print_excitation_energies('phRPA@GW@RHF','singlet',nS,Om) ! Compute correlation part of the self-energy - call GW_excitation_density(nBas_MOs, nC, nO, nR, nS, ERI_MO, XpY, rho) + call GW_excitation_density(nOrb, nC, nO, nR, nS, ERI_MO, XpY, rho) - if(regularize) call GW_regularization(nBas_MOs, nC, nO, nV, nR, nS, eGW, Om, rho) + if(regularize) call GW_regularization(nOrb, nC, nO, nV, nR, nS, eGW, Om, rho) - call GW_self_energy(eta, nBas_MOs, nC, nO, nV, nR, nS, eGW, Om, rho, EcGM, SigC, Z) + call GW_self_energy(eta, nOrb, nC, nO, nV, nR, nS, eGW, Om, rho, EcGM, SigC, Z) ! Make correlation self-energy Hermitian and transform it back to AO basis SigC = 0.5d0*(SigC + transpose(SigC)) - call MOtoAO(nBas_AOs, nBas_MOs, S, c, SigC, SigCp) + call MOtoAO(nBas, nOrb, S, c, SigC, SigCp) ! Solve the quasi-particle equation @@ -230,19 +230,19 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX ! Kinetic energy - ET = trace_matrix(nBas_AOs, matmul(P, T)) + ET = trace_matrix(nBas, matmul(P, T)) ! Potential energy - EV = trace_matrix(nBas_AOs, matmul(P, V)) + EV = trace_matrix(nBas, matmul(P, V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas_AOs, matmul(P, J)) + EJ = 0.5d0*trace_matrix(nBas, matmul(P, J)) ! Exchange energy - EK = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K)) + EK = 0.25d0*trace_matrix(nBas, matmul(P, K)) ! Total energy @@ -253,7 +253,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,err_diis,F_diis,err,F) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,F) end if @@ -261,9 +261,9 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX Fp = matmul(transpose(X), matmul(F, X)) cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas_MOs, cp, eGW) + call diagonalize_matrix(nOrb, cp, eGW) c = matmul(X, cp) - call AOtoMO(nBas_AOs, nBas_MOs, c, SigCp, SigC) + call AOtoMO(nBas, nOrb, c, SigCp, SigC) ! Density matrix @@ -271,8 +271,8 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX ! Print results - call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) - call print_qsRGW(nBas_AOs, nBas_MOs, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, Z, & + call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int_AO, dipole) + call print_qsRGW(nBas, nOrb, nO, nSCF, Conv, thresh, eHF, eGW, c, SigC, Z, & ENuc, ET, EV, EJ, EK, EcGM, EcRPA, EqsGW, dipole) end do @@ -304,7 +304,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX if(dophBSE) then call GW_phBSE(dophBSE2, TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, & - nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) + nOrb, nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eGW, eGW, EcBSE) if(exchange_kernel) then @@ -339,7 +339,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX end if call GW_phACFDT(exchange_kernel, doXBS, .true., TDA_W, TDA, dophBSE, singlet, triplet, & - eta, nBas_MOs, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) + eta, nOrb, nC, nO, nV, nR, nS, ERI_MO, eGW, eGW, EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -356,7 +356,7 @@ subroutine qsRGW(dotest, maxSCF, thresh, max_diis, doACFDT, exchange_kernel, doX if(doppBSE) then - call GW_ppBSE(TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nBas_MOs, & + call GW_ppBSE(TDA_W, TDA, dBSE, dTDA, singlet, triplet, eta, nOrb, & nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, eGW, EcBSE) EcBSE(2) = 3d0*EcBSE(2) diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index 5756744..1ae1b27 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -2,7 +2,7 @@ ! --- subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ZNuc, rNuc, ENuc, & - nBas_AOs, nBas_MOs, nO, S, T, V, Hc, ERI, dipole_int, X, ERHF, eHF, c, P) + nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, ERHF, eHF, c, P) ! Perform restricted Hartree-Fock calculation @@ -19,24 +19,24 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, double precision,intent(in) :: thresh double precision,intent(in) :: level_shift - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - 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) + 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,nOrb) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables integer :: nSCF - integer :: nBas_AOs_Sq + integer :: nBas_Sq integer :: n_diis double precision :: ET double precision :: EV @@ -59,9 +59,9 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Output variables double precision,intent(out) :: ERHF - double precision,intent(out) :: eHF(nBas_MOs) - double precision,intent(inout):: c(nBas_AOs,nBas_MOs) - double precision,intent(out) :: P(nBas_AOs,nBas_AOs) + double precision,intent(out) :: eHF(nOrb) + double precision,intent(inout):: c(nBas,nOrb) + double precision,intent(out) :: P(nBas,nBas) ! Hello world @@ -73,35 +73,37 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Useful quantities - nBas_AOs_Sq = nBas_AOs*nBas_AOs + nBas_Sq = nBas*nBas ! Memory allocation - allocate(J(nBas_AOs,nBas_AOs)) - allocate(K(nBas_AOs,nBas_AOs)) + allocate(J(nBas,nBas)) + allocate(K(nBas,nBas)) - allocate(err(nBas_AOs,nBas_AOs)) - allocate(F(nBas_AOs,nBas_AOs)) + allocate(err(nBas,nBas)) + allocate(F(nBas,nBas)) - allocate(cp(nBas_MOs,nBas_MOs)) - allocate(Fp(nBas_MOs,nBas_MOs)) + allocate(cp(nOrb,nOrb)) + allocate(Fp(nOrb,nOrb)) - allocate(err_diis(nBas_AOs_Sq,max_diis)) - allocate(F_diis(nBas_AOs_Sq,max_diis)) + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Guess coefficients and density matrix - call mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) + call mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) !P(:,:) = 2d0 * matmul(c(:,1:nO), transpose(c(:,1:nO))) - call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO, 2.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P, nBas_AOs) + call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & + c(1,1), nBas, c(1,1), nBas, & + 0.d0, P(1,1), nBas) ! Initialization n_diis = 0 F_diis(:,:) = 0d0 err_diis(:,:) = 0d0 - rcond = 0d0 + rcond = 0d0 Conv = 1d0 nSCF = 0 @@ -124,10 +126,14 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Build Fock matrix - call Hartree_matrix_AO_basis(nBas_AOs, P, ERI, J) - call exchange_matrix_AO_basis(nBas_AOs, P, ERI, K) + call Hartree_matrix_AO_basis(nBas, P, ERI, J) + call exchange_matrix_AO_basis(nBas, P, ERI, K) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) + if(nBas .ne. nOrb) then + call AOtoMO(nBas, nOrb, c(1,1), F(1,1), Fp(1,1)) + call MOtoAO(nBas, nOrb, S(1,1), c(1,1), Fp(1,1), F(1,1)) + endif ! Check convergence @@ -136,19 +142,19 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Kinetic energy - ET = trace_matrix(nBas_AOs, matmul(P, T)) + ET = trace_matrix(nBas, matmul(P, T)) ! Potential energy - EV = trace_matrix(nBas_AOs, matmul(P, V)) + EV = trace_matrix(nBas, matmul(P, V)) ! Hartree energy - EJ = 0.5d0*trace_matrix(nBas_AOs, matmul(P, J)) + EJ = 0.5d0*trace_matrix(nBas, matmul(P, J)) ! Exchange energy - EK = 0.25d0*trace_matrix(nBas_AOs, matmul(P, K)) + EK = 0.25d0*trace_matrix(nBas, matmul(P, K)) ! Total energy @@ -159,27 +165,36 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, if(max_diis > 1) then n_diis = min(n_diis+1, max_diis) - call DIIS_extrapolation(rcond, nBas_AOs_Sq, nBas_AOs_Sq, n_diis, err_diis, F_diis, err, F) + call DIIS_extrapolation(rcond, nBas_Sq, nBas_Sq, n_diis, err_diis, F_diis, err, F) end if ! Level shift if(level_shift > 0d0 .and. Conv > thresh) then - call level_shifting(level_shift, nBas_AOs, nBas_MOs, nO, S, c, F) + call level_shifting(level_shift, nBas, nOrb, nO, S, c, F) endif ! Diagonalize Fock matrix - Fp = matmul(transpose(X), matmul(F, X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas_MOs, cp, eHF) - c = matmul(X, cp) + if(nBas .eq. nOrb) then + Fp = matmul(transpose(X), matmul(F, X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nOrb, cp, eHF) + c = matmul(X, cp) + else + Fp = matmul(transpose(c), matmul(F, c)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nOrb, cp, eHF) + c = matmul(c, cp) + endif ! Density matrix !P(:,:) = 2d0*matmul(c(:,1:nO), transpose(c(:,1:nO))) - call dgemm('N', 'T', nBas_AOs, nBas_AOs, nO, 2.d0, c, nBas_AOs, c, nBas_AOs, 0.d0, P, nBas_AOs) + call dgemm('N', 'T', nBas, nBas, nO, 2.d0, & + c(1,1), nBas, c(1,1), nBas, & + 0.d0, P(1,1), nBas) ! Dump results @@ -210,8 +225,8 @@ subroutine RHF(dotest, maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, ! Compute dipole moments - call dipole_moment(nBas_AOs, P, nNuc, ZNuc, rNuc, dipole_int, dipole) - call print_RHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, EK, ERHF, dipole) + call dipole_moment(nBas, P, nNuc, ZNuc, rNuc, dipole_int, dipole) + call print_RHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, ERHF, dipole) ! Testing zone diff --git a/src/HF/RHF_search.f90 b/src/HF/RHF_search.f90 index 714d11c..d200139 100644 --- a/src/HF/RHF_search.f90 +++ b/src/HF/RHF_search.f90 @@ -2,7 +2,7 @@ ! --- 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, & + nBas, nOrb, 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 @@ -16,7 +16,7 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z double precision,intent(in) :: thresh double precision,intent(in) :: level_shift - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -25,15 +25,15 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z double precision,intent(in) :: ZNuc(nNuc) double precision,intent(in) :: rNuc(nNuc,ncart) double precision,intent(in) :: ENuc - 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) + 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,nOrb) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(inout):: ERI_MO(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(inout):: dipole_int_MO(nOrb,nOrb,ncart) ! Local variables @@ -62,9 +62,9 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z ! Output variables double precision,intent(out) :: ERHF - 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) + double precision,intent(out) :: e(nOrb) + double precision,intent(inout):: c(nBas,nOrb) + double precision,intent(out) :: P(nBas,nBas) ! Memory allocation @@ -80,7 +80,7 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z nS = (nO - nC)*(nV - nR) allocate(Aph(nS,nS), Bph(nS,nS), AB(nS,nS), Om(nS)) - allocate(R(nBas_MOs,nBas_MOs), ExpR(nBas_MOs,nBas_MOs)) + allocate(R(nOrb,nOrb), ExpR(nOrb,nOrb)) !------------------! ! Search algorithm ! @@ -97,7 +97,7 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z call wall_time(start_HF) 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) + nBas, nOrb, 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 @@ -113,9 +113,9 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z write(*,*) 'AO to MO transformation... Please be patient' write(*,*) do ixyz = 1, ncart - call AOtoMO(nBas_AOs, nBas_MOs, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) + call AOtoMO(nBas, nOrb, c, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do - call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, c, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, c, ERI_AO, ERI_MO) call wall_time(end_AOtoMO) t_AOtoMO = end_AOtoMO - start_AOtoMO @@ -128,8 +128,8 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z ispin = 1 - 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) + call phLR_A(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,e,ERI_MO,Aph) + call phLR_B(ispin,.false.,nOrb,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) AB(:,:) = Aph(:,:) + Bph(:,:) @@ -169,14 +169,14 @@ subroutine RHF_search(maxSCF, thresh, max_diis, guess_type, level_shift, nNuc, Z R(:,:) = 0d0 ia = 0 do i=nC+1,nO - do a=nO+1,nBas_MOs-nR + do a=nO+1,nOrb-nR ia = ia + 1 R(a,i) = +AB(ia,eig) R(i,a) = -AB(ia,eig) end do end do - call matrix_exponential(nBas_MOs, R, ExpR) + call matrix_exponential(nOrb, R, ExpR) c = matmul(c, ExpR) else diff --git a/src/HF/ROHF.f90 b/src/HF/ROHF.f90 index d67fc0b..5a0ff02 100644 --- a/src/HF/ROHF.f90 +++ b/src/HF/ROHF.f90 @@ -2,7 +2,7 @@ ! --- 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) + nBas, nOrb, nO, S, T, V, Hc, ERI, dipole_int, X, EROHF, eHF, c, Ptot) ! Perform restricted open-shell Hartree-Fock calculation @@ -19,7 +19,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, double precision,intent(in) :: mix double precision,intent(in) :: level_shift double precision,intent(in) :: thresh - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nNuc double precision,intent(in) :: ZNuc(nNuc) @@ -27,18 +27,18 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, double precision,intent(in) :: ENuc integer,intent(in) :: nO(nspin) - 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) + 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,nOrb) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int(nBas,nBas,ncart) ! Local variables integer :: nSCF - integer :: nBas_AOs_Sq + integer :: nBas_Sq integer :: n_diis double precision :: Conv double precision :: rcond @@ -65,9 +65,9 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Output variables double precision,intent(out) :: EROHF - 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) + double precision,intent(out) :: eHF(nOrb) + double precision,intent(inout):: c(nBas,nOrb) + double precision,intent(out) :: Ptot(nBas,nBas) ! Hello world @@ -79,30 +79,30 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Useful stuff - nBas_AOs_Sq = nBas_AOs*nBas_AOs + nBas_Sq = nBas*nBas ! Memory allocation - 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(J(nBas,nBas,nspin)) + allocate(K(nBas,nBas,nspin)) + allocate(F(nBas,nBas,nspin)) + allocate(Ftot(nBas,nBas)) + allocate(P(nBas,nBas,nspin)) + allocate(err(nBas,nBas)) - allocate(Fp(nBas_MOs,nBas_MOs)) - allocate(cp(nBas_MOs,nBas_MOs)) + allocate(Fp(nOrb,nOrb)) + allocate(cp(nOrb,nOrb)) - allocate(err_diis(nBas_AOs_Sq,max_diis)) - allocate(F_diis(nBas_AOs_Sq,max_diis)) + allocate(err_diis(nBas_Sq,max_diis)) + allocate(F_diis(nBas_Sq,max_diis)) ! Guess coefficients and demsity matrices - call mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) + call mo_guess(nBas, nOrb, 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) + call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas) end do Ptot(:,:) = P(:,:,1) + P(:,:,2) @@ -135,13 +135,13 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Build Hartree repulsion do ispin = 1, nspin - call Hartree_matrix_AO_basis(nBas_AOs, P(:,:,ispin), ERI(:,:,:,:), J(:,:,ispin)) + call Hartree_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), J(:,:,ispin)) end do ! Compute exchange potential do ispin = 1, nspin - call exchange_matrix_AO_basis(nBas_AOs, P(:,:,ispin), ERI(:,:,:,:), K(:,:,ispin)) + call exchange_matrix_AO_basis(nBas, P(:,:,ispin), ERI(:,:,:,:), K(:,:,ispin)) end do ! Build Fock operator @@ -150,7 +150,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, F(:,:,ispin) = Hc(:,:) + J(:,:,ispin) + J(:,:,mod(ispin,2)+1) + K(:,:,ispin) end do - call ROHF_fock_matrix(nBas_AOs, nBas_MOs, nO(1), nO(2), S, c, F(:,:,1), F(:,:,2), Ftot) + call ROHF_fock_matrix(nBas, nOrb, nO(1), nO(2), S, c, F(:,:,1), F(:,:,2), Ftot) ! Check convergence @@ -160,25 +160,25 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Kinetic energy do ispin = 1, nspin - ET(ispin) = trace_matrix(nBas_AOs, matmul(P(:,:,ispin), T(:,:))) + ET(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), T(:,:))) end do ! Potential energy do ispin = 1, nspin - EV(ispin) = trace_matrix(nBas_AOs, matmul(P(:,:,ispin), V(:,:))) + EV(ispin) = trace_matrix(nBas, matmul(P(:,:,ispin), V(:,:))) end do ! Hartree energy - 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))) + 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))) ! Exchange energy do ispin = 1, nspin - EK(ispin) = 0.5d0*trace_matrix(nBas_AOs, matmul(P(:,:,ispin), K(:,:,ispin))) + EK(ispin) = 0.5d0*trace_matrix(nBas, matmul(P(:,:,ispin), K(:,:,ispin))) end do ! Total energy @@ -190,7 +190,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, if(max_diis > 1) then n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBas_AOs_Sq,nBas_AOs_Sq,n_diis,err_diis,F_diis,err,Ftot) + call DIIS_extrapolation(rcond,nBas_Sq,nBas_Sq,n_diis,err_diis,F_diis,err,Ftot) end if @@ -199,7 +199,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, if(level_shift > 0d0 .and. Conv > thresh) then do ispin=1,nspin - call level_shifting(level_shift, nBas_AOs, nBas_MOs, maxval(nO), S, c, Ftot) + call level_shifting(level_shift, nBas, nOrb, maxval(nO), S, c, Ftot) end do end if @@ -211,7 +211,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Diagonalize Fock matrix to get eigenvectors and eigenvalues cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas_MOs, cp, eHF) + call diagonalize_matrix(nOrb, cp, eHF) ! Back-transform eigenvectors in non-orthogonal basis @@ -221,7 +221,7 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, 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) + call dgemm('N', 'T', nBas, nBas, nO(ispin), 1.d0, c, nBas, c, nBas, 0.d0, P(1,1,ispin), nBas) end do Ptot(:,:) = P(:,:,1) + P(:,:,2) @@ -254,8 +254,8 @@ subroutine ROHF(dotest, maxSCF, thresh, max_diis, guess_type, mix, level_shift, ! Compute final UHF energy - 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) + call dipole_moment(nBas,Ptot,nNuc,ZNuc,rNuc,dipole_int,dipole) + call print_ROHF(nBas, nOrb, nO, eHF, c, ENuc, ET, EV, EJ, EK, EROHF, dipole) ! Print test values diff --git a/src/HF/ROHF_fock_matrix.f90 b/src/HF/ROHF_fock_matrix.f90 index a3c5aad..a88c57a 100644 --- a/src/HF/ROHF_fock_matrix.f90 +++ b/src/HF/ROHF_fock_matrix.f90 @@ -1,7 +1,7 @@ ! --- -subroutine ROHF_fock_matrix(nBas_AOs, nBas_MOs, nOa, nOb, S, c, FaAO, FbAO, FAO) +subroutine ROHF_fock_matrix(nBas, nOrb, 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 @@ -20,14 +20,14 @@ subroutine ROHF_fock_matrix(nBas_AOs, nBas_MOs, nOa, nOb, S, c, FaAO, FbAO, FAO) ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nOa integer,intent(in) :: nOb - 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) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: c(nBas,nOrb) + double precision,intent(inout):: FaAO(nBas,nBas) + double precision,intent(inout):: FbAO(nBas,nBas) ! Local variables @@ -45,11 +45,11 @@ subroutine ROHF_fock_matrix(nBas_AOs, nBas_MOs, nOa, nOb, S, c, FaAO, FbAO, FAO) ! Output variables - double precision,intent(out) :: FAO(nBas_AOs,nBas_AOs) + double precision,intent(out) :: FAO(nBas,nBas) ! Memory allocation - allocate(F(nBas_MOs,nBas_MOs), Fa(nBas_MOs,nBas_MOs), Fb(nBas_MOs,nBas_MOs)) + allocate(F(nOrb,nOrb), Fa(nOrb,nOrb), Fb(nOrb,nOrb)) ! Roothan canonicalization parameters @@ -66,12 +66,12 @@ subroutine ROHF_fock_matrix(nBas_AOs, nBas_MOs, nOa, nOb, S, c, FaAO, FbAO, FAO) nC = min(nOa, nOb) nO = abs(nOa - nOb) - nV = nBas_AOs - nC - nO + nV = nBas - nC - nO ! Block-by-block Fock matrix - call AOtoMO(nBas_AOs, nBas_MOs, c, FaAO, Fa) - call AOtoMO(nBas_AOs, nBas_MOs, c, FbAO, Fb) + call AOtoMO(nBas, nOrb, c, FaAO, Fa) + call AOtoMO(nBas, nOrb, 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 ) @@ -85,9 +85,9 @@ subroutine ROHF_fock_matrix(nBas_AOs, nBas_MOs, 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_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) + call MOtoAO(nBas, nOrb, S, c, F, FAO) + call MOtoAO(nBas, nOrb, S, c, Fa, FaAO) + call MOtoAO(nBas, nOrb, S, c, Fb, FbAO) deallocate(F, Fa, Fb) diff --git a/src/HF/core_guess.f90 b/src/HF/core_guess.f90 index c48c3ae..8d34444 100644 --- a/src/HF/core_guess.f90 +++ b/src/HF/core_guess.f90 @@ -1,4 +1,4 @@ -subroutine core_guess(nBas_AOs, nBas_MOs, Hc, X, c) +subroutine core_guess(nBas, nOrb, Hc, X, c) ! Core guess of the molecular orbitals for HF calculation @@ -6,9 +6,9 @@ subroutine core_guess(nBas_AOs, nBas_MOs, Hc, X, c) ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs - double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) - double precision,intent(in) :: X(nBas_AOs,nBas_MOs) + integer,intent(in) :: nBas, nOrb + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) ! Local variables @@ -18,17 +18,17 @@ subroutine core_guess(nBas_AOs, nBas_MOs, Hc, X, c) ! Output variables - double precision,intent(out) :: c(nBas_AOs,nBas_MOs) + double precision,intent(out) :: c(nBas,nOrb) ! Memory allocation - allocate(cp(nBas_MOs,nBas_MOs), e(nBas_MOs)) + allocate(cp(nOrb,nOrb), e(nOrb)) ! Core guess cp(:,:) = matmul(transpose(X(:,:)), matmul(Hc(:,:), X(:,:))) - call diagonalize_matrix(nBas_MOs, cp, e) + call diagonalize_matrix(nOrb, cp, e) c(:,:) = matmul(X(:,:), cp(:,:)) deallocate(cp, e) diff --git a/src/HF/huckel_guess.f90 b/src/HF/huckel_guess.f90 index a8e7e52..7afd0a0 100644 --- a/src/HF/huckel_guess.f90 +++ b/src/HF/huckel_guess.f90 @@ -1,4 +1,4 @@ -subroutine huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) +subroutine huckel_guess(nBas, nOrb, S, Hc, X, c) ! Hickel guess @@ -6,10 +6,10 @@ subroutine huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs - double precision,intent(in) :: S(nBas_AOs,nBas_AOs) - double precision,intent(in) :: Hc(nBas_AOs,nBas_AOs) - double precision,intent(in) :: X(nBas_AOs,nBas_MOs) + integer,intent(in) :: nBas, nOrb + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) ! Local variables @@ -20,11 +20,11 @@ subroutine huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) ! Output variables - double precision,intent(out) :: c(nBas_AOs,nBas_MOs) + double precision,intent(out) :: c(nBas,nOrb) ! Memory allocation - allocate(F(nBas_AOs,nBas_AOs)) + allocate(F(nBas,nBas)) ! Extended Huckel parameter @@ -32,9 +32,9 @@ subroutine huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) ! GWH approximation - do mu = 1, nBas_AOs + do mu = 1, nBas F(mu,mu) = Hc(mu,mu) - do nu = mu+1, nBas_AOs + do nu = mu+1, nBas F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu)) F(nu,mu) = F(mu,nu) @@ -42,7 +42,7 @@ subroutine huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) end do end do - call core_guess(nBas_AOs, nBas_MOs, F, X, c) + call core_guess(nBas, nOrb, F, X, c) deallocate(F) diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index fd1f20d..2046569 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -1,7 +1,7 @@ ! --- -subroutine mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) +subroutine mo_guess(nBas, nOrb, guess_type, S, Hc, X, c) ! Guess of the molecular orbitals for HF calculation @@ -9,15 +9,15 @@ subroutine mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: guess_type - double precision,intent(in) :: S(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) :: S(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nOrb) ! Output variables - double precision,intent(inout) :: c(nBas_AOs,nBas_MOs) + double precision,intent(inout) :: c(nBas,nOrb) if(guess_type == 0) then @@ -27,12 +27,12 @@ subroutine mo_guess(nBas_AOs, nBas_MOs, guess_type, S, Hc, X, c) elseif(guess_type == 1) then write(*,*) 'Core guess...' - call core_guess(nBas_AOs, nBas_MOs, Hc, X, c) + call core_guess(nBas, nOrb, Hc, X, c) elseif(guess_type == 2) then write(*,*) 'Huckel guess...' - call huckel_guess(nBas_AOs, nBas_MOs, S, Hc, X, c) + call huckel_guess(nBas, nOrb, S, Hc, X, c) elseif(guess_type == 3) then diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 6790fcf..277db55 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -1,7 +1,7 @@ ! --- -subroutine print_RHF(nBas_AOs, nBas_MOs, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole) +subroutine print_RHF(nBas, nOrb, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERHF, dipole) ! Print one-electron energies and other stuff for G0W0 @@ -10,10 +10,10 @@ subroutine print_RHF(nBas_AOs, nBas_MOs, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERH ! Input variables - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: cHF(nBas_AOs,nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: cHF(nBas,nOrb) double precision,intent(in) :: ENuc double precision,intent(in) :: ET double precision,intent(in) :: EV @@ -78,13 +78,13 @@ subroutine print_RHF(nBas_AOs, nBas_MOs, nO, eHF, cHF, ENuc, ET, EV, EJ, EK, ERH write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' RHF orbital coefficients ' write(*,'(A50)') '---------------------------------------' - call matout(nBas_AOs, nBas_MOs, cHF) + call matout(nBas, nOrb, cHF) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' RHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas_MOs, eHF) + call vecout(nOrb, eHF) write(*,*) end subroutine diff --git a/src/HF/print_ROHF.f90 b/src/HF/print_ROHF.f90 index 09939f6..0297bb4 100644 --- a/src/HF/print_ROHF.f90 +++ b/src/HF/print_ROHF.f90 @@ -1,17 +1,17 @@ ! --- -subroutine print_ROHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, Ex, EROHF, dipole) +subroutine print_ROHF(nBas, nOrb, 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_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO(nspin) - double precision,intent(in) :: eHF(nBas_MOs) - double precision,intent(in) :: c(nBas_AOs,nBas_MOs) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: c(nBas,nOrb) double precision,intent(in) :: ENuc double precision,intent(in) :: ET(nspin) double precision,intent(in) :: EV(nspin) @@ -34,7 +34,7 @@ subroutine print_ROHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, Ex, EROH do ispin=1,nspin if(nO(ispin) > 0) then HOMO(ispin) = eHF(nO(ispin)) - if(nO(ispin) < nBas_MOs) then + if(nO(ispin) < nOrb) then LUMO(ispin) = eHF(nO(ispin)+1) else LUMO(ispin) = 0d0 @@ -105,13 +105,13 @@ subroutine print_ROHF(nBas_AOs, nBas_MOs, nO, eHF, c, ENuc, ET, EV, EJ, Ex, EROH write(*,'(A50)') '-----------------------------------------' write(*,'(A50)') 'ROHF orbital coefficients ' write(*,'(A50)') '-----------------------------------------' - call matout(nBas_AOs, nBas_MOs, c) + call matout(nBas, nOrb, c) write(*,*) end if write(*,'(A50)') '---------------------------------------' write(*,'(A50)') ' ROHF orbital energies (au) ' write(*,'(A50)') '---------------------------------------' - call vecout(nBas_MOs, eHF) + call vecout(nOrb, eHF) write(*,*) end subroutine diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 75514f4..c6a77c9 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -15,7 +15,7 @@ program QuAcK logical :: doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW logical :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh - integer :: nNuc, nBas_AOs, nBas_MOs + integer :: nNuc, nBas, nOrb integer :: nC(nspin) integer :: nO(nspin) integer :: nV(nspin) @@ -69,7 +69,9 @@ program QuAcK logical :: dotest,doRtest,doUtest,doGtest - integer :: i, j + integer :: i, j, j0 + double precision :: acc_d, acc_nd + double precision, allocatable :: tmp1(:,:), tmp2(:,:) !-------------! ! Hello World ! @@ -130,8 +132,8 @@ program QuAcK ! nO = number of occupied orbitals ! ! nV = number of virtual orbitals (see below) ! ! nR = number of Rydberg orbitals ! -! nBas_AOs = number of basis functions in AOs ! -! nBas_MOs = number of basis functions in MOs ! +! nBas = number of basis functions in AOs ! +! nOrb = number of basis functions in MOs ! !---------------------------------------------------! call read_molecule(nNuc,nO,nC,nR) @@ -145,7 +147,7 @@ program QuAcK ! Read basis set information from PySCF ! !---------------------------------------! - call read_basis_pyscf(nBas_AOs, nO, nV) + call read_basis_pyscf(nBas, nO, nV) !--------------------------------------! ! Read one- and two-electron integrals ! @@ -153,19 +155,19 @@ program QuAcK ! Memory allocation for one- and two-electron integrals - allocate(S(nBas_AOs,nBas_AOs)) - allocate(T(nBas_AOs,nBas_AOs)) - allocate(V(nBas_AOs,nBas_AOs)) - allocate(Hc(nBas_AOs,nBas_AOs)) - allocate(ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs)) - allocate(dipole_int_AO(nBas_AOs,nBas_AOs,ncart)) + allocate(S(nBas,nBas)) + allocate(T(nBas,nBas)) + allocate(V(nBas,nBas)) + allocate(Hc(nBas,nBas)) + allocate(ERI_AO(nBas,nBas,nBas,nBas)) + allocate(dipole_int_AO(nBas,nBas,ncart)) ! Read integrals call wall_time(start_int) - call read_integrals(nBas_AOs, S(1,1), T(1,1), V(1,1), Hc(1,1), ERI_AO(1,1,1,1)) - call read_dipole_integrals(nBas_AOs, dipole_int_AO) + call read_integrals(nBas, S(1,1), T(1,1), V(1,1), Hc(1,1), ERI_AO(1,1,1,1)) + call read_dipole_integrals(nBas, dipole_int_AO) call wall_time(end_int) @@ -176,39 +178,61 @@ program QuAcK ! Compute orthogonalization matrix - !call orthogonalization_matrix(nBas_AOs, S, X) + !call orthogonalization_matrix(nBas, S, X) - allocate(Uvec(nBas_AOs,nBas_AOs), Uval(nBas_AOs)) + allocate(Uvec(nBas,nBas), Uval(nBas)) - Uvec(1:nBas_AOs,1:nBas_AOs) = S(1:nBas_AOs,1:nBas_AOs) - call diagonalize_matrix(nBas_AOs, Uvec, Uval) + Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas) + call diagonalize_matrix(nBas, Uvec, Uval) - nBas_MOs = 0 - do i = 1, nBas_AOs + nOrb = 0 + do i = 1, nBas if(Uval(i) > 1d-6) then Uval(i) = 1d0 / dsqrt(Uval(i)) - nBas_MOs = nBas_MOs + 1 + nOrb = nOrb + 1 else write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization' end if end do write(*,'(A38)') '--------------------------------------' - write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas_AOs - write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nBas_MOs - write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nBas_MOs)/dble(nBas_AOs)) + write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas + write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nOrb + write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nOrb)/dble(nBas)) write(*,'(A38)') '--------------------------------------' write(*,*) - allocate(X(nBas_AOs,nBas_MOs)) - do j = 1, nBas_MOs - do i = 1, nBas_AOs - X(i,j) = Uvec(i,j) * Uval(j) + j0 = nBas - nOrb + allocate(X(nBas,nOrb)) + do j = j0+1, nBas + do i = 1, nBas + X(i,j-j0) = Uvec(i,j) * Uval(j) enddo enddo deallocate(Uvec, Uval) + !! check if X.T S X = 1_(nOrb,nOrb) + !allocate(tmp1(nOrb,nBas), tmp2(nOrb,nOrb)) + !call dgemm("T", "N", nOrb, nBas, nBas, 1.d0, & + ! X(1,1), nBas, S(1,1), nBas, & + ! 0.d0, tmp1(1,1), nOrb) + !call dgemm("N", "N", nOrb, nOrb, nBas, 1.d0, & + ! tmp1(1,1), nOrb, X(1,1), nBas, & + ! 0.d0, tmp2(1,1), nOrb) + !acc_d = 0.d0 + !acc_nd = 0.d0 + !do i = 1, nOrb + ! !write(*,'(1000(F15.7,2X))') (tmp2(i,j), j = 1, nOrb) + ! acc_d = acc_d + tmp2(i,i) + ! do j = 1, nOrb + ! if(j == i) cycle + ! acc_nd = acc_nd + dabs(tmp2(j,i)) + ! enddo + !enddo + !print*, ' diag part: ', dabs(acc_d - dble(nOrb)) / dble(nOrb) + !print*, ' non-diag part: ', acc_nd + !deallocate(tmp1, tmp2) !---------------------! ! Choose QuAcK branch ! @@ -240,7 +264,7 @@ program QuAcK dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas_AOs,nBas_MOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & @@ -256,7 +280,7 @@ program QuAcK dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas_AOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,spin_conserved,spin_flip,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & @@ -271,7 +295,7 @@ program QuAcK call GQuAcK(doGtest,doGHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,doDCD,doCCSD,doCCSDT, & dodrCCD,dorCCD,docrCCD,dolCCD,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0W0,doevGW,doqsGW,doG0F2,doevGF2,doqsGF2, & - nNuc,nBas_AOs,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & + nNuc,nBas,sum(nC),sum(nO),sum(nV),sum(nR),ENuc,ZNuc,rNuc,S,T,V,Hc,X,dipole_int_AO,ERI_AO, & maxSCF_HF,max_diis_HF,thresh_HF,level_shift,guess_type,mix,reg_MP, & maxSCF_CC,max_diis_CC,thresh_CC,TDA,maxSCF_GF,max_diis_GF,thresh_GF,lin_GF,reg_GF,eta_GF, & maxSCF_GW,max_diis_GW,thresh_GW,TDA_W,lin_GW,reg_GW,eta_GW, & diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index db16eab..1390366 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -2,7 +2,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d dodrCCD,dorCCD,docrCCD,dolCCD,doCIS,doCIS_D,doCID,doCISD,doFCI,dophRPA,dophRPAx,docrRPA,doppRPA, & doG0F2,doevGF2,doqsGF2,doufG0F02,doG0F3,doevGF3,doG0W0,doevGW,doqsGW,doufG0W0,doufGW,doSRGqsGW, & doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp,doG0T0eh,doevGTeh,doqsGTeh, & - nNuc,nBas_AOs,nBas_MOs,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & + nNuc,nBas,nOrb,nC,nO,nV,nR,ENuc,ZNuc,rNuc, & S,T,V,Hc,X,dipole_int_AO,ERI_AO,maxSCF_HF,max_diis_HF,thresh_HF,level_shift, & guess_type,mix,reg_MP,maxSCF_CC,max_diis_CC,thresh_CC,singlet,triplet,TDA, & maxSCF_GF,max_diis_GF,renorm_GF,thresh_GF,lin_GF,reg_GF,eta_GF,maxSCF_GW,max_diis_GW,thresh_GW, & @@ -29,7 +29,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d logical,intent(in) :: doG0T0pp,doevGTpp,doqsGTpp,doufG0T0pp logical,intent(in) :: doG0T0eh,doevGTeh,doqsGTeh - integer,intent(in) :: nNuc,nBas_AOs,nBas_MOs + integer,intent(in) :: nNuc,nBas,nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -38,13 +38,13 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d double precision,intent(in) :: ZNuc(nNuc),rNuc(nNuc,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) :: dipole_int_AO(nBas_AOs,nBas_AOs,ncart) - double precision,intent(in) :: ERI_AO(nBas_AOs,nBas_AOs,nBas_AOs,nBas_AOs) + 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,nOrb) + double precision,intent(in) :: dipole_int_AO(nBas,nBas,ncart) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) integer,intent(in) :: maxSCF_HF,max_diis_HF double precision,intent(in) :: thresh_HF,level_shift,mix @@ -109,11 +109,11 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Memory allocation ! !-------------------! - allocate(cHF(nBas_AOs,nBas_MOs)) - allocate(eHF(nBas_MOs)) - allocate(PHF(nBas_AOs,nBas_AOs)) - allocate(dipole_int_MO(nBas_MOs,nBas_MOs,ncart)) - allocate(ERI_MO(nBas_MOs,nBas_MOs,nBas_MOs,nBas_MOs)) + allocate(cHF(nBas,nOrb)) + allocate(eHF(nOrb)) + allocate(PHF(nBas,nBas)) + allocate(dipole_int_MO(nOrb,nOrb,ncart)) + allocate(ERI_MO(nOrb,nOrb,nOrb,nOrb)) !---------------------! ! Hartree-Fock module ! @@ -123,7 +123,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_HF) 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) + nBas, nOrb, 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 @@ -136,7 +136,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_HF) 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) + nBas, nOrb, 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 @@ -158,12 +158,12 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d ! Read and transform dipole-related integrals do ixyz = 1, ncart - call AOtoMO(nBas_AOs, nBas_MOs, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) + call AOtoMO(nBas, nOrb, cHF, dipole_int_AO(1,1,ixyz), dipole_int_MO(1,1,ixyz)) end do ! 4-index transform - call AOtoMO_ERI_RHF(nBas_AOs, nBas_MOs, cHF, ERI_AO, ERI_MO) + call AOtoMO_ERI_RHF(nBas, nOrb, cHF, ERI_AO, ERI_MO) call wall_time(end_AOtoMO) @@ -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_MOs, nC, nO, nV, nR, nS, eHF, ERI_MO) + call RHF_stability(nOrb, nC, nO, nV, nR, nS, eHF, ERI_MO) call wall_time(end_stab) t_stab = end_stab - start_stab @@ -193,7 +193,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_stab) 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, & + nBas, nOrb, 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) @@ -212,7 +212,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doMP) then call wall_time(start_MP) - call RMP(dotest, doMP2, doMP3, reg_MP, nBas_MOs, nBas_MOs, nC, nO, nV, nR, ERI_MO, ENuc, ERHF, eHF) + call RMP(dotest, doMP2, doMP3, reg_MP, nOrb, nOrb, nC, nO, nV, nR, ERI_MO, ENuc, ERHF, eHF) call wall_time(end_MP) t_MP = end_MP - start_MP @@ -232,7 +232,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_CC) call RCC(dotest, doCCD, dopCCD, doDCD, doCCSD, doCCSDT, dodrCCD, dorCCD, docrCCD, dolCCD, & - maxSCF_CC, thresh_CC, max_diis_CC, nBas_AOs, nBas_MOs, nC, nO, nV, nR, Hc, ERI_MO, ENuc, ERHF, eHF, cHF) + maxSCF_CC, thresh_CC, max_diis_CC, nBas, nOrb, nC, nO, nV, nR, Hc, ERI_MO, ENuc, ERHF, eHF, cHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -250,7 +250,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d if(doCI) then call wall_time(start_CI) - call RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nBas_MOs, & + call RCI(dotest, doCIS, doCIS_D, doCID, doCISD, doFCI, singlet, triplet, nOrb, & nC, nO, nV, nR, nS, ERI_MO, dipole_int_MO, eHF, ERHF) call wall_time(end_CI) @@ -270,7 +270,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_RPA) call RRPA(dotest, dophRPA, dophRPAx, docrRPA, doppRPA, TDA, doACFDT, exchange_kernel, singlet, triplet, & - nBas_MOs, nC, nO, nV, nR, nS, ENuc, ERHF, ERI_MO, dipole_int_MO, eHF) + nOrb, nC, nO, nV, nR, nS, ENuc, ERHF, ERI_MO, dipole_int_MO, eHF) call wall_time(end_RPA) t_RPA = end_RPA - start_RPA @@ -290,7 +290,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_GF) call RGF(dotest, doG0F2, doevGF2, doqsGF2, doufG0F02, doG0F3, doevGF3, renorm_GF, maxSCF_GF, & thresh_GF, max_diis_GF, dophBSE, doppBSE, TDA, dBSE, dTDA, singlet, triplet, lin_GF, & - eta_GF, reg_GF, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, & + eta_GF, reg_GF, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, & S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GF) @@ -311,7 +311,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_GW) call RGW(dotest, doG0W0, doevGW, doqsGW, doufG0W0, doufGW, doSRGqsGW, maxSCF_GW, thresh_GW, max_diis_GW, & doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, TDA_W, TDA, dBSE, dTDA, singlet, triplet, & - lin_GW, eta_GW, reg_GW, nNuc, ZNuc, rNuc, ENuc, nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, & + lin_GW, eta_GW, reg_GW, nNuc, ZNuc, rNuc, ENuc, nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, & V, Hc, ERI_AO, ERI_MO, dipole_int_AO, dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GW) @@ -333,7 +333,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call RGT(dotest, doG0T0pp, doevGTpp, doqsGTpp, doufG0T0pp, doG0T0eh, doevGTeh, doqsGTeh, & maxSCF_GT, thresh_GT, max_diis_GT, doACFDT, exchange_kernel, doXBS, dophBSE, dophBSE2, doppBSE, & TDA_T, TDA, dBSE, dTDA, singlet, triplet, lin_GT, eta_GT, reg_GT, nNuc, ZNuc, rNuc, ENuc, & - nBas_AOs, nBas_MOs, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, & + nBas, nOrb, nC, nO, nV, nR, nS, ERHF, S, X, T, V, Hc, ERI_AO, ERI_MO, dipole_int_AO, & dipole_int_MO, PHF, cHF, eHF) call wall_time(end_GT) diff --git a/src/utils/level_shifting.f90 b/src/utils/level_shifting.f90 index a006622..8b627fa 100644 --- a/src/utils/level_shifting.f90 +++ b/src/utils/level_shifting.f90 @@ -1,4 +1,4 @@ -subroutine level_shifting(level_shift, nBas_AOs, nBas_MOs, nO, S, c, F) +subroutine level_shifting(level_shift, nBas, nOrb, nO, S, c, F) ! Perform level-shifting on the Fock matrix @@ -7,10 +7,10 @@ subroutine level_shifting(level_shift, nBas_AOs, nBas_MOs, nO, S, c, F) ! Input variables double precision,intent(in) :: level_shift - integer,intent(in) :: nBas_AOs, nBas_MOs + integer,intent(in) :: nBas, nOrb integer,intent(in) :: nO - double precision,intent(in) :: S(nBas_AOs,nBas_AOs) - double precision,intent(in) :: c(nBas_AOs,nBas_MOs) + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: c(nBas,nOrb) ! Local variables @@ -21,13 +21,13 @@ subroutine level_shifting(level_shift, nBas_AOs, nBas_MOs, nO, S, c, F) ! Output variables - double precision,intent(inout):: F(nBas_AOs,nBas_AOs) + double precision,intent(inout):: F(nBas,nBas) - allocate(F_MO(nBas_MOs,nBas_MOs), Sc(nBas_AOs,nBas_MOs)) + allocate(F_MO(nOrb,nOrb), Sc(nBas,nOrb)) F_MO(:,:) = matmul(transpose(c), matmul(F, c)) - do a = nO+1, nBas_MOs + do a = nO+1, nOrb F_MO(a,a) = F_MO(a,a) + level_shift end do