mirror of
https://github.com/pfloos/quack
synced 2025-01-03 10:05:49 +01:00
new files
This commit is contained in:
parent
847160ad09
commit
22d83c824f
3
GoDuck
3
GoDuck
@ -11,8 +11,9 @@ if [ $# = 2 ]
|
|||||||
then
|
then
|
||||||
cp examples/molecule."$1" input/molecule
|
cp examples/molecule."$1" input/molecule
|
||||||
cp examples/basis."$1"."$2" input/basis
|
cp examples/basis."$1"."$2" input/basis
|
||||||
|
cp basis/"$2" input/basis.qcaml
|
||||||
cp examples/basis."$1"."$2" input/weight
|
cp examples/basis."$1"."$2" input/weight
|
||||||
./bin/IntPak
|
# ./bin/IntPak
|
||||||
./bin/QuAcK
|
./bin/QuAcK
|
||||||
fi
|
fi
|
||||||
|
|
||||||
|
|||||||
4
GoQCaml
Executable file
4
GoQCaml
Executable file
@ -0,0 +1,4 @@
|
|||||||
|
#! /bin/bash
|
||||||
|
|
||||||
|
cd int
|
||||||
|
../utils/QCaml/run_integrals -b ../input/basis.qcaml -x ../input/molecule.xyz
|
||||||
@ -2,4 +2,4 @@
|
|||||||
2 7 7 0 0
|
2 7 7 0 0
|
||||||
# Znuc x y z
|
# Znuc x y z
|
||||||
N 0. 0. 0.
|
N 0. 0. 0.
|
||||||
N 0. 0. 3.4
|
N 0. 0. 2.0
|
||||||
|
|||||||
@ -2,4 +2,4 @@
|
|||||||
2 7 7 0 0
|
2 7 7 0 0
|
||||||
# Znuc x y z
|
# Znuc x y z
|
||||||
N 0. 0. 0.
|
N 0. 0. 0.
|
||||||
N 0. 0. 3.4
|
N 0. 0. 2.0
|
||||||
|
|||||||
148
src/QuAcK/ACFDT.f90
Normal file
148
src/QuAcK/ACFDT.f90
Normal file
@ -0,0 +1,148 @@
|
|||||||
|
subroutine ACFDT(exchange_kernel,doXBS,dRPA,TDA,BSE,singlet_manifold,triplet_manifold, &
|
||||||
|
nBas,nC,nO,nV,nR,nS,ERI,e,Omega,XpY,XmY,rho,EcAC)
|
||||||
|
|
||||||
|
! Compute the correlation energy via the adiabatic connection fluctuation dissipation theorem
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
include 'quadrature.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
logical,intent(in) :: doXBS
|
||||||
|
logical,intent(in) :: exchange_kernel
|
||||||
|
logical,intent(in) :: dRPA
|
||||||
|
logical,intent(in) :: TDA
|
||||||
|
logical,intent(in) :: BSE
|
||||||
|
logical,intent(in) :: singlet_manifold
|
||||||
|
logical,intent(in) :: triplet_manifold
|
||||||
|
|
||||||
|
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
|
||||||
|
double precision,intent(in) :: e(nBas)
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
|
||||||
|
double precision :: Omega(nS,nspin)
|
||||||
|
double precision :: XpY(nS,nS,nspin)
|
||||||
|
double precision :: XmY(nS,nS,nspin)
|
||||||
|
double precision :: rho(nBas,nBas,nS,nspin)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: ispin
|
||||||
|
integer :: iAC
|
||||||
|
double precision :: lambda
|
||||||
|
double precision,allocatable :: Ec(:,:)
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: EcAC(nspin)
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(Ec(nAC,nspin))
|
||||||
|
|
||||||
|
! Antisymmetrized kernel version
|
||||||
|
|
||||||
|
if(exchange_kernel) then
|
||||||
|
|
||||||
|
write(*,*)
|
||||||
|
write(*,*) '*** Exchange kernel version ***'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Singlet manifold
|
||||||
|
|
||||||
|
if(singlet_manifold) then
|
||||||
|
|
||||||
|
ispin = 1
|
||||||
|
EcAC(ispin) = 0d0
|
||||||
|
Ec(:,ispin) = 0d0
|
||||||
|
|
||||||
|
write(*,*) '--------------'
|
||||||
|
write(*,*) 'Singlet states'
|
||||||
|
write(*,*) '--------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
|
||||||
|
do iAC=1,nAC
|
||||||
|
|
||||||
|
lambda = rAC(iAC)
|
||||||
|
|
||||||
|
if(doXBS) then
|
||||||
|
|
||||||
|
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
|
||||||
|
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
|
||||||
|
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
|
||||||
|
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
|
||||||
|
|
||||||
|
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
|
||||||
|
|
||||||
|
end do
|
||||||
|
|
||||||
|
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
|
||||||
|
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
! Triplet manifold
|
||||||
|
|
||||||
|
if(triplet_manifold) then
|
||||||
|
|
||||||
|
ispin = 2
|
||||||
|
EcAC(ispin) = 0d0
|
||||||
|
Ec(:,ispin) = 0d0
|
||||||
|
|
||||||
|
write(*,*) '--------------'
|
||||||
|
write(*,*) 'Triplet states'
|
||||||
|
write(*,*) '--------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A15,1X,A30,1X,A30)') 'lambda','Ec(lambda)','Tr(K x P_lambda)'
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
|
||||||
|
do iAC=1,nAC
|
||||||
|
|
||||||
|
lambda = rAC(iAC)
|
||||||
|
|
||||||
|
if(doXBS) then
|
||||||
|
|
||||||
|
call linear_response(ispin,dRPA,TDA,.false.,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
|
||||||
|
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
call excitation_density(nBas,nC,nO,nR,nS,ERI,XpY(:,:,ispin),rho(:,:,:,ispin))
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
call linear_response(ispin,dRPA,TDA,BSE,nBas,nC,nO,nV,nR,nS,lambda,e,ERI, &
|
||||||
|
rho(:,:,:,ispin),EcAC(ispin),Omega(:,ispin),XpY(:,:,ispin),XmY(:,:,ispin))
|
||||||
|
|
||||||
|
call ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY(:,:,ispin),XmY(:,:,ispin),Ec(iAC,ispin))
|
||||||
|
|
||||||
|
write(*,'(2X,F15.6,1X,F30.15,1X,F30.15)') lambda,EcAC(ispin),Ec(iAC,ispin)
|
||||||
|
|
||||||
|
end do
|
||||||
|
|
||||||
|
EcAC(ispin) = 0.5d0*dot_product(wAC,Ec(:,ispin))
|
||||||
|
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
write(*,'(2X,A50,1X,F15.6)') ' Ec(AC) via Gauss-Legendre quadrature:',EcAC(ispin)
|
||||||
|
write(*,*) '-----------------------------------------------------------------------------------'
|
||||||
|
write(*,*)
|
||||||
|
|
||||||
|
end if
|
||||||
|
|
||||||
|
end subroutine ACFDT
|
||||||
93
src/QuAcK/ACFDT_correlation_energy.f90
Normal file
93
src/QuAcK/ACFDT_correlation_energy.f90
Normal file
@ -0,0 +1,93 @@
|
|||||||
|
subroutine ACFDT_correlation_energy(ispin,exchange_kernel,nBas,nC,nO,nV,nR,nS,ERI,XpY,XmY,EcAC)
|
||||||
|
|
||||||
|
! Compute the correlation energy via the adiabatic connection formula
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
include 'parameters.h'
|
||||||
|
|
||||||
|
! Input variables
|
||||||
|
|
||||||
|
integer,intent(in) :: ispin
|
||||||
|
logical,intent(in) :: exchange_kernel
|
||||||
|
integer,intent(in) :: nBas,nC,nO,nV,nR,nS
|
||||||
|
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
|
||||||
|
double precision,intent(in) :: XpY(nS,nS)
|
||||||
|
double precision,intent(in) :: XmY(nS,nS)
|
||||||
|
|
||||||
|
! Local variables
|
||||||
|
|
||||||
|
integer :: i,j,a,b
|
||||||
|
integer :: ia,jb,kc
|
||||||
|
double precision :: delta_spin
|
||||||
|
double precision :: delta_Kx
|
||||||
|
double precision,allocatable :: P(:,:)
|
||||||
|
double precision,allocatable :: Ap(:,:)
|
||||||
|
double precision,allocatable :: Bp(:,:)
|
||||||
|
double precision,allocatable :: X(:,:)
|
||||||
|
double precision,allocatable :: Y(:,:)
|
||||||
|
double precision,external :: trace_matrix
|
||||||
|
|
||||||
|
! Output variables
|
||||||
|
|
||||||
|
double precision,intent(out) :: EcAC
|
||||||
|
|
||||||
|
! Singlet or triplet manifold?
|
||||||
|
|
||||||
|
delta_spin = 0d0
|
||||||
|
if(ispin == 1) delta_spin = +1d0
|
||||||
|
if(ispin == 2) delta_spin = -1d0
|
||||||
|
|
||||||
|
! Exchange kernel
|
||||||
|
|
||||||
|
delta_Kx = 0d0
|
||||||
|
if(exchange_kernel) delta_Kx = 1d0
|
||||||
|
|
||||||
|
! Memory allocation
|
||||||
|
|
||||||
|
allocate(P(nS,nS),Ap(nS,nS),Bp(nS,nS),X(nS,nS),Y(nS,nS))
|
||||||
|
|
||||||
|
! Compute P = (X+Y)^T.(X+Y) - 1
|
||||||
|
|
||||||
|
P(:,:) = matmul(transpose(XpY),XpY)
|
||||||
|
|
||||||
|
do ia=1,nS
|
||||||
|
P(ia,ia) = P(ia,ia) - 1d0
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! Compute Aiajb = (ia|bj) and Biajb = (ia|jb)
|
||||||
|
|
||||||
|
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
|
||||||
|
|
||||||
|
Ap(ia,jb) = (1d0 + delta_spin)*ERI(i,b,a,j) &
|
||||||
|
- delta_Kx*ERI(i,b,j,a)
|
||||||
|
|
||||||
|
Bp(ia,jb) = (1d0 + delta_spin)*ERI(i,j,a,b) &
|
||||||
|
- delta_Kx*ERI(i,j,b,a)
|
||||||
|
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
! Compute Tr(K x P_lambda)
|
||||||
|
|
||||||
|
! EcAC = trace_matrix(nS,matmul(Ap,P))
|
||||||
|
|
||||||
|
X(:,:) = 0.5d0*(XpY(:,:) + XmY(:,:))
|
||||||
|
Y(:,:) = 0.5d0*(XpY(:,:) - XmY(:,:))
|
||||||
|
|
||||||
|
EcAC = trace_matrix(nS,matmul(X,matmul(Bp,transpose(Y))) + matmul(Y,matmul(Bp,transpose(X)))) &
|
||||||
|
+ trace_matrix(nS,matmul(X,matmul(Ap,transpose(X))) + matmul(Y,matmul(Ap,transpose(Y)))) &
|
||||||
|
- trace_matrix(nS,Ap)
|
||||||
|
|
||||||
|
! call matout(nS,nS,matmul(transpose(X),X) - matmul(transpose(Y),Y))
|
||||||
|
|
||||||
|
end subroutine ACFDT_correlation_energy
|
||||||
|
|
||||||
@ -205,9 +205,8 @@ program QuAcK
|
|||||||
|
|
||||||
else
|
else
|
||||||
|
|
||||||
! call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
|
call system('./GoQCaml')
|
||||||
call
|
call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
|
||||||
! call read_integrals(nEl(:),nBas,S,T,V,Hc,ERI_AO_basis)
|
|
||||||
|
|
||||||
end if
|
end if
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user