From b96ea676fdebd3ba6357e458b9406cf20bb77b79 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 13 Apr 2020 11:33:48 +0200 Subject: [PATCH] b4 removing soG0T0 --- examples/molecule.H2 | 2 +- input/basis | 78 +++++++------ input/dft | 6 +- input/molecule | 6 +- input/molecule.xyz | 4 +- input/options | 2 +- input/weight | 78 +++++++------ src/QuAcK/G0T0.f90 | 109 ++++++++---------- src/QuAcK/QuAcK.f90 | 2 +- src/QuAcK/excitation_density_Tmatrix.f90 | 62 +++++++++- src/QuAcK/linear_response_B_pp.f90 | 14 +-- src/QuAcK/linear_response_C_pp.f90 | 12 +- src/QuAcK/linear_response_D_pp.f90 | 12 +- src/QuAcK/print_excitation.f90 | 3 +- src/QuAcK/renormalization_factor_Tmatrix.f90 | 65 +++-------- src/QuAcK/soG0T0.f90 | 2 +- ..._lda_exchange_derivative_discontinuity.f90 | 12 +- src/eDFT/RGIC_lda_exchange_energy.f90 | 12 +- .../RGIC_lda_exchange_individual_energy.f90 | 13 ++- src/eDFT/RGIC_lda_exchange_potential.f90 | 13 ++- 20 files changed, 276 insertions(+), 231 deletions(-) diff --git a/examples/molecule.H2 b/examples/molecule.H2 index cd43fec..8076140 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0.0 0.0 0.0 - H 0.0 0.0 3.7 + H 0.0 0.0 1.4 diff --git a/input/basis b/input/basis index dc1936c..20bc6a7 100644 --- a/input/basis +++ b/input/basis @@ -1,43 +1,49 @@ 1 9 +S 8 + 1 1469.0000000 0.0007660 + 2 220.5000000 0.0058920 + 3 50.2600000 0.0296710 + 4 14.2400000 0.1091800 + 5 4.5810000 0.2827890 + 6 1.5800000 0.4531230 + 7 0.5640000 0.2747740 + 8 0.0734500 0.0097510 +S 8 + 1 1469.0000000 -0.0001200 + 2 220.5000000 -0.0009230 + 3 50.2600000 -0.0046890 + 4 14.2400000 -0.0176820 + 5 4.5810000 -0.0489020 + 6 1.5800000 -0.0960090 + 7 0.5640000 -0.1363800 + 8 0.0734500 0.5751020 +S 1 + 1 0.0280500 1.0000000 +S 1 + 1 0.0086400 1.0000000 +P 3 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 +P 1 + 1 0.0240300 1.0000000 +P 1 + 1 0.0057900 1.0000000 +D 1 + 1 0.1239000 1.0000000 +D 1 + 1 0.0725000 1.0000000 +2 5 S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.3258000 1.0000000 + 1 0.1220000 1.0000000 S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 + 1 0.0297400 1.0000000 P 1 - 1 1.4070000 1.0000000 + 1 0.7270000 1.0000000 P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 -D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 -2 9 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 -S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 -P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 -D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 + 1 0.1410000 1.0000000 diff --git a/input/dft b/input/dft index 883e724..1c8aaa6 100644 --- a/input/dft +++ b/input/dft @@ -1,12 +1,12 @@ # Restricted or unrestricted KS calculation - LIM-RKS + GOK-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 # GGA = 2: RB88 # Hybrid = 4 # Hartree-Fock = 666 - 666 HF + 1 RS51 # correlation rung: # Hartree = 0 # LDA = 1: RVWN5,RMFL20 @@ -19,6 +19,6 @@ # Number of states in ensemble (nEns) 2 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 1.0 + 0.0 # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.00001 T 5 1 1 diff --git a/input/molecule b/input/molecule index cd43fec..e76247c 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd - 2 1 1 0 0 + 2 2 2 0 0 # Znuc x y z - H 0.0 0.0 0.0 - H 0.0 0.0 3.7 + Li 0. 0. 0. + H 0. 0. 3.099 diff --git a/input/molecule.xyz b/input/molecule.xyz index 45af661..fb94244 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,4 @@ 2 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 1.9579558213 + Li 0.0000000000 0.0000000000 0.0000000000 + H 0.0000000000 0.0000000000 1.6399202947 diff --git a/input/options b/input/options index f20ae20..b9af96e 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 F F F F F 0.000 + 256 0.00001 T 5 F F F F F F T 0.000 # ACFDT: AC Kx XBS T T T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/input/weight b/input/weight index dc1936c..20bc6a7 100644 --- a/input/weight +++ b/input/weight @@ -1,43 +1,49 @@ 1 9 +S 8 + 1 1469.0000000 0.0007660 + 2 220.5000000 0.0058920 + 3 50.2600000 0.0296710 + 4 14.2400000 0.1091800 + 5 4.5810000 0.2827890 + 6 1.5800000 0.4531230 + 7 0.5640000 0.2747740 + 8 0.0734500 0.0097510 +S 8 + 1 1469.0000000 -0.0001200 + 2 220.5000000 -0.0009230 + 3 50.2600000 -0.0046890 + 4 14.2400000 -0.0176820 + 5 4.5810000 -0.0489020 + 6 1.5800000 -0.0960090 + 7 0.5640000 -0.1363800 + 8 0.0734500 0.5751020 +S 1 + 1 0.0280500 1.0000000 +S 1 + 1 0.0086400 1.0000000 +P 3 + 1 1.5340000 0.0227840 + 2 0.2749000 0.1391070 + 3 0.0736200 0.5003750 +P 1 + 1 0.0240300 1.0000000 +P 1 + 1 0.0057900 1.0000000 +D 1 + 1 0.1239000 1.0000000 +D 1 + 1 0.0725000 1.0000000 +2 5 S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.3258000 1.0000000 + 1 0.1220000 1.0000000 S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 + 1 0.0297400 1.0000000 P 1 - 1 1.4070000 1.0000000 + 1 0.7270000 1.0000000 P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 -D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 -2 9 -S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 -S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 -S 1 - 1 0.0252600 1.0000000 -P 1 - 1 1.4070000 1.0000000 -P 1 - 1 0.3880000 1.0000000 -P 1 - 1 0.1020000 1.0000000 -D 1 - 1 1.0570000 1.0000000 -D 1 - 1 0.2470000 1.0000000 + 1 0.1410000 1.0000000 diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index 334908a..d3c471c 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -32,10 +32,12 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Local variables integer :: ispin + integer :: iblock integer :: nOOs,nOOt integer :: nVVs,nVVt double precision :: dERI double precision :: xERI + double precision :: alpha double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcAC(nspin) @@ -69,8 +71,11 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Dimensions of the pp-RPA linear reponse matrices - nOOs = nO*(nO + 1)/2 - nVVs = nV*(nV + 1)/2 + nOOs = nO*nO + nVVs = nV*nV + +! nOOs = nO*(nO + 1)/2 +! nVVs = nV*(nV + 1)/2 nOOt = nO*(nO - 1)/2 nVVt = nV*(nV - 1)/2 @@ -86,100 +91,76 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m SigT(nBas),Z(nBas),eG0T0(nBas)) !---------------------------------------------- -! Singlet manifold +! alpha-beta block !---------------------------------------------- - ispin = 1 + ispin = 1 + iblock = 3 ! Compute linear response - call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:),ERI(:,:,:,:), & - Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:), & - EcRPA(ispin)) + call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:),ERI(:,:,:,:), & + Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) EcRPA(ispin) = 1d0*EcRPA(ispin) - call print_excitation('pp-RPA (N+2)',ispin,nVVs,Omega1s(:)) - call print_excitation('pp-RPA (N-2)',ispin,nOOs,Omega2s(:)) + call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:)) + call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:)) !---------------------------------------------- -! Triplet manifold +! alpha-alpha block !---------------------------------------------- - ispin = 2 + ispin = 2 + iblock = 4 ! Compute linear response - call linear_response_pp(ispin,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:),ERI(:,:,:,:), & - Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:), & - EcRPA(ispin)) + call linear_response_pp(iblock,.true.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:),ERI(:,:,:,:), & + Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) - EcRPA(ispin) = 3d0*EcRPA(ispin) + EcRPA(ispin) = 2d0*EcRPA(ispin) +! EcRPA(ispin) = 3d0*EcRPA(ispin) - call print_excitation('pp-RPA (N+2)',ispin,nVVt,Omega1t(:)) - call print_excitation('pp-RPA (N-2)',ispin,nOOt,Omega2t(:)) + call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:)) + call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- SigT(:) = 0d0 + Z(:) = 0d0 - ispin = 2 - dERI = +1d0 - xERI = -1d0 + iblock = 3 + dERI = +1d0 + xERI = +0d0 + alpha = +1d0 - call excitation_density_Tmatrix(ispin,xERI,dERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & - X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) - - call self_energy_Tmatrix_diag(1d0,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:), & - Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),SigT(:)) - - ispin = 2 - dERI = +1d0 - xERI = +0d0 - - call excitation_density_Tmatrix(ispin,xERI,dERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & - X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) - - call self_energy_Tmatrix_diag(0.5d0,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:), & - Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),SigT(:)) - - ispin = 2 - dERI = +0d0 - xERI = +1d0 - - call excitation_density_Tmatrix(ispin,xERI,dERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & - X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) - - call self_energy_Tmatrix_diag(0.5d0,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:), & - Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),SigT(:)) - - ispin = 1 - dERI = +1d0 - xERI = +0d0 - call excitation_density_Tmatrix(ispin,xERI,dERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & + call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) - call self_energy_Tmatrix_diag(0.5d0,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:), & + call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:), & Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),SigT(:)) - ispin = 1 - dERI = +0d0 - xERI = +1d0 - call excitation_density_Tmatrix(ispin,xERI,dERI,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI(:,:,:,:), & - X1s(:,:),Y1s(:,:),rho1s(:,:,:),X2s(:,:),Y2s(:,:),rho2s(:,:,:)) + call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:), & + Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),Z(:)) - call self_energy_Tmatrix_diag(0.5d0,eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF(:), & - Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:),SigT(:)) + iblock = 4 + dERI = +1d0 + xERI = -1d0 + alpha = +1d0 -! Compute renormalization factor for T-matrix self-energy + call excitation_density_Tmatrix(iblock,dERI,xERI,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI(:,:,:,:), & + X1t(:,:),Y1t(:,:),rho1t(:,:,:),X2t(:,:),Y2t(:,:),rho2t(:,:,:)) -! call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,eHF(:), & -! Omega1s(:),rho1s(:,:,:),Omega2s(:),rho2s(:,:,:), & -! Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:), & -! Z(:)) - Z(:) = 1d0 + call self_energy_Tmatrix_diag(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:), & + Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),SigT(:)) + + call renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF(:), & + Omega1t(:),rho1t(:,:,:),Omega2t(:),rho2t(:,:,:),Z(:)) + + Z(:) = 1d0/(1d0 - Z(:)) !---------------------------------------------- ! Solve the quasi-particle equation diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index c2041c3..aa5b86b 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -685,7 +685,7 @@ program QuAcK if(doG0T0) then call cpu_time(start_G0T0) - call soG0T0(eta,nBas,nC(1),nO(1),nV(1),nR(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) diff --git a/src/QuAcK/excitation_density_Tmatrix.f90 b/src/QuAcK/excitation_density_Tmatrix.f90 index 1c225ec..ffc8c4e 100644 --- a/src/QuAcK/excitation_density_Tmatrix.f90 +++ b/src/QuAcK/excitation_density_Tmatrix.f90 @@ -104,7 +104,7 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E ! Triplet manifold !---------------------------------------------- - if(ispin == 2) then + if(ispin == 2 .or. ispin == 4) then do p=nC+1,nBas-nR @@ -160,4 +160,64 @@ subroutine excitation_density_Tmatrix(ispin,dERI,xERI,nBas,nC,nO,nV,nR,nOO,nVV,E end if +!---------------------------------------------- +! alpha-beta block +!---------------------------------------------- + + if(ispin == 3) then + + do p=nC+1,nBas-nR + + do i=nC+1,nO + do ab=1,nVV + + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + cd = cd + 1 + rho1(p,i,ab) = rho1(p,i,ab) & + + (dERI*ERI(p,i,c,d) + xERI*ERI(p,i,d,c))*X1(cd,ab) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=nC+1,nO + kl = kl + 1 + rho1(p,i,ab) = rho1(p,i,ab) & + + (dERI*ERI(p,i,k,l) + xERI*ERI(p,i,l,k))*Y1(kl,ab) + end do + end do + + end do + end do + + do a=1,nV-nR + do ij=1,nOO + + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + cd = cd + 1 + rho2(p,a,ij) = rho2(p,a,ij) & + + (dERI*ERI(p,nO+a,c,d) + xERI*ERI(p,nO+a,d,c))*X2(cd,ij) + end do + end do + + kl = 0 + do k=nC+1,nO + do l=nC+1,nO + kl = kl + 1 + rho2(p,a,ij) = rho2(p,a,ij) & + + (dERI*ERI(p,nO+a,k,l) + xERI*ERI(p,nO+a,l,k))*Y2(kl,ij) + end do + end do + + end do + end do + + end do + + end if + end subroutine excitation_density_Tmatrix diff --git a/src/QuAcK/linear_response_B_pp.f90 b/src/QuAcK/linear_response_B_pp.f90 index baa2079..fe2515d 100644 --- a/src/QuAcK/linear_response_B_pp.f90 +++ b/src/QuAcK/linear_response_B_pp.f90 @@ -42,20 +42,20 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) end if -! Build the B matrix for the triplet manifold +! Build the alpha-beta block of the B matrix - if(ispin == 2) then + if(ispin == 3) then ab = 0 do a=nO+1,nBas-nR - do b=a+1,nBas-nR + do b=nO+1,nBas-nR ab = ab + 1 ij = 0 do i=nC+1,nO - do j=i+1,nO + do j=nC+1,nO ij = ij + 1 - B_pp(ab,ij) = ERI(a,b,i,j) - ERI(a,b,j,i) + B_pp(ab,ij) = ERI(a,b,i,j) end do end do @@ -64,9 +64,9 @@ subroutine linear_response_B_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) end if -! Build the B matrix for the spinorbital basis +! Build the B matrix for the triplet manifold, or alpha-alpha, or in the spin-orbital basis - if(ispin == 3) then + if(ispin == 2 .or. ispin == 4) then ab = 0 do a=nO+1,nBas-nR diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90 index ae7d399..b1ba603 100644 --- a/src/QuAcK/linear_response_C_pp.f90 +++ b/src/QuAcK/linear_response_C_pp.f90 @@ -50,9 +50,9 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) end if -! Build C matrix for the triplet manifold +! Build C matrix for the triplet manifold, or alpha-alpha block, or in the spin-orbital basis - if(ispin == 2) then + if(ispin == 2 .or. ispin == 4) then ab = 0 do a=nO+1,nBas-nR @@ -73,21 +73,21 @@ subroutine linear_response_C_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) end if -! Build C matrix for the spinorbital basis +! Build the alpha-beta block of the C matrix if(ispin == 3) then ab = 0 do a=nO+1,nBas-nR - do b=a+1,nBas-nR + do b=nO+1,nBas-nR ab = ab + 1 cd = 0 do c=nO+1,nBas-nR - do d=c+1,nBas-nR + do d=nO+1,nBas-nR cd = cd + 1 C_pp(ab,cd) = + (e(a) + e(b) - eF)*Kronecker_delta(a,c)*Kronecker_delta(b,d) & - + ERI(a,b,c,d) - ERI(a,b,d,c) + + ERI(a,b,c,d) end do end do diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90 index 002a6bc..0a85d64 100644 --- a/src/QuAcK/linear_response_D_pp.f90 +++ b/src/QuAcK/linear_response_D_pp.f90 @@ -50,9 +50,9 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) end if -! Build the D matrix for the triplet manifold +! Build the D matrix for the triplet manifold, the alpha-alpha block, or in the spin-orbital basis - if(ispin == 2) then + if(ispin == 2 .or. ispin == 4) then ij = 0 do i=nC+1,nO @@ -73,21 +73,21 @@ subroutine linear_response_D_pp(ispin,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) end if -! Build the D matrix for the spinorbital basis +! Build the alpha-beta block of the D matrix if(ispin == 3) then ij = 0 do i=nC+1,nO - do j=i+1,nO + do j=nC+1,nO ij = ij + 1 kl = 0 do k=nC+1,nO - do l=k+1,nO + do l=nC+1,nO kl = kl + 1 D_pp(ij,kl) = - (e(i) + e(j) - eF)*Kronecker_delta(i,k)*Kronecker_delta(j,l) & - + ERI(i,j,k,l) - ERI(i,j,l,k) + + ERI(i,j,k,l) end do end do diff --git a/src/QuAcK/print_excitation.f90 b/src/QuAcK/print_excitation.f90 index fca9912..ea8fa8b 100644 --- a/src/QuAcK/print_excitation.f90 +++ b/src/QuAcK/print_excitation.f90 @@ -19,7 +19,8 @@ subroutine print_excitation(method,ispin,nS,Omega) if(ispin == 1) spin_manifold = 'singlet' if(ispin == 2) spin_manifold = 'triplet' - if(ispin == 3) spin_manifold = 'spinorb' + if(ispin == 3) spin_manifold = 'alp-bet' + if(ispin == 4) spin_manifold = 'alp-alp' write(*,*) write(*,*)'-------------------------------------------------------------' diff --git a/src/QuAcK/renormalization_factor_Tmatrix.f90 b/src/QuAcK/renormalization_factor_Tmatrix.f90 index ce4be02..b09bce8 100644 --- a/src/QuAcK/renormalization_factor_Tmatrix.f90 +++ b/src/QuAcK/renormalization_factor_Tmatrix.f90 @@ -1,6 +1,4 @@ -subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nVVt,e, & - Omega1s,rho1s,Omega2s,rho2s,Omega1t,rho1t,Omega2t,rho2t, & - Z) +subroutine renormalization_factor_Tmatrix(alpha,eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,Z) ! Compute renormalization factor of the T-matrix self-energy @@ -9,15 +7,16 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV ! Input variables + double precision,intent(in) :: alpha double precision,intent(in) :: eta integer,intent(in) :: nBas,nC,nO,nV,nR - integer,intent(in) :: nOOs,nOOt - integer,intent(in) :: nVVs,nVVt + integer,intent(in) :: nOO + integer,intent(in) :: nVV double precision,intent(in) :: e(nBas) - double precision,intent(in) :: Omega1s(nVVs),Omega1t(nVVt) - double precision,intent(in) :: rho1s(nBas,nO,nVVs),rho1t(nBas,nO,nVVt) - double precision,intent(in) :: Omega2s(nOOs),Omega2t(nOOt) - double precision,intent(in) :: rho2s(nBas,nV,nOOs),rho2t(nBas,nV,nOOt) + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: rho1(nBas,nO,nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho2(nBas,nV,nOO) ! Local variables @@ -28,21 +27,13 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV double precision,intent(out) :: Z(nBas) -! Initialize - - Z(:) = 0d0 - -!---------------------------------------------- -! Singlet part of the T-matrix self-energy -!---------------------------------------------- - ! Occupied part of the T-matrix self-energy do p=nC+1,nBas-nR do i=nC+1,nO - do cd=1,nVVs - eps = e(p) + e(i) - Omega1s(cd) - Z(p) = Z(p) + (rho1s(p,i,cd)/eps)**2 + do cd=1,nVV + eps = e(p) + e(i) - Omega1(cd) + Z(p) = Z(p) - (rho1(p,i,cd)/eps)**2 enddo enddo enddo @@ -51,41 +42,15 @@ subroutine renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,nOOt,nV do p=nC+1,nBas-nR do a=1,nV-nR - do kl=1,nOOs - eps = e(p) + e(nO+a) - Omega2s(kl) - Z(p) = Z(p) + (rho2s(p,a,kl)/eps)**2 - enddo - enddo - enddo - -!---------------------------------------------- -! Triplet part of the T-matrix self-energy -!---------------------------------------------- - -! Occupied part of the T-matrix self-energy - - do p=nC+1,nBas-nR - do i=nC+1,nO - do cd=1,nVVt - eps = e(p) + e(i) - Omega1t(cd) - Z(p) = Z(p) + (rho1t(p,i,cd)/eps)**2 - enddo - enddo - enddo - -! Virtual part of the T-matrix self-energy - - do p=nC+1,nBas-nR - do a=1,nV-nR - do kl=1,nOOt - eps = e(p) + e(nO+a) - Omega2t(kl) - Z(p) = Z(p) + (rho2t(p,a,kl)/eps)**2 + do kl=1,nOO + eps = e(p) + e(nO+a) - Omega2(kl) + Z(p) = Z(p) - (rho2(p,a,kl)/eps)**2 enddo enddo enddo ! Compute renormalization factor from derivative of SigT - Z(:) = 1d0/(1d0 + Z(:)) +! Z(:) = 1d0/(1d0 + Z(:)) end subroutine renormalization_factor_Tmatrix diff --git a/src/QuAcK/soG0T0.f90 b/src/QuAcK/soG0T0.f90 index 77ebc35..90baa42 100644 --- a/src/QuAcK/soG0T0.f90 +++ b/src/QuAcK/soG0T0.f90 @@ -75,7 +75,7 @@ subroutine soG0T0(eta,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,eHF) ! Spinorbital basis !---------------------------------------------- - ispin = 3 + ispin = 4 ! Compute linear response diff --git a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 index 5f8e758..479879a 100644 --- a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 @@ -49,9 +49,15 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho ! Parameters for He - a = 1.9125735895875828d0 - b = 2.715266992840757d0 - c = 2.1634223380633086d0 +! a = 1.9125735895875828d0 +! b = 2.715266992840757d0 +! c = 2.1634223380633086d0 + +! Parameters for HNO + + a = 0.0061158387543040335d0 + b = -0.00005968703047293955d0 + c = -0.00001692245714408755d0 w = wEns(2) dCxGICdw = (0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0))) diff --git a/src/eDFT/RGIC_lda_exchange_energy.f90 b/src/eDFT/RGIC_lda_exchange_energy.f90 index fb60b5d..bee986b 100644 --- a/src/eDFT/RGIC_lda_exchange_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_energy.f90 @@ -41,9 +41,15 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Parameters for He - a = 1.9125735895875828d0 - b = 2.715266992840757d0 - c = 2.1634223380633086d0 +! a = 1.9125735895875828d0 +! b = 2.715266992840757d0 +! c = 2.1634223380633086d0 + +! Parameters for HNO + + a = 0.0061158387543040335d0 + b = -0.00005968703047293955d0 + c = -0.00001692245714408755d0 w = wEns(2) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) diff --git a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 index e2b7b42..ecb5849 100644 --- a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 @@ -44,9 +44,16 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E ! Parameters for He - a = 1.9125735895875828d0 - b = 2.715266992840757d0 - c = 2.1634223380633086d0 +! a = 1.9125735895875828d0 +! b = 2.715266992840757d0 +! c = 2.1634223380633086d0 + + +! Parameters for HNO + + a = 0.0061158387543040335d0 + b = -0.00005968703047293955d0 + c = -0.00001692245714408755d0 w = wEns(2) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) diff --git a/src/eDFT/RGIC_lda_exchange_potential.f90 b/src/eDFT/RGIC_lda_exchange_potential.f90 index ad6e66f..b200c05 100644 --- a/src/eDFT/RGIC_lda_exchange_potential.f90 +++ b/src/eDFT/RGIC_lda_exchange_potential.f90 @@ -44,9 +44,16 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Parameters for He - a = 1.9125735895875828d0 - b = 2.715266992840757d0 - c = 2.1634223380633086d0 +! a = 1.9125735895875828d0 +! b = 2.715266992840757d0 +! c = 2.1634223380633086d0 + + +! Parameters for HNO + + a = 0.0061158387543040335d0 + b = -0.00005968703047293955d0 + c = -0.00001692245714408755d0 w = wEns(2) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2)