From f700c431d65dccfeb1ad29a778bbeca37ccf5bcb Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Wed, 25 Mar 2020 21:10:21 +0100 Subject: [PATCH] sort ppRPA --- input/basis | 33 ++------- input/molecule | 5 +- input/molecule.xyz | 5 +- input/weight | 33 ++------- src/QuAcK/linear_response_pp.f90 | 2 +- src/QuAcK/orthogonalization_matrix.f90 | 6 +- src/QuAcK/sort_ppRPA.f90 | 95 +++++++++++++++++++++++--- src/utils/wrap_lapack.f90 | 56 +-------------- 8 files changed, 108 insertions(+), 127 deletions(-) diff --git a/input/basis b/input/basis index 69b24c8..6796e3b 100644 --- a/input/basis +++ b/input/basis @@ -1,30 +1,9 @@ -1 6 +1 3 S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 + 1 0.2976000 1.0000000 P 1 - 1 1.4070000 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 + 1 1.2750000 1.0000000 diff --git a/input/molecule b/input/molecule index 81c624a..c78e87e 100644 --- a/input/molecule +++ b/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. - H 0. 0. 1.4 + He 0.0 0.0 0.0 diff --git a/input/molecule.xyz b/input/molecule.xyz index 6edc99d..797b5fc 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,4 +1,3 @@ - 2 + 1 - H 0.0000000000 0.0000000000 0.0000000000 - H 0.0000000000 0.0000000000 0.7408481486 + He 0.0000000000 0.0000000000 0.0000000000 diff --git a/input/weight b/input/weight index 69b24c8..6796e3b 100644 --- a/input/weight +++ b/input/weight @@ -1,30 +1,9 @@ -1 6 +1 3 S 3 - 1 33.8700000 0.0060680 - 2 5.0950000 0.0453080 - 3 1.1590000 0.2028220 + 1 38.3600000 0.0238090 + 2 5.7700000 0.1548910 + 3 1.2400000 0.4699870 S 1 - 1 0.3258000 1.0000000 -S 1 - 1 0.1027000 1.0000000 + 1 0.2976000 1.0000000 P 1 - 1 1.4070000 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 + 1 1.2750000 1.0000000 diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index 8f35f2a..f7a1ece 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/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' diff --git a/src/QuAcK/orthogonalization_matrix.f90 b/src/QuAcK/orthogonalization_matrix.f90 index 15ea4ac..84b7f00 100644 --- a/src/QuAcK/orthogonalization_matrix.f90 +++ b/src/QuAcK/orthogonalization_matrix.f90 @@ -34,9 +34,9 @@ subroutine orthogonalization_matrix(ortho_type,nBas,S,X) if(ortho_type == 1) then - write(*,*) - write(*,*) ' Lowdin orthogonalization' - write(*,*) +! write(*,*) +! write(*,*) ' Lowdin orthogonalization' +! write(*,*) Uvec = S call diagonalize_matrix(nBas,Uvec,Uval) diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 index a59df23..1e10b96 100644 --- a/src/QuAcK/sort_ppRPA.f90 +++ b/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), & - O1(nVV,nVV),O2(nOO,nOO)) + order1(nVV),order2(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)) - S2 = - matmul(transpose(Z2),matmul(M,Z2)) + ! Find degenerate eigenvalues - if(nVV > 0) call orthogonalization_matrix(1,nVV,S1,O1) - if(nOO > 0) call orthogonalization_matrix(1,nOO,S2,O2) + deg1 = 1 + ab_start = 1 - Z1 = matmul(Z1,O1) - Z2 = matmul(Z2,O2) + do ab=1,nVV + + 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 diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 724e248..2a38dfa 100644 --- a/src/utils/wrap_lapack.f90 +++ b/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)