diff --git a/src/GT/RG0T0pp.f90 b/src/GT/RG0T0pp.f90 index 07ab335..aaab3ff 100644 --- a/src/GT/RG0T0pp.f90 +++ b/src/GT/RG0T0pp.f90 @@ -1,5 +1,5 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,dTDA,doppBSE,singlet,triplet, & - linearize,eta,regularize,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) + linearize,eta,regularize,nOrb,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,dipole_int,eHF) ! Perform one-shot calculation with a T-matrix self-energy (G0T0) @@ -25,7 +25,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d double precision,intent(in) :: eta logical,intent(in) :: regularize - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -33,9 +33,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d 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) - double precision,intent(in) :: dipole_int(nBas,nBas,ncart) + double precision,intent(in) :: eHF(nOrb) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision,intent(in) :: dipole_int(nOrb,nOrb,ncart) ! Local variables @@ -45,6 +45,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d integer :: iblock integer :: nOOs,nOOt integer :: nVVs,nVVt + integer :: n_states, n_states_diag double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcGM @@ -63,6 +64,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d double precision,allocatable :: Z(:) double precision,allocatable :: eGT(:) double precision,allocatable :: eGTlin(:) + double precision, allocatable :: Om(:), R(:,:) ! Output variables @@ -92,8 +94,8 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! Dimensions of the pp-RPA linear reponse matrices -! nOOs = nO*(nO + 1)/2 -! nVVs = nV*(nV + 1)/2 + !nOOs = nO*(nO + 1)/2 + !nVVs = nV*(nV + 1)/2 nOOs = nO*nO nVVs = nV*nV @@ -105,30 +107,39 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Om1s(nVVs),X1s(nVVs,nVVs),Y1s(nOOs,nVVs), & Om2s(nOOs),X2s(nVVs,nOOs),Y2s(nOOs,nOOs), & - rho1s(nBas,nBas,nVVs),rho2s(nBas,nBas,nOOs), & + rho1s(nOrb,nOrb,nVVs),rho2s(nOrb,nOrb,nOOs), & Om1t(nVVt),X1t(nVVt,nVVt),Y1t(nOOt,nVVt), & Om2t(nOOt),X2t(nVVt,nOOt),Y2t(nOOt,nOOt), & - rho1t(nBas,nBas,nVVt),rho2t(nBas,nBas,nOOt), & - Sig(nBas),Z(nBas),eGT(nBas),eGTlin(nBas)) + rho1t(nOrb,nOrb,nVVt),rho2t(nOrb,nOrb,nOOt), & + Sig(nOrb),Z(nOrb),eGT(nOrb),eGTlin(nOrb)) !---------------------------------------------- ! alpha-beta block !---------------------------------------------- ispin = 1 -! iblock = 1 + !iblock = 1 iblock = 3 ! Compute linear response allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eHF,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eHF,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) + !n_states = nOOs + 5 + !n_states_diag = n_states + 4 + !allocate(Om(nOOs+nVVs), R(nOOs+nVVs,n_states_diag)) + !call ppLR_RG0T0pp_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOs, nVVs, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag) + !print*, 'LAPACK:' + !print*, Om2s + !print*, Om1s + !stop + deallocate(Bpp,Cpp,Dpp) if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-beta)',nVVs,Om1s) @@ -146,12 +157,21 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eHF,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eHF,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) + n_states = nOOt + 5 + n_states_diag = n_states + 4 + allocate(Om(nOOt+nVVt), R(nOOt+nVVt,n_states_diag)) + call ppLR_RG0T0pp_davidson(iblock, TDA_T, nC, nO, nR, nOrb, nOOt, nVVt, 1.d0, eHF, 0.d0, ERI, Om, R, n_states, n_states_diag) + print*, 'LAPACK:' + print*, Om2t + print*, Om1t + stop + deallocate(Bpp,Cpp,Dpp) if(print_T) call print_excitation_energies('ppRPA@RHF','2p (alpha-alpha)',nVVt,Om1t) @@ -164,23 +184,23 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d ! iblock = 1 iblock = 3 - call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) + call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) ! iblock = 2 iblock = 4 - call RGTpp_excitation_density(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) + call RGTpp_excitation_density(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- if(regularize) then - call GTpp_regularization(nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Om1s,rho1s,Om2s,rho2s) - call GTpp_regularization(nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Om1t,rho1t,Om2t,rho2t) + call GTpp_regularization(nOrb,nC,nO,nV,nR,nOOs,nVVs,eHF,Om1s,rho1s,Om2s,rho2s) + call GTpp_regularization(nOrb,nC,nO,nV,nR,nOOt,nVVt,eHF,Om1t,rho1t,Om2t,rho2t) end if - call RGTpp_self_energy_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & + call RGTpp_self_energy_diag(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & Om1t,rho1t,Om2t,rho2t,EcGM,Sig,Z) !---------------------------------------------- @@ -201,12 +221,12 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d write(*,*) ' *** Quasiparticle energies obtained by root search *** ' write(*,*) - call RGTpp_QP_graph(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & + call RGTpp_QP_graph(eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,Om1s,rho1s,Om2s,rho2s, & Om1t,rho1t,Om2t,rho2t,eGTlin,eHF,eGT,Z) end if -! call RGTpp_plot_self_energy(nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, & +! call RGTpp_plot_self_energy(nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF,eGT,Om1s,rho1s,Om2s,rho2s, & ! Om1t,rho1t,Om2t,rho2t) !---------------------------------------------- @@ -221,9 +241,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Bpp(nVVs,nOOs),Cpp(nVVs,nVVs),Dpp(nOOs,nOOs)) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVs,1d0,eGT,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOs,1d0,eGT,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOs,nVVs,1d0,ERI,Bpp) call ppLR(TDA_T,nOOs,nVVs,Bpp,Cpp,Dpp,Om1s,X1s,Y1s,Om2s,X2s,Y2s,EcRPA(ispin)) @@ -235,9 +255,9 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d allocate(Bpp(nVVt,nOOt),Cpp(nVVt,nVVt),Dpp(nOOt,nOOt)) - call ppLR_C(iblock,nBas,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) - call ppLR_D(iblock,nBas,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) - if(.not.TDA_T) call ppLR_B(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) + call ppLR_C(iblock,nOrb,nC,nO,nV,nR,nVVt,1d0,eGT,ERI,Cpp) + call ppLR_D(iblock,nOrb,nC,nO,nV,nR,nOOt,1d0,eGT,ERI,Dpp) + if(.not.TDA_T) call ppLR_B(iblock,nOrb,nC,nO,nV,nR,nOOt,nVVt,1d0,ERI,Bpp) call ppLR(TDA_T,nOOt,nVVt,Bpp,Cpp,Dpp,Om1t,X1t,Y1t,Om2t,X2t,Y2t,EcRPA(ispin)) @@ -246,13 +266,13 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) - call print_RG0T0pp(nBas,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA) + call print_RG0T0pp(nOrb,nO,eHF,ENuc,ERHF,Sig,Z,eGT,EcGM,EcRPA) ! Perform BSE calculation if(dophBSE) then - call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + call RGTpp_phBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI,dipole_int,eHF,eGT,EcBSE) @@ -288,7 +308,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d end if - call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS, & + call RGTpp_phACFDT(exchange_kernel,doXBS,.false.,TDA_T,TDA,dophBSE,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nS, & nOOs,nVVs,nOOt,nVVt,Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t, & Om2t,X2t,Y2t,rho1t,rho2t,ERI,eHF,eGT,EcBSE) @@ -314,7 +334,7 @@ subroutine RG0T0pp(dotest,doACFDT,exchange_kernel,doXBS,dophBSE,TDA_T,TDA,dBSE,d if(doppBSE) then - call RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, & + call RGTpp_ppBSE(TDA_T,TDA,dBSE,dTDA,singlet,triplet,eta,nOrb,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt, & Om1s,X1s,Y1s,Om2s,X2s,Y2s,rho1s,rho2s,Om1t,X1t,Y1t,Om2t,X2t,Y2t,rho1t,rho2t, & ERI,dipole_int,eHF,eGT,EcBSE) diff --git a/src/GW/RGW_ppBSE.f90 b/src/GW/RGW_ppBSE.f90 index 0df6c35..85b6620 100644 --- a/src/GW/RGW_ppBSE.f90 +++ b/src/GW/RGW_ppBSE.f90 @@ -125,6 +125,9 @@ subroutine RGW_ppBSE(TDA_W,TDA,dBSE,dTDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS call ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + ! TODO + !call ppLR_RGW_ppBSE_davidson(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcBSE(ispin)) + call ppLR_transition_vectors(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Om1,X1,Y1,Om2,X2,Y2) !----------------------------------------------------! diff --git a/src/LR/ppLR_B.f90 b/src/LR/ppLR_B.f90 index 0bff801..a4f1fab 100644 --- a/src/LR/ppLR_B.f90 +++ b/src/LR/ppLR_B.f90 @@ -1,4 +1,4 @@ -subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) +subroutine ppLR_B(ispin,nOrb,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) ! Compute the B matrix of the pp channel @@ -8,7 +8,7 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV @@ -16,7 +16,7 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) ! Local variables @@ -32,8 +32,8 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) if(ispin == 1) then ab = 0 - do a=nO+1,nBas-nR - do b=a,nBas-nR + do a=nO+1,nOrb-nR + do b=a,nOrb-nR ab = ab + 1 ij = 0 do i=nC+1,nO @@ -54,8 +54,8 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) if(ispin == 2 .or. ispin == 4) then ab = 0 - do a=nO+1,nBas-nR - do b=a+1,nBas-nR + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR ab = ab + 1 ij = 0 do i=nC+1,nO @@ -76,8 +76,8 @@ subroutine ppLR_B(ispin,nBas,nC,nO,nV,nR,nOO,nVV,lambda,ERI,Bpp) if(ispin == 3) then ab = 0 - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR ab = ab + 1 ij = 0 do i=nC+1,nO diff --git a/src/LR/ppLR_C.f90 b/src/LR/ppLR_C.f90 index 72f15f0..6170797 100644 --- a/src/LR/ppLR_C.f90 +++ b/src/LR/ppLR_C.f90 @@ -1,4 +1,4 @@ -subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) +subroutine ppLR_C(ispin,nOrb,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) ! Compute the C matrix of the pp channel @@ -8,14 +8,14 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nVV double precision,intent(in) :: lambda - double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: e(nOrb),ERI(nOrb,nOrb,nOrb,nOrb) ! Local variables @@ -39,15 +39,15 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) if(ispin == 1) then - a0 = nBas - nR - nO + a0 = nOrb - nR - nO !$OMP PARALLEL DEFAULT(NONE) & !$OMP PRIVATE(a, b, aa, ab, c, d, cd, e_ab, tmp_ab, delta_ac, tmp_cd) & - !$OMP SHARED(nO, nBas, nR, a0, eF, lambda, e, ERI, Cpp) + !$OMP SHARED(nO, nOrb, nR, a0, eF, lambda, e, ERI, Cpp) !$OMP DO - do a = nO+1, nBas-nR + do a = nO+1, nOrb-nR aa = a0 * (a - nO - 1) - (a - nO - 1) * (a - nO) / 2 - nO - do b = a, nBas-nR + do b = a, nOrb-nR ab = aa + b e_ab = e(a) + e(b) - eF @@ -58,14 +58,14 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) endif cd = 0 - do c = nO+1, nBas-nR + do c = nO+1, nOrb-nR delta_ac = 0.d0 if(a .eq. c) then delta_ac = 1.d0 endif - do d = c, nBas-nR + do d = c, nOrb-nR cd = cd + 1 tmp_cd = tmp_ab @@ -88,12 +88,12 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) ! ab = 0 -! do a=nO+1,nBas-nR -! do b=a,nBas-nR +! do a=nO+1,nOrb-nR +! do b=a,nOrb-nR ! ab = ab + 1 ! cd = 0 -! do c=nO+1,nBas-nR -! do d=c,nBas-nR +! do c=nO+1,nOrb-nR +! do d=c,nOrb-nR ! cd = cd + 1 ! ! Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & @@ -110,16 +110,16 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) if(ispin == 2 .or. ispin == 4) then !$OMP PARALLEL & - !$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nBas,nR) & + !$OMP SHARED(Cpp,lambda,ERI,e,eF,nC,nO,nOrb,nR) & !$OMP PRIVATE(c,d,a,b,ab,cd) & !$OMP DEFAULT(NONE) !$OMP DO - do c=nO+1,nBas-nR - do d=c+1,nBas-nR - cd = (c-(nO+1))*(nBas-nR-(nO+1)) - (c-1-(nO+1))*(c-(nO+1))/2 + d - c - do a=nO+1,nBas-nR - do b=a+1,nBas-nR - ab = (a-(nO+1))*(nBas-nR-(nO+1)) - (a-1-(nO+1))*(a-(nO+1))/2 + b - a + do c=nO+1,nOrb-nR + do d=c+1,nOrb-nR + cd = (c-(nO+1))*(nOrb-nR-(nO+1)) - (c-1-(nO+1))*(c-(nO+1))/2 + d - c + do a=nO+1,nOrb-nR + do b=a+1,nOrb-nR + ab = (a-(nO+1))*(nOrb-nR-(nO+1)) - (a-1-(nO+1))*(a-(nO+1))/2 + b - a Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + lambda*(ERI(a,b,c,d) - ERI(a,b,d,c)) @@ -138,12 +138,12 @@ subroutine ppLR_C(ispin,nBas,nC,nO,nV,nR,nVV,lambda,e,ERI,Cpp) if(ispin == 3) then ab = 0 - do a=nO+1,nBas-nR - do b=nO+1,nBas-nR + do a=nO+1,nOrb-nR + do b=nO+1,nOrb-nR ab = ab + 1 cd = 0 - do c=nO+1,nBas-nR - do d=nO+1,nBas-nR + do c=nO+1,nOrb-nR + do d=nO+1,nOrb-nR cd = cd + 1 Cpp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & diff --git a/src/LR/ppLR_D.f90 b/src/LR/ppLR_D.f90 index 4d19ea5..f55b66c 100644 --- a/src/LR/ppLR_D.f90 +++ b/src/LR/ppLR_D.f90 @@ -1,4 +1,4 @@ -subroutine ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) +subroutine ppLR_D(ispin,nOrb,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) ! Compute the D matrix of the pp channel @@ -8,14 +8,14 @@ subroutine ppLR_D(ispin,nBas,nC,nO,nV,nR,nOO,lambda,e,ERI,Dpp) ! Input variables integer,intent(in) :: ispin - integer,intent(in) :: nBas + integer,intent(in) :: nOrb integer,intent(in) :: nC integer,intent(in) :: nO integer,intent(in) :: nV integer,intent(in) :: nR integer,intent(in) :: nOO double precision,intent(in) :: lambda - double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: e(nOrb),ERI(nOrb,nOrb,nOrb,nOrb) ! Local variables diff --git a/src/LR/ppLR_RG0T0pp_davidson.f90 b/src/LR/ppLR_RG0T0pp_davidson.f90 new file mode 100644 index 0000000..6cf336a --- /dev/null +++ b/src/LR/ppLR_RG0T0pp_davidson.f90 @@ -0,0 +1,982 @@ + +! --- + +subroutine ppLR_RG0T0pp_davidson(ispin, TDA, nC, nO, nR, nOrb, nOO, nVV, lambda, e, eF, ERI, Om, R, n_states, n_states_diag) + + ! + ! Extract the low n_states + ! Om(i) (eigenvalues) and + ! R(:,i) (right-eigenvectors) + ! of the pp-RPA matrix + ! + ! (+C +B) + ! ( ) + ! (-B.T -D) + ! + + implicit none + + logical, intent(in) :: TDA + integer, intent(in) :: ispin + integer, intent(in) :: nC, nO, nR, nOrb, nOO, nVV + integer, intent(in) :: n_states ! nb of physical states + integer, intent(in) :: n_states_diag ! nb of states used to get n_states + double precision, intent(in) :: lambda, eF + double precision, intent(in) :: e(nOrb) + double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision, intent(out) :: Om(n_states) + double precision, intent(out) :: R(nOO+nVV,n_states_diag) + + integer :: N, M + integer :: iter, itermax, itertot + integer :: shift1, shift2 + integer :: i, j, ij, k, l, kl, a, b, ab, c, d, cd + integer :: i_omax(n_states) + logical :: converged + character*(16384) :: write_buffer + double precision :: r1, r2, dtwo_pi + double precision :: lambda_tmp + double precision :: to_print(2,n_states) + double precision :: mem + double precision, allocatable :: H_diag(:) + double precision, allocatable :: W(:,:) + double precision, allocatable :: U(:,:) + double precision, allocatable :: h(:,:), h_vec(:,:), h_val(:) + double precision, allocatable :: residual_norm(:) + double precision, allocatable :: overlap(:) + double precision, allocatable :: S_check(:,:) + + double precision, external :: u_dot_u + + dtwo_pi = 6.283185307179586d0 + + N = nOO + nVV + itermax = 8 + M = n_states_diag * itermax + + if(M .ge. N) then + print*, 'N = ', N + print*, 'M = ', M + print*, ' use Lapack or decrease n_states and/or itermax ' + stop + endif + + write(6,'(A)') '' + write(6,'(A)') 'Davidson Diagonalization' + write(6,'(A)') '------------------------' + write(6,'(A)') '' + + write(*,'(A40, I12)') 'Number of states = ', n_states + write(*,'(A40, I12)') 'Number of states in diagonalization = ', n_states_diag + write(*,'(A40, I12)') 'Number of basis functions = ', N + + + + + allocate(H_diag(N)) + allocate(U(N,M)) + allocate(W(N,M)) + allocate(h(M,M), h_vec(M,M), h_val(M)) + allocate(overlap(n_states_diag)) + allocate(residual_norm(n_states_diag)) + + mem = 8.d0 * dble(nOrb + nOrb**4 + N*n_states) / 1d6 + write(*,'(A40, F12.4)') 'I/O mem (MB) = ', mem + + mem = 8.d0 * dble(N + N*M + N*M + M*M + M*M + M + n_states_diag + n_states_diag) / 1d6 + write(*,'(A40, F12.4)') 'tmp mem (MB) = ', mem + + + call ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag) + + !print*, "H_diag:" + !do ab = 1, N + ! print*, ab, H_diag(ab) + !enddo + + ! initialize guess + R = 0.d0 + do k = 1, n_states + R(k,k) = 1.d0 + enddo + do k = n_states+1, n_states_diag + do i = 1, N + call random_number(r1) + call random_number(r2) + r1 = dsqrt(-2.d0*dlog(r1)) + r2 = dtwo_pi*r2 + R(i,k) = r1*dcos(r2) + enddo + R(k,k) = R(k,k) + 10.d0 + enddo + + do k = 1, n_states_diag + call normalize(R(1,k), N) + enddo + + !print*, 'guess' + !do k = 1, N + ! write(*,'(100(F15.7,2X))') (R(k,i), i = 1, n_states_diag) + !enddo + + ! working vectors + do k = 1, n_states_diag + do i = 1, N + U(i,k) = R(i,k) + enddo + enddo + + !print*, 'working vectors' + !do k = 1, N + ! write(*,'(100(F15.7,2X))') (U(k,i), i = 1, n_states_diag) + !enddo + + write(6,'(A)') '' + write_buffer = '=====' + do i = 1, n_states + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*n_states) + write_buffer = 'Iter' + do i = 1, n_states + write_buffer = trim(write_buffer)//' Energy Residual ' + enddo + write(6,'(A)') write_buffer(1:6+41*n_states) + write_buffer = '=====' + do i = 1, n_states + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') write_buffer(1:6+41*n_states) + + + converged = .False. + itertot = 0 + + do while (.not.converged) + + itertot = itertot + 1 + if(itertot == itermax) then + print*, 'exit before convergence !' + print*, 'itertot == itermax', itertot + exit + endif + + do iter = 1, itermax-1 + + shift1 = n_states_diag * (iter - 1) + shift2 = shift1 + n_states_diag + !print*, iter, shift1, shift2 + + if((iter > 1) .or. (itertot == 1)) then + + call ortho_qr(U(1,1), size(U, 1), N, shift2) + !call ortho_qr(U(1,1), size(U, 1), N, shift2) + + !print*, 'working vectors after qr' + !do k = 1, N + ! write(*,'(100(F15.7,2X))') (U(k,i), i = 1, n_states_diag) + !enddo + !allocate(S_check(shift2,shift2)) + !call dgemm("T", "N", shift2, shift2, N, 1.d0, U(1,1), size(U, 1), U(1,1), size(U, 1), 0.d0, S_check(1,1), size(S_check, 1)) + !do k = 1, shift2 + ! write(*,'(100(F15.7,2X))') (S_check(k,i), i = 1, shift2) + !enddo + !deallocate(S_check) + + call ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI(1,1,1,1), U(1,shift1+1), W(1,shift1+1)) + + else + + ! computed below + continue + endif + + ! h = U.T H U + call dgemm('T', 'N', shift2, shift2, N, 1.d0, & + U(1,1), size(U, 1), W(1,1), size(W, 1), & + 0.d0, h(1,1), size(h, 1)) + + ! h h_vec = h_val h_vec + call diag_nonsym_right(shift2, h(1,1), size(h, 1), h_vec(1,1), size(h_vec, 1), h_val(1), size(h_val, 1)) + !print*, 'h_val', h_val(1:shift2) + + ! U1 = U0 h_vec + call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0, & + U(1,1), size(U, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, U(1,shift2+1), size(U, 1)) + + do k = 1, n_states_diag + call normalize(U(1,shift2+k), N) + enddo + + do l = 1, n_states + !do k = 1, n_states_diag + ! overlap(k) = 0.d0 + ! do i = 1, N + ! overlap(k) = overlap(k) + U(i,shift2+k) * R(i,l) + ! enddo + ! overlap(k) = dabs(overlap(k)) + ! !print *, ' overlap =', k, overlap(k) + !enddo + !lambda_tmp = 0.d0 + !do k = 1, n_states_diag + ! if(overlap(k) .gt. lambda_tmp) then + ! i_omax(l) = k + ! lambda_tmp = overlap(k) + ! endif + !enddo + !if(lambda_tmp .lt. 0.7d0) then + ! print *, ' small overlap ...', l, i_omax(l) + ! print *, ' max overlap =', lambda_tmp + ! !stop + !endif + !if(i_omax(l) .ne. l) then + ! print *, ' !!! WARNING !!!' + ! print *, ' index of state', l, i_omax(l) + !endif + enddo + + ! W1 = W0 h_vec + call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0, & + W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, W(1,shift2+1), size(W, 1)) + + ! check if W1 = U1 h_val + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i, k) & + !$OMP SHARED(n_states, n_states_diag, N, shift2, U, h_val, W, H_diag, residual_norm, to_print) + !$OMP DO + do k = 1, n_states_diag + do i = 1, N + U(i,shift2+k) = (h_val(k) * U(i,shift2+k) - W(i,shift2+k)) / max(H_diag(i) - h_val(k), 1.d-2) + enddo + if(k <= n_states) then + residual_norm(k) = u_dot_u(U(1,shift2+k), N) + to_print(1,k) = h_val(k) + to_print(2,k) = residual_norm(k) + endif + enddo + !$OMP END DO + !$OMP END PARALLEL + !print*, " to_print", to_print + + if((itertot > 1) .and. (iter == 1)) then + continue + else + write(*,'(1X, I3, 1X, 100(1X, F16.10, 1X, F12.6))') iter-1, to_print(1:2,1:n_states) + !write(*, '(1X, I3, 1X, 100(1X, F16.10, 1X, F16.10, 1X, F16.10))') iter-1, to_print(1:2,1:n_states) + endif + + !print*, 'iter = ', iter + if(iter > 1) then + converged = dabs(maxval(residual_norm(1:n_states))) < 1d-15 + endif + + do k = 1, n_states + if(residual_norm(k) > 1.d8) then + print *, 'Davidson failed' + stop -1 + endif + enddo + + if(converged) exit + + enddo ! loop over iter + + + ! Re-contract U and update W + ! -------------------------------- + + call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0, & + W(1,1), size(W, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, R, size(R, 1)) + do k = 1, n_states_diag + do i = 1, N + W(i,k) = R(i,k) + enddo + enddo + + call dgemm('N', 'N', N, n_states_diag, shift2, 1.d0, & + U(1,1), size(U, 1), h_vec(1,1), size(h_vec, 1), & + 0.d0, R(1,1), size(R, 1)) + + do k = 1, n_states_diag + do i = 1, N + U(i,k) = R(i,k) + enddo + enddo + + call ortho_qr(U(1,1), size(U, 1), N, n_states_diag) + !call ortho_qr(U(1,1), size(U, 1), N, n_states_diag) + + do j = 1, n_states_diag + k = 1 + do while((k < N) .and. (U(k,j) == 0.d0)) + k = k+1 + enddo + if(U(k,j) * R(k,j) < 0.d0) then + do i = 1, N + W(i,j) = -W(i,j) + enddo + endif + enddo + + enddo ! loop over while + + ! --- + + write_buffer = '=====' + do i = 1, n_states + write_buffer = trim(write_buffer)//' ================ ===========' + enddo + write(6,'(A)') trim(write_buffer) + write(6,'(A)') '' + + + print*, " Davidson eigenvalues" + do k = 1, n_states + Om(k) = h_val(k) + print*, k, Om(k) + enddo + + deallocate(H_diag) + deallocate(U) + deallocate(W) + deallocate(h) + deallocate(h_vec) + deallocate(h_val) + deallocate(overlap) + deallocate(residual_norm) + + return +end + +! --- + +subroutine ppLR_HR_calc(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, n_states_diag, ERI, U, W) + + implicit none + + integer, intent(in) :: ispin + integer, intent(in) :: n_states_diag + integer, intent(in) :: nOO, nVV, nOrb, nC, nO, nR + double precision, intent(in) :: lambda, eF + double precision, intent(in) :: e(nOrb) + double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision, intent(in) :: U(nOO+nVV,n_states_diag) + double precision, intent(out) :: W(nOO+nVV,n_states_diag) + + integer :: i, j, ij, k, l, kl + integer :: a, b, c, d, ab, cd + integer :: state + double precision :: mat_tmp + double precision :: diff_loc, diff_tot + double precision, allocatable :: M_ref(:,:), W_ref(:,:) + double precision, allocatable :: Cpp_ref(:,:), Dpp_ref(:,:), Bpp_ref(:,:) + + double precision, external :: Kronecker_delta + + if(ispin .eq. 1) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a, nOrb-nR + ab = ab + 1 + + do state = 1, n_states_diag + + W(ab,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c, nOrb-nR + cd = cd + 1 + + mat_tmp = (e(a) + e(b) - eF) * Kronecker_delta(a, c) * Kronecker_delta(b, d) & + + lambda * (ERI(a,b,c,d) + ERI(a,b,d,c)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(c, d))) + + W(ab,state) = W(ab,state) + mat_tmp * U(cd,state) + enddo + enddo + + ij = nVV + do i = nC+1, nO + do j = i, nO + ij = ij + 1 + + mat_tmp = lambda * (ERI(a,b,i,j) + ERI(a,b,j,i)) / dsqrt( (1.d0 + Kronecker_delta(a, b)) & + * (1.d0 + Kronecker_delta(i, j))) + + W(ab,state) = W(ab,state) - mat_tmp * U(ij,state) + enddo + enddo + + enddo ! state + enddo ! b + enddo ! a + + ! --- + + ij = nVV + do i = nC+1, nO + do j = i, nO + ij = ij + 1 + + do state = 1, n_states_diag + + W(ij,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c, nOrb-nR + cd = cd + 1 + + mat_tmp = lambda * (ERI(c,d,i,j) + ERI(c,d,j,i)) / dsqrt( (1.d0 + Kronecker_delta(c, d)) & + * (1.d0 + Kronecker_delta(i, j))) + + W(ij,state) = W(ij,state) + mat_tmp * U(cd,state) + enddo + enddo + + kl = nVV + do k = nC+1, nO + do l = k, nO + kl = kl + 1 + + mat_tmp = - (e(i) + e(j) - eF) * Kronecker_delta(i, k) * Kronecker_delta(j, l) & + + lambda * (ERI(i,j,k,l) + ERI(i,j,l,k)) / dsqrt( (1.d0 + Kronecker_delta(i, j)) & + * (1.d0 + Kronecker_delta(k, l))) + + W(ij,state) = W(ij,state) - mat_tmp * U(kl,state) + enddo + enddo + + enddo ! state + enddo ! j + enddo ! i + + elseif((ispin .eq. 2) .or. (ispin .eq. 4)) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 + + do state = 1, n_states_diag + + W(ab,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + + mat_tmp = (e(a) + e(b) - eF) * Kronecker_delta(a, c) * Kronecker_delta(b, d) & + + lambda * (ERI(a,b,c,d) - ERI(a,b,d,c)) + + W(ab,state) = W(ab,state) + mat_tmp * U(cd,state) + enddo + enddo + + ij = nVV + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 + + mat_tmp = lambda * (ERI(a,b,i,j) - ERI(a,b,j,i)) + + W(ab,state) = W(ab,state) - mat_tmp * U(ij,state) + enddo + enddo + + enddo ! state + enddo ! b + enddo ! a + + ! --- + + ij = nVV + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 + + do state = 1, n_states_diag + + W(ij,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + + mat_tmp = lambda * (ERI(c,d,i,j) - ERI(c,d,j,i)) + + W(ij,state) = W(ij,state) + mat_tmp * U(cd,state) + enddo + enddo + + kl = nVV + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 + + mat_tmp = - (e(i) + e(j) - eF) * Kronecker_delta(i, k) * Kronecker_delta(j, l) & + + lambda * (ERI(i,j,k,l) - ERI(i,j,l,k)) + + W(ij,state) = W(ij,state) - mat_tmp * U(kl,state) + enddo + enddo + + enddo ! state + enddo ! j + enddo ! i + + elseif(ispin .eq. 3) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = nO+1, nOrb-nR + ab = ab + 1 + + do state = 1, n_states_diag + + W(ab,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = nO+1, nOrb-nR + cd = cd + 1 + + mat_tmp = (e(a) + e(b) - eF) * Kronecker_delta(a, c) * Kronecker_delta(b, d) & + + lambda * ERI(a,b,c,d) + + W(ab,state) = W(ab,state) + mat_tmp * U(cd,state) + enddo + enddo + + ij = nVV + do i = nC+1, nO + do j = nC+1, nO + ij = ij + 1 + + mat_tmp = lambda * ERI(a,b,i,j) + + W(ab,state) = W(ab,state) - mat_tmp * U(ij,state) + enddo + enddo + + enddo ! state + enddo ! b + enddo ! a + + ! --- + + ij = nVV + do i = nC+1, nO + do j = nC+1, nO + ij = ij + 1 + + do state = 1, n_states_diag + + W(ij,state) = 0.d0 + + cd = 0 + do c = nO+1, nOrb-nR + do d = nO+1, nOrb-nR + cd = cd + 1 + + mat_tmp = lambda * ERI(c,d,i,j) + + W(ij,state) = W(ij,state) + mat_tmp * U(cd,state) + enddo + enddo + + kl = nVV + do k = nC+1, nO + do l = nC+1, nO + kl = kl + 1 + + mat_tmp = - (e(i) + e(j) - eF) * Kronecker_delta(i, k) * Kronecker_delta(j, l) & + + lambda * ERI(i,j,k,l) + + W(ij,state) = W(ij,state) - mat_tmp * U(kl,state) + enddo + enddo + + enddo ! state + enddo ! j + enddo ! i + + else + + print*, ' ispin is not supported' + print*, ' ispin = ', ispin + stop + + endif + + +! print*, ' debug ppLR_HR_calc:' +! print*, ispin, nOO, nVV +! allocate(M_ref(nOO+nVV,nOO+nVV)) +! allocate(Bpp_ref(nVV,nOO), Cpp_ref(nVV,nVV), Dpp_ref(nOO,nOO)) +! allocate(W_ref(nOO+nVV,n_states_diag)) +! +! call ppLR_C(ispin, nOrb, nC, nO, nOrb-nO, nR, nVV, 1d0, e, ERI, Cpp_ref) +! call ppLR_D(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, 1d0, e, ERI, Dpp_ref) +! call ppLR_B(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, nVV, 1d0, ERI, Bpp_ref) +! M_ref = 0.d0 +! M_ref( 1:nVV , 1:nVV) = + Cpp_ref(1:nVV,1:nVV) +! M_ref(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp_ref(1:nOO,1:nOO) +! M_ref( 1:nVV ,nVV+1:nOO+nVV) = - Bpp_ref(1:nVV,1:nOO) +! M_ref(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp_ref(1:nVV,1:nOO)) +! +! call dgemm('N', 'N', nOO+nVV, n_states_diag, nOO+nVV, 1.d0, & +! M_ref(1,1), size(M_ref, 1), U(1,1), size(U, 1), & +! 0.d0, W_ref(1,1), size(W_ref, 1)) +! +! diff_tot = 0.d0 +! do state = 1, n_states_diag +! do ab = 1, nOO +! diff_loc = dabs(W(ab,state) - W_ref(ab,state)) +! if(diff_loc .gt. 1d-12) then +! print*, ' important diff on:', ab, state +! print*, W(ab,state), W_ref(ab,state) +! stop +! endif +! diff_tot = diff_tot + diff_loc +! enddo +! do ij = nVV+1, nVV+nOO +! diff_loc = dabs(W(ij,state) - W_ref(ij,state)) +! if(diff_loc .gt. 1d-12) then +! print*, ' important diff on:', ij, state +! print*, W(ij,state), W_ref(ij,state) +! stop +! endif +! diff_tot = diff_tot + diff_loc +! enddo +! enddo +! print*, 'diff_tot = ', diff_tot +! +! deallocate(M_ref) +! deallocate(Bpp_ref, Cpp_ref, Dpp_ref) +! deallocate(W_ref) + + return +end + +! --- + +subroutine ppLR_H_diag(ispin, nOrb, nC, nO, nR, nOO, nVV, lambda, e, eF, ERI, H_diag) + + implicit none + + integer, intent(in) :: ispin + integer, intent(in) :: nOO, nVV, nOrb, nC, nO, nR + double precision, intent(in) :: lambda, eF + double precision, intent(in) :: e(nOrb) + double precision, intent(in) :: ERI(nOrb,nOrb,nOrb,nOrb) + double precision, intent(out) :: H_diag(nOO+nVV) + + integer :: i, j, ij, k, l, kl + integer :: a, b, c, d, ab, cd + double precision :: diff_loc, diff_tot + double precision, allocatable :: M_ref(:,:) + double precision, allocatable :: Cpp_ref(:,:), Dpp_ref(:,:), Bpp_ref(:,:) + + double precision, external :: Kronecker_delta + + + if(ispin .eq. 1) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a, nOrb-nR + ab = ab + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = c, nOrb-nR + cd = cd + 1 + if(a .ne. c) cycle + if(b .ne. d) cycle + H_diag(ab) = e(a) + e(b) - eF & + + lambda * (ERI(a,b,c,d) + ERI(a,b,d,c)) / dsqrt((1.d0 + Kronecker_delta(a, b)) * (1.d0 + Kronecker_delta(c, d))) + enddo + enddo + enddo ! b + enddo ! a + + ij = nVV + do i = nC+1, nO + do j = i, nO + ij = ij + 1 + kl = 0 + do k = nC+1, nO + do l = k, nO + kl = kl + 1 + if(i .ne. k) cycle + if(j .ne. l) cycle + H_diag(ij) = (e(i) + e(j) - eF) & + - lambda * (ERI(i,j,k,l) + ERI(i,j,l,k)) / dsqrt((1.d0 + Kronecker_delta(i, j)) * (1.d0 + Kronecker_delta(k, l))) + enddo + enddo + enddo ! j + enddo ! i + + elseif((ispin .eq. 2) .or. (ispin .eq. 4)) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = a+1, nOrb-nR + ab = ab + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = c+1, nOrb-nR + cd = cd + 1 + if(a .ne. c) cycle + if(b .ne. d) cycle + H_diag(ab) = e(a) + e(b) - eF + lambda * (ERI(a,b,c,d) - ERI(a,b,d,c)) + enddo + enddo + enddo ! b + enddo ! a + + ij = nVV + do i = nC+1, nO + do j = i+1, nO + ij = ij + 1 + kl = 0 + do k = nC+1, nO + do l = k+1, nO + kl = kl + 1 + if(i .ne. k) cycle + if(j .ne. l) cycle + H_diag(ij) = e(i) + e(j) - eF - lambda * (ERI(i,j,k,l) - ERI(i,j,l,k)) + enddo + enddo + enddo ! j + enddo ! i + + elseif(ispin .eq. 3) then + + ab = 0 + do a = nO+1, nOrb-nR + do b = nO+1, nOrb-nR + ab = ab + 1 + cd = 0 + do c = nO+1, nOrb-nR + do d = nO+1, nOrb-nR + cd = cd + 1 + if(a .ne. c) cycle + if(b .ne. d) cycle + H_diag(ab) = (e(a) + e(b) - eF) + lambda * ERI(a,b,c,d) + enddo + enddo + enddo ! b + enddo ! a + + ij = nVV + do i = nC+1, nO + do j = nC+1, nO + ij = ij + 1 + kl = 0 + do k = nC+1, nO + do l = nC+1, nO + kl = kl + 1 + if(i .ne. k) cycle + if(j .ne. l) cycle + H_diag(ij) = (e(i) + e(j) - eF) - lambda * ERI(i,j,k,l) + enddo + enddo + enddo ! j + enddo ! i + + else + + print*, ' ispin is not supported' + print*, ' ispin = ', ispin + stop + + endif + + +! print*, ' debug ppLR_H_diag:' +! print*, ispin, nOO, nVV +! allocate(M_ref(nOO+nVV,nOO+nVV)) +! allocate(Bpp_ref(nVV,nOO), Cpp_ref(nVV,nVV), Dpp_ref(nOO,nOO)) +! +! call ppLR_C(ispin, nOrb, nC, nO, nOrb-nO, nR, nVV, 1d0, e, ERI, Cpp_ref) +! call ppLR_D(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, 1d0, e, ERI, Dpp_ref) +! call ppLR_B(ispin, nOrb, nC, nO, nOrb-nO, nR, nOO, nVV, 1d0, ERI, Bpp_ref) +! M_ref = 0.d0 +! M_ref( 1:nVV , 1:nVV) = + Cpp_ref(1:nVV,1:nVV) +! M_ref(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp_ref(1:nOO,1:nOO) +! M_ref( 1:nVV ,nVV+1:nOO+nVV) = - Bpp_ref(1:nVV,1:nOO) +! M_ref(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp_ref(1:nVV,1:nOO)) +! +! diff_tot = 0.d0 +! do ab = 1, nOO +! diff_loc = dabs(H_diag(ab) - M_ref(ab,ab)) +! if(diff_loc .gt. 1d-12) then +! print*, ' important diff on:', ab +! print*, H_diag(ab), M_ref(ab,ab) +! stop +! endif +! diff_tot = diff_tot + diff_loc +! enddo +! do ij = nVV+1, nVV+nOO +! diff_loc = dabs(H_diag(ij) - M_ref(ij,ij)) +! if(diff_loc .gt. 1d-12) then +! print*, ' important diff on:', ij +! print*, H_diag(ij), M_ref(ij,ij) +! stop +! endif +! diff_tot = diff_tot + diff_loc +! enddo +! print*, 'diff_tot = ', diff_tot +! +! deallocate(M_ref) +! deallocate(Bpp_ref, Cpp_ref, Dpp_ref) + + return +end + +! --- + +subroutine diag_nonsym_right(n, A, A_ldim, V, V_ldim, energy, E_ldim) + + implicit none + + integer, intent(in) :: n, A_ldim, V_ldim, E_ldim + double precision, intent(in) :: A(A_ldim,n) + double precision, intent(out) :: energy(E_ldim), V(V_ldim,n) + + character*1 :: JOBVL, JOBVR, BALANC, SENSE + integer :: i, j + integer :: ILO, IHI, lda, ldvl, ldvr, LWORK, INFO + double precision :: ABNRM + integer, allocatable :: iorder(:), IWORK(:) + double precision, allocatable :: WORK(:), SCALE_array(:), RCONDE(:), RCONDV(:) + double precision, allocatable :: Atmp(:,:), WR(:), WI(:), VL(:,:), VR(:,:), Vtmp(:) + double precision, allocatable :: energy_loc(:), V_loc(:,:) + + !print*, " n = ", n + allocate(Atmp(n,n)) + allocate(WI(n)) + allocate(WR(n)) + allocate(VR(n,n)) + allocate(VL(1,1)) + + do i = 1, n + do j = 1, n + Atmp(j,i) = A(j,i) + enddo + enddo + + JOBVL = "N" ! computes the left eigenvectors + JOBVR = "V" ! computes the right eigenvectors + BALANC = "B" ! Diagonal scaling and Permutation for optimization + SENSE = "V" ! Determines which reciprocal condition numbers are computed + lda = n + ldvr = n + ldvl = 1 + + allocate(WORK(1), SCALE_array(n), RCONDE(n), RCONDV(n), IWORK(2*n-2)) + + LWORK = -1 ! to ask for the optimal size of WORK + call dgeevx(BALANC, JOBVL, JOBVR, SENSE, & ! CHARACTERS + n, Atmp, lda, & ! MATRIX TO DIAGONALIZE + WR, WI, & ! REAL AND IMAGINARY PART OF EIGENVALUES + VL, ldvl, VR, ldvr, & ! LEFT AND RIGHT EIGENVECTORS + ILO, IHI, SCALE_array, ABNRM, RCONDE, RCONDV, & ! OUTPUTS OF OPTIMIZATION + WORK, LWORK, IWORK, INFO) + + if(INFO .ne. 0) then + print*, 'dgeevx failed !!', INFO + stop + endif + + LWORK = max(int(work(1)), 1) ! this is the optimal size of WORK + deallocate(WORK) + allocate(WORK(LWORK)) + call dgeevx(BALANC, JOBVL, JOBVR, SENSE, & + n, Atmp, lda, & + WR, WI, & + VL, ldvl, VR, ldvr, & + ILO, IHI, SCALE_array, ABNRM, RCONDE, RCONDV, & + WORK, LWORK, IWORK, INFO) + + if(INFO .ne. 0) then + print*, 'dgeevx failed !!', INFO + stop + endif + + deallocate(WORK, SCALE_array, RCONDE, RCONDV, IWORK) + deallocate(VL, Atmp) + + + allocate(energy_loc(n), V_loc(n,n)) + energy_loc = 0.d0 + V_loc = 0.d0 + + i = 1 + do while(i .le. n) + + !print*, i, WR(i), WI(i) + + if(dabs(WI(i)) .gt. 1d-7) then + + print*, ' Found an imaginary component to eigenvalue' + print*, ' Re(i) + Im(i)', i, WR(i), WI(i) + + energy_loc(i) = WR(i) + do j = 1, n + V_loc(j,i) = WR(i) * VR(j,i) - WI(i) * VR(j,i+1) + enddo + energy_loc(i+1) = WI(i) + do j = 1, n + V_loc(j,i+1) = WR(i) * VR(j,i+1) + WI(i) * VR(j,i) + enddo + i = i + 2 + + else + + energy_loc(i) = WR(i) + do j = 1, n + V_loc(j,i) = VR(j,i) + enddo + i = i + 1 + + endif + + enddo + + deallocate(WR, WI, VR) + + + ! ordering +! do j = 1, n +! write(444, '(100(1X, F16.10))') (V_loc(j,i), i=1,5) +! enddo + allocate(iorder(n)) + do i = 1, n + iorder(i) = i + enddo + call quick_sort(energy_loc, iorder, n) + do i = 1, n + energy(i) = energy_loc(i) + do j = 1, n + V(j,i) = V_loc(j,iorder(i)) + enddo + enddo + deallocate(iorder) +! do j = 1, n +! write(445, '(100(1X, F16.10))') (V_loc(j,i), i=1,5) +! enddo + deallocate(V_loc, energy_loc) + +end + +! --- + diff --git a/src/utils/non_sym_diag.f90 b/src/utils/non_sym_diag.f90 index 1659533..0bb19a6 100644 --- a/src/utils/non_sym_diag.f90 +++ b/src/utils/non_sym_diag.f90 @@ -519,7 +519,7 @@ subroutine impose_biorthog_svd(n, m, L, R) threshold = 1.d-6 num_linear_dependencies = 0 do i = 1, m - if(abs(D(i)) <= threshold) then + if(dabs(D(i)) <= threshold) then D(i) = 0.d0 num_linear_dependencies = num_linear_dependencies + 1 else diff --git a/src/utils/utils.f90 b/src/utils/utils.f90 index 316a87c..dfe511d 100644 --- a/src/utils/utils.f90 +++ b/src/utils/utils.f90 @@ -642,3 +642,113 @@ subroutine wall_time(t) t = dble(c)/dble(rate) end subroutine +! --- +! Compute +double precision function u_dot_u(u, sze) + + implicit none + + integer, intent(in) :: sze + double precision, intent(in) :: u(sze) + + double precision, external :: ddot + + !DIR$ FORCEINLINE + u_dot_u = ddot(sze,u,1,u,1) + +end + +! --- + +! Orthogonalization using Q.R factorization +subroutine ortho_qr(A, LDA, m, n) + + ! A : matrix to orthogonalize + ! LDA : leftmost dimension of A + ! m : Number of rows of A + ! n : Number of columns of A + ! /!\ int(WORK(1)) becomes negative when WORK(1) > 2147483648 + + integer, intent(in) :: m, n, LDA + double precision, intent(inout) :: A(LDA,n) + + integer :: LWORK, INFO + double precision, allocatable :: TAU(:), WORK(:) + + allocate(TAU(min(m, n))) + + allocate(WORK(1)) + + LWORK = -1 + + call dgeqrf(m, n, A(1,1), LDA, TAU, WORK, LWORK, INFO) + + if(INFO .lt. 0) then + print*, ' dgeqrf failed !! ', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + + deallocate(WORK) + allocate(WORK(LWORK)) + + call dgeqrf(m, n, A(1,1), LDA, TAU, WORK, LWORK, INFO) + + if(INFO .lt. 0) then + print*, ' dgeqrf failed !! ', INFO + stop + endif + + LWORK = -1 + call dorgqr(m, n, n, A(1,1), LDA, TAU, WORK, LWORK, INFO) + + if(INFO .lt. 0) then + print*, ' dorgqr failed !! ', INFO + stop + endif + + LWORK = max(n, int(WORK(1))) + + deallocate(WORK) + allocate(WORK(LWORK)) + + call dorgqr(m, n, n, A(1,1), LDA, TAU, WORK, LWORK, INFO) + + if(INFO .lt. 0) then + print*, ' dorgqr failed !! ', INFO + stop + endif + + deallocate(WORK, TAU) + +end + +! --- + +! Normalizes vector u +subroutine normalize(u, sze) + + implicit none + + integer, intent(in) :: sze + double precision, intent(inout) :: u(sze) + + integer :: i + double precision :: d + double precision, external :: dnrm2 + + !DIR$ FORCEINLINE + d = dnrm2(sze, u, 1) + if (d /= 0.d0) then + d = 1.d0/d + endif + if (d /= 1.d0) then + !DIR$ FORCEINLINE + call dscal(sze, d, u, 1) + endif + +end + +! --- +