diff --git a/input/methods b/input/methods index f5a20e0..c1d2dda 100644 --- a/input/methods +++ b/input/methods @@ -13,6 +13,6 @@ # G0W0 evGW qsGW F F F # G0T0 evGT qsGT - T T F + T F F # MCMP2 F diff --git a/input/options b/input/options index b9af96e..8efc658 100644 --- a/input/options +++ b/input/options @@ -11,6 +11,6 @@ # GW: maxSCF thresh DIIS n_diis COHSEX SOSEX BSE TDA G0W GW0 lin eta 256 0.00001 T 5 F F F F F F T 0.000 # ACFDT: AC Kx XBS - T T T + T F T # MCMP2: nMC nEq nWalk dt nPrint iSeed doDrift 1000000 100000 10 0.3 10000 1234 T diff --git a/src/QuAcK/G0T0.f90 b/src/QuAcK/G0T0.f90 index d6648d9..2b54b18 100644 --- a/src/QuAcK/G0T0.f90 +++ b/src/QuAcK/G0T0.f90 @@ -184,16 +184,16 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA,singlet_manifold,triplet_m ! Compute the ppRPA correlation energy - ispin = 1 - iblock = 3 - call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0(:),ERI(:,:,:,:), & - Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) - ispin = 2 - iblock = 4 - call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0(:),ERI(:,:,:,:), & - Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) - EcRPA(1) = EcRPA(1) - EcRPA(2) - EcRPA(2) = 3d0*EcRPA(2) +! ispin = 1 +! iblock = 3 +! call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOs,nVVs,eG0T0(:),ERI(:,:,:,:), & +! Omega1s(:),X1s(:,:),Y1s(:,:),Omega2s(:),X2s(:,:),Y2s(:,:),EcRPA(ispin)) +! ispin = 2 +! iblock = 4 +! call linear_response_pp(iblock,.false.,.false.,nBas,nC,nO,nV,nR,nOOt,nVVt,eG0T0(:),ERI(:,:,:,:), & +! Omega1t(:),X1t(:,:),Y1t(:,:),Omega2t(:),X2t(:,:),Y2t(:,:),EcRPA(ispin)) +! EcRPA(1) = EcRPA(1) - EcRPA(2) +! EcRPA(2) = 3d0*EcRPA(2) write(*,*) write(*,*)'-------------------------------------------------------------------------------' diff --git a/src/QuAcK/sort_ppRPA.f90 b/src/QuAcK/sort_ppRPA.f90 index 258dd9c..6db5120 100644 --- a/src/QuAcK/sort_ppRPA.f90 +++ b/src/QuAcK/sort_ppRPA.f90 @@ -129,73 +129,105 @@ subroutine sort_ppRPA(ortho_eigvec,nOO,nVV,Omega,Z,Omega1,X1,Y1,Omega2,X2,Y2) ! ! Find degenerate eigenvalues -! deg1 = 1 -! ab_start = 1 + deg1 = 1 + ab_start = 1 -! do ab=1,nVV + do ab=2,nVV -! if(ab < nVV .and. abs(Omega1(ab) - Omega1(ab+1)) < 1d-6) then + if(abs(Omega1(ab-1) - Omega1(ab)) < 1d-6) then -! if(deg1 == 1) ab_start = ab -! deg1 = deg1 + 1 + deg1 = deg1 + 1 -! else + if(ab == nVV) then + + 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) + + end if + + else + + ab_end = ab - 1 -! ab_end = ab ! print*,'deg = ',deg1,ab_start,ab_end -! allocate(S1(deg1,deg1),O1(deg1,deg1)) + 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) + 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) + deallocate(S1,O1) -! deg1 = 1 -! ab_start = ab + 1 + deg1 = 1 + ab_start = ab -! end if -! end do + end if + end do -! deg2 = 1 -! ij_start = 1 + deg2 = 1 + ij_start = 1 -! do ij=1,nOO + do ij=2,nOO -! if(ij < nOO .and. abs(Omega2(ij) - Omega2(ij+1)) < 1d-6) then + if(abs(Omega2(ij-1) - Omega2(ij)) < 1d-6) then -! if(deg2 == 1) ij_start = ij -! deg2 = deg2 + 1 + deg2 = deg2 + 1 -! else + if(ij == nOO) then + + 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) + + end if + + else + + ij_end = ij - 1 -! ij_end = ij ! print*,'deg = ',deg2,ij_start,ij_end -! allocate(S2(deg2,deg2),O2(deg2,deg2)) + 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) + 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) + deallocate(S2,O2) -! deg2 = 1 -! ij_start = ij + 1 + deg2 = 1 + ij_start = ij -! end if -! end do + 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)) +! 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) +! 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) +! Z1 = matmul(Z1,O1) +! Z2 = matmul(Z2,O2) end if