Skip to content

Commit aaaadf5

Browse files
author
lburth
committed
Fix LR at MOM refernce by staying (mostly) real
1 parent d93b237 commit aaaadf5

9 files changed

+225
-75
lines changed

‎src/LR/CVS_phLR_transition_vectors.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ subroutine CVS_phLR_transition_vectors(spin_allowed,nOrb,nC,nO,nV,nR,nS,dipole_i
2222
! Local variables
2323

2424
integer :: ia,jb,j,b
25-
integer :: maxS = 10
25+
integer :: maxS = 30
2626
double precision :: S2
2727
double precision,parameter :: thres_vec = 0.1d0
2828
double precision,allocatable :: X(:)

‎src/LR/CVS_phRLR.f90‎

Lines changed: 54 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -15,8 +15,11 @@ subroutine CVS_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
1515
! Local variables
1616

1717
double precision,external :: trace_matrix
18-
complex*16,allocatable :: RPA_matrix(:,:)
19-
complex*16,allocatable :: OmOmminus(:)
18+
complex*16,allocatable :: complex_RPA_matrix(:,:)
19+
double precision,allocatable :: RPA_matrix(:,:)
20+
double precision,allocatable :: vectors(:,:)
21+
double precision,allocatable :: OmOmminus(:)
22+
integer :: i
2023

2124
! Output variables
2225

@@ -38,32 +41,62 @@ subroutine CVS_phRLR(TDA,nS,Aph,Bph,EcRPA,Om,XpY,XmY)
3841

3942
else
4043

41-
allocate(RPA_matrix(2*nS,2*nS),OmOmminus(2*nS))
42-
RPA_matrix(1:nS,1:nS) = cmplx(Aph(:,:),0d0,kind=8)
43-
RPA_matrix(1:nS,nS+1:2*nS) = cmplx(Bph(:,:),0d0,kind=8)
44-
RPA_matrix(nS+1:2*nS,1:nS) = cmplx(-Bph(:,:),0d0,kind=8)
45-
RPA_matrix(nS+1:2*nS,nS+1:2*nS) = cmplx(-Aph(:,:),0d0,kind=8)
46-
call complex_diagonalize_matrix_without_sort(2*nS,RPA_matrix,OmOmminus)
47-
call complex_sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix)
48-
call complex_normalize_RPA(nS,RPA_matrix)
49-
call rotate_vectors_to_real_axis(2*nS,nS,RPA_matrix(1:2*nS,1:nS))
50-
51-
Om(:) = real(OmOmminus(1:nS))
44+
allocate(RPA_matrix(2*nS,2*nS),vectors(2*nS,2*nS),OmOmminus(2*nS))
45+
46+
RPA_matrix(1:nS,1:nS) = Aph(:,:)
47+
RPA_matrix(1:nS,nS+1:2*nS) = Bph(:,:)
48+
RPA_matrix(nS+1:2*nS,1:nS) = -Bph(:,:)
49+
RPA_matrix(nS+1:2*nS,nS+1:2*nS) = -Aph(:,:)
50+
51+
call diagonalize_general_matrix(2*nS,RPA_matrix,OmOmminus,vectors)
52+
53+
RPA_matrix(:,:) = vectors(:,:)
54+
55+
deallocate(vectors)
56+
57+
call sort_eigenvalues_RPA(2*nS,OmOmminus,RPA_matrix)
58+
59+
allocate(complex_RPA_matrix(2*nS,2*nS))
60+
61+
complex_RPA_matrix = cmplx(0.0d0, 0.0d0, kind=8)
62+
complex_RPA_matrix(1:nS, 1:nS) = cmplx(RPA_matrix(:,:), 0.0d0, kind=8)
63+
complex_RPA_matrix = cmplx(RPA_matrix(:,:),RPA_matrix(:,:)*0d0,kind=8)
64+
65+
deallocate(RPA_matrix)
66+
67+
call complex_normalize_RPA(nS,complex_RPA_matrix)
68+
5269
if(maxval(abs(OmOmminus(1:nS)+OmOmminus(nS+1:2*nS))) > 1e-8) then
5370
call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.')
5471
write(*,*) "Maximal difference :", maxval(abs(OmOmminus(1:nS)+OmOmminus(nS+1:2*nS)))
5572
end if
56-
if(maxval(abs(aimag(OmOmminus(:))))>1d-8)&
57-
call print_warning('You may have instabilities in linear response: complex excitation eigenvalue !')
58-
if(maxval(abs(aimag(RPA_matrix(1:2*nS,1:nS))))>1d-8) then
59-
print *, "Max imag value X+Y:",maxval(aimag(RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,1:nS)))
60-
print *, "Max imag value X-Y:",maxval(aimag(RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,1:nS)))
73+
74+
Om(:) = OmOmminus(1:nS)
75+
76+
deallocate(OmOmminus)
77+
78+
! Set transition vectors of zero modes to 0
79+
do i=1,nS
80+
if(abs(Om(i))<1d-8) then
81+
complex_RPA_matrix(:,i) = (0d0,0d0)
82+
complex_RPA_matrix(:,i+nS) = (0d0,0d0)
83+
end if
84+
end do
85+
86+
87+
if(maxval(aimag(complex_RPA_matrix(1:2*nS,1:nS)))>1d-8) then
88+
call print_warning('You may have instabilities in linear response: complex transition vectors !')
89+
print *, "Max imag value X+Y:",maxval(aimag(transpose(complex_RPA_matrix(1:nS,1:nS) + complex_RPA_matrix(nS+1:2*nS,1:nS)))),&
90+
"at" ,maxloc(aimag(transpose(complex_RPA_matrix(1:nS,1:nS) + complex_RPA_matrix(nS+1:2*nS,1:nS))))
91+
92+
print *, "Max imag value X-Y:",maxval(aimag(transpose(complex_RPA_matrix(1:nS,1:nS) - complex_RPA_matrix(nS+1:2*nS,1:nS)))),&
93+
"at" ,maxloc(aimag(transpose(complex_RPA_matrix(1:nS,1:nS) - complex_RPA_matrix(nS+1:2*nS,1:nS))))
6194
end if
6295

