diff --git a/input/dft b/input/dft index b16601a..aaee5fd 100644 --- a/input/dft +++ b/input/dft @@ -31,7 +31,7 @@ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 1 0.0 0.0 + 0.0 0.0 0.0 # Ncentered ? F # Parameters for CC weight-dependent exchange functional diff --git a/input/methods b/input/methods index 979b54d..9b32026 100644 --- a/input/methods +++ b/input/methods @@ -1,9 +1,9 @@ # RHF UHF KS MOM - F F T F + T F F F # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) - F F F F F + F F F T F # drCCD rCCD crCCD lCCD F F F F # CIS* CIS(D) CID CISD FCI @@ -13,9 +13,9 @@ # G0F2* evGF2* qsGF2* G0F3 evGF3 F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F F F F F + T F F F F # G0T0 evGT qsGT - F F F + T F F # MCMP2 F # * unrestricted version available diff --git a/input/options b/input/options index f4a2779..2dc2665 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ -# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability - 1024 0.00001 T 2 1 1 F F +# HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability + 256 0.0000001 T 5 2 1 T 1.0 F # MP: # CC: maxSCF thresh DIIS n_diis @@ -7,14 +7,14 @@ # spin: TDA singlet triplet spin_conserved spin_flip F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg - 256 0.00001 T 5 T 0.0 3 F + 256 0.00001 T 5 T 0.0 3 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg - 256 0.00001 T 5 T 0.0 F F F F F F + 256 0.00001 T 5 T 0.0 F F F F F F # GT: maxSCF thresh DIIS n_diis lin eta TDA_T reg 256 0.00001 T 5 T 0.0 F F # ACFDT: AC Kx XBS - F F F + F T T # BSE: BSE dBSE dTDA evDyn - T T T F + T F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/mol/h2.xyz b/mol/h2.xyz index 6a4e902..9e53116 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 -H 0.0 0.0 0.0 -H 0.0 0.0 0.7 +H 0. 0. 0. +H 0. 0. 0.741 diff --git a/src/CI/CISD.f90 b/src/CI/CISD.f90 index 079521c..02f7fd3 100644 --- a/src/CI/CISD.f90 +++ b/src/CI/CISD.f90 @@ -86,7 +86,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI write(*,*) 'nH = ',nH write(*,*) - maxH = min(nH,41) + maxH = min(nH,51) ! Memory allocation diff --git a/src/GF/BSE2.f90 b/src/GF/BSE2.f90 index 25bc56c..e971ceb 100644 --- a/src/GF/BSE2.f90 +++ b/src/GF/BSE2.f90 @@ -1,6 +1,6 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,EcBSE) -! Compute the Bethe-Salpeter excitation energies +! Compute the second-order Bethe-Salpeter excitation energies implicit none include 'parameters.h' @@ -33,6 +33,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) double precision :: rho + double precision,allocatable :: A_sta(:,:,:) + double precision,allocatable :: B_sta(:,:,:) ! Output variables @@ -40,7 +42,8 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ! Memory allocation - allocate(OmBSE(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(OmBSE(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin), & + A_sta(nS,nS,nspin),B_sta(nS,nS,nspin)) !------------------- ! Singlet manifold @@ -51,13 +54,17 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ispin = 1 EcBSE(ispin) = 0d0 + ! Compute static kernel + + call BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin)) + if(.not.TDA) call BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin)) + ! Compute BSE2 excitation energies - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI, & - OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta(:,:,ispin),-B_sta(:,:,ispin), & + EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) - call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -67,11 +74,11 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, if(evDyn) then call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) else call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if @@ -88,13 +95,17 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, ispin = 2 EcBSE(ispin) = 0d0 + ! Compute static kernel + + call BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,A_sta(:,:,ispin)) + if(.not.TDA) call BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,eGF,B_sta(:,:,ispin)) + ! Compute BSE2 excitation energies - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF(:),ERI(:,:,:,:), & - OmBSE(:,ispin),rho,EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGF,ERI,-A_sta(:,:,ispin),-B_sta(:,:,ispin), & + EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('BSE2 ',ispin,nS,OmBSE(:,ispin)) - call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) ! Compute dynamic correction for BSE via perturbation theory @@ -103,11 +114,11 @@ subroutine BSE2(TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI, if(evDyn) then call BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) else call BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF, & - OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + A_sta(:,:,ispin),B_sta(:,:,ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) end if diff --git a/src/GF/BSE2_A_matrix_static.f90 b/src/GF/BSE2_A_matrix_static.f90 new file mode 100644 index 0000000..dee7202 --- /dev/null +++ b/src/GF/BSE2_A_matrix_static.f90 @@ -0,0 +1,157 @@ +subroutine BSE2_A_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,A_sta) + +! Compute the resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: A_sta(nS,nS) + +! Initialization + + A_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the singlet manifold + + if(ispin == 1) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = - (eGF(c) - eGF(k)) + num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) & + - ERI(j,k,c,i)*ERI(a,c,b,k) + 2d0*ERI(j,k,c,i)*ERI(a,c,k,b) + + A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = + (eGF(c) - eGF(k)) + num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) & + - ERI(j,c,k,i)*ERI(a,k,b,c) + 2d0*ERI(j,c,k,i)*ERI(a,k,c,b) + + A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - (eGF(c) + eGF(d)) + num = 2d0*ERI(a,j,c,d)*ERI(c,d,i,b) - ERI(a,j,c,d)*ERI(c,d,b,i) & + - ERI(a,j,d,c)*ERI(c,d,i,b) + 2d0*ERI(a,j,d,c)*ERI(c,d,b,i) + + A_sta(ia,jb) = A_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = - (eGF(k) + eGF(l)) + num = 2d0*ERI(a,j,k,l)*ERI(k,l,i,b) - ERI(a,j,k,l)*ERI(k,l,b,i) & + - ERI(a,j,l,k)*ERI(k,l,i,b) + 2d0*ERI(a,j,l,k)*ERI(k,l,b,i) + + A_sta(ia,jb) = A_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +! Second-order correlation kernel for the block A of the triplet manifold + + if(ispin == 2) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = - (eGF(c) - eGF(k)) + num = 2d0*ERI(j,k,i,c)*ERI(a,c,b,k) - ERI(j,k,i,c)*ERI(a,c,k,b) - ERI(j,k,c,i)*ERI(a,c,b,k) + + A_sta(ia,jb) = A_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = + (eGF(c) - eGF(k)) + num = 2d0*ERI(j,c,i,k)*ERI(a,k,b,c) - ERI(j,c,i,k)*ERI(a,k,c,b) - ERI(j,c,k,i)*ERI(a,k,b,c) + + A_sta(ia,jb) = A_sta(ia,jb) + num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - (eGF(c) + eGF(d)) + num = ERI(a,j,c,d)*ERI(c,d,b,i) + ERI(a,j,d,c)*ERI(c,d,i,b) + + A_sta(ia,jb) = A_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = - (eGF(k) + eGF(l)) + num = ERI(a,j,k,l)*ERI(k,l,b,i) + ERI(a,j,l,k)*ERI(k,l,i,b) + + A_sta(ia,jb) = A_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + + +end subroutine BSE2_A_matrix_static diff --git a/src/GF/BSE2_B_matrix_static.f90 b/src/GF/BSE2_B_matrix_static.f90 new file mode 100644 index 0000000..e93abba --- /dev/null +++ b/src/GF/BSE2_B_matrix_static.f90 @@ -0,0 +1,157 @@ +subroutine BSE2_B_matrix_static(ispin,eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,eGF,B_sta) + +! Compute the anti-resonant part of the static BSE2 matrix + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: eGF(nBas) + +! Local variables + + double precision :: dem,num + integer :: i,j,k,l + integer :: a,b,c,d + integer :: ia,jb + +! Output variables + + double precision,intent(out) :: B_sta(nS,nS) + +! Initialization + + B_sta(:,:) = 0d0 + +! Second-order correlation kernel for the block A of the singlet manifold + + if(ispin == 1) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = + eGF(k) - eGF(c) + num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) & + - ERI(b,k,c,i)*ERI(a,c,j,k) + 2d0*ERI(b,k,c,i)*ERI(a,c,k,j) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = - eGF(c) + eGF(k) + num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) & + - ERI(b,c,k,i)*ERI(a,k,j,c) + 2d0*ERI(b,c,k,i)*ERI(a,k,c,j) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - eGF(c) - eGF(d) + num = 2d0*ERI(a,b,c,d)*ERI(c,d,i,j) - ERI(a,b,c,d)*ERI(c,d,j,i) & + - ERI(a,b,d,c)*ERI(c,d,i,j) + 2d0*ERI(a,b,d,c)*ERI(c,d,j,i) + + B_sta(ia,jb) = B_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = + eGF(k) + eGF(l) + num = 2d0*ERI(a,b,k,l)*ERI(k,l,i,j) - ERI(a,b,k,l)*ERI(k,l,j,i) & + - ERI(a,b,l,k)*ERI(k,l,i,j) + 2d0*ERI(a,b,l,k)*ERI(k,l,j,i) + + B_sta(ia,jb) = B_sta(ia,jb) + 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + +! Second-order correlation kernel for the block A of the triplet manifold + + if(ispin == 2) then + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + do k=nC+1,nO + do c=nO+1,nBas-nR + + dem = + eGF(k) - eGF(c) + num = 2d0*ERI(b,k,i,c)*ERI(a,c,j,k) - ERI(b,k,i,c)*ERI(a,c,k,j) - ERI(b,k,c,i)*ERI(a,c,j,k) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + dem = - eGF(c) + eGF(k) + num = 2d0*ERI(b,c,i,k)*ERI(a,k,j,c) - ERI(b,c,i,k)*ERI(a,k,c,j) - ERI(b,c,k,i)*ERI(a,k,j,c) + + B_sta(ia,jb) = B_sta(ia,jb) - num*dem/(dem**2 + eta**2) + + end do + end do + + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + dem = - eGF(c) - eGF(d) + num = ERI(a,b,c,d)*ERI(c,d,j,i) + ERI(a,b,d,c)*ERI(c,d,i,j) + + B_sta(ia,jb) = B_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + do k=nC+1,nO + do l=nC+1,nO + + dem = + eGF(k) + eGF(l) + num = ERI(a,b,k,l)*ERI(k,l,j,i) + ERI(a,b,l,k)*ERI(k,l,i,j) + + B_sta(ia,jb) = B_sta(ia,jb) - 0.5d0*num*dem/(dem**2 + eta**2) + + end do + end do + + end do + end do + + end do + end do + + end if + + +end subroutine BSE2_B_matrix_static diff --git a/src/GF/BSE2_dynamic_perturbation.f90 b/src/GF/BSE2_dynamic_perturbation.f90 index 94c7c59..50077f1 100644 --- a/src/GF/BSE2_dynamic_perturbation.f90 +++ b/src/GF/BSE2_dynamic_perturbation.f90 @@ -1,4 +1,4 @@ -subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,OmBSE,XpY,XmY) +subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY) ! Compute dynamical effects via perturbation theory for BSE @@ -21,6 +21,8 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: A_sta(nS,nS) + double precision,intent(in) :: B_sta(nS,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -81,7 +83,7 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo if(dTDA) then ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) else @@ -94,10 +96,10 @@ subroutine BSE2_dynamic_perturbation(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipo ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + dot_product(Y,matmul(ZAm_dyn,Y)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & - - dot_product(Y,matmul(Am_dyn,Y)) & - + dot_product(X,matmul(B_dyn,Y)) & - - dot_product(Y,matmul(B_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) & + - dot_product(Y,matmul(Am_dyn - A_sta,Y)) & + + dot_product(X,matmul(B_dyn - B_sta,Y)) & + - dot_product(Y,matmul(B_dyn - B_sta,X)) end if diff --git a/src/GF/BSE2_dynamic_perturbation_iterative.f90 b/src/GF/BSE2_dynamic_perturbation_iterative.f90 index eedfad1..9147bcd 100644 --- a/src/GF/BSE2_dynamic_perturbation_iterative.f90 +++ b/src/GF/BSE2_dynamic_perturbation_iterative.f90 @@ -1,5 +1,5 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int, & - eHF,eGF,OmBSE,XpY,XmY) + eHF,eGF,A_sta,B_sta,OmBSE,XpY,XmY) ! Compute self-consistently the dynamical effects via perturbation theory for BSE2 @@ -22,6 +22,8 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n double precision,intent(in) :: dipole_int(nBas,nBas,ncart) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: eGF(nBas) + double precision,intent(in) :: A_sta(nS,nS) + double precision,intent(in) :: B_sta(nS,nS) double precision,intent(in) :: OmBSE(nS) double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) @@ -102,7 +104,7 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n if(dTDA) then - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) else @@ -116,10 +118,10 @@ subroutine BSE2_dynamic_perturbation_iterative(dTDA,ispin,eta,nBas,nC,nO,nV,nR,n ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & + dot_product(Y,matmul(ZAm_dyn,Y)) - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & - - dot_product(Y,matmul(Am_dyn,Y)) & - + dot_product(X,matmul(B_dyn,Y)) & - - dot_product(Y,matmul(B_dyn,X)) + OmDyn(ia) = dot_product(X,matmul(Ap_dyn - A_sta,X)) & + - dot_product(Y,matmul(Am_dyn - A_sta,Y)) & + + dot_product(X,matmul(B_dyn - B_sta,Y)) & + - dot_product(Y,matmul(B_dyn - B_sta,X)) end if diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90 index b038740..ae8bd95 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix.f90 @@ -124,8 +124,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE singlet excitation energies - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt+TAs,TBt+TBs, & - EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt+TAs,TBt+TBs, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & @@ -163,8 +163,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE triplet excitation energies - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt-TAs,TBt-TBs, & - EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt-TAs,TBt-TBs, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) diff --git a/src/GT/Bethe_Salpeter_Tmatrix_so.f90 b/src/GT/Bethe_Salpeter_Tmatrix_so.f90 new file mode 100644 index 0000000..d990fdb --- /dev/null +++ b/src/GT/Bethe_Salpeter_Tmatrix_so.f90 @@ -0,0 +1,79 @@ +subroutine Bethe_Salpeter_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,Omega1,X1,Y1,Omega2,X2,Y2,rho1,rho2, & + ERI,eT,eGT,EcBSE) + +! Compute the Bethe-Salpeter excitation energies with the T-matrix kernel + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + integer,intent(in) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: eT(nBas) + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: ispin + + double precision :: EcRPA + double precision,allocatable :: TA(:,:),TB(:,:) + double precision,allocatable :: OmBSE(:) + double precision,allocatable :: XpY_BSE(:,:) + double precision,allocatable :: XmY_BSE(:,:) + +! Output variables + + double precision,intent(out) :: EcBSE + +! Memory allocation + + allocate(TA(nS,nS),TB(nS,nS),OmBSE(nS),XpY_BSE(nS,nS),XmY_BSE(nS,nS)) + +!------------------! +! Compute T-matrix ! +!------------------! + + ispin = 4 + + call linear_response_pp(ispin,.false.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eT,ERI, & + Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) + + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Omega1,rho1,Omega2,rho2,TA) + call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,Omega1,rho1,Omega2,rho2,TB) + +!------------------! +! Singlet manifold ! +!------------------! + + ispin = 3 + EcBSE = 0d0 + + ! Compute BSE singlet excitation energies + + call linear_response_BSE(ispin,.false.,.false.,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & + EcBSE,OmBSE,XpY_BSE,XmY_BSE) + + call print_excitation('BSE@GT ',ispin,nS,OmBSE) + +end subroutine Bethe_Salpeter_Tmatrix_so diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index 11f4415..e75ba9d 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -238,11 +238,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! Free memory - - deallocate(Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,rho1t,rho2t) - ! Compute the BSE correlation energy via the adiabatic connection if(doACFDT) then @@ -260,7 +255,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing end if call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI_MO,eHF,eG0T0,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eHF,eG0T0,EcAC) if(exchange_kernel) then diff --git a/src/GT/evGT.f90 b/src/GT/evGT.f90 index 1c5e518..2af3bdd 100644 --- a/src/GT/evGT.f90 +++ b/src/GT/evGT.f90 @@ -265,7 +265,8 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & end if call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI_MO,eGT,eGT,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) if(exchange_kernel) then diff --git a/src/GT/excitation_density_Tmatrix_so.f90 b/src/GT/excitation_density_Tmatrix_so.f90 new file mode 100644 index 0000000..f0749f6 --- /dev/null +++ b/src/GT/excitation_density_Tmatrix_so.f90 @@ -0,0 +1,80 @@ +subroutine excitation_density_Tmatrix_so(nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,rho1,X2,Y2,rho2) + +! Compute excitation densities for T-matrix self-energy + + implicit none + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: X1(nVV,nVV) + double precision,intent(in) :: Y1(nOO,nVV) + double precision,intent(in) :: X2(nVV,nOO) + double precision,intent(in) :: Y2(nOO,nOO) + +! Local variables + + integer :: j,k,l + integer :: b,c,d + integer :: p,q + integer :: ab,cd,ij,kl + double precision,external :: Kronecker_delta + +! Output variables + + double precision,intent(out) :: rho1(nBas,nBas,nVV) + double precision,intent(out) :: rho2(nBas,nBAs,nOO) + +! Initialization + + rho1(:,:,:) = 0d0 + rho2(:,:,:) = 0d0 + + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + + do ab=1,nVV + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X1(cd,ab) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho1(p,q,ab) = rho1(p,q,ab) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y1(kl,ab) + end do + end do + + end do + + do ij=1,nOO + + cd = 0 + do c=nO+1,nBas-nR + do d=c+1,nBas-nR + cd = cd + 1 + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,c,d) - ERI(p,q,d,c))*X2(cd,ij) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=k+1,nO + kl = kl + 1 + rho2(p,q,ij) = rho2(p,q,ij) + (ERI(p,q,k,l) - ERI(p,q,l,k))*Y2(kl,ij) + end do + end do + + end do + end do + + end do + +end subroutine excitation_density_Tmatrix_so diff --git a/src/GT/qsGT.f90 b/src/GT/qsGT.f90 index 1e4c349..4344878 100644 --- a/src/GT/qsGT.f90 +++ b/src/GT/qsGT.f90 @@ -389,7 +389,8 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T end if call ACFDT_Tmatrix(exchange_kernel,doXBS,.false.,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI_MO,eGT,eGT,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI_MO,eGT,eGT,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/GT/renormalization_factor_Tmatrix_so.f90 b/src/GT/renormalization_factor_Tmatrix_so.f90 new file mode 100644 index 0000000..ca98217 --- /dev/null +++ b/src/GT/renormalization_factor_Tmatrix_so.f90 @@ -0,0 +1,63 @@ +subroutine renormalization_factor_Tmatrix_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z) + +! Compute renormalization factor of the T-matrix self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,j,k,l,a,b,c,d,p,cd,kl + double precision :: eps + +! Output variables + + double precision,intent(out) :: Z(nBas) + +! Initialize + + Z(:) = 0d0 + +!---------------------------------------------- +! T-matrix renormalization factor in the spinorbital basis +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do cd=1,nVV + eps = e(p) + e(i) - Omega1(cd) + Z(p) = Z(p) + (rho1(p,i,cd)/eps)**2 + enddo + enddo + enddo + +! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do kl=1,nOO + eps = e(p) + e(a) - Omega2(kl) + Z(p) = Z(p) + (rho2(p,a,kl)/eps)**2 + enddo + enddo + enddo + +! Compute renormalization factor from derivative of SigT + + Z(:) = 1d0/(1d0 + Z(:)) + +end subroutine renormalization_factor_Tmatrix_so diff --git a/src/GT/self_energy_Tmatrix_diag_so.f90 b/src/GT/self_energy_Tmatrix_diag_so.f90 new file mode 100644 index 0000000..9c7834b --- /dev/null +++ b/src/GT/self_energy_Tmatrix_diag_so.f90 @@ -0,0 +1,63 @@ +subroutine self_energy_Tmatrix_diag_so(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT) + +! Compute diagonal of the correlation part of the T-matrix self-energy + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + integer :: i,j,k,l,a,b,c,d,p,cd,kl + double precision :: eps + +! Output variables + + double precision,intent(out) :: SigT(nBas) + +! Initialize + + SigT(:) = 0d0 + +!---------------------------------------------- +! T-matrix self-energy in the spinorbital basis +!---------------------------------------------- + +! Occupied part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do i=nC+1,nO + do cd=1,nVV + eps = e(p) + e(i) - Omega1(cd) + SigT(p) = SigT(p) + rho1(p,i,cd)**2/eps + enddo + enddo + enddo + +! Virtual part of the T-matrix self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do kl=1,nOO + eps = e(p) + e(a) - Omega2(kl) + SigT(p) = SigT(p) + rho2(p,a,kl)**2/eps + enddo + enddo + enddo + +end subroutine self_energy_Tmatrix_diag_so diff --git a/src/GT/soG0T0.f90 b/src/GT/soG0T0.f90 new file mode 100644 index 0000000..540ed47 --- /dev/null +++ b/src/GT/soG0T0.f90 @@ -0,0 +1,124 @@ +subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) + +! Perform G0W0 calculation with a T-matrix self-energy (G0T0) in the spinorbital basis + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + + integer,intent(in) :: nBas,nC,nO,nV,nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: ispin + integer :: nOO + integer :: nVV + double precision :: EcRPA + double precision :: EcGM + double precision :: EcBSE + integer :: nBas2,nC2,nO2,nV2,nR2,nS2 + double precision,allocatable :: Omega1(:) + double precision,allocatable :: X1(:,:) + double precision,allocatable :: Y1(:,:) + double precision,allocatable :: rho1(:,:,:) + double precision,allocatable :: Omega2(:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: Y2(:,:) + double precision,allocatable :: rho2(:,:,:) + double precision,allocatable :: SigT(:) + double precision,allocatable :: Z(:) + double precision,allocatable :: eG0T0(:) + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'| One-shot soG0T0 calculation |' + write(*,*)'************************************************' + write(*,*) + +! Define occupied and virtual spaces + + nBas2 = 2*nBas + nO2 = 2*nO + nV2 = 2*nV + nC2 = 2*nC + nR2 = 2*nR + nS2 = nO2*nV2 + + +! Spatial to spin orbitals + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Dimensions of the rr-RPA linear reponse matrices + + nOO = nO2*(nO2 - 1)/2 + nVV = nV2*(nV2 - 1)/2 + +! Memory allocation + + allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & + Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & + rho1(nBas2,nBas2,nVV),rho2(nBas2,nBas2,nOO), & + eG0T0(nBas2),SigT(nBas2),Z(nBas2)) + +!---------------------------------------------- +! Spinorbital basis +!---------------------------------------------- + + ispin = 4 + +! Compute linear response + + call linear_response_pp(ispin,.false.,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,1d0,seHF,sERI, & + Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) + + call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) + call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) + +! Compute excitation densities for the T-matrix + + call excitation_density_Tmatrix_so(nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI,X1,Y1,rho1,X2,Y2,rho2) + +!---------------------------------------------- +! Compute T-matrix version of the self-energy +!---------------------------------------------- + + call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF, & + Omega1,rho1,Omega2,rho2,SigT) + +! Compute renormalization factor for T-matrix self-energy + + call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF, & + Omega1,rho1,Omega2,rho2,Z) + +!---------------------------------------------- +! Solve the quasi-particle equation +!---------------------------------------------- + + eG0T0(:) = seHF(:) + Z(:)*SigT(:) + +!---------------------------------------------- +! Dump results +!---------------------------------------------- + + call print_G0T0(nBas2,nO2,seHF,ENuc,ERHF,SigT,Z,eG0T0,EcGM,EcRPA) + + + call Bethe_Salpeter_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nS2,nOO,nVV,Omega1,X1,Y1,Omega2,X2,Y2,rho1,rho2, & + sERI,seHF,eG0T0,EcBSE) + +end subroutine soG0T0 diff --git a/src/GW/Bethe_Salpeter.f90 b/src/GW/Bethe_Salpeter.f90 index a71949a..2440585 100644 --- a/src/GW/Bethe_Salpeter.f90 +++ b/src/GW/Bethe_Salpeter.f90 @@ -42,14 +42,17 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) + double precision,allocatable :: WA_sta(:,:) + double precision,allocatable :: WB_sta(:,:) + ! Output variables double precision,intent(out) :: EcBSE(nspin) ! Memory allocation - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - allocate(OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + WA_sta(nS,nS),WB_sta(nS,nS),OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) !--------------------------------- ! Compute (singlet) RPA screening @@ -58,10 +61,13 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, isp_W = 1 EcRPA = 0d0 - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA_sta) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB_sta) + !------------------- ! Singlet manifold !------------------- @@ -73,8 +79,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ! Compute BSE excitation energies - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & - rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) @@ -112,8 +118,8 @@ subroutine Bethe_Salpeter(TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC, ! Compute BSE excitation energies - call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,OmRPA, & - rho_RPA,EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + call linear_response_BSE(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI,WA_sta,WB_sta, & + EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GW ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) diff --git a/src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 b/src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 index 2d82750..69b026b 100644 --- a/src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 +++ b/src/GW/Bethe_Salpeter_AB_matrix_dynamic.f90 @@ -90,16 +90,16 @@ subroutine Bethe_Salpeter_AB_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,eGW,O eps_Am = - OmBSE - OmRPA(kc) - (eGW(b) - eGW(i)) chi_Am = chi_Am + rho_RPA(i,j,kc)*rho_RPA(a,b,kc)*eps_Am/(eps_Am**2 + eta**2) - eps_Bp = + OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) + eps_Bp = - OmRPA(kc) - (eGW(a) - eGW(b)) chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) - eps_Bp = + OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) + eps_Bp = - OmRPA(kc) - (eGW(j) - eGW(i)) chi_Bp = chi_Bp + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bp/(eps_Bp**2 + eta**2) - eps_Bm = - OmBSE - OmRPA(kc) - (eGW(a) - eGW(b)) + eps_Bm = - OmRPA(kc) - (eGW(a) - eGW(b)) chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) - eps_Bm = - OmBSE - OmRPA(kc) - (eGW(j) - eGW(i)) + eps_Bm = - OmRPA(kc) - (eGW(j) - eGW(i)) chi_Bm = chi_Bm + rho_RPA(i,b,kc)*rho_RPA(a,j,kc)*eps_Bm/(eps_Bm**2 + eta**2) enddo diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index aa35dd4..e6f49bc 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -109,8 +109,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & ! Compute screening ! !-------------------! - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0, & - eHF,ERI_MO,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0, & + eHF,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) @@ -168,8 +168,8 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, & ! Compute the RPA correlation energy - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI_MO,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eG0W0,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) !--------------! ! Dump results ! diff --git a/src/GW/G0W0_SOSEX.f90 b/src/GW/G0W0_SOSEX.f90 index fa1eeb5..696fa3b 100644 --- a/src/GW/G0W0_SOSEX.f90 +++ b/src/GW/G0W0_SOSEX.f90 @@ -93,8 +93,8 @@ subroutine G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDy !-------------------! do ispin=1,nspin - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO, & - OmRPA(:,ispin),rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO, & + EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) end do @@ -126,8 +126,8 @@ subroutine G0W0_SOSEX(doACFDT,exchange_kernel,doXBS,BSE,TDA_W,TDA,dBSE,dTDA,evDy ! Compute the RPA correlation energy do ispin=1,nspin - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eSOSEX,ERI_MO,OmRPA(:,ispin), & - rho_RPA(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eSOSEX,ERI_MO, & + EcRPA(ispin),OmRPA(:,ispin),XpY_RPA(:,:,ispin),XmY_RPA(:,:,ispin)) end do !--------------! diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index f510975..c691616 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -151,8 +151,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) end if diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 806ae6d..9f1d33b 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -181,8 +181,8 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, if(.not. GW0 .or. nSCF == 0) then - call linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & - OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,OmRPA) endif diff --git a/src/GW/self_energy_exchange_diag.f90 b/src/GW/self_energy_exchange_diag.f90 index 0e305b7..eb0cebc 100644 --- a/src/GW/self_energy_exchange_diag.f90 +++ b/src/GW/self_energy_exchange_diag.f90 @@ -14,7 +14,8 @@ subroutine self_energy_exchange_diag(nBas,c,P,ERI,SigX) ! Local variables - integer :: q,mu,nu + integer :: mu,nu + integer :: q double precision,allocatable :: Fx(:,:) ! Output variables diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index a66f2aa..0e3989e 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -1,4 +1,5 @@ -subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx) +subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,F,ERI,dipole_int,X,ERHF,e,c,P,Vx) ! Perform restricted Hartree-Fock calculation @@ -7,8 +8,11 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Input variables - integer,intent(in) :: maxSCF,max_diis,guess_type + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + integer,intent(in) :: guess_type double precision,intent(in) :: thresh + double precision,intent(in) :: level_shift integer,intent(in) :: nBas integer,intent(in) :: nO @@ -46,7 +50,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T double precision,allocatable :: K(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: Fp(:,:) - double precision,allocatable :: ON(:) ! Output variables @@ -71,30 +74,22 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Memory allocation - allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas), & - cp(nBas,nBas),Fp(nBas,nBas),ON(nBas), & + allocate(J(nBas,nBas),K(nBas,nBas),error(nBas,nBas),cp(nBas,nBas),Fp(nBas,nBas), & error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) -! Guess coefficients and eigenvalues +! Guess coefficients and density matrix - call mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) - -! ON(:) = 0d0 -! do i=1,nO -! ON(i) = 1d0 -! ON(i) = dble(2*i-1) -! end do - -! call density_matrix(nBas,ON,c,P) + call mo_guess(nBas,guess_type,S,Hc,X,c) + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) ! Initialization - n_diis = 0 F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - Conv = 1d0 - nSCF = 0 - rcond = 0d0 + Conv = 1d0 + n_diis = 0 + nSCF = 0 + rcond = 0d0 !------------------------------------------------------------------------ ! Main SCF loop @@ -123,17 +118,25 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Check convergence error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) - Conv = maxval(abs(error)) + Conv = maxval(abs(error)) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - if(abs(rcond) > 1d-7) then + if(abs(rcond) > 1d-15) then + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,error_diis,F_diis,error,F) + else + n_diis = 0 + end if +! Level-shifting + + if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F) + ! Diagonalize Fock matrix Fp = matmul(transpose(X),matmul(F,X)) @@ -144,14 +147,12 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Density matrix P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - -! call density_matrix(nBas,ON,c,P) ! Compute HF energy - ERHF = trace_matrix(nBas,matmul(P,Hc)) & - + 0.5d0*trace_matrix(nBas,matmul(P,J)) & - + 0.25d0*trace_matrix(nBas,matmul(P,K)) + ERHF = trace_matrix(nBas,matmul(P,Hc)) & + + 0.5d0*trace_matrix(nBas,matmul(P,J)) & + + 0.25d0*trace_matrix(nBas,matmul(P,K)) ! Compute HOMO-LUMO gap @@ -205,6 +206,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Compute Vx for post-HF calculations - call mo_fock_exchange_potential(nBas,c,K,Vx) + call mo_fock_exchange_potential(nBas,c,P,ERI,Vx) end subroutine RHF diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index afa1fad..39bbaf6 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -1,4 +1,5 @@ -subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx) +subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,nO,S,T,V,Hc,ERI,dipole_int,X,EUHF,e,c,P,Vx) ! Perform unrestricted Hartree-Fock calculation @@ -11,6 +12,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO integer,intent(in) :: max_diis integer,intent(in) :: guess_type logical,intent(in) :: mix + double precision,intent(in) :: level_shift double precision,intent(in) :: thresh integer,intent(in) :: nBas @@ -79,22 +81,12 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO K(nBas,nBas,nspin),err(nBas,nBas,nspin),cp(nBas,nBas,nspin), & err_diis(nBasSq,max_diis,nspin),F_diis(nBasSq,max_diis,nspin)) -! Guess coefficients and eigenvalues +! Guess coefficients and demsity matrices - if(guess_type == 1) then - - F(:,:,:) = 0d0 - do ispin=1,nspin - F(:,:,ispin) = Hc(:,:) - end do - - else if(guess_type == 2) then - - do ispin=1,nspin - call random_number(F(:,:,ispin)) - end do - - end if + do ispin=1,nspin + call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin)) + P(:,:,ispin) = matmul(c(:,1:nO(ispin),ispin),transpose(c(:,1:nO(ispin),ispin))) + end do ! Initialization @@ -179,13 +171,28 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - if(minval(rcond(:)) > 1d-7) then + + if(minval(rcond(:)) > 1d-15) then + do ispin=1,nspin if(nO(ispin) > 1) call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis,err_diis(:,1:n_diis,ispin), & F_diis(:,1:n_diis,ispin),err(:,:,ispin),F(:,:,ispin)) end do + else + n_diis = 0 + + end if + +! Level-shifting + + if(level_shift > 0d0 .and. Conv > thresh) then + + do ispin=1,nspin + call level_shifting(level_shift,nBas,nO(ispin),S,c(:,:,ispin),F(:,:,ispin)) + end do + end if !------------------------------------------------------------------------ @@ -253,7 +260,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO ! Compute Vx for post-HF calculations do ispin=1,nspin - call mo_fock_exchange_potential(nBas,c(:,:,ispin),K(:,:,ispin),Vx(:,ispin)) + call mo_fock_exchange_potential(nBas,c(:,:,ispin),P(:,:,ispin),ERI,Vx(:,ispin)) end do end subroutine UHF diff --git a/src/HF/core_guess.f90 b/src/HF/core_guess.f90 new file mode 100644 index 0000000..7e45e8e --- /dev/null +++ b/src/HF/core_guess.f90 @@ -0,0 +1,33 @@ +subroutine core_guess(nBas,Hc,X,c) + +! Core guess of the molecular orbitals for HF calculation + + implicit none + +! Input variables + + integer,intent(in) :: nBas + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + +! Local variables + + double precision,allocatable :: cp(:,:) + double precision,allocatable :: e(:) + + +! Output variables + + double precision,intent(out) :: c(nBas,nBas) + +! Memory allocation + + allocate(cp(nBas,nBas),e(nBas)) + +! Core guess + + cp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) + call diagonalize_matrix(nBas,cp,e) + c(:,:) = matmul(X(:,:),cp(:,:)) + +end subroutine diff --git a/src/HF/huckel_guess.f90 b/src/HF/huckel_guess.f90 index 91da7cd..80efaa8 100644 --- a/src/HF/huckel_guess.f90 +++ b/src/HF/huckel_guess.f90 @@ -1,47 +1,39 @@ -subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) +subroutine huckel_guess(nBas,S,Hc,X,c) -! Hickel guess of the molecular orbitals for HF calculation +! Hickel guess implicit none ! Input variables integer,intent(in) :: nBas - integer,intent(in) :: nO double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(inout):: J(nBas,nBas) - double precision,intent(inout):: K(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) - double precision,intent(inout):: cp(nBas,nBas) - double precision,intent(inout):: F(nBas,nBas) - double precision,intent(inout):: Fp(nBas,nBas) - double precision,intent(inout):: e(nBas) - double precision,intent(inout):: P(nBas,nBas) ! Local variables integer :: mu,nu double precision :: a + double precision,allocatable :: F(:,:) + ! Output variables double precision,intent(out) :: c(nBas,nBas) +! Memory allocation + + allocate(F(nBas,nBas)) + +! Extended Huckel parameter + a = 1.75d0 - Fp(:,:) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) - c(:,:) = matmul(X(:,:),cp(:,:)) - - call Coulomb_matrix_AO_basis(nBas,P,ERI,J) - call exchange_matrix_AO_basis(nBas,P,ERI,K) - - F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) +! GWH approximation do mu=1,nBas + F(mu,mu) = Hc(mu,mu) do nu=mu+1,nBas F(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu)) @@ -50,9 +42,6 @@ subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) enddo enddo - Fp(:,:) = matmul(transpose(X(:,:)),matmul(F(:,:),X(:,:))) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) - c(:,:) = matmul(X(:,:),cp(:,:)) + call core_guess(nBas,F,X,c) end subroutine diff --git a/src/HF/mo_fock_exchange_potential.f90 b/src/HF/mo_fock_exchange_potential.f90 index 56203c3..d07fcce 100644 --- a/src/HF/mo_fock_exchange_potential.f90 +++ b/src/HF/mo_fock_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx) +subroutine mo_fock_exchange_potential(nBas,c,P,ERI,Vx) ! Compute the exchange potential in the MO basis @@ -9,12 +9,14 @@ subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx) integer,intent(in) :: nBas double precision,intent(in) :: c(nBas,nBas) - double precision,intent(in) :: Fx(nBas,nBas) + double precision,intent(in) :: P(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) ! Local variables integer :: mu,nu - integer :: p + integer :: q + double precision,allocatable :: Fx(:,:) ! Output variables @@ -22,13 +24,18 @@ subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx) ! Compute Vx + allocate(Fx(nBas,nBas)) + call exchange_matrix_AO_basis(nBas,P,ERI,Fx) + Vx(:) = 0d0 - do p=1,nBas + do q=1,nBas do mu=1,nBas do nu=1,nBas - Vx(p) = Vx(p) + c(mu,p)*Fx(mu,nu)*c(nu,p) + Vx(q) = Vx(q) + c(mu,q)*Fx(mu,nu)*c(nu,q) end do end do end do + deallocate(Fx) + end subroutine mo_fock_exchange_potential diff --git a/src/HF/mo_guess.f90 b/src/HF/mo_guess.f90 index d478ab6..387f371 100644 --- a/src/HF/mo_guess.f90 +++ b/src/HF/mo_guess.f90 @@ -1,4 +1,4 @@ -subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) +subroutine mo_guess(nBas,guess_type,S,Hc,X,c) ! Guess of the molecular orbitals for HF calculation @@ -7,19 +7,10 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) ! Input variables integer,intent(in) :: nBas - integer,intent(in) :: nO integer,intent(in) :: guess_type double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(inout):: J(nBas,nBas) - double precision,intent(inout):: K(nBas,nBas) double precision,intent(in) :: X(nBas,nBas) - double precision,intent(inout):: cp(nBas,nBas) - double precision,intent(inout):: F(nBas,nBas) - double precision,intent(inout):: Fp(nBas,nBas) - double precision,intent(inout):: e(nBas) - double precision,intent(inout):: P(nBas,nBas) ! Local variables @@ -31,14 +22,11 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) if(guess_type == 1) then - Fp = matmul(transpose(X),matmul(Hc,X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,e) - c = matmul(X,cp) + call core_guess(nBas,Hc,X,c) elseif(guess_type == 2) then - call huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) + call huckel_guess(nBas,S,Hc,X,c) elseif(guess_type == 3) then @@ -51,6 +39,4 @@ subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,F,Fp,e,c,P) endif - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - end subroutine diff --git a/src/LR/linear_response.f90 b/src/LR/linear_response.f90 index ee80c72..e239f54 100644 --- a/src/LR/linear_response.f90 +++ b/src/LR/linear_response.f90 @@ -1,4 +1,4 @@ -subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Omega_RPA,rho_RPA,EcRPA,Omega,XpY,XmY) +subroutine linear_response(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,Ec,Omega,XpY,XmY) ! Compute linear response @@ -9,7 +9,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E logical,intent(in) :: dRPA logical,intent(in) :: TDA - logical,intent(in) :: BSE double precision,intent(in) :: eta integer,intent(in) :: ispin integer,intent(in) :: nBas @@ -22,9 +21,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: Omega_RPA(nS) - double precision,intent(in) :: rho_RPA(nBas,nBas,nS) - ! Local variables double precision :: trace_matrix @@ -38,7 +34,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E ! Output variables - double precision,intent(out) :: EcRPA + double precision,intent(out) :: Ec double precision,intent(out) :: Omega(nS) double precision,intent(out) :: XpY(nS,nS) double precision,intent(out) :: XmY(nS,nS) @@ -51,8 +47,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) - if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,A) - ! Tamm-Dancoff approximation if(TDA) then @@ -67,8 +61,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) - if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega_RPA,rho_RPA,B) - ! Build A + B and A - B matrices ApB = A + B @@ -111,6 +103,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E ! Compute the RPA correlation energy - EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) + Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) end subroutine linear_response diff --git a/src/LR/linear_response_Tmatrix.f90 b/src/LR/linear_response_BSE.f90 similarity index 75% rename from src/LR/linear_response_Tmatrix.f90 rename to src/LR/linear_response_BSE.f90 index 89b9207..cd02e93 100644 --- a/src/LR/linear_response_Tmatrix.f90 +++ b/src/LR/linear_response_BSE.f90 @@ -1,4 +1,4 @@ -subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_BSE,B_BSE,EcRPA,Omega,XpY,XmY) +subroutine linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A_BSE,B_BSE,Ec,Omega,XpY,XmY) ! Compute linear response @@ -7,9 +7,17 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! Input variables - logical,intent(in) :: dRPA,TDA + integer,intent(in) :: ispin + logical,intent(in) :: dRPA + logical,intent(in) :: TDA + logical,intent(in) :: BSE double precision,intent(in) :: eta - integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS double precision,intent(in) :: lambda double precision,intent(in) :: e(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -29,7 +37,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! Output variables - double precision,intent(out) :: EcRPA + double precision,intent(out) :: Ec double precision,intent(out) :: Omega(nS) double precision,intent(out) :: XpY(nS,nS) double precision,intent(out) :: XmY(nS,nS) @@ -42,12 +50,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) -! print*,'A' -! call matout(nS,nS,A) -! print*,'TA' -! call matout(nS,nS,A_BSE) - - A(:,:) = A(:,:) - A_BSE(:,:) + if(BSE) A(:,:) = A(:,:) - A_BSE(:,:) ! Tamm-Dancoff approximation @@ -63,12 +66,7 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) -! print*,'B' -! call matout(nS,nS,B) -! print*,'TB' -! call matout(nS,nS,B_BSE) - - B(:,:) = B(:,:) - B_BSE(:,:) + if(BSE) B(:,:) = B(:,:) - B_BSE(:,:) ! Build A + B and A - B matrices @@ -82,10 +80,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda if(minval(Omega) < 0d0) & call print_warning('You may have instabilities in linear response: A-B is not positive definite!!') -! do ia=1,nS -! if(Omega(ia) < 0d0) Omega(ia) = 0d0 -! end do - call ADAt(nS,AmB,1d0*sqrt(Omega),AmBSq) call ADAt(nS,AmB,1d0/sqrt(Omega),AmBIv) @@ -95,10 +89,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda if(minval(Omega) < 0d0) & call print_warning('You may have instabilities in linear response: negative excitations!!') - - ! do ia=1,nS - ! if(Omega(ia) < 0d0) Omega(ia) = 0d0 - ! end do Omega = sqrt(Omega) @@ -112,6 +102,6 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda ! Compute the RPA correlation energy - EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) + Ec = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) -end subroutine linear_response_Tmatrix +end subroutine linear_response_BSE diff --git a/src/LR/unrestricted_linear_response_C_pp.f90 b/src/LR/unrestricted_linear_response_C_pp.f90 index 45e8643..7e37f78 100644 --- a/src/LR/unrestricted_linear_response_C_pp.f90 +++ b/src/LR/unrestricted_linear_response_C_pp.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,& - e,ERI_aaaa,ERI_aabb,ERI_bbbb,C_pp) +subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,& + lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C_pp) ! Compute linear response @@ -54,7 +54,7 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do d=nO(2)+1,nBas-nR(2) cd = cd + 1 C_pp(ab,cd) = (e(a,1) + e(b,2))*Kronecker_delta(a,c) & -*Kronecker_delta(b,d) + lambda*ERI_aabb(a,b,c,d) + *Kronecker_delta(b,d) + lambda*ERI_aabb(a,b,c,d) end do end do end do @@ -79,9 +79,10 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP do d=c+1,nBas-nR(1) cd = cd + 1 - C_pp(ab,cd) = (e(a,1) + e(b,1) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & - + lambda*(ERI_aaaa(a,b,c,d) - ERI_aaaa(a,b,d,c)) -!write(*,*) C_pp(ab,cd) + C_pp(ab,cd) = (e(a,1) + e(b,1) - eF)*Kronecker_delta(a,c)& + *Kronecker_delta(b,d) + lambda*(ERI_aaaa(a,b,c,d) & + - ERI_aaaa(a,b,d,c)) + end do end do end do @@ -103,12 +104,13 @@ subroutine unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nP cd = cd + 1 C_pp(ab,cd) = (e(a,2) + e(b,2) - eF)*Kronecker_delta(a,c) & - *Kronecker_delta(b,d) + lambda*(ERI_bbbb(a,b,c,d) - ERI_bbbb(a,b,d,c)) + *Kronecker_delta(b,d) + lambda*(ERI_bbbb(a,b,c,d) & + - ERI_bbbb(a,b,d,c)) + end do end do - end do end do - end do + end do end if diff --git a/src/LR/unrestricted_linear_response_D_pp.f90 b/src/LR/unrestricted_linear_response_D_pp.f90 index 02ab1d0..3d7c10e 100644 --- a/src/LR/unrestricted_linear_response_D_pp.f90 +++ b/src/LR/unrestricted_linear_response_D_pp.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,& - e,ERI_aaaa,ERI_aabb,ERI_bbbb,D_pp) +subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,& + lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,D_pp) ! Compute linear response @@ -55,7 +55,7 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH do l=nC(2)+1,nO(2) kl = kl + 1 D_pp(ij,kl) = -(e(i,1) + e(j,2))*Kronecker_delta(i,k)& - *Kronecker_delta(j,l) +lambda*ERI_aabb(i,j,k,l) + *Kronecker_delta(j,l) +lambda*ERI_aabb(i,j,k,l) end do end do end do @@ -81,14 +81,15 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH do l=k+1,nO(1) kl = kl + 1 - D_pp(ij,kl) = -(e(i,1) + e(j,1) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + lambda*(ERI_aaaa(i,j,k,l) - ERI_aaaa(i,j,l,k)) + D_pp(ij,kl) = -(e(i,1) + e(j,1) - eF)*Kronecker_delta(i,k)& + *Kronecker_delta(j,l) + lambda*(ERI_aaaa(i,j,k,l) & + - ERI_aaaa(i,j,l,k)) end do end do end do end do - end if + end if if (ispin == 3) then @@ -104,7 +105,8 @@ subroutine unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nH kl = kl + 1 D_pp(ij,kl) = -(e(i,2) + e(j,2) - eF)*Kronecker_delta(i,k) & - *Kronecker_delta(j,l) + lambda*(ERI_bbbb(i,j,k,l) - ERI_bbbb(i,j,l,k)) + *Kronecker_delta(j,l) + lambda*(ERI_bbbb(i,j,k,l) & + - ERI_bbbb(i,j,l,k)) end do end do diff --git a/src/LR/unrestricted_linear_response_pp.f90 b/src/LR/unrestricted_linear_response_pp.f90 index 3111e89..89acb7f 100644 --- a/src/LR/unrestricted_linear_response_pp.f90 +++ b/src/LR/unrestricted_linear_response_pp.f90 @@ -1,6 +1,7 @@ -subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt, & -nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1,X1,Y1,Omega2,X2,Y2,& -EcRPA) +subroutine unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,& + nPt,nHaa,nHab,nHbb,nHt,lambda,e,ERI_aaaa,& + ERI_aabb,ERI_bbbb,Omega1,X1,Y1,Omega2,X2,Y2,& + EcRPA) ! Compute linear response for unrestricted formalism @@ -56,53 +57,47 @@ EcRPA) ! Memory allocation - allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt)& -,Omega(nPt+nHt)) -!write(*,*) 'Hello' + allocate(C(nPt,nPt),B(nPt,nHt),D(nHt,nHt),M(nPt+nHt,nPt+nHt),Z(nPt+nHt,nPt+nHt),& + Omega(nPt+nHt)) + ! Build C, B and D matrices for the pp channel - call unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,lambda,& -e,ERI_aaaa,ERI_aabb,ERI_bbbb,C) + call unrestricted_linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,& + lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,C) - call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,& -nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,ERI_bbbb,B) + call unrestricted_linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nPt,nHaa,& + nHab,nHbb,nHt,lambda,ERI_aaaa,ERI_aabb,& + ERI_bbbb,B) -!call matout(nPt,nHt,B) -!write(*,*) 'Hello' -call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,lambda,& -e,ERI_aaaa,ERI_aabb,ERI_bbbb,D) + call unrestricted_linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nHaa,nHab,nHbb,nHt,& + lambda,e,ERI_aaaa,ERI_aabb,ERI_bbbb,D) -!call matout(nHt,nHt,D) -!write(*,*) 'Hello' ! Diagonal blocks - M( 1:nPt , 1:nPt) = + C(1:nPt,1:nPt) - M(nPt+1:nPt+nHt,nPt+1:nPt+nHt) = - D(1:nHt,1:nHt) + M( 1:nPt , 1:nPt) = + C(1:nPt,1:nPt) + M(nPt+1:nPt+nHt,nPt+1:nPt+nHt) = - D(1:nHt,1:nHt) - ! Off-diagonal blocks +! Off-diagonal blocks - M( 1:nPt ,nPt+1:nHt+nPt) = - B(1:nPt,1:nHt) - M(nPt+1:nHt+nPt, 1:nPt) = + transpose(B(1:nPt,1:nHt)) - -!call matout(nPt+nHt,nPt+nHt,M) + M( 1:nPt ,nPt+1:nHt+nPt) = - B(1:nPt,1:nHt) + M(nPt+1:nHt+nPt, 1:nPt) = + transpose(B(1:nPt,1:nHt)) ! Diagonalize the p-h matrix if(nHt+nPt > 0) call diagonalize_general_matrix(nHt+nPt,M,Omega,Z) - ! Split the various quantities in p-p and h-h parts +! Split the various quantities in p-p and h-h parts - call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),& -Y2(:,:)) + call sort_ppRPA(nHt,nPt,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),& + Y2(:,:)) ! Compute the RPA correlation energy - EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) - trace_matrix(nHt,D(:,:)) ) + EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nPt,C(:,:)) & + - trace_matrix(nHt,D(:,:)) ) EcRPA1 = +sum(Omega1(:)) - trace_matrix(nPt,C(:,:)) EcRPA2 = -sum(Omega2(:)) - trace_matrix(nHt,D(:,:)) if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' - - end subroutine unrestricted_linear_response_pp diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 2959cec..40d1415 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -111,7 +111,7 @@ program QuAcK double precision :: start_Bas ,end_Bas ,t_Bas integer :: maxSCF_HF,n_diis_HF - double precision :: thresh_HF + double precision :: thresh_HF,level_shift logical :: DIIS_HF,guess_type,ortho_type,mix integer :: maxSCF_CC,n_diis_CC @@ -180,15 +180,15 @@ program QuAcK ! Read options for methods - call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, & - maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - TDA,singlet,triplet,spin_conserved,spin_flip, & - maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, & - COHSEX,SOSEX,TDA_W,G0W,GW0, & - maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & - doACFDT,exchange_kernel,doXBS, & - BSE,dBSE,dTDA,evDyn, & + call read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & + maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & + TDA,singlet,triplet,spin_conserved,spin_flip, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, & + COHSEX,SOSEX,TDA_W,G0W,GW0, & + maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & + doACFDT,exchange_kernel,doXBS, & + BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Weird stuff @@ -292,7 +292,7 @@ program QuAcK end if call cpu_time(start_HF) - call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nNuc,ZNuc,rNuc,ENuc, & + call RHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc, & nBas,nO,S,T,V,Hc,F_AO,ERI_AO,dipole_int_AO,X,ERHF,eHF,cHF,PHF,Vxc) call cpu_time(end_HF) @@ -312,7 +312,7 @@ program QuAcK unrestricted = .true. call cpu_time(start_HF) - call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc, & + call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & nBas,nO,S,T,V,Hc,ERI_AO,dipole_int_AO,X,EUHF,eHF,cHF,PHF,Vxc) call cpu_time(end_HF) @@ -332,7 +332,7 @@ program QuAcK unrestricted = .true. call cpu_time(start_KS) - call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, & + call eDFT(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC, & nO,nV,nR,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI_AO,dipole_int_AO,EUHF,eHF,cHF,PHF,Vxc) @@ -1159,6 +1159,7 @@ program QuAcK else +! call soG0T0(eta_GT,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI_MO,eHF) call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet, & linGT,eta_GT,regGT,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int_MO, & PHF,cHF,eHF,Vxc,eG0T0) diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 383aaee..8a18c34 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -1,12 +1,12 @@ -subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,dostab, & - maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & - TDA,singlet,triplet,spin_conserved,spin_flip, & - maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & - maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, & - COHSEX,SOSEX,TDA_W,G0W,GW0, & - maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & - doACFDT,exchange_kernel,doXBS, & - BSE,dBSE,dTDA,evDyn, & +subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_type,mix,level_shift,dostab, & + maxSCF_CC,thresh_CC,DIIS_CC,n_diis_CC, & + TDA,singlet,triplet,spin_conserved,spin_flip, & + maxSCF_GF,thresh_GF,DIIS_GF,n_diis_GF,linGF,eta_GF,renormGF,regGF, & + maxSCF_GW,thresh_GW,DIIS_GW,n_diis_GW,linGW,eta_GW,regGW, & + COHSEX,SOSEX,TDA_W,G0W,GW0, & + maxSCF_GT,thresh_GT,DIIS_GT,n_diis_GT,linGT,eta_GT,regGT,TDA_T, & + doACFDT,exchange_kernel,doXBS, & + BSE,dBSE,dTDA,evDyn, & nMC,nEq,nWalk,dt,nPrint,iSeed,doDrift) ! Read desired methods @@ -22,6 +22,7 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t integer,intent(out) :: guess_type integer,intent(out) :: ortho_type logical,intent(out) :: mix + double precision,intent(out) :: level_shift logical,intent(out) :: dostab integer,intent(out) :: maxSCF_CC @@ -100,14 +101,15 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t guess_type = 1 ortho_type = 1 mix = .false. + level_shift = 0d0 dostab = .false. read(1,*) - read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,answer3 + read(1,*) maxSCF_HF,thresh_HF,answer1,n_diis_HF,guess_type,ortho_type,answer2,level_shift,answer3 - if(answer1 == 'T') DIIS_HF = .true. - if(answer2 == 'T') mix = .true. - if(answer3 == 'T') dostab = .true. + if(answer1 == 'T') DIIS_HF = .true. + if(answer2 == 'T') mix = .true. + if(answer3 == 'T') dostab = .true. if(.not.DIIS_HF) n_diis_HF = 1 diff --git a/src/RPA/ACFDT.f90 b/src/RPA/ACFDT.f90 index 24d9de0..4a7f33a 100644 --- a/src/RPA/ACFDT.f90 +++ b/src/RPA/ACFDT.f90 @@ -32,6 +32,8 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB double precision,allocatable :: Ec(:,:) double precision :: EcRPA + double precision,allocatable :: WA(:,:) + double precision,allocatable :: WB(:,:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) @@ -48,7 +50,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + allocate(WA(nS,nS),WB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) ! Antisymmetrized kernel version @@ -69,10 +71,13 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB isp_W = 1 EcRPA = 0d0 - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) + ! Singlet manifold if(singlet) then @@ -94,15 +99,18 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) ! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -143,14 +151,17 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 index 9b25e29..051deab 100644 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ b/src/RPA/ACFDT_Tmatrix.f90 @@ -1,5 +1,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - ERI,eT,eGT,EcAC) + nOOs,nVVs,nOOt,nVVt,Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,rho1s,rho2s,Omega1t,X1t,Y1t, & + Omega2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC) ! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix @@ -26,6 +27,28 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple integer,intent(in) :: nR integer,intent(in) :: nS + integer,intent(in) :: nOOs + integer,intent(in) :: nOOt + integer,intent(in) :: nVVs + integer,intent(in) :: nVVt + + double precision,intent(in) :: Omega1s(nVVs) + double precision,intent(in) :: X1s(nVVs,nVVs) + double precision,intent(in) :: Y1s(nOOs,nVVs) + double precision,intent(in) :: Omega2s(nOOs) + double precision,intent(in) :: X2s(nVVs,nOOs) + double precision,intent(in) :: Y2s(nOOs,nOOs) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs) + double precision,intent(in) :: Omega1t(nVVt) + double precision,intent(in) :: X1t(nVVt,nVVt) + double precision,intent(in) :: Y1t(nOOt,nVVt) + double precision,intent(in) :: Omega2t(nOOt) + double precision,intent(in) :: X2t(nVVt,nOOt) + double precision,intent(in) :: Y2t(nOOt,nOOt) + double precision,intent(in) :: rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: rho2t(nBas,nBas,nOOt) + double precision,intent(in) :: eT(nBas) double precision,intent(in) :: eGT(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -39,46 +62,23 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple double precision :: lambda double precision,allocatable :: Ec(:,:) - integer :: nOOs,nOOt - integer :: nVVs,nVVt - double precision :: EcRPA(nspin) - double precision,allocatable :: TA(:,:) - double precision,allocatable :: TB(:,:) + double precision,allocatable :: TAs(:,:) + double precision,allocatable :: TBs(:,:) + double precision,allocatable :: TAt(:,:) + double precision,allocatable :: TBt(:,:) double precision,allocatable :: Omega(:,:) double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: Omega1s(:),Omega1t(:) - double precision,allocatable :: X1s(:,:),X1t(:,:) - double precision,allocatable :: Y1s(:,:),Y1t(:,:) - double precision,allocatable :: rho1s(:,:,:),rho1t(:,:,:) - double precision,allocatable :: Omega2s(:),Omega2t(:) - double precision,allocatable :: X2s(:,:),X2t(:,:) - double precision,allocatable :: Y2s(:,:),Y2t(:,:) - double precision,allocatable :: rho2s(:,:,:),rho2t(:,:,:) - ! Output variables double precision,intent(out) :: EcAC(nspin) -! Useful quantities - - nOOs = nO*nO - nVVs = nV*nV - - nOOt = nO*(nO-1)/2 - nVVt = nV*(nV-1)/2 - ! Memory allocation - allocate(Omega1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & - Omega2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & - Omega1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & - Omega2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & - rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt)) - allocate(TA(nS,nS),TB(nS,nS),Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) + allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), & + Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) allocate(Ec(nAC,nspin)) ! Antisymmetrized kernel version @@ -113,11 +113,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple lambda = rAC(iAC) - ! Initialize T matrix - - TA(:,:) = 0d0 - TB(:,:) = 0d0 - if(doXBS) then isp_T = 1 @@ -128,8 +123,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs) isp_T = 2 iblock = 4 @@ -139,13 +134,13 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt) end if - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -183,11 +178,6 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple lambda = rAC(iAC) - ! Initialize T matrix - - TA(:,:) = 0d0 - TB(:,:) = 0d0 - if(doXBS) then isp_T = 1 @@ -198,8 +188,8 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TAs) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Omega1s,rho1s,Omega2s,rho2s,TBs) isp_T = 2 iblock = 4 @@ -209,13 +199,13 @@ subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triple call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TAt) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Omega1t,rho1t,Omega2t,rho2t,TBt) end if - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TA,TB, & - EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) diff --git a/src/RPA/ACFDT_cr.f90 b/src/RPA/ACFDT_cr.f90 index d82c77e..7d128a6 100644 --- a/src/RPA/ACFDT_cr.f90 +++ b/src/RPA/ACFDT_cr.f90 @@ -33,6 +33,8 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta double precision,allocatable :: Ec(:,:) double precision :: EcRPA + double precision,allocatable :: WA(:,:) + double precision,allocatable :: WB(:,:) double precision,allocatable :: OmRPA(:) double precision,allocatable :: XpY_RPA(:,:) double precision,allocatable :: XmY_RPA(:,:) @@ -49,7 +51,7 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta ! Memory allocation allocate(Ec(nAC,nspin)) - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) + allocate(WA(nS,nS),WB(nS,nS),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) allocate(Omega(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) ! Antisymmetrized kernel version @@ -70,10 +72,13 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta isp_W = 1 EcRPA = 0d0 - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,1d0,ERI,OmRPA,rho_RPA,WB) + ! Singlet manifold if(singlet) then @@ -95,15 +100,18 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) ! call print_excitation('W^lambda: ',isp_W,nS,OmRPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS, & ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) @@ -144,14 +152,17 @@ subroutine ACFDT_cr(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta if(doXBS) then - call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & - rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call linear_response(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI, & + EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) + call static_screening_WA(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WA) + call static_screening_WB(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,OmRPA,rho_RPA,WB) + end if - call linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,OmRPA, & - rho_RPA,EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call linear_response_BSE(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,WA,WB, & + EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) diff --git a/src/RPA/RPA.f90 b/src/RPA/RPA.f90 index 0dce1a7..e02423a 100644 --- a/src/RPA/RPA.f90 +++ b/src/RPA/RPA.f90 @@ -33,7 +33,6 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision :: rho double precision :: EcRPA(nspin) double precision :: EcAC(nspin) @@ -67,7 +66,7 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, ispin = 1 - call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & + call linear_response(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -80,7 +79,7 @@ subroutine RPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR, ispin = 2 - call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,rho,Omega(:,ispin), & + call linear_response(ispin,.true.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/RPA/RPAx.f90 b/src/RPA/RPAx.f90 index f47351d..422d6ec 100644 --- a/src/RPA/RPAx.f90 +++ b/src/RPA/RPAx.f90 @@ -34,7 +34,6 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision :: rho double precision :: EcRPAx(nspin) double precision :: EcAC(nspin) @@ -69,7 +68,7 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR ispin = 1 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -82,7 +81,7 @@ subroutine RPAx(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,nR ispin = 2 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('RPAx@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/RPA/crRPA.f90 b/src/RPA/crRPA.f90 index 7671595..3312092 100644 --- a/src/RPA/crRPA.f90 +++ b/src/RPA/crRPA.f90 @@ -34,7 +34,6 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n double precision,allocatable :: XpY(:,:,:) double precision,allocatable :: XmY(:,:,:) - double precision :: rho double precision :: EcRPAx(nspin) double precision :: EcAC(nspin) @@ -68,7 +67,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 1 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) @@ -81,7 +80,7 @@ subroutine crRPA(TDA,doACFDT,exchange_kernel,singlet,triplet,eta,nBas,nC,nO,nV,n ispin = 2 - call linear_response(ispin,.false.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI,Omega(:,ispin),rho, & + call linear_response(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,-1d0,eHF,ERI, & EcRPAx(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call print_excitation('crRPA@HF ',ispin,nS,Omega(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) diff --git a/src/RPA/ppURPA.f90 b/src/RPA/ppURPA.f90 index 73f960f..75b0dd4 100644 --- a/src/RPA/ppURPA.f90 +++ b/src/RPA/ppURPA.f90 @@ -25,7 +25,7 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH ! Local variables - integer :: ispin + integer :: ispin,iblock integer :: nPaa,nPbb,nPab,nP_sc,nP_sf integer :: nHaa,nHbb,nHab,nH_sc,nH_sf double precision,allocatable :: Omega1sc(:),Omega1sf(:) @@ -56,25 +56,28 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH if(spin_conserved) then ispin = 1 + iblock = 1 -!spin-conserved quantities +!Spin-conserved quantities - nPab = nV(1)*nV(2) - nHab = nO(1)*nO(2) + nPab = nV(1)*nV(2) + nHab = nO(1)*nO(2) - nP_sc = nPab - nH_sc = nHab + nP_sc = nPab + nH_sc = nHab ! Memory allocation - allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), & - Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc)) + allocate(Omega1sc(nP_sc),X1sc(nP_sc,nP_sc),Y1sc(nH_sc,nP_sc), & + Omega2sc(nH_sc),X2sc(nP_sc,nH_sc),Y2sc(nH_sc,nH_sc)) - call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sc,nHaa,nHab,nHbb,nH_sc,1d0, & - e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc,Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin)) - - call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc) - call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc) + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nP_sc,nHaa,nHab,nHbb,nH_sc,1d0,e,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Omega1sc,X1sc,Y1sc, & + Omega2sc,X2sc,Y2sc,Ec_ppURPA(ispin)) + + call print_excitation('pp-RPA (N+2)',5,nP_sc,Omega1sc) + call print_excitation('pp-RPA (N-2)',5,nH_sc,Omega2sc) endif @@ -83,40 +86,45 @@ subroutine ppURPA(TDA,doACFDT,spin_conserved,spin_flip,nBas,nC,nO,nV,nR,ENuc,EUH if(spin_flip) then ispin = 2 + iblock = 2 -!spin-flip quantities +!Spin-flip quantities - nPaa = nV(1)*(nV(1)-1)/2 - nPbb = nV(2)*(nV(2)-1)/2 + nPaa = nV(1)*(nV(1)-1)/2 + nPbb = nV(2)*(nV(2)-1)/2 - nP_sf = nPaa + nP_sf = nPaa - nHaa = nO(1)*(nO(1)-1)/2 - nHbb = nO(2)*(nO(2)-1)/2 + nHaa = nO(1)*(nO(1)-1)/2 + nHbb = nO(2)*(nO(2)-1)/2 - nH_sf = nHaa + nH_sf = nHaa -allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & - Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) + allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & + Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) -call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sf, & -nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf,& -Ec_ppURPA(ispin)) + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb, & + nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa, & + ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf, & + Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) - ispin = 3 + deallocate(Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf) - nP_sf = nPbb - nH_sf = nHbb + iblock = 3 -!allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & -! Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) + nP_sf = nPbb + nH_sf = nHbb -call unrestricted_linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,nP_sf, & -nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,Omega2sf,X2sf,Y2sf,& -Ec_ppURPA(ispin)) + allocate(Omega1sf(nP_sf),X1sf(nP_sf,nP_sf),Y1sf(nH_sf,nP_sf), & + Omega2sf(nH_sf),X2sf(nP_sf,nH_sf),Y2sf(nH_sf,nH_sf)) - call print_excitation('pp-RPA (N+2)',6,nP_sf,Omega1sf) - call print_excitation('pp-RPA (N-2)',6,nH_sf,Omega2sf) + call unrestricted_linear_response_pp(iblock,TDA,nBas,nC,nO,nV,nR,nPaa,nPab,nPbb,& + nP_sf,nHaa,nHab,nHbb,nH_sf,1d0,e,ERI_aaaa,& + ERI_aabb,ERI_bbbb,Omega1sf,X1sf,Y1sf,& + Omega2sf,X2sf,Y2sf,Ec_ppURPA(ispin)) + + call print_excitation('pp-RPA (N+2)',6,nP_sf,Omega1sf) + call print_excitation('pp-RPA (N-2)',6,nH_sf,Omega2sf) endif diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/UKS.f90 similarity index 90% rename from src/eDFT/eDFT_UKS.f90 rename to src/eDFT/UKS.f90 index 9ed88dc..afa4914 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/UKS.f90 @@ -1,5 +1,5 @@ -subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & - nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,c,Pw,Vxc) +subroutine UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix,level_shift, & + nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,c,Pw,Vxc) ! Perform unrestricted Kohn-Sham calculation for ensembles @@ -20,6 +20,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max integer,intent(in) :: max_diis integer,intent(in) :: guess_type logical,intent(in) :: mix + double precision,intent(in) :: level_shift double precision,intent(in) :: thresh integer,intent(in) :: nBas double precision,intent(in) :: AO(nBas,nGrid) @@ -45,10 +46,11 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max integer :: xc_rung logical :: LDA_centered = .false. - integer :: nSCF,nBasSq + integer :: nSCF + integer :: nBasSq integer :: n_diis - integer :: nO(nspin) - double precision :: conv + integer :: nO(nspin,nEns) + double precision :: Conv double precision :: rcond(nspin) double precision :: ET(nspin) double precision :: EV(nspin) @@ -117,34 +119,22 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Guess coefficients and eigenvalues - nO(:) = 0 - do ispin=1,nspin - nO(ispin) = int(sum(occnum(:,ispin,1))) + nO(:,:) = 0 + do iEns=1,nEns + do ispin=1,nspin + nO(ispin,iEns) = int(sum(occnum(:,ispin,iEns))) + end do end do - - if(guess_type == 1) then - do ispin=1,nspin - cp(:,:,ispin) = matmul(transpose(X(:,:)),matmul(Hc(:,:),X(:,:))) - call diagonalize_matrix(nBas,cp(:,:,ispin),eKS(:,ispin)) - c(:,:,ispin) = matmul(X(:,:),cp(:,:,ispin)) - end do + do ispin=1,nspin + call mo_guess(nBas,guess_type,S,Hc,X,c(:,:,ispin)) + end do - ! Mix guess to enforce symmetry breaking - - if(mix) call mix_guess(nBas,nO,c) - - else if(guess_type == 2) then - - do ispin=1,nspin - call random_number(F(:,:,ispin)) - end do - - else - - print*,'Wrong guess option' - stop +! Mix guess for UHF solution in singlet states + if(mix) then + write(*,*) '!! guess mixing disabled in UKS !!' + write(*,*) end if ! Initialization @@ -175,7 +165,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max '|','#','|','E(KS)','|','Ex(KS)','|','Ec(KS)','|','Conv','|','nEl','|' write(*,*)'------------------------------------------------------------------------------------------' - do while(conv > thresh .and. nSCF < maxSCF) + do while(Conv > thresh .and. nSCF < maxSCF) ! Increment @@ -264,19 +254,34 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max err(:,:,ispin) = matmul(F(:,:,ispin),matmul(Pw(:,:,ispin),S(:,:))) - matmul(matmul(S(:,:),Pw(:,:,ispin)),F(:,:,ispin)) end do - if(nSCF > 1) conv = maxval(abs(err(:,:,:))) + if(nSCF > 1) Conv = maxval(abs(err(:,:,:))) ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - if(minval(rcond(:)) > 1d-15) then - do ispin=1,nspin + do ispin=1,nspin + + if(rcond(ispin) > 1d-15) then + call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, & err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) + else + + n_diis = 0 + + end if + + end do + +! Level-shifting + + if(level_shift > 0d0 .and. Conv > thresh) then + + do ispin=1,nspin + call level_shifting(level_shift,nBas,maxval(nO(ispin,:)),S,c,F(:,:,ispin)) end do - else - n_diis = 0 - end if + + end if ! Transform Fock matrix in orthogonal basis @@ -342,7 +347,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Dump results write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',conv,'|',sum(nEl(:)),'|' + '|',nSCF,'|',Ew + ENuc,'|',sum(Ex(:)),'|',sum(Ec(:)),'|',Conv,'|',sum(nEl(:)),'|' end do write(*,*)'------------------------------------------------------------------------------------------' @@ -384,4 +389,4 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max call individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & AO,dAO,T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew) -end subroutine eDFT_UKS +end subroutine UKS diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 0ce947b..ffc0632 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -1,4 +1,4 @@ -subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, & +subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc,nBas,nEl,nC,nO,nV,nR, & nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell, & max_ang_mom,min_exponent,max_exponent,S,T,V,Hc,X,ERI,dipole_int,Ew,eKS,cKS,PKS,Vxc) @@ -15,6 +15,7 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n integer,intent(in) :: max_diis integer,intent(in) :: guess_type logical,intent(in) :: mix + logical,intent(in) :: level_shift double precision,intent(in) :: thresh integer,intent(in) :: nNuc @@ -165,8 +166,9 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n end do call cpu_time(start_KS) - call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & - nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc) + call UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight, & + maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc) call cpu_time(end_KS) t_KS = end_KS - start_KS @@ -182,8 +184,9 @@ subroutine eDFT(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,n if(method == 'eDFT-UKS') then call cpu_time(start_KS) - call eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight,maxSCF,thresh,max_diis,guess_type,mix, & - nNuc,ZNuc,rNuc,ENuc,nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc) + call UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC(1:nCC,1:nEns-1),nGrid,weight, & + maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc,ENuc, & + nBas,AO,dAO,S,T,V,Hc,ERI,dipole_int,X,occnum,Cx_choice,doNcentered,Ew,eKS,cKS,PKS,Vxc) call cpu_time(end_KS) t_KS = end_KS - start_KS diff --git a/src/utils/level_shifting.f90 b/src/utils/level_shifting.f90 new file mode 100644 index 0000000..9e8e23c --- /dev/null +++ b/src/utils/level_shifting.f90 @@ -0,0 +1,37 @@ +subroutine level_shifting(level_shift,nBas,nO,S,c,F) + +! Perform level-shifting on the Fock matrix + + implicit none + +! Input variables + + double precision,intent(in) :: level_shift + integer,intent(in) :: nBas + integer,intent(in) :: nO + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: c(nBas,nBas) + +! Local variables + + double precision,allocatable :: F_MO(:,:) + double precision,allocatable :: Sc(:,:) + + integer :: a + +! Output variables + + double precision,intent(inout):: F(nBas,nBas) + + allocate(F_MO(nBas,nBas),Sc(nBas,nBas)) + + F_MO(:,:) = matmul(transpose(c),matmul(F,c)) + + do a=nO+1,nBas + F_MO(a,a) = F_MO(a,a) + level_shift + end do + + Sc(:,:) = matmul(S,c) + F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc))) + +end subroutine level_shifting