mirror of
https://github.com/pfloos/quack
synced 2024-12-23 12:56:38 +01:00
typos EOM
This commit is contained in:
parent
ca3a985d50
commit
0405a53fa8
@ -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
|
|
||||||
@ -1,170 +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 later on
|
|
||||||
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
|
|
||||||
|
|
||||||
@ -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,118 +0,0 @@
|
|||||||
subroutine read_basis(nNuc,rNuc,nBas,nO,nV,nShell,TotAngMomShell,CenterShell,KShell,DShell,ExpShell)
|
|
||||||
|
|
||||||
! Read basis set information
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: nNuc,nO
|
|
||||||
double precision,intent(in) :: rNuc(nNuc,ncart)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: nShAt,iNuc,iShell
|
|
||||||
integer :: i,j,k
|
|
||||||
character :: shelltype
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
integer,intent(out) :: nShell,nBas
|
|
||||||
double precision,intent(out) :: CenterShell(maxShell,ncart)
|
|
||||||
integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell)
|
|
||||||
double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
|
|
||||||
integer,intent(out) :: nV
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
! 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,nNuc
|
|
||||||
read(2,*) iNuc,nShAt
|
|
||||||
write(*,'(A28,1X,I16)') 'Atom n. ',iNuc
|
|
||||||
write(*,'(A28,1X,I16)') 'number of shells ',nShAt
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
|
|
||||||
! Basis function centers
|
|
||||||
|
|
||||||
do j=1,nShAt
|
|
||||||
nShell = nShell + 1
|
|
||||||
do k=1,ncart
|
|
||||||
CenterShell(nShell,k) = rNuc(iNuc,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
|
|
||||||
@ -1,68 +0,0 @@
|
|||||||
subroutine read_geometry(nAt,ZNuc,rA,ENuc)
|
|
||||||
|
|
||||||
! Read molecular geometry
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Ouput variables
|
|
||||||
|
|
||||||
integer,intent(in) :: nAt
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: i,j
|
|
||||||
double precision :: RAB
|
|
||||||
character(len=2) :: El
|
|
||||||
integer,external :: element_number
|
|
||||||
|
|
||||||
! Ouput variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: ZNuc(NAt),rA(nAt,ncart),ENuc
|
|
||||||
|
|
||||||
! Open file with geometry specification
|
|
||||||
|
|
||||||
open(unit=1,file='input/molecule')
|
|
||||||
|
|
||||||
! Read geometry
|
|
||||||
|
|
||||||
read(1,*)
|
|
||||||
read(1,*)
|
|
||||||
read(1,*)
|
|
||||||
|
|
||||||
do i=1,nAt
|
|
||||||
read(1,*) El,rA(i,1),rA(i,2),rA(i,3)
|
|
||||||
ZNuc(i) = element_number(El)
|
|
||||||
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,ncart)
|
|
||||||
enddo
|
|
||||||
write(*,*)
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
end subroutine read_geometry
|
|
||||||
@ -1,39 +0,0 @@
|
|||||||
subroutine read_molecule(nAt,nEl,nO,nC,nR)
|
|
||||||
|
|
||||||
! Read number of atoms nAt and number of electrons nEl
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(out) :: nAt,nEl,nO,nC,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
|
|
||||||
|
|
||||||
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
|
|
||||||
@ -1,385 +0,0 @@
|
|||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
function KroneckerDelta(i,j) result(delta)
|
|
||||||
|
|
||||||
! Kronecker Delta
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: i,j
|
|
||||||
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
integer :: delta
|
|
||||||
|
|
||||||
if(i == j) then
|
|
||||||
delta = 1
|
|
||||||
else
|
|
||||||
delta = 0
|
|
||||||
endif
|
|
||||||
|
|
||||||
end function KroneckerDelta
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
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(7X,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
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
subroutine CalcTrAB(n,A,B,Tr)
|
|
||||||
|
|
||||||
! Calculate the trace of the square matrix A.B
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: n
|
|
||||||
double precision,intent(in) :: A(n,n),B(n,n)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: i,j
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: Tr
|
|
||||||
|
|
||||||
Tr = 0d0
|
|
||||||
do i=1,n
|
|
||||||
do j=1,n
|
|
||||||
Tr = Tr + A(i,j)*B(j,i)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine CalcTrAB
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
function EpsilonSwitch(i,j) result(delta)
|
|
||||||
|
|
||||||
! Epsilon function
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: i,j
|
|
||||||
integer :: delta
|
|
||||||
|
|
||||||
if(i <= j) then
|
|
||||||
delta = 1
|
|
||||||
else
|
|
||||||
delta = -1
|
|
||||||
endif
|
|
||||||
|
|
||||||
end function EpsilonSwitch
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
function KappaCross(i,j,k) result(kappa)
|
|
||||||
|
|
||||||
! kappa(i,j,k) = eps(i,j) delta(i,k) - eps(k,i) delta(i,j)
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: i,j,k
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: EpsilonSwitch,KroneckerDelta
|
|
||||||
double precision :: kappa
|
|
||||||
|
|
||||||
kappa = dble(EpsilonSwitch(i,j)*KroneckerDelta(i,k) - EpsilonSwitch(k,i)*KroneckerDelta(i,j))
|
|
||||||
|
|
||||||
end function KappaCross
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
subroutine CalcInv3(a,det)
|
|
||||||
|
|
||||||
! Calculate the inverse and the determinant of a 3x3 matrix
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
double precision,intent(inout) :: a(3,3)
|
|
||||||
double precision, intent(inout) :: det
|
|
||||||
double precision :: b(3,3)
|
|
||||||
integer :: i,j
|
|
||||||
|
|
||||||
det = a(1,1)*(a(2,2)*a(3,3)-a(2,3)*a(3,2)) &
|
|
||||||
- a(1,2)*(a(2,1)*a(3,3)-a(2,3)*a(3,1)) &
|
|
||||||
+ a(1,3)*(a(2,1)*a(3,2)-a(2,2)*a(3,1))
|
|
||||||
|
|
||||||
do i=1,3
|
|
||||||
b(i,1) = a(i,1)
|
|
||||||
b(i,2) = a(i,2)
|
|
||||||
b(i,3) = a(i,3)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
a(1,1) = b(2,2)*b(3,3) - b(2,3)*b(3,2)
|
|
||||||
a(2,1) = b(2,3)*b(3,1) - b(2,1)*b(3,3)
|
|
||||||
a(3,1) = b(2,1)*b(3,2) - b(2,2)*b(3,1)
|
|
||||||
|
|
||||||
a(1,2) = b(1,3)*b(3,2) - b(1,2)*b(3,3)
|
|
||||||
a(2,2) = b(1,1)*b(3,3) - b(1,3)*b(3,1)
|
|
||||||
a(3,2) = b(1,2)*b(3,1) - b(1,1)*b(3,2)
|
|
||||||
|
|
||||||
a(1,3) = b(1,2)*b(2,3) - b(1,3)*b(2,2)
|
|
||||||
a(2,3) = b(1,3)*b(2,1) - b(1,1)*b(2,3)
|
|
||||||
a(3,3) = b(1,1)*b(2,2) - b(1,2)*b(2,1)
|
|
||||||
|
|
||||||
do i=1,3
|
|
||||||
do j=1,3
|
|
||||||
a(i,j) = a(i,j)/det
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine CalcInv3
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
subroutine CalcInv4(a,det)
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
double precision,intent(inout) :: a(4,4)
|
|
||||||
double precision,intent(inout) :: det
|
|
||||||
double precision :: b(4,4)
|
|
||||||
integer :: i,j
|
|
||||||
|
|
||||||
det = a(1,1)*(a(2,2)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
|
|
||||||
-a(2,3)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
|
|
||||||
+a(2,4)*(a(3,2)*a(4,3)-a(3,3)*a(4,2))) &
|
|
||||||
- a(1,2)*(a(2,1)*(a(3,3)*a(4,4)-a(3,4)*a(4,3)) &
|
|
||||||
-a(2,3)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) &
|
|
||||||
+a(2,4)*(a(3,1)*a(4,3)-a(3,3)*a(4,1))) &
|
|
||||||
+ a(1,3)*(a(2,1)*(a(3,2)*a(4,4)-a(3,4)*a(4,2)) &
|
|
||||||
-a(2,2)*(a(3,1)*a(4,4)-a(3,4)*a(4,1)) &
|
|
||||||
+a(2,4)*(a(3,1)*a(4,2)-a(3,2)*a(4,1))) &
|
|
||||||
- a(1,4)*(a(2,1)*(a(3,2)*a(4,3)-a(3,3)*a(4,2)) &
|
|
||||||
-a(2,2)*(a(3,1)*a(4,3)-a(3,3)*a(4,1)) &
|
|
||||||
+a(2,3)*(a(3,1)*a(4,2)-a(3,2)*a(4,1)))
|
|
||||||
do i=1,4
|
|
||||||
b(1,i) = a(1,i)
|
|
||||||
b(2,i) = a(2,i)
|
|
||||||
b(3,i) = a(3,i)
|
|
||||||
b(4,i) = a(4,i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
a(1,1) = b(2,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(2,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(2,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
|
|
||||||
a(2,1) = -b(2,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(2,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(2,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
|
|
||||||
a(3,1) = b(2,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(2,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(2,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
|
|
||||||
a(4,1) = -b(2,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))+b(2,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))-b(2,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
|
|
||||||
|
|
||||||
a(1,2) = -b(1,2)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))+b(1,3)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))-b(1,4)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))
|
|
||||||
a(2,2) = b(1,1)*(b(3,3)*b(4,4)-b(3,4)*b(4,3))-b(1,3)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))+b(1,4)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))
|
|
||||||
a(3,2) = -b(1,1)*(b(3,2)*b(4,4)-b(3,4)*b(4,2))+b(1,2)*(b(3,1)*b(4,4)-b(3,4)*b(4,1))-b(1,4)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
|
|
||||||
a(4,2) = b(1,1)*(b(3,2)*b(4,3)-b(3,3)*b(4,2))-b(1,2)*(b(3,1)*b(4,3)-b(3,3)*b(4,1))+b(1,3)*(b(3,1)*b(4,2)-b(3,2)*b(4,1))
|
|
||||||
|
|
||||||
a(1,3) = b(1,2)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))-b(1,3)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))+b(1,4)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))
|
|
||||||
a(2,3) = -b(1,1)*(b(2,3)*b(4,4)-b(2,4)*b(4,3))+b(1,3)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))-b(1,4)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))
|
|
||||||
a(3,3) = b(1,1)*(b(2,2)*b(4,4)-b(2,4)*b(4,2))-b(1,2)*(b(2,1)*b(4,4)-b(2,4)*b(4,1))+b(1,4)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))
|
|
||||||
a(4,3) = -b(1,1)*(b(2,2)*b(4,3)-b(2,3)*b(4,2))+b(1,2)*(b(2,1)*b(4,3)-b(2,3)*b(4,1))-b(1,3)*(b(2,1)*b(4,2)-b(2,2)*b(4,1))
|
|
||||||
|
|
||||||
a(1,4) = -b(1,2)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))+b(1,3)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))-b(1,4)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))
|
|
||||||
a(2,4) = b(1,1)*(b(2,3)*b(3,4)-b(2,4)*b(3,3))-b(1,3)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))+b(1,4)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))
|
|
||||||
a(3,4) = -b(1,1)*(b(2,2)*b(3,4)-b(2,4)*b(3,2))+b(1,2)*(b(2,1)*b(3,4)-b(2,4)*b(3,1))-b(1,4)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
|
|
||||||
a(4,4) = b(1,1)*(b(2,2)*b(3,3)-b(2,3)*b(3,2))-b(1,2)*(b(2,1)*b(3,3)-b(2,3)*b(3,1))+b(1,3)*(b(2,1)*b(3,2)-b(2,2)*b(3,1))
|
|
||||||
|
|
||||||
do i=1,4
|
|
||||||
do j=1,4
|
|
||||||
a(i,j) = a(i,j)/det
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine CalcInv4
|
|
||||||
|
|
||||||
|
|
||||||
!double precision function binom(i,j)
|
|
||||||
! implicit none
|
|
||||||
! integer,intent(in) :: i,j
|
|
||||||
! double precision :: logfact
|
|
||||||
! integer, save :: ifirst
|
|
||||||
! double precision, save :: memo(0:15,0:15)
|
|
||||||
! integer :: k,l
|
|
||||||
! if (ifirst == 0) then
|
|
||||||
! ifirst = 1
|
|
||||||
! do k=0,15
|
|
||||||
! do l=0,15
|
|
||||||
! memo(k,l) = dexp( logfact(k)-logfact(l)-logfact(k-l) )
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! endif
|
|
||||||
! if ( (i<=15).and.(j<=15) ) then
|
|
||||||
! binom = memo(i,j)
|
|
||||||
! else
|
|
||||||
! binom = dexp( logfact(i)-logfact(j)-logfact(i-j) )
|
|
||||||
! endif
|
|
||||||
!end
|
|
||||||
!
|
|
||||||
!double precision function fact(n)
|
|
||||||
! implicit none
|
|
||||||
! integer :: n
|
|
||||||
! double precision, save :: memo(1:100)
|
|
||||||
! integer, save :: memomax = 1
|
|
||||||
!
|
|
||||||
! if (n<=memomax) then
|
|
||||||
! if (n<2) then
|
|
||||||
! fact = 1.d0
|
|
||||||
! else
|
|
||||||
! fact = memo(n)
|
|
||||||
! endif
|
|
||||||
! return
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
! integer :: i
|
|
||||||
! memo(1) = 1.d0
|
|
||||||
! do i=memomax+1,min(n,100)
|
|
||||||
! memo(i) = memo(i-1)*dble(i)
|
|
||||||
! enddo
|
|
||||||
! memomax = min(n,100)
|
|
||||||
! double precision :: logfact
|
|
||||||
! fact = dexp(logfact(n))
|
|
||||||
!end function
|
|
||||||
!
|
|
||||||
!double precision function logfact(n)
|
|
||||||
! implicit none
|
|
||||||
! integer :: n
|
|
||||||
! double precision, save :: memo(1:100)
|
|
||||||
! integer, save :: memomax = 1
|
|
||||||
!
|
|
||||||
! if (n<=memomax) then
|
|
||||||
! if (n<2) then
|
|
||||||
! logfact = 0.d0
|
|
||||||
! else
|
|
||||||
! logfact = memo(n)
|
|
||||||
! endif
|
|
||||||
! return
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
! integer :: i
|
|
||||||
! memo(1) = 0.d0
|
|
||||||
! do i=memomax+1,min(n,100)
|
|
||||||
! memo(i) = memo(i-1)+dlog(dble(i))
|
|
||||||
! enddo
|
|
||||||
! memomax = min(n,100)
|
|
||||||
! logfact = memo(memomax)
|
|
||||||
! do i=101,n
|
|
||||||
! logfact += dlog(dble(i))
|
|
||||||
! enddo
|
|
||||||
!end function
|
|
||||||
!
|
|
||||||
!double precision function dble_fact(n)
|
|
||||||
! implicit none
|
|
||||||
! integer :: n
|
|
||||||
! double precision :: dble_fact_even, dble_fact_odd
|
|
||||||
!
|
|
||||||
! dble_fact = 1.d0
|
|
||||||
!
|
|
||||||
! if(n.lt.0) return
|
|
||||||
!
|
|
||||||
! if(iand(n,1).eq.0)then
|
|
||||||
! dble_fact = dble_fact_even(n)
|
|
||||||
! else
|
|
||||||
! dble_fact= dble_fact_odd(n)
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
!end function
|
|
||||||
!
|
|
||||||
!double precision function dble_fact_even(n) result(fact2)
|
|
||||||
! implicit none
|
|
||||||
! integer :: n,k
|
|
||||||
! double precision, save :: memo(0:100)
|
|
||||||
! integer, save :: memomax = 0
|
|
||||||
! double precision :: prod
|
|
||||||
!
|
|
||||||
!
|
|
||||||
! if (n <= memomax) then
|
|
||||||
! if (n < 2) then
|
|
||||||
! fact2 = 1.d0
|
|
||||||
! else
|
|
||||||
! fact2 = memo(n)
|
|
||||||
! endif
|
|
||||||
! return
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
! integer :: i
|
|
||||||
! memo(0)=1.d0
|
|
||||||
! memo(1)=1.d0
|
|
||||||
! do i=memomax+2,min(n,100),2
|
|
||||||
! memo(i) = memo(i-2)* dble(i)
|
|
||||||
! enddo
|
|
||||||
! memomax = min(n,100)
|
|
||||||
! fact2 = memo(memomax)
|
|
||||||
!
|
|
||||||
! if (n > 100) then
|
|
||||||
! double precision :: dble_logfact
|
|
||||||
! fact2 = dexp(dble_logfact(n))
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
!end function
|
|
||||||
!
|
|
||||||
!double precision function dble_fact_odd(n) result(fact2)
|
|
||||||
! implicit none
|
|
||||||
! integer :: n
|
|
||||||
! double precision, save :: memo(1:100)
|
|
||||||
! integer, save :: memomax = 1
|
|
||||||
!
|
|
||||||
! if (n<=memomax) then
|
|
||||||
! if (n<3) then
|
|
||||||
! fact2 = 1.d0
|
|
||||||
! else
|
|
||||||
! fact2 = memo(n)
|
|
||||||
! endif
|
|
||||||
! return
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
! integer :: i
|
|
||||||
! memo(1) = 1.d0
|
|
||||||
! do i=memomax+2,min(n,99),2
|
|
||||||
! memo(i) = memo(i-2)* dble(i)
|
|
||||||
! enddo
|
|
||||||
! memomax = min(n,99)
|
|
||||||
! fact2 = memo(memomax)
|
|
||||||
!
|
|
||||||
! if (n > 99) then
|
|
||||||
! double precision :: dble_logfact
|
|
||||||
! fact2 = dexp(dble_logfact(n))
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
!end function
|
|
||||||
|
|
||||||
@ -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
|
|
||||||
@ -1,170 +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 later on
|
|
||||||
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
|
|
||||||
|
|
||||||
@ -1,118 +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,ncart)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: nShAt,iAt,iShell
|
|
||||||
integer :: i,j,k
|
|
||||||
character :: shelltype
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
integer,intent(out) :: nShell,nBas
|
|
||||||
double precision,intent(out) :: CenterShell(maxShell,ncart)
|
|
||||||
integer,intent(out) :: TotAngMomShell(maxShell),KShell(maxShell)
|
|
||||||
double precision,intent(out) :: DShell(maxShell,maxK),ExpShell(maxShell,maxK)
|
|
||||||
integer,intent(out) :: nV
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
! 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,ncart
|
|
||||||
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
|
|
||||||
@ -1,68 +0,0 @@
|
|||||||
subroutine read_geometry(nNuc,ZNuc,rNuc,ENuc)
|
|
||||||
|
|
||||||
! Read molecular geometry
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Ouput variables
|
|
||||||
|
|
||||||
integer,intent(in) :: nNuc
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: i,j
|
|
||||||
double precision :: RAB
|
|
||||||
character(len=2) :: El
|
|
||||||
integer,external :: element_number
|
|
||||||
|
|
||||||
! Ouput variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: ZNuc(nNuc),rNuc(nNuc,ncart),ENuc
|
|
||||||
|
|
||||||
! Open file with geometry specification
|
|
||||||
|
|
||||||
open(unit=1,file='input/molecule')
|
|
||||||
|
|
||||||
! Read geometry
|
|
||||||
|
|
||||||
read(1,*)
|
|
||||||
read(1,*)
|
|
||||||
read(1,*)
|
|
||||||
|
|
||||||
do i=1,nNuc
|
|
||||||
read(1,*) El,rNuc(i,1),rNuc(i,2),rNuc(i,3)
|
|
||||||
ZNuc(i) = element_number(El)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
! Compute nuclear repulsion energy
|
|
||||||
|
|
||||||
ENuc = 0
|
|
||||||
|
|
||||||
do i=1,nNuc-1
|
|
||||||
do j=i+1,nNuc
|
|
||||||
RAB = (rNuc(i,1)-rNuc(j,1))**2 + (rNuc(i,2)-rNuc(j,2))**2 + (rNuc(i,3)-rNuc(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,nNuc
|
|
||||||
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:',(rNuc(i,j),j=1,ncart)
|
|
||||||
enddo
|
|
||||||
write(*,*)
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,'(A28,1X,F16.10)') 'Nuclear repulsion energy = ',ENuc
|
|
||||||
write(*,'(A28)') '------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
end subroutine read_geometry
|
|
||||||
@ -1,119 +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
|
|
||||||
double precision :: scale
|
|
||||||
|
|
||||||
! 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.
|
|
||||||
|
|
||||||
scale = 1d0
|
|
||||||
|
|
||||||
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/scale**2
|
|
||||||
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
|
|
||||||
|
|
||||||
ERI = ERI/scale
|
|
||||||
! <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
|
|
||||||
@ -1,64 +0,0 @@
|
|||||||
subroutine read_molecule(nNuc,nEl,nO,nC,nR)
|
|
||||||
|
|
||||||
! Read number of atoms and number of electrons
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
include 'parameters.h'
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(out) :: nNuc
|
|
||||||
integer,intent(out) :: nEl(nspin)
|
|
||||||
integer,intent(out) :: nO(nspin)
|
|
||||||
integer,intent(out) :: nC(nspin)
|
|
||||||
integer,intent(out) :: nR(nspin)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: nCore
|
|
||||||
integer :: nRyd
|
|
||||||
|
|
||||||
! Open file with geometry specification
|
|
||||||
|
|
||||||
open(unit=1,file='input/molecule')
|
|
||||||
|
|
||||||
! Read number of atoms and number of electrons
|
|
||||||
|
|
||||||
read(1,*)
|
|
||||||
read(1,*) nNuc,nEl(1),nEl(2),nCore,nRyd
|
|
||||||
|
|
||||||
if(mod(nCore,2) /= 0 .or. mod(nRyd,2) /= 0) then
|
|
||||||
|
|
||||||
print*, 'The number of core and Rydberg electrons must be even!'
|
|
||||||
stop
|
|
||||||
|
|
||||||
end if
|
|
||||||
|
|
||||||
nO(:) = nEl(:)
|
|
||||||
nC(:) = nCore/2
|
|
||||||
nR(:) = nRyd/2
|
|
||||||
|
|
||||||
! Print results
|
|
||||||
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of atoms',nNuc
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,*)
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of spin-up electrons',nEl(1)
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of spin-down electrons',nEl(2)
|
|
||||||
write(*,'(A28,1X,I16)') ' Total number of electrons',sum(nEl(:))
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,*)
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of core electrons',sum(nC(:))
|
|
||||||
write(*,'(A28,1X,I16)') 'Number of Rydberg electrons',sum(nR(:))
|
|
||||||
write(*,'(A28)') '----------------------'
|
|
||||||
write(*,*)
|
|
||||||
|
|
||||||
! Close file with geometry specification
|
|
||||||
|
|
||||||
close(unit=1)
|
|
||||||
|
|
||||||
end subroutine read_molecule
|
|
||||||
@ -1,339 +0,0 @@
|
|||||||
!------------------------------------------------------------------------
|
|
||||||
function Kronecker_delta(i,j) result(delta)
|
|
||||||
|
|
||||||
! Kronecker Delta
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: i,j
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision :: delta
|
|
||||||
|
|
||||||
if(i == j) then
|
|
||||||
delta = 1d0
|
|
||||||
else
|
|
||||||
delta = 0d0
|
|
||||||
endif
|
|
||||||
|
|
||||||
end function Kronecker_delta
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
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
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
subroutine trace_vector(n,v,Tr)
|
|
||||||
|
|
||||||
! Calculate the trace of the vector v of length n
|
|
||||||
!!! Please use the intrinsic fortran sum() !!!
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: n
|
|
||||||
double precision,intent(in) :: v(n)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: i
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: Tr
|
|
||||||
|
|
||||||
Tr = 0d0
|
|
||||||
do i=1,n
|
|
||||||
Tr = Tr + v(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine trace_vector
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
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 compute_error(nData,Mean,Var,Error)
|
|
||||||
|
|
||||||
! Calculate the statistical error
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
double precision,intent(in) :: nData,Mean(3)
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: Error(3)
|
|
||||||
double precision,intent(inout):: Var(3)
|
|
||||||
|
|
||||||
Error = sqrt((Var-Mean**2/nData)/nData/(nData-1d0))
|
|
||||||
|
|
||||||
end subroutine compute_error
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
subroutine identity_matrix(N,A)
|
|
||||||
|
|
||||||
! Set the matrix A to the identity matrix
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: N
|
|
||||||
|
|
||||||
! Local viaruabkes
|
|
||||||
|
|
||||||
integer :: i
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: A(N,N)
|
|
||||||
|
|
||||||
A = 0d0
|
|
||||||
|
|
||||||
do i=1,N
|
|
||||||
A(i,i) = 1d0
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine identity_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,k
|
|
||||||
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,k
|
|
||||||
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
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
subroutine print_warning(message)
|
|
||||||
|
|
||||||
! Print warning
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
character(len=*),intent(in) :: message
|
|
||||||
|
|
||||||
write(*,*) message
|
|
||||||
|
|
||||||
end subroutine print_warning
|
|
||||||
|
|
||||||
|
|
||||||
@ -1,207 +0,0 @@
|
|||||||
!subroutine eigenvalues_non_symmetric_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(eRe(N),eIm(N),work(3*N))
|
|
||||||
! lwork = size(work)
|
|
||||||
!
|
|
||||||
! call DGEEV('N','N',N,A,N, eRe, eIm, 0d0,1, VR,LDVR, WORK, LWORK, INFO )
|
|
||||||
!
|
|
||||||
! if(info /= 0) then
|
|
||||||
! print*,'Problem in diagonalize_matrix (dseev)!!'
|
|
||||||
! stop
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
!end subroutine eigenvalues_non_symmetric_matrix
|
|
||||||
|
|
||||||
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)!!'
|
|
||||||
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,rcond)
|
|
||||||
|
|
||||||
! 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),rcond
|
|
||||||
double precision,intent(out) :: x(N)
|
|
||||||
|
|
||||||
integer :: info,lwork
|
|
||||||
double precision :: ferr,berr
|
|
||||||
integer,allocatable :: ipiv(:),iwork(:)
|
|
||||||
double precision,allocatable :: AF(:,:),work(:)
|
|
||||||
|
|
||||||
lwork = 3*N
|
|
||||||
allocate(AF(N,N),ipiv(N),work(lwork),iwork(N))
|
|
||||||
|
|
||||||
call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info)
|
|
||||||
|
|
||||||
! if (info /= 0) then
|
|
||||||
|
|
||||||
! print *, info
|
|
||||||
! stop 'error in linear_solve (dsysvx)!!'
|
|
||||||
|
|
||||||
! endif
|
|
||||||
|
|
||||||
end subroutine linear_solve
|
|
||||||
|
|
||||||
subroutine easy_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 easy_linear_solve
|
|
||||||
|
|
||||||
@ -1,400 +0,0 @@
|
|||||||
!------------------------------------------------------------------------
|
|
||||||
function Kronecker_delta(i,j) result(delta)
|
|
||||||
|
|
||||||
! Kronecker Delta
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: i,j
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision :: delta
|
|
||||||
|
|
||||||
if(i == j) then
|
|
||||||
delta = 1d0
|
|
||||||
else
|
|
||||||
delta = 0d0
|
|
||||||
endif
|
|
||||||
|
|
||||||
end function Kronecker_delta
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
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
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
subroutine trace_vector(n,v,Tr)
|
|
||||||
|
|
||||||
! Calculate the trace of the vector v of length n
|
|
||||||
!!! Please use the intrinsic fortran sum() !!!
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: n
|
|
||||||
double precision,intent(in) :: v(n)
|
|
||||||
|
|
||||||
! Local variables
|
|
||||||
|
|
||||||
integer :: i
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: Tr
|
|
||||||
|
|
||||||
Tr = 0d0
|
|
||||||
do i=1,n
|
|
||||||
Tr = Tr + v(i)
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine trace_vector
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
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 compute_error(nData,Mean,Var,Error)
|
|
||||||
|
|
||||||
! Calculate the statistical error
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
double precision,intent(in) :: nData,Mean(3)
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: Error(3)
|
|
||||||
double precision,intent(inout):: Var(3)
|
|
||||||
|
|
||||||
Error = sqrt((Var-Mean**2/nData)/nData/(nData-1d0))
|
|
||||||
|
|
||||||
end subroutine compute_error
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
subroutine identity_matrix(N,A)
|
|
||||||
|
|
||||||
! Set the matrix A to the identity matrix
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
! Input variables
|
|
||||||
|
|
||||||
integer,intent(in) :: N
|
|
||||||
|
|
||||||
! Local viaruabkes
|
|
||||||
|
|
||||||
integer :: i
|
|
||||||
|
|
||||||
! Output variables
|
|
||||||
|
|
||||||
double precision,intent(out) :: A(N,N)
|
|
||||||
|
|
||||||
A = 0d0
|
|
||||||
|
|
||||||
do i=1,N
|
|
||||||
A(i,i) = 1d0
|
|
||||||
enddo
|
|
||||||
|
|
||||||
end subroutine identity_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,k
|
|
||||||
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,k
|
|
||||||
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
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
subroutine print_warning(message)
|
|
||||||
|
|
||||||
! Print warning
|
|
||||||
|
|
||||||
implicit none
|
|
||||||
|
|
||||||
character(len=*),intent(in) :: message
|
|
||||||
|
|
||||||
write(*,*) message
|
|
||||||
|
|
||||||
end subroutine print_warning
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
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
|
|
||||||
|
|
||||||
!------------------------------------------------------------------------
|
|
||||||
|
|
||||||
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
|
|
||||||
@ -1,207 +0,0 @@
|
|||||||
!subroutine eigenvalues_non_symmetric_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(eRe(N),eIm(N),work(3*N))
|
|
||||||
! lwork = size(work)
|
|
||||||
!
|
|
||||||
! call DGEEV('N','N',N,A,N, eRe, eIm, 0d0,1, VR,LDVR, WORK, LWORK, INFO )
|
|
||||||
!
|
|
||||||
! if(info /= 0) then
|
|
||||||
! print*,'Problem in diagonalize_matrix (dseev)!!'
|
|
||||||
! stop
|
|
||||||
! endif
|
|
||||||
!
|
|
||||||
!end subroutine eigenvalues_non_symmetric_matrix
|
|
||||||
|
|
||||||
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)!!'
|
|
||||||
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,rcond)
|
|
||||||
|
|
||||||
! 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),rcond
|
|
||||||
double precision,intent(out) :: x(N)
|
|
||||||
|
|
||||||
integer :: info,lwork
|
|
||||||
double precision :: ferr,berr
|
|
||||||
integer,allocatable :: ipiv(:),iwork(:)
|
|
||||||
double precision,allocatable :: AF(:,:),work(:)
|
|
||||||
|
|
||||||
lwork = 3*N
|
|
||||||
allocate(AF(N,N),ipiv(N),work(lwork),iwork(N))
|
|
||||||
|
|
||||||
call dsysvx('N','U',N,1,A,N,AF,N,ipiv,b,N,x,N,rcond,ferr,berr,work,lwork,iwork,info)
|
|
||||||
|
|
||||||
! if (info /= 0) then
|
|
||||||
|
|
||||||
! print *, info
|
|
||||||
! stop 'error in linear_solve (dsysvx)!!'
|
|
||||||
|
|
||||||
! endif
|
|
||||||
|
|
||||||
end subroutine linear_solve
|
|
||||||
|
|
||||||
subroutine easy_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 easy_linear_solve
|
|
||||||
|
|
||||||
Loading…
Reference in New Issue
Block a user