63-
XpY(:,:) = transpose(real(RPA_matrix(1:nS,1:nS) + RPA_matrix(nS+1:2*nS,1:nS)))
64-
XmY(:,:) = transpose(real(RPA_matrix(1:nS,1:nS) - RPA_matrix(nS+1:2*nS,1:nS)))
96+
XpY(:,:) = transpose(real(complex_RPA_matrix(1:nS,1:nS) + complex_RPA_matrix(nS+1:2*nS,1:nS)))
97+
XmY(:,:) = transpose(real(complex_RPA_matrix(1:nS,1:nS) - complex_RPA_matrix(nS+1:2*nS,1:nS)))
6598

66-
deallocate(RPA_matrix,OmOmminus)
99+
deallocate(complex_RPA_matrix)
67100

68101
end if
69102

‎src/LR/CVS_phULR.f90‎

Lines changed: 53 additions & 21 deletions
Original file line numberDiff line numberDiff line change
@@ -17,8 +17,11 @@ subroutine CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
1717
! Local variables
1818

1919
double precision,external :: trace_matrix
20-
complex*16,allocatable :: RPA_matrix(:,:)
21-
complex*16,allocatable :: OmOmminus(:)
20+
double precision,allocatable :: RPA_matrix(:,:)
21+
complex*16,allocatable :: complex_RPA_matrix(:,:)
22+
double precision,allocatable :: vectors(:,:)
23+
double precision,allocatable :: OmOmminus(:)
24+
integer :: i
2225

2326
! Output variables
2427

@@ -38,34 +41,63 @@ subroutine CVS_phULR(TDA,nSa,nSb,nSt,Aph,Bph,EcRPA,Om,XpY,XmY)
3841

3942
else
4043

