PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-12-22 20:34:46 +01:00
Code Releases Activity

faster orthogonalization in ppRPA

Browse Source
This commit is contained in:
Pierre-Francois Loos 2020-04-13 20:13:38 +02:00
parent 2f835304f7
commit 57003acb75
4 changed files with 85 additions and 53 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
input/methods
View File

@ -13,6 +13,6 @@
# G0W0 evGW qsGW
F F F
# G0T0 evGT qsGT
T T F
T F F
# MCMP2
F

2
input/options
View File

@ -11,6 +11,6 @@
# GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta
256 0.00001 T 5 F F F F F F T 0.000
# ACFDT: AC Kx XBS
T T T
T F T
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T

20
src/QuAcK/G0T0.f90
View File

@ -184,16 +184,16 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m
! Compute the ppRPA correlation energy
ispin = 1
iblock = 3
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0(:),ERI(:,:,:,:), &
Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin))
ispin = 2
iblock = 4
call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0(:),ERI(:,:,:,:), &
Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin))
EcRPA(1) = EcRPA(1) - EcRPA(2)
EcRPA(2) = 3d0*EcRPA(2)
! ispin = 1
! iblock = 3
! call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0(:),ERI(:,:,:,:), &
! Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin))
! ispin = 2
! iblock = 4
! call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0(:),ERI(:,:,:,:), &
! Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin))
! EcRPA(1) = EcRPA(1) - EcRPA(2)
! EcRPA(2) = 3d0*EcRPA(2)
write(*,*)
write(*,*)'-------------------------------------------------------------------------------'

114
src/QuAcK/sort_ppRPA.f90
View File

@ -129,73 +129,105 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! ! Find degenerate eigenvalues
! deg1 = 1
! ab_start = 1
deg1 = 1
ab_start = 1
! do ab=1,nVV
do ab=2,nVV
! if(ab < nVV .and. abs(Omega1(ab) - Omega1(ab+1)) < 1d-6) then
if(abs(Omega1(ab-1) - Omega1(ab)) < 1d-6) then
! if(deg1 == 1) ab_start = ab
! deg1 = deg1 + 1
deg1 = deg1 + 1
! else
if(ab == nVV) then
ab_end = ab
! ab_end = ab
! print*,'deg = ',deg1,ab_start,ab_end
! allocate(S1(deg1,deg1),O1(deg1,deg1))
allocate(S1(deg1,deg1),O1(deg1,deg1))
! S1 = matmul(transpose(Z1(:,ab_start:ab_end)),matmul(M,Z1(:,ab_start:ab_end)))
! call orthogonalization_matrix(1,deg1,S1,O1)
! Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1)
S1 = matmul(transpose(Z1(:,ab_start:ab_end)),matmul(M,Z1(:,ab_start:ab_end)))
call orthogonalization_matrix(1,deg1,S1,O1)
Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1)
! deallocate(S1,O1)
deallocate(S1,O1)
! deg1 = 1
! ab_start = ab + 1
end if
! end if
! end do
else
! deg2 = 1
! ij_start = 1
ab_end = ab - 1
! do ij=1,nOO
! print*,'deg = ',deg1,ab_start,ab_end
! if(ij < nOO .and. abs(Omega2(ij) - Omega2(ij+1)) < 1d-6) then
allocate(S1(deg1,deg1),O1(deg1,deg1))
! if(deg2 == 1) ij_start = ij
! deg2 = deg2 + 1
S1 = matmul(transpose(Z1(:,ab_start:ab_end)),matmul(M,Z1(:,ab_start:ab_end)))
call orthogonalization_matrix(1,deg1,S1,O1)
Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1)
! else
deallocate(S1,O1)
deg1 = 1
ab_start = ab
end if
end do
deg2 = 1
ij_start = 1
do ij=2,nOO
if(abs(Omega2(ij-1) - Omega2(ij)) < 1d-6) then
deg2 = deg2 + 1
if(ij == nOO) then
ij_end = ij
! ij_end = ij
! print*,'deg = ',deg2,ij_start,ij_end
! allocate(S2(deg2,deg2),O2(deg2,deg2))
allocate(S2(deg2,deg2),O2(deg2,deg2))
! S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end)))
! call orthogonalization_matrix(1,deg2,S2,O2)
! Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2)
S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end)))
call orthogonalization_matrix(1,deg2,S2,O2)
Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2)
! deallocate(S2,O2)
deallocate(S2,O2)
! deg2 = 1
! ij_start = ij + 1
end if
! end if
! end do
else
allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO))
S1 = + matmul(transpose(Z1),matmul(M,Z1))
S2 = - matmul(transpose(Z2),matmul(M,Z2))
ij_end = ij - 1
if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
! print*,'deg = ',deg2,ij_start,ij_end
Z1 = matmul(Z1,O1)
Z2 = matmul(Z2,O2)
allocate(S2(deg2,deg2),O2(deg2,deg2))
S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end)))
call orthogonalization_matrix(1,deg2,S2,O2)
Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2)
deallocate(S2,O2)
deg2 = 1
ij_start = ij
end if
end do
! allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO))
! S1 = + matmul(transpose(Z1),matmul(M,Z1))
! S2 = - matmul(transpose(Z2),matmul(M,Z2))
! if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
! if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
! Z1 = matmul(Z1,O1)
! Z2 = matmul(Z2,O2)
end if