diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 new file mode 100644 index 0000000..a811843 --- /dev/null +++ b/src/MBPT/unrestricted_Bethe_Salpeter_ZA_matrix_dynamic.f90 @@ -0,0 +1,176 @@ +subroutine unrestricted_Bethe_Salpeter_ZA_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,lambda,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,OmBSE,ZA_dyn) + +! Compute the extra term for dynamical Bethe-Salpeter equation for linear response in the unrestricted formalism + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nSa + integer,intent(in) :: nSb + integer,intent(in) :: nSt + integer,intent(in) :: nS_sc + double precision,intent(in) :: eta + double precision,intent(in) :: lambda + double precision,intent(in) :: eGW(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: OmRPA(nS_sc) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin) + double precision,intent(in) :: OmBSE + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,kc + +! Output variables + + double precision,intent(out) :: ZA_dyn(nSt,nSt) + +!--------------------------------------------------! +! Build BSE matrix for spin-conserving transitions ! +!--------------------------------------------------! + + ZA_dyn(:,:) = 0d0 + + if(ispin == 1) then + + ! aaaa block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,1)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,1)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi + + enddo + enddo + enddo + enddo + + ! bbbb block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,2)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,2)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi + + enddo + enddo + enddo + enddo + + end if + +!--------------------------------------------! +! Build BSE matrix for spin-flip transitions ! +!--------------------------------------------! + + if(ispin == 2) then + + ! abab block + + ia = 0 + do i=nC(1)+1,nO(1) + do a=nO(2)+1,nBas-nR(2) + ia = ia + 1 + jb = 0 + do j=nC(1)+1,nO(1) + do b=nO(2)+1,nBas-nR(2) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,1) - eGW(j,2)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,1) - eGW(i,2)) + chi = chi + rho_RPA(i,j,kc,1)*rho_RPA(a,b,kc,2)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(ia,jb) = ZA_dyn(ia,jb) + lambda*chi + + end do + end do + end do + end do + + ! baba block + + ia = 0 + do i=nC(2)+1,nO(2) + do a=nO(1)+1,nBas-nR(1) + ia = ia + 1 + jb = 0 + do j=nC(2)+1,nO(2) + do b=nO(1)+1,nBas-nR(1) + jb = jb + 1 + + chi = 0d0 + do kc=1,nS_sc + + eps = + OmBSE - OmRPA(kc) - (eGW(a,2) - eGW(j,1)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + eps = + OmBSE - OmRPA(kc) - (eGW(b,2) - eGW(i,1)) + chi = chi + rho_RPA(i,j,kc,2)*rho_RPA(a,b,kc,1)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + + enddo + + ZA_dyn(nSa+ia,nSa+jb) = ZA_dyn(nSa+ia,nSa+jb) + lambda*chi + + end do + end do + end do + end do + + end if + +end subroutine unrestricted_Bethe_Salpeter_ZA_matrix_dynamic diff --git a/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 b/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 new file mode 100644 index 0000000..1057eee --- /dev/null +++ b/src/MBPT/unrestricted_Bethe_Salpeter_dynamic_perturbation.f90 @@ -0,0 +1,110 @@ +subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,nS_sc,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,dipole_int_aa,dipole_int_bb, & + OmRPA,rho_RPA,OmBSE,XpY_BSE,XmY_BSE) + +! Compute dynamical effects via perturbation theory for BSE + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: ispin + logical,intent(in) :: dTDA + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC(nspin) + integer,intent(in) :: nO(nspin) + integer,intent(in) :: nV(nspin) + integer,intent(in) :: nR(nspin) + integer,intent(in) :: nS(nspin) + integer,intent(in) :: nSa + integer,intent(in) :: nSb + integer,intent(in) :: nSt + integer,intent(in) :: nS_sc + + double precision,intent(in) :: eGW(nBas,nspin) + double precision,intent(in) :: ERI_aaaa(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_aabb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_bbbb(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_abab(nBas,nBas,nBas,nBas) + double precision,intent(in) :: dipole_int_aa(nBas,nBas,ncart) + double precision,intent(in) :: dipole_int_bb(nBas,nBas,ncart) + double precision,intent(in) :: OmRPA(nS_sc) + double precision,intent(in) :: rho_RPA(nBas,nBas,nS_sc,nspin) + double precision,intent(in) :: OmBSE(nSt) + double precision,intent(in) :: XpY_BSE(nSt,nSt) + double precision,intent(in) :: XmY_BSE(nSt,nSt) + +! Local variables + + integer :: ia + + integer,parameter :: maxS = 10 + double precision :: gapGW + + double precision,allocatable :: OmDyn(:) + double precision,allocatable :: ZDyn(:) + double precision,allocatable :: X(:) + double precision,allocatable :: Y(:) + + double precision,allocatable :: A_dyn(:,:) + double precision,allocatable :: ZA_dyn(:,:) + +! Memory allocation + + allocate(OmDyn(nSt),ZDyn(nSt),X(nSt),Y(nSt),A_dyn(nSt,nSt),ZA_dyn(nSt,nSt)) + +! Print main components of transition vectors + + if(dTDA) then + write(*,*) + write(*,*) '*** dynamical TDA activated ***' + write(*,*) + end if + + gapGW = min(eGW(nO(1)+1,1),eGW(nO(2)+1,2)) - max(eGW(nO(1),1),eGW(nO(2),2)) + + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) ' First-order dynamical correction to static Bethe-Salpeter excitation energies ' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(A57,F10.6,A3)') ' BSE neutral excitation must be lower than the GW gap = ',gapGW*HaToeV,' eV' + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' + write(*,*) '---------------------------------------------------------------------------------------------------' + + do ia=1,min(nSt,maxS) + + X(:) = 0.5d0*(XpY_BSE(ia,:) + XmY_BSE(ia,:)) + Y(:) = 0.5d0*(XpY_BSE(ia,:) - XmY_BSE(ia,:)) + + ! First-order correction + + if(dTDA) then + + ! Resonant part of the BSE correction for dynamical TDA + + call unrestricted_Bethe_Salpeter_A_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,1d0,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,OmBSE(ia),A_dyn) + + ! Renormalization factor of the resonant parts for dynamical TDA + + call unrestricted_Bethe_Salpeter_ZA_matrix_dynamic(ispin,eta,nBas,nC,nO,nV,nR,nSa,nSb,nSt,nS_sc,1d0,eGW, & + ERI_aaaa,ERI_aabb,ERI_bbbb,ERI_abab,OmRPA,rho_RPA,OmBSE(ia),ZA_dyn) + + ZDyn(ia) = dot_product(X,matmul(ZA_dyn,X)) + OmDyn(ia) = dot_product(X,matmul( A_dyn,X)) + + end if + + ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) + OmDyn(ia) = ZDyn(ia)*OmDyn(ia) + + write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & + ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia) + + end do + write(*,*) '---------------------------------------------------------------------------------------------------' + write(*,*) + +end subroutine unrestricted_Bethe_Salpeter_dynamic_perturbation diff --git a/src/Makefile.common b/src/Makefile.common new file mode 100644 index 0000000..2f08f04 --- /dev/null +++ b/src/Makefile.common @@ -0,0 +1,12 @@ +IDIR =../../include +MODDIR=../modules + +FC = gfortran -I$(IDIR) -J$(MODDIR) -I$(MODDIR) +#FC = ifort -I$(IDIR) -module $(MODDIR) +AR = libtool +ifeq ($(DEBUG),1) +FFLAGS = -Wall -g -msse4.2 -fcheck=all -Waliasing -Wampersand -Wconversion -Wsurprising -Wintrinsics-std -Wno-tabs -Wintrinsic-shadow -Wline-truncation -Wreal-q-constant +else +FFLAGS = -Wall -Wno-unused -Wno-unused-dummy-argument -O3 +endif + diff --git a/src/modules/.gitignore b/src/modules/.gitignore new file mode 100644 index 0000000..e69de29