From 7d63abee672c0e2ebf06f746beb9360f31832c0f Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 21 Jan 2025 15:01:59 +0100 Subject: [PATCH 1/2] rCCD for Petros --- src/CC/EE_EOM_CCD_1h1p.f90 | 14 +++++++--- src/CC/rCCD.f90 | 3 +-- src/utils/sort.f90 | 32 +++++++++++++++++++++++ src/utils/wrap_lapack.f90 | 52 +++++++++++++++++++++++++++++++++++++- 4 files changed, 94 insertions(+), 7 deletions(-) diff --git a/src/CC/EE_EOM_CCD_1h1p.f90 b/src/CC/EE_EOM_CCD_1h1p.f90 index d50ed49..f93c1c1 100644 --- a/src/CC/EE_EOM_CCD_1h1p.f90 +++ b/src/CC/EE_EOM_CCD_1h1p.f90 @@ -28,7 +28,8 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) double precision,allocatable :: Wovvo(:,:,:,:) double precision,allocatable :: H(:,:) double precision,allocatable :: Om(:) - double precision,allocatable :: Z(:,:) + double precision,allocatable :: VL(:,:) + double precision,allocatable :: VR(:,:) integer,allocatable :: order(:) @@ -46,7 +47,7 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) ! Memory allocation - allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS),Z(nS,nS)) + allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS),VL(nS,nS),VR(nS,nS)) allocate(order(nS)) @@ -134,17 +135,22 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) if(nS > 0) then - call diagonalize_general_matrix(nS,H,Om,Z) + call diagonalize_general_matrix_LR(nS,H,Om,VL,VR) do ia=1,nS order(ia) = ia end do call quick_sort(Om,order,nS) - call set_order(Z,order,nS,nS) + call set_order_LR(VL,VR,order,nS,nS) call print_excitation_energies('EE-EOM-CCD','spinorbital',nS,Om) + write(*,*) 'Right Eigenvectors' + call matout(nS,nS,VR) + + call matout(nS,3,VR(:,1:3)) + end if end subroutine diff --git a/src/CC/rCCD.f90 b/src/CC/rCCD.f90 index aa4c80b..debb801 100644 --- a/src/CC/rCCD.f90 +++ b/src/CC/rCCD.f90 @@ -220,8 +220,7 @@ subroutine rCCD(dotest,maxSCF,thresh,max_diis,nBasin,nCin,nOin,nVin,nRin,ERI,ENu ! EE-EOM-CCD (1h1p) -! if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) - + if(do_EE_EOM_CC_1h1p) call EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) if(dotest) then diff --git a/src/utils/sort.f90 b/src/utils/sort.f90 index 36793ae..93ce1a3 100644 --- a/src/utils/sort.f90 +++ b/src/utils/sort.f90 @@ -60,6 +60,38 @@ recursive subroutine rec_quicksort(x,iorder,isize,first,last,level) end if end subroutine +subroutine set_order_LR(x,y,iorder,isize,jsize) + + implicit none + + integer :: isize,jsize + double precision :: x(isize,jsize) + double precision :: y(isize,jsize) + double precision,allocatable :: xtmp(:,:) + double precision,allocatable :: ytmp(:,:) + integer :: iorder(*) + integer :: i,j + + allocate(xtmp(isize,jsize),ytmp(isize,jsize)) + + do i=1,isize + do j=1,jsize + xtmp(i,j) = x(i,iorder(j)) + ytmp(i,j) = y(i,iorder(j)) + end do + end do + + do i=1,isize + do j=1,jsize + x(i,j) = xtmp(i,j) + y(i,j) = ytmp(i,j) + end do + end do + + deallocate(xtmp,ytmp) + +end subroutine + subroutine set_order(x,iorder,isize,jsize) implicit none diff --git a/src/utils/wrap_lapack.f90 b/src/utils/wrap_lapack.f90 index 86280fe..9cb4fff 100644 --- a/src/utils/wrap_lapack.f90 +++ b/src/utils/wrap_lapack.f90 @@ -1,3 +1,53 @@ +subroutine diagonalize_general_matrix_LR(N,A,WR,VL,VR) + +! Diagonalize a non-symmetric square matrix + + implicit none + +! Input variables + + integer,intent(in) :: N + double precision,intent(inout):: A(N,N) + double precision,intent(out) :: VL(N,N) + double precision,intent(out) :: VR(N,N) + double precision,intent(out) :: WR(N) + +! Local variables + + integer :: i + double precision :: tmp + integer :: lwork,info + double precision,allocatable :: work(:),WI(:) + +! Memory allocation + + allocate(work(1),WI(N)) + + lwork = -1 + call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info) + lwork = int(work(1)) + + deallocate(work) + + allocate(work(lwork)) + + call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info) + + do i=1,N + tmp = dot_product(vl(:,i),vr(:,i)) + vl(:,i) = vl(:,i)/tmp + end do + + call matout(N,N,matmul(transpose(VL),VR)) + + deallocate(work,WI) + + if(info /= 0) then + print*,'Problem in diagonalize_general_matrix (dgeev)!!' + end if + +end subroutine + subroutine diagonalize_general_matrix(N,A,WR,VR) ! Diagonalize a non-symmetric square matrix @@ -31,7 +81,7 @@ subroutine diagonalize_general_matrix(N,A,WR,VR) call dgeev('V','V',N,A,N,WR,WI,VL,N,VR,N,work,lwork,info) - deallocate(work, WI, VL) + deallocate(work,WI,VL) if(info /= 0) then print*,'Problem in diagonalize_general_matrix (dgeev)!!' From 251ff8fdfe3b08d9b3f2c5372fdf2cffa5230aa9 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Tue, 21 Jan 2025 16:31:13 +0100 Subject: [PATCH 2/2] debug petros --- src/CC/EE_EOM_CCD_1h1p.f90 | 138 +++++++++++++++++++++++++++++++++++-- 1 file changed, 133 insertions(+), 5 deletions(-) diff --git a/src/CC/EE_EOM_CCD_1h1p.f90 b/src/CC/EE_EOM_CCD_1h1p.f90 index f93c1c1..ef3f7f9 100644 --- a/src/CC/EE_EOM_CCD_1h1p.f90 +++ b/src/CC/EE_EOM_CCD_1h1p.f90 @@ -30,9 +30,20 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) double precision,allocatable :: Om(:) double precision,allocatable :: VL(:,:) double precision,allocatable :: VR(:,:) + double precision,allocatable :: Leom(:,:,:) + double precision,allocatable :: Reom(:,:,:) + + integer :: nstate,m + double precision :: Ex,tmp integer,allocatable :: order(:) + double precision,allocatable :: rdm1_oo(:,:) + double precision,allocatable :: rdm1_vv(:,:) + + double precision,allocatable :: rdm2_oovv(:,:,:,:) + double precision,allocatable :: rdm2_ovvo(:,:,:,:) + ! Hello world write(*,*) @@ -47,10 +58,9 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) ! Memory allocation - allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS),VL(nS,nS),VR(nS,nS)) + allocate(Foo(nO,nO),Fvv(nV,nV),Wovvo(nO,nV,nV,nO),H(nS,nS),Om(nS)) allocate(order(nS)) - ! Form one-body terms do a=1,nV-nR @@ -133,6 +143,8 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) ! Diagonalize EOM Hamiltonian + allocate(VL(nS,nS),VR(nS,nS)) + if(nS > 0) then call diagonalize_general_matrix_LR(nS,H,Om,VL,VR) @@ -146,11 +158,127 @@ subroutine EE_EOM_CCD_1h1p(nC,nO,nV,nR,eO,eV,OOVV,OVVO,t) call print_excitation_energies('EE-EOM-CCD','spinorbital',nS,Om) - write(*,*) 'Right Eigenvectors' - call matout(nS,nS,VR) +! write(*,*) 'Right Eigenvectors' +! call matout(nS,nS,VR) - call matout(nS,3,VR(:,1:3)) +! call matout(nS,3,VR(:,1:3)) end if + allocate(Leom(nO,nV,nS),Reom(nO,nV,nS)) + + do m=1,nS + ia = 0 + do i=1,nO + do a=1,nV + ia = ia + 1 + Leom(i,a,m) = VL(ia,m) + Reom(i,a,m) = VR(ia,m) + end do + end do + end do + + deallocate(VL,VR) + +!------------------------------------------------------------------------ +! EOM section +!------------------------------------------------------------------------ + + allocate(rdm1_oo(nO,nO),rdm1_vv(nV,nV)) + allocate(rdm2_oovv(nO,nO,nV,nV),rdm2_ovvo(nO,nV,nV,nO)) + + nstate = 1 + + tmp = 0d0 + do i=1,nO + do a=1,nV + tmp = tmp + Leom(i,a,nstate)*Reom(i,a,nstate) + end do + end do + print*,tmp + + rdm1_oo(:,:) = 0d0 + do i=1,nO + do j=1,nO + do c=1,nV + + rdm1_oo(i,j) = rdm1_oo(i,j) - Reom(i,c,nstate)*Leom(j,c,nstate) + + end do + end do + end do + + rdm1_vv(:,:) = 0d0 + do a=1,nV + do b=1,nV + do k=1,nO + + rdm1_vv(a,b) = rdm1_vv(a,b) + Reom(k,b,nstate)*Leom(k,a,nstate) + + end do + end do + end do + + rdm2_ovvo(:,:,:,:) = 0d0 + do i=1,nO + do a=1,nV + do b=1,nV + do j=1,nO + + rdm2_ovvo(i,a,b,j) = Reom(i,b,nstate)*Leom(j,a,nstate) + + end do + end do + end do + end do + + rdm2_oovv(:,:,:,:) = 0d0 + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + do k=1,nO + do c=1,nV + + rdm2_oovv(i,j,a,b) = rdm2_oovv(i,j,a,b) & + + Reom(j,b,nstate)*t(k,i,c,a)*Leom(k,c,nstate) & + - Reom(i,b,nstate)*t(k,j,c,a)*Leom(k,c,nstate) & + - Reom(j,a,nstate)*t(k,i,c,b)*Leom(k,c,nstate) & + + Reom(i,a,nstate)*t(k,j,c,b)*Leom(k,c,nstate) + + end do + end do + + end do + end do + end do + end do + + Ex = 0d0 + + do i=1,nO + Ex = Ex + rdm1_oo(i,i)*eO(i) + end do + + do a=1,nV + Ex = Ex + rdm1_vv(a,a)*eV(a) + end do + + do i=1,nO + do a=1,nV + do b=1,nV + do j=1,nO + + Ex = Ex + rdm2_ovvo(i,a,b,j)*OVVO(i,a,b,j) + 0.25d0*rdm2_oovv(i,j,a,b)*OOVV(i,j,a,b) + + end do + end do + end do + end do + + print*,'Ex = ',Ex + print*,'Om = ',Om(nstate) + + end subroutine