mirror of
https://github.com/pfloos/quack
synced 2024-12-22 20:34:46 +01:00
Merge branch 'master' of https://github.com/pfloos/QuAcK
This commit is contained in:
Abdallah Ammar
commit
0deffa8638
@ -1,85 +0,0 @@
|
||||
subroutine AOtoMO_oooa(nBas,nO,nA,cO,cA,O,ooOoa)
|
||||
|
||||
! AO to MO transformation of two-electron integrals for the block oooa
|
||||
! Semi-direct O(N^5) algorithm
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas,nO,nA
|
||||
double precision,intent(in) :: cO(nBas,nO),cA(nBas,nA),O(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:)
|
||||
integer :: mu,nu,la,si,i,j,k,x
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: ooOoa(nO,nO,nO,nA)
|
||||
|
||||
! Memory allocation
|
||||
allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas))
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A42)') '----------------------------------------'
|
||||
write(*,'(A42)') ' AO to MO transformation for oooa block '
|
||||
write(*,'(A42)') '----------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
scr1 = 0d0
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do si=1,nBas
|
||||
do x=1,nA
|
||||
scr1(mu,nu,la,x) = scr1(mu,nu,la,x) + O(mu,nu,la,si)*cA(si,x)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
scr2 = 0d0
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do i=1,nO
|
||||
do x=1,nA
|
||||
scr2(i,nu,la,x) = scr2(i,nu,la,x) + cO(mu,i)*scr1(mu,nu,la,x)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
scr1 = 0d0
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do i=1,nO
|
||||
do k=1,nO
|
||||
do x=1,nA
|
||||
scr1(i,nu,k,x) = scr1(i,nu,k,x) + scr2(i,nu,la,x)*cO(la,k)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
ooOoa = 0d0
|
||||
do nu=1,nBas
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do k=1,nO
|
||||
do x=1,nA
|
||||
ooOoa(i,j,k,x) = ooOoa(i,j,k,x) + cO(nu,j)*scr1(i,nu,k,x)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
deallocate(scr1,scr2)
|
||||
|
||||
end subroutine
|
||||
@ -1,85 +0,0 @@
|
||||
subroutine AOtoMO_oooo(nBas,nO,cO,O,ooOoo)
|
||||
|
||||
! AO to MO transformation of two-electron integrals for the block oooo
|
||||
! Semi-direct O(N^5) algorithm
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas,nO
|
||||
double precision,intent(in) :: cO(nBas,nO),O(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:)
|
||||
integer :: mu,nu,la,si,i,j,k,l
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: ooOoo(nO,nO,nO,nO)
|
||||
|
||||
! Memory allocation
|
||||
allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas))
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A42)') '----------------------------------------'
|
||||
write(*,'(A42)') ' AO to MO transformation for oooo block '
|
||||
write(*,'(A42)') '----------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
scr1 = 0d0
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do si=1,nBas
|
||||
do l=1,nO
|
||||
scr1(mu,nu,la,l) = scr1(mu,nu,la,l) + O(mu,nu,la,si)*cO(si,l)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
scr2 = 0d0
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do i=1,nO
|
||||
do l=1,nO
|
||||
scr2(i,nu,la,l) = scr2(i,nu,la,l) + cO(mu,i)*scr1(mu,nu,la,l)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
scr1 = 0d0
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do i=1,nO
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
scr1(i,nu,k,l) = scr1(i,nu,k,l) + scr2(i,nu,la,l)*cO(la,k)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
ooOoo = 0d0
|
||||
do nu=1,nBas
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do k=1,nO
|
||||
do l=1,nO
|
||||
ooOoo(i,j,k,l) = ooOoo(i,j,k,l) + cO(nu,j)*scr1(i,nu,k,l)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
deallocate(scr1,scr2)
|
||||
|
||||
end subroutine
|
||||
@ -1,77 +0,0 @@
|
||||
subroutine AOtoMO_oovv(nBas,nO,nV,cO,cV,O,ooOvv)
|
||||
|
||||
! AO to MO transformation of two-electron integrals for the block oovv
|
||||
! Semi-direct O(N^5) algorithm
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas,nO,nV
|
||||
double precision,intent(in) :: cO(nBas,nO),cV(nBas,nV),O(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision,allocatable :: scr1(:,:,:,:),scr2(:,:,:,:)
|
||||
integer :: mu,nu,la,si,i,j,a,b
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: ooOvv(nO,nO,nV,nV)
|
||||
|
||||
! Memory allocation
|
||||
allocate(scr1(nBas,nBas,nBas,nBas),scr2(nBas,nBas,nBas,nBas))
|
||||
|
||||
scr1 = 0d0
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do si=1,nBas
|
||||
do b=1,nV
|
||||
scr1(mu,nu,la,b) = scr1(mu,nu,la,b) + O(mu,nu,la,si)*cV(si,b)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
scr2 = 0d0
|
||||
do mu=1,nBas
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do i=1,nO
|
||||
do b=1,nV
|
||||
scr2(i,nu,la,b) = scr2(i,nu,la,b) + cO(mu,i)*scr1(mu,nu,la,b)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
scr1 = 0d0
|
||||
do nu=1,nBas
|
||||
do la=1,nBas
|
||||
do i=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
scr1(i,nu,a,b) = scr1(i,nu,a,b) + scr2(i,nu,la,b)*cV(la,a)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
ooOvv = 0d0
|
||||
do nu=1,nBas
|
||||
do i=1,nO
|
||||
do j=1,nO
|
||||
do a=1,nV
|
||||
do b=1,nV
|
||||
ooOvv(i,j,a,b) = ooOvv(i,j,a,b) + cO(nu,j)*scr1(i,nu,a,b)
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
end subroutine
|
||||
@ -1,26 +0,0 @@
|
||||
subroutine Hartree_matrix_MO_basis(nBas,c,P,G,H)
|
||||
|
||||
! Compute Hartree matrix in the MO basis
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas)
|
||||
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: H(nBas,nBas)
|
||||
|
||||
! Compute Hartree matrix in the AO basis
|
||||
|
||||
call Hartree_matrix_AO_basis(nBas,P,G,H)
|
||||
|
||||
! Transform Hartree matrix in the MO basis
|
||||
|
||||
H = matmul(transpose(c),matmul(H,c))
|
||||
|
||||
end subroutine
|
||||
@ -1,26 +0,0 @@
|
||||
subroutine exchange_matrix_MO_basis(nBas,c,P,G,K)
|
||||
|
||||
! Compute exchange matrix in the MO basis
|
||||
|
||||
implicit none
|
||||
include 'parameters.h'
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBas
|
||||
double precision,intent(in) :: c(nBas,nBas),P(nBas,nBas)
|
||||
double precision,intent(in) :: G(nBas,nBas,nBas,nBas)
|
||||
|
||||
! Output variables
|
||||
|
||||
double precision,intent(out) :: K(nBas,nBas)
|
||||
|
||||
! Compute Hartree Hamiltonian in the AO basis
|
||||
|
||||
call exchange_matrix_AO_basis(nBas,P,G,K)
|
||||
|
||||
! Transform Coulomb matrix in the MO basis
|
||||
|
||||
K = matmul(transpose(c),matmul(K,c))
|
||||
|
||||
end subroutine exchange_matrix_MO_basis
|
||||
@ -31,7 +31,7 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
! Output variables
|
||||
|
||||
!------------------------------------------------------------------------
|
||||
! Compute MP3 energy
|
||||
! Compute MP2 energy
|
||||
!------------------------------------------------------------------------
|
||||
|
||||
if(doMP2) then
|
||||
@ -53,11 +53,12 @@ subroutine RMP(dotest,doMP2,doMP3,regularize,nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
if(doMP3) then
|
||||
|
||||
call wall_time(start_MP)
|
||||
call RMP3(nOrb,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF)
|
||||
write(*,*) 'Restricted MP3 NYI... Sorry'
|
||||
write(*,*)
|
||||
call wall_time(end_MP)
|
||||
|
||||
t_MP = end_MP - start_MP
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP2 = ',t_MP,' seconds'
|
||||
write(*,'(A65,1X,F9.3,A8)') 'Total wall time for MP3 = ',t_MP,' seconds'
|
||||
write(*,*)
|
||||
|
||||
end if
|
||||
|
||||
194
src/MP/RMP3.f90
194
src/MP/RMP3.f90
@ -1,194 +0,0 @@
|
||||
subroutine RMP3(nBasin,nCin,nOin,nVin,nRin,ERI,ENuc,EHF,e)
|
||||
|
||||
! Perform third-order Moller-Plesset calculation
|
||||
|
||||
implicit none
|
||||
|
||||
! Input variables
|
||||
|
||||
integer,intent(in) :: nBasin
|
||||
integer,intent(in) :: nCin
|
||||
integer,intent(in) :: nOin
|
||||
integer,intent(in) :: nVin
|
||||
integer,intent(in) :: nRin
|
||||
double precision,intent(in) :: ENuc,EHF
|
||||
double precision,intent(in) :: e(nBasin)
|
||||
double precision,intent(in) :: ERI(nBasin,nBasin,nBasin,nBasin)
|
||||
|
||||
! Local variables
|
||||
|
||||
double precision :: eps,E2,EcMP2
|
||||
double precision :: eps1,eps2,E3a,E3b,E3c
|
||||
double precision :: EcMP3
|
||||
|
||||
integer :: nBas
|
||||
integer :: nC
|
||||
integer :: nO
|
||||
integer :: nV
|
||||
integer :: nR
|
||||
integer :: i,j,k,l,a,b,c,d
|
||||
|
||||
double precision,allocatable :: se(:)
|
||||
double precision,allocatable :: eO(:)
|
||||
double precision,allocatable :: eV(:)
|
||||
|
||||
double precision,allocatable :: sERI(:,:,:,:)
|
||||
double precision,allocatable :: dbERI(:,:,:,:)
|
||||
|
||||
double precision,allocatable :: OOOO(:,:,:,:)
|
||||
double precision,allocatable :: OOVV(:,:,:,:)
|
||||
double precision,allocatable :: OVVO(:,:,:,:)
|
||||
double precision,allocatable :: VVOO(:,:,:,:)
|
||||
double precision,allocatable :: VVVV(:,:,:,:)
|
||||
|
||||
|
||||
! Hello world
|
||||
|
||||
write(*,*)
|
||||
write(*,*)'******************************'
|
||||
write(*,*)'* Restricted MP3 Calculation *'
|
||||
write(*,*)'******************************'
|
||||
write(*,*)
|
||||
|
||||
! Spatial to spin orbitals
|
||||
|
||||
nBas = 2*nBasin
|
||||
nC = 2*nCin
|
||||
nO = 2*nOin
|
||||
nV = 2*nVin
|
||||
nR = 2*nRin
|
||||
|
||||
allocate(se(nBas),sERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call spatial_to_spin_MO_energy(nBasin,e,nBas,se)
|
||||
call spatial_to_spin_ERI(nBasin,ERI,nBas,sERI)
|
||||
|
||||
! Antysymmetrize ERIs
|
||||
|
||||
allocate(dbERI(nBas,nBas,nBas,nBas))
|
||||
|
||||
call antisymmetrize_ERI(2,nBas,sERI,dbERI)
|
||||
|
||||
deallocate(sERI)
|
||||
|
||||
! Form energy denominator
|
||||
|
||||
allocate(eO(nO),eV(nV))
|
||||
|
||||
eO(:) = se(1:nO)
|
||||
eV(:) = se(nO+1:nBas)
|
||||
|
||||
deallocate(se)
|
||||
|
||||
! Create integral batches
|
||||
|
||||
allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVOO(nV,nV,nO,nO),VVVV(nV,nV,nV,nV))
|
||||
|
||||
OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO )
|
||||
OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas,nO+1:nBas)
|
||||
OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas,nO+1:nBas, 1:nO )
|
||||
VVOO(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas, 1:nO , 1:nO )
|
||||
VVVV(:,:,:,:) = dbERI(nO+1:nBas,nO+1:nBas,nO+1:nBas,nO+1:nBas)
|
||||
|
||||
deallocate(dbERI)
|
||||
|
||||
! Compute MP2 energy
|
||||
|
||||
E2 = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
eps = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
|
||||
E2 = E2 + OOVV(i,j,a,b)*OOVV(i,j,a,b)/eps
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
EcMP2 = 0.25d0*E2
|
||||
|
||||
! Compute MP3 energy
|
||||
|
||||
E3a = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do l=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
|
||||
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
eps2 = eO(k) + eO(l) - eV(a) - eV(b)
|
||||
|
||||
E3a = E3a + OOVV(i,j,a,b)*OOOO(k,l,i,j)*VVOO(a,b,k,l)/(eps1*eps2)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
E3b = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
do d=1,nV-nR
|
||||
|
||||
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
eps2 = eO(i) + eO(j) - eV(c) - eV(d)
|
||||
|
||||
E3b = E3b + OOVV(i,j,a,b)*VVVV(a,b,c,d)*VVOO(c,d,i,j)/(eps1*eps2)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
E3c = 0d0
|
||||
|
||||
do i=nC+1,nO
|
||||
do j=nC+1,nO
|
||||
do k=nC+1,nO
|
||||
do a=1,nV-nR
|
||||
do b=1,nV-nR
|
||||
do c=1,nV-nR
|
||||
|
||||
eps1 = eO(i) + eO(j) - eV(a) - eV(b)
|
||||
eps2 = eO(i) + eO(k) - eV(a) - eV(c)
|
||||
|
||||
E3c = E3c + OOVV(i,j,a,b)*OVVO(k,b,c,j)*VVOO(a,c,i,k)/(eps1*eps2)
|
||||
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
end do
|
||||
|
||||
EcMP3 = 0.125d0*E3a + 0.125d0*E3b + E3c
|
||||
|
||||
write(*,*)
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,'(A32)') ' MP3 calculation '
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' MP2 contribution ',EcMP2
|
||||
write(*,'(A32,1X,F16.10)') ' MP3 contribution ',EcMP3
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,'(A32,1X,F16.10)') ' MP3 correlation energy', EcMP2 + EcMP3
|
||||
write(*,'(A32,1X,F16.10)') ' MP3 total energy',ENuc + EHF + EcMP2 + EcMP3
|
||||
write(*,'(A32)') '-----------------------'
|
||||
write(*,*)
|
||||
|
||||
end subroutine
|
||||
@ -30,7 +30,7 @@ program QuAcK
|
||||
double precision,allocatable :: T(:,:)
|
||||
double precision,allocatable :: V(:,:)
|
||||
double precision,allocatable :: Hc(:,:)
|
||||
double precision,allocatable :: X(:,:)
|
||||
double precision,allocatable :: X(:,:),X_tmp(:,:)
|
||||
double precision,allocatable :: dipole_int_AO(:,:,:)
|
||||
double precision,allocatable :: ERI_AO(:,:,:,:)
|
||||
double precision,allocatable :: Uvec(:,:), Uval(:)
|
||||
@ -71,10 +71,6 @@ program QuAcK
|
||||
|
||||
logical :: dotest,doRtest,doUtest,doGtest
|
||||
|
||||
integer :: i, j, j0
|
||||
double precision :: acc_d, acc_nd
|
||||
double precision, allocatable :: tmp1(:,:), tmp2(:,:)
|
||||
|
||||
!-------------!
|
||||
! Hello World !
|
||||
!-------------!
|
||||
@ -180,61 +176,11 @@ program QuAcK
|
||||
|
||||
! Compute orthogonalization matrix
|
||||
|
||||
!call orthogonalization_matrix(nBas, S, X)
|
||||
|
||||
allocate(Uvec(nBas,nBas), Uval(nBas))
|
||||
|
||||
Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas)
|
||||
call diagonalize_matrix(nBas, Uvec, Uval)
|
||||
|
||||
nOrb = 0
|
||||
do i = 1, nBas
|
||||
if(Uval(i) > 1d-6) then
|
||||
Uval(i) = 1d0 / dsqrt(Uval(i))
|
||||
nOrb = nOrb + 1
|
||||
else
|
||||
write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization'
|
||||
end if
|
||||
end do
|
||||
|
||||
write(*,'(A38)') '--------------------------------------'
|
||||
write(*,'(A38,1X,I16)') 'Number of basis functions (AOs)', nBas
|
||||
write(*,'(A38,1X,I16)') 'Number of basis functions (MOs)', nOrb
|
||||
write(*,'(A38,1X,F9.3)') ' % of discarded orbitals = ', 100.d0 * (1.d0 - dble(nOrb)/dble(nBas))
|
||||
write(*,'(A38)') '--------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
j0 = nBas - nOrb
|
||||
allocate(X_tmp(nBas,nBas))
|
||||
call orthogonalization_matrix(nBas,nOrb,S,X_tmp)
|
||||
allocate(X(nBas,nOrb))
|
||||
do j = j0+1, nBas
|
||||
do i = 1, nBas
|
||||
X(i,j-j0) = Uvec(i,j) * Uval(j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
deallocate(Uvec, Uval)
|
||||
|
||||
!! check if X.T S X = 1_(nOrb,nOrb)
|
||||
!allocate(tmp1(nOrb,nBas), tmp2(nOrb,nOrb))
|
||||
!call dgemm("T", "N", nOrb, nBas, nBas, 1.d0, &
|
||||
! X(1,1), nBas, S(1,1), nBas, &
|
||||
! 0.d0, tmp1(1,1), nOrb)
|
||||
!call dgemm("N", "N", nOrb, nOrb, nBas, 1.d0, &
|
||||
! tmp1(1,1), nOrb, X(1,1), nBas, &
|
||||
! 0.d0, tmp2(1,1), nOrb)
|
||||
!acc_d = 0.d0
|
||||
!acc_nd = 0.d0
|
||||
!do i = 1, nOrb
|
||||
! !write(*,'(1000(F15.7,2X))') (tmp2(i,j), j = 1, nOrb)
|
||||
! acc_d = acc_d + tmp2(i,i)
|
||||
! do j = 1, nOrb
|
||||
! if(j == i) cycle
|
||||
! acc_nd = acc_nd + dabs(tmp2(j,i))
|
||||
! enddo
|
||||
!enddo
|
||||
!print*, ' diag part: ', dabs(acc_d - dble(nOrb)) / dble(nOrb)
|
||||
!print*, ' non-diag part: ', acc_nd
|
||||
!deallocate(tmp1, tmp2)
|
||||
X(1:nBas,1:nOrb) = X_tmp(1:nBas,1:nOrb)
|
||||
deallocate(X_tmp)
|
||||
|
||||
!---------------------!
|
||||
! Choose QuAcK branch !
|
||||
|
||||
@ -1,7 +1,4 @@
|
||||
|
||||
! ---
|
||||
|
||||
subroutine orthogonalization_matrix(nBas, S, X)
|
||||
subroutine orthogonalization_matrix(nBas,nOrb,S,X)
|
||||
|
||||
! Compute the orthogonalization matrix X
|
||||
|
||||
@ -17,113 +14,47 @@ subroutine orthogonalization_matrix(nBas, S, X)
|
||||
logical :: debug
|
||||
double precision,allocatable :: UVec(:,:),Uval(:)
|
||||
double precision,parameter :: thresh = 1d-6
|
||||
integer,parameter :: ortho_type = 1
|
||||
|
||||
integer :: i
|
||||
integer :: i,j,j0
|
||||
|
||||
! Output variables
|
||||
|
||||
integer :: nOrb
|
||||
double precision,intent(out) :: X(nBas,nBas)
|
||||
|
||||
debug = .false.
|
||||
|
||||
! Type of orthogonalization ortho_type
|
||||
!
|
||||
! 1 = Lowdin
|
||||
! 2 = Canonical
|
||||
! 3 = SVD
|
||||
!
|
||||
|
||||
allocate(Uvec(nBas,nBas),Uval(nBas))
|
||||
|
||||
if(ortho_type == 1) then
|
||||
|
||||
!
|
||||
! S V = V s where
|
||||
!
|
||||
! V.T V = 1 and s > 0 (S is positive def)
|
||||
!
|
||||
! S = V s V.T
|
||||
! = V s^0.5 s^0.5 V.T
|
||||
! = V s^0.5 V.T V s^0.5 V.T
|
||||
! = S^0.5 S^0.5
|
||||
!
|
||||
! where
|
||||
!
|
||||
! S^0.5 = V s^0.5 V.T
|
||||
!
|
||||
! X = S^(-0.5)
|
||||
! = V s^(-0.5) V.T
|
||||
!
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*) ' Lowdin orthogonalization'
|
||||
! write(*,*)
|
||||
|
||||
Uvec = S
|
||||
call diagonalize_matrix(nBas, Uvec, Uval)
|
||||
|
||||
do i = 1, nBas
|
||||
|
||||
if(Uval(i) < thresh) then
|
||||
|
||||
write(*,*) 'Eigenvalue',i,' is very small in Lowdin orthogonalization = ',Uval(i)
|
||||
|
||||
end if
|
||||
|
||||
Uval(i) = 1d0 / dsqrt(Uval(i))
|
||||
|
||||
end do
|
||||
|
||||
call ADAt(nBas, Uvec(1,1), Uval(1), X(1,1))
|
||||
|
||||
elseif(ortho_type == 2) then
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*) 'Canonical orthogonalization'
|
||||
! write(*,*)
|
||||
|
||||
Uvec = S
|
||||
call diagonalize_matrix(nBas, Uvec, Uval)
|
||||
|
||||
do i = 1, nBas
|
||||
Uvec(1:nBas,1:nBas) = S(1:nBas,1:nBas)
|
||||
call diagonalize_matrix(nBas,Uvec,Uval)
|
||||
|
||||
nOrb = 0
|
||||
do i = 1,nBas
|
||||
if(Uval(i) > thresh) then
|
||||
|
||||
Uval(i) = 1d0 / dsqrt(Uval(i))
|
||||
|
||||
nOrb = nOrb + 1
|
||||
else
|
||||
|
||||
write(*,*) ' Eigenvalue',i,'too small for canonical orthogonalization'
|
||||
|
||||
end if
|
||||
|
||||
end do
|
||||
|
||||
call AD(nBas, Uvec, Uval)
|
||||
X = Uvec
|
||||
write(*,'(A50)') '------------------------------------------------'
|
||||
write(*,'(A40,1X,I5)') 'Number of basis functions = ',nBas
|
||||
write(*,'(A40,1X,I5)') 'Number of spatial orbitals = ',nOrb
|
||||
write(*,'(A40,1X,I5)') 'Number of discarded functions = ',nBas - nOrb
|
||||
write(*,'(A50)') '------------------------------------------------'
|
||||
write(*,*)
|
||||
|
||||
elseif(ortho_type == 3) then
|
||||
j0 = nBas - nOrb
|
||||
|
||||
! write(*,*)
|
||||
! write(*,*) ' SVD-based orthogonalization NYI'
|
||||
! write(*,*)
|
||||
do j = j0+1,nBas
|
||||
do i = 1,nBas
|
||||
X(i,j-j0) = Uvec(i,j) * Uval(j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Uvec = S
|
||||
! call diagonalize_matrix(nBas,Uvec,Uval)
|
||||
|
||||
! do i=1,nBas
|
||||
! if(Uval(i) > thresh) then
|
||||
! Uval(i) = 1d0/sqrt(Uval(i))
|
||||
! else
|
||||
! write(*,*) 'Eigenvalue',i,'too small for canonical orthogonalization'
|
||||
! end if
|
||||
! end do
|
||||
!
|
||||
! call AD(nBas,Uvec,Uval)
|
||||
! X = Uvec
|
||||
|
||||
end if
|
||||
deallocate(Uvec,Uval)
|
||||
|
||||
! Print results
|
||||
|
||||
|
||||
Loading…
Reference in New Issue
Block a user