LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-26 14:23:43 +01:00
Code Issues Releases Wiki Activity
f313e5b740
quantum_package/plugins/loos/NOFT/NOFT_JKfunc.irp.f
Anthony Scemama 5b6067ea2f Cleaning plugins
2018-10-19 11:32:58 +02:00

261 lines
6.1 KiB
Fortran
Raw Blame History

subroutine NOFT_JKfunc(nMO,FL,ET,EV,n)
implicit none
BEGIN_DOC
! JK-only functionals for NOFT
END_DOC
PROVIDE mo_bielec_integrals_in_map
! Input variables
integer,intent(in) :: nMO,FL
double precision,intent(in) :: ET,EV
double precision,intent(in) :: n(nMO)
! Local variables
integer :: i,j
double precision :: EJ_SD
double precision :: EK_SD,EK_MBB,EK_POWER,EK_BBC1,EK_BBC2,EK_CA,EK_CGA,EK_GU,EK_ML,EK_MLSIC
double precision :: E_SD,E_MBB,E_POWER,E_BBC1,E_BBC2,E_CA,E_CGA,E_GU,E_ML,E_MLSIC
double precision :: alpha,a0,a1,b1
double precision :: f_ij,get_mo_bielec_integral
double precision,allocatable :: Jint(:,:),Kint(:,:)
! Memory allocation
allocate(Jint(nMO,nMO),Kint(nMO,nMO))
! Coulomb, exchange and time-inversion integrals
do i=1,nMO
do j=1,nMO
Jint(i,j) = get_mo_bielec_integral(i,j,i,j,mo_integrals_map)
Kint(i,j) = get_mo_bielec_integral(i,j,j,i,mo_integrals_map)
enddo
enddo
! Compute SD Coulomb energy
EJ_SD = 2d0*dot_product(n,matmul(Jint,n))
! Compute SD exchange energy
EK_SD = dot_product(n,matmul(Kint,n))
! Compute MBB exchange energy
EK_MBB = 0d0
do i=1,nMO
do j=1,nMO
EK_MBB = EK_MBB + sqrt(n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute BBC1 exchange energy
EK_BBC1 = 0d0
do i=1,nMO
do j=1,nMO
if(i /= j .and. i > FL .and. j > FL) then
f_ij = - sqrt(n(i)*n(j))
else
f_ij = sqrt(n(i)*n(j))
endif
EK_BBC1 = EK_BBC1 + f_ij*Kint(i,j)
enddo
enddo
! Compute BBC2 exchange energy
EK_BBC2 = 0d0
do i=1,nMO
do j=1,i-1
if(i > FL .and. j > FL) then
f_ij = - sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
f_ij = n(i)*n(j)
else
f_ij = sqrt(n(i)*n(j))
endif
EK_BBC2 = EK_BBC2 + f_ij*Kint(i,j)
enddo
EK_BBC2 = EK_BBC2 + n(i)*Kint(i,i)
do j=i+1,nMO
if(i > FL .and. j > FL) then
f_ij = - sqrt(n(i)*n(j))
elseif(i <= FL .and. j <= FL) then
f_ij = n(i)*n(j)
else
f_ij = sqrt(n(i)*n(j))
endif
EK_BBC2 = EK_BBC2 + f_ij*Kint(i,j)
enddo
enddo
! Compute CA exchange energy
EK_CA = 0d0
do i=1,nMO
do j=1,nMO
EK_CA = EK_CA + (sqrt(n(i)*n(j)*(1d0 - n(i))*(1d0 - n(j))) + n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute CGA exchange energy
EK_CGA = 0d0
do i=1,nMO
do j=1,nMO
EK_CGA = EK_CGA + 0.5d0*(sqrt(n(i)*n(j)*(2d0 - n(i))*(2d0 - n(j))) + n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute ML exchange energy
EK_ML = 0d0
a0 = 126.3101d0
a1 = 2213.33d0
b1 = 2338.64d0
do i=1,nMO
do j=1,nMO
EK_ML = EK_ML + n(i)*n(j)*(a0 + a1*n(i)*n(j))/(1d0 + b1*n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute MLSIC exchange energy
EK_MLSIC = 0d0
a0 = 1298.78d0
a1 = 35114.4d0
b1 = 36412.2d0
do i=1,nMO
do j=1,i-1
EK_MLSIC = EK_MLSIC + n(i)*n(j)*(a0 + a1*n(i)*n(j))/(1d0 + b1*n(i)*n(j))*Kint(i,j)
enddo
EK_MLSIC = EK_MLSIC + n(i)*n(i)*Kint(i,i)
do j=i+1,nMO
EK_MLSIC = EK_MLSIC + n(i)*n(j)*(a0 + a1*n(i)*n(j))/(1d0 + b1*n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute GU exchange energy
EK_GU = 0d0
do i=1,nMO
do j=1,i-1
EK_GU = EK_GU + sqrt(n(i)*n(j))*Kint(i,j)
enddo
EK_GU = EK_GU + n(i)*n(i)*Kint(i,j)
do j=i+1,nMO
EK_GU = EK_GU + sqrt(n(i)*n(j))*Kint(i,j)
enddo
enddo
! Compute POWER exchange energy
EK_POWER = 0d0
alpha = 1d0/3d0
do i=1,nMO
do j=1,nMO
EK_POWER = EK_POWER + (n(i)*n(j))**alpha*Kint(i,j)
enddo
enddo
! Compute total energies
E_SD = ET + EV + EJ_SD - EK_SD
E_MBB = ET + EV + EJ_SD - EK_MBB
E_BBC1 = ET + EV + EJ_SD - EK_BBC1
E_BBC2 = ET + EV + EJ_SD - EK_BBC2
E_CA = ET + EV + EJ_SD - EK_CA
E_CGA = ET + EV + EJ_SD - EK_CGA
E_ML = ET + EV + EJ_SD - EK_ML
E_MLSIC = ET + EV + EJ_SD - EK_MLSIC
E_GU = ET + EV + EJ_SD - EK_GU
E_POWER = ET + EV + EJ_SD - EK_POWER
! Dump energies
print*, '*******************************'
print*, '*** JK NOFT functionals ***'
print*, '*******************************'
print*, ''
print*, '*** Coulomb energies ***'
print*, 'Coulomb SD energy = ',EJ_SD
print*, ''
print*, '*** Exchange energies ***'
print*, 'Exchange SD energy = ',-EK_SD
print*, 'Exchange MBB energy = ',-EK_MBB
print*, 'Exchange BBC1 energy = ',-EK_BBC1
print*, 'Exchange BBC2 energy = ',-EK_BBC2
print*, 'Exchange CA energy = ',-EK_CA
print*, 'Exchange CGA energy = ',-EK_CGA
print*, 'Exchange ML energy = ',-EK_ML
print*, 'Exchange MLSIC energy = ',-EK_MLSIC
print*, 'Exchange GU energy = ',-EK_GU
print*, 'Exchange POWER energy = ',-EK_POWER
print*, ''
print*, ''
print*, '*** Two-electron energies ***'
print*, 'J+K SD energy = ',EJ_SD - EK_SD
print*, 'J+K MBB energy = ',EJ_SD - EK_MBB
print*, 'J+K BBC1 energy = ',EJ_SD - EK_BBC1
print*, 'J+K BBC2 energy = ',EJ_SD - EK_BBC2
print*, 'J+K CA energy = ',EJ_SD - EK_CA
print*, 'J+K CGA energy = ',EJ_SD - EK_CGA
print*, 'J+K ML energy = ',EJ_SD - EK_ML
print*, 'J+K MLSIC energy = ',EJ_SD - EK_MLSIC
print*, 'J+K GU energy = ',EJ_SD - EK_GU
print*, 'J+K POWER energy = ',EJ_SD - EK_POWER
print*, ''
print*, '*** Total energies ***'
print*, 'Total SD energy = ',E_SD
print*, 'Total MBB energy = ',E_MBB
print*, 'Total BBC1 energy = ',E_BBC1
print*, 'Total BBC2 energy = ',E_BBC2
print*, 'Total CA energy = ',E_CA
print*, 'Total CGA energy = ',E_CGA
print*, 'Total ML energy = ',E_ML
print*, 'Total MLSIC energy = ',E_MLSIC
print*, 'Total GU energy = ',E_GU
print*, 'Total POWER energy = ',E_POWER
print*, ''
end subroutine NOFT_JKfunc
Reference in New Issue View Git Blame Copy Permalink