LCPQ/quack
4
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-06-18 11:15:30 +02:00
Code Releases Activity

remove xcDFT

Browse Source
This commit is contained in:
Pierre-Francois Loos 2019-02-07 22:54:43 +01:00
parent 86964672d7
commit 905a4821cd
67 changed files with 0 additions and 5418 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

101
src/xcDFT/AO_values_grid.f90
View File

@ -1,101 +0,0 @@
subroutine AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
nGrid,root,AO,dAO)
! Compute values of the AOs and their derivatives with respect to the cartesian coordinates
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nBas,nShell
double precision,intent(in) :: CenterShell(maxShell,3)
integer,intent(in) :: TotAngMomShell(maxShell)
integer,intent(in) :: KShell(maxShell)
double precision,intent(in) :: DShell(maxShell,maxK)
double precision,intent(in) :: ExpShell(maxShell,maxK)
double precision,intent(in) :: root(3,nGrid)
integer,intent(in) :: nGrid
! Local variables
integer :: atot,nShellFunction,a(3)
integer,allocatable :: ShellFunction(:,:)
double precision :: rASq,xA,yA,zA,NormCoeff,prim
integer :: iSh,iShF,iK,iG,iBas
! Output variables
double precision,intent(out) :: AO(nBas,nGrid)
double precision,intent(out) :: dAO(3,nBas,nGrid)
! Initialization
iBas = 0
AO(:,:) = 0d0
dAO(:,:,:) = 0d0
!------------------------------------------------------------------------
! Loops over shells
!------------------------------------------------------------------------
do iSh=1,nShell
atot = TotAngMomShell(iSh)
nShellFunction = (atot*atot + 3*atot + 2)/2
allocate(ShellFunction(1:nShellFunction,1:3))
call generate_shell(atot,nShellFunction,ShellFunction)
do iShF=1,nShellFunction
iBas = iBas + 1
a(:) = ShellFunction(iShF,:)
do iG=1,nGrid
xA = root(1,iG) - CenterShell(iSh,1)
yA = root(2,iG) - CenterShell(iSh,2)
zA = root(3,iG) - CenterShell(iSh,3)
! Calculate distance for exponential
rASq = xA**2 + yA**2 + zA**2
!------------------------------------------------------------------------
! Loops over contraction degrees
!-------------------------------------------------------------------------
do iK=1,KShell(iSh)
! Calculate the exponential part
prim = DShell(iSh,iK)*NormCoeff(ExpShell(iSh,iK),a)*exp(-ExpShell(iSh,iK)*rASq)
AO(iBas,iG) = AO(iBas,iG) + prim
prim = -2d0*ExpShell(iSh,iK)*prim
dAO(:,iBas,iG) = dAO(:,iBas,iG) + prim
enddo
dAO(1,iBas,iG) = xA**(a(1)+1)*yA**a(2)*zA**a(3)*dAO(1,iBas,iG)
if(a(1) > 0) dAO(1,iBas,iG) = dAO(1,iBas,iG) + dble(a(1))*xA**(a(1)-1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
dAO(2,iBas,iG) = xA**a(1)*yA**(a(2)+1)*zA**a(3)*dAO(2,iBas,iG)
if(a(2) > 0) dAO(2,iBas,iG) = dAO(2,iBas,iG) + dble(a(2))*xA**a(1)*yA**(a(2)-1)*zA**a(3)*AO(iBas,iG)
dAO(3,iBas,iG) = xA**a(1)*yA**a(2)*zA**(a(3)+1)*dAO(3,iBas,iG)
if(a(3) > 0) dAO(3,iBas,iG) = dAO(3,iBas,iG) + dble(a(3))*xA**a(1)*yA**a(2)*zA**(a(3)-1)*AO(iBas,iG)
! Calculate polynmial part
AO(iBas,iG) = xA**a(1)*yA**a(2)*zA**a(3)*AO(iBas,iG)
enddo
enddo
deallocate(ShellFunction)
enddo
!------------------------------------------------------------------------
! End loops over shells
!------------------------------------------------------------------------
end subroutine AO_values_grid

52
src/xcDFT/DIIS_extrapolation.f90
View File

@ -1,52 +0,0 @@
subroutine DIIS_extrapolation(n,n_diis,error,e,error_in,e_inout)
! Perform DIIS extrapolation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: n,n_diis
double precision,intent(in) :: error(n,n_diis),e(n,n_diis),error_in(n)
! Local variables
double precision,allocatable :: A(:,:),b(:),w(:)
! Output variables
double precision,intent(inout):: e_inout(n)
! Memory allocaiton
allocate(A(n_diis+1,n_diis+1),b(n_diis+1),w(n_diis+1))
! Update DIIS "history"
call prepend(n,n_diis,error,error_in)
call prepend(n,n_diis,e,e_inout)
! Build A matrix
A(1:n_diis,1:n_diis) = matmul(transpose(error),error)
A(1:n_diis,n_diis+1) = -1d0
A(n_diis+1,1:n_diis) = -1d0
A(n_diis+1,n_diis+1) = +0d0
! Build x matrix
b(1:n_diis) = +0d0
b(n_diis+1) = -1d0
! Solve linear system
call linear_solve(n_diis+1,A,b,w)
! Extrapolate
e_inout(:) = matmul(w(1:n_diis),transpose(e(:,1:n_diis)))
end subroutine DIIS_extrapolation

34
src/xcDFT/Makefile
View File

@ -1,34 +0,0 @@
IDIR =../../include
BDIR =../../bin
ODIR = obj
SDIR =.
FC = gfortran -I$(IDIR)
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 -O2
endif
LIBS = ~/Dropbox/quack/lib/*.a
#LIBS = -lblas -llapack
SRCF90 = $(wildcard *.f90)
SRC = $(wildcard *.f)
OBJ = $(patsubst %.f90,$(ODIR)/%.o,$(SRCF90)) $(patsubst %.f,$(ODIR)/%.o,$(SRC))
$(ODIR)/%.o: %.f90
$(FC) -c -o $@ $< $(FFLAGS)
$(ODIR)/%.o: %.f
$(FC) -c -o $@ $< $(FFLAGS)
$(BDIR)/xcDFT: $(OBJ)
$(FC) -o $@ $^ $(FFLAGS) $(LIBS)
debug:
DEBUG=1 make clean $(BDIR)/xcDFT
clean:
rm -f $(ODIR)/*.o $(BDIR)/xcDFT $(BDIR)/debug

31
src/xcDFT/NormCoeff.f90
View File

@ -1,31 +0,0 @@
function NormCoeff(alpha,a)
! Compute normalization coefficients for cartesian gaussians
implicit none
! Input variables
double precision,intent(in) :: alpha
integer,intent(in) :: a(3)
! local variable
double precision :: pi,dfa(3),dfac
integer :: atot
! Output variable
double precision NormCoeff
pi = 4d0*atan(1d0)
atot = a(1) + a(2) + a(3)
dfa(1) = dfac(2*a(1))/(2d0**a(1)*dfac(a(1)))
dfa(2) = dfac(2*a(2))/(2d0**a(2)*dfac(a(2)))
dfa(3) = dfac(2*a(3))/(2d0**a(3)*dfac(a(3)))
NormCoeff = (2d0*alpha/pi)**(3d0/2d0)*(4d0*alpha)**atot
NormCoeff = NormCoeff/(dfa(1)*dfa(2)*dfa(3))
NormCoeff = sqrt(NormCoeff)
end function NormCoeff

221
src/xcDFT/RKS.f90
View File

@ -1,221 +0,0 @@
subroutine RKS(rung,nGrid,weight,nBas,AO,dAO,nO,S,T,V,Hc,ERI,X,ENuc,EKS)
! Perform a restricted Kohn-Sham calculation
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
integer,intent(in) :: nO
double precision,intent(in) :: S(nBas,nBas)
double precision,intent(in) :: T(nBas,nBas)
double precision,intent(in) :: V(nBas,nBas)
double precision,intent(in) :: Hc(nBas,nBas)
double precision,intent(in) :: X(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: ENuc
! Local variables
integer,parameter :: maxSCF = 64
double precision,parameter :: thresh = 1d-5
integer,parameter :: n_diis = 1
integer :: nSCF
double precision :: Conv
double precision :: ET,EV,EJ
double precision :: Ex
double precision :: Ec
double precision,allocatable :: e(:)
double precision,allocatable :: c(:,:),cp(:,:)
double precision,allocatable :: P(:,:),Pa(:,:)
double precision,allocatable :: J(:,:)
double precision,allocatable :: F(:,:),Fp(:,:)
double precision,allocatable :: Fx(:,:),FxHF(:,:)
double precision,allocatable :: Fc(:,:)
double precision,allocatable :: error(:,:)
double precision,allocatable :: error_diis(:,:),F_diis(:,:)
double precision,external :: trace_matrix
double precision,external :: exchange_energy
double precision,external :: electron_number
double precision,allocatable :: rhoa(:)
double precision,allocatable :: drhoa(:,:)
double precision :: nEl
! Output variables
double precision,intent(out) :: EKS
! Hello world
write(*,*)
write(*,*)'************************************************'
write(*,*)'| Restricted Kohn-Sham calculation |'
write(*,*)'************************************************'
write(*,*)
!------------------------------------------------------------------------
! Rung of Jacob's ladder
!------------------------------------------------------------------------
call select_rung(rung)
! Memory allocation
allocate(e(nBas),c(nBas,nBas),cp(nBas,nBas),P(nBas,nBas),Pa(nBas,nBas), &
J(nBas,nBas),F(nBas,nBas),Fp(nBas,nBas),Fx(nBas,nBas),FxHF(nBas,nBas), &
Fc(nBas,nBas),error(nBas,nBas),rhoa(nGrid),drhoa(3,nGrid), &
error_diis(nBas*nBas,n_diis),F_diis(nBas*nBas,n_diis))
! Guess coefficients and eigenvalues
F(:,:) = Hc(:,:)
! Initialization
nSCF = 0
Conv = 1d0
nEl = 0d0
Ex = 0d0
Ec = 0d0
Fx(:,:) = 0d0
FxHF(:,:) = 0d0
Fc(:,:) = 0d0
F_diis(:,:) = 0d0
error_diis(:,:) = 0d0
!------------------------------------------------------------------------
! Main SCF loop
!------------------------------------------------------------------------
write(*,*)
write(*,*)'------------------------------------------------------------------------------------------'
write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') &
'|','#','|','EKS','|','ExKS','|','EcKS','|','Conv','|','nEl','|'
write(*,*)'------------------------------------------------------------------------------------------'
do while(Conv > thresh .and. nSCF < maxSCF)
! Increment
nSCF = nSCF + 1
! Transform Fock matrix in orthogonal basis
Fp = matmul(transpose(X),matmul(F,X))
! Diagonalize Fock matrix to get eigenvectors and eigenvalues
cp(:,:) = Fp(:,:)
call diagonalize_matrix(nBas,cp,e)
! Back-transform eigenvectors in non-orthogonal basis
c = matmul(X,cp)
! Compute density matrix
Pa(:,:) = matmul(c(:,1:nO),transpose(c(:,1:nO)))
P(:,:) = 2d0*Pa(:,:)
! Compute one-electron density and its gradient if necessary
call density(nGrid,nBas,Pa,AO,rhoa)
if(rung > 1) call gradient_density(nGrid,nBas,Pa,AO,dAO,drhoa)
! Build Coulomb repulsion
call hartree_coulomb(nBas,P,ERI,J)
! Compute exchange potential
call exchange_potential(rung,nGrid,weight,nBas,Pa,ERI,AO,dAO,rhoa,drhoa,Fx,FxHF)
! Compute correlation potential
! call correlation_potential(rung,nGrid,weight,nBas,Pa,ERI,AO,dAO,rhoa,drhoa,Fc)
! Build Fock operator
F(:,:) = Hc(:,:) + J(:,:) + Fx(:,:) + Fc(:,:)
! Check convergence
error = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F)
Conv = maxval(abs(error))
! DIIS extrapolation
call DIIS_extrapolation(nBas*nBas,min(n_diis,nSCF),error_diis,F_diis,error,F)
!------------------------------------------------------------------------
! Compute KS energy
!------------------------------------------------------------------------
! Kinetic energy
ET = trace_matrix(nBas,matmul(P,T))
! Potential energy
EV = trace_matrix(nBas,matmul(P,V))
! Coulomb energy
EJ = 0.5d0*trace_matrix(nBas,matmul(P,J))
! Exchange energy
Ex = exchange_energy(rung,nGrid,weight,nBas,Pa,FxHF,rhoa,drhoa)
! Correlation energy
! call correlation_energy(rung,nGrid,weight,nBas,Pa,rhoa,drhoa,Ec)
EKS = ET + EV + EJ + Ex + Ec
! Check the grid accuracy by computing the number of electrons
nEl = electron_number(nGrid,weight,rhoa)
! Dump results
write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') &
'|',nSCF,'|',EKS+ENuc,'|',Ex,'|',Ec,'|',Conv,'|',nEl,'|'
enddo
write(*,*)'------------------------------------------------------------------------------------------'
!------------------------------------------------------------------------
! End of SCF loop
!------------------------------------------------------------------------
! Did it actually converge?
if(nSCF == maxSCF) then
write(*,*)
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)' Convergence failed '
write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!'
write(*,*)
stop
endif
! Compute final KS energy
call print_RKS(nBas,nO,e,C,ENuc,ET,EV,EJ,Ex,Ec,EKS)
end subroutine RKS

38
src/xcDFT/density.f90
View File

@ -1,38 +0,0 @@
subroutine density(nGrid,nBas,P,AO,rho)
! Calculate one-electron density
implicit none
include 'parameters.h'
! Input variables
double precision,parameter :: thresh = 1d-15
integer,intent(in) :: nGrid
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
! Local variables
integer :: iG,mu,nu
! Output variables
double precision,intent(out) :: rho(nGrid)
rho(:) = 0d0
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
enddo
enddo
enddo
! do iG=1,nGrid
! rho(iG) = max(rho(iG),thresh)
! enddo
end subroutine density

3017
src/xcDFT/dft_grid.f
View File

File diff suppressed because it is too large Load Diff

20
src/xcDFT/electron_number.f90
View File

@ -1,20 +0,0 @@
function electron_number(nGrid,w,rho) result(nEl)
! Compute the number of electrons via integration of the one-electron density
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: w(nGrid)
double precision,intent(in) :: rho(nGrid)
! Output variables
double precision :: nEl
nEl = 2d0*dot_product(w,rho)
end function electron_number

171
src/xcDFT/elements.f90
View File

@ -1,171 +0,0 @@
function element_number(element_name)
implicit none
integer,parameter :: nelement_max = 103
character(len=2),intent(in) :: element_name
integer :: element_number
character(len=2),parameter :: element_list(nelement_max) = &
(/' H', 'He', & ! 2
'Li','Be', ' B',' C',' N',' O',' F','Ne', & ! 10
'Na','Mg', 'Al','Si',' P',' S','Cl','Ar', & ! 18
' K','Ca','Sc','Ti',' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se','Br','Kr', & ! 36
'Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd','Ag','Cd','In','Sn','Sb','Te',' I','Xe', & ! 54
'Cs','Ba', & ! 56
'La','Ce','Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb', & ! 70
'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb','Bi','Po','At','Rn', & ! 86
'Fr','Ra', & ! 88
'Ac','Th','Pa',' U','Np','Pu','Am','Cm','Bk','Cf','Es','Fm','Md','No', & ! 102
'Lr' & ! 103
/)
!=====
integer :: ielement
!=====
ielement=1
do while( ADJUSTL(element_name) /= ADJUSTL(element_list(ielement)) )
if( ielement == nelement_max ) then
write(*,'(a,a)') ' Input symbol ',element_name
write(*,'(a,i3,a)') ' Element symbol is not one of first ',nelement_max,' elements'
write(*,*) '!!! element symbol not understood !!!'
stop
endif
ielement = ielement + 1
enddo
element_number = ielement
end function element_number
function element_core(zval,zatom)
implicit none
double precision,intent(in) :: zval
double precision,intent(in) :: zatom
integer :: element_core
!=====
!
! If zval /= zatom, this is certainly an effective core potential
! and no core states should be frozen.
if( ABS(zval - zatom) > 1d0-3 ) then
element_core = 0
else
if( zval <= 4.00001d0 ) then ! up to Be
element_core = 0
else if( zval <= 12.00001d0 ) then ! up to Mg
element_core = 1
else if( zval <= 30.00001d0 ) then ! up to Ca
element_core = 5
else if( zval <= 48.00001d0 ) then ! up to Sr
element_core = 9
else
write(*,*) '!!! not imlemented in element_core !!!'
stop
endif
endif
end function element_core
function element_covalent_radius(zatom)
! Return covalent radius of an atom
implicit none
include 'parameters.h'
integer,intent(in) :: zatom
double precision :: element_covalent_radius
!
! Data from Cambridge Structural Database
! http://en.wikipedia.org/wiki/Covalent_radius
!
! Values are first given in picometer
! They will be converted in bohr just after
select case(zatom)
case( 1)
element_covalent_radius = 31.
case( 2)
element_covalent_radius = 28.
case( 3)
element_covalent_radius = 128.
case( 4)
element_covalent_radius = 96.
case( 5)
element_covalent_radius = 84.
case( 6)
element_covalent_radius = 73.
case( 7)
element_covalent_radius = 71.
case( 8)
element_covalent_radius = 66.
case( 9)
element_covalent_radius = 57.
case(10) ! Ne.
element_covalent_radius = 58.
case(11)
element_covalent_radius = 166.
case(12)
element_covalent_radius = 141.
case(13)
element_covalent_radius = 121.
case(14)
element_covalent_radius = 111.
case(15)
element_covalent_radius = 107.
case(16)
element_covalent_radius = 105.
case(17)
element_covalent_radius = 102.
case(18) ! Ar.
element_covalent_radius = 106.
case(19)
element_covalent_radius = 203.
case(20)
element_covalent_radius = 176.
case(21)
element_covalent_radius = 170.
case(22)
element_covalent_radius = 160.
case(23)
element_covalent_radius = 153.
case(24)
element_covalent_radius = 139.
case(25)
element_covalent_radius = 145.
case(26)
element_covalent_radius = 145.
case(27)
element_covalent_radius = 140.
case(28)
element_covalent_radius = 124.
case(29)
element_covalent_radius = 132.
case(30)
element_covalent_radius = 122.
case(31)
element_covalent_radius = 120.
case(32)
element_covalent_radius = 119.
case(34)
element_covalent_radius = 120.
case(35)
element_covalent_radius = 120.
case(36) ! Kr.
element_covalent_radius = 116.
case default
write(*,*) '!!! covalent radius not available !!!'
stop
end select
! pm to bohr conversion
element_covalent_radius = element_covalent_radius*pmtoau
end function element_covalent_radius

79
src/xcDFT/exchange_energy.f90
View File

@ -1,79 +0,0 @@
function exchange_energy(rung,nGrid,weight,nBas,P,FxHF,rho,drho) result(Ex)
! Compute the exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: FxHF(nBas,nBas)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
double precision :: ExLDA,ExGGA,ExHF
double precision :: cX,aX,aC
double precision :: Ex
! Output variables
! Memory allocation
Ex = 0d0
ExLDA = 0d0
ExGGA = 0d0
ExHF = 0d0
select case (rung)
! Hartree calculation
case(0)
Ex = 0d0
! LDA functionals
case(1)
call lda_exchange_energy(nGrid,weight,rho,ExLDA)
Ex = ExLDA
! GGA functionals
case(2)
call gga_exchange_energy(nGrid,weight,rho,drho,ExGGA)
Ex = ExGGA
! Hybrid functionals
case(4)
cX = 0.20d0
aX = 0.72d0
aC = 0.81d0
call lda_exchange_energy(nGrid,weight,rho,ExLDA)
call gga_exchange_energy(nGrid,weight,rho,drho,ExGGA)
call fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = ExLDA &
+ cX*(ExHF - ExLDA) &
+ aX*(ExGGA - ExLDA)
! Hartree-Fock calculation
case(666)
call fock_exchange_energy(nBas,P,FxHF,ExHF)
Ex = ExHF
end select
end function exchange_energy

82
src/xcDFT/exchange_potential.f90
View File

@ -1,82 +0,0 @@
subroutine exchange_potential(rung,nGrid,weight,nBas,P,ERI,AO,dAO,rho,drho,Fx,FxHF)
! Compute the exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
double precision,allocatable :: FxLDA(:,:),FxGGA(:,:)
double precision :: cX,aX,aC
! Output variables
double precision,intent(out) :: Fx(nBas,nBas),FxHF(nBas,nBas)
! Memory allocation
allocate(FxLDA(nBas,nBas),FxGGA(nBas,nBas))
FxLDA(:,:) = 0d0
FxGGA(:,:) = 0d0
select case (rung)
! Hartree calculation
case(0)
Fx(:,:) = 0d0
! LDA functionals
case(1)
call lda_exchange_potential(nGrid,weight,nBas,AO,rho,FxLDA)
Fx(:,:) = FxLDA(:,:)
! GGA functionals
case(2)
call gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
Fx(:,:) = FxGGA(:,:)
! Hybrid functionals
case(4)
cX = 0.20d0
aX = 0.72d0
aC = 0.81d0
call lda_exchange_potential(nGrid,weight,nBas,AO,rho,FxLDA)
call gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA)
call fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = FxLDA(:,:) &
+ cX*(FxHF(:,:) - FxLDA(:,:)) &
+ aX*(FxGGA(:,:) - FxLDA(:,:))
! Hartree-Fock calculation
case(666)
call fock_exchange_potential(nBas,P,ERI,FxHF)
Fx(:,:) = FxHF(:,:)
end select
end subroutine exchange_potential

25
src/xcDFT/fock_exchange_energy.f90
View File

@ -1,25 +0,0 @@
subroutine fock_exchange_energy(nBas,P,Fx,Ex)
! Compute the (exact) Fock exchange energy
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: Fx(nBas,nBas)
! Local variables
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: Ex
! Compute HF exchange energy
Ex = trace_matrix(nBas,matmul(P,Fx))
end subroutine fock_exchange_energy

34
src/xcDFT/fock_exchange_potential.f90
View File

@ -1,34 +0,0 @@
subroutine fock_exchange_potential(nBas,P,ERI,Fx)
! Compute the Fock exchange potential
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Compute HF exchange matrix
Fx(:,:) = 0d0
do nu=1,nBas
do si=1,nBas
do la=1,nBas
do mu=1,nBas
Fx(mu,nu) = Fx(mu,nu) - P(la,si)*ERI(mu,la,si,nu)
enddo
enddo
enddo
enddo
end subroutine fock_exchange_potential

32
src/xcDFT/generate_shell.f90
View File

@ -1,32 +0,0 @@
subroutine generate_shell(atot,nShellFunction,ShellFunction)
! Generate shells for a given total angular momemtum
implicit none
! Input variables
integer,intent(in) :: atot,nShellFunction
! Local variables
integer :: ax,ay,az,ia
! Output variables
integer,intent(out) :: ShellFunction(nShellFunction,3)
ia = 0
do ax=atot,0,-1
do az=0,atot
ay = atot - ax - az
if(ay >= 0) then
ia = ia + 1
ShellFunction(ia,1) = ax
ShellFunction(ia,2) = ay
ShellFunction(ia,3) = az
endif
enddo
enddo
end subroutine generate_shell

44
src/xcDFT/gga_exchange_energy.f90
View File

@ -1,44 +0,0 @@
subroutine gga_exchange_energy(nGrid,weight,rho,drho,Ex)
! Compute GGA exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
integer :: iG
double precision :: alpha,beta
double precision :: r,g
! Output variables
double precision :: Ex
! Coefficients for G96 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = 1d0/137d0
! Compute GGA exchange energy
Ex = 0d0
do iG=1,nGrid
r = rho(iG)
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
Ex = Ex + weight(iG)*r**(4d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2)
enddo
Ex = 2d0*Ex
end subroutine gga_exchange_energy

62
src/xcDFT/gga_exchange_potential.f90
View File

@ -1,62 +0,0 @@
subroutine gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx)
! Compute GGA exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
double precision,intent(in) :: drho(3,nGrid)
! Local variables
double precision,parameter :: thresh = 1d-15
integer :: mu,nu,iG
double precision :: alpha,beta
double precision :: r,g,vAO,gAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Coefficients for G96 GGA exchange functional
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
beta = +1d0/137d0
beta = 0d0
! Compute GGA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = rho(iG)
g = drho(1,iG)**2 + drho(2,iG)**2 + drho(3,iG)**2
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*(4d0/3d0*r**(1d0/3d0)*(alpha - beta*g**(3d0/4d0)/r**2) &
+ 2d0*beta*g**(3d0/4d0)/r**(5d0/3d0))
gAO = drho(1,iG)*(dAO(1,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(1,nu,iG)) &
+ drho(2,iG)*(dAO(2,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(2,nu,iG)) &
+ drho(3,iG)*(dAO(3,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(3,nu,iG))
gAO = weight(iG)*gAO
Fx(mu,nu) = Fx(mu,nu) - 2d0*gAO*3d0/4d0*beta*g**(-1d0/4d0)/r**(2d0/3d0)
enddo
enddo
enddo
end subroutine gga_exchange_potential

45
src/xcDFT/gradient_density.f90
View File

@ -1,45 +0,0 @@
subroutine gradient_density(nGrid,nBas,P,AO,dAO,drho)
! Calculate gradient of the one-electron density
implicit none
include 'parameters.h'
! Input variables
double precision,parameter :: thresh = 1d-15
integer,intent(in) :: nGrid
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
! Local variables
integer :: ixyz,iG,mu,nu
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: drho(3,nGrid)
drho(:,:) = 0d0
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
do ixyz=1,3
drho(ixyz,iG) = drho(ixyz,iG) &
+ P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG))
enddo
enddo
enddo
enddo
do iG=1,nGrid
do ixyz=1,3
if(abs(drho(ixyz,iG)) < thresh) drho(ixyz,iG) = thresh
enddo
enddo
end subroutine gradient_density

33
src/xcDFT/hartree_coulomb.f90
View File

@ -1,33 +0,0 @@
subroutine hartree_coulomb(nBas,P,ERI,J)
! Compute Coulomb matrix
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas)
! Local variables
integer :: mu,nu,la,si
! Output variables
double precision,intent(out) :: J(nBas,nBas)
J = 0d0
do mu=1,nBas
do nu=1,nBas
do la=1,nBas
do si=1,nBas
J(mu,nu) = J(mu,nu) + P(la,si)*ERI(mu,la,nu,si)
enddo
enddo
enddo
enddo
end subroutine hartree_coulomb

36
src/xcDFT/lda_exchange_energy.f90
View File

@ -1,36 +0,0 @@
subroutine lda_exchange_energy(nGrid,weight,rho,Ex)
! Compute LDA exchange energy
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: iG
double precision :: alpha
! Output variables
double precision :: Ex
! Cx coefficient for Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
! Compute LDA exchange energy
Ex = 0d0
do iG=1,nGrid
Ex = Ex + weight(iG)*alpha*rho(iG)**(4d0/3d0)
enddo
Ex = 2d0*Ex
end subroutine lda_exchange_energy

46
src/xcDFT/lda_exchange_potential.f90
View File

@ -1,46 +0,0 @@
subroutine lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx)
! Compute LDA exchange potential
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
double precision,intent(in) :: weight(nGrid)
integer,intent(in) :: nBas
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: rho(nGrid)
! Local variables
integer :: mu,nu,iG
double precision :: alpha
double precision :: r,vAO
! Output variables
double precision,intent(out) :: Fx(nBas,nBas)
! Cx coefficient for Slater LDA exchange
alpha = -(3d0/2d0)*(3d0/(4d0*pi))**(1d0/3d0)
! Compute LDA exchange matrix in the AO basis
Fx(:,:) = 0d0
do mu=1,nBas
do nu=1,nBas
do iG=1,nGrid
r = rho(iG)
vAO = weight(iG)*AO(mu,iG)*AO(nu,iG)
Fx(mu,nu) = Fx(mu,nu) &
+ vAO*4d0/3d0*alpha*r**(1d0/3d0)
enddo
enddo
enddo
end subroutine lda_exchange_potential

BIN
src/xcDFT/obj/AO_values_grid.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/DIIS_extrapolation.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/NormCoeff.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/RKS.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/density.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/dft_grid.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/electron_number.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/elements.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/exchange_energy.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/exchange_potential.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/fock_exchange_energy.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/fock_exchange_potential.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/generate_shell.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/gga_exchange_energy.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/gga_exchange_potential.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/gradient_density.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/hartree_coulomb.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/lda_exchange_energy.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/lda_exchange_potential.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/one_electron_density.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/orthogonalization_matrix.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/print_RKS.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/quadrature_grid.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/read_basis.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/read_geometry.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/read_grid.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/read_integrals.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/read_molecule.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/read_options.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/select_rung.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/utils.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/wrap_lapack.o
View File

Binary file not shown.

BIN
src/xcDFT/obj/xcDFT.o
View File

Binary file not shown.

47
src/xcDFT/one_electron_density.f90
View File

@ -1,47 +0,0 @@
subroutine one_electron_density(nGrid,nBas,P,AO,dAO,rho,drho)
! Calculate one-electron density
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nGrid
integer,intent(in) :: nBas
double precision,intent(in) :: P(nBas,nBas)
double precision,intent(in) :: AO(nBas,nGrid)
double precision,intent(in) :: dAO(3,nBas,nGrid)
! Local variables
integer :: ixyz,iG,mu,nu
double precision,external :: trace_matrix
! Output variables
double precision,intent(out) :: rho(nGrid)
double precision,intent(out) :: drho(3,nGrid)
rho(:) = 0d0
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
rho(iG) = rho(iG) + AO(mu,iG)*P(mu,nu)*AO(nu,iG)
enddo
enddo
enddo
drho(:,:) = 0d0
do ixyz=1,3
do iG=1,nGrid
do mu=1,nBas
do nu=1,nBas
drho(ixyz,iG) = drho(ixyz,iG) &
+ P(mu,nu)*(dAO(ixyz,mu,iG)*AO(nu,iG) + AO(mu,iG)*dAO(ixyz,nu,iG))
enddo
enddo
enddo
enddo
end subroutine one_electron_density

63
src/xcDFT/orthogonalization_matrix.f90
View File

@ -1,63 +0,0 @@
subroutine orthogonalization_matrix(nBas,S,X)
! Compute the orthogonalization matrix X = S^(-1/2)
implicit none
! Input variables
integer,intent(in) :: nBas
double precision,intent(in) :: S(nBas,nBas)
! Local variables
logical :: debug
double precision,allocatable :: UVec(:,:),Uval(:)
double precision,parameter :: thresh = 1d-6
integer :: i
! Output variables
double precision,intent(out) :: X(nBas,nBas)
debug = .false.
allocate(Uvec(nBas,nBas),Uval(nBas))
write(*,*)
write(*,*) ' *** Lowdin orthogonalization X = S^(-1/2) *** '
write(*,*)
Uvec = S
call diagonalize_matrix(nBas,Uvec,Uval)
do i=1,nBas
if(Uval(i) > thresh) then
Uval(i) = 1d0/sqrt(Uval(i))
else
write(*,*) 'Eigenvalue',i,'too small for Lowdin orthogonalization'
endif
enddo
call ADAt(nBas,Uvec,Uval,X)
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Orthogonalization matrix'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,X)
write(*,*)
endif
end subroutine orthogonalization_matrix

61
src/xcDFT/print_RKS.f90
View File

@ -1,61 +0,0 @@
subroutine print_RKS(nBas,nO,e,C,ENuc,ET,EV,EJ,Ex,Ec,EKS)
! Print one- and two-electron energies and other stuff for RKS calculation
implicit none
include 'parameters.h'
integer,intent(in) :: nBas,nO
double precision,intent(in) :: e(nBas),c(nBas,nBas),ENuc,ET,EV,EJ,Ex,Ec,EKS
integer :: HOMO,LUMO
double precision :: Gap
! HOMO and LUMO
HOMO = nO
LUMO = HOMO + 1
Gap = e(LUMO) - e(HOMO)
! Dump results
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A32)') ' Summary '
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10)') ' One-electron energy ',ET + EV
write(*,'(A32,1X,F16.10)') ' Kinetic energy ',ET
write(*,'(A32,1X,F16.10)') ' Potential energy ',EV
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10)') ' Two-electron energy ',EJ + Ex + Ec
write(*,'(A32,1X,F16.10)') ' Coulomb energy ',EJ
write(*,'(A32,1X,F16.10)') ' Exchange energy ',Ex
write(*,'(A32,1X,F16.10)') ' Correlation energy ',Ec
write(*,'(A50)') '---------------------------------------'
write(*,'(A32,1X,F16.10)') ' Electronic energy ',EKS
write(*,'(A32,1X,F16.10)') ' Nuclear repulsion ',ENuc
write(*,'(A32,1X,F16.10)') ' Kohn-Sham energy ',EKS + ENuc
write(*,'(A50)') '---------------------------------------'
write(*,'(A36,F13.6)') ' KS HOMO energy (eV):',e(HOMO)*HatoeV
write(*,'(A36,F13.6)') ' KS LUMO energy (eV):',e(LUMO)*Hatoev
write(*,'(A36,F13.6)') ' KS HOMO-LUMO gap (eV):',Gap*Hatoev
write(*,'(A50)') '---------------------------------------'
write(*,*)
! Print results
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') 'Kohn-Sham orbital coefficients '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,nBas,C)
write(*,*)
write(*,'(A50)') '---------------------------------------'
write(*,'(A50)') ' Kohn-Sham orbital energies '
write(*,'(A50)') '---------------------------------------'
call matout(nBas,1,e)
write(*,*)
end subroutine print_RKS

77
src/xcDFT/quadrature_grid.f90
View File

@ -1,77 +0,0 @@
subroutine quadrature_grid(nRad,nAng,nGrid,root,weight)
! Build roots and weights of quadrature grid
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nRad,nAng,nGrid
! Local variables
integer :: i,j,k
double precision :: scale
double precision,allocatable :: Radius(:)
double precision,allocatable :: RadWeight(:)
double precision,allocatable :: XYZ(:,:)
double precision,allocatable :: XYZWeight(:)
! Output variables
double precision,intent(out) :: root(3,nGrid)
double precision,intent(out) :: weight(nGrid)
! Memory allocation
allocate(Radius(nRad),RadWeight(nRad),XYZ(3,nAng),XYZWeight(nAng))
! Findthe radial grid
scale = 1d0
call EulMac(Radius,RadWeight,nRad,scale)
write(*,20)
write(*,30)
write(*,20)
do i = 1,nRad
write(*,40) i,Radius(i),RadWeight(i)
end do
write(*,20)
write(*,*)
! Find the angular grid
call Lebdev(XYZ,XYZWeight,nAng)
write(*,20)
write(*,50)
write(*,20)
do j = 1,nAng
write(*,60) j,(XYZ(k,j),k=1,3), XYZWeight(j)
end do
write(*,20)
! Form the roots and weights
k = 0
do i=1,nRad
do j=1,nAng
k = k + 1
root(:,k) = Radius(i)*XYZ(:,j)
weight(k) = RadWeight(i)*XYZWeight(j)
enddo
enddo
! Compute values of the basis functions (and the its gradient if required) at each grid point
20 format(T2,58('-'))
30 format(T20,'Radial Quadrature',/, &
T6,'I',T26,'Radius',T50,'Weight')
40 format(T3,I4,T18,F17.10,T35,F25.10)
50 format(T20,'Angular Quadrature',/, &
T6,'I',T19,'X',T29,'Y',T39,'Z',T54,'Weight')
60 format(T3,I4,T13,3F10.5,T50,F10.5)
end subroutine quadrature_grid

117
src/xcDFT/read_basis.f90
View File

@ -1,117 +0,0 @@
subroutine read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
! Read basis set information
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: nAt,nO
double precision,intent(in) :: rAt(nAt,3)
! Local variables
integer :: nShAt,iAt,iShell
integer :: i,j,k
character :: shelltype
! Output variables
integer,intent(out) :: nShell,nBas,nV
double precision,intent(out) :: CenterShell(maxShell,3)
integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell)
double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
!------------------------------------------------------------------------
! Primary basis set information
!------------------------------------------------------------------------
! Open file with basis set specification
open(unit=2,file='input/basis')
! Read basis information
write(*,'(A28)') 'Gaussian basis set'
write(*,'(A28)') '------------------'
nShell = 0
do i=1,nAt
read(2,*) iAt,nShAt
write(*,'(A28,1X,I16)') 'Atom n. ',iAt
write(*,'(A28,1X,I16)') 'number of shells ',nShAt
write(*,'(A28)') '------------------'
! Basis function centers
do j=1,nShAt
nShell = nShell + 1
do k=1,3
CenterShell(nShell,k) = rAt(iAt,k)
enddo
! Shell type and contraction degree
read(2,*) shelltype,KShell(nShell)
if(shelltype == "S") then
TotAngMomShell(nShell) = 0
write(*,'(A28,1X,I16)') 's-type shell with K = ',KShell(nShell)
elseif(shelltype == "P") then
TotAngMomShell(nShell) = 1
write(*,'(A28,1X,I16)') 'p-type shell with K = ',KShell(nShell)
elseif(shelltype == "D") then
TotAngMomShell(nShell) = 2
write(*,'(A28,1X,I16)') 'd-type shell with K = ',KShell(nShell)
elseif(shelltype == "F") then
TotAngMomShell(nShell) = 3
write(*,'(A28,1X,I16)') 'f-type shell with K = ',KShell(nShell)
elseif(shelltype == "G") then
TotAngMomShell(nShell) = 4
write(*,'(A28,1X,I16)') 'g-type shell with K = ',KShell(nShell)
elseif(shelltype == "H") then
TotAngMomShell(nShell) = 5
write(*,'(A28,1X,I16)') 'h-type shell with K = ',KShell(nShell)
elseif(shelltype == "I") then
TotAngMomShell(nShell) = 6
write(*,'(A28,1X,I16)') 'i-type shell with K = ',KShell(nShell)
endif
! Read exponents and contraction coefficients
write(*,'(A28,1X,A16,A16)') '','Exponents','Contraction'
do k=1,Kshell(nShell)
read(2,*) ExpShell(nShell,k),DShell(nShell,k)
write(*,'(A28,1X,F16.10,F16.10)') '',ExpShell(nShell,k),DShell(nShell,k)
enddo
enddo
write(*,'(A28)') '------------------'
enddo
! Total number of shells
write(*,'(A28,1X,I16)') 'Number of shells',nShell
write(*,'(A28)') '------------------'
write(*,*)
! Close file with basis set specification
close(unit=2)
! Calculate number of basis functions
nBas = 0
do iShell=1,nShell
nBas = nBas + (TotAngMomShell(iShell)*TotAngMomShell(iShell) + 3*TotAngMomShell(iShell) + 2)/2
enddo
write(*,'(A28)') '------------------'
write(*,'(A28,1X,I16)') 'Number of basis functions',NBas
write(*,'(A28)') '------------------'
write(*,*)
! Number of virtual orbitals
nV = nBas - nO
end subroutine read_basis

58
src/xcDFT/read_geometry.f90
View File

@ -1,58 +0,0 @@
subroutine read_geometry(nAt,ZNuc,rA,ENuc)
! Read molecular geometry
implicit none
! Ouput variables
integer,intent(in) :: nAt
! Local variables
integer :: i,j
double precision :: RAB
! Ouput variables
double precision,intent(out) :: ZNuc(NAt),rA(nAt,3),ENuc
! Open file with geometry specification
open(unit=1,file='input/molecule')
! Read number of atoms
read(1,*)
read(1,*)
read(1,*)
do i=1,nAt
read(1,*) ZNuc(i),rA(i,1),rA(i,2),rA(i,3)
enddo
! Compute nuclear repulsion energy
ENuc = 0
do i=1,nAt-1
do j=i+1,nAt
RAB = (rA(i,1)-rA(j,1))**2 + (rA(i,2)-rA(j,2))**2 + (rA(i,3)-rA(j,3))**2
ENuc = ENuc + ZNuc(i)*ZNuc(j)/sqrt(RAB)
enddo
enddo
! Close file with geometry specification
close(unit=1)
! Print geometry
write(*,'(A28)') '------------------'
write(*,'(A28)') 'Molecular geometry'
write(*,'(A28)') '------------------'
do i=1,NAt
write(*,'(A28,1X,I16)') 'Atom n. ',i
write(*,'(A28,1X,F16.10)') 'Z = ',ZNuc(i)
write(*,'(A28,1X,F16.10,F16.10,F16.10)') 'Atom coordinates:',(rA(i,j),j=1,3)
enddo
write(*,*)
write(*,'(A28)') '------------------'
write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc
write(*,'(A28)') '------------------'
write(*,*)
end subroutine read_geometry

47
src/xcDFT/read_grid.f90
View File

@ -1,47 +0,0 @@
subroutine read_grid(SGn,nRad,nAng,nGrid)
! Read grid type
implicit none
! Input variables
integer,intent(in) :: SGn
! Output variables
integer,intent(out) :: nRad
integer,intent(out) :: nAng
integer,intent(out) :: nGrid
write(*,*)'----------------------------------------------------------'
write(*,'(A22,I1)')' Quadrature grid: SG-',SGn
write(*,*)'----------------------------------------------------------'
select case (SGn)
case(0)
nRad = 23
nAng = 170
case(1)
nRad = 50
nAng = 194
case(2)
nRad = 75
nAng = 302
case(3)
nRad = 99
nAng = 590
case default
write(*,*) '!!! Quadrature grid not available !!!'
stop
end select
nGrid = nRad*nAng
end subroutine read_grid

114
src/xcDFT/read_integrals.f90
View File

@ -1,114 +0,0 @@
subroutine read_integrals(nBas,S,T,V,Hc,G)
! Read one- and two-electron integrals from files
implicit none
! Input variables
integer,intent(in) :: nBas
! Local variables
logical :: debug
integer :: mu,nu,la,si
double precision :: Ov,Kin,Nuc,ERI
! Output variables
double precision,intent(out) :: S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),G(nBas,nBas,nBas,nBas)
! Open file with integrals
debug = .false.
open(unit=8 ,file='int/Ov.dat')
open(unit=9 ,file='int/Kin.dat')
open(unit=10,file='int/Nuc.dat')
open(unit=11,file='int/ERI.dat')
! Read overlap integrals
S = 0d0
do
read(8,*,end=8) mu,nu,Ov
S(mu,nu) = Ov
enddo
8 close(unit=8)
! Read kinetic integrals
T = 0d0
do
read(9,*,end=9) mu,nu,Kin
T(mu,nu) = Kin
enddo
9 close(unit=9)
! Read nuclear integrals
V = 0d0
do
read(10,*,end=10) mu,nu,Nuc
V(mu,nu) = Nuc
enddo
10 close(unit=10)
! Define core Hamiltonian
Hc = T + V
! Read nuclear integrals
G = 0d0
do
read(11,*,end=11) mu,nu,la,si,ERI
! <12|34>
G(mu,nu,la,si) = ERI
! <32|14>
G(la,nu,mu,si) = ERI
! <14|32>
G(mu,si,la,nu) = ERI
! <34|12>
G(la,si,mu,nu) = ERI
! <41|23>
G(si,mu,nu,la) = ERI
! <23|41>
G(nu,la,si,mu) = ERI
! <21|43>
G(nu,mu,si,la) = ERI
! <43|21>
G(si,la,nu,mu) = ERI
enddo
11 close(unit=11)
! Print results
if(debug) then
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Overlap integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,S)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Kinetic integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,T)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Nuclear integrals'
write(*,'(A28)') '----------------------'
call matout(nBas,nBas,V)
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28)') 'Electron repulsion integrals'
write(*,'(A28)') '----------------------'
do la=1,nBas
do si=1,nBas
call matout(nBas,nBas,G(1,1,la,si))
enddo
enddo
write(*,*)
endif
end subroutine read_integrals

42
src/xcDFT/read_molecule.f90
View File

@ -1,42 +0,0 @@
subroutine read_molecule(nAt,nEl,nO)
! Read number of atoms nAt and number of electrons nEl
implicit none
! Input variables
integer,intent(out) :: nAt,nEl,nO
! Open file with geometry specification
open(unit=1,file='input/molecule')
! Read number of atoms and number of electrons
read(1,*)
read(1,*) nAt,nEl
! Number of occupied orbitals
if(mod(nEl,2) /= 0) then
write(*,*) 'closed-shell system required!'
stop
endif
nO = nEl/2
! Print results
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of atoms',nAt
write(*,'(A28)') '----------------------'
write(*,*)
write(*,'(A28)') '----------------------'
write(*,'(A28,1X,I16)') 'Number of electrons',nEl
write(*,'(A28)') '----------------------'
write(*,*)
! Close file with geometry specification
close(unit=1)
end subroutine read_molecule

31
src/xcDFT/read_options.f90
View File

@ -1,31 +0,0 @@
subroutine read_options(rung,SGn)
! Read DFT options
implicit none
! Input variables
integer,intent(out) :: rung
integer,intent(out) :: SGn
! Open file with method specification
open(unit=1,file='input/options')
! Default values
rung = 1
SGn = 0
! Read rung of Jacob's ladder
read(1,*)
read(1,*) rung
! Read SG-n grid
read(1,*)
read(1,*) SGn
end subroutine read_options

45
src/xcDFT/select_rung.f90
View File

@ -1,45 +0,0 @@
subroutine select_rung(rung)
! Select rung of Jacob's ladder
implicit none
include 'parameters.h'
! Input variables
integer,intent(in) :: rung
select case (rung)
! Hartree calculation
case(0)
write(*,*) " *** 0th rung of Jacob's ladder: Hartree calculation *** "
! LDA functionals
case(1)
write(*,*) " *** 1st rung of Jacob's ladder: local-density approximation (LDA) *** "
! GGA functionals
case(2)
write(*,*) " *** 2nd rung of Jacob's ladder: generalized gradient approximation (GGA) *** "
! meta-GGA functionals
case(3)
write(*,*) " *** 3rd rung of Jacob's ladder: meta-GGA functional (MGGA) *** "
! Hybrid functionals
case(4)
write(*,*) " *** 4th rung of Jacob's ladder: hybrid functional *** "
! Hartree-Fock calculation
case(666)
write(*,*) " *** rung 666: Hartree-Fock calculation *** "
! Default
case default
write(*,*) "!!! rung not available !!!"
stop
end select
end subroutine select_rung

246
src/xcDFT/utils.f90
View File

@ -1,246 +0,0 @@
!------------------------------------------------------------------------
subroutine matout(m,n,A)
! Print the MxN array A
implicit none
integer,parameter :: ncol = 5
double precision,parameter :: small = 1d-10
integer,intent(in) :: m,n
double precision,intent(in) :: A(m,n)
double precision :: B(ncol)
integer :: ilower,iupper,num,i,j
do ilower=1,n,ncol
iupper = min(ilower + ncol - 1,n)
num = iupper - ilower + 1
write(*,'(3X,10(9X,I6))') (j,j=ilower,iupper)
do i=1,m
do j=ilower,iupper
B(j-ilower+1) = A(i,j)
enddo
do j=1,num
if(abs(B(j)) < small) B(j) = 0d0
enddo
write(*,'(I7,10F15.8)') i,(B(j),j=1,num)
enddo
enddo
end subroutine matout
!------------------------------------------------------------------------
function trace_matrix(n,A) result(Tr)
! Calculate the trace of the square matrix A
implicit none
! Input variables
integer,intent(in) :: n
double precision,intent(in) :: A(n,n)
! Local variables
integer :: i
! Output variables
double precision :: Tr
Tr = 0d0
do i=1,n
Tr = Tr + A(i,i)
enddo
end function trace_matrix
!------------------------------------------------------------------------
subroutine prepend(N,M,A,b)
! Prepend the vector b of size N into the matrix A of size NxM
implicit none
! Input variables
integer,intent(in) :: N,M
double precision,intent(in) :: b(N)
! Local viaruabkes
integer :: i,j
! Output variables
double precision,intent(out) :: A(N,M)
! print*,'b in append'
! call matout(N,1,b)
do i=1,N
do j=M-1,1,-1
A(i,j+1) = A(i,j)
enddo
A(i,1) = b(i)
enddo
end subroutine prepend
!------------------------------------------------------------------------
subroutine append(N,M,A,b)
! Append the vector b of size N into the matrix A of size NxM
implicit none
! Input variables
integer,intent(in) :: N,M
double precision,intent(in) :: b(N)
! Local viaruabkes
integer :: i,j
! Output variables
double precision,intent(out) :: A(N,M)
do i=1,N
do j=2,M
A(i,j-1) = A(i,j)
enddo
A(i,M) = b(i)
enddo
end subroutine append
!------------------------------------------------------------------------
subroutine AtDA(N,A,D,B)
! Perform B = At.D.A where A is a NxN matrix and D is a diagonal matrix given
! as a vector of length N
implicit none
! Input variables
integer,intent(in) :: N
double precision,intent(in) :: A(N,N),D(N)
! Local viaruabkes
integer :: i,j,k
! Output variables
double precision,intent(out) :: B(N,N)
B = 0d0
do i=1,N
do j=1,N
do k=1,N
B(i,k) = B(i,k) + A(j,i)*D(j)*A(j,k)
enddo
enddo
enddo
end subroutine AtDA
!------------------------------------------------------------------------
subroutine ADAt(N,A,D,B)
! Perform B = A.D.At where A is a NxN matrix and D is a diagonal matrix given
! as a vector of length N
implicit none
! Input variables
integer,intent(in) :: N
double precision,intent(in) :: A(N,N),D(N)
! Local viaruabkes
integer :: i,j,k
! Output variables
double precision,intent(out) :: B(N,N)
B = 0d0
do i=1,N
do j=1,N
do k=1,N
B(i,k) = B(i,k) + A(i,j)*D(j)*A(k,j)
enddo
enddo
enddo
end subroutine ADAt
!------------------------------------------------------------------------
subroutine DA(N,D,A)
! Perform A <- D.A where A is a NxN matrix and D is a diagonal matrix given
! as a vector of length N
implicit none
integer,intent(in) :: N
integer :: i,j
double precision,intent(in) :: D(N)
double precision,intent(inout):: A(N,N)
do i=1,N
do j=1,N
A(i,j) = D(i)*A(i,j)
enddo
enddo
end subroutine DA
!------------------------------------------------------------------------
subroutine AD(N,A,D)
! Perform A <- A.D where A is a NxN matrix and D is a diagonal matrix given
! as a vector of length N
implicit none
integer,intent(in) :: N
integer :: i,j
double precision,intent(in) :: D(N)
double precision,intent(inout):: A(N,N)
do i=1,N
do j=1,N
A(i,j) = A(i,j)*D(j)
enddo
enddo
end subroutine AD
!------------------------------------------------------------------------
recursive function fac(n) result(fact)
implicit none
integer :: fact
integer, intent(in) :: n
if (n == 0) then
fact = 1
else
fact = n * fac(n-1)
end if
end function fac
function dfac(n) result(fact)
implicit none
double precision :: fact
integer, intent(in) :: n
integer :: fac
fact = dble(fac(n))
end function dfac

147
src/xcDFT/wrap_lapack.f90
View File

@ -1,147 +0,0 @@
subroutine diagonalize_matrix(N,A,e)
! Diagonalize a square matrix
implicit none
! Input variables
integer,intent(in) :: N
double precision,intent(inout):: A(N,N)
double precision,intent(out) :: e(N)
! Local variables
integer :: lwork,info
double precision,allocatable :: work(:)
! Memory allocation
allocate(work(3*N))
lwork = size(work)
call dsyev('V','U',N,A,N,e,work,lwork,info)
if(info /= 0) then
print*,'Problem in diagonalize_matrix (dsyev)!!'
stop
endif
end subroutine diagonalize_matrix
subroutine svd(N,A,U,D,Vt)
! Compute A = U.D.Vt
! Dimension of A is NxN
implicit none
integer, intent(in) :: N
double precision,intent(in) :: A(N,N)
double precision,intent(out) :: U(N,N)
double precision,intent(out) :: Vt(N,N)
double precision,intent(out) :: D(N)
double precision,allocatable :: work(:)
integer :: info,lwork
double precision,allocatable :: scr(:,:)
allocate (scr(N,N))
scr(:,:) = A(:,:)
! Find optimal size for temporary arrays
allocate(work(1))
lwork = -1
call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info)
lwork = int(work(1))
deallocate(work)
allocate(work(lwork))
call dgesvd('A','A',N,N,scr,N,D,U,N,Vt,N,work,lwork,info)
deallocate(work,scr)
if (info /= 0) then
print *, info, ': SVD failed'
stop
endif
end
subroutine inverse_matrix(N,A,B)
! Returns the inverse of the square matrix A in B
implicit none
integer,intent(in) :: N
double precision, intent(in) :: A(N,N)
double precision, intent(out) :: B(N,N)
integer :: info,lwork
integer, allocatable :: ipiv(:)
double precision,allocatable :: work(:)
allocate (ipiv(N),work(N*N))
lwork = size(work)
B(1:N,1:N) = A(1:N,1:N)
call dgetrf(N,N,B,N,ipiv,info)
if (info /= 0) then
print*,info
stop 'error in inverse (dgetrf)!!'
endif
call dgetri(N,B,N,ipiv,work,lwork,info)
if (info /= 0) then
print *, info
stop 'error in inverse (dgetri)!!'
endif
deallocate(ipiv,work)
end subroutine inverse_matrix
subroutine linear_solve(N,A,b,x)
! Solve the linear system A.x = b where A is a NxN matrix
! and x and x are vectors of size N
implicit none
integer,intent(in) :: N
double precision,intent(in) :: A(N,N),b(N)
double precision,intent(out) :: x(N)
integer :: info,lwork
integer,allocatable :: ipiv(:)
double precision,allocatable :: work(:)
allocate(ipiv(N),work(N*N))
lwork = size(work)
x = b
call dsysv('U',N,1,A,N,ipiv,x,N,work,lwork,info)
if (info /= 0) then
print *, info
stop 'error in linear_solve (dsysv)!!'
endif
end subroutine linear_solve

120
src/xcDFT/xcDFT.f90
View File

@ -1,120 +0,0 @@
program xcDFT
! exchange-correlation density-functional theory calculations
include 'parameters.h'
integer :: nAt,nBas,nEl,nO,nV
double precision :: ENuc,EKS
double precision,allocatable :: ZNuc(:),rAt(:,:)
integer :: nShell
integer,allocatable :: TotAngMomShell(:)
integer,allocatable :: KShell(:)
double precision,allocatable :: CenterShell(:,:)
double precision,allocatable :: DShell(:,:)
double precision,allocatable :: ExpShell(:,:)
double precision,allocatable :: S(:,:),T(:,:),V(:,:),Hc(:,:),X(:,:)
double precision,allocatable :: ERI(:,:,:,:)
integer :: rung
integer :: SGn
integer :: nRad,nAng,nGrid
double precision,allocatable :: root(:,:)
double precision,allocatable :: weight(:)
double precision,allocatable :: AO(:,:)
double precision,allocatable :: dAO(:,:,:)
double precision :: start_KS,end_KS,t_KS
! Hello World
write(*,*)
write(*,*) '********************************'
write(*,*) '* TCCM winter school 2008: DFT *'
write(*,*) '********************************'
write(*,*)
!------------------------------------------------------------------------
! Read input information
!------------------------------------------------------------------------
! Read number of atoms, number of electrons of the system
! nO = number of occupied orbitals
! nV = number of virtual orbitals (see below)
! nBas = number of basis functions (see below)
! = nO + nV
call read_molecule(nAt,nEl,nO)
allocate(ZNuc(nAt),rAt(nAt,3))
! Read geometry
call read_geometry(nAt,ZNuc,rAt,ENuc)
allocate(CenterShell(maxShell,3),TotAngMomShell(maxShell),KShell(maxShell), &
DShell(maxShell,maxK),ExpShell(maxShell,maxK))
!------------------------------------------------------------------------
! Read basis set information
!------------------------------------------------------------------------
call read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
!------------------------------------------------------------------------
! Read one- and two-electron integrals
!------------------------------------------------------------------------
! Memory allocation for one- and two-electron integrals
allocate(S(nBas,nBas),T(nBas,nBas),V(nBas,nBas),Hc(nBas,nBas),X(nBas,nBas), &
ERI(nBas,nBas,nBas,nBas))
! Read integrals
call read_integrals(nBas,S,T,V,Hc,ERI)
! Orthogonalization X = S^(-1/2)
call orthogonalization_matrix(nBas,S,X)
!------------------------------------------------------------------------
! DFT options
!------------------------------------------------------------------------
call read_options(rung,SGn)
!------------------------------------------------------------------------
! Construct quadrature grid
!------------------------------------------------------------------------
call read_grid(SGn,nRad,nAng,nGrid)
allocate(root(3,nGrid),weight(nGrid))
call quadrature_grid(nRad,nAng,nGrid,root,weight)
!------------------------------------------------------------------------
! Calculate AO values at grid points
!------------------------------------------------------------------------
allocate(AO(nBas,nGrid),dAO(3,nBas,nGrid))
call AO_values_grid(nBas,nShell,CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
nGrid,root,AO,dAO)
!------------------------------------------------------------------------
! Compute KS energy
!------------------------------------------------------------------------
call cpu_time(start_KS)
call RKS(rung,nGrid,weight,nBas,AO,dAO,nO,S,T,V,Hc,ERI,X,ENuc,EKS)
call cpu_time(end_KS)
t_KS = end_KS - start_KS
write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for KS = ',t_KS,' seconds'
write(*,*)
!------------------------------------------------------------------------
! End of xcDFT
!------------------------------------------------------------------------
end program xcDFT