LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:35:36 +01:00
Code Releases Activity

UMP2

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-09-18 13:52:35 +02:00
parent eb82dd2d67
commit 26cec28a31
12 changed files with 419 additions and 104 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
examples/molecule.C2
View File

@ -1,5 +1,5 @@
# nAt nEla nElb nCore nRyd
2 7 7 0 0
2 6 6 0 0
# Znuc x y z
C 0. 0. 0.
C 0. 0. 2.0

56
input/basis
View File

@ -1,30 +1,32 @@
1 6
S 8
1 1469.0000000 0.0007660
2 220.5000000 0.0058920
3 50.2600000 0.0296710
4 14.2400000 0.1091800
5 4.5810000 0.2827890
6 1.5800000 0.4531230
7 0.5640000 0.2747740
8 0.0734500 0.0097510
S 8
1 1469.0000000 -0.0001200
2 220.5000000 -0.0009230
3 50.2600000 -0.0046890
4 14.2400000 -0.0176820
5 4.5810000 -0.0489020
6 1.5800000 -0.0960090
7 0.5640000 -0.1363800
8 0.0734500 0.5751020
1 6
S 9
1 6.665000E+03 6.920000E-04
2 1.000000E+03 5.329000E-03
3 2.280000E+02 2.707700E-02
4 6.471000E+01 1.017180E-01
5 2.106000E+01 2.747400E-01
6 7.495000E+00 4.485640E-01
7 2.797000E+00 2.850740E-01
8 5.215000E-01 1.520400E-02
9 1.596000E-01 -3.191000E-03
S 9
1 6.665000E+03 -1.460000E-04
2 1.000000E+03 -1.154000E-03
3 2.280000E+02 -5.725000E-03
4 6.471000E+01 -2.331200E-02
5 2.106000E+01 -6.395500E-02
6 7.495000E+00 -1.499810E-01
7 2.797000E+00 -1.272620E-01
8 5.215000E-01 5.445290E-01
9 1.596000E-01 5.804960E-01
S 1
1 0.0280500 1.0000000
P 3
1 1.5340000 0.0227840
2 0.2749000 0.1391070
3 0.0736200 0.5003750
1 1.596000E-01 1.000000E+00
P 4
1 9.439000E+00 3.810900E-02
2 2.002000E+00 2.094800E-01
3 5.456000E-01 5.085570E-01
4 1.517000E-01 4.688420E-01
P 1
1 0.0240300 1.0000000
1 1.517000E-01 1.000000E+00
D 1
1 0.1239000 1.0000000
1 5.500000E-01 1.0000000

6
input/methods
View File

@ -1,7 +1,7 @@
# RHF UHF MOM
T F F
F T F
# MP2 MP3 MP2-F12
F F F
T F F
# CCD CCSD CCSD(T)
F F F
# drCCD rCCD lCCD pCCD
@ -13,7 +13,7 @@
# G0F2 evGF2 G0F3 evGF3
F F F F
# G0W0 evGW qsGW
T F F
F F F
# G0T0 evGT qsGT
F F F
# MCMP2

4
input/molecule
View File

@ -1,4 +1,4 @@
# nAt nEla nElb nCore nRyd
1 2 1 0 0
1 5 1 0 0
# Znuc x y z
Li 0.0 0.0 0.0
C 0. 0. 0.

2
input/molecule.xyz
View File

@ -1,3 +1,3 @@
1
Li 0.0000000000 0.0000000000 0.0000000000
C 0.0000000000 0.0000000000 0.0000000000

4
input/options
View File

@ -7,9 +7,9 @@
# spin: singlet triplet TDA
T T F
# GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.00367493 3
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.00367493 F F F F F
256 0.00001 T 5 T 0.0 F F F F F
# ACFDT: AC Kx XBS
F F T
# BSE: BSE dBSE dTDA evDyn

22
src/QuAcK/AOtoMO_integral_transform.f90
View File

