diff --git a/src/AOtoMO/AOtoMO_integral_transform.f90 b/src/AOtoMO/AOtoMO_integral_transform.f90 index 7463dc4..0b005fd 100644 --- a/src/AOtoMO/AOtoMO_integral_transform.f90 +++ b/src/AOtoMO/AOtoMO_integral_transform.f90 @@ -1,4 +1,4 @@ -subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI_MO_basis) +subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO,ERI_MO) ! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm ! bra and ket are the spin of (bra1 bra2|ket1 ket2) @@ -11,7 +11,7 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI integer,intent(in) :: bra1,bra2 integer,intent(in) :: ket1,ket2 integer,intent(in) :: nBas - double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin) ! Local variables @@ -20,7 +20,7 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI ! Output variables - double precision,intent(out) :: ERI_MO_basis(nBas,nBas,nBas,nBas) + double precision,intent(out) :: ERI_MO(nBas,nBas,nBas,nBas) ! Memory allocation @@ -28,76 +28,12 @@ subroutine AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,c,ERI_AO_basis,ERI ! Four-index transform via semi-direct O(N^5) algorithm -! scr(:,:,:,:) = 0d0 - -! ! do l=1,nBas -! ! do si=1,nBas -! ! do la=1,nBas -! ! do nu=1,nBas -! ! do mu=1,nBas -! ! scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket2) -! ! enddo -! ! enddo -! ! enddo -! ! enddo -! ! enddo - -! call dgemm ('N', 'N', nBas*nBas*nBas, nBas, nBas, 1d0, ERI_AO_basis(1,1,1,1), & -! size(ERI_AO_basis,1)*size(ERI_AO_basis,2)*size(ERI_AO_basis,3), c(1,1,ket2), & -! size(c,1), 0d0, scr, size(scr,1)*size(scr,2)*size(scr,3) ) - -! ERI_MO_basis(:,:,:,:) = 0d0 - -! ! do l=1,nBas -! ! do la=1,nBas -! ! do nu=1,nBas -! ! do i=1,nBas -! ! do mu=1,nBas -! ! ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra1)*scr(mu,nu,la,l) -! ! enddo -! ! enddo -! ! enddo -! ! enddo -! ! enddo - -! call dgemm ('T', 'N', nBas, nBas*nBas*nBas, nBas, 1d0, c(1,1,bra1), size(c,1), & -! scr(1,1,1,1), size(scr,1), 0d0, ERI_MO_basis, size(ERI_MO_basis,1)) - -! scr(:,:,:,:) = 0d0 - -! do l=1,nBas -! do k=1,nBas -! do la=1,nBas -! do nu=1,nBas -! do i=1,nBas -! scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra2) -! enddo -! enddo -! enddo -! enddo -! enddo - -! ERI_MO_basis(:,:,:,:) = 0d0 - -! do l=1,nBas -! do k=1,nBas -! do j=1,nBas -! do i=1,nBas -! do nu=1,nBas -! ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket1)*scr(i,nu,k,l) -! enddo -! ! write(11,'(I5,I5,I5,I5,F16.10)') i,j,k,l,ERI_MO_basis(i,j,k,l) -! enddo -! enddo -! enddo - ! enddo - - call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_AO_basis, nBas, c(1,1,bra2), size(c,1), 0d0, scr, nBas**3) + call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_AO, nBas, c(1,1,bra2), size(c,1), 0d0, scr, nBas**3) - call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,bra1), size(c,1), 0d0, ERI_MO_basis, nBas**3) + call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,bra1), size(c,1), 0d0, ERI_MO, nBas**3) - call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_MO_basis, nBas, c(1,1,ket1), size(c,1), 0d0, scr, nBas**3) + call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, ERI_MO, nBas, c(1,1,ket1), size(c,1), 0d0, scr, nBas**3) - call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,ket2), size(c,1), 0d0, ERI_MO_basis, nBas**3) + call dgemm ('T', 'N', nBas**3, nBas, nBas, 1d0, scr, nBas, c(1,1,ket2), size(c,1), 0d0, ERI_MO, nBas**3) -end subroutine AOtoMO_integral_transform +end subroutine diff --git a/src/AOtoMO/AOtoMO_oooa.f90 b/src/AOtoMO/AOtoMO_oooa.f90 index fc474e0..ce1a559 100644 --- a/src/AOtoMO/AOtoMO_oooa.f90 +++ b/src/AOtoMO/AOtoMO_oooa.f90 @@ -82,4 +82,4 @@ subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa) deallocate(scr1,scr2) -end subroutine AOtoMO_oooa +end subroutine diff --git a/src/AOtoMO/AOtoMO_oooo.f90 b/src/AOtoMO/AOtoMO_oooo.f90 index d9ebe47..e177290 100644 --- a/src/AOtoMO/AOtoMO_oooo.f90 +++ b/src/AOtoMO/AOtoMO_oooo.f90 @@ -82,4 +82,4 @@ subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo) deallocate(scr1,scr2) -end subroutine AOtoMO_oooo +end subroutine diff --git a/src/AOtoMO/AOtoMO_oooooo.f90 b/src/AOtoMO/AOtoMO_oooooo.f90 deleted file mode 100644 index e8bba04..0000000 --- a/src/AOtoMO/AOtoMO_oooooo.f90 +++ /dev/null @@ -1,135 +0,0 @@ -subroutine AOtoMO_oooooo(nBas,nO,cO,O,oooOooo) - -! AO to MO transformation of three-electron integrals for the block oooooo -! Semi-direct O(N^7) algorithm - - implicit none - -! Input variables - - integer,intent(in) :: nBas,nO - double precision,intent(in) :: cO(nBas,nO),O(nBas,nBas,nBas,nBas,nBas,nBas) - -! Local variables - - double precision,allocatable :: scr1(:,:,:,:,:,:),scr2(:,:,:,:,:,:) - integer :: mu,nu,la,si,ka,ta,i,j,k,l,m,n - -! Output variables - - double precision,intent(out) :: oooOooo(nO,nO,nO,nO,nO,nO) - -! Memory allocation - allocate(scr1(nBas,nBas,nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas,nBas,nBas)) - - write(*,*) - write(*,'(A42)') '------------------------------------------' - write(*,'(A42)') ' AO to MO transformation for oooooo block ' - write(*,'(A42)') '------------------------------------------' - write(*,*) - - scr1 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do si=1,nBas - do ka=1,nBas - do ta=1,nBas - do n=1,nO - scr1(mu,nu,la,si,ka,n) = scr1(mu,nu,la,si,ka,n) + O(mu,nu,la,si,ka,ta)*cO(ta,n) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - scr2 = 0d0 - do mu=1,nBas - do nu=1,nBas - do la=1,nBas - do si=1,nBas - do ka=1,nBas - do i=1,nO - do n=1,nO - scr2(i,nu,la,si,ka,n) = scr2(i,nu,la,si,ka,n) + cO(mu,i)*scr1(mu,nu,la,si,ka,n) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - scr1 = 0d0 - do nu=1,nBas - do la=1,nBas - do si=1,nBas - do ka=1,nBas - do i=1,nO - do m=1,nO - do n=1,nO - scr1(i,nu,la,si,m,n) = scr1(i,nu,la,si,m,n) + scr2(i,nu,la,si,m,n)*cO(ka,m) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - scr2 = 0d0 - do nu=1,nBas - do la=1,nBas - do si=1,nBas - do i=1,nO - do j=1,nO - do m=1,nO - do n=1,nO - scr2(i,j,la,si,m,n) = scr2(i,j,la,si,m,n) + cO(nu,j)*scr1(i,nu,la,si,m,n) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - scr1 = 0d0 - do la=1,nBas - do si=1,nBas - do i=1,nO - do j=1,nO - do l=1,nO - do m=1,nO - do n=1,nO - scr1(i,j,la,l,m,n) = scr1(i,j,la,l,m,n) + scr2(i,j,la,si,m,n)*cO(si,l) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - oooOooo = 0d0 - do si=1,nBas - do i=1,nO - do j=1,nO - do k=1,nO - do l=1,nO - do m=1,nO - do n=1,nO - oooOooo(i,j,k,l,m,n) = oooOooo(i,j,k,l,m,n) + cO(la,k)*scr1(i,j,la,k,m,n) - enddo - enddo - enddo - enddo - enddo - enddo - enddo - - deallocate(scr1,scr2) - -end subroutine AOtoMO_oooooo diff --git a/src/AOtoMO/AOtoMO_oovv.f90 b/src/AOtoMO/AOtoMO_oovv.f90 index 05365c1..e82a19f 100644 --- a/src/AOtoMO/AOtoMO_oovv.f90 +++ b/src/AOtoMO/AOtoMO_oovv.f90 @@ -74,4 +74,4 @@ subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv) enddo enddo -end subroutine AOtoMO_oovv +end subroutine diff --git a/src/AOtoMO/AOtoMO_transform.f90 b/src/AOtoMO/AOtoMO_transform.f90 index 7919084..410d3ab 100644 --- a/src/AOtoMO/AOtoMO_transform.f90 +++ b/src/AOtoMO/AOtoMO_transform.f90 @@ -15,4 +15,4 @@ subroutine AOtoMO_transform(nBas,c,A) A = matmul(transpose(c),matmul(A,c)) -end subroutine AOtoMO_transform +end subroutine diff --git a/src/AOtoMO/Coulomb_matrix_AO_basis.f90 b/src/AOtoMO/Coulomb_matrix_AO_basis.f90 index eaf6a3e..7268acb 100644 --- a/src/AOtoMO/Coulomb_matrix_AO_basis.f90 +++ b/src/AOtoMO/Coulomb_matrix_AO_basis.f90 @@ -31,4 +31,4 @@ subroutine Coulomb_matrix_AO_basis(nBas,P,G,J) enddo -end subroutine Coulomb_matrix_AO_basis +end subroutine diff --git a/src/AOtoMO/Coulomb_matrix_MO_basis.f90 b/src/AOtoMO/Coulomb_matrix_MO_basis.f90 index 1fea11e..d64bfe1 100644 --- a/src/AOtoMO/Coulomb_matrix_MO_basis.f90 +++ b/src/AOtoMO/Coulomb_matrix_MO_basis.f90 @@ -23,4 +23,4 @@ subroutine Coulomb_matrix_MO_basis(nBas,c,P,G,J) J = matmul(transpose(c),matmul(J,c)) -end subroutine Coulomb_matrix_MO_basis +end subroutine diff --git a/src/AOtoMO/Hartree_matrix_AO_basis.f90 b/src/AOtoMO/Hartree_matrix_AO_basis.f90 index 3ee7368..0d0e5a7 100644 --- a/src/AOtoMO/Hartree_matrix_AO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_AO_basis.f90 @@ -30,4 +30,4 @@ subroutine Hartree_matrix_AO_basis(nBas,P,Hc,G,H) enddo enddo -end subroutine Hartree_matrix_AO_basis +end subroutine diff --git a/src/AOtoMO/Hartree_matrix_MO_basis.f90 b/src/AOtoMO/Hartree_matrix_MO_basis.f90 index 6cf85bd..829c0e9 100644 --- a/src/AOtoMO/Hartree_matrix_MO_basis.f90 +++ b/src/AOtoMO/Hartree_matrix_MO_basis.f90 @@ -23,4 +23,4 @@ subroutine Hartree_matrix_MO_basis(nBas,c,P,Hc,G,H) H = matmul(transpose(c),matmul(H,c)) -end subroutine Hartree_matrix_MO_basis +end subroutine diff --git a/src/AOtoMO/MOtoAO_transform.f90 b/src/AOtoMO/MOtoAO_transform.f90 index b4a8b4f..005c84f 100644 --- a/src/AOtoMO/MOtoAO_transform.f90 +++ b/src/AOtoMO/MOtoAO_transform.f90 @@ -24,4 +24,4 @@ subroutine MOtoAO_transform(nBas,S,c,A) Sc = matmul(S,c) A = matmul(Sc,matmul(A,transpose(Sc))) -end subroutine MOtoAO_transform +end subroutine diff --git a/src/AOtoMO/exchange_matrix_AO_basis.f90 b/src/AOtoMO/exchange_matrix_AO_basis.f90 index 8a38c0e..6568030 100644 --- a/src/AOtoMO/exchange_matrix_AO_basis.f90 +++ b/src/AOtoMO/exchange_matrix_AO_basis.f90 @@ -30,4 +30,4 @@ subroutine exchange_matrix_AO_basis(nBas,P,ERI,K) enddo enddo -end subroutine exchange_matrix_AO_basis +end subroutine diff --git a/src/CI/CID.f90 b/src/CI/CID.f90 index b8f52d6..a640cb8 100644 --- a/src/CI/CID.f90 +++ b/src/CI/CID.f90 @@ -215,4 +215,4 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi write(*,*) endif -end subroutine CID +end subroutine diff --git a/src/CI/CIS.f90 b/src/CI/CIS.f90 index 172d966..bdaff1d 100644 --- a/src/CI/CIS.f90 +++ b/src/CI/CIS.f90 @@ -19,10 +19,11 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) logical :: dump_matrix = .false. logical :: dump_trans = .false. + logical :: dRPA = .false. integer :: ispin integer :: maxS = 10 double precision :: lambda - double precision,allocatable :: A(:,:),Omega(:) + double precision,allocatable :: A(:,:),Om(:) ! Hello world @@ -38,14 +39,14 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) ! Memory allocation - allocate(A(nS,nS),Omega(nS)) + allocate(A(nS,nS),Om(nS)) ! Compute CIS matrix if(singlet) then ispin = 1 - call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A) if(dump_matrix) then print*,'CIS matrix (singlet state)' @@ -53,9 +54,9 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) write(*,*) endif - call diagonalize_matrix(nS,A,Omega) - call print_excitation('CIS ',ispin,nS,Omega) - call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A)) + call diagonalize_matrix(nS,A,Om) + call print_excitation('CIS ',ispin,nS,Om) + call print_transition_vectors_ph(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,transpose(A),transpose(A)) if(dump_trans) then print*,'Singlet CIS transition vectors' @@ -66,14 +67,14 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) ! Compute CIS(D) correction maxS = min(maxS,nS) - if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS)) + if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Om(1:maxS),A(:,1:maxS)) endif if(triplet) then ispin = 2 - call phLR_A(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A) if(dump_matrix) then print*,'CIS matrix (triplet state)' @@ -81,9 +82,9 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) write(*,*) endif - call diagonalize_matrix(nS,A,Omega) - call print_excitation('CIS ',ispin,nS,Omega) - call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A)) + call diagonalize_matrix(nS,A,Om) + call print_excitation('CIS ',ispin,nS,Om) + call print_transition_vectors_ph(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Om,transpose(A),transpose(A)) if(dump_trans) then print*,'Triplet CIS transition vectors' @@ -94,7 +95,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) ! Compute CIS(D) correction maxS = min(maxS,nS) - if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS)) + if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Om(1:maxS),A(:,1:maxS)) endif diff --git a/src/CI/CISD.f90 b/src/CI/CISD.f90 index 02f7fd3..039d6b4 100644 --- a/src/CI/CISD.f90 +++ b/src/CI/CISD.f90 @@ -307,4 +307,4 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI write(*,*) endif -end subroutine CISD +end subroutine diff --git a/src/CI/D_correction.f90 b/src/CI/D_correction.f90 index ee2a7d2..63551d8 100644 --- a/src/CI/D_correction.f90 +++ b/src/CI/D_correction.f90 @@ -274,4 +274,4 @@ subroutine D_correction(ispin,nBasin,nCin,nOin,nVin,nRin,nSin,maxS,eHF,ERI,w,X) !------------------------------------------------------------------------ ! End of loop over single excitations !------------------------------------------------------------------------ -end subroutine D_correction +end subroutine diff --git a/src/CI/FCI.f90 b/src/CI/FCI.f90 index 1fcabbd..d2c44b5 100644 --- a/src/CI/FCI.f90 +++ b/src/CI/FCI.f90 @@ -28,4 +28,4 @@ subroutine FCI(nBas,nC,nO,nV,nR,ERI,e) ! Diagonalize FCI matrix -end subroutine FCI +end subroutine diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index 19bace6..6c51ea2 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -134,4 +134,4 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E endif -end subroutine UCIS +end subroutine diff --git a/src/GF/QP_graph_GF2.f90 b/src/GF/QP_graph_GF2.f90 index 1017505..2e273cf 100644 --- a/src/GF/QP_graph_GF2.f90 +++ b/src/GF/QP_graph_GF2.f90 @@ -73,4 +73,4 @@ subroutine QP_graph_GF2(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2lin,ERI,eGF2) end do -end subroutine QP_graph_GF2 +end subroutine diff --git a/src/GF/regularized_self_energy_GF2_diag.f90 b/src/GF/regularized_self_energy_GF2_diag.f90 index 634eb9d..ce424cb 100644 --- a/src/GF/regularized_self_energy_GF2_diag.f90 +++ b/src/GF/regularized_self_energy_GF2_diag.f90 @@ -83,4 +83,4 @@ subroutine regularized_self_energy_GF2_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,eGF2,ERI Z(:) = 1d0/(1d0 - Z(:)) -end subroutine regularized_self_energy_GF2_diag +end subroutine diff --git a/src/GW/Bethe_Salpeter_pp_so.f90 b/src/GW/Bethe_Salpeter_pp_so.f90 deleted file mode 100644 index b38c38b..0000000 --- a/src/GW/Bethe_Salpeter_pp_so.f90 +++ /dev/null @@ -1,98 +0,0 @@ -subroutine Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eW,eGW,EcBSE) - -! Compute the Bethe-Salpeter excitation energies at the pp level - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: TDA_W - logical,intent(in) :: TDA - logical,intent(in) :: singlet - logical,intent(in) :: triplet - - 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 - double precision,intent(in) :: eW(nBas) - double precision,intent(in) :: eGW(nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - -! Local variables - - integer :: ispin - integer :: isp_W - - integer :: nOO - integer :: nVV - - double precision :: EcRPA - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) - - double precision,allocatable :: Omega1(:) - double precision,allocatable :: X1(:,:) - double precision,allocatable :: Y1(:,:) - - double precision,allocatable :: Omega2(:) - double precision,allocatable :: X2(:,:) - double precision,allocatable :: Y2(:,:) - - double precision,allocatable :: WB(:,:) - double precision,allocatable :: WC(:,:) - double precision,allocatable :: WD(:,:) - -! Output variables - - double precision,intent(out) :: EcBSE(nspin) - -!--------------------------------- -! Compute RPA screening -!--------------------------------- - - isp_W = 3 - EcRPA = 0d0 - - ! Memory allocation - - allocate(OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS)) - - call phLR(isp_W,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) - - write(*,*) '****************' - write(*,*) '*** SpinOrbs ***' - write(*,*) '****************' - write(*,*) - - ispin = 1 - EcBSE(:) = 0d0 - - nOO = nO*(nO-1)/2 - nVV = nV*(nV-1)/2 - - allocate(Omega1(nVV),X1(nVV,nVV),Y1(nOO,nVV), & - Omega2(nOO),X2(nVV,nOO),Y2(nOO,nOO), & - WB(nVV,nOO),WC(nVV,nVV),WD(nOO,nOO)) - - if(.not.TDA) call static_screening_WB_pp(4,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WB) - call static_screening_WC_pp(4,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WC) - call static_screening_WD_pp(4,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,ERI,OmRPA,rho_RPA,WD) - - ! Compute BSE excitation energies - - call linear_response_pp_BSE(4,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & - Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin)) - - call print_excitation('pp-BSE (N+2)',isp_W,nVV,Omega1) - call print_excitation('pp-BSE (N-2)',isp_W,nOO,Omega2) - -end subroutine Bethe_Salpeter_pp_so diff --git a/src/GW/G0W0.f90 b/src/GW/G0W0.f90 index 1ddfbe2..9e00a42 100644 --- a/src/GW/G0W0.f90 +++ b/src/GW/G0W0.f90 @@ -1,6 +1,6 @@ -subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & - ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc) +subroutine G0W0(doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, & + singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO, & + dipole_int,PHF,cHF,eHF,Vxc) ! Perform G0W0 calculation @@ -13,10 +13,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS - logical,intent(in) :: COHSEX - logical,intent(in) :: BSE - logical,intent(in) :: BSE2 - logical,intent(in) :: ppBSE + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 + logical,intent(in) :: doppBSE logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE @@ -47,30 +46,25 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD ! Local variables logical :: print_W = .true. + logical :: dRPA integer :: ispin double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) - double precision :: EcppBSE(nspin) double precision :: EcGM + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: SigX(:) double precision,allocatable :: SigC(:) double precision,allocatable :: Z(:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) double precision,allocatable :: eGW(:) double precision,allocatable :: eGWlin(:) - integer :: nBas2 - integer :: nC2 - integer :: nO2 - integer :: nV2 - integer :: nR2 - integer :: nS2 - ! Output variables ! Hello world @@ -83,15 +77,9 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD ! Initialization + dRPA = .true. EcRPA = 0d0 -! COHSEX approximation - - if(COHSEX) then - write(*,*) 'COHSEX approximation activated!' - write(*,*) - end if - ! TDA for W if(TDA_W) then @@ -112,21 +100,25 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD ! Memory allocation - allocate(SigC(nBas),SigX(nBas),Z(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS),eGW(nBas),eGWlin(nBas)) + allocate(Aph(nS,nS),Bph(nS,nS),SigC(nBas),SigX(nBas),Z(nBas),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & + eGW(nBas),eGWlin(nBas)) !-------------------! ! Compute screening ! !-------------------! - call phLR(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) - if(print_W) call print_excitation('RPA@HF ',ispin,nS,OmRPA) + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + + if(print_W) call print_excitation('RPA@HF ',ispin,nS,Om) !--------------------------! ! Compute spectral weights ! !--------------------------! - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) !------------------------! ! Compute GW self-energy ! @@ -136,12 +128,12 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD if(regularize) then - call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) + call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC) + call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,Z) else - call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC,Z) + call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eHF,Om,rho,EcGM,SigC,Z) end if @@ -165,17 +157,20 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW,regularize) + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGWlin,eGW,regularize) ! Find all the roots of the QP equation if necessary - ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,OmRPA,rho_RPA,eGWlin) + ! call QP_roots(nBas,nC,nO,nV,nR,nS,eta,eHF,Om,rho,eGWlin) end if ! Compute the RPA correlation energy - call phLR(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) !--------------! ! Dump results ! @@ -185,17 +180,13 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD ! Deallocate memory -! deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) - -! Plot stuff - -! call plot_GW(nBas,nC,nO,nV,nR,nS,eHF,eGW,OmRPA,rho_RPA) + deallocate(SigC,Z,Om,XpY,XmY,rho,eGWlin) ! Perform BSE calculation - if(BSE) then + if(dophBSE) then - call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) + call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) if(exchange_kernel) then @@ -229,14 +220,14 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eHF,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (singlet) =',EcAC(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy (triplet) =',EcAC(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 correlation energy =',EcAC(1) + EcAC(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'AC@BSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (singlet) =',EcAC(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy (triplet) =',EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 correlation energy =',EcAC(1) + EcAC(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'AC@phBSE@G0W0 total energy =',ENuc + ERHF + EcAC(1) + EcAC(2),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) @@ -244,40 +235,19 @@ subroutine G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTD end if - if(ppBSE) then + if(doppBSE) then - call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcppBSE) + call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcBSE(1),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 correlation energy =',EcBSE(1) + 3d0*EcBSE(2),' au' + write(*,'(2X,A50,F20.10,A3)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2),' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! nBas2 = 2*nBas -! nO2 = 2*nO -! nV2 = 2*nV -! nC2 = 2*nC -! nR2 = 2*nR -! nS2 = nO2*nV2 -! -! allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) -! -! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) -! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW) -! call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI) -! -! call GW_ppBSE_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE) - end if -! if(BSE) call ufBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,eGW) -! if(BSE) call ufXBSE(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_MO,eHF,OmRPA,rho_RPA) - - if(BSE) call XBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) - - end subroutine diff --git a/src/GW/GW_self_energy.f90 b/src/GW/GW_self_energy.f90 index 687385f..9e7f852 100644 --- a/src/GW/GW_self_energy.f90 +++ b/src/GW/GW_self_energy.f90 @@ -1,4 +1,4 @@ -subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) +subroutine GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) ! Compute correlation part of the self-energy and the renormalization factor @@ -7,7 +7,6 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) ! Input variables - logical,intent(in) :: COHSEX double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -37,113 +36,72 @@ subroutine GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) Sig(:,:) = 0d0 Z(:) = 0d0 -!-----------------------------! -! COHSEX static approximation ! -!-----------------------------! - - if(COHSEX) then - - ! COHSEX: SEX of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - Sig(p,q) = Sig(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb) - end do - end do - end do - end do - - ! COHSEX: COH part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do r=nC+1,nBas-nR - do jb=1,nS - Sig(p,q) = Sig(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb) - end do - end do - end do - end do - - EcGM = 0d0 - do i=nC+1,nO - EcGM = EcGM + 0.5d0*Sig(i,i) - end do - - Z(:) = 0d0 - - else - !----------------! ! GW self-energy ! !----------------! - ! Occupied part of the correlation self-energy - - !$OMP PARALLEL & - !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & - !$OMP PRIVATE(jb,i,q,p,eps,num) & - !$OMP DEFAULT(NONE) - !$OMP DO - do q=nC+1,nBas-nR - do p=nC+1,nBas-nR - do jb=1,nS - do i=nC+1,nO - - eps = e(p) - e(i) + Omega(jb) - num = 2d0*rho(p,i,jb)*rho(q,i,jb) - Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) - if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - end do - end do - end do - end do - !$OMP END DO - !$OMP END PARALLEL - - ! Virtual part of the correlation self-energy - - !$OMP PARALLEL & - !$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & - !$OMP PRIVATE(jb,a,q,p,eps,num) & - !$OMP DEFAULT(NONE) - !$OMP DO - do q=nC+1,nBas-nR - do p=nC+1,nBas-nR - do jb=1,nS - do a=nO+1,nBas-nR - - eps = e(p) - e(a) - Omega(jb) - num = 2d0*rho(p,a,jb)*rho(q,a,jb) - Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) - if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - end do - end do - end do - end do - !$OMP END DO - !$OMP END PARALLEL - - ! Galitskii-Migdal correlation energy - - EcGM = 0d0 - do jb=1,nS - do a=nO+1,nBas-nR - do i=nC+1,nO - - eps = e(a) - e(i) + Omega(jb) - num = 4d0*rho(a,i,jb)*rho(a,i,jb) - EcGM = EcGM - num*eps/(eps**2 + eta**2) +! Occupied part of the correlation self-energy +!$OMP PARALLEL & +!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & +!$OMP PRIVATE(jb,i,q,p,eps,num) & +!$OMP DEFAULT(NONE) +!$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + do jb=1,nS + do i=nC+1,nO + + eps = e(p) - e(i) + Omega(jb) + num = 2d0*rho(p,i,jb)*rho(q,i,jb) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do end do + end do + end do + !$OMP END DO + !$OMP END PARALLEL + +! Virtual part of the correlation self-energy + +!$OMP PARALLEL & +!$OMP SHARED(Sig,rho,eta,nS,nC,nO,nBas,nR,e,Omega) & +!$OMP PRIVATE(jb,a,q,p,eps,num) & +!$OMP DEFAULT(NONE) +!$OMP DO + do q=nC+1,nBas-nR + do p=nC+1,nBas-nR + do jb=1,nS + do a=nO+1,nBas-nR + + eps = e(p) - e(a) - Omega(jb) + num = 2d0*rho(p,a,jb)*rho(q,a,jb) + Sig(p,q) = Sig(p,q) + num*eps/(eps**2 + eta**2) + if(p == q) Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + end do + end do + end do + end do +!$OMP END DO +!$OMP END PARALLEL + +! Galitskii-Migdal correlation energy + + EcGM = 0d0 + do jb=1,nS + do a=nO+1,nBas-nR + do i=nC+1,nO + + eps = e(a) - e(i) + Omega(jb) + num = 4d0*rho(a,i,jb)*rho(a,i,jb) + EcGM = EcGM - num*eps/(eps**2 + eta**2) + end do end do - - end if + end do ! Compute renormalization factor from derivative diff --git a/src/GW/GW_self_energy_diag.f90 b/src/GW/GW_self_energy_diag.f90 index 735591e..c8c0272 100644 --- a/src/GW/GW_self_energy_diag.f90 +++ b/src/GW/GW_self_energy_diag.f90 @@ -1,4 +1,4 @@ -subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) +subroutine GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,Sig,Z) ! Compute diagonal of the correlation part of the self-energy and the renormalization factor @@ -7,7 +7,6 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S ! Input variables - logical,intent(in) :: COHSEX double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -35,91 +34,54 @@ subroutine GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,S Sig(:) = 0d0 Z(:) = 0d0 -!----------------------------- -! COHSEX static self-energy -!----------------------------- - - if(COHSEX) then - - ! COHSEX: SEX part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - Sig(p) = Sig(p) + 4d0*rho(p,i,jb)**2/Omega(jb) - end do - end do - end do +!----------------! +! GW self-energy ! +!----------------! - ! COHSEX: COH part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do jb=1,nS - Sig(p) = Sig(p) - 2d0*rho(p,q,jb)**2/Omega(jb) - end do - end do - end do +! Occupied part of the correlation self-energy - ! GM correlation energy - - EcGM = 0d0 + do p=nC+1,nBas-nR do i=nC+1,nO - EcGM = EcGM - Sig(i) - end do + do jb=1,nS -!----------------------------- -! GW self-energy -!----------------------------- + eps = e(p) - e(i) + Omega(jb) + num = 2d0*rho(p,i,jb)**2 + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - else - - ! Occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - - eps = e(p) - e(i) + Omega(jb) - num = 2d0*rho(p,i,jb)**2 - Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) - Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - - end do end do end do - - ! Virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS + end do - eps = e(p) - e(a) - Omega(jb) - num = 2d0*rho(p,a,jb)**2 - Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) - Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 +! Virtual part of the correlation self-energy + + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do jb=1,nS + + eps = e(p) - e(a) - Omega(jb) + num = 2d0*rho(p,a,jb)**2 + Sig(p) = Sig(p) + num*eps/(eps**2 + eta**2) + Z(p) = Z(p) - num*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do end do end do + end do - ! GM correlation energy +! Galitskii-Migdal correlation energy - EcGM = 0d0 - do i=nC+1,nO - do a=nO+1,nBas-nR - do jb=1,nS + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do jb=1,nS - eps = e(a) - e(i) + Omega(jb) - num = 4d0*rho(a,i,jb)**2 - EcGM = EcGM - num*eps/(eps**2 + eta**2) + eps = e(a) - e(i) + Omega(jb) + num = 4d0*rho(a,i,jb)**2 + EcGM = EcGM - num*eps/(eps**2 + eta**2) - end do end do end do - - end if + end do ! Compute renormalization factor from derivative diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 658c255..425474f 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -1,4 +1,4 @@ -subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & +subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, & singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI_AO,ERI_MO,dipole_int,PHF, & cHF,eHF,Vxc) @@ -17,15 +17,14 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS - logical,intent(in) :: COHSEX - logical,intent(in) :: BSE - logical,intent(in) :: BSE2 + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn - logical,intent(in) :: ppBSE + logical,intent(in) :: doppBSE logical,intent(in) :: singlet logical,intent(in) :: triplet logical,intent(in) :: linearize @@ -49,6 +48,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Local variables logical :: linear_mixing + logical :: dRPA = .true. integer :: ispin integer :: nSCF integer :: n_diis @@ -58,10 +58,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, double precision :: EcRPA double precision :: EcBSE(nspin) double precision :: EcAC(nspin) - double precision :: EcppBSE(nspin) double precision :: EcGM double precision :: alpha - double precision :: Dpijb,Dpajb + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) double precision,allocatable :: error_diis(:,:) double precision,allocatable :: e_diis(:,:) double precision,allocatable :: eGW(:) @@ -69,21 +69,10 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, double precision,allocatable :: Z(:) double precision,allocatable :: SigX(:) double precision,allocatable :: SigC(:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) - - double precision,allocatable :: eGWlin(:) - -! integer :: nBas2 -! integer :: nC2 -! integer :: nO2 -! integer :: nV2 -! integer :: nR2 -! integer :: nS2 - -! double precision,allocatable :: seHF(:),seGW(:),sERI(:,:,:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) ! Hello world @@ -93,13 +82,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, write(*,*)'************************************************' write(*,*) -! COHSEX approximation - - if(COHSEX) then - write(*,*) 'COHSEX approximation activated!' - write(*,*) - end if - ! TDA for W if(TDA_W) then @@ -121,8 +103,8 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Memory allocation - allocate(eGW(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas),OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS), & - rho_RPA(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis),eGWlin(nBas)) + allocate(Aph(nS,nS),Bph(nS,nS),eGW(nBas),eOld(nBas),Z(nBas),SigX(nBas),SigC(nBas), & + Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS),error_diis(nBas,max_diis),e_diis(nBas,max_diis)) ! Compute the exchange part of the self-energy @@ -140,8 +122,6 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, eOld(:) = eGW(:) Z(:) = 1d0 rcond = 0d0 - - !------------------------------------------------------------------------ ! Main loop @@ -151,28 +131,31 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Compute screening - call phLR(ispin,.true.,TDA_W,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,EcRPA,OmRPA,XpY_RPA,XmY_RPA) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) ! Compute spectral weights - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) ! Compute correlation part of the self-energy if(regularize) then - call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + call regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC) + call renormalization_factor_SRG(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,Z) else - call GW_self_energy_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z) + call GW_self_energy_diag(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) end if ! Solve the quasi-particle equation - eGWlin(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) + eGW(:) = eHF(:) + SigX(:) + SigC(:) - Vxc(:) ! Linearized or graphical solution? @@ -181,14 +164,14 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' write(*,*) - eGW(:) = eGWlin(:) + eGW(:) = eGW(:) else write(*,*) ' *** Quasiparticle energies obtained by root search (experimental) *** ' write(*,*) - call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,OmRPA,rho_RPA,eGWlin,eGW,regularize) + call QP_graph(nBas,nC,nO,nV,nR,nS,eta,eHF,SigX,Vxc,Om,rho,eGW,eGW,regularize) end if @@ -250,13 +233,13 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Deallocate memory - deallocate(eOld,Z,SigC,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error_diis,e_diis) + deallocate(Aph,Bph,eOld,Z,SigC,Om,XpY,XmY,rho,error_diis,e_diis) ! Perform BSE calculation - if(BSE) then + if(dophBSE) then - call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) + call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eGW,eGW,EcBSE) if(exchange_kernel) then @@ -290,7 +273,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -305,34 +288,19 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, end if - if(ppBSE) then + if(doppBSE) then - call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcppBSE) + call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF,eGW,EcBSE) write(*,*) write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcppBSE(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcppBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy (triplet) =',3d0*EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@evGW total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) -! nBas2 = 2*nBas -! nO2 = 2*nO -! nV2 = 2*nV -! nC2 = 2*nC -! nR2 = 2*nR -! nS2 = nO2*nV2 -! -! allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) -! -! call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) -! call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW) -! call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI) -! -! call GW_ppBSE_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE) - end if end subroutine diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 6faa700..2da86c9 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -1,4 +1,4 @@ -subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE,BSE2,TDA_W,TDA,dBSE,dTDA,evDyn, & +subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dophBSE,dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,doppBSE, & singlet,triplet,eta,regularize,nNuc,ZNuc,rNuc,ENuc,nBas,nC,nO,nV,nR,nS,ERHF,S,X,T,V,Hc,ERI_AO, & ERI_MO,dipole_int_AO,dipole_int_MO,PHF,cHF,eHF) @@ -15,14 +15,14 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, logical,intent(in) :: doACFDT logical,intent(in) :: exchange_kernel logical,intent(in) :: doXBS - logical,intent(in) :: COHSEX - logical,intent(in) :: BSE - logical,intent(in) :: BSE2 + logical,intent(in) :: dophBSE + logical,intent(in) :: dophBSE2 logical,intent(in) :: TDA_W logical,intent(in) :: TDA logical,intent(in) :: dBSE logical,intent(in) :: dTDA logical,intent(in) :: evDyn + logical,intent(in) :: doppBSE logical,intent(in) :: singlet logical,intent(in) :: triplet double precision,intent(in) :: eta @@ -73,13 +73,16 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, double precision,external :: trace_matrix double precision :: dipole(ncart) + logical :: dRPA = .true. logical :: print_W = .true. double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) + double precision,allocatable :: Aph(:,:) + double precision,allocatable :: Bph(:,:) + double precision,allocatable :: Om(:) + double precision,allocatable :: XpY(:,:) + double precision,allocatable :: XmY(:,:) + double precision,allocatable :: rho(:,:,:) double precision,allocatable :: c(:,:) double precision,allocatable :: cp(:,:) double precision,allocatable :: eGW(:) @@ -112,13 +115,6 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, nBasSq = nBas*nBas -! COHSEX approximation - - if(COHSEX) then - write(*,*) 'COHSEX approximation activated!' - write(*,*) - end if - ! TDA for W if(TDA_W) then @@ -137,7 +133,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, allocate(eGW(nBas),eOld(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas), & J(nBas,nBas),K(nBas,nBas),SigC(nBas,nBas),SigCp(nBas,nBas),SigCm(nBas,nBas),Z(nBas), & - OmRPA(nS),XpY_RPA(nS,nS),XmY_RPA(nS,nS),rho_RPA(nBas,nBas,nS), & + Aph(nS,nS),Bph(nS,nS),Om(nS),XpY(nS,nS),XmY(nS,nS),rho(nBas,nBas,nS), & error(nBas,nBas),error_diis(nBasSq,max_diis),F_diis(nBasSq,max_diis)) ! Initialization @@ -178,21 +174,24 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Compute linear response - call phLR(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) + call phLR_A(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI_MO,Aph) + if(.not.TDA_W) call phLR_B(ispin,dRPA,nBas,nC,nO,nV,nR,nS,1d0,ERI_MO,Bph) + + call phLR(TDA_W,nS,Aph,Bph,EcRPA,Om,XpY,XmY) + if(print_W) call print_excitation('RPA@qsGW ',ispin,nS,Om) ! Compute correlation part of the self-energy - call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY_RPA,rho_RPA) + call GW_excitation_density(nBas,nC,nO,nR,nS,ERI_MO,XpY,rho) if(regularize) then - call regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC) - call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,Z) + call regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC) + call regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,Z) else - call GW_self_energy(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eGW,OmRPA,rho_RPA,EcGM,SigC,Z) + call GW_self_energy(eta,nBas,nC,nO,nV,nR,nS,eGW,Om,rho,EcGM,SigC,Z) endif @@ -291,13 +290,13 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, ! Deallocate memory - deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,error,error_diis,F_diis) + deallocate(c,cp,P,F,Fp,J,K,SigC,SigCp,SigCm,Z,Om,XpY,XmY,rho,error,error_diis,F_diis) ! Perform BSE calculation - if(BSE) then + if(dophBSE) then - call GW_phBSE(BSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) + call GW_phBSE(dophBSE2,TDA_W,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eGW,eGW,EcBSE) if(exchange_kernel) then @@ -331,7 +330,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, end if - call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) + call GW_phACFDT(exchange_kernel,doXBS,.true.,TDA_W,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,eGW,eGW,EcAC) write(*,*) write(*,*)'-------------------------------------------------------------------------------' @@ -346,4 +345,19 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, end if + if(doppBSE) then + + call GW_ppBSE(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int_MO,eHF,eGW,EcBSE) + + write(*,*) + write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (singlet) =',EcBSE(1) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy (triplet) =',3d0*EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW correlation energy =',EcBSE(1) + 3d0*EcBSE(2) + write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@qsGW total energy =',ENuc + ERHF + EcBSE(1) + 3d0*EcBSE(2) + write(*,*)'-------------------------------------------------------------------------------' + write(*,*) + + end if + end subroutine diff --git a/src/GW/regularized_renormalization_factor.f90 b/src/GW/regularized_renormalization_factor.f90 index f1ddc05..5e141a9 100644 --- a/src/GW/regularized_renormalization_factor.f90 +++ b/src/GW/regularized_renormalization_factor.f90 @@ -1,4 +1,4 @@ -subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) +subroutine regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) ! Compute the regularized version of the GW renormalization factor @@ -7,7 +7,6 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O ! Input variables - logical,intent(in) :: COHSEX double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -41,45 +40,34 @@ subroutine regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,O kappa = 1d0 -! static COHSEX approximation +! Occupied part of the correlation self-energy - if(COHSEX) then - - Z(:) = 1d0 - return - - else - - ! Occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(p) - e(i) + Omega(jb) - fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) - Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk - end do + do p=nC+1,nBas-nR + do i=nC+1,nO + do jb=1,nS + eps = e(p) - e(i) + Omega(jb) + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) + Z(p) = Z(p) - 2d0*rho(p,i,jb)**2*dfk end do end do + end do - ! Virtual part of the correlation self-energy +! Virtual part of the correlation self-energy - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(p) - e(a) - Omega(jb) - fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) - Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk - end do + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do jb=1,nS + eps = e(p) - e(a) - Omega(jb) + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + dfk = - fk/eps + 4d0*kappa**2*exp(-2d0*eps**2/kappa**2) + Z(p) = Z(p) - 2d0*rho(p,a,jb)**2*dfk end do end do - - end if + end do ! Compute renormalization factor from derivative of SigC Z(:) = 1d0/(1d0 - Z(:)) -end subroutine regularized_renormalization_factor +end subroutine diff --git a/src/GW/regularized_self_energy_correlation.f90 b/src/GW/regularized_self_energy_correlation.f90 index 16808a9..bd3cc18 100644 --- a/src/GW/regularized_self_energy_correlation.f90 +++ b/src/GW/regularized_self_energy_correlation.f90 @@ -1,4 +1,4 @@ -subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) +subroutine regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) ! Compute correlation part of the regularized self-energy @@ -7,7 +7,6 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e, ! Input variables - logical,intent(in) :: COHSEX double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -44,88 +43,49 @@ subroutine regularized_self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e, kappa = 1d0 -!-----------------------------! -! COHSEX static approximation ! -!-----------------------------! - - if(COHSEX) then - - ! COHSEX: SEX of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - SigC(p,q) = SigC(p,q) + 4d0*rho(p,i,jb)*rho(q,i,jb)/Omega(jb) - end do - end do - end do - end do - - ! COHSEX: COH part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do r=nC+1,nBas-nR - do jb=1,nS - SigC(p,q) = SigC(p,q) - 2d0*rho(p,r,jb)*rho(q,r,jb)/Omega(jb) - end do - end do - end do - end do - - EcGM = 0d0 - do i=nC+1,nO - EcGM = EcGM + 0.5d0*SigC(i,i) - end do - - else - !----------------! ! GW self-energy ! !----------------! - ! Occupied part of the correlation self-energy +! Occupied part of the correlation self-energy - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(p) - e(i) + Omega(jb) - fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk - end do - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(p) - e(a) - Omega(jb) - fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk - end do + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR + do i=nC+1,nO + do jb=1,nS + eps = e(p) - e(i) + Omega(jb) + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + SigC(p,q) = SigC(p,q) + 2d0*rho(p,i,jb)*rho(q,i,jb)*fk end do end do end do + end do - ! GM correlation energy +! Virtual part of the correlation self-energy - EcGM = 0d0 - do i=nC+1,nO + do p=nC+1,nBas-nR + do q=nC+1,nBas-nR do a=nO+1,nBas-nR do jb=1,nS - eps = e(a) - e(i) + Omega(jb) + eps = e(p) - e(a) - Omega(jb) fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps - EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk + SigC(p,q) = SigC(p,q) + 2d0*rho(p,a,jb)*rho(q,a,jb)*fk end do end do end do + end do - end if +! GM correlation energy -end subroutine regularized_self_energy_correlation + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do jb=1,nS + eps = e(a) - e(i) + Omega(jb) + fk = (1d0 - exp(-2d0*eps**2/kappa**2))/eps + EcGM = EcGM - 4d0*rho(a,i,jb)**2*fk + end do + end do + end do + +end subroutine diff --git a/src/GW/regularized_self_energy_correlation_diag.f90 b/src/GW/regularized_self_energy_correlation_diag.f90 index 545c73a..9c4fb7d 100644 --- a/src/GW/regularized_self_energy_correlation_diag.f90 +++ b/src/GW/regularized_self_energy_correlation_diag.f90 @@ -1,4 +1,4 @@ -subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) +subroutine regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) ! Compute diagonal of the correlation part of the regularized self-energy @@ -7,7 +7,6 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, ! Input variables - logical,intent(in) :: COHSEX double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -33,78 +32,41 @@ subroutine regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR, SigC(:) = 0d0 -!----------------------------- -! COHSEX static self-energy -!----------------------------- - - if(COHSEX) then - - ! COHSEX: SEX part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - SigC(p) = SigC(p) + 4d0*rho(p,i,jb)**2/Omega(jb) - end do - end do - end do - - ! COHSEX: COH part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do jb=1,nS - SigC(p) = SigC(p) - 2d0*rho(p,q,jb)**2/Omega(jb) - end do - end do - end do - - ! GM correlation energy - - EcGM = 0d0 - do i=nC+1,nO - EcGM = EcGM - SigC(i) - end do - !----------------------------- ! GW self-energy !----------------------------- + +! Occupied part of the correlation self-energy - else - - ! Occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - Dpijb = e(p) - e(i) + Omega(jb) - SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*(1d0 - exp(-2d0*eta*Dpijb*Dpijb))/Dpijb - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - Dpajb = e(p) - e(a) - Omega(jb) - SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*(1d0 - exp(-2d0*eta*Dpajb*Dpajb))/Dpajb - end do - end do - end do - - ! GM correlation energy - - EcGM = 0d0 + do p=nC+1,nBas-nR do i=nC+1,nO - do a=nO+1,nBas-nR - do jb=1,nS - EcGM = EcGM - 4d0*rho(a,i,jb)**2 - end do + do jb=1,nS + Dpijb = e(p) - e(i) + Omega(jb) + SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2*(1d0 - exp(-2d0*eta*Dpijb*Dpijb))/Dpijb end do end do + end do - end if +! Virtual part of the correlation self-energy -end subroutine regularized_self_energy_correlation_diag + do p=nC+1,nBas-nR + do a=nO+1,nBas-nR + do jb=1,nS + Dpajb = e(p) - e(a) - Omega(jb) + SigC(p) = SigC(p) + 2d0*rho(p,a,jb)**2*(1d0 - exp(-2d0*eta*Dpajb*Dpajb))/Dpajb + end do + end do + end do + +! Galitskii-Migdal correlation energy + + EcGM = 0d0 + do i=nC+1,nO + do a=nO+1,nBas-nR + do jb=1,nS + EcGM = EcGM - 4d0*rho(a,i,jb)**2 + end do + end do + end do + +end subroutine diff --git a/src/GW/renormalization_factor_so.f90 b/src/GW/renormalization_factor_so.f90 deleted file mode 100644 index 1484f17..0000000 --- a/src/GW/renormalization_factor_so.f90 +++ /dev/null @@ -1,72 +0,0 @@ -subroutine renormalization_factor_so(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Z) - -! Compute renormalization factor for GW - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: COHSEX - 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 - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: p,i,a,jb - double precision :: eps - -! Output variables - - double precision,intent(out) :: Z(nBas) - -! Initialize - - Z(:) = 0d0 - -! static COHSEX approximation - - if(COHSEX) then - - Z(:) = 1d0 - return - - else - - ! Occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(p) - e(i) + Omega(jb) - Z(p) = Z(p) - 1d0*rho(p,i,jb)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(p) - e(a) - Omega(jb) - Z(p) = Z(p) - 1d0*rho(p,a,jb)**2*(eps/(eps**2 + eta**2))**2 - end do - end do - end do - - end if - -! Compute renormalization factor from derivative of SigC - - Z(:) = 1d0/(1d0 - Z(:)) - -end subroutine renormalization_factor_so diff --git a/src/GW/self_energy_correlation_diag_so.f90 b/src/GW/self_energy_correlation_diag_so.f90 deleted file mode 100644 index b26c2d6..0000000 --- a/src/GW/self_energy_correlation_diag_so.f90 +++ /dev/null @@ -1,111 +0,0 @@ -subroutine self_energy_correlation_diag_so(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,EcGM,SigC) - -! Compute diagonal of the correlation part of the self-energy - - implicit none - include 'parameters.h' - -! Input variables - - logical,intent(in) :: COHSEX - 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 - double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega(nS) - double precision,intent(in) :: rho(nBas,nBas,nS) - -! Local variables - - integer :: i,a,p,q,jb - double precision :: eps - -! Output variables - - double precision,intent(out) :: SigC(nBas) - double precision,intent(out) :: EcGM - -! Initialize - - SigC(:) = 0d0 - -!----------------------------- -! COHSEX static self-energy -!----------------------------- - - if(COHSEX) then - - ! COHSEX: SEX part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - SigC(p) = SigC(p) + 2d0*rho(p,i,jb)**2/Omega(jb) - end do - end do - end do - - ! COHSEX: COH part of the COHSEX correlation self-energy - - do p=nC+1,nBas-nR - do q=nC+1,nBas-nR - do jb=1,nS - SigC(p) = SigC(p) - rho(p,q,jb)**2/Omega(jb) - end do - end do - end do - - ! GM correlation energy - - EcGM = 0d0 - do i=nC+1,nO - EcGM = EcGM - SigC(i) - end do - -!----------------------------- -! GW self-energy -!----------------------------- - - else - - ! Occupied part of the correlation self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do jb=1,nS - eps = e(p) - e(i) + Omega(jb) - SigC(p) = SigC(p) + rho(p,i,jb)**2*eps/(eps**2 + eta**2) - end do - end do - end do - - ! Virtual part of the correlation self-energy - - do p=nC+1,nBas-nR - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(p) - e(a) - Omega(jb) - SigC(p) = SigC(p) + rho(p,a,jb)**2*eps/(eps**2 + eta**2) - end do - end do - end do - - ! GM correlation energy - - EcGM = 0d0 - do i=nC+1,nO - do a=nO+1,nBas-nR - do jb=1,nS - eps = e(a) - e(i) + Omega(jb) - EcGM = EcGM - 2d0*rho(a,i,jb)*rho(a,i,jb)*eps/(eps**2 + eta**2) - end do - end do - end do - - end if - -end subroutine self_energy_correlation_diag_so diff --git a/src/GW/soG0W0.f90 b/src/GW/soG0W0.f90 deleted file mode 100644 index 77f58e3..0000000 --- a/src/GW/soG0W0.f90 +++ /dev/null @@ -1,211 +0,0 @@ -subroutine soG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,evDyn,ppBSE, & - singlet,triplet,linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, & - ERI_AO,ERI_MO,dipole_int,PHF,cHF,eHF,Vxc,eGW) - -! Perform G0W0 calculation in the spin-orbital basis - - implicit none - include 'parameters.h' - include 'quadrature.h' - -! Input variables - - logical,intent(in) :: doACFDT - logical,intent(in) :: exchange_kernel - logical,intent(in) :: doXBS - logical,intent(in) :: COHSEX - logical,intent(in) :: BSE - logical,intent(in) :: ppBSE - logical,intent(in) :: TDA_W - logical,intent(in) :: TDA - logical,intent(in) :: dBSE - logical,intent(in) :: dTDA - logical,intent(in) :: evDyn - logical,intent(in) :: singlet - logical,intent(in) :: triplet - logical,intent(in) :: linearize - double precision,intent(in) :: eta - logical,intent(in) :: regularize - - 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) :: ENuc - double precision,intent(in) :: ERHF - double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) - double precision,intent(in) :: Vxc(nBas) - double precision,intent(in) :: eHF(nBas) - double precision,intent(in) :: cHF(nBas,nBas) - double precision,intent(in) :: PHF(nBas,nBas) - -! Local variables - - logical :: print_W = .true. - integer :: ispin - double precision :: EcRPA - double precision :: EcBSE(nspin) - double precision :: EcAC(nspin) - double precision :: EcppBSE(nspin) - double precision :: EcGM - double precision,allocatable :: SigC(:) - double precision,allocatable :: Z(:) - double precision,allocatable :: OmRPA(:) - double precision,allocatable :: XpY_RPA(:,:) - double precision,allocatable :: XmY_RPA(:,:) - double precision,allocatable :: rho_RPA(:,:,:) - - double precision,allocatable :: eGWlin(:) - - integer :: nBas2 - integer :: nC2 - integer :: nO2 - integer :: nV2 - integer :: nR2 - integer :: nS2 - - double precision,allocatable :: seHF(:),seGW(:),sERI(:,:,:,:) - -! Output variables - - double precision :: eGW(nBas) - -! Hello world - - write(*,*) - write(*,*)'************************************************' - write(*,*)'| One-shot soG0W0 calculation |' - write(*,*)'************************************************' - write(*,*) - -! Initialization - - EcRPA = 0d0 - -! COHSEX approximation - - if(COHSEX) then - write(*,*) 'COHSEX approximation activated!' - write(*,*) - end if - -! TDA for W - - if(TDA_W) then - write(*,*) 'Tamm-Dancoff approximation for dynamic screening!' - write(*,*) - end if - -! TDA - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) - end if - -! spatial to spin transformation - - nBas2 = 2*nBas - nO2 = 2*nO - nV2 = 2*nV - nC2 = 2*nC - nR2 = 2*nR - nS2 = nO2*nV2 - - allocate(seHF(nBas2),seGW(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) - - call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) - call spatial_to_spin_MO_energy(nBas,eGW,nBas2,seGW) - call spatial_to_spin_ERI(nBas,ERI_MO,nBas2,sERI) - -! Spin manifold - - ispin = 3 - -! Memory allocation - - allocate(SigC(nBas2),Z(nBas2),OmRPA(nS2),XpY_RPA(nS2,nS2),XmY_RPA(nS2,nS2),rho_RPA(nBas2,nBas2,nS2),eGWlin(nBas2)) - -!-------------------! -! Compute screening ! -!-------------------! - - call phLR(ispin,.true.,TDA_W,eta,nBas2,nC2,nO2,nV2,nR2,nS2,1d0,seHF,sERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - - if(print_W) call print_excitation('RPA@HF ',ispin,nS2,OmRPA) - -!--------------------------! -! Compute spectral weights ! -!--------------------------! - - call GW_excitation_density(nBas2,nC2,nO2,nR2,nS2,sERI,XpY_RPA,rho_RPA) - -!------------------------! -! Compute GW self-energy ! -!------------------------! - - if(regularize) then - - call regularized_self_energy_correlation_diag(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,EcGM,SigC) - call regularized_renormalization_factor(COHSEX,eta,nBas,nC,nO,nV,nR,nS,eHF,OmRPA,rho_RPA,Z) - - else - - call self_energy_correlation_diag_so(COHSEX,eta,nBas2,nC2,nO2,nV2,nR2,nS2,seHF,OmRPA,rho_RPA,EcGM,SigC) - call renormalization_factor_so(COHSEX,eta,nBas2,nC2,nO2,nV2,nR2,nS2,seHF,OmRPA,rho_RPA,Z) - - end if - -!-----------------------------------! -! Solve the quasi-particle equation ! -!-----------------------------------! - - eGWlin(:) = seHF(:) + Z(:)*SigC(:) - - ! Linearized or graphical solution? - - if(linearize) then - - write(*,*) ' *** Quasiparticle energies obtained by linearization *** ' - write(*,*) - - seGW(:) = eGWlin(:) - - end if - -! Compute the RPA correlation energy - - call phLR(ispin,.true.,TDA_W,eta,nBas2,nC2,nO2,nV2,nR2,nS2,1d0,seGW,sERI,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - -!--------------! -! Dump results ! -!--------------! - - call print_G0W0(nBas2,nO2,seHF,ENuc,ERHF,SigC,Z,seGW,EcRPA,EcGM) - -! Deallocate memory - - deallocate(SigC,Z,OmRPA,XpY_RPA,XmY_RPA,rho_RPA,eGWlin) - -! Perform BSE calculation - - if(ppBSE) then - - call Bethe_Salpeter_pp_so(TDA_W,TDA,singlet,triplet,eta,nBas2,nC2,nO2,nV2,nR2,nS2,sERI,dipole_int,seHF,seGW,EcppBSE) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (singlet) =',EcppBSE(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy (triplet) =',3d0*EcppBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 correlation energy =',EcppBSE(1) + 3d0*EcppBSE(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppBSE@G0W0 total energy =',ENuc + ERHF + EcppBSE(1) + 3d0*EcppBSE(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - - end if - -end subroutine diff --git a/src/HF/MOM_overlap.f90 b/src/HF/MOM_overlap.f90 index 3ca91a0..0e97fb4 100644 --- a/src/HF/MOM_overlap.f90 +++ b/src/HF/MOM_overlap.f90 @@ -49,4 +49,4 @@ subroutine MOM_overlap(nBas,nO,S,cG,c,ON) -end subroutine MOM_overlap +end subroutine diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index e54d3f2..19e9b60 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -204,4 +204,4 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rNuc,ENuc call mo_fock_exchange_potential(nBas,c,P,ERI,Vx) -end subroutine RHF +end subroutine diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index d3912eb..2b3d949 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -258,4 +258,4 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,level_shift,nNuc,ZNuc,rNuc, call mo_fock_exchange_potential(nBas,c(:,:,ispin),P(:,:,ispin),ERI,Vx(:,ispin)) end do -end subroutine UHF +end subroutine diff --git a/src/HF/UHF_stability.f90 b/src/HF/UHF_stability.f90 index db77387..8550855 100644 --- a/src/HF/UHF_stability.f90 +++ b/src/HF/UHF_stability.f90 @@ -171,4 +171,4 @@ subroutine UHF_stability(nBas,nC,nO,nV,nR,nS,eHF,ERI_aaaa,ERI_aabb,ERI_bbbb) write(*,*)'-------------------------------------------------------------' write(*,*) -end subroutine UHF_stability +end subroutine diff --git a/src/HF/density.f90 b/src/HF/density.f90 index 8e6a574..a365405 100644 --- a/src/HF/density.f90 +++ b/src/HF/density.f90 @@ -32,4 +32,4 @@ subroutine density(nGrid,nBas,P,AO,rho) enddo enddo -end subroutine density +end subroutine diff --git a/src/HF/density_matrix.f90 b/src/HF/density_matrix.f90 index 6f1e6a7..5d40950 100644 --- a/src/HF/density_matrix.f90 +++ b/src/HF/density_matrix.f90 @@ -27,4 +27,4 @@ subroutine density_matrix(nBas,ON,c,P) enddo enddo -end subroutine density_matrix +end subroutine diff --git a/src/HF/dipole_moment.f90 b/src/HF/dipole_moment.f90 index 2419f08..656a39b 100644 --- a/src/HF/dipole_moment.f90 +++ b/src/HF/dipole_moment.f90 @@ -50,4 +50,4 @@ subroutine dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int,dipole) enddo -end subroutine dipole_moment +end subroutine diff --git a/src/HF/mix_guess.f90 b/src/HF/mix_guess.f90 index 3cca61a..9fc2ec0 100644 --- a/src/HF/mix_guess.f90 +++ b/src/HF/mix_guess.f90 @@ -45,4 +45,4 @@ subroutine mix_guess(nBas,nO,c) c(:,nO(2),2) = HOMOb(:) c(:,nO(2)+1,2) = LUMOb(:) -end subroutine mix_guess +end subroutine diff --git a/src/HF/mo_fock_exchange_potential.f90 b/src/HF/mo_fock_exchange_potential.f90 index d07fcce..c30357e 100644 --- a/src/HF/mo_fock_exchange_potential.f90 +++ b/src/HF/mo_fock_exchange_potential.f90 @@ -38,4 +38,4 @@ subroutine mo_fock_exchange_potential(nBas,c,P,ERI,Vx) deallocate(Fx) -end subroutine mo_fock_exchange_potential +end subroutine diff --git a/src/HF/print_RHF.f90 b/src/HF/print_RHF.f90 index 9b88837..4ca1945 100644 --- a/src/HF/print_RHF.f90 +++ b/src/HF/print_RHF.f90 @@ -74,6 +74,4 @@ subroutine print_RHF(nBas,nO,eHF,cHF,ENuc,ET,EV,EJ,EK,ERHF,dipole) call matout(nBas,1,eHF) write(*,*) -end subroutine print_RHF - - +end subroutine diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index 17834a6..46a0551 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -124,4 +124,4 @@ subroutine print_UHF(nBas,nO,Ov,e,c,ENuc,ET,EV,EJ,Ex,EUHF,dipole) call matout(nBas,1,e(:,2)) write(*,*) -end subroutine print_UHF +end subroutine diff --git a/src/LR/print_excitation.f90 b/src/LR/print_excitation.f90 index acf7045..c28cb4b 100644 --- a/src/LR/print_excitation.f90 +++ b/src/LR/print_excitation.f90 @@ -41,6 +41,4 @@ subroutine print_excitation(method,ispin,nS,Omega) write(*,*)'-------------------------------------------------------------' write(*,*) -end subroutine print_excitation - - +end subroutine diff --git a/src/LR/print_transition_vectors_pp.f90 b/src/LR/print_transition_vectors_pp.f90 index 4d686fe..a339dcd 100644 --- a/src/LR/print_transition_vectors_pp.f90 +++ b/src/LR/print_transition_vectors_pp.f90 @@ -191,4 +191,4 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip if(nOO > 0) write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:)) write(*,*) -end subroutine print_transition_vectors_pp +end subroutine diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index e36695d..30e7a79 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -178,4 +178,4 @@ subroutine read_methods(doRHF,doUHF,doRMOM,doUMOM,doKS, & close(unit=1) -end subroutine read_methods +end subroutine diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 523e8e7..fbcde29 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -243,4 +243,4 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t close(unit=1) -end subroutine read_options +end subroutine diff --git a/src/RPA/ACFDT_Tmatrix.f90 b/src/RPA/ACFDT_Tmatrix.f90 deleted file mode 100644 index 5953979..0000000 --- a/src/RPA/ACFDT_Tmatrix.f90 +++ /dev/null @@ -1,223 +0,0 @@ -subroutine ACFDT_Tmatrix(exchange_kernel,doXBS,dRPA,TDA_T,TDA,BSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & - nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & - Om2t,X2t,Y2t,rho1t,rho2t,ERI,eT,eGT,EcAC) - -! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem for the T-matrix - - implicit none - include 'parameters.h' - include 'quadrature.h' - -! Input variables - - logical,intent(in) :: doXBS - logical,intent(in) :: exchange_kernel - logical,intent(in) :: dRPA - logical,intent(in) :: TDA_T - logical,intent(in) :: TDA - logical,intent(in) :: BSE - logical,intent(in) :: singlet - logical,intent(in) :: triplet - - 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) :: nOOs - integer,intent(in) :: nOOt - integer,intent(in) :: nVVs - integer,intent(in) :: nVVt - - double precision,intent(in) :: Om1s(nVVs) - double precision,intent(in) :: X1s(nVVs,nVVs) - double precision,intent(in) :: Y1s(nOOs,nVVs) - double precision,intent(in) :: Om2s(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) :: Om1t(nVVt) - double precision,intent(in) :: X1t(nVVt,nVVt) - double precision,intent(in) :: Y1t(nOOt,nVVt) - double precision,intent(in) :: Om2t(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) - -! Local variables - - integer :: ispin - integer :: isp_T - integer :: iblock - integer :: iAC - double precision :: lambda - double precision,allocatable :: Ec(:,:) - - double precision :: EcRPA(nspin) - double precision,allocatable :: TAs(:,:) - double precision,allocatable :: TBs(:,:) - double precision,allocatable :: TAt(:,:) - double precision,allocatable :: TBt(:,:) - double precision,allocatable :: Om(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - -! Output variables - - double precision,intent(out) :: EcAC(nspin) - -! Memory allocation - - allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), & - Om(nS,nspin),XpY(nS,nS,nspin),XmY(nS,nS,nspin)) - allocate(Ec(nAC,nspin)) - -! Antisymmetrized kernel version - - if(exchange_kernel) then - - write(*,*) - write(*,*) '*** Exchange kernel version ***' - write(*,*) - - end if - - EcAC(:) = 0d0 - Ec(:,:) = 0d0 - -! Singlet manifold - - if(singlet) then - - ispin = 1 - - write(*,*) '--------------' - write(*,*) 'Singlet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - if(doXBS) then - - isp_T = 1 - iblock = 3 - - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) - - call GTpp_excitation_density(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,Om1s,rho1s,Om2s,rho2s,TAs) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs) - - isp_T = 2 - iblock = 4 - - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) - - call GTpp_excitation_density(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,Om1t,rho1t,Om2t,rho2t,TAt) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt) - - end if - - call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt+TAs,TBt+TBs, & - EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) - - end do - - EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) - - if(exchange_kernel) EcAC(ispin) = 0.5d0*EcAC(ispin) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - -! Triplet manifold - - if(triplet) then - - ispin = 2 - - write(*,*) '--------------' - write(*,*) 'Triplet states' - write(*,*) '--------------' - write(*,*) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)' - write(*,*) '-----------------------------------------------------------------------------------' - - do iAC=1,nAC - - lambda = rAC(iAC) - - if(doXBS) then - - isp_T = 1 - iblock = 3 - - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,lambda,eT,ERI,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(isp_T)) - - call GTpp_excitation_density(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,Om1s,rho1s,Om2s,rho2s,TAs) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,lambda,Om1s,rho1s,Om2s,rho2s,TBs) - - isp_T = 2 - iblock = 4 - - call ppLR(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,lambda,eT,ERI,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(isp_T)) - - call GTpp_excitation_density(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,Om1t,rho1t,Om2t,rho2t,TAt) - if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,lambda,Om1t,rho1t,Om2t,rho2t,TBt) - - end if - - call linear_response_BSE(ispin,.false.,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,eGT,ERI,TAt-TAs,TBt-TBs, & - EcAC(ispin),Om(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - - call phACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin)) - - write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin) - - end do - - EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin)) - - if(exchange_kernel) EcAC(ispin) = 1.5d0*EcAC(ispin) - - write(*,*) '-----------------------------------------------------------------------------------' - write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin) - write(*,*) '-----------------------------------------------------------------------------------' - write(*,*) - - end if - -end subroutine diff --git a/src/RPA/soRPAx.f90 b/src/RPA/soRPAx.f90 deleted file mode 100644 index 840b711..0000000 --- a/src/RPA/soRPAx.f90 +++ /dev/null @@ -1,113 +0,0 @@ -subroutine soRPAx(TDA,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) - -! Perform random phase approximation calculation with exchange (aka TDHF) in the -! spinorbital basis - - implicit none - include 'parameters.h' - include 'quadrature.h' - -! Input variables - - logical,intent(in) :: TDA - 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) :: 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 - double precision,allocatable :: Omega(:) - double precision,allocatable :: XpY(:,:) - double precision,allocatable :: XmY(:,:) - double precision,allocatable :: X(:,:) - double precision,allocatable :: Y(:,:) - double precision,allocatable :: Xinv(:,:) - double precision,allocatable :: t(:,:,:,:) - - double precision :: EcRPAx - - integer ::i,a,j,b,k,c,ia,jb,kc - -! Hello world - - write(*,*) - write(*,*)'***********************************************************' - write(*,*)'| Random phase approximation calculation with exchange |' - write(*,*)'***********************************************************' - write(*,*) - -! TDA - - if(TDA) then - write(*,*) 'Tamm-Dancoff approximation activated!' - write(*,*) ' => RPAx + TDA = CIS ' - write(*,*) - end if - -! Initialization - - EcRPAx = 0d0 - -! Memory allocation - - allocate(Omega(nS),XpY(nS,nS),XmY(nS,nS)) - - ispin = 3 - - call pphLR(ispin,.false.,TDA,0d0,nBas,nC,nO,nV,nR,nS,1d0,eHF,ERI,EcRPAx,Omega,XpY,XmY) - call print_excitation('soRPAx@HF ',ispin,nS,Omega) - - EcRPAx = 0.5d0*EcRPAx - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@RPAx correlation energy =',EcRPAx - write(*,'(2X,A50,F20.10)') 'Tr@RPAx total energy =',ENuc + ERHF + EcRPAx - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - -! allocate(X(nS,nS),Y(nS,nS),Xinv(nS,nS),t(nO,nO,nV,nV)) - -! X(:,:) = transpose(0.5d0*(XpY(:,:) + XmY(:,:))) -! Y(:,:) = transpose(0.5d0*(XpY(:,:) - XmY(:,:))) - -! call matout(nS,nS,matmul(transpose(X),X)-matmul(transpose(Y),Y)) - -! call inverse_matrix(nS,X,Xinv) - -! t = 0d0 -! ia = 0 -! do i=1,nO -! do a=1,nV -! ia = ia + 1 - -! jb = 0 -! do j=1,nO -! do b=1,nV -! jb = jb + 1 - -! kc = 0 -! do k=1,nO -! do c=1,nV -! kc = kc + 1 - -! t(i,j,a,b) = t(i,j,a,b) + Y(ia,kc)*Xinv(kc,jb) - -! end do -! end do -! end do -! end do -! end do -! end do - -! call matout(nO*nO,nV*nV,t) - -end subroutine soRPAx diff --git a/src/RPA/sort_ppRPA.f90 b/src/RPA/sort_ppRPA.f90 index 993017f..c5e2ca6 100644 --- a/src/RPA/sort_ppRPA.f90 +++ b/src/RPA/sort_ppRPA.f90 @@ -234,4 +234,4 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! call matout(nVV,nOO,X2) ! call matout(nOO,nOO,Y2) -end subroutine sort_ppRPA +end subroutine diff --git a/src/utils/DIIS_extrapolation.f90 b/src/utils/DIIS_extrapolation.f90 index ed5d728..fdf8019 100644 --- a/src/utils/DIIS_extrapolation.f90 +++ b/src/utils/DIIS_extrapolation.f90 @@ -52,4 +52,4 @@ subroutine DIIS_extrapolation(rcond,n_err,n_e,n_diis,error,e,error_in,e_inout) e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis))) -end subroutine DIIS_extrapolation +end subroutine diff --git a/src/utils/chem_to_phys_ERI.f90 b/src/utils/chem_to_phys_ERI.f90 index 9764bcb..34c7084 100644 --- a/src/utils/chem_to_phys_ERI.f90 +++ b/src/utils/chem_to_phys_ERI.f90 @@ -30,4 +30,4 @@ subroutine chem_to_phys_ERI(nBas,ERI) ERI(:,:,:,:) = pERI(:,:,:,:) -end subroutine chem_to_phys_ERI +end subroutine diff --git a/src/utils/level_shifting.f90 b/src/utils/level_shifting.f90 index 9e8e23c..8435b32 100644 --- a/src/utils/level_shifting.f90 +++ b/src/utils/level_shifting.f90 @@ -34,4 +34,4 @@ subroutine level_shifting(level_shift,nBas,nO,S,c,F) Sc(:,:) = matmul(S,c) F(:,:) = matmul(Sc,matmul(F_MO,transpose(Sc))) -end subroutine level_shifting +end subroutine diff --git a/src/utils/read_integrals_sph.f90 b/src/utils/read_integrals_sph.f90 deleted file mode 100644 index 0ee515e..0000000 --- a/src/utils/read_integrals_sph.f90 +++ /dev/null @@ -1,128 +0,0 @@ -subroutine read_integrals_sph(nEl,nBas,S,T,V,Hc,G) - -! Read one- and two-electron integrals from files - - implicit none - include 'parameters.h' - -! Input variables - - integer,intent(in) :: nBas - -! Local variables - - logical :: debug - integer :: nEl(nspin) - integer :: mu,nu,la,si - double precision :: Ov,Kin,Nuc,ERI - double precision :: rs,R,Rinv - -! Output variables - - double precision,intent(out) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas) - -! Open file with integrals - - debug = .false. - - open(unit=1,file='input/sph') - read(1,*) - read(1,*) rs - - R = sqrt(dble(sum(nEl(:))))/2d0*rs - Rinv = 1d0/R - - print*, 'Scaling integrals by ',R - - open(unit=8 ,file='/Users/loos/Integrals/QuAcK_Sph/Ov.dat') - open(unit=9 ,file='/Users/loos/Integrals/QuAcK_Sph/Kin.dat') - open(unit=10,file='/Users/loos/Integrals/QuAcK_Sph/Nuc.dat') - open(unit=11,file='/Users/loos/Integrals/QuAcK_Sph/ERI.dat') - -! Read overlap integrals - - S(:,:) = 0d0 - do - read(8,*,end=8) mu,nu,Ov - S(mu,nu) = Ov - enddo - 8 close(unit=8) - -! Read kinetic integrals - - T(:,:) = 0d0 - do - read(9,*,end=9) mu,nu,Kin - T(mu,nu) = Rinv**2*Kin - enddo - 9 close(unit=9) - -! Read nuclear integrals - - V(:,:) = 0d0 - do - read(10,*,end=10) mu,nu,Nuc - V(mu,nu) = Nuc - enddo - 10 close(unit=10) - -! Define core Hamiltonian - - Hc(:,:) = T(:,:) + V(:,:) - -! Read nuclear integrals - - G(:,:,:,:) = 0d0 - do - read(11,*,end=11) mu,nu,la,si,ERI - - ERI = Rinv*ERI -! <12|34> - G(mu,nu,la,si) = ERI -! <32|14> - G(la,nu,mu,si) = ERI -! <14|32> - G(mu,si,la,nu) = ERI -! <34|12> - G(la,si,mu,nu) = ERI -! <41|23> - G(si,mu,nu,la) = ERI -! <23|41> - G(nu,la,si,mu) = ERI -! <21|43> - G(nu,mu,si,la) = ERI -! <43|21> - G(si,la,nu,mu) = ERI - enddo - 11 close(unit=11) - - -! Print results - if(debug) then - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Overlap integrals' - write(*,'(A28)') '----------------------' - call matout(nBas,nBas,S) - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Kinetic integrals' - write(*,'(A28)') '----------------------' - call matout(nBas,nBas,T) - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Nuclear integrals' - write(*,'(A28)') '----------------------' - call matout(nBas,nBas,V) - write(*,*) - write(*,'(A28)') '----------------------' - write(*,'(A28)') 'Electron repulsion integrals' - write(*,'(A28)') '----------------------' - do la=1,nBas - do si=1,nBas - call matout(nBas,nBas,G(1,1,la,si)) - enddo - enddo - write(*,*) - endif - -end subroutine read_integrals_sph diff --git a/src/utils/rij.f90 b/src/utils/rij.f90 deleted file mode 100644 index ff9b6f8..0000000 --- a/src/utils/rij.f90 +++ /dev/null @@ -1,24 +0,0 @@ -subroutine rij(nWalk,r,r12) - -! Compute the interelectronic distances - - implicit none - -! Input variables - - integer,intent(in) :: nWalk - double precision,intent(in) :: r(nWalk,1:2,1:3) - -! Output variables - - double precision,intent(out) :: r12(nWalk) - -! Compute - - r12(1:nWalk) = (r(1:nWalk,1,1)-r(1:nWalk,2,1))**2 & - + (r(1:nWalk,1,2)-r(1:nWalk,2,2))**2 & - + (r(1:nWalk,1,3)-r(1:nWalk,2,3))**2 - - r12 = sqrt(r12) - -end subroutine rij diff --git a/src/utils/sort.f90 b/src/utils/sort.f90 index 91d4902..87d3239 100644 --- a/src/utils/sort.f90 +++ b/src/utils/sort.f90 @@ -8,7 +8,7 @@ call rec_quicksort(x,iorder,isize,1,isize,1) -end subroutine quick_sort +end subroutine recursive subroutine rec_quicksort(x,iorder,isize,first,last,level) @@ -58,7 +58,7 @@ recursive subroutine rec_quicksort(x,iorder,isize,first,last,level) call rec_quicksort(x, iorder, isize, j+1, last,level/2) endif endif -end subroutine rec_quicksort +end subroutine subroutine set_order(x,iorder,isize,jsize) @@ -86,4 +86,4 @@ subroutine set_order(x,iorder,isize,jsize) deallocate(xtmp) -end subroutine set_order +end subroutine diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index 9192cbc..6832245 100644 --- a/src/utils/utils.f90 +++ b/src/utils/utils.f90 @@ -19,7 +19,7 @@ function Kronecker_delta(i,j) result(delta) delta = 0d0 endif -end function Kronecker_delta +end function function KroneckerDelta(i,j) result(delta) @@ -42,7 +42,7 @@ function KroneckerDelta(i,j) result(delta) delta = 0 endif -end function KroneckerDelta +end function !------------------------------------------------------------------------ subroutine matout(m,n,A) @@ -73,7 +73,7 @@ subroutine matout(m,n,A) enddo enddo -end subroutine matout +end subroutine !------------------------------------------------------------------------ subroutine trace_vector(n,v,Tr) @@ -101,7 +101,7 @@ subroutine trace_vector(n,v,Tr) Tr = Tr + v(i) enddo -end subroutine trace_vector +end subroutine !------------------------------------------------------------------------ function trace_matrix(n,A) result(Tr) @@ -128,7 +128,7 @@ function trace_matrix(n,A) result(Tr) Tr = Tr + A(i,i) enddo -end function trace_matrix +end function !------------------------------------------------------------------------ subroutine compute_error(nData,Mean,Var,Error) @@ -148,7 +148,7 @@ subroutine compute_error(nData,Mean,Var,Error) Error = sqrt((Var-Mean**2/nData)/nData/(nData-1d0)) -end subroutine compute_error +end subroutine !------------------------------------------------------------------------ subroutine identity_matrix(N,A) @@ -175,7 +175,7 @@ subroutine identity_matrix(N,A) A(i,i) = 1d0 enddo -end subroutine identity_matrix +end subroutine !------------------------------------------------------------------------ subroutine prepend(N,M,A,b) @@ -208,7 +208,7 @@ subroutine prepend(N,M,A,b) A(i,1) = b(i) enddo -end subroutine prepend +end subroutine !------------------------------------------------------------------------ subroutine append(N,M,A,b) @@ -237,7 +237,7 @@ subroutine append(N,M,A,b) A(i,M) = b(i) enddo -end subroutine append +end subroutine !------------------------------------------------------------------------ subroutine AtDA(N,A,D,B) @@ -270,7 +270,7 @@ subroutine AtDA(N,A,D,B) enddo enddo -end subroutine AtDA +end subroutine !------------------------------------------------------------------------ subroutine ADAt(N,A,D,B) @@ -303,7 +303,7 @@ subroutine ADAt(N,A,D,B) enddo enddo -end subroutine ADAt +end subroutine !------------------------------------------------------------------------ subroutine DA(N,D,A) @@ -323,7 +323,7 @@ subroutine DA(N,D,A) enddo enddo -end subroutine DA +end subroutine !------------------------------------------------------------------------ subroutine AD(N,A,D) @@ -344,7 +344,7 @@ subroutine AD(N,A,D) enddo enddo -end subroutine AD +end subroutine !------------------------------------------------------------------------ subroutine print_warning(message) @@ -357,7 +357,7 @@ subroutine print_warning(message) write(*,*) message -end subroutine print_warning +end subroutine !------------------------------------------------------------------------ @@ -387,7 +387,7 @@ subroutine CalcTrAB(n,A,B,Tr) enddo enddo -end subroutine CalcTrAB +end subroutine !------------------------------------------------------------------------ @@ -408,7 +408,7 @@ function EpsilonSwitch(i,j) result(delta) delta = -1 endif -end function EpsilonSwitch +end function !------------------------------------------------------------------------ @@ -429,7 +429,7 @@ function KappaCross(i,j,k) result(kappa) kappa = dble(EpsilonSwitch(i,j)*KroneckerDelta(i,k) - EpsilonSwitch(k,i)*KroneckerDelta(i,j)) -end function KappaCross +end function !------------------------------------------------------------------------ @@ -472,7 +472,7 @@ subroutine CalcInv3(a,det) enddo enddo -end subroutine CalcInv3 +end subroutine !------------------------------------------------------------------------ @@ -530,7 +530,7 @@ subroutine CalcInv4(a,det) enddo enddo -end subroutine CalcInv4 +end subroutine subroutine wall_time(t) implicit none @@ -542,5 +542,5 @@ subroutine wall_time(t) endif CALL SYSTEM_CLOCK(count=c) t = dble(c)/dble(rate) -end subroutine wall_time +end subroutine diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 7aa870c..84b8332 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -35,7 +35,7 @@ subroutine diagonalize_general_matrix(N,A,WR,VR) print*,'Problem in diagonalize_general_matrix (dgeev)!!' endif -end subroutine diagonalize_general_matrix +end subroutine subroutine diagonalize_matrix(N,A,e) @@ -72,7 +72,7 @@ subroutine diagonalize_matrix(N,A,e) print*,'Problem in diagonalize_matrix (dsyev)!!' endif -end subroutine diagonalize_matrix +end subroutine subroutine svd(N,A,U,D,Vt) @@ -157,7 +157,7 @@ subroutine inverse_matrix(N,A,B) deallocate(ipiv,work) -end subroutine inverse_matrix +end subroutine subroutine linear_solve(N,A,b,x,rcond) @@ -187,7 +187,7 @@ subroutine linear_solve(N,A,b,x,rcond) ! endif -end subroutine linear_solve +end subroutine subroutine easy_linear_solve(N,A,b,x) @@ -218,5 +218,5 @@ subroutine easy_linear_solve(N,A,b,x) endif -end subroutine easy_linear_solve +end subroutine