diff --git a/mol/MgO.xyz b/mol/MgO.xyz index 091c39f..78b1186 100644 --- a/mol/MgO.xyz +++ b/mol/MgO.xyz @@ -1,4 +1,4 @@ 2 Mg 0.0000 0.0000 0.0000 -O 0.0000 0.0000 1.728 +O 0.0000 0.0000 1.749 diff --git a/src/Parquet/GParquet.f90 b/src/Parquet/GParquet.f90 index ddda574..cf5ff1c 100644 --- a/src/Parquet/GParquet.f90 +++ b/src/Parquet/GParquet.f90 @@ -7,8 +7,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ! Hard-coded parameters - logical :: linearize = .true. - logical :: TDA = .true. + logical :: TDA = .true. + logical :: linearize = .true. logical :: print_phLR = .false. logical :: print_ppLR = .false. @@ -25,10 +25,10 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, integer :: n_it_1b,n_it_2b double precision :: err_1b,err_2b double precision :: err_eh, err_pp - double precision :: err_eig_eh, err_eig_hh, err_eig_ee - double precision :: start_t, end_t, t - double precision :: start_1b, end_1b, t_1b - double precision :: start_2b, end_2b, t_2b + double precision :: err_eig_eh,err_eig_pp,err_eig_hh,err_eig_ee + double precision :: start_t,end_t,t + double precision :: start_1b,end_1b,t_1b + double precision :: start_2b,end_2b,t_2b integer :: nOO,nVV @@ -100,11 +100,11 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, write(*,*)'---------------------------------------------------------------' write(*,*)' Parquet parameters for one-body and two-body self-consistency ' write(*,*)'---------------------------------------------------------------' - write(*,'(1X,A50,1X,I5)') 'Maximum number for one-body self-consistency:', max_it_1b - write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for one-body energies:', conv_1b + write(*,'(1X,A50,1X,I5)') 'Maximum number of one-body iteration:',max_it_1b + write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for one-body energies:',conv_1b write(*,*)'---------------------------------------------------------------' - write(*,'(1X,A50,1X,I5)') 'Maximum number for two-body self-consistency:', max_it_2b - write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for two-body energies:', conv_2b + write(*,'(1X,A50,1X,I5)') 'Maximum number of two-body iteration:',max_it_2b + write(*,'(1X,A50,1X,E10.5)') 'Convergence threshold for two-body energies:',conv_2b write(*,*)'---------------------------------------------------------------' write(*,*) @@ -177,14 +177,10 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ! eh channel ! !-----------------! - write(*,*)' ------------------------------' - write(*,*)' | Diagonalizing ehBSE |' - write(*,*)' ------------------------------' - write(*,*) + write(*,*) 'Diagonalizing eh BSE problem...' allocate(Aph(nS,nS),Bph(nS,nS),eh_Om(nS),XpY(nS,nS),XmY(nS,nS),eh_Gam_A(nS,nS),eh_Gam_B(nS,nS)) - Aph(:,:) = 0d0 Bph(:,:) = 0d0 @@ -208,14 +204,12 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, Aph(:,:) = Aph(:,:) + eh_Gam_A(:,:) Bph(:,:) = Bph(:,:) + eh_Gam_B(:,:) - - call phGLR(TDA,nS,Aph,Bph,EcRPA,eh_Om,XpY,XmY) call wall_time(end_t) t = end_t - start_t - write(*,'(A50,1X,F9.3,A8)') 'Wall time for phBSE =',t,' seconds' + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for phBSE problem =',t,' seconds' write(*,*) if(print_phLR) call print_excitation_energies('phBSE@Parquet','eh generalized',nS,eh_Om) @@ -228,10 +222,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, ! pp channel ! !-----------------! - write(*,*)' ------------------------------' - write(*,*)' | Diagonalizing ppBSE |' - write(*,*)' ------------------------------' - write(*,*) + write(*,*) 'Diagonalizing pp BSE problem...' + allocate(Bpp(nVV,nOO),Cpp(nVV,nVV),Dpp(nOO,nOO), & ee_Om(nVV),X1(nVV,nVV),Y1(nOO,nVV), & @@ -256,8 +248,8 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, else if(.not.TDA) call G_pp_Gamma_B(nOrb,nC,nO,nR,nOO,nVV,old_eh_Phi,pp_Gam_B) - call G_pp_Gamma_C(nOrb,nO,nR,nVV,old_eh_Phi,pp_Gam_C) - call G_pp_Gamma_D(nOrb,nC,nO,nOO,old_eh_Phi,pp_Gam_D) + call G_pp_Gamma_C(nOrb,nO,nR,nVV,old_eh_Phi,pp_Gam_C) + call G_pp_Gamma_D(nOrb,nC,nO,nOO,old_eh_Phi,pp_Gam_D) end if @@ -269,7 +261,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, call wall_time(end_t) t = end_t - start_t - write(*,'(A50,1X,F9.3,A8)') 'Wall time for ppBSE =',t,' seconds' + write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for ppBSE problem =',t,' seconds' write(*,*) if(print_ppLR) call print_excitation_energies('ppBSE@Parquet','2p generalized',nVV,ee_Om) @@ -277,18 +269,10 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, err_eig_ee = maxval(abs(old_ee_Om - ee_Om)) err_eig_hh = maxval(abs(old_hh_Om - hh_Om)) + err_eig_pp = max(err_eig_ee,err_eig_hh) deallocate(Bpp,Cpp,Dpp,pp_Gam_B,pp_Gam_C,pp_Gam_D) - - write(*,*) '----------------------------------------' - write(*,*) ' Two-body (eigenvalue) convergence ' - write(*,*) '----------------------------------------' - write(*,'(1X,A30,F10.6)')'Error for eh channel = ',err_eig_eh - write(*,'(1X,A30,F10.6)')'Error for pp channel = ',max(err_eig_ee,err_eig_hh) - write(*,*) '----------------------------------------' - write(*,*) - !----------! ! Updating ! !----------! @@ -344,7 +328,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, allocate(pp_Phi(nOrb,nOrb,nOrb,nOrb)) ! Build eh reducible kernels - write(*,*) 'Computing eh reducible kernel ...' + write(*,*) 'Computing eh reducible kernel...' call wall_time(start_t) call G_eh_Phi(nOrb,nC,nR,nS,old_eh_Om,eh_rho,eh_Phi) @@ -354,7 +338,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, write(*,*) ! Build pp reducible kernels - write(*,*) 'Computing pp reducible kernel ...' + write(*,*) 'Computing pp reducible kernel...' call wall_time(start_t) call G_pp_Phi(nOrb,nC,nR,nOO,nVV,old_ee_Om,ee_rho,old_hh_Om,hh_rho,pp_Phi) @@ -363,15 +347,15 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, write(*,'(1X,A50,1X,F9.3,A8)') 'Wall time for pp reducible kernel =',t,' seconds' write(*,*) - ! alpha = 0.05d0 - ! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:) - ! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:) - err_eh = maxval(abs(eh_Phi - old_eh_Phi)) err_pp = maxval(abs(pp_Phi - old_pp_Phi)) - call matout(nOrb**2,nOrb**2,eh_Phi - old_eh_Phi) - call matout(nOrb**2,nOrb**2,pp_Phi - old_pp_Phi) +! alpha = 0.05d0 +! eh_Phi(:,:,:,:) = alpha * eh_Phi(:,:,:,:) + (1d0 - alpha) * old_eh_Phi(:,:,:,:) +! pp_Phi(:,:,:,:) = alpha * pp_Phi(:,:,:,:) + (1d0 - alpha) * old_pp_Phi(:,:,:,:) + +! call matout(nOrb**2,nOrb**2,eh_Phi - old_eh_Phi) +! call matout(nOrb**2,nOrb**2,pp_Phi - old_pp_Phi) !--------------------! ! DIIS extrapolation ! @@ -425,12 +409,12 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, deallocate(eh_Phi,pp_Phi) - write(*,*) '----------------------------------------' - write(*,*) ' Two-body (kernel) convergence ' - write(*,*) '----------------------------------------' - write(*,'(1X,A30,F10.6)')'Error for eh channel = ',err_eh - write(*,'(1X,A30,F10.6)')'Error for pp channel = ',err_pp - write(*,*) '----------------------------------------' + write(*,*) '----------------------------------------------' + write(*,*) ' Two-body (frequency/kernel) convergence ' + write(*,*) '----------------------------------------------' + write(*,'(1X,A24,F10.6,1X,F10.6)')'Error for eh channel = ',err_eig_eh,err_eh + write(*,'(1X,A24,F10.6,1X,F10.6)')'Error for pp channel = ',err_eig_pp,err_pp + write(*,*) '----------------------------------------------' write(*,*) ! Convergence criteria @@ -438,7 +422,7 @@ subroutine GParquet(max_it_1b,conv_1b,max_it_2b,conv_2b,nOrb,nC,nO,nV,nR,nS,eHF, call wall_time(end_2b) t_2b = end_2b - start_2b - write(*,'(A50,1X,F9.3,A8)') 'Wall time for two-body iteration =',t_2b,' seconds' + write(*,'(A50,1X,I4,A2,F9.3,A8)') 'Wall time for two-body iteration #',n_it_2b,' =',t_2b,' seconds' write(*,*) end do diff --git a/src/Parquet/G_screened_integrals.f90 b/src/Parquet/G_screened_integrals.f90 index 9d40221..e921bbb 100644 --- a/src/Parquet/G_screened_integrals.f90 +++ b/src/Parquet/G_screened_integrals.f90 @@ -37,9 +37,9 @@ subroutine G_eh_screened_integral(nOrb,nC,nO,nR,nS,ERI,eh_Phi,pp_Phi,XpY,XmY,rho Y = 0.5d0*(XpY(ia,jb) - XmY(ia,jb)) rho(p,q,ia) = rho(p,q,ia) & - + (ERI(q,j,p,b) - ERI(q,j,b,p) & + + (1d0*ERI(q,j,p,b) - 1d0*ERI(q,j,b,p) & - 1d0*eh_Phi(q,j,b,p) + 1d0*pp_Phi(q,j,p,b)) * X & - + (ERI(q,b,p,j) - ERI(q,b,j,p) & + + (1d0*ERI(q,b,p,j) - 1d0*ERI(q,b,j,p) & - 1d0*eh_Phi(q,b,j,p) + 1d0*pp_Phi(q,b,p,j)) * Y end do @@ -107,7 +107,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do c=nO+1,nOrb-nR do d=c+1,nOrb-nR cd = cd + 1 - rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & + rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) & + 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X1(cd,ab) end do end do @@ -116,7 +116,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do k=nC+1,nO do l=k+1,nO kl = kl + 1 - rho1(p,q,ab) = rho1(p,q,ab) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & + rho1(p,q,ab) = rho1(p,q,ab) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) & + 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y1(kl,ab) end do end do @@ -132,7 +132,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do c=nO+1,nOrb-nR do d=c+1,nOrb-nR cd = cd + 1 - rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,c,d) - ERI(p,q,d,c) & + rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,c,d) - 1d0*ERI(p,q,d,c) & + 1d0*eh_Phi(p,q,c,d) - 1d0*eh_Phi(p,q,d,c) ) * X2(cd,ij) end do end do @@ -141,7 +141,7 @@ subroutine G_pp_screened_integral(nOrb,nC,nO,nR,nOO,nVV,ERI,eh_Phi,X1,Y1,rho1,X2 do k=nC+1,nO do l=k+1,nO kl = kl + 1 - rho2(p,q,ij) = rho2(p,q,ij) + ( ERI(p,q,k,l) - ERI(p,q,l,k) & + rho2(p,q,ij) = rho2(p,q,ij) + ( 1d0*ERI(p,q,k,l) - 1d0*ERI(p,q,l,k) & + 1d0*eh_Phi(p,q,k,l) - 1d0*eh_Phi(p,q,l,k) ) * Y2(kl,ij) end do end do