@ -1,14 +1,17 @@
subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
subroutine AOtoMO_integral_transform(bra,ket,nBas,c,ERI_AO_basis,ERI_MO_basis)
! AO to MO transformation of two-electron integrals
! Semi-direct O(N^5) algorithm
! AO to MO transformation of two-electron integrals via the semi-direct O(N^5) algorithm
! bra and ket are the spin of (bra|ket)
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: bra
integer,intent(in) :: ket
integer,intent(in) :: nBas
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas)
double precision,intent(in) :: ERI_AO_basis(nBas,nBas,nBas,nBas),c(nBas,nBas,nspin)
! Local variables
@ -20,8 +23,11 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
double precision,intent(out) :: ERI_MO_basis(nBas,nBas,nBas,nBas)
! Memory allocation
allocate(scr(nBas,nBas,nBas,nBas))
! Four-index transform via semi-direct O(N^5) algorithm
scr(:,:,:,:) = 0d0
do l=1,nBas
@ -29,7 +35,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
do la=1,nBas
do nu=1,nBas
do mu=1,nBas
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l)
scr(mu,nu,la,l) = scr(mu,nu,la,l) + ERI_AO_basis(mu,nu,la,si)*c(si,l,ket)
enddo
enddo
enddo
@ -43,7 +49,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
do nu=1,nBas
do i=1,nBas
do mu=1,nBas
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i)*scr(mu,nu,la,l)
ERI_MO_basis(i,nu,la,l) = ERI_MO_basis(i,nu,la,l) + c(mu,i,bra)*scr(mu,nu,la,l)
enddo
enddo
enddo
@ -57,7 +63,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
do la=1,nBas
do nu=1,nBas
do i=1,nBas
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k)
scr(i,nu,k,l) = scr(i,nu,k,l) + ERI_MO_basis(i,nu,la,l)*c(la,k,bra)
enddo
enddo
enddo
@ -71,7 +77,7 @@ subroutine AOtoMO_integral_transform(nBas,c,ERI_AO_basis,ERI_MO_basis)
do j=1,nBas
do i=1,nBas
do nu=1,nBas
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j)*scr(i,nu,k,l)
ERI_MO_basis(i,j,k,l) = ERI_MO_basis(i,j,k,l) + c(nu,j,ket)*scr(i,nu,k,l)
enddo
! print*,i,k,j,l,ERI_MO_basis(i,j,k,l)
enddo

2
src/QuAcK/MP2.f90
View File

@ -1,6 +1,6 @@
subroutine MP2(nBas,nC,nO,nV,nR,ERI,ENuc,EHF,e,EcMP2)
! Perform third-order Moller-Plesset calculation
! Perform second-order Moller-Plesset calculation
implicit none

93
src/QuAcK/QuAcK.f90
View File

