diff --git a/examples/molecule.Be b/examples/molecule.Be index 7f85101..6a6f6d1 100644 --- a/examples/molecule.Be +++ b/examples/molecule.Be @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 3 1 0 0 + 1 2 2 0 0 # Znuc x y z Be 0.0 0.0 0.0 diff --git a/examples/molecule.He b/examples/molecule.He index 1ac370f..c78e87e 100644 --- a/examples/molecule.He +++ b/examples/molecule.He @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 2 0 0 0 + 1 1 1 0 0 # Znuc x y z He 0.0 0.0 0.0 diff --git a/input/methods b/input/methods index bb07f30..76ead5a 100644 --- a/input/methods +++ b/input/methods @@ -1,13 +1,13 @@ # RHF UHF MOM - F T F + T F F # MP2* MP3 MP2-F12 F F F # CCD CCSD CCSD(T) F F F # drCCD rCCD lCCD pCCD F F F F -# CIS* CID CISD - T F F +# CIS* CIS(D) CID CISD + T T F F # RPA* RPAx* ppRPA F F F # G0F2 evGF2 G0F3 evGF3 diff --git a/src/CI/CIS.f90 b/src/CI/CIS.f90 index ac4923c..80bf16c 100644 --- a/src/CI/CIS.f90 +++ b/src/CI/CIS.f90 @@ -1,5 +1,4 @@ -subroutine CIS(singlet_manifold,triplet_manifold, & - nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) +subroutine CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI,dipole_int,eHF) ! Perform configuration interaction single calculation` @@ -8,8 +7,9 @@ subroutine CIS(singlet_manifold,triplet_manifold, & ! Input variables - logical,intent(in) :: singlet_manifold - logical,intent(in) :: triplet_manifold + logical,intent(in) :: singlet + logical,intent(in) :: triplet + logical,intent(in) :: doCIS_D integer,intent(in) :: nBas,nC,nO,nV,nR,nS double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) @@ -20,6 +20,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, & logical :: dump_matrix = .false. logical :: dump_trans = .false. integer :: ispin + integer :: maxS = 10 double precision :: lambda double precision,allocatable :: A(:,:),Omega(:) @@ -41,7 +42,7 @@ subroutine CIS(singlet_manifold,triplet_manifold, & ! Compute CIS matrix - if(singlet_manifold) then + if(singlet) then ispin = 1 call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A) @@ -61,9 +62,14 @@ subroutine CIS(singlet_manifold,triplet_manifold, & write(*,*) endif + ! Compute CIS(D) correction + + maxS = min(maxS,nS) + if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS)) + endif - if(triplet_manifold) then + if(triplet) then ispin = 2 call linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS,lambda,eHF,ERI,A) @@ -83,6 +89,11 @@ subroutine CIS(singlet_manifold,triplet_manifold, & write(*,*) endif + ! Compute CIS(D) correction + + maxS = min(maxS,nS) + if(doCIS_D) call D_correction(ispin,nBas,nC,nO,nV,nR,nS,maxS,eHF,ERI,Omega(1:maxS),A(:,1:maxS)) + endif end subroutine CIS diff --git a/src/CI/D_correction.f90 b/src/CI/D_correction.f90 new file mode 100644 index 0000000..9fb2b4b --- /dev/null +++ b/src/CI/D_correction.f90 @@ -0,0 +1,224 @@ +subroutine D_correction(ispin,nBasin,nCin,nOin,nVin,nRin,nSin,maxS,eHF,ERI,w,X) + +! Compute the D correction of CIS(D) + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBasin + integer,intent(in) :: nCin + integer,intent(in) :: nOin + integer,intent(in) :: nVin + integer,intent(in) :: nRin + integer,intent(in) :: nSin + integer,intent(in) :: maxS + double precision,intent(in) :: eHF(nBasin) + double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin) + double precision,intent(in) :: w(maxS) + double precision,intent(in) :: X(nSin,maxS) + +! Local variables + + integer :: i,j,k + integer :: a,b,c + integer :: m,ia + + integer :: nBas + integer :: nC + integer :: nO + integer :: nV + integer :: nR + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta(:,:,:,:) + + double precision,allocatable :: OOOV(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVVV(:,:,:,:) + + double precision,allocatable :: u(:,:,:,:) + double precision,allocatable :: v(:,:) + double precision,allocatable :: t(:,:,:,:) + double precision,allocatable :: rr(:,:),r(:,:) + double precision :: wD + + double precision,external :: Kronecker_delta + +! Spatial to spin orbitals + + nBas = 2*nBasin + nC = 2*nCin + nO = 2*nOin + nV = 2*nVin + nR = 2*nRin + + allocate(seHF(nBas),sERI(nBas,nBas,nBas,nBas)) + + call spatial_to_spin_MO_energy(nBasin,eHF,nBas,seHF) + call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas,nBas,nBas,nBas)) + + call antisymmetrize_ERI(2,nBas,sERI,dbERI) + + deallocate(sERI) + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta(nO,nO,nV,nV)) + + eO(:) = seHF(1:nO) + eV(:) = seHF(nO+1:nBas) + + call form_delta_OOVV(nC,nO,nV,nR,eO,eV,delta) + + deallocate(seHF,eO,eV) + +! Create integral batches + + allocate(OOOV(nO,nO,nO,nV),OOVV(nO,nO,nV,nV),OVVV(nO,nV,nV,nV)) + + OOOV(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO ,nO+1:nBas) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas) + OVVV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas,nO+1:nBas) + + deallocate(dbERI) + +! Memory allocation + + allocate(t(nO,nO,nV,nV),r(nO,nV),u(nO,nO,nV,nV),v(nO,nV)) + +! MP2 guess amplitudes + + t(:,:,:,:) = -OOVV(:,:,:,:)/delta(:,:,:,:) + +!------------------------------------------------------------------------ +! Loop over single excitations +!------------------------------------------------------------------------ + + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) ' CIS(D) correction ' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A20,1X,A20,1X,A20)') '#','CIS (eV)','CIS(D) (eV)','Correction (eV)' + write(*,*) '---------------------------------------------------------------------------------------------------' + + do m=1,maxS + + ! Unfold r + + allocate(rr(nOin,nVin)) + + ia = 0 + do i=nCin+1,nOin + do a=1,nVin-nRin + ia = ia + 1 + rr(i,a) = x(ia,m) + end do + end do + + if(ispin == 1) then + + do i=nC+1,nO + do a=1,nV-nR + r(i,a) = rr((i+1)/2,(a+1)/2)*Kronecker_delta(mod(i,2),mod(a,2)) + end do + end do + + elseif(ispin == 2) then + + do i=nC+1,nO + do a=1,nV-nR + r(i,a) = rr((i+1)/2,(a+1)/2)*Kronecker_delta(mod(i,2),mod(a+1,2)) + end do + end do + + else + + print*,'!!! CIS(D) must be for singlet or triplet !!!' + stop + + end if + + deallocate(rr) + + ! Compute u array + + u(:,:,:,:) = 0d0 + + do i=nC+1,nO + do a=1,nV-nR + do j=nC+1,nO + do b=1,nV-nR + + do c=1,nV-nR + u(i,j,a,b) = u(i,j,a,b) + OVVV(i,c,a,b)*r(j,c) - OVVV(j,c,a,b)*r(i,c) + end do + + do k=nC+1,nO + u(i,j,a,b) = u(i,j,a,b) + OOOV(i,j,k,a)*r(k,b) - OOOV(i,j,k,b)*r(k,a) + end do + + end do + end do + end do + end do + + ! Compute v array + + v(:,:) = 0d0 + + do i=nC+1,nO + do j=nC+1,nO + do k=nC+1,nO + do a=1,nV-nR + do b=1,nV-nR + do c=1,nV-nR + v(i,a) = v(i,a) + OOVV(j,k,b,c)*(r(i,b)*t(j,k,c,a) + r(j,a)*t(i,k,c,b) + 2d0*r(j,b)*t(i,k,a,c)) + end do + end do + end do + end do + end do + end do + + v(:,:) = 0.5d0*v(:,:) + + ! Compute CIS(D) correction to CIS excitation energies + + wD = 0d0 + + do i=nC+1,nO + do a=1,nV-nR + do j=nC+1,nO + do b=1,nV-nR + + wD = wD - 0.25d0*u(i,j,a,b)**2/(delta(i,j,a,b) - w(m)) + + enddo + enddo + wD = wD + r(i,a)*v(i,a) + enddo + enddo + wD = 0.5d0*wD + + ! Flush results + + write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6)') m,w(m)*HaToeV,(w(m)+wD)*HaToeV,wD*HaToeV + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) +!------------------------------------------------------------------------ +! End of loop over single excitations +!------------------------------------------------------------------------ +end subroutine D_correction diff --git a/src/CI/UCIS.f90 b/src/CI/UCIS.f90 index 5a8e7bd..e4e01d9 100644 --- a/src/CI/UCIS.f90 +++ b/src/CI/UCIS.f90 @@ -37,7 +37,9 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E integer :: nS_ab,nS_ba,nS_sf double precision,allocatable :: A_sf(:,:) + double precision,allocatable :: Z_sf(:,:) double precision,allocatable :: Omega_sf(:) + integer ,allocatable :: order(:) ! Hello world @@ -103,7 +105,7 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E nS_ba = (nO(2) - nC(2))*(nV(1) - nR(1)) nS_sf = nS_ab + nS_ba - allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf)) + allocate(A_sf(nS_sf,nS_sf),Omega_sf(nS_sf),Z_sf(nS_sf,nS_sf)) call unrestricted_linear_response_A_matrix(ispin,.false.,nBas,nC,nO,nV,nR,nS_ab,nS_ba,nS_sf,lambda,eHF, & ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,A_sf) @@ -114,7 +116,11 @@ subroutine UCIS(spin_conserved,spin_flip,nBas,nC,nO,nV,nR,nS,ERI_aaaa,ERI_aabb,E write(*,*) endif - call diagonalize_matrix(nS_sf,A_sf,Omega_sf) +! call diagonalize_matrix(nS_sf,A_sf,Omega_sf) + allocate(order(nS_sf)) + call diagonalize_general_matrix(nS_sf,A_sf,Omega_sf,Z_sf) + call quick_sort(Omega_sf,order(:),nS_sf) + call print_excitation('UCIS ',6,nS_sf,Omega_sf) if(dump_trans) then diff --git a/src/HF/print_UHF.f90 b/src/HF/print_UHF.f90 index 640d023..e70b5b6 100644 --- a/src/HF/print_UHF.f90 +++ b/src/HF/print_UHF.f90 @@ -26,7 +26,11 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) do ispin=1,nspin if(nO(ispin) > 0) then HOMO(ispin) = e(nO(ispin),ispin) - LUMO(ispin) = e(nO(ispin)+1,ispin) + if(nO(ispin) < nBas) then + LUMO(ispin) = e(nO(ispin)+1,ispin) + else + LUMO(ispin) = 0d0 + end if Gap(ispin) = LUMO(ispin) - HOMO(ispin) else HOMO(ispin) = 0d0 @@ -37,7 +41,6 @@ subroutine print_UHF(nBas,nO,e,c,ENuc,ET,EV,EJ,Ex,EUHF) ! Dump results - write(*,*) write(*,'(A60)') '-------------------------------------------------' write(*,'(A40)') ' Summary ' diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index f16eab9..ed41f7a 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -9,7 +9,7 @@ program QuAcK logical :: doMP2,doMP3,doMP2F12 logical :: doCCD,doCCSD,doCCSDT logical :: do_drCCD,do_rCCD,do_lCCD,do_pCCD - logical :: doCIS,doCID,doCISD + logical :: doCIS,doCIS_D,doCID,doCISD logical :: doRPA,doRPAx,doppRPA logical :: doADC logical :: doG0F2,doevGF2,doG0F3,doevGF3 @@ -158,7 +158,7 @@ program QuAcK doMP2,doMP3,doMP2F12, & doCCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doCID,doCISD, & + doCIS,doCIS_D,doCID,doCISD, & doRPA,doRPAx,doppRPA, & doG0F2,doevGF2,doG0F3,doevGF3, & doG0W0,doevGW,doqsGW, & @@ -621,7 +621,7 @@ program QuAcK else - call CIS(singlet,triplet,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF) + call CIS(singlet,triplet,doCIS_D,nBas,nC,nO,nV,nR,nS,ERI_MO,dipole_int,eHF) end if call cpu_time(end_CIS) diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index ec9afc1..3b4770b 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -2,7 +2,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & doMP2,doMP3,doMP2F12, & doCCD,doCCSD,doCCSDT, & do_drCCD,do_rCCD,do_lCCD,do_pCCD, & - doCIS,doCID,doCISD, & + doCIS,doCIS_D,doCID,doCISD, & doRPA,doRPAx,doppRPA, & doG0F2,doevGF2,doG0F3,doevGF3, & doG0W0,doevGW,doqsGW, & @@ -19,7 +19,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & logical,intent(out) :: doMP2,doMP3,doMP2F12 logical,intent(out) :: doCCD,doCCSD,doCCSDT logical,intent(out) :: do_drCCD,do_rCCD,do_lCCD,do_pCCD - logical,intent(out) :: doCIS,doCID,doCISD + logical,intent(out) :: doCIS,doCIS_D,doCID,doCISD logical,intent(out) :: doRPA,doRPAx,doppRPA logical,intent(out) :: doG0F2,doevGF2,doG0F3,doevGF3 logical,intent(out) :: doG0W0,doevGW,doqsGW @@ -54,6 +54,7 @@ subroutine read_methods(doRHF,doUHF,doMOM, & do_pCCD = .false. doCIS = .false. + doCIS_D = .false. doCID = .false. doCISD = .false. @@ -111,10 +112,12 @@ subroutine read_methods(doRHF,doUHF,doMOM, & ! Read excited state methods read(1,*) - read(1,*) answer1,answer2,answer3 + read(1,*) answer1,answer2,answer3,answer4 if(answer1 == 'T') doCIS = .true. - if(answer2 == 'T') doCID = .true. - if(answer3 == 'T') doCISD = .true. + if(answer2 == 'T') doCIS_D = .true. + if(answer3 == 'T') doCID = .true. + if(answer4 == 'T') doCISD = .true. + if(doCIS_D) doCIS = .true. read(1,*) read(1,*) answer1,answer2,answer3 diff --git a/src/RPA/unrestricted_linear_response_A_matrix.f90 b/src/RPA/unrestricted_linear_response_A_matrix.f90 index 8754667..79a9a26 100644 --- a/src/RPA/unrestricted_linear_response_A_matrix.f90 +++ b/src/RPA/unrestricted_linear_response_A_matrix.f90 @@ -121,6 +121,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa end do end do + print*,'nSa,nSb,nSt',nSa,nSb,nSt + call matout(nSt,nSt,A_lr) + end if !----------------------------------------------- @@ -141,9 +144,9 @@ subroutine unrestricted_linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nSa do j=nC(1)+1,nO(1) do b=nO(2)+1,nBas-nR(2) jb = jb + 1 - - A_lr(ia,jb) = (e(a,2) - e(i,1))*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(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 @@ -161,15 +164,18 @@ 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(j,a,i,b) -! print*,'ia,jb,A(ia,jb) = ',i,a,j,b,ia,jb,A_lr(ia,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 end do end do + print*,'nSa,nSb,nSt',nSa,nSb,nSt + call matout(nSt,nSt,A_lr) + end if diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 6fe4fc5..11cda87 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -107,8 +107,8 @@ program eDFT ! Allocate ensemble weights allocate(wEns(maxEns),occnum(nBas,nspin,maxEns)) - call read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & - maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice) + call read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & + maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice) !------------------------------------------------------------------------ ! Read one- and two-electron integrals diff --git a/src/eDFT/read_options.f90 b/src/eDFT/read_options_dft.f90 similarity index 95% rename from src/eDFT/read_options.f90 rename to src/eDFT/read_options_dft.f90 index cb1d05b..1448aff 100644 --- a/src/eDFT/read_options.f90 +++ b/src/eDFT/read_options_dft.f90 @@ -1,5 +1,5 @@ -subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & - maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice) +subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_w1,aCC_w2, & + maxSCF,thresh,DIIS,max_diis,guess_type,ortho_type,doNcentered,occnum,Cx_choice) ! Read DFT options @@ -201,4 +201,4 @@ subroutine read_options(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,aCC_ close(unit=1) -end subroutine read_options +end subroutine read_options_dft