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

safer but slower orthogonalization in pp-RPA

Browse Source
This commit is contained in:
Pierre-Francois Loos 2021-10-21 22:40:20 +02:00
parent 52a7e632c4
commit c4ac7f69e6
6 changed files with 76 additions and 68 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
input/methods
View File

@ -13,9 +13,9 @@
# G0F2* evGF2* qsGF2* G0F3 evGF3 # G0F2* evGF2* qsGF2* G0F3 evGF3
F F F F F F F F F F
# G0W0* evGW* qsGW* ufG0W0 ufGW # G0W0* evGW* qsGW* ufG0W0 ufGW
T F F F F F F F F F
# G0T0 evGT qsGT # G0T0 evGT qsGT
F F F T F F
# MCMP2 # MCMP2
F F
# * unrestricted version available # * unrestricted version available

4
input/options
View File

@ -5,7 +5,7 @@
# CC: maxSCF thresh DIIS n_diis # CC: maxSCF thresh DIIS n_diis
64 0.0000000001 T 5 64 0.0000000001 T 5
# spin: TDA singlet triplet spin_conserved spin_flip # spin: TDA singlet triplet spin_conserved spin_flip
T T F T T F T T T T
# GF: maxSCF thresh DIIS n_diis lin eta renorm # GF: maxSCF thresh DIIS n_diis lin eta renorm
256 0.00001 T 5 T 0.0 3 256 0.00001 T 5 T 0.0 3
# GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 # GW/GT: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0
@ -13,6 +13,6 @@
# ACFDT: AC Kx XBS # ACFDT: AC Kx XBS
F F T F F T
# BSE: BSE dBSE dTDA evDyn # BSE: BSE dBSE dTDA evDyn
T F T F T T T F
# MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift
1000000 100000 10 0.3 10000 1234 T 1000000 100000 10 0.3 10000 1234 T

2
mol/carbon_monoxide.xyz
View File

@ -1,4 +1,4 @@
2 2
Carbon monoxide,^1\Sigma^+,CC3,aug-cc-pVTZ
C 0.00000000 0.00000000 -0.66116488 C 0.00000000 0.00000000 -0.66116488
O 0.00000000 0.00000000 0.47237899 O 0.00000000 0.00000000 0.47237899

2
mol/hydrogen_chloride.xyz
View File

@ -1,4 +1,4 @@
2 2
Hydrogen chloride,^1\Sigma^+,CC3,aug-cc-pVTZ
Cl 0.00000000 0.00000000 -0.01317536 Cl 0.00000000 0.00000000 -0.01317536
H 0.00000000 0.00000000 1.26199843 H 0.00000000 0.00000000 1.26199843

10
src/LR/linear_response_pp.f90
View File

@ -95,6 +95,16 @@ subroutine linear_response_pp(ispin,TDA,nBas,nC,nO,nV,nR,nOO,nVV,lambda,e,ERI,Om
call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:)) call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:))
! print*,'Omega1'
! call matout(nVV,1,Omega1)
! print*,'X1t.X1 - Y1t.Y1'
! call matout(nVV,nVV,matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1))
! print*,'Omega2'
! call matout(nOO,1,Omega2)
! print*,'Y2t.Y2 - X2t.X2'
! call matout(nOO,nOO,matmul(transpose(Y2),Y2) - matmul(transpose(X2),X2))
end if end if
! Compute the RPA correlation energy ! Compute the RPA correlation energy

118
src/RPA/sort_ppRPA.f90
View File

@ -116,99 +116,97 @@ subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
! Orthogonalize eigenvectors ! Orthogonalize eigenvectors
! ! Find degenerate eigenvalues ! deg1 = 1
! ab_start = 1
deg1 = 1 ! do ab=2,nVV
ab_start = 1
do ab=2,nVV ! if(abs(Omega1(ab-1) - Omega1(ab)) < 1d-6) then
if(abs(Omega1(ab-1) - Omega1(ab)) < 1d-6) then ! deg1 = deg1 + 1
deg1 = deg1 + 1 ! if(ab == nVV) then
if(ab == nVV) then ! ab_end = ab
ab_end = ab ! 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))) ! deallocate(S1,O1)
call orthogonalization_matrix(1,deg1,S1,O1)
Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1)
deallocate(S1,O1) ! end if
end if ! else
else ! ab_end = ab - 1
ab_end = ab - 1 ! 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))) ! deallocate(S1,O1)
call orthogonalization_matrix(1,deg1,S1,O1)
Z1(:,ab_start:ab_end) = matmul(Z1(:,ab_start:ab_end),O1)
deallocate(S1,O1) ! deg1 = 1
! ab_start = ab
deg1 = 1 ! end if
ab_start = ab ! end do
end if ! deg2 = 1
end do ! ij_start = 1
deg2 = 1 ! do ij=2,nOO
ij_start = 1
do ij=2,nOO ! if(abs(Omega2(ij-1) - Omega2(ij)) < 1d-6) then
if(abs(Omega2(ij-1) - Omega2(ij)) < 1d-6) then ! deg2 = deg2 + 1
deg2 = deg2 + 1 ! if(ij == nOO) then
if(ij == nOO) then ! ij_end = ij
!
! 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)
ij_end = ij ! end if
allocate(S2(deg2,deg2),O2(deg2,deg2)) ! else
S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end))) ! ij_end = ij - 1
call orthogonalization_matrix(1,deg2,S2,O2)
Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2)
deallocate(S2,O2) ! allocate(S2(deg2,deg2),O2(deg2,deg2))
end if ! 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)
else ! deallocate(S2,O2)
ij_end = ij - 1 ! deg2 = 1
! ij_start = ij
allocate(S2(deg2,deg2),O2(deg2,deg2)) ! end if
! end do
S2 = - matmul(transpose(Z2(:,ij_start:ij_end)),matmul(M,Z2(:,ij_start:ij_end))) allocate(S1(nVV,nVV),S2(nOO,nOO),O1(nVV,nVV),O2(nOO,nOO))
call orthogonalization_matrix(1,deg2,S2,O2) S1 = + matmul(transpose(Z1),matmul(M,Z1))
Z2(:,ij_start:ij_end) = matmul(Z2(:,ij_start:ij_end),O2) S2 = - matmul(transpose(Z2),matmul(M,Z2))
deallocate(S2,O2) if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
deg2 = 1 Z1 = matmul(Z1,O1)
ij_start = ij Z2 = matmul(Z2,O2)
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)
! Define submatrices X1, Y1, X2, & Y2 ! Define submatrices X1, Y1, X2, & Y2