@ -4,6 +4,7 @@ program QuAcK
include 'parameters.h'
logical :: doSph
logical :: unrestricted
logical :: doRHF,doUHF,doMOM
logical :: doMP2,doMP3,doMP2F12
logical :: doCCD,doCCSD,doCCSDT
@ -27,8 +28,8 @@ program QuAcK
double precision,allocatable :: ZNuc(:),rNuc(:,:)
double precision,allocatable :: cHF(:,:,:),eHF(:,:),PHF(:,:,:)
double precision,allocatable :: eG0W0(:)
double precision,allocatable :: eG0T0(:)
double precision,allocatable :: eG0W0(:,:)
double precision,allocatable :: eG0T0(:,:)
logical :: doACFDT
logical :: exchange_kernel
@ -50,6 +51,11 @@ program QuAcK
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),H(:,:),X(:,:)
double precision,allocatable :: ERI_AO(:,:,:,:)
double precision,allocatable :: ERI_MO(:,:,:,:)
integer :: bra
integer :: ket
double precision,allocatable :: ERI_MO_aa(:,:,:,:)
double precision,allocatable :: ERI_MO_ab(:,:,:,:)
double precision,allocatable :: ERI_MO_bb(:,:,:,:)
double precision,allocatable :: ERI_ERF_AO(:,:,:,:)
double precision,allocatable :: ERI_ERF_MO(:,:,:,:)
double precision,allocatable :: F12(:,:,:,:),Yuk(:,:,:,:),FC(:,:,:,:,:,:)
@ -211,9 +217,8 @@ program QuAcK
! Memory allocation for one- and two-electron integrals
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas),eG0T0(nBas),PHF(nBas,nBas,nspin), &
S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas), &
ERI_AO(nBas,nBas,nBas,nBas),ERI_MO(nBas,nBas,nBas,nBas))
allocate(cHF(nBas,nBas,nspin),eHF(nBas,nspin),eG0W0(nBas,nspin),eG0T0(nBas,nspin),PHF(nBas,nBas,nspin), &
S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),H(nBas,nBas),X(nBas,nBas),ERI_AO(nBas,nBas,nBas,nBas))
! Read integrals
@ -263,6 +268,9 @@ program QuAcK
if(doUHF) then
! Switch on the unrestricted flag
unrestricted = .true.
call cpu_time(start_HF)
call UHF(maxSCF_HF,thresh_HF,n_diis_HF,guess_type,nBas,nO,S,T,V,Hc,ERI_AO,X,ENuc,EUHF,eHF,cHF,PHF)
call cpu_time(end_HF)
@ -307,13 +315,50 @@ program QuAcK
if(doSph) then
allocate(ERI_MO(nBas,nBas,nBas,nBas))
ERI_MO(:,:,:,:) = ERI_AO(:,:,:,:)
print*,'!!! MO = AO !!!'
deallocate(ERI_AO)
else
call AOtoMO_integral_transform(nBas,cHF,ERI_AO,ERI_MO)
if(unrestricted) then
! Memory allocation
allocate(ERI_MO_aa(nBas,nBas,nBas,nBas),ERI_MO_ab(nBas,nBas,nBas,nBas),ERI_MO_bb(nBas,nBas,nBas,nBas))
! 4-index transform for (aa|aa) block
bra = 1
ket = 1
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_aa)
! 4-index transform for (bb|bb) block
bra = 1
ket = 2
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_ab)
! 4-index transform for (aa|bb) block
bra = 2
ket = 2
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO_bb)
else
! Memory allocation
allocate(ERI_MO(nBas,nBas,nBas,nBas))
! 4-index transform
bra = 1
ket = 1
call AOtoMO_integral_transform(bra,ket,nBas,cHF,ERI_AO,ERI_MO)
end if
end if
@ -330,7 +375,17 @@ program QuAcK
if(doMP2) then
call cpu_time(start_MP2)
call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2)
if(unrestricted) then
call UMP2(nBas,nC,nO,nV,nR,ERI_MO_aa,ERI_MO_ab,ERI_MO_bb,ENuc,EUHF,eHF,EcMP2)
else
call MP2(nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF,EcMP2)
end if
call cpu_time(end_MP2)
t_MP2 = end_MP2 - start_MP2
@ -670,14 +725,24 @@ program QuAcK
! Perform G0W0 calculatiom
!------------------------------------------------------------------------
eG0W0(:) = eHF(:,1)
eG0W0(:,:) = eHF(:,:)
if(doG0W0) then
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_MO,PHF,cHF,eHF,eG0W0)
if(unrestricted) then
call UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0)
else
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_MO,PHF,cHF,eHF,eG0W0)
end if
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
@ -726,7 +791,7 @@ program QuAcK
! Perform G0T0 calculatiom
!------------------------------------------------------------------------
eG0T0(:) = eHF(:,1)
eG0T0(:,:) = eHF(:,:)
if(doG0T0) then
@ -868,7 +933,7 @@ program QuAcK
call cpu_time(start_G0W0)
call G0W0(doACFDT,exchange_kernel,doXBS,COHSEX,SOSEX,BSE,TDA_W,TDA, &
dBSE,dTDA,evDyn,singlet_manifold,triplet_manifold,linGW,eta_GW, &
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,H,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
nBas,nC,nO,nV,nR,nS,ENuc,ERHF,Hc,ERI_ERF_MO,PHF,cHF,eHF,eG0W0)
call cpu_time(end_G0W0)
t_G0W0 = end_G0W0 - start_G0W0
@ -889,7 +954,7 @@ program QuAcK
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for G0T0 = ',t_G0T0,' seconds'
write(*,*)
call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV)
! call matout(nBas,1,(eG0W0+eG0T0-eHF(:,1))*HaToeV)
end if

