diff --git a/GoQCaml b/GoQCaml index 57b0d2c..5ab2190 100755 --- a/GoQCaml +++ b/GoQCaml @@ -1,5 +1,5 @@ #! /bin/bash cd int -###../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -m 0.5 +../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz +###../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz -m 0.5 diff --git a/input/basis b/input/basis index 6796e3b..7d8f45c 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,83 @@ -1 3 +1 5 S 3 - 1 38.3600000 0.0238090 - 2 5.7700000 0.1548910 - 3 1.2400000 0.4699870 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.2976000 1.0000000 + 1 0.1220000 1.0000000 +S 1 + 1 0.0297400 1.0000000 P 1 - 1 1.2750000 1.0000000 + 1 0.7270000 1.0000000 +P 1 + 1 0.1410000 1.0000000 +2 9 +S 8 + 1 9046.0000000 0.0007000 + 2 1357.0000000 0.0053890 + 3 309.3000000 0.0274060 + 4 87.7300000 0.1032070 + 5 28.5600000 0.2787230 + 6 10.2100000 0.4485400 + 7 3.8380000 0.2782380 + 8 0.7466000 0.0154400 +S 8 + 1 9046.0000000 -0.0001530 + 2 1357.0000000 -0.0012080 + 3 309.3000000 -0.0059920 + 4 87.7300000 -0.0245440 + 5 28.5600000 -0.0674590 + 6 10.2100000 -0.1580780 + 7 3.8380000 -0.1218310 + 8 0.7466000 0.5490030 +S 1 + 1 0.2248000 1.0000000 +S 1 + 1 0.0612400 1.0000000 +P 3 + 1 13.5500000 0.0399190 + 2 2.9170000 0.2171690 + 3 0.7973000 0.5103190 +P 1 + 1 0.2185000 1.0000000 +P 1 + 1 0.0561100 1.0000000 +D 1 + 1 0.8170000 1.0000000 +D 1 + 1 0.2300000 1.0000000 +3 9 +S 8 + 1 11720.0000000 0.0007100 + 2 1759.0000000 0.0054700 + 3 400.8000000 0.0278370 + 4 113.7000000 0.1048000 + 5 37.0300000 0.2830620 + 6 13.2700000 0.4487190 + 7 5.0250000 0.2709520 + 8 1.0130000 0.0154580 +S 8 + 1 11720.0000000 -0.0001600 + 2 1759.0000000 -0.0012630 + 3 400.8000000 -0.0062670 + 4 113.7000000 -0.0257160 + 5 37.0300000 -0.0709240 + 6 13.2700000 -0.1654110 + 7 5.0250000 -0.1169550 + 8 1.0130000 0.5573680 +S 1 + 1 0.3023000 1.0000000 +S 1 + 1 0.0789600 1.0000000 +P 3 + 1 17.7000000 0.0430180 + 2 3.8540000 0.2289130 + 3 1.0460000 0.5087280 +P 1 + 1 0.2753000 1.0000000 +P 1 + 1 0.0685600 1.0000000 +D 1 + 1 1.1850000 1.0000000 +D 1 + 1 0.3320000 1.0000000 diff --git a/input/methods b/input/methods index ff1d101..63c03d4 100644 --- a/input/methods +++ b/input/methods @@ -7,13 +7,13 @@ # drCCD rCCD lCCD pCCD F F F F # CIS CID CISD - F F T + F F F # RPA RPAx ppRPA F F F # G0F2 evGF2 G0F3 evGF3 F F F F # G0W0 evGW qsGW - F F F + T F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/molecule b/input/molecule index c78e87e..634523f 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,6 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 3 8 8 0 0 # Znuc x y z - He 0.0 0.0 0.0 +H 1.18163475 0.00000000 -1.17386890 +N -0.44776863 0.00000000 -0.03589263 +O 0.21099695 0.00000000 2.15462460 diff --git a/input/molecule.xyz b/input/molecule.xyz index 797b5fc..38e47cd 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,5 @@ - 1 + 3 - He 0.0000000000 0.0000000000 0.0000000000 + H 0.6252942263 0.0000000000 -0.6211847152 + N -0.2369489718 0.0000000000 -0.0189935632 + O 0.1116547855 0.0000000000 1.1401783185 diff --git a/input/options b/input/options index ee86c35..8771279 100644 --- a/input/options +++ b/input/options @@ -4,12 +4,12 @@ # CC: maxSCF thresh DIIS n_diis 64 0.0000001 T 5 -# CIS/TDHF/BSE: singlet triplet - T F +# spin: singlet triplet + T T # 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 T 0.000 + 256 0.00001 T 5 F F T F F F T 0.000 # ACFDT: AC Kx XBS F F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift diff --git a/input/weight b/input/weight index 6796e3b..7d8f45c 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,83 @@ -1 3 +1 5 S 3 - 1 38.3600000 0.0238090 - 2 5.7700000 0.1548910 - 3 1.2400000 0.4699870 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.2976000 1.0000000 + 1 0.1220000 1.0000000 +S 1 + 1 0.0297400 1.0000000 P 1 - 1 1.2750000 1.0000000 + 1 0.7270000 1.0000000 +P 1 + 1 0.1410000 1.0000000 +2 9 +S 8 + 1 9046.0000000 0.0007000 + 2 1357.0000000 0.0053890 + 3 309.3000000 0.0274060 + 4 87.7300000 0.1032070 + 5 28.5600000 0.2787230 + 6 10.2100000 0.4485400 + 7 3.8380000 0.2782380 + 8 0.7466000 0.0154400 +S 8 + 1 9046.0000000 -0.0001530 + 2 1357.0000000 -0.0012080 + 3 309.3000000 -0.0059920 + 4 87.7300000 -0.0245440 + 5 28.5600000 -0.0674590 + 6 10.2100000 -0.1580780 + 7 3.8380000 -0.1218310 + 8 0.7466000 0.5490030 +S 1 + 1 0.2248000 1.0000000 +S 1 + 1 0.0612400 1.0000000 +P 3 + 1 13.5500000 0.0399190 + 2 2.9170000 0.2171690 + 3 0.7973000 0.5103190 +P 1 + 1 0.2185000 1.0000000 +P 1 + 1 0.0561100 1.0000000 +D 1 + 1 0.8170000 1.0000000 +D 1 + 1 0.2300000 1.0000000 +3 9 +S 8 + 1 11720.0000000 0.0007100 + 2 1759.0000000 0.0054700 + 3 400.8000000 0.0278370 + 4 113.7000000 0.1048000 + 5 37.0300000 0.2830620 + 6 13.2700000 0.4487190 + 7 5.0250000 0.2709520 + 8 1.0130000 0.0154580 +S 8 + 1 11720.0000000 -0.0001600 + 2 1759.0000000 -0.0012630 + 3 400.8000000 -0.0062670 + 4 113.7000000 -0.0257160 + 5 37.0300000 -0.0709240 + 6 13.2700000 -0.1654110 + 7 5.0250000 -0.1169550 + 8 1.0130000 0.5573680 +S 1 + 1 0.3023000 1.0000000 +S 1 + 1 0.0789600 1.0000000 +P 3 + 1 17.7000000 0.0430180 + 2 3.8540000 0.2289130 + 3 1.0460000 0.5087280 +P 1 + 1 0.2753000 1.0000000 +P 1 + 1 0.0685600 1.0000000 +D 1 + 1 1.1850000 1.0000000 +D 1 + 1 0.3320000 1.0000000 diff --git a/src/QuAcK/BSE_dynamic_perturbation.f90 b/src/QuAcK/BSE_dynamic_perturbation.f90 new file mode 100644 index 0000000..63dc325 --- /dev/null +++ b/src/QuAcK/BSE_dynamic_perturbation.f90 @@ -0,0 +1,61 @@ +subroutine BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA,OmBSE,XpY,XmY,rho) + +! Compute dynamical effects via perturbation theory for BSE + + implicit none + include 'parameters.h' + +! Input variables + + logical,intent(in) :: TDA + double precision,intent(in) :: eta + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + + double precision :: OmRPA(nS) + double precision :: OmBSE(nS) + double precision :: XpY(nS,nS) + double precision :: XmY(nS,nS) + double precision :: rho(nBas,nBas,nS) + +! Local variables + + integer :: ia + integer,parameter :: maxS = 10 + + double precision,allocatable :: OmDyn(:) + double precision,allocatable :: X(:) + double precision,allocatable :: Y(:) + double precision,allocatable :: A_dyn(:,:) + double precision,allocatable :: B_dyn(:,:) + +! Memory allocation + + allocate(OmDyn(nS),X(nS),Y(nS),A_dyn(nS,nS),B_dyn(nS,nS)) + + write(*,*) '-----------------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A30,1X,A30,1X,A30)') '#','Static excitation (eV)','Dynamic correction (eV)','Dynamic excitation (eV)' + write(*,*) '-----------------------------------------------------------------------------------------------------------' + + do ia=1,maxS + + X(:) = 0.5d0*(XpY(:,ia) + XmY(:,ia)) + Y(:) = 0.5d0*(XpY(:,ia) - XmY(:,ia)) + + call Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(:),rho(:,:,:),A_dyn(:,:)) + + if(TDA) then + B_dyn(:,:) = 0d0 + else + call Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,1d0,OmRPA(:),OmBSE(:),rho(:,:,:),B_dyn(:,:)) + end if + + OmDyn(ia) = dot_product(X(:),matmul(A_dyn(:,:),X(:))) - dot_product(Y(:),matmul(A_dyn(:,:),Y(:))) & + + dot_product(X(:),matmul(B_dyn(:,:),Y(:))) - dot_product(Y(:),matmul(B_dyn(:,:),X(:))) + + write(*,'(2X,I5,15X,F15.6,15X,F15.6,15X,F15.6)') ia,OmBSE(ia)*HaToeV,OmDyn(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV + + end do + write(*,*) '-----------------------------------------------------------------------------------------------------------' + write(*,*) + +end subroutine BSE_dynamic_perturbation diff --git a/src/QuAcK/Bethe_Salpeter.f90 b/src/QuAcK/Bethe_Salpeter.f90 index b0932ac..ba77a3a 100644 --- a/src/QuAcK/Bethe_Salpeter.f90 +++ b/src/QuAcK/Bethe_Salpeter.f90 @@ -1,5 +1,5 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & - nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,Omega,XpY,XmY,rho,EcRPA,EcBSE) + nBas,nC,nO,nV,nR,nS,ERI,eW,eGW,OmRPA,XpY,XmY,rho,EcRPA,EcBSE) ! Compute the Bethe-Salpeter excitation energies @@ -18,7 +18,7 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & double precision,intent(in) :: eGW(nBas) double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) - double precision :: Omega(nS,nspin) + double precision :: OmRPA(nS,nspin) double precision :: XpY(nS,nS,nspin) double precision :: XmY(nS,nS,nspin) double precision :: rho(nBas,nBas,nS,nspin) @@ -26,13 +26,20 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & ! Local variables integer :: ispin + double precision,allocatable :: OmBSE(:,:) ! Output variables double precision,intent(out) :: EcRPA(nspin) double precision,intent(out) :: EcBSE(nspin) +! Memory allocation + + allocate(OmBSE(nS,nspin)) + +!------------------- ! Singlet manifold +!------------------- if(singlet_manifold) then @@ -40,16 +47,24 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & EcBSE(ispin) = 0d0 call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + OmBSE(:,ispin) = OmRPA(:,ispin) call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + + ! Compute dynamic correction for BSE via perturbation theory + + call BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & + XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin)) end if +!------------------- ! Triplet manifold +!------------------- if(triplet_manifold) then @@ -57,12 +72,19 @@ subroutine Bethe_Salpeter(TDA,singlet_manifold,triplet_manifold,eta, & EcBSE(ispin) = 0d0 call linear_response(ispin,.true.,TDA,.false.,eta,nBas,nC,nO,nV,nR,nS,1d0,eW,ERI, & - rho(:,:,:,ispin),EcRPA(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + rho(:,:,:,ispin),EcRPA(ispin),OmRPA(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin)) + OmBSE(:,ispin) = OmRPA(:,ispin) call linear_response(ispin,.true.,TDA,.true.,eta,nBas,nC,nO,nV,nR,nS,1d0,eGW,ERI, & - rho(:,:,:,ispin),EcBSE(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) - call print_excitation('BSE ',ispin,nS,Omega(:,ispin)) + rho(:,:,:,ispin),EcBSE(ispin),OmBSE(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin)) + call print_excitation('BSE ',ispin,nS,OmBSE(:,ispin)) + + ! Compute dynamic correction for BSE via perturbation theory + + call BSE_dynamic_perturbation(TDA,eta,nBas,nC,nO,nV,nR,nS,OmRPA(:,ispin),OmBSE(:,ispin), & + XpY(:,:,ispin),XmY(:,:,ispin),rho(:,:,:,ispin)) + end if diff --git a/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 new file mode 100644 index 0000000..b456b66 --- /dev/null +++ b/src/QuAcK/Bethe_Salpeter_A_matrix_dynamic.f90 @@ -0,0 +1,66 @@ +subroutine Bethe_Salpeter_A_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,OmBSE,rho,A_lr) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: OmBSE(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: A_lr(nS,nS) + + A_lr(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi = 0d0 + do kc=1,nS + + eps = OmRPA(kc)**2 + eta**2 + chi = chi + rho(i,j,kc)*rho(a,b,kc)*OmRPA(kc)/eps + + enddo + + A_lr(ia,jb) = A_lr(ia,jb) - 4d0*lambda*chi + + chi = 0d0 + do kc=1,nS + + eps = (OmBSE(kc) - OmRPA(kc))**2 + eta**2 + chi = chi + rho(i,j,kc)*rho(a,b,kc)*(OmBSE(kc) - OmRPA(kc))/eps + + eps = (OmBSE(kc) + OmRPA(kc))**2 + eta**2 + chi = chi - rho(i,j,kc)*rho(a,b,kc)*(OmBSE(kc) + OmRPA(kc))/eps + + enddo + + A_lr(ia,jb) = A_lr(ia,jb) - 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine Bethe_Salpeter_A_matrix_dynamic diff --git a/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 b/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 new file mode 100644 index 0000000..8b46d84 --- /dev/null +++ b/src/QuAcK/Bethe_Salpeter_B_matrix_dynamic.f90 @@ -0,0 +1,67 @@ +subroutine Bethe_Salpeter_B_matrix_dynamic(eta,nBas,nC,nO,nV,nR,nS,lambda,OmRPA,OmBSE,rho,B_lr) + +! Compute the dynamic part of the Bethe-Salpeter equation matrices + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas,nC,nO,nV,nR,nS + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: OmRPA(nS) + double precision,intent(in) :: OmBSE(nS) + double precision,intent(in) :: rho(nBas,nBas,nS) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: B_lr(nS,nS) + + B_lr(:,:) = 0d0 + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi = 0d0 + do kc=1,nS + + eps = OmRPA(kc)**2 + eta**2 + chi = chi + rho(i,b,kc)*rho(a,j,kc)*OmRPA(kc)/eps + + enddo + + B_lr(ia,jb) = B_lr(ia,jb) - 4d0*lambda*chi + + chi = 0d0 + do kc=1,nS + + eps = (OmBSE(kc) - OmRPA(kc))**2 + eta**2 + chi = chi + rho(i,b,kc)*rho(a,j,kc)*(OmBSE(kc) - OmRPA(kc))/eps + + eps = (OmBSE(kc) + OmRPA(kc))**2 + eta**2 + chi = chi - rho(i,b,kc)*rho(a,j,kc)*(OmBSE(kc) + OmRPA(kc))/eps + + + enddo + + B_lr(ia,jb) = B_lr(ia,jb) - 2d0*lambda*chi + + enddo + enddo + enddo + enddo + +end subroutine Bethe_Salpeter_B_matrix_dynamic diff --git a/src/QuAcK/linear_response.f90 b/src/QuAcK/linear_response.f90 index bf64262..c3ea157 100644 --- a/src/QuAcK/linear_response.f90 +++ b/src/QuAcK/linear_response.f90 @@ -41,6 +41,7 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E ! Build A and B matrices call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) + if(BSE) call Bethe_Salpeter_A_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,A) ! Tamm-Dancoff approximation @@ -49,9 +50,10 @@ subroutine linear_response(ispin,dRPA,TDA,BSE,eta,nBas,nC,nO,nV,nR,nS,lambda,e,E if(.not. TDA) then call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) + if(BSE) call Bethe_Salpeter_B_matrix(eta,nBas,nC,nO,nV,nR,nS,lambda,ERI,Omega,rho,B) - endif + end if ! Build A + B and A - B matrices