41-
allocate(RPA_matrix(2*nSt,2*nSt),OmOmminus(2*nSt))
44+
allocate(RPA_matrix(2*nSt,2*nSt),vectors(2*nSt,2*nSt),OmOmminus(2*nSt))
4245

43-
RPA_matrix(1:nSt,1:nSt) = cmplx( Aph(:,:),0d0,kind=8)
44-
RPA_matrix(1:nSt,nSt+1:2*nSt) = cmplx( Bph(:,:),0d0,kind=8)
45-
RPA_matrix(nSt+1:2*nSt,1:nSt) = cmplx(-Bph(:,:),0d0,kind=8)
46-
RPA_matrix(nSt+1:2*nSt,nSt+1:2*nSt) = cmplx(-Aph(:,:),0d0,kind=8)
46+
RPA_matrix(1:nSt,1:nSt) = Aph(:,:)
47+
RPA_matrix(1:nSt,nSt+1:2*nSt) = Bph(:,:)
48+
RPA_matrix(nSt+1:2*nSt,1:nSt) = -Bph(:,:)
49+
RPA_matrix(nSt+1:2*nSt,nSt+1:2*nSt) = -Aph(:,:)
50+
51+
call diagonalize_general_matrix(2*nSt,RPA_matrix,OmOmminus,vectors)
52+
53+
RPA_matrix(:,:) = vectors(:,:)
54+
55+
deallocate(vectors)
56+
57+
call sort_eigenvalues_RPA(2*nSt,OmOmminus,RPA_matrix)
58+
59+
allocate(complex_RPA_matrix(2*nSt,2*nSt))
60+
61+
complex_RPA_matrix = cmplx(0.0d0, 0.0d0, kind=8)
62+
complex_RPA_matrix(1:nSt, 1:nSt) = cmplx(RPA_matrix(:,:), 0.0d0, kind=8)
63+
complex_RPA_matrix = cmplx(RPA_matrix(:,:),RPA_matrix(:,:)*0d0,kind=8)
64+
65+
deallocate(RPA_matrix)
66+
67+
call complex_normalize_RPA(nSt,complex_RPA_matrix)
4768

48-
call complex_diagonalize_matrix_without_sort(2*nSt,RPA_matrix,OmOmminus)
49-
call complex_sort_eigenvalues_RPA(2*nSt,OmOmminus,RPA_matrix)
50-
call complex_normalize_RPA(nSt,RPA_matrix)
51-
call rotate_vectors_to_real_axis(2*nSt,nSt,RPA_matrix(1:2*nSt,1:nSt))
52-
53-
Om(:) = real(OmOmminus(1:nSt))
5469
if(maxval(abs(OmOmminus(1:nSt)+OmOmminus(nSt+1:2*nSt))) > 1e-8) then
5570
call print_warning('We dont find a Om and -Om structure as solution of the RPA. There might be a problem somewhere.')
5671
write(*,*) "Maximal difference :", maxval(abs(OmOmminus(1:nSt)+OmOmminus(nSt+1:2*nSt)))
5772
end if
58-
if(maxval(abs(aimag(OmOmminus(:))))>1d-8)&
59-
call print_warning('You may have instabilities in linear response: complex excitation eigenvalue !')
60-
if(maxval(aimag(RPA_matrix(1:2*nSt,1:nSt)))>1d-8) then
73+
74+
Om(:) = OmOmminus(1:nSt)
75+
76+
deallocate(OmOmminus)
77+
78+
! Set transition vectors of zero modes to 0
79+
do i=1,nSt
80+
if(abs(Om(i))<1d-8) then
81+
complex_RPA_matrix(:,i) = (0d0,0d0)
82+
complex_RPA_matrix(:,i+nSt) = (0d0,0d0)
83+
end if
84+
end do
85+
86+
87+
if(maxval(aimag(complex_RPA_matrix(1:2*nSt,1:nSt)))>1d-8) then
6188
call print_warning('You may have instabilities in linear response: complex transition vectors !')
62-
print *, "Max imag value X+Y:",maxval(aimag(transpose(RPA_matrix(1:nSt,1:nSt) + RPA_matrix(nSt+1:2*nSt,1:nSt))))
63-
print *, "Max imag value X-Y:",maxval(aimag(transpose(RPA_matrix(1:nSt,1:nSt) - RPA_matrix(nSt+1:2*nSt,1:nSt))))
89+
print *, "Max imag value X+Y:",maxval(aimag(transpose(complex_RPA_matrix(1:nSt,1:nSt) + complex_RPA_matrix(nSt+1:2*nSt,1:nSt)))),&
90+
"at" ,maxloc(aimag(transpose(complex_RPA_matrix(1:nSt,1:nSt) + complex_RPA_matrix(nSt+1:2*nSt,1:nSt))))
91+
92+
print *, "Max imag value X-Y:",maxval(aimag(transpose(complex_RPA_matrix(1:nSt,1:nSt) - complex_RPA_matrix(nSt+1:2*nSt,1:nSt)))),&
93+
"at" ,maxloc(aimag(transpose(complex_RPA_matrix(1:nSt,1:nSt) - complex_RPA_matrix(nSt+1:2*nSt,1:nSt))))
6494
end if
65-
XpY(:,:) = transpose(real(RPA_matrix(1:nSt,1:nSt) + RPA_matrix(nSt+1:2*nSt,1:nSt)))
66-
XmY(:,:) = transpose(real(RPA_matrix(1:nSt,1:nSt) - RPA_matrix(nSt+1:2*nSt,1:nSt)))
67-
deallocate(RPA_matrix,OmOmminus)
6895

