diff --git a/input/basis b/input/basis index 4c61da7..27724b4 100644 --- a/input/basis +++ b/input/basis @@ -1,16 +1,37 @@ -1 6 -S 4 1.00 - 234.0000000 0.0025870 - 35.1600000 0.0195330 - 7.9890000 0.0909980 - 2.2120000 0.2720500 +1 10 +S 8 1.00 + 24350.0000000 0.0005020 + 3650.0000000 0.0038810 + 829.6000000 0.0199970 + 234.0000000 0.0784180 + 75.6100000 0.2296760 + 26.7300000 0.4327220 + 9.9270000 0.3506420 + 1.1020000 -0.0076450 +S 8 1.00 + 24350.0000000 -0.0001180 + 3650.0000000 -0.0009150 + 829.6000000 -0.0047370 + 234.0000000 -0.0192330 + 75.6100000 -0.0603690 + 26.7300000 -0.1425080 + 9.9270000 -0.1777100 + 1.1020000 0.6058360 S 1 1.00 - 0.6669000 1.0000000 + 2.8360000 1.0000000 S 1 1.00 - 0.2089000 1.0000000 + 0.3782000 1.0000000 +P 3 1.00 + 54.7000000 0.0171510 + 12.4300000 0.1076560 + 3.6790000 0.3216810 P 1 1.00 - 3.0440000 1.0000000 + 1.1430000 1.0000000 P 1 1.00 - 0.7580000 1.0000000 + 0.3300000 1.0000000 D 1 1.00 - 1.9650000 1.0000000 + 4.0140000 1.0000000 +D 1 1.00 + 1.0960000 1.0000000 +F 1 1.00 + 2.5440000 1.0000000 diff --git a/input/molecule b/input/molecule index c78e87e..edeba31 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,4 @@ # nAt nEla nElb nCore nRyd - 1 1 1 0 0 + 1 5 5 0 0 # Znuc x y z - He 0.0 0.0 0.0 + Ne 0.0 0.0 0.0 diff --git a/input/options b/input/options index 721bbaa..bc5c578 100644 --- a/input/options +++ b/input/options @@ -3,7 +3,7 @@ # MP: # CC: maxSCF thresh DIIS n_diis - 64 0.00001 F 1 + 64 0.0000001 F 1 # CIS/TDHF/BSE: singlet triplet T T # GF: maxSCF thresh DIIS n_diis renormalization diff --git a/input/weight b/input/weight index 4c61da7..27724b4 100644 --- a/input/weight +++ b/input/weight @@ -1,16 +1,37 @@ -1 6 -S 4 1.00 - 234.0000000 0.0025870 - 35.1600000 0.0195330 - 7.9890000 0.0909980 - 2.2120000 0.2720500 +1 10 +S 8 1.00 + 24350.0000000 0.0005020 + 3650.0000000 0.0038810 + 829.6000000 0.0199970 + 234.0000000 0.0784180 + 75.6100000 0.2296760 + 26.7300000 0.4327220 + 9.9270000 0.3506420 + 1.1020000 -0.0076450 +S 8 1.00 + 24350.0000000 -0.0001180 + 3650.0000000 -0.0009150 + 829.6000000 -0.0047370 + 234.0000000 -0.0192330 + 75.6100000 -0.0603690 + 26.7300000 -0.1425080 + 9.9270000 -0.1777100 + 1.1020000 0.6058360 S 1 1.00 - 0.6669000 1.0000000 + 2.8360000 1.0000000 S 1 1.00 - 0.2089000 1.0000000 + 0.3782000 1.0000000 +P 3 1.00 + 54.7000000 0.0171510 + 12.4300000 0.1076560 + 3.6790000 0.3216810 P 1 1.00 - 3.0440000 1.0000000 + 1.1430000 1.0000000 P 1 1.00 - 0.7580000 1.0000000 + 0.3300000 1.0000000 D 1 1.00 - 1.9650000 1.0000000 + 4.0140000 1.0000000 +D 1 1.00 + 1.0960000 1.0000000 +F 1 1.00 + 2.5440000 1.0000000 diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 22c8012..33cce74 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -345,8 +345,8 @@ program QuAcK call cpu_time(start_CCD) ! call ring_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) - call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) -! call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) +! call ladder_CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) + call CCD(maxSCF_CC,thresh_CC,n_diis_CC,nBas,nEl,ERI_MO_basis,ENuc,ERHF,eHF) call cpu_time(end_CCD) t_CCD = end_CCD - start_CCD diff --git a/src/QuAcK/form_ladder_r.f90 b/src/QuAcK/form_ladder_r.f90 new file mode 100644 index 0000000..b66e1bd --- /dev/null +++ b/src/QuAcK/form_ladder_r.f90 @@ -0,0 +1,59 @@ +subroutine form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2) + +! Form residuals for ladder CCD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: OOOO(nO,nO,nO,nO) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + double precision,intent(in) :: VVVV(nV,nV,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: r2(nO,nO,nV,nV) + + r2(:,:,:,:) = 0d0 + + do i=1,nO + do j=1,nO + do a=1,nV + do b=1,nV + + do k=1,nO + do l=1,nO + r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*t2(k,l,a,b)*OOOO(i,j,k,l) + end do + end do + + do c=1,nV + do d=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + 0.5d0*VVVV(c,d,a,b)*t2(i,j,c,d) + end do + end do + + do k=1,nO + do l=1,nO + do c=1,nV + do d=1,nV + r2(i,j,a,b) = r2(i,j,a,b) + 0.25d0*t2(i,j,c,d)*OOVV(k,l,c,d)*t2(k,l,a,b) + + end do + end do + end do + end do + + end do + end do + end do + end do + +end subroutine form_ladder_r diff --git a/src/QuAcK/form_ring_r.f90 b/src/QuAcK/form_ring_r.f90 new file mode 100644 index 0000000..099af66 --- /dev/null +++ b/src/QuAcK/form_ring_r.f90 @@ -0,0 +1,55 @@ +subroutine form_ring_r(nO,nV,OVVO,OOVV,t2,r2) + +! Form residuals for ring CCD + + implicit none + +! Input variables + + integer,intent(in) :: nO,nV + double precision,intent(in) :: t2(nO,nO,nV,nV) + double precision,intent(in) :: OVVO(nO,nV,nV,nO) + double precision,intent(in) :: OOVV(nO,nO,nV,nV) + +! Local variables + + integer :: i,j,k,l + integer :: a,b,c,d + +! Output variables + + double precision,intent(out) :: r2(nO,nO,nV,nV) + + r2(:,:,:,:) = 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 + + r2(i,j,a,b) = r2(i,j,a,b) + OVVO(i,c,a,k)*t2(k,j,c,b) + OVVO(j,c,b,k)*t2(i,k,a,c) + + end do + end do + + do k=1,nO + do l=1,nO + do c=1,nV + do d=1,nV + + r2(i,j,a,b) = r2(i,j,a,b) + t2(i,k,a,c)*OOVV(k,l,c,d)*t2(l,j,d,b) + + end do + end do + end do + end do + + end do + end do + end do + end do + +end subroutine form_ring_r diff --git a/src/QuAcK/ladder_CCD.f90 b/src/QuAcK/ladder_CCD.f90 new file mode 100644 index 0000000..fd6da4a --- /dev/null +++ b/src/QuAcK/ladder_CCD.f90 @@ -0,0 +1,193 @@ +subroutine ladder_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) + +! Ladder CCD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas,nEl + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nBas2 + integer :: nO + integer :: nV + integer :: nSCF + double precision :: Conv + double precision :: EcMP2,EcMP3,EcMP4 + double precision :: ECCD,EcCCD + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVOV(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: X1(:,:,:,:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: X3(:,:) + double precision,allocatable :: X4(:,:,:,:) + + double precision,allocatable :: u(:,:,:,:) + double precision,allocatable :: v(:,:,:,:) + + double precision,allocatable :: r2(:,:,:,:) + double precision,allocatable :: t2(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************' + write(*,*)'| ladder-CCD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas2 = 2*nBas + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas2,nBas2,nBas2,nBas2)) + + call antisymmetrize_ERI(2,nBas2,sERI,dbERI) + + deallocate(sERI) + +! Define occupied and virtual spaces + + nO = 2*nEl + nV = nBas2 - nO + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = seHF(1:nO) + eV(:) = seHF(nO+1:nBas2) + + call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV) + + deallocate(seHF) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVOV(nO,nV,nO,nV),VVVV(nV,nV,nV,nV)) + + OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2) + OVOV(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2, 1:nO ,nO+1:nBas2) + VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + EcMP4 = 0d0 + +! Initialization + + allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) + allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| ladder-CCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Compute residual + + call form_ladder_r(nO,nV,OOOO,OOVV,VVVV,t2,r2) + + r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) + +! Check convergence + + Conv = maxval(abs(r2(:,:,:,:))) + +! Update amplitudes + + t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + +! Compute correlation energy + + EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + + if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.)) + +! Dump results + + ECCD = ERHF + EcCCD + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' + + enddo + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Moller-Plesset energies + + write(*,*) + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3 + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4 + write(*,*) + +end subroutine ladder_CCD diff --git a/src/QuAcK/linear_response_pp.f90 b/src/QuAcK/linear_response_pp.f90 index cb9df48..c0b796d 100644 --- a/src/QuAcK/linear_response_pp.f90 +++ b/src/QuAcK/linear_response_pp.f90 @@ -70,23 +70,25 @@ subroutine linear_response_pp(ispin,BSE,nBas,nC,nO,nV,nR,nOO,nVV,e,ERI,Omega1,X1 Z(:,:) = M(:,:) call diagonalize_general_matrix(nOO+nVV,M(:,:),Omega(:),Z(:,:)) +! call diagonalize_matrix(nOO+nVV,Z(:,:),Omega(:)) - write(*,*) 'pp-RPA excitation energies' - call matout(nOO+nVV,1,Omega(:)) - write(*,*) +! write(*,*) 'pp-RPA excitation energies' +! call matout(nOO+nVV,1,Omega(:)) +! write(*,*) + +! call matout(nOO+nVV,nOO+nVV,matmul(transpose(Z(:,:)),matmul(W(:,:),Z(:,:)))) +! write(*,*) ! Split the various quantities in p-p and h-h parts - Omega1(:) = Omega(nOO+1:nOO+nVV) - Omega2(:) = Omega(1:nOO) + call sort_ppRPA(nOO,nVV,Omega(:),Z(:,:),Omega1(:),X1(:,:),Y1(:,:),Omega2(:),X2(:,:),Y2(:,:)) +! Omega1(:) = Omega(nOO+1:nOO+nVV) +! Omega2(:) = Omega(1:nOO) - X1(:,:) = Z(nOO+1:nOO+nVV,nOO+1:nOO+nVV) - Y1(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV) - X2(:,:) = Z(nOO+1:nOO+nVV, 1:nOO ) - Y2(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV) - - if(minval(Omega1(:)) < 0d0) call print_warning('You may have instabilities in pp-RPA!!') - if(maxval(Omega2(:)) > 0d0) call print_warning('You may have instabilities in pp-RPA!!') +! X1(:,:) = Z(nOO+1:nOO+nVV,nOO+1:nOO+nVV) +! Y1(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV) +! X2(:,:) = Z(nOO+1:nOO+nVV, 1:nOO ) +! Y2(:,:) = Z( 1:nOO ,nOO+1:nOO+nVV) ! Compute the RPA correlation energy diff --git a/src/QuAcK/ppRPA.f90 b/src/QuAcK/ppRPA.f90 index 9fec1b2..7fd4536 100644 --- a/src/QuAcK/ppRPA.f90 +++ b/src/QuAcK/ppRPA.f90 @@ -110,9 +110,9 @@ subroutine ppRPA(singlet_manifold,triplet_manifold,nBas,nC,nO,nV,nR,ENuc,ERHF,ER write(*,*) write(*,*)'-------------------------------------------------------------------------------' write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (singlet) =',Ec_ppRPA(1) - write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (triplet) =',Ec_ppRPA(2) - write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy =',Ec_ppRPA(1) + Ec_ppRPA(2) - write(*,'(2X,A40,F15.6)') 'pp-RPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy (triplet) =',3d0*Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'pp-RPA correlation energy =',Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) + write(*,'(2X,A40,F15.6)') 'pp-RPA total energy =',ENuc + ERHF + Ec_ppRPA(1) + 3d0*Ec_ppRPA(2) write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/QuAcK/ring_CCD.f90 b/src/QuAcK/ring_CCD.f90 new file mode 100644 index 0000000..4cef290 --- /dev/null +++ b/src/QuAcK/ring_CCD.f90 @@ -0,0 +1,193 @@ +subroutine ring_CCD(maxSCF,thresh,max_diis,nBas,nEl,ERI,ENuc,ERHF,eHF) + +! Ring CCD module + + implicit none + +! Input variables + + integer,intent(in) :: maxSCF + integer,intent(in) :: max_diis + double precision,intent(in) :: thresh + + integer,intent(in) :: nBas,nEl + double precision,intent(in) :: ENuc,ERHF + double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + +! Local variables + + integer :: nBas2 + integer :: nO + integer :: nV + integer :: nSCF + double precision :: Conv + double precision :: EcMP2,EcMP3,EcMP4 + double precision :: ECCD,EcCCD + double precision,allocatable :: seHF(:) + double precision,allocatable :: sERI(:,:,:,:) + double precision,allocatable :: dbERI(:,:,:,:) + + double precision,allocatable :: eO(:) + double precision,allocatable :: eV(:) + double precision,allocatable :: delta_OOVV(:,:,:,:) + + double precision,allocatable :: OOOO(:,:,:,:) + double precision,allocatable :: OOVV(:,:,:,:) + double precision,allocatable :: OVVO(:,:,:,:) + double precision,allocatable :: VVVV(:,:,:,:) + + double precision,allocatable :: X1(:,:,:,:) + double precision,allocatable :: X2(:,:) + double precision,allocatable :: X3(:,:) + double precision,allocatable :: X4(:,:,:,:) + + double precision,allocatable :: u(:,:,:,:) + double precision,allocatable :: v(:,:,:,:) + + double precision,allocatable :: r2(:,:,:,:) + double precision,allocatable :: t2(:,:,:,:) + +! Hello world + + write(*,*) + write(*,*)'**************************************' + write(*,*)'| ring-CCD calculation |' + write(*,*)'**************************************' + write(*,*) + +! Spatial to spin orbitals + + nBas2 = 2*nBas + + allocate(seHF(nBas2),sERI(nBas2,nBas2,nBas2,nBas2)) + + call spatial_to_spin_MO_energy(nBas,eHF,nBas2,seHF) + call spatial_to_spin_ERI(nBas,ERI,nBas2,sERI) + +! Antysymmetrize ERIs + + allocate(dbERI(nBas2,nBas2,nBas2,nBas2)) + + call antisymmetrize_ERI(2,nBas2,sERI,dbERI) + + deallocate(sERI) + +! Define occupied and virtual spaces + + nO = 2*nEl + nV = nBas2 - nO + +! Form energy denominator + + allocate(eO(nO),eV(nV)) + allocate(delta_OOVV(nO,nO,nV,nV)) + + eO(:) = seHF(1:nO) + eV(:) = seHF(nO+1:nBas2) + + call form_delta_OOVV(nO,nV,eO,eV,delta_OOVV) + + deallocate(seHF) + +! Create integral batches + + allocate(OOOO(nO,nO,nO,nO),OOVV(nO,nO,nV,nV),OVVO(nO,nV,nV,nO),VVVV(nV,nV,nV,nV)) + + OOOO(:,:,:,:) = dbERI( 1:nO , 1:nO , 1:nO , 1:nO ) + OOVV(:,:,:,:) = dbERI( 1:nO , 1:nO ,nO+1:nBas2,nO+1:nBas2) + OVVO(:,:,:,:) = dbERI( 1:nO ,nO+1:nBas2,nO+1:nBas2, 1:nO ) + VVVV(:,:,:,:) = dbERI(nO+1:nBas2,nO+1:nBas2,nO+1:nBas2,nO+1:nBas2) + + deallocate(dbERI) + +! MP2 guess amplitudes + + allocate(t2(nO,nO,nV,nV)) + + t2(:,:,:,:) = -OOVV(:,:,:,:)/delta_OOVV(:,:,:,:) + + EcMP2 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + EcMP4 = 0d0 + +! Initialization + + allocate(r2(nO,nO,nV,nV),u(nO,nO,nV,nV),v(nO,nO,nV,nV)) + allocate(X1(nO,nO,nO,nO),X2(nV,nV),X3(nO,nO),X4(nO,nO,nV,nV)) + + Conv = 1d0 + nSCF = 0 + +!------------------------------------------------------------------------ +! Main SCF loop +!------------------------------------------------------------------------ + write(*,*) + write(*,*)'----------------------------------------------------' + write(*,*)'| ring-CCD calculation |' + write(*,*)'----------------------------------------------------' + write(*,'(1X,A1,1X,A3,1X,A1,1X,A16,1X,A1,1X,A10,1X,A1,1X,A10,1X,A1,1X)') & + '|','#','|','E(CCD)','|','Ec(CCD)','|','Conv','|' + write(*,*)'----------------------------------------------------' + + do while(Conv > thresh .and. nSCF < maxSCF) + +! Increment + + nSCF = nSCF + 1 + +! Compute residual + + call form_ring_r(nO,nV,OVVO,OOVV,t2,r2) + + r2(:,:,:,:) = OOVV(:,:,:,:) + delta_OOVV(:,:,:,:)*t2(:,:,:,:) + r2(:,:,:,:) + +! Check convergence + + Conv = maxval(abs(r2(:,:,:,:))) + +! Update amplitudes + + t2(:,:,:,:) = t2(:,:,:,:) - r2(:,:,:,:)/delta_OOVV(:,:,:,:) + +! Compute correlation energy + + EcCCD = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2,.true.)) + + if(nSCF == 1) EcMP3 = 0.25d0*dot_product(pack(OOVV,.true.),pack(t2 + v/delta_OOVV,.true.)) + +! Dump results + + ECCD = ERHF + EcCCD + + write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F10.6,1X,A1,1X,F10.6,1X,A1,1X)') & + '|',nSCF,'|',ECCD+ENuc,'|',EcCCD,'|',Conv,'|' + + enddo + write(*,*)'----------------------------------------------------' +!------------------------------------------------------------------------ +! End of SCF loop +!------------------------------------------------------------------------ + +! Did it actually converge? + + if(nSCF == maxSCF) then + + write(*,*) + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*)' Convergence failed ' + write(*,*)'!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!' + write(*,*) + + stop + + endif + +! Moller-Plesset energies + + write(*,*) + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP2) = ',EcMP2 + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP3) = ',EcMP3 + write(*,'(1X,A15,1X,F10.6)') 'Ec(MP4-SDQ) = ',EcMP4 + write(*,*) + +end subroutine ring_CCD diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 new file mode 100644 index 0000000..d63ff12 --- /dev/null +++ b/src/QuAcK/sort_ppRPA.f90 @@ -0,0 +1,88 @@ +subroutine sort_ppRPA(nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) + +! Compute the metric matrix for pp-RPA + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nOO + integer,intent(in) :: nVV + double precision,intent(in) :: Omega(nOO+nVV) + double precision,intent(in) :: Z(nOO+nVV,nOO+nVV) + +! Local variables + + integer :: pq,ab,ij + double precision,allocatable :: s(:,:) + +! Output variables + + double precision,intent(out) :: Omega1(nVV) + double precision,intent(out) :: X1(nVV,nVV) + double precision,intent(out) :: Y1(nOO,nVV) + double precision,intent(out) :: Omega2(nOO) + double precision,intent(out) :: X2(nVV,nOO) + double precision,intent(out) :: Y2(nOO,nOO) + +! Memory allocation + + allocate(s(nOO+nVV,nOO+nVV)) + +! Initializatiom + + Omega1(:) = 0d0 + Omega2(:) = 0d0 + + ab = 0 + ij = 0 + do pq=1,nOO+nVV + + if(Omega(pq) > 0d0) then + + ab = ab + 1 + Omega1(ab) = Omega(pq) + X1(1:nVV,ab) = Z( 1: nVV,pq) + Y1(1:nOO,ab) = Z(nVV+1:nOO+nVV,pq) + + else + + ij = ij + 1 + Omega2(ij) = Omega(pq) + X2(1:nVV,ij) = Z( 1: nVV,pq) + Y2(1:nOO,ij) = Z(nVV+1:nOO+nVV,pq) + + end if + + end do + + 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!!') + + write(*,*) 'pp-RPA positive excitation energies' + call matout(nVV,1,Omega1(:)) + write(*,*) + + write(*,*) 'pp-RPA negative excitation energies' + call matout(nOO,1,Omega2(:)) + write(*,*) + +! Check eigenvector signatures + +! s( 1: nVV, 1: nVV) = matmul(transpose(X1),X1) - matmul(transpose(Y1),Y1) +! s(nVV+1:nOO+nVV,nVV+1:nOO+nVV) = matmul(transpose(X2),X2) - matmul(transpose(Y2),Y2) + +! write(*,*) 'Signatures pp' +! do ab=1,nVV +! write(*,'(I6,F10.6)') ab,s(ab,ab) +! end do +! write(*,*) + +! write(*,*) 'Signatures hh' +! do ij=1,nOO +! write(*,'(I6,F10.6)') ij,s(ij,ij) +! end do +! write(*,*) + +end subroutine sort_ppRPA