diff --git a/examples/molecule.Be b/examples/molecule.Be index 6a6f6d1..7f85101 100644 --- a/examples/molecule.Be +++ b/examples/molecule.Be @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 2 2 0 0 + 1 3 1 0 0 # Znuc x y z Be 0.0 0.0 0.0 diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 7225285..d335f5d 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -1,5 +1,5 @@ # nAt nEla nElb nCore nRyd 2 1 1 0 0 # Znuc x y z - H 0. 0. -0.7 - H 0. 0. 0.7 + H 0. 0. 0.0 + H 0. 0. 1.4 diff --git a/examples/molecule.He b/examples/molecule.He index c78e87e..1ac370f 100644 --- a/examples/molecule.He +++ b/examples/molecule.He @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 2 0 0 0 # Znuc x y z He 0.0 0.0 0.0 diff --git a/input/methods b/input/methods index 76ead5a..3a1c5a4 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF MOM - T F F + F T F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) @@ -7,7 +7,7 @@ # drCCD rCCD lCCD pCCD F F F F # CIS* CIS(D) CID CISD - T T F F + F F F F # RPA* RPAx* ppRPA F F F # G0F2 evGF2 G0F3 evGF3 diff --git a/input/options b/input/options index 488a820..4b6f510 100644 --- a/input/options +++ b/input/options @@ -5,14 +5,14 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 # spin: singlet triplet spin_conserved spin_flip TDA - T T T T T + T T T F T # GF: maxSCF thresh DIIS n_diis lin eta renorm 256 0.00001 T 5 T 0.0 3 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 256 0.00001 T 5 T 0.0 F F F F F # ACFDT: AC Kx XBS - T F T + F F T # BSE: BSE dBSE dTDA evDyn - T F T F + F F T F # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/CI/CIS.f90 b/src/CI/CIS.f90 index 80bf16c..b4bf0fa 100644 --- a/src/CI/CIS.f90 +++ b/src/CI/CIS.f90 @@ -55,6 +55,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) call diagonalize_matrix(nS,A,Omega) call print_excitation('CIS ',ispin,nS,Omega) + call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A)) if(dump_trans) then print*,'Singlet CIS transition vectors' @@ -82,6 +83,7 @@ subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) call diagonalize_matrix(nS,A,Omega) call print_excitation('CIS ',ispin,nS,Omega) + call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int,Omega,transpose(A),transpose(A)) if(dump_trans) then print*,'Triplet CIS transition vectors' diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index e4e01d9..cc06f85 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -80,6 +80,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E call diagonalize_matrix(nS_sc,A_sc,Omega_sc) call print_excitation('UCIS ',5,nS_sc,Omega_sc) + call print_unrestricted_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,nS_aa,nS_bb,nS_sc,dipole_int_aa,dipole_int_bb, & + Omega_sc,transpose(A_sc),transpose(A_sc)) if(dump_trans) then print*,'Spin-conserved CIS transition vectors' @@ -122,6 +124,8 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E call quick_sort(Omega_sf,order(:),nS_sf) call print_excitation('UCIS ',6,nS_sf,Omega_sf) + call print_unrestricted_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,nS_ab,nS_ba,nS_sf,dipole_int_aa,dipole_int_bb, & + Omega_sf,transpose(A_sf),transpose(A_sf)) if(dump_trans) then print*,'Spin-flip CIS transition vectors' diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 34eb9ef..ac0e9c5 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -233,6 +233,6 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,nBas,nO,S,T,V,Hc,ERI,X,ENuc,EUH ! Compute final UHF energy - call print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) + call print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF) end subroutine UHF diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index e70b5b6..f5e9eb5 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -1,4 +1,4 @@ -subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) +subroutine print_UHF(nBas,nO,S,e,c,ENuc,ET,EV,EJ,Ex,EUHF) ! Print one- and two-electron energies and other stuff for UHF calculation @@ -7,6 +7,7 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) integer,intent(in) :: nBas integer,intent(in) :: nO(nspin) + double precision,intent(in) :: S(nBas,nBas) double precision,intent(in) :: e(nBas,nspin) double precision,intent(in) :: c(nBas,nBas,nspin) double precision,intent(in) :: ENuc @@ -16,10 +17,14 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) double precision,intent(in) :: Ex(nspin) double precision,intent(in) :: EUHF + integer :: i,j integer :: ispin double precision :: HOMO(nspin) double precision :: LUMO(nspin) double precision :: Gap(nspin) + double precision :: S2_exact + double precision :: S2 + integer :: spin_state ! HOMO and LUMO @@ -39,6 +44,11 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) end if end do + S2_exact = dble(nO(1) - nO(2))/2d0*(dble(nO(1) - nO(2))/2d0 + 1d0) + S2 = S2_exact + nO(2) - sum(matmul(transpose(c(:,1:nO(1),1)),matmul(S,c(:,1:nO(2),2)))**2) + + spin_state = nO(1) - nO(2) + 1 + ! Dump results write(*,*) @@ -79,6 +89,10 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) write(*,'(A40,F13.6,A3)') ' UHF LUMO b energy:',LUMO(2)*HatoeV,' eV' write(*,'(A40,F13.6,A3)') ' UHF HOMOb-LUMOb gap :',Gap(2)*HatoeV,' eV' write(*,'(A60)') '-------------------------------------------------' + write(*,'(A40,I6)') ' 2S+1 :',spin_state + write(*,'(A40,F13.6)') ' (exact) :',S2_exact + write(*,'(A40,F13.6)') ' :',S2 + write(*,'(A60)') '-------------------------------------------------' write(*,*) ! Print results diff --git a/src/RPA/ACFDT.f90 b/src/RPA/ACFDT.f90 index cd8843b..9853403 100644 --- a/src/RPA/ACFDT.f90 +++ b/src/RPA/ACFDT.f90 @@ -94,7 +94,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB if(doXBS) then - call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) @@ -143,7 +143,7 @@ subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA_W,TDA,BSE,singlet,triplet,eta,nB if(doXBS) then - call linear_response(isp_W,dRPA,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & + call linear_response(isp_W,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS,lambda,eW,ERI,OmRPA, & rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY_RPA,rho_RPA) diff --git a/src/RPA/print_transition_vectors.f90 b/src/RPA/print_transition_vectors.f90 index 605f439..b7f6e13 100644 --- a/src/RPA/print_transition_vectors.f90 +++ b/src/RPA/print_transition_vectors.f90 @@ -108,8 +108,12 @@ subroutine print_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,nS,dipole_int, end do write(*,*) + print*,2d0*sum(X(:)**2 + Y(:)**2) + end do +! Thomas-Reiche-Kuhn sum rule + write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:)) write(*,*) diff --git a/src/RPA/print_unrestricted_transition_vectors.f90 b/src/RPA/print_unrestricted_transition_vectors.f90 index 862b5b9..ab9edfc 100644 --- a/src/RPA/print_unrestricted_transition_vectors.f90 +++ b/src/RPA/print_unrestricted_transition_vectors.f90 @@ -31,7 +31,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n integer :: ixyz integer :: ispin integer,parameter :: maxS = 10 - double precision :: norm + double precision :: S2 double precision,parameter :: thres_vec = 0.1d0 double precision,allocatable :: X(:) double precision,allocatable :: Y(:) @@ -103,9 +103,20 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - print*,'---------------------------------------------' - write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A1)') ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' (f = ',os(ia),')' - print*,'---------------------------------------------' + S2 = (nO(1) - nO(2))/2d0 + S2 = 2d0*S2+1d0 + S2 = 0.0d0 + do jb=1,nSa + S2 = S2 + 4d0*(X(jb)**2 + Y(jb)**2) + end do + do jb=1,nSb + S2 = S2 - 4d0*(X(nSa+jb)**2 + Y(nSa+jb)**2) + end do + + print*,'-------------------------------------------------------------' + write(*,'(A15,I3,A2,F10.6,A3,A6,F6.4,A11,F6.4)') & + ' Excitation n. ',ia,': ',Omega(ia)*HaToeV,' eV',' f = ',os(ia),' = ',S2 + print*,'-------------------------------------------------------------' ! Spin-up transitions @@ -113,7 +124,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(1)+1,nO(1) do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb)/sqrt(2d0) + if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A -> ',b,'A = ',X(jb) end do end do @@ -121,7 +132,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(1)+1,nO(1) do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb)/sqrt(2d0) + if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'A <- ',b,'A = ',Y(jb) end do end do @@ -131,7 +142,7 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(2)+1,nO(2) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - if(abs(X(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(jb)/sqrt(2d0) + if(abs(X(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B -> ',b,'B = ',X(nSa+jb) end do end do @@ -139,13 +150,15 @@ subroutine print_unrestricted_transition_vectors(spin_allowed,nBas,nC,nO,nV,nR,n do j=nC(2)+1,nO(2) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - if(abs(Y(jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(jb)/sqrt(2d0) + if(abs(Y(nSa+jb)) > thres_vec) write(*,'(I3,A5,I3,A4,F10.6)') j,'B <- ',b,'B = ',Y(nSa+jb) end do end do write(*,*) end do +! Thomas-Reiche-Kuhn sum rule + write(*,'(A30,F10.6)') 'Thomas-Reiche-Kuhn sum rule = ',sum(os(:)) write(*,*) diff --git a/src/RPA/unrestricted_linear_response_A_matrix.f90 b/src/RPA/unrestricted_linear_response_A_matrix.f90 index 79a9a26..6a49cc9 100644 --- a/src/RPA/unrestricted_linear_response_A_matrix.f90 +++ b/src/RPA/unrestricted_linear_response_A_matrix.f90 @@ -145,8 +145,8 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(2)+1,nBas-nR(2) jb = jb + 1 ! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B)' - A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) !& -! - (1d0 - delta_dRPA)*lambda*ERI_abab(a,j,b,i) + A_lr(ia,jb) = (e(a,2) - e(i,1))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + - (1d0 - delta_dRPA)*lambda*ERI_abab(a,j,b,i) end do end do @@ -164,9 +164,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do b=nO(1)+1,nBas-nR(1) jb = jb + 1 - A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) !& -! - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) - print*,'(',i,'A',a,'B) -> (',j,'A',b,'B) -> ',A_lr(nSa+ia,nSa+jb) + A_lr(nSa+ia,nSa+jb) = (e(a,1) - e(i,2))*Kronecker_delta(i,j)*Kronecker_delta(a,b) & + - (1d0 - delta_dRPA)*lambda*ERI_abab(i,b,j,a) +! print*,'(',i,'A',a,'B) -> (',j,'A',b,'B) -> ',A_lr(nSa+ia,nSa+jb) end do end do