sort ppRPA

input/basis

@@ -1,30 +1,9 @@
-1   6
+1   3
 S   3
-  1     33.8700000              0.0060680        
+  1     38.3600000              0.0238090        
-  2      5.0950000              0.0453080        
+  2      5.7700000              0.1548910        
-  3      1.1590000              0.2028220        
+  3      1.2400000              0.4699870        
 S   1
-  1      0.3258000              1.0000000        
+  1      0.2976000              1.0000000        
-S   1
-  1      0.1027000              1.0000000        
 P   1
-  1      1.4070000              1.0000000        
+  1      1.2750000              1.0000000
-P   1
-  1      0.3880000              1.0000000        
-D   1
-  1      1.0570000              1.0000000
-2   6
-S   3
-  1     33.8700000              0.0060680        
-  2      5.0950000              0.0453080        
-  3      1.1590000              0.2028220        
-S   1
-  1      0.3258000              1.0000000        
-S   1
-  1      0.1027000              1.0000000        
-P   1
-  1      1.4070000              1.0000000        
-P   1
-  1      0.3880000              1.0000000        
-D   1
-  1      1.0570000              1.0000000

input/molecule

@@ -1,5 +1,4 @@
 # nAt nEla nElb nCore nRyd
- 2   1     1   0     0
+  1   1     1   0     0
 # Znuc   x            y           z
-  H     0.           0.         0.
+  He      0.0          0.0         0.0
-  H     0.           0.         1.4

input/molecule.xyz

@@ -1,4 +1,3 @@
-  2
+  1
 
- H      0.0000000000    0.0000000000    0.0000000000
+ He     0.0000000000    0.0000000000    0.0000000000
- H      0.0000000000    0.0000000000    0.7408481486

input/weight

@@ -1,30 +1,9 @@
-1   6
+1   3
 S   3
-  1     33.8700000              0.0060680        
+  1     38.3600000              0.0238090        
-  2      5.0950000              0.0453080        
+  2      5.7700000              0.1548910        
-  3      1.1590000              0.2028220        
+  3      1.2400000              0.4699870        
 S   1
-  1      0.3258000              1.0000000        
+  1      0.2976000              1.0000000        
-S   1
-  1      0.1027000              1.0000000        
 P   1
-  1      1.4070000              1.0000000        
+  1      1.2750000              1.0000000
-P   1
-  1      0.3880000              1.0000000        
-D   1
-  1      1.0570000              1.0000000
-2   6
-S   3
-  1     33.8700000              0.0060680        
-  2      5.0950000              0.0453080        
-  3      1.1590000              0.2028220        
-S   1
-  1      0.3258000              1.0000000        
-S   1
-  1      0.1027000              1.0000000        
-P   1
-  1      1.4070000              1.0000000        
-P   1
-  1      0.3880000              1.0000000        
-D   1
-  1      1.0570000              1.0000000

src/QuAcK/linear_response_pp.f90

@@ -138,7 +138,7 @@ subroutine linear_response_pp(ispin,ortho_eigvec,BSE,nBas,nC,nO,nV,nR,nOO,nVV, &
   EcRPA = 0.5d0*( sum(Omega1(:)) - sum(Omega2(:)) - trace_matrix(nVV,C(:,:)) - trace_matrix(nOO,D(:,:)) )
   EcRPA1 = +sum(Omega1(:)) - trace_matrix(nVV,C(:,:))
   EcRPA2 = -sum(Omega2(:)) - trace_matrix(nOO,D(:,:))
-  if(abs(EcRPA - EcRPA1) > 1d-10 .or. abs(EcRPA - EcRPA2) > 1d-10) & 
+  if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & 
     print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!'
 
 ! write(*,*)'X1'

src/QuAcK/orthogonalization_matrix.f90

@@ -34,9 +34,9 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X)
 
   if(ortho_type == 1) then
 
-    write(*,*)
+!   write(*,*)
-    write(*,*) ' Lowdin orthogonalization'
+!   write(*,*) ' Lowdin orthogonalization'
-    write(*,*)
+!   write(*,*)
 
     Uvec = S
     call diagonalize_matrix(nBas,Uvec,Uval)

src/QuAcK/sort_ppRPA.f90

@@ -16,6 +16,8 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 ! Local variables
 
   integer                       :: pq,ab,ij
+  integer                       :: deg1,ab_start,ab_end
+  integer                       :: deg2,ij_start,ij_end
   double precision,allocatable  :: M(:,:)
   double precision,allocatable  :: Z1(:,:)
   double precision,allocatable  :: Z2(:,:)
@@ -24,6 +26,9 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
   double precision,allocatable  :: O1(:,:)
   double precision,allocatable  :: O2(:,:)
 
+  integer,allocatable           :: order1(:)
+  integer,allocatable           :: order2(:)
+
 ! Output variables
 
   double precision,intent(out)  :: Omega1(nVV)
