diff --git a/input/basis b/input/basis index 120ad98..b9ca7b5 100644 --- a/input/basis +++ b/input/basis @@ -1,58 +1,9 @@ -1 6 -S 8 1.00 - 9046.0000000 0.0007000 - 1357.0000000 0.0053890 - 309.3000000 0.0274060 - 87.7300000 0.1032070 - 28.5600000 0.2787230 - 10.2100000 0.4485400 - 3.8380000 0.2782380 - 0.7466000 0.0154400 -S 8 1.00 - 9046.0000000 -0.0001530 - 1357.0000000 -0.0012080 - 309.3000000 -0.0059920 - 87.7300000 -0.0245440 - 28.5600000 -0.0674590 - 10.2100000 -0.1580780 - 3.8380000 -0.1218310 - 0.7466000 0.5490030 +1 3 +S 3 1.00 + 38.3600000 0.0238090 + 5.7700000 0.1548910 + 1.2400000 0.4699870 S 1 1.00 - 0.2248000 1.0000000 -P 3 1.00 - 13.5500000 0.0399190 - 2.9170000 0.2171690 - 0.7973000 0.5103190 + 0.2976000 1.0000000 P 1 1.00 - 0.2185000 1.0000000 -D 1 1.00 - 0.8170000 1.0000000 -2 6 -S 8 1.00 - 9046.0000000 0.0007000 - 1357.0000000 0.0053890 - 309.3000000 0.0274060 - 87.7300000 0.1032070 - 28.5600000 0.2787230 - 10.2100000 0.4485400 - 3.8380000 0.2782380 - 0.7466000 0.0154400 -S 8 1.00 - 9046.0000000 -0.0001530 - 1357.0000000 -0.0012080 - 309.3000000 -0.0059920 - 87.7300000 -0.0245440 - 28.5600000 -0.0674590 - 10.2100000 -0.1580780 - 3.8380000 -0.1218310 - 0.7466000 0.5490030 -S 1 1.00 - 0.2248000 1.0000000 -P 3 1.00 - 13.5500000 0.0399190 - 2.9170000 0.2171690 - 0.7973000 0.5103190 -P 1 1.00 - 0.2185000 1.0000000 -D 1 1.00 - 0.8170000 1.0000000 + 1.2750000 1.0000000 diff --git a/input/methods b/input/methods index 46c8861..2b692b1 100644 --- a/input/methods +++ b/input/methods @@ -5,10 +5,10 @@ # CCD CCSD CCSD(T) F F F # CIS TDHF ADC - T T F + F T F # GF2 GF3 F F # G0W0 evGW qsGW - T T F + F F F # MCMP2 F diff --git a/input/molecule b/input/molecule index 93b448f..c78e87e 100644 --- a/input/molecule +++ b/input/molecule @@ -1,5 +1,4 @@ # nAt nEla nElb nCore nRyd - 2 7 7 0 0 + 1 1 1 0 0 # Znuc x y z - N 0. 0. 0. - N 0. 0. 3.4 + He 0.0 0.0 0.0 diff --git a/input/options b/input/options index 23aacb6..1b1eb4a 100644 --- a/input/options +++ b/input/options @@ -4,8 +4,8 @@ # CC: maxSCF thresh DIIS n_diis 64 0.00001 F 1 -# CIS/TDHF: singlet triplet - T T +# CIS/TDHF/BSE: singlet triplet + T F # GF: maxSCF thresh DIIS n_diis renormalization 64 0.00001 T 10 3 # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 linearize eta diff --git a/input/weight b/input/weight index 120ad98..b9ca7b5 100644 --- a/input/weight +++ b/input/weight @@ -1,58 +1,9 @@ -1 6 -S 8 1.00 - 9046.0000000 0.0007000 - 1357.0000000 0.0053890 - 309.3000000 0.0274060 - 87.7300000 0.1032070 - 28.5600000 0.2787230 - 10.2100000 0.4485400 - 3.8380000 0.2782380 - 0.7466000 0.0154400 -S 8 1.00 - 9046.0000000 -0.0001530 - 1357.0000000 -0.0012080 - 309.3000000 -0.0059920 - 87.7300000 -0.0245440 - 28.5600000 -0.0674590 - 10.2100000 -0.1580780 - 3.8380000 -0.1218310 - 0.7466000 0.5490030 +1 3 +S 3 1.00 + 38.3600000 0.0238090 + 5.7700000 0.1548910 + 1.2400000 0.4699870 S 1 1.00 - 0.2248000 1.0000000 -P 3 1.00 - 13.5500000 0.0399190 - 2.9170000 0.2171690 - 0.7973000 0.5103190 + 0.2976000 1.0000000 P 1 1.00 - 0.2185000 1.0000000 -D 1 1.00 - 0.8170000 1.0000000 -2 6 -S 8 1.00 - 9046.0000000 0.0007000 - 1357.0000000 0.0053890 - 309.3000000 0.0274060 - 87.7300000 0.1032070 - 28.5600000 0.2787230 - 10.2100000 0.4485400 - 3.8380000 0.2782380 - 0.7466000 0.0154400 -S 8 1.00 - 9046.0000000 -0.0001530 - 1357.0000000 -0.0012080 - 309.3000000 -0.0059920 - 87.7300000 -0.0245440 - 28.5600000 -0.0674590 - 10.2100000 -0.1580780 - 3.8380000 -0.1218310 - 0.7466000 0.5490030 -S 1 1.00 - 0.2248000 1.0000000 -P 3 1.00 - 13.5500000 0.0399190 - 2.9170000 0.2171690 - 0.7973000 0.5103190 -P 1 1.00 - 0.2185000 1.0000000 -D 1 1.00 - 0.8170000 1.0000000 + 1.2750000 1.0000000 diff --git a/src/QuAcK/TDHF.f90 b/src/QuAcK/TDHF.f90 index de05458..adf4d3a 100644 --- a/src/QuAcK/TDHF.f90 +++ b/src/QuAcK/TDHF.f90 @@ -66,10 +66,16 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, ispin = 1 - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, & - rho,EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & + EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) + call linear_response_ph(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & + EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + call print_excitation('p-h TDHF ',ispin,nS,Omega(:,ispin)) + + call linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho) + endif ! Triplet manifold @@ -78,8 +84,8 @@ subroutine TDHF(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,nS,ENuc,ERHF, ispin = 2 - call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI, & - rho,EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) + call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho, & + EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin)) call print_excitation('TDHF ',ispin,nS,Omega(:,ispin)) endif diff --git a/src/QuAcK/huckel_guess.f90 b/src/QuAcK/huckel_guess.f90 new file mode 100644 index 0000000..07fea2a --- /dev/null +++ b/src/QuAcK/huckel_guess.f90 @@ -0,0 +1,61 @@ +subroutine huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P) + +! Hickel guess of the molecular orbitals for HF calculation + + implicit none + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(inout):: J(nBas,nBas) + double precision,intent(inout):: K(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(inout):: cp(nBas,nBas) + double precision,intent(inout):: Fp(nBas,nBas) + double precision,intent(inout):: e(nBas) + double precision,intent(inout):: P(nBas,nBas) + +! Local variables + + integer :: mu,nu + double precision :: a + +! Output variables + + double precision,intent(out) :: c(nBas,nBas) + + a = 1.75d0 + + Fp = matmul(transpose(X),matmul(Hc,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,e) + c = matmul(X,cp) + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + + call Coulomb_matrix_AO_basis(nBas,P,ERI,J) + call exchange_matrix_AO_basis(nBas,P,ERI,K) + + do mu=1,nBas + + Fp(mu,mu) = Hc(mu,mu) + J(mu,mu) + 0.5d0*K(mu,mu) + + do nu=mu+1,nBas + + Fp(mu,nu) = 0.5d0*a*S(mu,nu)*(Hc(mu,mu) + Hc(nu,nu)) + Fp(nu,mu) = Fp(mu,nu) + + enddo + + enddo + + Fp = matmul(transpose(X),matmul(Fp,X)) + + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,e) + c = matmul(X,cp) + +end subroutine diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90 index 3d86ddc..a4313e3 100644 --- a/src/QuAcK/linear_response.f90 +++ b/src/QuAcK/linear_response.f90 @@ -90,6 +90,6 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho,EcRP EcRPA = 0.5d0*(sum(Omega) - trace_matrix(nS,A)) -! print*,'EcRPA = ',EcRPA + print*,'EcRPA = ',EcRPA end subroutine linear_response diff --git a/src/QuAcK/linear_response_B_pp.f90 b/src/QuAcK/linear_response_B_pp.f90 new file mode 100644 index 0000000..8ab931a --- /dev/null +++ b/src/QuAcK/linear_response_B_pp.f90 @@ -0,0 +1,56 @@ +subroutine linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B_pp) + +! Compute the B matrix of the pp channel + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: dRPA + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: delta_spin + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: a,b,i,j,ab,ij + +! Output variables + + double precision,intent(out) :: B_pp(nVV,nOO) + +! Singlet or triplet manifold? + + delta_spin = 0d0 + if(ispin == 1) delta_spin = +1d0 + if(ispin == 2) delta_spin = -1d0 + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +! Build A matrix + + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + ab = ab + 1 + ij = 0 + do i=nC+1,nO + do j=nC+1,nO + ij = ij + 1 + + B_pp(ab,ij) = (1d0 + delta_spin)*ERI(a,b,i,j) - (1d0 - delta_dRPA)*ERI(a,b,j,j) + + enddo + enddo + enddo + enddo + +end subroutine linear_response_B_pp diff --git a/src/QuAcK/linear_response_C_pp.f90 b/src/QuAcK/linear_response_C_pp.f90 new file mode 100644 index 0000000..ef3b3ee --- /dev/null +++ b/src/QuAcK/linear_response_C_pp.f90 @@ -0,0 +1,57 @@ +subroutine linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C_pp) + +! Compute the C matrix of the pp channel + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: dRPA + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: delta_spin + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: a,b,c,d,ab,cd + +! Output variables + + double precision,intent(out) :: C_pp(nVV,nVV) + +! Singlet or triplet manifold? + + delta_spin = 0d0 + if(ispin == 1) delta_spin = +1d0 + if(ispin == 2) delta_spin = -1d0 + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +! Build A matrix + + ab = 0 + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + ab = ab + 1 + cd = 0 + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + cd = cd + 1 + + C_pp(ab,cd) = (e(a) + e(b))*Kronecker_delta(a,c)*Kronecker_delta(b,d) & + + (1d0 + delta_spin)*ERI(a,b,c,d) - (1d0 - delta_dRPA)*ERI(a,b,d,c) + + enddo + enddo + enddo + enddo + +end subroutine linear_response_C_pp diff --git a/src/QuAcK/linear_response_D_pp.f90 b/src/QuAcK/linear_response_D_pp.f90 new file mode 100644 index 0000000..28455b6 --- /dev/null +++ b/src/QuAcK/linear_response_D_pp.f90 @@ -0,0 +1,57 @@ +subroutine linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D_pp) + +! Compute the D matrix of the pp channel + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: dRPA + integer,intent(in) :: nBas,nC,nO,nV,nR,nOO,nVV + double precision,intent(in) :: e(nBas),ERI(nBas,nBas,nBas,nBas) + +! Local variables + + double precision :: delta_spin + double precision :: delta_dRPA + double precision,external :: Kronecker_delta + + integer :: i,j,k,l,ij,kl + +! Output variables + + double precision,intent(out) :: D_pp(nOO,nOO) + +! Singlet or triplet manifold? + + delta_spin = 0d0 + if(ispin == 1) delta_spin = +1d0 + if(ispin == 2) delta_spin = -1d0 + +! Direct RPA + + delta_dRPA = 0d0 + if(dRPA) delta_dRPA = 1d0 + +! Build A matrix + + ij = 0 + 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 + + D_pp(ij,kl) = - (e(i) + e(j))*Kronecker_delta(i,k)*Kronecker_delta(j,l) & + + (1d0 + delta_spin)*ERI(k,l,i,j) - (1d0 - delta_dRPA)*ERI(k,l,j,i) + + enddo + enddo + enddo + enddo + +end subroutine linear_response_D_pp diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 new file mode 100644 index 0000000..437a3a3 --- /dev/null +++ b/src/QuAcK/linear_response_pp.f90 @@ -0,0 +1,89 @@ +subroutine linear_response_pp(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,e,ERI,rho) + +! Compute the p-p channel of the linear response + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: dRPA + logical,intent(in) :: TDA + logical,intent(in) :: BSE + integer,intent(in) :: ispin,nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: e(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + integer :: nOO + integer :: nVV + double precision :: trace_matrix + double precision,allocatable :: B(:,:) + double precision,allocatable :: C(:,:) + double precision,allocatable :: D(:,:) + double precision,allocatable :: M(:,:) + double precision,allocatable :: w(:) + double precision :: Ec_ppRPA + +! Output variables + +! Useful quantities + + nOO = nO*nO + nVV = nV*nV + +! Memory allocation + + allocate(B(nVV,nOO),C(nVV,nVV),D(nOO,nOO),M(nOO+nVV,nOO+nVV),w(nOO+nVV)) + +! Build B, C and D matrices for the pp channel + + call linear_response_B_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,B) + call linear_response_C_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,C) + call linear_response_D_pp(ispin,dRPA,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,D) + +!------------------------------------------------------------------------ +! Solve the p-p eigenproblem +!------------------------------------------------------------------------ +! +! | C -B | | X1 X2 | | w1 0 | | X1 X2 | +! | | | | = | | | | +! | Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 | +! + +! Diagonal blocks + + M( 1:nVV , 1:nVV) = +C(1:nVV,1:nVV) + M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = -D(1:nOO,1:nOO) + +! Off-diagonal blocks + + M(1:nOO,nOO+1:nOO+nVV) = +transpose(B(1:nVV,1:nOO)) + M(nOO+1:nOO+nVV,1:nOO) = - B(1:nVV,1:nOO) + + print*, 'pp-RPA matrix' + call matout(nOO+nVV,nOO+nVV,M(:,:)) + +! Diagonalize the p-h matrix + + call diagonalize_matrix(nOO+nVV,M(:,:),w(:)) + +! Build X+Y + +! XpY(1:nS,1:nS) = M(nS+1:2*nS,1:nS) + M(nS+1:2*nS,nS+1:2*nS) + +! print*,'X+Y' +! call matout(nS,nS,XpY) + + print*,'pp-RPA excitation energies' + call matout(2*nS,1,w) + +! Compute the RPA correlation energy + + Ec_ppRPA = 0.5d0*(sum(abs(w)) - trace_matrix(nVV,C) - trace_matrix(nOO,D)) + + print*,'Ec(pp-RPA) = ',Ec_ppRPA + +end subroutine linear_response_pp diff --git a/src/QuAcK/mo_guess.f90 b/src/QuAcK/mo_guess.f90 new file mode 100644 index 0000000..9321d13 --- /dev/null +++ b/src/QuAcK/mo_guess.f90 @@ -0,0 +1,51 @@ +subroutine mo_guess(nBas,nO,guess_type,S,Hc,ERI,J,K,X,cp,Fp,e,c,P) + +! Guess of the molecular orbitals for HF calculation + + implicit none + +! Input variables + + integer,intent(in) :: nBas + integer,intent(in) :: nO + integer,intent(in) :: guess_type + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(inout):: J(nBas,nBas) + double precision,intent(inout):: K(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(inout):: cp(nBas,nBas) + double precision,intent(inout):: Fp(nBas,nBas) + double precision,intent(inout):: e(nBas) + double precision,intent(inout):: P(nBas,nBas) + +! Local variables + + integer :: nSCF + +! Output variables + + double precision,intent(out) :: c(nBas,nBas) + + if(guess_type == 1) then + + Fp = matmul(transpose(X),matmul(Hc,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,e) + c = matmul(X,cp) + + elseif(guess_type == 2) then + + call huckel_guess(nBas,nO,S,Hc,ERI,J,K,X,cp,Fp,e,c,P) + + elseif(guess_type == 3) then + + call random_number(c) + + endif + + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + + +end subroutine