96+
XpY(:,:) = transpose(real(complex_RPA_matrix(1:nSt,1:nSt) + complex_RPA_matrix(nSt+1:2*nSt,1:nSt)))
97+
XmY(:,:) = transpose(real(complex_RPA_matrix(1:nSt,1:nSt) - complex_RPA_matrix(nSt+1:2*nSt,1:nSt)))
98+
99+
deallocate(complex_RPA_matrix)
100+
69101
end if
70102

71103
! Compute the RPA correlation energy

‎src/LR/CVS_phULR_transition_vectors.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -32,7 +32,7 @@ subroutine CVS_phULR_transition_vectors(ispin,nBas,nC,nO,nV,nR,nS,nSa,nSb,nSt,&
3232
! Local variables
3333

3434
integer :: ia,jb,j,b
35-
integer :: maxS = 10
35+
integer :: maxS = 30
3636
double precision,parameter :: thres_vec = 0.1d0
3737
double precision,allocatable :: X(:)
3838
double precision,allocatable :: Y(:)

‎src/LR/complex_print_excitation_energies.f90‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -14,7 +14,7 @@ subroutine complex_print_excitation_energies(method,manifold,nS,Om)
1414

1515
! Local variables
1616

17-
integer,parameter :: maxS = 300
17+
integer,parameter :: maxS = 50
1818
integer :: m
1919

2020
write(*,*)

‎src/RPA/CVS_phURPAx.f90‎

Lines changed: 0 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -76,8 +76,6 @@ subroutine CVS_phURPAx(dotest,TDA,doACFDT,exchange_kernel,spin_conserved,spin_fl
7676

7777
nFC(1) = MERGE(1,0,FC(1)/=0)
7878
nFC(2) = MERGE(1,0,FC(2)/=0)
79-
print *, nFC
80-
print *, FC
8179
allocate(occupations_fc(maxval(nO-nFC),nspin))
8280
! remove FC from occupations
8381
do ispin=1,nspin

‎src/utils/rotate_vectors_to_real_axis.f90‎

Lines changed: 0 additions & 27 deletions
This file was deleted.

‎src/utils/sort_eigenvalues_RPA.f90‎

Lines changed: 112 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,112 @@
1+
subroutine sort_eigenvalues_RPA(n, evals, evecs)
2+
implicit none
3+
integer, intent(in) :: n
4+
double precision, intent(inout) :: evals(n)
5+
double precision, intent(inout) :: evecs(n, n)
6+
integer :: i, j, k, nhalf, counter
7+
double precision :: temp_val
8+
double precision :: temp_vec(n)
9+
double precision :: etanorm
10+
double precision :: threshold = 1d-8
11+
12+
! -------------------------------
13+
! Step 1: initial bubble sort
14+
! Sort in descending order by real part
15+
! -------------------------------
16+
do i = 1, n-1
17+
do j = i+1, n
18+
19+
if (evals(i) < evals(j)) then
20+
! swap
21+
temp_val = evals(i)
22+
evals(i) = evals(j)
23+
evals(j) = temp_val
24+
25+
temp_vec = evecs(:, i)
26+
evecs(:, i) = evecs(:, j)
27+
evecs(:, j) = temp_vec
28+
end if
29+
30+
end do
31+
end do
32+
33+
! -------------------------------
34+
! Step 2: classify first half
35+
! - eta-norm > 0 -> excitation
36+
! - eta-norm < 0 -> swap
37+
! -------------------------------
38+
nhalf = n/2
39+
counter = 0
40+
41+
do i = 1, nhalf
42+
43+
etanorm = sum(abs(evecs(1:nhalf,i))**2) - &
44+
sum(abs(evecs(nhalf+1:n,i))**2)
45+
46+
if(abs(etanorm) < threshold) then
47+
print *, 'Mode with vanishing eta-norm has been found !'
48+
print *, 'Omega:', evals(i)
49+
print *, 'eta-norm', etanorm
50+
end if
51+
52+
! normal case: classify by eta-norm
53+
if (etanorm < -threshold) then
54+
55+
j = n + 1 - i
56+
57+
temp_val = evals(i)
58+
evals(i) = evals(j)
59+
evals(j) = temp_val
60+
61+
temp_vec = evecs(:, i)
62+
evecs(:, i) = evecs(:, j)
63+
evecs(:, j) = temp_vec
64+
65+
counter = counter + 1
66+
67+
end if
68+
69+
end do
70+
71+
if(counter /= 0) then
72+
print *, "Attention:", counter, "excitations with negative excitation energies have been found!"
73+
end if
74+
75+
! -------------------------------
76+
! Step 3: separate positives and negatives
77+
! -------------------------------
78+
! Assuming after step 2, first nhalf are excitations (positive eta norm) and last nhalf are de-excitations (negative eta norm)
79+
80+
! Sort positive excitations ascending
81+
do i = 1, nhalf-1
82+
do j = i+1, nhalf
83+
84+
if (real(evals(i)) > real(evals(j))) then
85+
temp_val = evals(i)
86+
evals(i) = evals(j)
87+
evals(j) = temp_val
88+
temp_vec = evecs(:, i)
89+
evecs(:, i) = evecs(:, j)
90+
evecs(:, j) = temp_vec
91+
end if
92+
93+
end do
94+
end do
95+
96+
! Sort negative de-excitations descending (real value), and descding in imaginary value
97+
do i = nhalf+1, n-1
98+
do j = i+1, n
99+
100+
if (real(evals(i)) < real(evals(j))) then
101+
temp_val = evals(i)
102+
evals(i) = evals(j)
103+
evals(j) = temp_val
104+
temp_vec = evecs(:, i)
105+
evecs(:, i) = evecs(:, j)
106+
evecs(:, j) = temp_vec
107+
end if
108+
109+
end do
110+
end do
111+
112+
end subroutine

‎src/utils/wrap_lapack.f90‎

Lines changed: 3 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -56,7 +56,6 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
5656

5757
! Input variables
5858

59-
integer :: i,j,k
6059
integer,intent(in) :: N
6160
double precision,intent(inout):: A(N,N)
6261
double precision,intent(out) :: VR(N,N)
@@ -80,6 +79,9 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
8079
allocate(work(lwork))
8180

8281
call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
82+
83+
if(any(abs(WI)>1d-8))&
84+
print*, 'Found eigenvalue with imaginary part > 1e-8 (dgeev)!!'
8385

8486
deallocate(work,WI,VL)
8587

0 commit comments

Comments
 (0)