158
src/QuAcK/UMP2.f90 Normal file
View File

@ -0,0 +1,158 @@
subroutine UMP2(nBas,nC,nO,nV,nR,ERI_aa,ERI_ab,ERI_bb,ENuc,EHF,e,Ec)
! Perform unrestricted second-order Moller-Plesset calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas
integer,intent(in) :: nC(nspin)
integer,intent(in) :: nO(nspin)
integer,intent(in) :: nV(nspin)
integer,intent(in) :: nR(nspin)
double precision,intent(in) :: ENuc
double precision,intent(in) :: EHF
double precision,intent(in) :: ERI_aa(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_ab(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ERI_bb(nBas,nBas,nBas,nBas)
double precision,intent(in) :: e(nBas,nspin)
! Local variables
integer :: bra,ket
integer :: i,j,a,b
double precision :: eps
double precision :: Edaa,Exaa,Ecaa
double precision :: Edab,Exab,Ecab
double precision :: Edbb,Exbb,Ecbb
double precision :: Ed,Ex
! Output variables
double precision,intent(out) :: Ec
! Hello world
write(*,*)
write(*,*)'********************************************************'
write(*,*)'| Unrestricted second-order Moller-Plesset calculation |'
write(*,*)'********************************************************'
write(*,*)
!---------------------!
! Compute UMP2 energy |
!---------------------!
! aaaa block
bra = 1
ket = 1
Edaa = 0d0
Exaa = 0d0
do i=nC(bra)+1,nO(bra)
do a=nO(bra)+1,nBas-nR(bra)
do j=nC(ket)+1,nO(ket)
do b=nO(ket)+1,nBas-nR(ket)
eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket)
Edaa = Edaa + 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,a,b)/eps
Exaa = Exaa - 0.5d0*ERI_aa(i,j,a,b)*ERI_aa(i,j,b,a)/eps
enddo
enddo
enddo
enddo
Ecaa = Edaa + Exaa
! aabb block
bra = 1
ket = 2
Edab = 0d0
Exab = 0d0
do i=nC(bra)+1,nO(bra)
do a=nO(bra)+1,nBas-nR(bra)
do j=nC(ket)+1,nO(ket)
do b=nO(ket)+1,nBas-nR(ket)
eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket)
Edab = Edab + ERI_ab(i,j,a,b)*ERI_ab(i,j,a,b)/eps
enddo
enddo
enddo
enddo
Ecab = Edab + Exab
! bbbb block
bra = 2
ket = 2
Edbb = 0d0
Exbb = 0d0
do i=nC(bra)+1,nO(bra)
do a=nO(bra)+1,nBas-nR(bra)
do j=nC(ket)+1,nO(ket)
do b=nO(ket)+1,nBas-nR(ket)
eps = e(i,bra) + e(j,ket) - e(a,bra) - e(b,ket)
Edbb = Edbb + 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,a,b)/eps
Exbb = Exbb - 0.5d0*ERI_bb(i,j,a,b)*ERI_bb(i,j,b,a)/eps
enddo
enddo
enddo
enddo
Ecbb = Edbb + Exbb
! Final flush
Ed = Edaa + Edab + Edbb
Ex = Exaa + Exab + Exbb
Ec = Ed + Ex
write(*,*)
write(*,'(A32)') '--------------------------'
write(*,'(A32)') ' MP2 calculation '
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 correlation energy = ',Ec
write(*,'(A32,1X,F16.10)') ' alpha-alpha = ',Ecaa
write(*,'(A32,1X,F16.10)') ' alpha-beta = ',Ecab
write(*,'(A32,1X,F16.10)') ' beta-beta = ',Ecbb
write(*,*)
write(*,'(A32,1X,F16.10)') ' Direct part = ',Ed
write(*,'(A32,1X,F16.10)') ' alpha-alpha = ',Edaa
write(*,'(A32,1X,F16.10)') ' alpha-beta = ',Edab
write(*,'(A32,1X,F16.10)') ' beta-beta = ',Edbb
write(*,*)
write(*,'(A32,1X,F16.10)') ' Exchange part = ',Ex
write(*,'(A32,1X,F16.10)') ' alpha-alpha = ',Exaa
write(*,'(A32,1X,F16.10)') ' alpha-beta = ',Exab
write(*,'(A32,1X,F16.10)') ' beta-beta = ',Exbb
write(*,'(A32)') '--------------------------'
write(*,'(A32,1X,F16.10)') ' MP2 electronic energy = ', EHF + Ec
write(*,'(A32,1X,F16.10)') ' MP2 total energy = ',ENuc + EHF + Ec
write(*,'(A32)') '--------------------------'
write(*,*)
end subroutine UMP2

