diff --git a/input/options b/input/options index 62bd5b2..53c481b 100644 --- a/input/options +++ b/input/options @@ -1,5 +1,5 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess level_shift stability - 512 0.0000001 T 5 2 1 F 0.0 F + 512 0.0000001 T 5 2 1 F 1.0 F # MP: # CC: maxSCF thresh DIIS n_diis diff --git a/mol/formaldehyde.xyz b/mol/formaldehyde.xyz index 16d1880..5fdafd5 100644 --- a/mol/formaldehyde.xyz +++ b/mol/formaldehyde.xyz @@ -1,6 +1,6 @@ 4 -Formaldehyde,^1A_1,CC3,aug-cc-pVTZ + C 0.00000000 0.00000000 -0.60298508 O 0.00000000 0.00000000 0.60539399 H 0.00000000 0.93467313 -1.18217476 -H 0.00000000 -0.93467313 -1.18217476 \ No newline at end of file +H 0.00000000 -0.93467313 -1.18217476 diff --git a/mol/h2.xyz b/mol/h2.xyz index fc52308..2a6dce9 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 H 0. 0. 0. -H 0. 0. 1.0 +H 0. 0. 0.7 diff --git a/src/GW/Bethe_Salpeter_pp.f90 b/src/GW/Bethe_Salpeter_pp.f90 index 5e46305..c592269 100644 --- a/src/GW/Bethe_Salpeter_pp.f90 +++ b/src/GW/Bethe_Salpeter_pp.f90 @@ -94,8 +94,8 @@ subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,E call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin)) - call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) - call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) +! call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) +! call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) @@ -128,8 +128,8 @@ subroutine Bethe_Salpeter_pp(TDA_W,TDA,singlet,triplet,eta,nBas,nC,nO,nV,nR,nS,E call linear_response_pp_BSE(ispin,TDA,.true.,nBas,nC,nO,nV,nR,nOO,nVV,1d0,eGW,ERI,WB,WC,WD, & Omega1,X1,Y1,Omega2,X2,Y2,EcBSE(ispin)) - call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) - call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) +! call print_excitation('pp-BSE (N+2)',ispin,nVV,Omega1) +! call print_excitation('pp-BSE (N-2)',ispin,nOO,Omega2) call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) diff --git a/src/GW/static_screening_WB_pp.f90 b/src/GW/static_screening_WB_pp.f90 index 07047c8..9cb9b89 100644 --- a/src/GW/static_screening_WB_pp.f90 +++ b/src/GW/static_screening_WB_pp.f90 @@ -24,6 +24,7 @@ subroutine static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E ! Local variables + double precision,external :: Kronecker_delta double precision :: chi double precision :: eps integer :: a,b,i,j,ab,ij,m @@ -57,7 +58,7 @@ subroutine static_screening_WB_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E chi = chi + rho(a,j,m)*rho(b,i,m)*Omega(m)/eps enddo - WB(ab,ij) = WB(ab,ij) - 4d0*lambda*chi + WB(ab,ij) = WB(ab,ij) + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(i,j))) end do end do diff --git a/src/GW/static_screening_WC_pp.f90 b/src/GW/static_screening_WC_pp.f90 index 7422f76..da10b25 100644 --- a/src/GW/static_screening_WC_pp.f90 +++ b/src/GW/static_screening_WC_pp.f90 @@ -24,6 +24,7 @@ subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E ! Local variables + double precision,external :: Kronecker_delta double precision :: chi double precision :: eps integer :: a,b,c,d,ab,cd,m @@ -57,8 +58,7 @@ subroutine static_screening_WC_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E chi = chi + rho(a,d,m)*rho(b,c,m)*Omega(m)/eps enddo - WC(ab,cd) = WC(ab,cd) - 4d0*lambda*chi - + WC(ab,cd) = WC(ab,cd) + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do diff --git a/src/GW/static_screening_WD_pp.f90 b/src/GW/static_screening_WD_pp.f90 index ae4ffef..80d4178 100644 --- a/src/GW/static_screening_WD_pp.f90 +++ b/src/GW/static_screening_WD_pp.f90 @@ -24,6 +24,7 @@ subroutine static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E ! Local variables + double precision,external :: Kronecker_delta double precision :: chi double precision :: eps integer :: i,j,k,l,ij,kl,m @@ -57,7 +58,7 @@ subroutine static_screening_WD_pp(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,E chi = chi + rho(i,l,m)*rho(j,k,m)*Omega(m)/eps enddo - WD(ij,kl) = WD(ij,kl) - 4d0*lambda*chi + WD(ij,kl) = WD(ij,kl) + 4d0*lambda*chi/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do diff --git a/src/LR/print_transition_vectors_pp.f90 b/src/LR/print_transition_vectors_pp.f90 index 620ed28..0186c6e 100644 --- a/src/LR/print_transition_vectors_pp.f90 +++ b/src/LR/print_transition_vectors_pp.f90 @@ -27,24 +27,24 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip integer :: a,b,c,d,ab,cd integer :: i,j,k,l,ij,kl - integer :: maxOO + integer :: maxOO = 10 integer :: maxVV = 10 double precision :: S2 double precision,parameter :: thres_vec = 0.1d0 - double precision,allocatable :: osOO(:) - double precision,allocatable :: osVV(:) + double precision,allocatable :: os1(:) + double precision,allocatable :: os2(:) ! Memory allocation maxOO = min(nOO,maxOO) maxVV = min(nVV,maxVV) - allocate(osOO(maxOO),osVV(maxVV)) + allocate(os1(nVV),os2(nOO)) ! Compute oscillator strengths - osOO(:) = 0d0 - osVV(:) = 0d0 + os1(:) = 0d0 + os2(:) = 0d0 ! if(spin_allowed) call oscillator_strength(nBas,nC,nO,nV,nR,nS,maxS,dipole_int,Omega,XpY,XmY,os) @@ -63,8 +63,8 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip end if print*,'-------------------------------------------------------------' - write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' p-p excitation n. ',ab,': ',Omega1(ab)*HaToeV,' eV',' f = ',osVV(ab),' = ',S2 + write(*,'(A20,I3,A2,F10.4,A3,A6,F6.4,A11,F6.4)') & + ' p-p excitation n. ',ab,': ',Omega1(ab)*HaToeV,' eV',' f = ',os1(ab),' = ',S2 print*,'-------------------------------------------------------------' if(spin_allowed) then @@ -111,14 +111,14 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip ! Thomas-Reiche-Kuhn sum rule - write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(osVV(:)) + write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for p-p sector = ',sum(os1(:)) write(*,*) !-----------------------------------------------! ! Print details about excitations for hh sector ! !-----------------------------------------------! - do ij=1,maxOO + do ij=nOO,nOO-maxOO+1,-1 ! values @@ -129,8 +129,8 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip end if print*,'-------------------------------------------------------------' - write(*,'(A20,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & - ' h-h excitation n. ',ij,': ',Omega2(ij)*HaToeV,' eV',' f = ',osOO(ij),' = ',S2 + write(*,'(A20,I3,A2,F10.4,A3,A6,F6.4,A11,F6.4)') & + ' h-h excitation n. ',ij,': ',Omega2(ij)*HaToeV,' eV',' f = ',os2(ij),' = ',S2 print*,'-------------------------------------------------------------' if(spin_allowed) then @@ -177,7 +177,7 @@ subroutine print_transition_vectors_pp(spin_allowed,nBas,nC,nO,nV,nR,nOO,nVV,dip ! Thomas-Reiche-Kuhn sum rule - write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(osOO(:)) + write(*,'(A50,F10.6)') 'Thomas-Reiche-Kuhn sum rule for h-h sector = ',sum(os2(:)) write(*,*) end subroutine print_transition_vectors_pp diff --git a/src/RPA/ppRPA.f90 b/src/RPA/ppRPA.f90 index 2e6894f..d2d420b 100644 --- a/src/RPA/ppRPA.f90 +++ b/src/RPA/ppRPA.f90 @@ -69,8 +69,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin)) - call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) - call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) +! call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) +! call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) call print_transition_vectors_pp(.true.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2) @@ -92,8 +92,8 @@ subroutine ppRPA(TDA,doACFDT,singlet,triplet,nBas,nC,nO,nV,nR,ENuc,ERHF,ERI,dipo call linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,1d0,e,ERI, & Omega1,X1,Y1,Omega2,X2,Y2,Ec_ppRPA(ispin)) - call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) - call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) +! call print_excitation('pp-RPA (N+2)',ispin,nVV,Omega1) +! call print_excitation('pp-RPA (N-2)',ispin,nOO,Omega2) call print_transition_vectors_pp(.false.,nBas,nC,nO,nV,nR,nOO,nVV,dipole_int,Omega1,X1,Y1,Omega2,X2,Y2)