mirror of
https://github.com/pfloos/quack
synced 2025-01-03 10:05:59 +01:00
Fixing IntPak
This commit is contained in:
parent
f4329480ba
commit
a3b5d63ca2
@ -1,4 +1,7 @@
|
|||||||
|
integer,parameter :: ncart = 3
|
||||||
integer,parameter :: nspin = 2
|
integer,parameter :: nspin = 2
|
||||||
|
integer,parameter :: nsp = 3
|
||||||
|
integer,parameter :: maxEns = 10
|
||||||
integer,parameter :: maxShell = 50
|
integer,parameter :: maxShell = 50
|
||||||
integer,parameter :: n1eInt = 3
|
integer,parameter :: n1eInt = 3
|
||||||
integer,parameter :: n2eInt = 4
|
integer,parameter :: n2eInt = 4
|
||||||
@ -6,6 +9,7 @@
|
|||||||
integer,parameter :: n4eInt = 3
|
integer,parameter :: n4eInt = 3
|
||||||
integer,parameter :: maxK = 20
|
integer,parameter :: maxK = 20
|
||||||
|
|
||||||
|
double precision,parameter :: threshold = 1d-15
|
||||||
double precision,parameter :: pi = acos(-1d0)
|
double precision,parameter :: pi = acos(-1d0)
|
||||||
double precision,parameter :: HaToeV = 27.21138602d0
|
double precision,parameter :: HaToeV = 27.21138602d0
|
||||||
double precision,parameter :: pmtoau = 0.0188973d0
|
double precision,parameter :: pmtoau = 0.0188973d0
|
||||||
|
|||||||
@ -2,6 +2,8 @@
|
|||||||
F
|
F
|
||||||
# Chemist notation for two-electron integral?
|
# Chemist notation for two-electron integral?
|
||||||
T
|
T
|
||||||
|
# Exposant of the Slater geminal
|
||||||
|
1.0
|
||||||
# One-electron integrals: Ov Kin Nuc
|
# One-electron integrals: Ov Kin Nuc
|
||||||
T T T
|
T T T
|
||||||
# Two-electron integrals: ERI F12 Yuk Erf
|
# Two-electron integrals: ERI F12 Yuk Erf
|
||||||
|
|||||||
@ -1,4 +1,4 @@
|
|||||||
# nAt nEl nCore nRyd
|
# nAt nEl nEla nElb
|
||||||
1 2 0 0
|
1 2 1 1
|
||||||
# Znuc x y z
|
# Znuc x y z
|
||||||
2.0 0.0 0.0 0.0
|
He 0.0 0.0 0.0
|
||||||
|
|||||||
@ -1,28 +0,0 @@
|
|||||||
subroutine CalcNBasis(nShell,atot,NBasis)
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: nShell
|
|
||||||
integer,intent(in) :: atot(nShell)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: iShell
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
integer,intent(out) :: NBasis
|
|
||||||
|
|
||||||
NBasis = 0
|
|
||||||
do iShell=1,nShell
|
|
||||||
NBasis = NBasis + (atot(iShell)*atot(iShell) + 3*atot(iShell) + 2)/2
|
|
||||||
enddo
|
|
||||||
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of basis functions',NBasis
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
end subroutine CalcNBasis
|
|
||||||
@ -18,11 +18,13 @@ program IntPak
|
|||||||
logical :: do3eInt(n3eInt)
|
logical :: do3eInt(n3eInt)
|
||||||
logical :: do4eInt(n4eInt)
|
logical :: do4eInt(n4eInt)
|
||||||
|
|
||||||
integer :: NAtoms,NBasis,iType
|
integer :: nNuc,nBas,iType
|
||||||
|
integer :: nEl,nO,nV
|
||||||
double precision :: ExpS
|
double precision :: ExpS
|
||||||
|
double precision :: ENuc
|
||||||
integer :: KG
|
integer :: KG
|
||||||
double precision,allocatable :: DG(:),ExpG(:)
|
double precision,allocatable :: DG(:),ExpG(:)
|
||||||
double precision,allocatable :: ZNuc(:),XYZAtoms(:,:)
|
double precision,allocatable :: ZNuc(:),rNuc(:,:)
|
||||||
|
|
||||||
integer :: nShell
|
integer :: nShell
|
||||||
integer,allocatable :: TotAngMomShell(:),KShell(:)
|
integer,allocatable :: TotAngMomShell(:),KShell(:)
|
||||||
@ -52,7 +54,7 @@ program IntPak
|
|||||||
|
|
||||||
! Read options for integral calculations
|
! Read options for integral calculations
|
||||||
|
|
||||||
call read_options(debug,chemist_notation,doOv,doKin,doNuc,doERI,doF12,doYuk,doErf,do3eInt,do4eInt)
|
call read_options(debug,chemist_notation,ExpS,doOv,doKin,doNuc,doERI,doF12,doYuk,doErf,do3eInt,do4eInt)
|
||||||
|
|
||||||
! Which integrals do you want?
|
! Which integrals do you want?
|
||||||
|
|
||||||
@ -61,26 +63,29 @@ program IntPak
|
|||||||
! Read input information
|
! Read input information
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
call ReadNAtoms(NAtoms)
|
! 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
|
||||||
|
|
||||||
allocate(ZNuc(1:NAtoms),XYZAtoms(1:NAtoms,1:3))
|
call read_molecule(nNuc,nEl,nO)
|
||||||
|
|
||||||
call ReadGeometry(NAtoms,ZNuc,XYZAtoms)
|
allocate(ZNuc(1:nNuc),rNuc(1:nNuc,1:3))
|
||||||
|
|
||||||
|
! Read geometry
|
||||||
|
|
||||||
|
call read_geometry(nNuc,ZNuc,rNuc,ENuc)
|
||||||
|
|
||||||
allocate(CenterShell(1:maxShell,1:3),TotAngMomShell(1:maxShell),KShell(1:maxShell), &
|
allocate(CenterShell(1:maxShell,1:3),TotAngMomShell(1:maxShell),KShell(1:maxShell), &
|
||||||
DShell(1:maxShell,1:maxK),ExpShell(1:maxShell,1:maxK))
|
DShell(1:maxShell,1:maxK),ExpShell(1:maxShell,1:maxK))
|
||||||
|
|
||||||
call ReadBasis(NAtoms,XYZAtoms,nShell,CenterShell, &
|
call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
|
||||||
TotAngMomShell,KShell,DShell,ExpShell)
|
|
||||||
|
|
||||||
call CalcNBasis(nShell,TotAngMomShell,NBasis)
|
|
||||||
|
|
||||||
call ReadGeminal(ExpS)
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
! Memory allocation
|
! Memory allocation
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
allocate(S(1:NBasis,1:NBasis))
|
allocate(S(1:nBas,1:nBas))
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
! Compute one-electron overlap integrals
|
! Compute one-electron overlap integrals
|
||||||
@ -90,7 +95,7 @@ program IntPak
|
|||||||
iType = 1
|
iType = 1
|
||||||
|
|
||||||
call cpu_time(start_1eInt(iType))
|
call cpu_time(start_1eInt(iType))
|
||||||
call ComputeOv(debug,NBasis,nShell, &
|
call ComputeOv(debug,nBas,nShell, &
|
||||||
CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
|
CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
|
||||||
np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType),S)
|
np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType),S)
|
||||||
call cpu_time(end_1eInt(iType))
|
call cpu_time(end_1eInt(iType))
|
||||||
@ -148,7 +153,7 @@ program IntPak
|
|||||||
call cpu_time(start_1eInt(iType))
|
call cpu_time(start_1eInt(iType))
|
||||||
call ComputeNuc(debug,nShell, &
|
call ComputeNuc(debug,nShell, &
|
||||||
CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
|
CenterShell,TotAngMomShell,KShell,DShell,ExpShell, &
|
||||||
NAtoms,ZNuc,XYZAtoms, &
|
nNuc,ZNuc,rNuc, &
|
||||||
np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType))
|
np1eInt(iType),nSigp1eInt(iType),nc1eInt(iType),nSigc1eInt(iType))
|
||||||
call cpu_time(end_1eInt(iType))
|
call cpu_time(end_1eInt(iType))
|
||||||
|
|
||||||
|
|||||||
@ -1,176 +0,0 @@
|
|||||||
subroutine ReadBasis(NAtoms,XYZAtoms,nShell,CenterShell, &
|
|
||||||
TotAngMomShell,KShell,DShell,ExpShell)
|
|
||||||
|
|
||||||
! Read basis set information
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: NAtoms
|
|
||||||
double precision,intent(in) :: XYZAtoms(NAtoms,3)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: nShAt,iAt
|
|
||||||
integer :: i,j,k
|
|
||||||
character :: shelltype
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
integer,intent(out) :: nShell
|
|
||||||
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,NAtoms
|
|
||||||
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) = XYZAtoms(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 in OBS',nShell
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
! Close file with basis set specification
|
|
||||||
|
|
||||||
close(unit=2)
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
! Auxiliary basis set information
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
! Open file with auxilairy basis specification
|
|
||||||
|
|
||||||
open(unit=3,file='input/auxbasis')
|
|
||||||
|
|
||||||
! Read basis information
|
|
||||||
|
|
||||||
write(*,'(A28)') 'Auxiliary basis set'
|
|
||||||
write(*,'(A28)') '-------------------'
|
|
||||||
|
|
||||||
do i=1,NAtoms
|
|
||||||
read(3,*) 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) = XYZAtoms(iAt,k)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! Shell type and contraction degree
|
|
||||||
|
|
||||||
read(3,*) 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(3,*) 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 in ABS',nShell
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
! Close file with basis set specification
|
|
||||||
|
|
||||||
close(unit=3)
|
|
||||||
|
|
||||||
end subroutine ReadBasis
|
|
||||||
@ -1,20 +0,0 @@
|
|||||||
subroutine ReadNAtoms(NAtoms)
|
|
||||||
|
|
||||||
! Read number of atoms
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
integer,intent(out) :: NAtoms
|
|
||||||
|
|
||||||
! Open file with geometry specification
|
|
||||||
open(unit=1,file='input/molecule')
|
|
||||||
|
|
||||||
! Read number of atoms
|
|
||||||
read(1,*)
|
|
||||||
read(1,*) NAtoms
|
|
||||||
|
|
||||||
! Close file with geometry specification
|
|
||||||
close(unit=1)
|
|
||||||
|
|
||||||
end subroutine ReadNAtoms
|
|
||||||
@ -65,13 +65,13 @@ program MCQC
|
|||||||
! Hello World
|
! Hello World
|
||||||
|
|
||||||
write(*,*)
|
write(*,*)
|
||||||
write(*,*) '********************************'
|
write(*,*) '******************************************************************************************'
|
||||||
write(*,*) '* Quack *'
|
write(*,*) '* Quack Quack Quack *'
|
||||||
write(*,*) '* __ __ __ *'
|
write(*,*) '* __ __ __ __ __ __ __ __ __ *'
|
||||||
write(*,*) '* <(o )___ <(o )___ <(o )___ *'
|
write(*,*) '* <(o )___ <(o )___ <(o )___ <(o )___ <(o )___ <(o )___ <(o )___ <(o )___ <(o )___ *'
|
||||||
write(*,*) '* ( ._> / ( ._> / ( ._> / *'
|
write(*,*) '* ( ._> / ( ._> / ( ._> / ( ._> / ( ._> / ( ._> / ( ._> / ( ._> / ( ._> / *'
|
||||||
write(*,*) '*|----------------------------|*'
|
write(*,*) '*|--------------------------------------------------------------------------------------|*'
|
||||||
write(*,*) '********************************'
|
write(*,*) '******************************************************************************************'
|
||||||
write(*,*)
|
write(*,*)
|
||||||
|
|
||||||
! Which calculations do you want to do?
|
! Which calculations do you want to do?
|
||||||
|
|||||||
@ -1,29 +0,0 @@
|
|||||||
function NormCoeff(alpha,a)
|
|
||||||
|
|
||||||
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
|
|
||||||
@ -16,7 +16,6 @@ subroutine qsGW_PT(nBas,nC,nO,nV,nR,nS,e0,SigCm)
|
|||||||
integer :: x,y,z,t
|
integer :: x,y,z,t
|
||||||
double precision :: eps
|
double precision :: eps
|
||||||
double precision,allocatable :: e1(:),e2(:),e3(:),e4(:)
|
double precision,allocatable :: e1(:),e2(:),e3(:),e4(:)
|
||||||
double precision,parameter :: threshold = 1d-15
|
|
||||||
|
|
||||||
! Allocation
|
! Allocation
|
||||||
|
|
||||||
|
|||||||
@ -1,176 +0,0 @@
|
|||||||
subroutine read_auxiliary_basis(NAtoms,XYZAtoms,nShell,CenterShell, &
|
|
||||||
TotAngMomShell,KShell,DShell,ExpShell)
|
|
||||||
|
|
||||||
! Read auxiliary basis set information
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: NAtoms
|
|
||||||
double precision,intent(in) :: XYZAtoms(NAtoms,3)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: nShAt,iAt
|
|
||||||
integer :: i,j,k
|
|
||||||
character :: shelltype
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
integer,intent(out) :: nShell
|
|
||||||
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,NAtoms
|
|
||||||
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) = XYZAtoms(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 in OBS',nShell
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
! Close file with basis set specification
|
|
||||||
|
|
||||||
close(unit=2)
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
! Auxiliary basis set information
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
! Open file with auxilairy basis specification
|
|
||||||
|
|
||||||
open(unit=3,file='input/auxbasis')
|
|
||||||
|
|
||||||
! Read basis information
|
|
||||||
|
|
||||||
write(*,'(A28)') 'Auxiliary basis set'
|
|
||||||
write(*,'(A28)') '-------------------'
|
|
||||||
|
|
||||||
do i=1,NAtoms
|
|
||||||
read(3,*) 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) = XYZAtoms(iAt,k)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! Shell type and contraction degree
|
|
||||||
|
|
||||||
read(3,*) 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(3,*) 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 in ABS',nShell
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
! Close file with basis set specification
|
|
||||||
|
|
||||||
close(unit=3)
|
|
||||||
|
|
||||||
end subroutine read_auxiliary_basis
|
|
||||||
@ -1,128 +0,0 @@
|
|||||||
subroutine read_basis(nAt,rAt,nBas,nC,nO,nV,nR,nS, &
|
|
||||||
nShell,atot,CenterShell,TotAngMomShell,KShell,DShell,ExpShell)
|
|
||||||
|
|
||||||
! Read basis set information
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: nAt,nC,nO,nR,atot(maxShell)
|
|
||||||
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,nS
|
|
||||||
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 in OBS',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 + (atot(iShell)*atot(iShell) + 3*atot(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
|
|
||||||
|
|
||||||
if(nR > nV) then
|
|
||||||
write(*,*) 'Number of Rydberg orbitals greater than number of virtual orbitals!'
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
! Number of single excitation
|
|
||||||
|
|
||||||
nS = (nO - nC)*(nV - nR)
|
|
||||||
|
|
||||||
|
|
||||||
end subroutine read_basis
|
|
||||||
@ -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
|
|
||||||
@ -13,7 +13,6 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
|
|||||||
logical :: debug
|
logical :: debug
|
||||||
integer :: mu,nu,la,si
|
integer :: mu,nu,la,si
|
||||||
double precision :: Ov,Kin,Nuc,ERI
|
double precision :: Ov,Kin,Nuc,ERI
|
||||||
double precision :: scale
|
|
||||||
|
|
||||||
! Output variables
|
! Output variables
|
||||||
|
|
||||||
@ -23,8 +22,6 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
|
|||||||
|
|
||||||
debug = .false.
|
debug = .false.
|
||||||
|
|
||||||
scale = 1d0
|
|
||||||
|
|
||||||
open(unit=8 ,file='int/Ov.dat')
|
open(unit=8 ,file='int/Ov.dat')
|
||||||
open(unit=9 ,file='int/Kin.dat')
|
open(unit=9 ,file='int/Kin.dat')
|
||||||
open(unit=10,file='int/Nuc.dat')
|
open(unit=10,file='int/Nuc.dat')
|
||||||
@ -44,7 +41,7 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
|
|||||||
T = 0d0
|
T = 0d0
|
||||||
do
|
do
|
||||||
read(9,*,end=9) mu,nu,Kin
|
read(9,*,end=9) mu,nu,Kin
|
||||||
T(mu,nu) = Kin/scale**2
|
T(mu,nu) = Kin
|
||||||
enddo
|
enddo
|
||||||
9 close(unit=9)
|
9 close(unit=9)
|
||||||
|
|
||||||
@ -65,25 +62,22 @@ subroutine read_integrals(nBas,S,T,V,Hc,G)
|
|||||||
|
|
||||||
G = 0d0
|
G = 0d0
|
||||||
do
|
do
|
||||||
read(11,*,end=11) mu,la,nu,si,ERI
|
read(11,*,end=11) mu,nu,la,si,ERI
|
||||||
! read(11,*,end=11) ERI,mu,nu,la,si
|
! (12|34)
|
||||||
|
|
||||||
ERI = ERI/scale
|
|
||||||
! <12|34>
|
|
||||||
G(mu,nu,la,si) = ERI
|
G(mu,nu,la,si) = ERI
|
||||||
! <32|14>
|
! (21|34)
|
||||||
G(la,nu,mu,si) = ERI
|
G(nu,mu,la,si) = ERI
|
||||||
! <14|32>
|
! (12|43)
|
||||||
G(mu,si,la,nu) = ERI
|
G(mu,nu,si,la) = ERI
|
||||||
! <34|12>
|
! (21|43)
|
||||||
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
|
G(nu,mu,si,la) = ERI
|
||||||
! <43|21>
|
! (34|12)
|
||||||
|
G(la,si,mu,nu) = ERI
|
||||||
|
! (43|12)
|
||||||
|
G(si,la,mu,nu) = ERI
|
||||||
|
! (34|21)
|
||||||
|
G(la,si,nu,mu) = ERI
|
||||||
|
! (43|21)
|
||||||
G(si,la,nu,mu) = ERI
|
G(si,la,nu,mu) = ERI
|
||||||
enddo
|
enddo
|
||||||
11 close(unit=11)
|
11 close(unit=11)
|
||||||
|
|||||||
@ -1,67 +0,0 @@
|
|||||||
subroutine read_molecule(nAt,nEl,nC,nO,nR)
|
|
||||||
|
|
||||||
! Read number of atoms nAt,
|
|
||||||
! number of electrons nEl,
|
|
||||||
! number of core electrons nC,
|
|
||||||
! number of Rydberg orbitals nR
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
integer,intent(out) :: nAt,nEl,nC,nO,nR
|
|
||||||
|
|
||||||
! 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,nC,nR
|
|
||||||
|
|
||||||
! Number of occupied orbitals
|
|
||||||
|
|
||||||
if(mod(nEl,2) /= 0) then
|
|
||||||
write(*,*) 'closed-shell system required!'
|
|
||||||
! stop
|
|
||||||
endif
|
|
||||||
nO = nEl/2
|
|
||||||
|
|
||||||
! Number of core orbitals
|
|
||||||
|
|
||||||
if(mod(nC,2) /= 0) then
|
|
||||||
write(*,*) 'Number of core electrons not even!'
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
nC = nC/2
|
|
||||||
|
|
||||||
if(nC > nO) then
|
|
||||||
write(*,*) 'Number of core electrons greater than number of electrons!'
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
|
|
||||||
! 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(*,*)
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of core electrons',2*nC
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,*)
|
|
||||||
write(*,*)
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of Rydberg orbitals',nR
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
! Close file with geometry specification
|
|
||||||
|
|
||||||
close(unit=1)
|
|
||||||
|
|
||||||
end subroutine read_molecule
|
|
||||||
@ -4,10 +4,10 @@ program eDFT
|
|||||||
|
|
||||||
include 'parameters.h'
|
include 'parameters.h'
|
||||||
|
|
||||||
integer :: nAt,nBas,nEl(nspin),nO(nspin),nV(nspin)
|
integer :: nNuc,nBas,nEl(nspin),nO(nspin),nV(nspin)
|
||||||
double precision :: ENuc,EKS
|
double precision :: ENuc,EKS
|
||||||
|
|
||||||
double precision,allocatable :: ZNuc(:),rAt(:,:)
|
double precision,allocatable :: ZNuc(:),rNuc(:,:)
|
||||||
|
|
||||||
integer :: nShell
|
integer :: nShell
|
||||||
integer,allocatable :: TotAngMomShell(:)
|
integer,allocatable :: TotAngMomShell(:)
|
||||||
@ -55,12 +55,12 @@ program eDFT
|
|||||||
! nBas = number of basis functions (see below)
|
! nBas = number of basis functions (see below)
|
||||||
! = nO + nV
|
! = nO + nV
|
||||||
|
|
||||||
call read_molecule(nAt,nEl,nO)
|
call read_molecule(nNuc,nEl,nO)
|
||||||
allocate(ZNuc(nAt),rAt(nAt,ncart))
|
allocate(ZNuc(nNuc),rNuc(nNuc,ncart))
|
||||||
|
|
||||||
! Read geometry
|
! Read geometry
|
||||||
|
|
||||||
call read_geometry(nAt,ZNuc,rAt,ENuc)
|
call read_geometry(nNuc,ZNuc,rNuc,ENuc)
|
||||||
|
|
||||||
allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), &
|
allocate(CenterShell(maxShell,ncart),TotAngMomShell(maxShell),KShell(maxShell), &
|
||||||
DShell(maxShell,maxK),ExpShell(maxShell,maxK))
|
DShell(maxShell,maxK),ExpShell(maxShell,maxK))
|
||||||
@ -69,7 +69,7 @@ program eDFT
|
|||||||
! Read basis set information
|
! Read basis set information
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
|
|
||||||
call read_basis(nAt,rAt,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
|
call read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
!------------------------------------------------------------------------
|
||||||
! Read one- and two-electron integrals
|
! Read one- and two-electron integrals
|
||||||
|
|||||||
@ -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
|
|
||||||
! (21|34)
|
|
||||||
G(nu,mu,la,si) = ERI
|
|
||||||
! (12|43)
|
|
||||||
G(mu,nu,si,la) = ERI
|
|
||||||
! (21|43)
|
|
||||||
G(nu,mu,si,la) = ERI
|
|
||||||
! (34|12)
|
|
||||||
G(la,si,mu,nu) = ERI
|
|
||||||
! (43|12)
|
|
||||||
G(si,la,mu,nu) = ERI
|
|
||||||
! (34|21)
|
|
||||||
G(la,si,nu,mu) = 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
|
|
||||||
@ -31,7 +31,7 @@ subroutine read_options(x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns,maxSCF,thresh,DI
|
|||||||
|
|
||||||
! Open file with method specification
|
! Open file with method specification
|
||||||
|
|
||||||
open(unit=1,file='input/options')
|
open(unit=1,file='input/dft')
|
||||||
|
|
||||||
! Default values
|
! Default values
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user