@@ -37,8 +42,7 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 
  allocate(M(nOO+nVV,nOO+nVV),              &
           Z1(nOO+nVV,nVV),Z2(nOO+nVV,nOO), &
-          S1(nVV,nVV),S2(nOO,nOO),         &
+          order1(nVV),order2(nOO))
-          O1(nVV,nVV),O2(nOO,nOO))
 
 ! Initializatiom
 
@@ -88,6 +92,20 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
   if(minval(Omega1(:)) < 0d0 .or. ab /= nVV) call print_warning('You may have instabilities in pp-RPA!!')
   if(maxval(Omega2(:)) > 0d0 .or. ij /= nOO) call print_warning('You may have instabilities in pp-RPA!!')
 
+  do ab=1,nVV
+    order1(ab) = ab
+  end do
+
+  do ij=1,nOO
+    order2(ij) = ij
+  end do
+
+  call quick_sort(Omega1(:),order1(:),nVV)
+  call set_order(Z1(:,:),order1(:),nOO+nVV,nVV)
+
+  call quick_sort(Omega2(:),order2(:),nOO)
+  call set_order(Z2(:,:),order2(:),nOO+nVV,nOO)
+
 ! write(*,*) 'pp-RPA positive excitation energies'
 ! call matout(nVV,1,Omega1(:))
 ! write(*,*)
@@ -100,14 +118,75 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2)
 
   if(ortho_eigvec) then
 
-    S1 = + matmul(transpose(Z1),matmul(M,Z1))
+    ! Find degenerate eigenvalues
-    S2 = - matmul(transpose(Z2),matmul(M,Z2))
 
-    if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1)
+      deg1 = 1
-    if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2)
+      ab_start = 1
 
-    Z1 = matmul(Z1,O1)
+      do ab=1,nVV
-    Z2 = matmul(Z2,O2)
+
+        if(ab < nVV .and. abs(Omega1(ab) - Omega1(ab+1)) < 1d-10) then 
+
+          if(deg1 == 1) ab_start = ab
+          deg1 = deg1 + 1
+
+        else
+
+          ab_end = ab 
+!         print*,'deg = ',deg1,ab_start,ab_end
+
+          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)
+
+          deallocate(S1,O1)
+
+          deg1 = 1
+          ab_start = ab + 1
+
+        end if
+      end do
+
+      deg2 = 1
+      ij_start = 1
+
+      do ij=1,nOO
+
+        if(ij < nOO .and. abs(Omega2(ij) - Omega2(ij+1)) < 1d-10) then 
+
+          if(deg2 == 1) ij_start = ij 
+          deg2 = deg2 + 1
+
+        else
+
+          ij_end = ij 
+!         print*,'deg = ',deg2,ij_start,ij_end
+
+          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 + 1
+
+        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
 

src/utils/wrap_lapack.f90

@@ -1,35 +1,3 @@
-!subroutine diagonalize_matrix_lowest(N,M,A,e)
-!
-!! Diagonalize a square matrix but only provide the M lowest eigenvalues/eigenvectors
-!
-!  implicit none
-!
-!! Input variables
-!
-!  integer,intent(in)            :: N
-!  integer,intent(in)            :: M
-!  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)
-!  abstol = 1d-15
-!
-!  call dsyevx('V','I','U',N,A,N,VL,VU,1,M,abstol,M,e,work,lwork,info)
-! 
-!  if(info /= 0) then 
-!    print*,'Problem in diagonalize_matrix (dsyev)!!'
-!  endif
-!
-!end subroutine diagonalize_matrix_lowest
-
 subroutine diagonalize_general_matrix(N,A,WR,VR)
 
 ! Diagonalize a non-symmetric square matrix
@@ -53,35 +21,13 @@ subroutine diagonalize_general_matrix(N,A,WR,VR)
 
   lwork = 4*N
   allocate(WI(N),VL(N,N),work(lwork))
-! tmp1 = A
+
   call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info)
 
   if(info /= 0) then 
     print*,'Problem in diagonalize_general_matrix (dgeev)!!'
   endif
 
-! tmp2 = 0d0
-! do i=1,N
-!   do j=1,N
-!     do k=1,N
-!       tmp2(i,j) = tmp2(i,j) + tmp1(i,k)*vr(k,j)
-!     end do
-!   end do
-! end do
-
-! print*,'tmp2'
-! call matout(N,N,tmp2)
-
-! tmp1 = 0d0
-! do i=1,N
-!   do j=1,N
-!     tmp1(i,j) = wr(j)*vr(i,j)
-!   end do
-! end do
-
-! print*,'tmp1'
-! call matout(N,N,tmp1)
-
 end subroutine diagonalize_general_matrix
 
 subroutine diagonalize_matrix(N,A,e)