89
src/QuAcK/linear_response_A_matrix.f90
View File

@ -15,7 +15,7 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
! Local variables
double precision :: delta_spin,delta_dRPA
double precision :: delta_dRPA
double precision,external :: Kronecker_delta
integer :: i,j,a,b,ia,jb
@ -24,35 +24,78 @@ subroutine linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,
double precision,intent(out) :: A_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build A matrix
! Build A matrix for single manifold
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
if(ispin == 1) then
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ (1d0 + delta_spin)*lambda*ERI(i,b,a,j) &
- (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
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
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ 2d0*lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
enddo
enddo
enddo
enddo
end do
end do
end do
end do
end if
! Build A matrix for triplet manifold
if(ispin == 2) then
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
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
- (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
end do
end do
end do
end do
end if
! Build A matrix for spin orbitals
if(ispin == 3) then
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
A_lr(ia,jb) = (e(a) - e(i))*Kronecker_delta(i,j)*Kronecker_delta(a,b) &
+ lambda*ERI(i,b,a,j) - (1d0 - delta_dRPA)*lambda*ERI(i,b,j,a)
end do
end do
end do
end do
end if
end subroutine linear_response_A_matrix

85
src/QuAcK/linear_response_B_matrix.f90
View File

@ -14,7 +14,7 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
! Local variables
double precision :: delta_spin,delta_dRPA
double precision :: delta_dRPA
integer :: i,j,a,b,ia,jb
@ -22,34 +22,75 @@ subroutine linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B_
double precision,intent(out) :: B_lr(nS,nS)
! Singlet or triplet manifold?
delta_spin = 0d0
if(ispin == 1) delta_spin = +1d0
if(ispin == 2) delta_spin = -1d0
! Direct RPA
delta_dRPA = 0d0
if(dRPA) delta_dRPA = 1d0
! Build B matrix
! Build B matrix for singlet manifold
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
if(ispin == 1) then
B_lr(ia,jb) = (1d0 + delta_spin)*lambda*ERI(i,j,a,b) &
- (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
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
B_lr(ia,jb) = 2d0*lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
end do
end do
end do
end do
enddo
enddo
enddo
enddo
end if
! Build B matrix for triplet manifold
if(ispin == 2) then
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
B_lr(ia,jb) = - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
end do
end do
end do
end do
end if
! Build B matrix for spin orbitals
if(ispin == 3) then
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
B_lr(ia,jb) = lambda*ERI(i,j,a,b) - (1d0 - delta_dRPA)*lambda*ERI(i,j,b,a)
end do
end do
end do
end do
end if
end subroutine linear_response_B_matrix