diff --git a/src/LR/phGLR_A.f90 b/src/LR/phGLR_A.f90 index 12d33d5..3b6f6c4 100644 --- a/src/LR/phGLR_A.f90 +++ b/src/LR/phGLR_A.f90 @@ -24,6 +24,9 @@ subroutine phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) double precision,external :: Kronecker_delta integer :: i,j,a,b,ia,jb + integer :: nn,jb0 + logical :: i_eq_j + double precision :: ct1,ct2 ! Output variables @@ -35,7 +38,36 @@ subroutine phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) if(dRPA) delta_dRPA = 1d0 ! Build A matrix for spin orbitals + ! nn = nOrb - nR - nO + ! ct1 = lambda + ! ct2 = - (1d0 - delta_dRPA) * lambda + + ! !$OMP PARALLEL DEFAULT(NONE) & + ! !$OMP PRIVATE (i, a, j, b, i_eq_j, ia, jb0, jb) & + ! !$OMP SHARED (nC, nO, nR, nOrb, nn, ct1, ct2, e, ERI, Aph) + ! !$OMP DO COLLAPSE(2) + ! do i = nC+1, nO + ! do a = nO+1, nOrb-nR + ! ia = a - nO + (i - nC - 1) * nn + ! do j = nC+1, nO + ! i_eq_j = i == j + ! jb0 = (j - nC - 1) * nn - nO + ! do b = nO+1, nOrb-nR + ! jb = b + jb0 + + ! Aph(ia,jb) = ct1 * ERI(b,i,j,a) + ct2 * ERI(b,j,a,i) + ! if(i_eq_j) then + ! if(a == b) Aph(ia,jb) = Aph(ia,jb) + e(a) - e(i) + ! endif + + ! enddo + ! enddo + + ! enddo + ! enddo + ! !$OMP END DO + ! !$OMP END PARALLEL ia = 0 do i=nC+1,nO do a=nO+1,nOrb-nR @@ -53,4 +85,4 @@ subroutine phGLR_A(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,e,ERI,Aph) end do end do -end subroutine +end subroutine diff --git a/src/LR/phGLR_B.f90 b/src/LR/phGLR_B.f90 index acf94ec..bfc1ce5 100644 --- a/src/LR/phGLR_B.f90 +++ b/src/LR/phGLR_B.f90 @@ -22,7 +22,9 @@ subroutine phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) double precision :: delta_dRPA integer :: i,j,a,b,ia,jb - + integer :: nn,jb0 + double precision :: ct1,ct2 + ! Output variables double precision,intent(out) :: Bph(nS,nS) @@ -33,7 +35,32 @@ subroutine phGLR_B(dRPA,nOrb,nC,nO,nV,nR,nS,lambda,ERI,Bph) if(dRPA) delta_dRPA = 1d0 ! Build B matrix for spin orbitals + ! nn = nOrb - nR - nO + ! ct1 = lambda + ! ct2 = - (1d0 - delta_dRPA) * lambda + + ! !$OMP PARALLEL DEFAULT(NONE) & + ! !$OMP PRIVATE (i, a, j, b, ia, jb0, jb) & + ! !$OMP SHARED (nC, nO, nR, nOrb, nn, ct1, ct2, ERI, Bph) + ! !$OMP DO COLLAPSE(2) + ! do i = nC+1, nO + ! do a = nO+1, nOrb-nR + ! ia = a - nO + (i - nC - 1) * nn + ! do j = nC+1, nO + ! jb0 = (j - nC - 1) * nn - nO + ! do b = nO+1, nOrb-nR + ! jb = b + jb0 + + ! Bph(ia,jb) = ct1 * ERI(i,j,a,b) + ct2 * ERI(i,j,b,a) + + ! enddo + ! enddo + + ! enddo + ! enddo + ! !$OMP END DO + ! !$OMP END PARALLEL ia = 0 do i=nC+1,nO do a=nO+1,nOrb-nR diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index b47a38a..3dc247a 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -7,8 +7,8 @@ subroutine GParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c ! Hard-coded parameters - logical :: print_phLR = .false. - logical :: print_ppLR = .false. + logical :: print_phLR = .true. + logical :: print_ppLR = .true. ! Input variables diff --git a/src/Parquet/G_irred_Parquet_self_energy.f90 b/src/Parquet/G_irred_Parquet_self_energy.f90 index ad78e8c..e091560 100644 --- a/src/Parquet/G_irred_Parquet_self_energy.f90 +++ b/src/Parquet/G_irred_Parquet_self_energy.f90 @@ -24,6 +24,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& integer :: p,n double precision :: eps,dem1,dem2,reg,reg1,reg2 double precision :: num + double precision :: start_t,end_t,t ! Output variables double precision,intent(out) :: SigC(nOrb) @@ -39,6 +40,7 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& !-----------------------------! ! GF2 part of the self-energy ! !-----------------------------! + call wall_time(start_t) do p=nC+1,nOrb-nR ! 2h1p sum do i=nC+1,nO @@ -71,11 +73,15 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& end do end do end do + call wall_time(end_t) + t = end_t - start_t + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building GF(2) self-energy =',t,' seconds' + write(*,*) !-----------------------------! ! eh part of the self-energy ! !-----------------------------! - + call wall_time(start_t) do p=nC+1,nOrb-nR do i=nC+1,nO @@ -182,11 +188,15 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& end do ! i end do ! p - + call wall_time(end_t) + t = end_t - start_t + + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building eh self-energy =',t,' seconds' + write(*,*) !-----------------------------! ! pp part of the self-energy ! !-----------------------------! - + call wall_time(start_t) do p=nC+1,nOrb-nR do i=nC+1,nO @@ -303,7 +313,11 @@ subroutine G_Parquet_self_energy(eta,nOrb,nC,nO,nV,nR,nS,nOO,nVV,eQP,ERI,& end do ! a end do ! p - + call wall_time(end_t) + t = end_t - start_t + + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for building pp self-energy =',t,' seconds' + write(*,*) !-----------------------------! ! Renormalization factor ! !-----------------------------! diff --git a/src/Parquet/G_screened_integrals.f90 b/src/Parquet/G_screened_integrals.f90 index aa66cb2..31cca2c 100644 --- a/src/Parquet/G_screened_integrals.f90 +++ b/src/Parquet/G_screened_integrals.f90 @@ -38,9 +38,9 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho rho(p,q,ia) = rho(p,q,ia) & - + (1d0*ERI(q,j,p,b) - 1d0*ERI(q,j,b,p) & + + (ERI(q,j,p,b) - ERI(q,j,b,p) & - 0d0*eh_Phi(q,j,b,p) + 0d0*pp_Phi(q,j,p,b)) * X & - + (1d0*ERI(q,b,p,j) - 1d0*ERI(q,b,j,p) & + + (ERI(q,b,p,j) - ERI(q,b,j,p) & - 0d0*eh_Phi(q,b,j,p) + 0d0*pp_Phi(q,b,p,j)) * Y @@ -110,7 +110,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do d=c+1,nOrb-nR cd = cd + 1 - rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) & + rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & + 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X1(cd,ab) end do @@ -121,7 +121,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do l=k+1,nO kl = kl + 1 - rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) & + rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & + 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab) end do @@ -134,12 +134,13 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do i=nC+1,nO do j=i+1,nO ij = ij + 1 + cd = 0 do c=nO+1,nOrb-nR do d=c+1,nOrb-nR cd = cd + 1 - rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) & + rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & + 0d0*eh_Phi(p,q,c,d) - 0d0*eh_Phi(p,q,d,c) ) * X2(cd,ij) end do @@ -150,7 +151,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do l=k+1,nO kl = kl + 1 - rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) & + rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & + 0d0*eh_Phi(p,q,k,l) - 0d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij) end do diff --git a/src/Parquet/RParquet.f90 b/src/Parquet/RParquet.f90 index a8634b6..f519c43 100644 --- a/src/Parquet/RParquet.f90 +++ b/src/Parquet/RParquet.f90 @@ -7,8 +7,8 @@ subroutine RParquet(TDAeh,TDApp,linearize,eta,ENuc,max_it_1b,conv_1b,max_it_2b,c ! Hard-coded parameters - logical :: print_phLR = .false. - logical :: print_ppLR = .false. + logical :: print_phLR = .true. + logical :: print_ppLR = .true. ! Input variables diff --git a/src/QuAcK/RQuAcK.f90 b/src/QuAcK/RQuAcK.f90 index ff64438..20c5642 100644 --- a/src/QuAcK/RQuAcK.f90 +++ b/src/QuAcK/RQuAcK.f90 @@ -380,7 +380,7 @@ subroutine RQuAcK(working_dir,use_gpu,dotest,doRHF,doROHF,dostab,dosearch,doMP2, if(doParquet) then call wall_time(start_Parquet) - call RParquet(TDAeh,TDApp,lin_parquet,reg_parquet,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,ERHF,eGW,ERI_MO) + call RParquet(TDAeh,TDApp,lin_parquet,reg_parquet,ENuc,max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,ERHF,eHF,ERI_MO) call wall_time(end_Parquet) t_Parquet = end_Parquet - start_Parquet