From 9910e6e88924b6f788171dff1386b3132a5d6f1d Mon Sep 17 00:00:00 2001 From: pfloos Date: Thu, 29 Aug 2024 19:21:25 +0200 Subject: [PATCH 1/3] working on pCCD --- src/CC/pCCD.f90 | 79 ++++++++++++++++++++++++++++++++++++++++++++----- 1 file changed, 71 insertions(+), 8 deletions(-) diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index 4f8b67d..570a679 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -57,10 +57,13 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, double precision :: E1,E2 double precision,allocatable :: h(:,:) + double precision,allocatable :: c(:,:) double precision,allocatable :: grad(:) double precision,allocatable :: tmp(:,:,:,:) double precision,allocatable :: hess(:,:) - double precision,allocatable :: eig(:) + double precision,allocatable :: hessInv(:,:) + double precision,allocatable :: kappa(:,:) + double precision,allocatable :: ekappa(:,:) integer :: O,V,N integer :: n_diis @@ -126,6 +129,13 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis)) +!------------------------------------------------------------------------ +! Start orbital optimization +!------------------------------------------------------------------------ + + allocate(c(N,N)) + c(:,:) = cHF(:,:) + !------------------------------------------------------------------------ ! Compute t ampltiudes !------------------------------------------------------------------------ @@ -540,6 +550,13 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, write(*,*) 'Orbital gradient at the pCCD level:' call matout(N,N,grad) + write(*,*) + +! Convergence + + Conv = maxval(abs(grad)) + write(*,*) ' Convergence of orbtial gradient = ',Conv + write(*,*) ! Compute Hessian @@ -550,7 +567,6 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, do p=1,N do q=1,N - rs = 0 do r=1,N do s=1,N @@ -614,7 +630,8 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, rs = rs + 1 - hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r) + hess(pq,rs) = tmp(p,r,q,s) - tmp(r,p,q,s) - tmp(p,r,s,q) + tmp(r,p,s,q) +! hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r) end do end do @@ -624,13 +641,59 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, call matout(N**2,N**2,hess) - deallocate(tmp) +! allocate(eig(N**2)) - allocate(eig(N**2)) - - call diagonalize_matrix(N**2,hess,eig) +! call diagonalize_matrix(N**2,hess,eig) - call vecout(N**2,eig) +! call vecout(N**2,eig) + + allocate(hessInv(N**2,N**2)) + + call inverse_matrix(N**2,hess,hessInv) + + allocate(kappa(N,N),ekappa(N,N)) + + kappa(:,:) = 0d0 + + pq = 0 + do p=1,N + do q=1,N + + pq = pq + 1 + + rs = 0 + do r=1,N + do s=1,N + + rs = rs + 1 + + kappa(p,q) = kappa(p,q) + hessInv(pq,rs)*grad(rs) + + end do + end do + + end do + end do + + write(*,*) 'kappa' + call matout(N,N,kappa) + write(*,*) + + call matrix_exponential(N,kappa,ekappa) + + write(*,*) 'e^kappa' + call matout(N,N,ekappa) + write(*,*) + + write(*,*) 'Old orbitals' + call matout(N,N,c) + write(*,*) + + c = matmul(c,ekappa) + + write(*,*) 'Rotated orbitals' + call matout(N,N,c) + write(*,*) ! Testing zone From 08bf6632dfbdd65f5aa18b9c1e8d67cf4d7d8049 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Fri, 30 Aug 2024 16:26:59 +0200 Subject: [PATCH 2/3] saving work in pCCD --- src/CC/RCC.f90 | 22 +- src/CC/pCCD.f90 | 1141 ++++++++++++++++++++++-------------------- src/QuAcK/RQuAcK.f90 | 2 +- 3 files changed, 603 insertions(+), 562 deletions(-) diff --git a/src/CC/RCC.f90 b/src/CC/RCC.f90 index 21f6866..aa8923b 100644 --- a/src/CC/RCC.f90 +++ b/src/CC/RCC.f90 @@ -1,5 +1,5 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF) + maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF) ! Coupled-cluster module @@ -34,7 +34,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_MO(nBas,nBas,nBas,nBas) ! Local variables @@ -47,7 +48,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCD) then call wall_time(start_CC) - call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call CCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -64,7 +65,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d call wall_time(start_CC) call DCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR, & - ERI,ENuc,ERHF,eHF) + ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -82,7 +83,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(doCCSD) then call wall_time(start_CC) - call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call CCSD(dotest,maxSCF,thresh,max_diis,doCCSDT,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -98,7 +99,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dodrCCD) then call wall_time(start_CC) - call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call drCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -114,7 +115,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dorCCD) then call wall_time(start_CC) - call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call rCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -130,7 +131,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(docrCCD) then call wall_time(start_CC) - call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call crCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -146,7 +147,7 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dolCCD) then call wall_time(start_CC) - call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI,ENuc,ERHF,eHF) + call lCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,ERI_MO,ENuc,ERHF,eHF) call wall_time(end_CC) t_CC = end_CC - start_CC @@ -162,7 +163,8 @@ subroutine RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,d if(dopCCD) then call wall_time(start_CC) - call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF) + call pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF) + call wall_time(end_CC) t_CC = end_CC - start_CC diff --git a/src/CC/pCCD.f90 b/src/CC/pCCD.f90 index 570a679..bcfd50b 100644 --- a/src/CC/pCCD.f90 +++ b/src/CC/pCCD.f90 @@ -1,4 +1,4 @@ -subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF,eHF,cHF) +subroutine pCCD(dotest,maxIt,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI_AO,ENuc,ERHF,eHF,cHF) ! pair CCD module @@ -8,25 +8,32 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, logical,intent(in) :: dotest - integer,intent(in) :: maxSCF + integer,intent(in) :: maxIt integer,intent(in) :: max_diis double precision,intent(in) :: thresh - integer,intent(in) :: nBas,nC,nO,nV,nR - double precision,intent(in) :: ENuc,ERHF + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: cHF(nBas,nBas) double precision,intent(in) :: Hc(nBas,nBas) - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ERI_AO(nBas,nBas,nBas,nBas) ! Local variables - integer :: p,q,r,s,t,u + integer :: p,q,r,s,t,u,w integer :: pq,rs integer :: i,j,a,b - integer :: nSCF - double precision :: Conv + integer :: nItAmp + integer :: nItOrb + double precision :: CvgAmp + double precision :: CvgOrb double precision :: ECC double precision :: EcCC @@ -56,14 +63,15 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, double precision :: tr_2rdm double precision :: E1,E2 - double precision,allocatable :: h(:,:) double precision,allocatable :: c(:,:) + double precision,allocatable :: h(:,:) + double precision,allocatable :: ERI_MO(:,:,:,:) double precision,allocatable :: grad(:) double precision,allocatable :: tmp(:,:,:,:) double precision,allocatable :: hess(:,:) double precision,allocatable :: hessInv(:,:) - double precision,allocatable :: kappa(:,:) - double precision,allocatable :: ekappa(:,:) + double precision,allocatable :: Kap(:,:) + double precision,allocatable :: ExpKap(:,:) integer :: O,V,N integer :: n_diis @@ -77,9 +85,9 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, ! Hello world write(*,*) - write(*,*)'**************************************' - write(*,*)'| pair CCD calculation |' - write(*,*)'**************************************' + write(*,*)'*******************************' + write(*,*)'* Restricted pCCD Calculation *' + write(*,*)'*******************************' write(*,*) ! Useful quantities @@ -88,612 +96,643 @@ subroutine pCCD(dotest,maxSCF,thresh,max_diis,nBas,nC,nO,nV,nR,Hc,ERI,ENuc,ERHF, V = nV - nR N = O + V -! Form energy denominator + !------------------------------------! + ! Star Loop for orbital optimization ! + !------------------------------------! + allocate(ERI_MO(N,N,N,N)) + allocate(c(N,N),h(N,N)) allocate(eO(O),eV(V),delta_OV(O,V)) - - eO(:) = eHF(nC+1:nO) - eV(:) = eHF(nO+1:nBas-nR) - - call form_delta_OV(nC,nO,nV,nR,eO,eV,delta_OV) - -! Create integral batches - allocate(OOOO(O,O),OOVV(O,V),OVOV(O,V),OVVO(O,V),VVVV(V,V)) - do i=1,O - do j=1,O - OOOO(i,j) = ERI(nC+i,nC+i,nC+j,nC+j) - end do - end do + c(:,:) = cHF(nC+1:nBas-nR,nC+1:nBas-nR) - do i=1,O - do a=1,V - OOVV(i,a) = ERI(nC+i,nC+i,nO+a,nO+a) - OVOV(i,a) = ERI(nC+i,nO+a,nC+i,nO+a) - OVVO(i,a) = ERI(nC+i,nO+a,nO+a,nC+i) - end do - end do + CvgOrb = 1d0 + nItOrb = 0 - do a=1,V - do b=1,V - VVVV(a,b) = ERI(nO+a,nO+a,nO+b,nO+b) - end do - end do - -! Initialization - - allocate(t2(O,V),r2(O,V),yO(O,O),yV(V,V)) - -! Memory allocation for DIIS - - allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis)) - -!------------------------------------------------------------------------ -! Start orbital optimization -!------------------------------------------------------------------------ - - allocate(c(N,N)) - c(:,:) = cHF(:,:) - -!------------------------------------------------------------------------ -! Compute t ampltiudes -!------------------------------------------------------------------------ - - Conv = 1d0 - nSCF = 0 - ECC = ERHF - EcCC = 0d0 - - n_diis = 0 - t2(:,:) = 0d0 - t2_diis(:,:) = 0d0 - err_diis(:,:) = 0d0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ write(*,*) write(*,*)'----------------------------------------------------' - write(*,*)'| pCCD calculation: t amplitudes |' - write(*,*)'----------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|' + write(*,*)'| Orbital Optimization for pCCD |' write(*,*)'----------------------------------------------------' - do while(Conv > thresh .and. nSCF < maxSCF) + do while(CvgOrb > thresh .and. nItOrb < 1) - ! Increment + nItOrb = nItOrb + 1 - nSCF = nSCF + 1 + ! Transform integrals - ! Form intermediate array - - yO(:,:) = matmul(t2,transpose(OOVV)) - - ! Compute residual + h = matmul(transpose(c),matmul(Hc(nC+1:nBas-nR,nC+1:nBas-nR),c)) + call AOtoMO_ERI_RHF(N,c,ERI_AO(nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR,nC+1:nBas-nR),ERI_MO) + ! Form energy denominator - r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) & - - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:) + eO(:) = eHF(nC+1:nO) + eV(:) = eHF(nO+1:nBas-nR) do i=1,O do a=1,V - - do j=1,O - r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*t2(i,a) + OOOO(j,i)*t2(j,a) + yO(i,j)*t2(j,a) - end do - - do b=1,V - r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*t2(i,a) + VVVV(a,b)*t2(i,b) - end do - + delta_OV(i,a) = eV(a) - eO(i) end do end do - ! Check convergence - - Conv = maxval(abs(r2(:,:))) - - ! Update amplitudes - - t2(:,:) = t2(:,:) - 0.5d0*r2(:,:)/delta_OV(:,:) - - ! Compute correlation energy - - EcCC = trace_matrix(V,matmul(transpose(OOVV),t2)) - - ! Dump results - - ECC = ERHF + EcCC - - ! DIIS extrapolation - - if(max_diis > 1) then - - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,err_diis,t2_diis,-0.5d0*r2/delta_OV,t2) - - end if - - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - - end do - write(*,*)'----------------------------------------------------' -!------------------------------------------------------------------------ -! End of SCF loop -!------------------------------------------------------------------------ - -! Did it actually converge? - - if(nSCF == maxSCF) then - - write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed for t ampitudes ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - - stop - - end if - -! Deallocate memory - - deallocate(err_diis,t2_diis) - -! Memory allocation - - allocate(z2(O,V)) - -! Memory allocation for DIIS - - allocate(err_diis(O*V,max_diis),z2_diis(O*V,max_diis)) - -!------------------------------------------------------------------------ -! Compute z ampltiudes -!------------------------------------------------------------------------ - - Conv = 1d0 - nSCF = 0 - - n_diis = 0 - z2_diis(:,:) = 0d0 - err_diis(:,:) = 0d0 - -!------------------------------------------------------------------------ -! Main SCF loop -!------------------------------------------------------------------------ - write(*,*) - write(*,*)'----------------------------------------------------' - write(*,*)'| pCCD calculation: z amplitudes |' - write(*,*)'----------------------------------------------------' - write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & - '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|' - write(*,*)'----------------------------------------------------' - - do while(Conv > thresh .and. nSCF < maxSCF) - - ! Increment - - nSCF = nSCF + 1 - - ! Form intermediate array - - yO(:,:) = matmul(OOVV,transpose(t2)) - yV(:,:) = matmul(transpose(OOVV),t2) - - ! Compute residual - - r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*z2(:,:) & - - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - 2d0*OOVV(:,:)*t2(:,:))*z2(:,:) + ! Create integral batches + do i=1,O + do j=1,O + OOOO(i,j) = ERI_MO(i,i,j,j) + end do + end do + do i=1,O do a=1,V - - do j=1,O - r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*z2(i,a) - 2d0*OOVV(i,a)*z2(j,a)*t2(j,a) & - + OOOO(i,j)*z2(j,a) + yO(i,j)*z2(j,a) - end do - - do b=1,V - r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*z2(i,a) - 2d0*OOVV(i,a)*z2(i,b)*t2(i,b) & - + VVVV(b,a)*z2(i,b) + yV(a,b)*z2(i,b) - end do - + OOVV(i,a) = ERI_MO(i,i,O+a,O+a) + OVOV(i,a) = ERI_MO(i,O+a,i,O+a) + OVVO(i,a) = ERI_MO(i,O+a,O+a,i) + end do + end do + + do a=1,V + do b=1,V + VVVV(a,b) = ERI_MO(O+a,O+a,O+b,O+b) end do end do - ! Check convergence + !----------------------------! + ! Star Loop for t amplitudes ! + !----------------------------! + + allocate(t2(O,V),r2(O,V),yO(O,O)) + allocate(err_diis(O*V,max_diis),t2_diis(O*V,max_diis)) - Conv = maxval(abs(r2(:,:))) + CvgAmp = 1d0 + nItAmp = 0 + ECC = ERHF + EcCC = 0d0 + + n_diis = 0 + t2(:,:) = 0d0 + t2_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 + + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| pCCD calculation: t amplitudes |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(CvgAmp > thresh .and. nItAmp < maxIt) + + ! Increment + + nItAmp = nItAmp + 1 + + ! Form intermediate array + + yO(:,:) = matmul(t2,transpose(OOVV)) + + ! Compute residual + + r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*t2(:,:) & + - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - OOVV(:,:)*t2(:,:))*t2(:,:) + + do i=1,O + do a=1,V + + do j=1,O + r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*t2(i,a) + OOOO(j,i)*t2(j,a) + yO(i,j)*t2(j,a) + end do + + do b=1,V + r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*t2(i,a) + VVVV(a,b)*t2(i,b) + end do + + end do + end do + + ! Check convergence + + CvgAmp = maxval(abs(r2(:,:))) + + ! Update amplitudes + + t2(:,:) = t2(:,:) - 0.5d0*r2(:,:)/delta_OV(:,:) + + ! Compute correlation energy + + EcCC = 0d0 + do i=1,O + do a=1,V + EcCC = EcCC + OOVV(i,a)*t2(i,a) + end do + end do + + ! Dump results + + ECC = ERHF + EcCC + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nO*nV,nO*nV,n_diis,err_diis,t2_diis,-0.5d0*r2/delta_OV,t2) + + end if + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|' + + end do + write(*,*)'----------------------------------------------------' + + !---------------------------! + ! End Loop for t amplitudes ! + !---------------------------! + + deallocate(r2,yO) + deallocate(err_diis,t2_diis) + + ! Did it actually converge? + + if(nItAmp == maxIt) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)'! Convergence failed for t ampitudes !' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + + stop - ! Update amplitudes - - z2(:,:) = z2(:,:) - 0.5d0*r2(:,:)/delta_OV(:,:) - - ! DIIS extrapolation - - if(max_diis > 1) then - - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,O*V,O*V,n_diis,err_diis,z2_diis,-0.5d0*r2/delta_OV,z2) - end if + + !-----------------------------! + ! Start Loop for z amplitudes ! + !-----------------------------! + + allocate(z2(O,V),r2(O,V),yO(O,O),yV(V,V)) + allocate(err_diis(O*V,max_diis),z2_diis(O*V,max_diis)) - write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & - '|',nSCF,'|',ECC+ENuc,'|',EcCC,'|',Conv,'|' - - end do - write(*,*)'----------------------------------------------------' - write(*,*) -!------------------------------------------------------------------------ -! End of SCF loop -!------------------------------------------------------------------------ - -! Did it actually converge? - - if(nSCF == maxSCF) then + CvgAmp = 1d0 + nItAmp = 0 + + n_diis = 0 + z2_diis(:,:) = 0d0 + err_diis(:,:) = 0d0 write(*,*) - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - write(*,*)' Convergence failed ' - write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - - stop + write(*,*)'----------------------------------------------------' + write(*,*)'| pCCD calculation: z amplitudes |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(pCCD)','|','Ec(pCCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(CvgAmp > thresh .and. nItAmp < maxIt) + + ! Increment + + nItAmp = nItAmp + 1 + + ! Form intermediate array + + yO(:,:) = matmul(OOVV,transpose(t2)) + yV(:,:) = matmul(transpose(OOVV),t2) + + ! Compute residual + + r2(:,:) = OOVV(:,:) + 2d0*delta_OV(:,:)*z2(:,:) & + - 2d0*(2d0*OVOV(:,:) - OVVO(:,:) - 2d0*OOVV(:,:)*t2(:,:))*z2(:,:) + + do i=1,O + do a=1,V + + do j=1,O + r2(i,a) = r2(i,a) - 2d0*OOVV(j,a)*t2(j,a)*z2(i,a) - 2d0*OOVV(i,a)*z2(j,a)*t2(j,a) & + + OOOO(i,j)*z2(j,a) + yO(i,j)*z2(j,a) + end do + + do b=1,V + r2(i,a) = r2(i,a) - 2d0*OOVV(i,b)*t2(i,b)*z2(i,a) - 2d0*OOVV(i,a)*z2(i,b)*t2(i,b) & + + VVVV(b,a)*z2(i,b) + yV(a,b)*z2(i,b) + end do + + end do + end do + + ! Check convergence + + CvgAmp = maxval(abs(r2(:,:))) + + ! Update amplitudes + + z2(:,:) = z2(:,:) - 0.5d0*r2(:,:)/delta_OV(:,:) + + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,O*V,O*V,n_diis,err_diis,z2_diis,-0.5d0*r2/delta_OV,z2) + + end if + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nItAmp,'|',ECC+ENuc,'|',EcCC,'|',CvgAmp,'|' - end if - -! Deallocate memory - - deallocate(err_diis,z2_diis,r2) - -!--------------------------! -! Compute density matrices ! -!--------------------------! - - allocate(xOO(O,O),xVV(V,V),xOV(O,V)) - - xOO(:,:) = matmul(t2,transpose(z2)) - xVV(:,:) = matmul(transpose(z2),t2) - xOV(:,:) = matmul(t2,matmul(transpose(z2),t2)) - -! Form 1RDM - - allocate(rdm1(N,N)) - - rdm1(:,:) = 0d0 - - do i=1,O - rdm1(i,i) = 2d0*(1d0 - xOO(i,i)) - end do - - do a=1,V - rdm1(O+a,O+a) = 2d0*xVV(a,a) - end do - -! Check 1RDM - - tr_1rdm = trace_matrix(N,rdm1) - write(*,'(A25,F16.10)') ' --> Trace of the 1RDM = ',tr_1rdm - - if( abs(dble(2*O) - tr_1rdm) > thresh ) & - write(*,*) ' !!! Your 1RDM seems broken !!! ' - write(*,*) - -! write(*,*) '1RDM is diagonal at the pCCD level:' -! call matout(N,N,rdm1) - -! Form 2RM - - allocate(rdm2(N,N,N,N)) - - rdm2(:,:,:,:) = 0d0 - - ! iijj - - do i=1,O - do j=1,O - rdm2(i,i,j,j) = 2d0*xOO(i,j) end do - end do + write(*,*)'----------------------------------------------------' + write(*,*) - ! iiaa + !---------------------------! + ! End Loop for z ampltiudes ! + !---------------------------! - do i=1,O + deallocate(r2,yO,yV) + deallocate(err_diis,z2_diis) + + ! Did it actually converge? + + if(nItAmp == maxIt) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)'! Convergence failed for z ampltiudes !' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + + stop + + end if + + !--------------------------! + ! Compute density matrices ! + !--------------------------! + + allocate(rdm1(N,N),rdm2(N,N,N,N)) + allocate(xOO(O,O),xVV(V,V),xOV(O,V)) + + xOO(:,:) = matmul(t2,transpose(z2)) + xVV(:,:) = matmul(transpose(z2),t2) + xOV(:,:) = matmul(t2,matmul(transpose(z2),t2)) + + ! Form 1RDM + + rdm1(:,:) = 0d0 + + do i=1,O + rdm1(i,i) = 2d0*(1d0 - xOO(i,i)) + end do + do a=1,V - rdm2(i,i,O+a,O+a) = 2d0*(t2(i,a) + xOV(i,a) - 2d0*t2(i,a)*(xVV(a,a) + xOO(i,i) - t2(i,a)*z2(i,a))) + rdm1(O+a,O+a) = 2d0*xVV(a,a) end do - end do - - ! aaii - - do i=1,O + + ! Check 1RDM + + tr_1rdm = trace_matrix(N,rdm1) + write(*,'(A25,F16.10)') ' --> Trace of the 1RDM = ',tr_1rdm + + if( abs(dble(2*O) - tr_1rdm) > thresh ) & + write(*,*) ' !!! Your 1RDM seems broken !!! ' + write(*,*) + + ! write(*,*) '1RDM is diagonal at the pCCD level:' + ! call matout(N,N,rdm1) + + ! Form 2RM + + rdm2(:,:,:,:) = 0d0 + + ! iijj + + do i=1,O + do j=1,O + rdm2(i,i,j,j) = 2d0*xOO(i,j) + end do + end do + + ! iiaa + + do i=1,O + do a=1,V + rdm2(i,i,O+a,O+a) = 2d0*(t2(i,a) + xOV(i,a) - 2d0*t2(i,a)*(xVV(a,a) + xOO(i,i) - t2(i,a)*z2(i,a))) + end do + end do + + ! aaii + + do i=1,O + do a=1,V + rdm2(O+a,O+a,i,i) = 2d0*z2(i,a) + end do + end do + + ! aabb + do a=1,V - rdm2(O+a,O+a,i,i) = 2d0*z2(i,a) + do b=1,V + rdm2(O+a,O+a,O+b,O+b) = 2d0*xVV(a,b) + end do end do - end do - - ! aabb - - do a=1,V - do b=1,V - rdm2(O+a,O+a,O+b,O+b) = 2d0*xVV(a,b) + + ! ijij + + do i=1,O + do j=1,O + rdm2(i,j,i,j) = 4d0*(1d0 - xOO(i,i) - xOO(j,j)) + end do end do - end do - - ! ijij - - do i=1,O - do j=1,O - rdm2(i,j,i,j) = 4d0*(1d0 - xOO(i,i) - xOO(j,j)) + + ! ijji + + do i=1,O + do j=1,O + rdm2(i,j,j,i) = - 2d0*(1d0 - xOO(i,i) - xOO(j,j)) + end do end do - end do - - ! ijji - - do i=1,O - do j=1,O - rdm2(i,j,j,i) = - 2d0*(1d0 - xOO(i,i) - xOO(j,j)) + + ! iiii + + do i=1,O + rdm2(i,i,i,i) = 2d0*(1d0 - xOO(i,i)) end do - end do - - ! iiii - - do i=1,O - rdm2(i,i,i,i) = 2d0*(1d0 - xOO(i,i)) - end do - - ! iaia - - do i=1,O + + ! iaia + + do i=1,O + do a=1,V + rdm2(i,O+a,i,O+a) = 4d0*(xVV(a,a) - t2(i,a)*z2(i,a)) + end do + end do + + ! iaai + + do i=1,O + do a=1,V + rdm2(i,O+a,O+a,i) = - 2d0*(xVV(a,a) - t2(i,a)*z2(i,a)) + end do + end do + + ! aiai + + do i=1,O + do a=1,V + rdm2(O+a,i,O+a,i) = 4d0*(xVV(a,a) - t2(i,a)*z2(i,a)) + end do + end do + + ! aiia + + do i=1,O + do a=1,V + rdm2(O+a,i,i,O+a) = - 2d0*(xVV(a,a) - t2(i,a)*z2(i,a)) + end do + end do + + ! abab + do a=1,V - rdm2(i,O+a,i,O+a) = 4d0*(xVV(a,a) - t2(i,a)*z2(i,a)) + rdm2(O+a,O+a,O+a,O+a) = 2d0*xVV(a,a) end do - end do + + ! Check 2RDM + + tr_2rdm = trace_matrix(N**2,rdm2) + write(*,'(A25,F16.10)') ' --> Trace of the 2RDM = ',tr_2rdm + + if( abs(dble(2*O*(2*O-1)) - tr_2rdm) > thresh ) & + write(*,*) ' !!! Your 2RDM seems broken !!! ' + write(*,*) + + ! write(*,*) '2RDM is not diagonal at the pCCD level:' + ! call matout(N**2,N**2,rdm2) - ! iaai - - do i=1,O - do a=1,V - rdm2(i,O+a,O+a,i) = - 2d0*(xVV(a,a) - t2(i,a)*z2(i,a)) - end do - end do - - ! aiai - - do i=1,O - do a=1,V - rdm2(O+a,i,O+a,i) = 4d0*(xVV(a,a) - t2(i,a)*z2(i,a)) - end do - end do - - ! aiia - - do i=1,O - do a=1,V - rdm2(O+a,i,i,O+a) = - 2d0*(xVV(a,a) - t2(i,a)*z2(i,a)) - end do - end do - - ! abab - - do a=1,V - rdm2(O+a,O+a,O+a,O+a) = 2d0*xVV(a,a) - end do - -! Check 2RDM - - tr_2rdm = trace_matrix(N**2,rdm2) - write(*,'(A25,F16.10)') ' --> Trace of the 2RDM = ',tr_2rdm - - if( abs(dble(2*O*(2*O-1)) - tr_2rdm) > thresh ) & - write(*,*) ' !!! Your 2RDM seems broken !!! ' - write(*,*) - -! write(*,*) '2RDM is not diagonal at the pCCD level:' -! call matout(N**2,N**2,rdm2) - -! Compute electronic energy - - allocate(h(N,N)) - h = matmul(transpose(cHF),matmul(Hc,cHF)) - - E1 = 0d0 - E2 = 0d0 - - do p=1,N - do q=1,N - E1 = E1 + rdm1(p,q)*h(p,q) - do r=1,N - do s=1,N - E2 = E2 + rdm2(p,q,r,s)*ERI(p,q,r,s) + deallocate(xOO,xVV,xOV) + deallocate(t2,z2) + + ! Compute electronic energy + + E1 = 0d0 + E2 = 0d0 + + do p=1,N + do q=1,N + E1 = E1 + rdm1(p,q)*h(p,q) + do r=1,N + do s=1,N + E2 = E2 + rdm2(p,q,r,s)*ERI_MO(p,q,r,s) + end do end do end do end do - end do - - E2 = 0.5d0*E2 - - write(*,'(A25,F16.10)') ' One-electron energy = ',E1 - write(*,'(A25,F16.10)') ' Two-electron energy = ',E2 - write(*,'(A25,F16.10)') ' Electronic energy = ',E1 + E2 - write(*,'(A25,F16.10)') ' Total energy = ',E1 + E2 + ENuc - write(*,*) - -! Compute gradient - - allocate(grad(N**2)) - - grad(:) = 0d0 - - pq = 0 - do p=1,N - do q=1,N - - pq = pq + 1 - - do r=1,N - grad(pq) = grad(pq) + h(r,p)*rdm1(r,q) - h(q,r)*rdm1(p,r) - end do - - do r=1,N - do s=1,N - do t=1,N - grad(pq) = grad(pq) + (ERI(r,s,p,t)*rdm2(r,s,q,t) - ERI(q,t,r,s)*rdm2(p,t,r,s)) + + E2 = 0.5d0*E2 + + write(*,'(A25,F16.10)') ' One-electron energy = ',E1 + write(*,'(A25,F16.10)') ' Two-electron energy = ',E2 + write(*,'(A25,F16.10)') ' Electronic energy = ',E1 + E2 + write(*,'(A25,F16.10)') ' Total energy = ',E1 + E2 + ENuc + write(*,*) + + !--------------------------! + ! Compute orbital gradient ! + !--------------------------! + + allocate(grad(N**2)) + + grad(:) = 0d0 + + pq = 0 + do p=1,N + do q=1,N + + pq = pq + 1 + + do r=1,N + grad(pq) = grad(pq) + h(r,p)*rdm1(r,q) - h(q,r)*rdm1(p,r) + end do + + do r=1,N + do s=1,N + do t=1,N + grad(pq) = grad(pq) + (ERI_MO(r,s,p,t)*rdm2(r,s,q,t) - ERI_MO(q,t,r,s)*rdm2(p,t,r,s)) + end do end do end do - end do - - end do - end do - - write(*,*) 'Orbital gradient at the pCCD level:' - call matout(N,N,grad) - write(*,*) - -! Convergence - - Conv = maxval(abs(grad)) - write(*,*) ' Convergence of orbtial gradient = ',Conv - write(*,*) -! Compute Hessian - - allocate(hess(N**2,N**2),tmp(N,N,N,N)) - - tmp(:,:,:,:) = 0d0 - - do p=1,N - do q=1,N - - do r=1,N - do s=1,N - - tmp(p,q,r,s) = - h(s,p)*rdm1(r,q) - h(q,r)*rdm1(p,s) - - do u=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & - Kronecker_delta(q,r)*(h(u,p)*rdm1(u,s) + h(s,u)*rdm1(p,u)) & - + Kronecker_delta(p,s)*(h(u,r)*rdm1(u,q) + h(q,u)*rdm1(r,u)) ) - - end do - - do u=1,N - do v=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) + ERI(u,v,p,r)*rdm2(u,v,q,s) + ERI(q,s,u,v)*rdm2(p,r,u,v) - - end do - end do - - do t=1,N + end do + end do + + write(*,*) 'Orbital gradient at the pCCD level:' + call matout(N,N,grad) + write(*,*) + + ! Check convergence of orbital optimization + + CvgOrb = maxval(abs(grad)) + write(*,*) ' Iteration',nItOrb,'for pCCD orbital optimization' + write(*,*) ' Convergence of orbital gradient = ',CvgOrb + write(*,*) + + !-------------------------! + ! Compute orbital Hessian ! + !-------------------------! + + allocate(hess(N**2,N**2),tmp(N,N,N,N)) + + tmp(:,:,:,:) = 0d0 + + do p=1,N + do q=1,N + + do r=1,N + do s=1,N + + tmp(p,q,r,s) = - h(s,p)*rdm1(r,q) - h(q,r)*rdm1(p,s) + do u=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) - ( & - ERI(s,t,p,u)*rdm2(r,t,q,u) + ERI(t,s,p,u)*rdm2(t,r,q,u) & - + ERI(q,u,r,t)*rdm2(p,u,s,t) + ERI(q,u,t,r)*rdm2(p,u,t,s) ) - + + tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & + Kronecker_delta(q,r)*(h(u,p)*rdm1(u,s) + h(s,u)*rdm1(p,u)) & + + Kronecker_delta(p,s)*(h(u,r)*rdm1(u,q) + h(q,u)*rdm1(r,u)) ) + end do - end do - - do t=1,N + do u=1,N - do v=1,N - - tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & - Kronecker_delta(q,r)*(ERI(u,v,p,t)*rdm2(u,v,s,t) + ERI(s,t,u,v)*rdm2(p,t,u,v)) & - + Kronecker_delta(p,s)*(ERI(q,t,u,v)*rdm2(r,t,u,v) + ERI(u,v,r,t)*rdm2(u,v,q,t)) ) - + do w=1,N + + tmp(p,q,r,s) = tmp(p,q,r,s) + ERI_MO(u,w,p,r)*rdm2(u,w,q,s) + ERI_MO(q,s,u,w)*rdm2(p,r,u,w) + end do end do + + do t=1,N + do u=1,N + + tmp(p,q,r,s) = tmp(p,q,r,s) - ( & + ERI_MO(s,t,p,u)*rdm2(r,t,q,u) + ERI_MO(t,s,p,u)*rdm2(t,r,q,u) & + + ERI_MO(q,u,r,t)*rdm2(p,u,s,t) + ERI_MO(q,u,t,r)*rdm2(p,u,t,s) ) + + end do + end do + + do t=1,N + do u=1,N + do w=1,N + + tmp(p,q,r,s) = tmp(p,q,r,s) + 0.5d0*( & + Kronecker_delta(q,r)*(ERI_MO(u,w,p,t)*rdm2(u,w,s,t) + ERI_MO(s,t,u,w)*rdm2(p,t,u,w)) & + + Kronecker_delta(p,s)*(ERI_MO(q,t,u,w)*rdm2(r,t,u,w) + ERI_MO(u,w,r,t)*rdm2(u,w,q,t)) ) + + end do + end do + end do + end do - end do + end do - end do + + ! Flatten Hessian matrix and add permutations + + pq = 0 + do p=1,N + do q=1,N + + pq = pq + 1 + + rs = 0 + do r=1,N + do s=1,N + + rs = rs + 1 + + hess(pq,rs) = tmp(p,r,q,s) - tmp(r,p,q,s) - tmp(p,r,s,q) + tmp(r,p,s,q) +!! hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r) + + end do + end do + + end do + end do + + deallocate(rdm1,rdm2,tmp) + + allocate(hessInv(N**2,N**2)) + + call inverse_matrix(N**2,hess,hessInv) + + deallocate(hess) + + allocate(Kap(N,N)) + + Kap(:,:) = 0d0 + + pq = 0 + do p=1,N + do q=1,N + + pq = pq + 1 + + rs = 0 + do r=1,N + do s=1,N + + rs = rs + 1 + + Kap(p,q) = Kap(p,q) - hessInv(pq,rs)*grad(rs) + + end do + end do + + end do + end do + + deallocate(hessInv,grad) + + write(*,*) 'kappa' + call matout(N,N,Kap) + write(*,*) + + allocate(ExpKap(N,N)) + call matrix_exponential(N,Kap,ExpKap) + deallocate(Kap) + + write(*,*) 'e^kappa' + call matout(N,N,ExpKap) + write(*,*) + + write(*,*) 'Old orbitals' + call matout(N,N,c) + write(*,*) + + c = matmul(c,ExpKap) + deallocate(ExpKap) + + write(*,*) 'Rotated orbitals' + call matout(N,N,c) + write(*,*) + end do -! Flatten Hessian matrix and add permutations + !-----------------------------------! + ! End Loop for orbital optimization ! + !-----------------------------------! - pq = 0 - do p=1,N - do q=1,N + ! Did it actually converge? - pq = pq + 1 - - rs = 0 - do r=1,N - do s=1,N + if(nItOrb == maxIt) then - rs = rs + 1 + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)'! Convergence failed for orbital optimization !' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' - hess(pq,rs) = tmp(p,r,q,s) - tmp(r,p,q,s) - tmp(p,r,s,q) + tmp(r,p,s,q) -! hess(pq,rs) = tmp(p,q,r,s) - tmp(q,p,r,s) - tmp(p,q,s,r) + tmp(q,p,s,r) + stop - end do - end do - - end do - end do - - call matout(N**2,N**2,hess) - -! allocate(eig(N**2)) - -! call diagonalize_matrix(N**2,hess,eig) - -! call vecout(N**2,eig) - - allocate(hessInv(N**2,N**2)) - - call inverse_matrix(N**2,hess,hessInv) - - allocate(kappa(N,N),ekappa(N,N)) - - kappa(:,:) = 0d0 - - pq = 0 - do p=1,N - do q=1,N - - pq = pq + 1 - - rs = 0 - do r=1,N - do s=1,N - - rs = rs + 1 - - kappa(p,q) = kappa(p,q) + hessInv(pq,rs)*grad(rs) - - end do - end do - - end do - end do - - write(*,*) 'kappa' - call matout(N,N,kappa) - write(*,*) - - call matrix_exponential(N,kappa,ekappa) - - write(*,*) 'e^kappa' - call matout(N,N,ekappa) - write(*,*) - - write(*,*) 'Old orbitals' - call matout(N,N,c) - write(*,*) - - c = matmul(c,ekappa) - - write(*,*) 'Rotated orbitals' - call matout(N,N,c) - write(*,*) + end if ! Testing zone diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index 3897f37..9632c32 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -228,7 +228,7 @@ subroutine RQuAcK(dotest,doRHF,doROHF,dostab,dosearch,doMP2,doMP3,doCCD,dopCCD,d call wall_time(start_CC) call RCC(dotest,doCCD,dopCCD,doDCD,doCCSD,doCCSDT,dodrCCD,dorCCD,docrCCD,dolCCD, & - maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_MO,ENuc,ERHF,eHF,cHF) + maxSCF_CC,thresh_CC,max_diis_CC,nBas,nC,nO,nV,nR,Hc,ERI_AO,ERI_MO,ENuc,ERHF,eHF,cHF) call wall_time(end_CC) t_CC = end_CC - start_CC From eb6e0bcc02f9cc58fe2748f4e098e8c59a256b2a Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Fri, 30 Aug 2024 19:10:54 +0200 Subject: [PATCH 3/3] added non-sym diag with biorthog condition --- src/LR/ppLR.f90 | 138 +++++---- src/make_ninja.py | 4 +- src/utils/non_sym_diag.f90 | 571 +++++++++++++++++++++++++++++++++++++ 3 files changed, 657 insertions(+), 56 deletions(-) create mode 100644 src/utils/non_sym_diag.f90 diff --git a/src/LR/ppLR.f90 b/src/LR/ppLR.f90 index 2e54490..64d533c 100644 --- a/src/LR/ppLR.f90 +++ b/src/LR/ppLR.f90 @@ -1,93 +1,123 @@ -subroutine ppLR(TDA,nOO,nVV,Bpp,Cpp,Dpp,Om1,X1,Y1,Om2,X2,Y2,EcRPA) -! Solve the pp-RPA linear eigenvalue problem +! --- + +subroutine ppLR(TDA, nOO, nVV, Bpp, Cpp, Dpp, Om1, X1, Y1, Om2, X2, Y2, EcRPA) + + ! + ! Solve the pp-RPA linear eigenvalue problem + ! + ! right eigen-problem: H R = R w + ! left eigen-problem: H.T L = L w + ! + ! where L.T R = 1 + ! + ! + ! (+C +B) + ! H = ( ) where C = C.T and D = D.T + ! (-B.T -D) + ! + ! (w1 0) (X1 X2) (+X1 +X2) + ! w = ( ), R = ( ) and L = ( ) + ! (0 w2) (Y1 Y2) (-Y1 -Y2) + ! + ! + ! the normalisation condition reduces to + ! + ! X1.T X2 - Y1.T Y2 = 0 + ! X1.T X1 - Y1.T Y1 = 1 + ! X2.T X2 - Y2.T Y2 = 1 + ! implicit none include 'parameters.h' -! Input variables - - logical,intent(in) :: TDA - integer,intent(in) :: nOO - integer,intent(in) :: nVV - double precision,intent(in) :: Bpp(nVV,nOO) - double precision,intent(in) :: Cpp(nVV,nVV) - double precision,intent(in) :: Dpp(nOO,nOO) + logical, intent(in) :: TDA + integer, intent(in) :: nOO, nVV + double precision, intent(in) :: Bpp(nVV,nOO), Cpp(nVV,nVV), Dpp(nOO,nOO) + double precision, intent(out) :: Om1(nVV), X1(nVV,nVV), Y1(nOO,nVV) + double precision, intent(out) :: Om2(nOO), X2(nVV,nOO), Y2(nOO,nOO) + double precision, intent(out) :: EcRPA -! Local variables + logical :: imp_bio, verbose + integer :: i, j, N + double precision :: EcRPA1, EcRPA2 + double precision :: thr_d, thr_nd, thr_deg + double precision,allocatable :: M(:,:), Z(:,:), Om(:) - double precision :: trace_matrix - double precision :: EcRPA1 - double precision :: EcRPA2 - double precision,allocatable :: M(:,:) - double precision,allocatable :: Z(:,:) - double precision,allocatable :: Om(:) + double precision, external :: trace_matrix -! Output variables - double precision,intent(out) :: Om1(nVV) - double precision,intent(out) :: X1(nVV,nVV) - double precision,intent(out) :: Y1(nOO,nVV) - double precision,intent(out) :: Om2(nOO) - double precision,intent(out) :: X2(nVV,nOO) - double precision,intent(out) :: Y2(nOO,nOO) - double precision,intent(out) :: EcRPA -! Memory allocation + N = nOO + nVV - allocate(M(nOO+nVV,nOO+nVV),Z(nOO+nVV,nOO+nVV),Om(nOO+nVV)) - -!-------------------------------------------------! -! Solve the p-p eigenproblem ! -!-------------------------------------------------! -! ! -! | C B | | X1 X2 | | w1 0 | | X1 X2 | ! -! | | | | = | | | | ! -! | -Bt -D | | Y1 Y2 | | 0 w2 | | Y1 Y2 | ! -! ! -!-------------------------------------------------! + allocate(M(N,N), Z(N,N), Om(N)) if(TDA) then X1(:,:) = +Cpp(:,:) Y1(:,:) = 0d0 - if(nVV > 0) call diagonalize_matrix(nVV,X1,Om1) + if(nVV > 0) call diagonalize_matrix(nVV, X1, Om1) X2(:,:) = 0d0 Y2(:,:) = -Dpp(:,:) - if(nOO > 0) call diagonalize_matrix(nOO,Y2,Om2) + if(nOO > 0) call diagonalize_matrix(nOO, Y2, Om2) else - ! Diagonal blocks - + ! Diagonal blocks M( 1:nVV , 1:nVV) = + Cpp(1:nVV,1:nVV) M(nVV+1:nVV+nOO,nVV+1:nVV+nOO) = - Dpp(1:nOO,1:nOO) - ! Off-diagonal blocks - + ! Off-diagonal blocks M( 1:nVV ,nVV+1:nOO+nVV) = - Bpp(1:nVV,1:nOO) M(nVV+1:nOO+nVV, 1:nVV) = + transpose(Bpp(1:nVV,1:nOO)) -! call matout(nOO,nOO,Dpp) + !! Diagonalize the p-p matrix + !if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV, M, Om, Z) + !! Split the various quantities in p-p and h-h parts + !call sort_ppRPA(nOO, nVV, Om, Z, Om1, X1, Y1, Om2, X2, Y2) - ! Diagonalize the p-p matrix - if(nOO+nVV > 0) call diagonalize_general_matrix(nOO+nVV,M,Om,Z) + thr_d = 1d-6 ! to determine if diagonal elements of L.T x R are close enouph to 1 + thr_nd = 1d-6 ! to determine if non-diagonal elements of L.T x R are close enouph to 1 + thr_deg = 1d-8 ! to determine if two eigenvectors are degenerate or not + imp_bio = .True. ! impose bi-orthogonality + verbose = .False. + call diagonalize_nonsym_matrix(N, M, Z, Om, thr_d, thr_nd, thr_deg, imp_bio, verbose) - ! Split the various quantities in p-p and h-h parts + do i = 1, nOO + Om2(i) = Om(i) + do j = 1, nVV + X2(j,i) = Z(j,i) + enddo + do j = 1, nOO + Y2(j,i) = Z(nVV+j,i) + enddo + enddo - call sort_ppRPA(nOO,nVV,Om,Z,Om1,X1,Y1,Om2,X2,Y2) + do i = 1, nVV + Om1(i) = Om(nOO+i) + do j = 1, nVV + X1(j,i) = M(j,nOO+i) + enddo + do j = 1, nOO + Y1(j,i) = M(nVV+j,nOO+i) + enddo + enddo end if -! Compute the RPA correlation energy + ! Compute the RPA correlation energy + EcRPA = 0.5d0 * (sum(Om1) - sum(Om2) - trace_matrix(nVV, Cpp) - trace_matrix(nOO, Dpp)) + EcRPA1 = +sum(Om1) - trace_matrix(nVV, Cpp) + EcRPA2 = -sum(Om2) - trace_matrix(nOO, Dpp) - EcRPA = 0.5d0*( sum(Om1) - sum(Om2) - trace_matrix(nVV,Cpp) - trace_matrix(nOO,Dpp) ) - EcRPA1 = +sum(Om1) - trace_matrix(nVV,Cpp) - EcRPA2 = -sum(Om2) - trace_matrix(nOO,Dpp) - - if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) & + if(abs(EcRPA - EcRPA1) > 1d-6 .or. abs(EcRPA - EcRPA2) > 1d-6) then print*,'!!! Issue in pp-RPA linear reponse calculation RPA1 != RPA2 !!!' + endif + + deallocate(M, Z, Om) end subroutine + + diff --git a/src/make_ninja.py b/src/make_ninja.py index 1e72639..a5e6e84 100755 --- a/src/make_ninja.py +++ b/src/make_ninja.py @@ -89,8 +89,8 @@ FIX_ORDER_OF_LIBS=-Wl,--start-group if sys.platform in ["linux", "linux2"]: # compiler = compile_gfortran_linux -# compiler = compile_ifort_linux - compiler = compile_olympe + compiler = compile_ifort_linux +# compiler = compile_olympe elif sys.platform == "darwin": compiler = compile_gfortran_mac else: diff --git a/src/utils/non_sym_diag.f90 b/src/utils/non_sym_diag.f90 new file mode 100644 index 0000000..04279ac --- /dev/null +++ b/src/utils/non_sym_diag.f90 @@ -0,0 +1,571 @@ + +! --- + +subroutine diagonalize_nonsym_matrix(N, A, L, e_re, thr_d, thr_nd, thr_deg, imp_bio, verbose) + + ! Diagonalize a non-symmetric matrix + ! + ! Output + ! right-eigenvectors are saved in A + ! left-eigenvectors are saved in L + ! eigenvalues are saved in e = e_re + i e_im + + implicit none + + integer, intent(in) :: N + logical, intent(in) :: imp_bio, verbose + double precision, intent(in) :: thr_d, thr_nd, thr_deg + double precision, intent(inout) :: A(N,N) + double precision, intent(out) :: e_re(N), L(N,N) + + integer :: i, j, ii + integer :: lwork, info + double precision :: accu_d, accu_nd + integer, allocatable :: iorder(:), deg_num(:) + double precision, allocatable :: Atmp(:,:), Ltmp(:,:), work(:), e_im(:) + double precision, allocatable :: S(:,:) + + if(verbose) then + print*, ' Starting a non-Hermitian diagonalization ...' + print*, ' Good Luck ;)' + print*, ' imp_bio = ', imp_bio + endif + + ! --- + ! diagonalize + + allocate(Atmp(N,N), e_im(N)) + Atmp(1:N,1:N) = A(1:N,1:N) + + allocate(work(1)) + lwork = -1 + call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info) + if(info .gt. 0) then + print*,'dgeev failed !!', info + stop + endif + + lwork = max(int(work(1)), 1) + deallocate(work) + allocate(work(lwork)) + + call dgeev('V', 'V', N, Atmp, N, e_re, e_im, L, N, A, N, work, lwork, info) + if(info .ne. 0) then + print*,'dgeev failed !!', info + stop + endif + + deallocate(Atmp, WORK) + + + ! --- + ! check if eigenvalues are real + + i = 1 + ii = 0 + do while(i .le. N) + if(dabs(e_im(i)) .gt. 1.d-12) then + ii = ii + 1 + if(verbose) then + print*, ' Warning: complex eigenvalue !' + print*, i, e_re(i), e_im(i) + if(dabs(e_im(i)/e_re(i)) .lt. 1.d-6) then + print*, ' small enouph to be igored' + else + print*, ' IMAGINARY PART IS SIGNIFANT !!!' + endif + endif + endif + i = i + 1 + enddo + + if(verbose) then + if(ii .eq. 0) print*, ' congratulations :) eigenvalues are real-valued !!' + endif + + + ! --- + ! track & sort the real eigenvalues + + allocate(Atmp(N,N), Ltmp(N,N), iorder(N)) + + do i = 1, N + iorder(i) = i + enddo + call quick_sort(e_re, iorder, N) + + Atmp(:,:) = A(:,:) + Ltmp(:,:) = L(:,:) + do i = 1, N + do j = 1, N + A(j,i) = Atmp(j,iorder(i)) + L(j,i) = Ltmp(j,iorder(i)) + enddo + enddo + + deallocate(Atmp, Ltmp, iorder) + + + + + ! --- + ! check bi-orthog + + allocate(S(N,N)) + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose) + + if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then + + if(verbose) then + print *, ' lapack vectors are normalized and bi-orthogonalized' + endif + + elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then + + if(verbose) then + print *, ' lapack vectors are not normalized but bi-orthogonalized' + endif + + call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.) + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose) + + else + + if(verbose) then + print *, ' lapack vectors are not normalized neither bi-orthogonalized' + endif + + allocate(deg_num(N)) + call reorder_degen_eigvec(N, thr_deg, deg_num, e_re, L, A) + call impose_biorthog_degen_eigvec(N, deg_num, e_re, L, A) + deallocate(deg_num) + + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .false., verbose) + if((accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(N))/dble(N) .lt. thr_d)) then + if(verbose) then + print *, ' lapack vectors are now normalized and bi-orthogonalized' + endif + elseif((accu_nd .lt. thr_nd) .and. (dabs(accu_d - dble(N)) .gt. thr_d)) then + if(verbose) then + print *, ' lapack vectors are now not normalized but bi-orthogonalized' + endif + call check_biorthog_binormalize(N, N, L, A, thr_d, thr_nd, .true.) + call check_biorthog(N, N, L, A, accu_d, accu_nd, S, thr_d, thr_nd, .true., verbose) + else + if(verbose) then + print*, ' bi-orthogonalization failed !' + endif + if(imp_bio) then + print*, ' bi-orthogonalization failed !' + deallocate(S) + stop + endif + endif + + endif + + deallocate(S) + return + +end + +! --- + +subroutine check_biorthog(n, m, Vl, Vr, accu_d, accu_nd, S, thr_d, thr_nd, stop_ifnot, verbose) + + implicit none + + integer, intent(in) :: n, m + logical, intent(in) :: stop_ifnot, verbose + double precision, intent(in) :: Vl(n,m), Vr(n,m) + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(out) :: accu_d, accu_nd, S(m,m) + + integer :: i, j + double precision, allocatable :: SS(:,:) + + if(verbose) then + print *, ' check bi-orthogonality' + endif + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + dabs(S(i,i)) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + + if(verbose) then + if((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d) then + print *, ' non bi-orthogonal vectors !' + print *, ' accu_nd = ', accu_nd + print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + else + print *, ' vectors are bi-orthogonals' + endif + endif + + ! --- + + if(stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) then + print *, ' non bi-orthogonal vectors !' + print *, ' accu_nd = ', accu_nd + print *, ' accu_d = ', dabs(accu_d-dble(m))/dble(m) + stop + endif + +end + +! --- + +subroutine check_biorthog_binormalize(n, m, Vl, Vr, thr_d, thr_nd, stop_ifnot) + + implicit none + + integer, intent(in) :: n, m + logical, intent(in) :: stop_ifnot + double precision, intent(in) :: thr_d, thr_nd + double precision, intent(inout) :: Vl(n,m), Vr(n,m) + + integer :: i, j + double precision :: accu_d, accu_nd, s_tmp + double precision, allocatable :: S(:,:) + + ! --- + + allocate(S(m,m)) + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + do i = 1, m + if(S(i,i) .lt. 0.d0) then + do j = 1, n + Vl(j,i) = -1.d0 * Vl(j,i) + enddo + S(i,i) = -S(i,i) + endif + enddo + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + + ! --- + + if( (accu_nd .lt. thr_nd) .and. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d) ) then + + do i = 1, m + if(S(i,i) <= 0.d0) then + print *, ' overap negative' + print *, i, S(i,i) + exit + endif + if(dabs(S(i,i) - 1.d0) .gt. thr_d) then + s_tmp = 1.d0 / dsqrt(S(i,i)) + do j = 1, n + Vl(j,i) = Vl(j,i) * s_tmp + Vr(j,i) = Vr(j,i) * s_tmp + enddo + endif + + enddo + + endif + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , Vl, size(Vl, 1), Vr, size(Vr, 1) & + , 0.d0, S, size(S, 1) ) + + accu_d = 0.d0 + accu_nd = 0.d0 + do i = 1, m + do j = 1, m + if(i==j) then + accu_d = accu_d + S(i,i) + else + accu_nd = accu_nd + S(j,i) * S(j,i) + endif + enddo + enddo + accu_nd = dsqrt(accu_nd) / dble(m) + + deallocate(S) + + ! --- + + if( stop_ifnot .and. ((accu_nd .gt. thr_nd) .or. (dabs(accu_d-dble(m))/dble(m) .gt. thr_d)) ) then + print *, accu_nd, thr_nd + print *, dabs(accu_d-dble(m))/dble(m), thr_d + print *, ' biorthog_binormalize failed !' + stop + endif + +end + +! --- + +subroutine reorder_degen_eigvec(n, thr_deg, deg_num, e0, L0, R0) + + implicit none + + integer, intent(in) :: n + double precision, intent(in) :: thr_deg + double precision, intent(inout) :: e0(n), L0(n,n), R0(n,n) + integer, intent(out) :: deg_num(n) + + logical :: complex_root + integer :: i, j, k, m, ii, j_tmp + double precision :: ei, ej, de + double precision :: accu_d, accu_nd + double precision :: e0_tmp, L0_tmp(n), R0_tmp(n) + double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:) + + do i = 1, n + deg_num(i) = 1 + enddo + + do i = 1, n-1 + ei = e0(i) + + ! already considered in degen vectors + if(deg_num(i) .eq. 0) cycle + + ii = 0 + do j = i+1, n + ej = e0(j) + de = dabs(ei - ej) + + if(de .lt. thr_deg) then + ii = ii + 1 + + j_tmp = i + ii + + deg_num(j_tmp) = 0 + + e0_tmp = e0(j_tmp) + e0(j_tmp) = e0(j) + e0(j) = e0_tmp + + L0_tmp(1:n) = L0(1:n,j_tmp) + L0(1:n,j_tmp) = L0(1:n,j) + L0(1:n,j) = L0_tmp(1:n) + + R0_tmp(1:n) = R0(1:n,j_tmp) + R0(1:n,j_tmp) = R0(1:n,j) + R0(1:n,j) = R0_tmp(1:n) + endif + enddo + + deg_num(i) = ii + 1 + enddo + + ii = 0 + do i = 1, n + if(deg_num(i) .gt. 1) then + ii = ii + 1 + endif + enddo + + if(ii .eq. 0) then + print*, ' WARNING: bi-orthogonality is lost but there is no degeneracies' + print*, ' rotations may change energy' + stop + endif + +end + +! --- + +subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0) + + implicit none + + integer, intent(in) :: n, deg_num(n) + double precision, intent(in) :: e0(n) + double precision, intent(inout) :: L0(n,n), R0(n,n) + + logical :: complex_root + integer :: i, j, k, m + double precision :: ei, ej, de + double precision :: accu_d, accu_nd + double precision, allocatable :: L(:,:), R(:,:), S(:,:), S_inv_half(:,:) + + !do i = 1, n + ! if(deg_num(i) .gt. 1) then + ! print *, ' degen on', i, deg_num(i), e0(i) + ! endif + !enddo + + ! --- + + do i = 1, n + m = deg_num(i) + + if(m .gt. 1) then + + allocate(L(n,m), R(n,m), S(m,m)) + + do j = 1, m + L(1:n,j) = L0(1:n,i+j-1) + R(1:n,j) = R0(1:n,i+j-1) + enddo + + ! --- + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + accu_nd = 0.d0 + do j = 1, m + do k = 1, m + if(j==k) cycle + accu_nd = accu_nd + dabs(S(j,k)) + enddo + enddo + + if(accu_nd .lt. 1d-12) then + deallocate(S, L, R) + cycle + endif + + call impose_biorthog_svd(n, m, L, R) + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + accu_nd = 0.d0 + do j = 1, m + do k = 1, m + if(j==k) cycle + accu_nd = accu_nd + dabs(S(j,k)) + enddo + enddo + if(accu_nd .gt. 1d-12) then + print*, ' accu_nd =', accu_nd + print*, ' your strategy for degenerates orbitals failed !' + print*, m, 'deg on', i + stop + endif + + deallocate(S) + + ! --- + + do j = 1, m + L0(1:n,i+j-1) = L(1:n,j) + R0(1:n,i+j-1) = R(1:n,j) + enddo + + deallocate(L, R) + + endif + enddo + +end + +! --- + +subroutine impose_biorthog_svd(n, m, L, R) + + implicit none + + integer, intent(in) :: n, m + double precision, intent(inout) :: L(n,m), R(n,m) + + integer :: i, j, num_linear_dependencies + double precision :: threshold + double precision, allocatable :: S(:,:), tmp(:,:) + double precision, allocatable :: U(:,:), V(:,:), Vt(:,:), D(:) + + allocate(S(m,m)) + + call dgemm( 'T', 'N', m, m, n, 1.d0 & + , L, size(L, 1), R, size(R, 1) & + , 0.d0, S, size(S, 1) ) + + ! --- + + allocate(U(m,m), Vt(m,m), D(m)) + + call svd(S, m, U, m, D, Vt, m, m, m) + + deallocate(S) + + threshold = 1.d-6 + num_linear_dependencies = 0 + do i = 1, m + if(abs(D(i)) <= threshold) then + D(i) = 0.d0 + num_linear_dependencies = num_linear_dependencies + 1 + else + D(i) = 1.d0 / dsqrt(D(i)) + endif + enddo + if(num_linear_dependencies > 0) then + write(*,*) ' linear dependencies = ', num_linear_dependencies + write(*,*) ' m = ', m + stop + endif + + allocate(V(m,m)) + do i = 1, m + do j = 1, m + V(j,i) = Vt(i,j) + enddo + enddo + deallocate(Vt) + + ! --- + + ! R <-- R x V x D^{-0.5} + ! L <-- L x U x D^{-0.5} + + do i = 1, m + do j = 1, m + V(j,i) = V(j,i) * D(i) + U(j,i) = U(j,i) * D(i) + enddo + enddo + + allocate(tmp(n,m)) + tmp(:,:) = R(:,:) + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , tmp, size(tmp, 1), V, size(V, 1) & + , 0.d0, R, size(R, 1)) + + tmp(:,:) = L(:,:) + call dgemm( 'N', 'N', n, m, m, 1.d0 & + , tmp, size(tmp, 1), U, size(U, 1) & + , 0.d0, L, size(L, 1)) + + deallocate(tmp, U, V, D) + +end + +! --- +