From 5ef788cbc93d861f8d4e102cb87c54f2c47f638d Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 24 Mar 2020 12:28:56 +0100 Subject: [PATCH] pCCD and ppRPA --- examples/molecule.LiF | 2 +- examples/molecule.LiH | 2 +- input/basis | 46 +++++++++++++----------- input/methods | 4 +-- input/molecule | 4 +-- input/molecule.xyz | 2 +- input/options | 2 +- input/weight | 46 +++++++++++++----------- src/QuAcK/CCD.f90 | 2 +- src/QuAcK/G0T0.f90 | 3 ++ src/QuAcK/QuAcK.f90 | 7 ++-- src/QuAcK/excitation_density_Tmatrix.f90 | 36 +++++++------------ src/QuAcK/linear_response_C_pp.f90 | 2 +- src/QuAcK/linear_response_D_pp.f90 | 2 +- src/QuAcK/linear_response_pp.f90 | 27 ++++++-------- src/QuAcK/pCCD.f90 | 9 ++++- src/QuAcK/print_excitation.f90 | 2 +- src/QuAcK/self_energy_Tmatrix_diag.f90 | 4 +-- src/QuAcK/soG0T0.f90 | 13 +++---- src/QuAcK/sort_ppRPA.f90 | 14 ++++---- src/eDFT/LIM_RKS.f90 | 2 +- src/eDFT/quadrature_grid.f90 | 2 +- 22 files changed, 117 insertions(+), 116 deletions(-) diff --git a/examples/molecule.LiF b/examples/molecule.LiF index 07213b0..892b553 100644 --- a/examples/molecule.LiF +++ b/examples/molecule.LiF @@ -2,4 +2,4 @@ 2 6 6 0 0 # Znuc x y z Li 0. 0. 0. - F 0. 0. 3.5 + F 0. 0. 2.954767548731 diff --git a/examples/molecule.LiH b/examples/molecule.LiH index 7b1b846..e91dfd2 100644 --- a/examples/molecule.LiH +++ b/examples/molecule.LiH @@ -2,4 +2,4 @@ 2 2 2 0 0 # Znuc x y z Li 0. 0. 0. - H 0. 0. 3.5 + H 0. 0. 3.051018991463 diff --git a/input/basis b/input/basis index 1e6d894..f19a2d0 100644 --- a/input/basis +++ b/input/basis @@ -1,24 +1,30 @@ -1 9 -S 4 - 1 234.0000000 0.0025870 - 2 35.1600000 0.0195330 - 3 7.9890000 0.0909980 - 4 2.2120000 0.2720500 +1 6 +S 8 + 1 17880.0000000 0.0007380 + 2 2683.0000000 0.0056770 + 3 611.5000000 0.0288830 + 4 173.5000000 0.1085400 + 5 56.6400000 0.2909070 + 6 20.4200000 0.4483240 + 7 7.8100000 0.2580260 + 8 1.6530000 0.0150630 +S 8 + 1 17880.0000000 -0.0001720 + 2 2683.0000000 -0.0013570 + 3 611.5000000 -0.0067370 + 4 173.5000000 -0.0276630 + 5 56.6400000 -0.0762080 + 6 20.4200000 -0.1752270 + 7 7.8100000 -0.1070380 + 8 1.6530000 0.5670500 S 1 - 1 0.6669000 1.0000000 -S 1 - 1 0.2089000 1.0000000 -S 1 - 1 0.0513800 1.0000000 + 1 0.4869000 1.0000000 +P 3 + 1 28.3900000 0.0460870 + 2 6.2700000 0.2401810 + 3 1.6950000 0.5087440 P 1 - 1 3.0440000 1.0000000 -P 1 - 1 0.7580000 1.0000000 -P 1 - 1 0.1993000 1.0000000 + 1 0.4317000 1.0000000 D 1 - 1 1.9650000 1.0000000 -D 1 - 1 0.4592000 1.0000000 - + 1 2.2020000 1.0000000 diff --git a/input/methods b/input/methods index e126029..7895c92 100644 --- a/input/methods +++ b/input/methods @@ -3,7 +3,7 @@ # MP2 MP3 MP2-F12 T F F # CCD CCSD CCSD(T) drCCD rCCD lCCD pCCD - F F F F F F F + F F F F F F T # CIS RPA RPAx ppRPA ADC F F F F F # G0F2 evGF2 G0F3 evGF3 @@ -11,6 +11,6 @@ # G0W0 evGW qsGW F F F # G0T0 evGT qsGT - T F F + F F F # MCMP2 F diff --git a/input/molecule b/input/molecule index c78e87e..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 5 5 0 0 # Znuc x y z - He 0.0 0.0 0.0 + Ne 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 797b5fc..1c70680 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,3 @@ 1 - He 0.0000000000 0.0000000000 0.0000000000 + Ne 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/options b/input/options index 7881bfa..2de160f 100644 --- a/input/options +++ b/input/options @@ -9,7 +9,7 @@ # GF: maxSCF thresh DIIS n_diis lin renorm 256 0.00001 T 5 T 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta - 256 0.00001 T 5 F F T F F F T 0.000 + 256 0.00001 T 5 F F F F F F F 0.000 # ACFDT: AC Kx XBS T F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/input/weight b/input/weight index 1e6d894..f19a2d0 100644 --- a/input/weight +++ b/input/weight @@ -1,24 +1,30 @@ -1 9 -S 4 - 1 234.0000000 0.0025870 - 2 35.1600000 0.0195330 - 3 7.9890000 0.0909980 - 4 2.2120000 0.2720500 +1 6 +S 8 + 1 17880.0000000 0.0007380 + 2 2683.0000000 0.0056770 + 3 611.5000000 0.0288830 + 4 173.5000000 0.1085400 + 5 56.6400000 0.2909070 + 6 20.4200000 0.4483240 + 7 7.8100000 0.2580260 + 8 1.6530000 0.0150630 +S 8 + 1 17880.0000000 -0.0001720 + 2 2683.0000000 -0.0013570 + 3 611.5000000 -0.0067370 + 4 173.5000000 -0.0276630 + 5 56.6400000 -0.0762080 + 6 20.4200000 -0.1752270 + 7 7.8100000 -0.1070380 + 8 1.6530000 0.5670500 S 1 - 1 0.6669000 1.0000000 -S 1 - 1 0.2089000 1.0000000 -S 1 - 1 0.0513800 1.0000000 + 1 0.4869000 1.0000000 +P 3 + 1 28.3900000 0.0460870 + 2 6.2700000 0.2401810 + 3 1.6950000 0.5087440 P 1 - 1 3.0440000 1.0000000 -P 1 - 1 0.7580000 1.0000000 -P 1 - 1 0.1993000 1.0000000 + 1 0.4317000 1.0000000 D 1 - 1 1.9650000 1.0000000 -D 1 - 1 0.4592000 1.0000000 - + 1 2.2020000 1.0000000 diff --git a/src/QuAcK/CCD.f90 b/src/QuAcK/CCD.f90 index 20ba3a0..3929415 100644 --- a/src/QuAcK/CCD.f90 +++ b/src/QuAcK/CCD.f90 @@ -48,7 +48,7 @@ subroutine CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) double precision,allocatable :: r2(:,:,:,:) double precision,allocatable :: t2(:,:,:,:) - integer :: n_diis + integer :: n_diis,i,j,a,b double precision :: rcond double precision,allocatable :: error_diis(:,:) double precision,allocatable :: t_diis(:,:) diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 8c43f4b..b8b1dd7 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -137,6 +137,9 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute T-matrix version of the self-energy !---------------------------------------------- + rho1s = 0d0 + rho1t = 0d0 + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 0374372..7a8349f 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -367,6 +367,7 @@ program QuAcK call cpu_time(start_CCD) call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) + call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for CCD = ',t_CCD,' seconds' @@ -669,10 +670,10 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) -! call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, & -! singlet_manifold,triplet_manifold,linGW,eta, & -! nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF) call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(1),ENuc,ERHF,ERI_MO_basis,eHF) + ! call G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA, & + ! singlet_manifold,triplet_manifold,linGW,eta, & + ! nBas,nC(1),nO(1),nV(1),nR(1),nS(1),ENuc,ERHF,ERI_MO_basis,eHF) call cpu_time(end_G0T0) diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index c1de99e..5b176a7 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -46,9 +46,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do d=c,nBas-nR cd = cd + 1 rho1(p,i,ab) = rho1(p,i,ab) & -! + ERI(p,i,c,d)*X1(cd,ab) - + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) & - /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) + + ERI(p,i,c,d)*X1(cd,ab) +! + (ERI(p,i,c,d) + ERI(p,i,d,c))*X1(cd,ab) & +! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -58,9 +58,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do l=k,nO kl = kl + 1 rho1(p,i,ab) = rho1(p,i,ab) & -! + ERI(p,i,k,l)*Y1(kl,ab) - + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) & - /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) + + ERI(p,i,k,l)*Y1(kl,ab) +! + (ERI(p,i,k,l) + ERI(p,i,l,k))*Y1(kl,ab) & +! /sqrt((1d0 + Kronecker_delta(p,i))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -76,9 +76,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do d=c,nBas-nR cd = cd + 1 rho2(p,a,ij) = rho2(p,a,ij) & -! + ERI(p,nO+a,c,d)*X2(cd,ij) - + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) & - /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d))) + + ERI(p,nO+a,c,d)*X2(cd,ij) +! + (ERI(p,nO+a,c,d) + ERI(p,nO+a,d,c))*X2(cd,ij) & +! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(c,d))) end do end do @@ -88,9 +88,9 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do l=k,nO kl = kl + 1 rho2(p,a,ij) = rho2(p,a,ij) & -! + ERI(p,nO+a,k,l)*Y2(kl,ij) - + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) & - /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l))) + + ERI(p,nO+a,k,l)*Y2(kl,ij) +! + (ERI(p,nO+a,k,l) + ERI(p,nO+a,l,k))*Y2(kl,ij) & +! /sqrt((1d0 + Kronecker_delta(p,nO+a))*(1d0 + Kronecker_delta(k,l))) end do end do @@ -118,7 +118,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r cd = cd + 1 rho1(p,i,ab) = rho1(p,i,ab) & + (ERI(p,i,c,d) - ERI(p,i,d,c))*X1(cd,ab) - print*,rho1(p,i,ab),ERI(p,i,c,d),X1(cd,ab) end do end do @@ -129,7 +128,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r kl = kl + 1 rho1(p,i,ab) = rho1(p,i,ab) & + (ERI(p,i,k,l) - ERI(p,i,l,k))*Y1(kl,ab) - print*,rho1(p,i,ab),ERI(p,i,k,l),Y1(kl,ab) end do end do @@ -222,16 +220,6 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r end do -! do p=nC+1,nBas-nR -! do i=nC+1,nO -! do ab=1,nVV -! print*,p,i,ab,rho1(p,i,ab) -! end do -! end do -! end do - -! call matout(nVV,nOO,Y1) - end if end subroutine excitation_density_Tmatrix diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90 index 9c93daa..ae7d399 100644 --- a/src/QuAcK/linear_response_C_pp.f90 +++ b/src/QuAcK/linear_response_C_pp.f90 @@ -25,7 +25,7 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) ! Define the chemical potential eF = e(nO) + e(nO+1) -! eF = 0d0 + eF = 0d0 ! Build C matrix for the singlet manifold diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90 index 93167af..002a6bc 100644 --- a/src/QuAcK/linear_response_D_pp.f90 +++ b/src/QuAcK/linear_response_D_pp.f90 @@ -25,7 +25,7 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) ! Define the chemical potential eF = e(nO) + e(nO+1) -! eF = 0d0 + eF = 0d0 ! Build the D matrix for the singlet manifold diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index f178b7f..8f35f2a 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -1,5 +1,5 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & - e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcppRPA) + e,ERI,Omega1,X1,Y1,Omega2,X2,Y2,EcRPA) ! Compute the p-p channel of the linear response: see Scuseria et al. JCP 139, 104113 (2013) @@ -22,6 +22,8 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & integer :: ab,cd,ij,kl integer :: p,q,r,s double precision :: trace_matrix + double precision :: EcRPA1 + double precision :: EcRPA2 double precision,allocatable :: B(:,:) double precision,allocatable :: C(:,:) double precision,allocatable :: D(:,:) @@ -37,7 +39,7 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & double precision,intent(out) :: Omega2(nOO) double precision,intent(out) :: X2(nVV,nOO) double precision,intent(out) :: Y2(nOO,nOO) - double precision,intent(out) :: EcppRPA + double precision,intent(out) :: EcRPA ! Memory allocation @@ -117,10 +119,11 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & ! print*, 'pp-RPA matrix' ! call matout(nOO+nVV,nOO+nVV,M(:,:)) +! write(*,*) ! Diagonalize the p-h matrix - call diagonalize_general_matrix(nOO+nVV,M,Omega,Z) + if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Omega,Z) ! write(*,*) 'pp-RPA excitation energies' ! call matout(nOO+nVV,1,Omega(:)) @@ -132,11 +135,11 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & ! Compute the RPA correlation energy -! EcppRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) -! print*,+sum(Omega1(:)),- trace_matrix(nVV,C(:,:)) -! print*,-sum(Omega2(:)),- trace_matrix(nOO,D(:,:)) - EcppRPA = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:)) -! EcppRPA = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:)) + EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) ) + EcRPA1 = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:)) + EcRPA2 = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:)) + if(abs(EcRPA - EcRPA1) > 1d-10 .or. abs(EcRPA - EcRPA2) > 1d-10) & + print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' ! write(*,*)'X1' ! call matout(nVV,nVV,X1) @@ -147,12 +150,4 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, & ! write(*,*)'Y2' ! call matout(nOO,nOO,Y2) -! print*,'Ec(pp-RPA) = ',EcppRPA - -! print*,'Eigenvalues' -! call matout(nOO+nVV,1,Omega) - -! print*,'Eigenvectors' -! call matout(nOO+nVV,nOO+nVV,matmul(transpose(Z),Z)) - end subroutine linear_response_pp diff --git a/src/QuAcK/pCCD.f90 b/src/QuAcK/pCCD.f90 index 49e44f7..dcfe7ab 100644 --- a/src/QuAcK/pCCD.f90 +++ b/src/QuAcK/pCCD.f90 @@ -53,6 +53,8 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) allocate(delta_OOVV(nO,nV)) + delta_OOVV(:,:) = 0d0 + do i=nC+1,nO do a=1,nV-nR delta_OOVV(i,a) = 2d0*(eHF(nO+a) - eHF(i)) @@ -63,6 +65,11 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) allocate(OOOO(nO,nO),OOVV(nO,nV),OVOV(nO,nV),OVVO(nO,nV),VVVV(nV,nV)) + OOOO(:,:) = 0d0 + OOVV(:,:) = 0d0 + OVVO(:,:) = 0d0 + VVVV(:,:) = 0d0 + do i=nC+1,nO do j=nC+1,nO OOOO(i,j) = ERI(i,i,j,j) @@ -141,7 +148,7 @@ subroutine pCCD(maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) - 2d0*(2d0*OVOV(i,a) - OVVO(i,a) - OOVV(i,a)*t(i,a))*t(i,a) do j=nC+1,nO - r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + y(i,j)*t(j,a) + r(i,a) = r(i,a) - 2d0*OOVV(j,a)*t(j,a)*t(i,a) + OOOO(j,i)*t(j,a) + y(i,j)*t(j,a) end do do b=1,nV-nR diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index a32831e..fca9912 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -7,7 +7,7 @@ subroutine print_excitation(method,ispin,nS,Omega) ! Input variables - character*6,intent(in) :: method + character*12,intent(in) :: method integer,intent(in) :: ispin,nS double precision,intent(in) :: Omega(nS) diff --git a/src/QuAcK/self_energy_Tmatrix_diag.f90 b/src/QuAcK/self_energy_Tmatrix_diag.f90 index 5e509ba..f1a3e00 100644 --- a/src/QuAcK/self_energy_Tmatrix_diag.f90 +++ b/src/QuAcK/self_energy_Tmatrix_diag.f90 @@ -72,7 +72,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do i=nC+1,nO do cd=1,nVVt eps = e(p) + e(i) - Omega1t(cd) - SigT(p) = SigT(p) + rho1t(p,i,cd)*rho1t(p,i,cd)/eps + ! SigT(p) = SigT(p) + rho1t(p,i,cd)*rho1t(p,i,cd)/eps enddo enddo enddo @@ -83,7 +83,7 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, do a=1,nV-nR do kl=1,nOOt eps = e(p) + e(nO+a) - Omega2t(kl) - SigT(p) = SigT(p) + rho2t(p,a,kl)*rho2t(p,a,kl)/eps + ! SigT(p) = SigT(p) + rho2t(p,a,kl)*rho2t(p,a,kl)/eps enddo enddo enddo diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index 8fc0f0b..47e0998 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -89,31 +89,28 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) ! Compute excitation densities for the T-matrix call excitation_density_Tmatrix(ispin,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,sERI(:,:,:,:), & - X1(:,:),Y1(:,:),rho1(:,:,:), & - X2(:,:),Y2(:,:),rho2(:,:,:)) + X1(:,:),Y1(:,:),rho1(:,:,:),X2(:,:),Y2(:,:),rho2(:,:,:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- call self_energy_Tmatrix_diag_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & - Omega1(:),rho1(:,:,:), & - Omega2(:),rho2(:,:,:), & + Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & SigT(:)) ! Compute renormalization factor for T-matrix self-energy call renormalization_factor_Tmatrix_so(eta,nBas2,nC2,nO2,nV2,nR2,nOO,nVV,seHF(:), & - Omega1(:),rho1(:,:,:), & - Omega2(:),rho2(:,:,:), & + Omega1(:),rho1(:,:,:),Omega2(:),rho2(:,:,:), & Z(:)) !---------------------------------------------- ! Solve the quasi-particle equation !---------------------------------------------- - eG0T0(:) = seHF(:) + SigT(:) -! eG0T0(:) = seHF(:) + Z(:)*SigT(:) +! eG0T0(:) = seHF(:) + SigT(:) + eG0T0(:) = seHF(:) + Z(:)*SigT(:) !---------------------------------------------- ! Dump results diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 index c3e7162..a59df23 100644 --- a/src/QuAcK/sort_ppRPA.f90 +++ b/src/QuAcK/sort_ppRPA.f90 @@ -74,16 +74,12 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ab = ab + 1 Omega1(ab) = Omega(pq) Z1(1:nOO+nVV,ab) = Z(1:nOO+nVV,pq) -! X1(1:nVV,ab) = Z( 1: nVV,pq) -! Y1(1:nOO,ab) = Z(nVV+1:nOO+nVV,pq) else ij = ij + 1 Omega2(ij) = Omega(pq) Z2(1:nOO+nVV,ij) = Z(1:nOO+nVV,pq) -! X2(1:nVV,ij) = Z( 1: nVV,pq) -! Y2(1:nOO,ij) = Z(nVV+1:nOO+nVV,pq) end if @@ -115,11 +111,13 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) end if - X1(1:nVV,1:nVV) = Z1( 1: nVV,1:nVV) - Y1(1:nOO,1:nVV) = Z1(nVV+1:nOO+nVV,1:nVV) +! Define submatrices X1, Y1, X2, & Y2 - X2(1:nVV,1:nOO) = Z2( 1: nVV,1:nOO) - Y2(1:nOO,1:nOO) = Z2(nVV+1:nOO+nVV,1:nOO) + X1(1:nVV,1:nVV) = + Z1( 1: nVV,1:nVV) + Y1(1:nOO,1:nVV) = - Z1(nVV+1:nOO+nVV,1:nVV) + + X2(1:nVV,1:nOO) = + Z2( 1: nVV,1:nOO) + Y2(1:nOO,1:nOO) = - Z2(nVV+1:nOO+nVV,1:nOO) ! write(*,*) 'Z1t.M.Z1' ! call matout(nVV,nVV,matmul(matmul(transpose(Z1),M),Z1)) diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 index fb65711..6158264 100644 --- a/src/eDFT/LIM_RKS.f90 +++ b/src/eDFT/LIM_RKS.f90 @@ -73,7 +73,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nGrid,weight,maxSCF,thres !------------------------------------------------------------------------ write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' ZERO-WEIGHT CALCULATION ' + write(*,'(A40)') ' EQUI-WEIGHT CALCULATION ' write(*,'(A40)') '*************************************************' wLIM(1:nEns) = 1d0/dble(nEns) diff --git a/src/eDFT/quadrature_grid.f90 b/src/eDFT/quadrature_grid.f90 index 437b5f7..c1f12c2 100644 --- a/src/eDFT/quadrature_grid.f90 +++ b/src/eDFT/quadrature_grid.f90 @@ -66,7 +66,7 @@ subroutine quadrature_grid(nRad,nAng,nGrid,root,weight) enddo enddo -! Compute values of the basis functions (and the its gradient if required) at each grid point +! Print grids 20 format(T2,58('-')) 30 format(T20,'Radial Quadrature',/, &