From cdb92208729a9db29e9f49fb35dac268cd1b4bd1 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 27 Nov 2023 13:30:53 +0100 Subject: [PATCH] cleanup qsGW --- src/GW/qsRGW.f90 | 49 +++++++++++++++++++++-------------------------- src/GW/ufG0W0.f90 | 6 +++--- src/GW/ufGW.f90 | 6 +++--- src/HF/RHF.f90 | 26 ++++++++++++------------- 4 files changed, 41 insertions(+), 46 deletions(-) diff --git a/src/GW/qsRGW.f90 b/src/GW/qsRGW.f90 index b6a74fe..26fcc6d 100644 --- a/src/GW/qsRGW.f90 +++ b/src/GW/qsRGW.f90 @@ -158,12 +158,9 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop nSCF = nSCF + 1 - ! Buid Hartree matrix + ! Build Hartree-exchange matrix call Hartree_matrix_AO_basis(nBas,P,ERI_AO,J) - - ! Compute exchange part of the self-energy - call exchange_matrix_AO_basis(nBas,P,ERI_AO,K) ! AO to MO transformation of two-electron integrals @@ -206,29 +203,6 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop if(nSCF > 1) Conv = maxval(abs(err)) - ! DIIS extrapolation - - if(max_diis > 1) then - n_diis = min(n_diis+1,max_diis) - call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F) - end if - - ! Diagonalize Hamiltonian in AO basis - - Fp = matmul(transpose(X),matmul(F,X)) - cp(:,:) = Fp(:,:) - call diagonalize_matrix(nBas,cp,eGW) - c = matmul(X,cp) - call AOtoMO(nBas,c,SigCp,SigC) - - ! Compute new density matrix in the AO basis - - P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) - - !------------------------------------------------------------------------ - ! Compute total energy - !------------------------------------------------------------------------ - ! Kinetic energy ET = trace_matrix(nBas,matmul(P,T)) @@ -249,6 +223,27 @@ subroutine qsRGW(dotest,maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,dop EqsGW = ET + EV + EJ + EK + ! DIIS extrapolation + + if(max_diis > 1) then + + n_diis = min(n_diis+1,max_diis) + call DIIS_extrapolation(rcond,nBasSq,nBasSq,n_diis,err_diis,F_diis,err,F) + + end if + + ! Diagonalize Hamiltonian in AO basis + + Fp = matmul(transpose(X),matmul(F,X)) + cp(:,:) = Fp(:,:) + call diagonalize_matrix(nBas,cp,eGW) + c = matmul(X,cp) + call AOtoMO(nBas,c,SigCp,SigC) + + ! Density matrix + + P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) + ! Print results call dipole_moment(nBas,P,nNuc,ZNuc,rNuc,dipole_int_AO,dipole) diff --git a/src/GW/ufG0W0.f90 b/src/GW/ufG0W0.f90 index a7935ac..faa6a4f 100644 --- a/src/GW/ufG0W0.f90 +++ b/src/GW/ufG0W0.f90 @@ -209,7 +209,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) timing = end_timing - start_timing write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for construction of supermatrix = ',timing,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds' write(*,*) else @@ -324,7 +324,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) timing = end_timing - start_timing write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for construction of supermatrix = ',timing,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds' write(*,*) end if @@ -341,7 +341,7 @@ subroutine ufG0W0(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) timing = end_timing - start_timing write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for diagonalization of supermatrix = ',timing,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for diagonalization of supermatrix = ',timing,' seconds' write(*,*) !-----------------! diff --git a/src/GW/ufGW.f90 b/src/GW/ufGW.f90 index be31c7c..00a1aa6 100644 --- a/src/GW/ufGW.f90 +++ b/src/GW/ufGW.f90 @@ -213,7 +213,7 @@ subroutine ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) timing = end_timing - start_timing write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for construction of supermatrix = ',timing,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds' write(*,*) else @@ -338,7 +338,7 @@ subroutine ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) timing = end_timing - start_timing write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for construction of supermatrix = ',timing,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for construction of supermatrix = ',timing,' seconds' write(*,*) end if @@ -355,7 +355,7 @@ subroutine ufGW(dotest,TDA_W,nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) timing = end_timing - start_timing write(*,*) - write(*,'(A65,1X,F9.3,A8)') 'Total wall time for diagonalization of supermatrix = ',timing,' seconds' + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for diagonalization of supermatrix = ',timing,' seconds' write(*,*) !-----------------! diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index eb99670..8b1e441 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -104,43 +104,43 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN do while(Conv > thresh .and. nSCF < maxSCF) -! Increment + ! Increment nSCF = nSCF + 1 -! Build Fock matrix + ! Build Fock matrix call Hartree_matrix_AO_basis(nBas,P,ERI,J) call exchange_matrix_AO_basis(nBas,P,ERI,K) F(:,:) = Hc(:,:) + J(:,:) + 0.5d0*K(:,:) -! Check convergence + ! Check convergence err = matmul(F,matmul(P,S)) - matmul(matmul(S,P),F) if(nSCF > 1) Conv = maxval(abs(err)) -! Kinetic energy + ! Kinetic energy ET = trace_matrix(nBas,matmul(P,T)) -! Potential energy + ! Potential energy EV = trace_matrix(nBas,matmul(P,V)) -! Hartree energy + ! Hartree energy EJ = 0.5d0*trace_matrix(nBas,matmul(P,J)) -! Exchange energy + ! Exchange energy EK = 0.25d0*trace_matrix(nBas,matmul(P,K)) -! Total energy + ! Total energy ERHF = ET + EV + EJ + EK -! DIIS extrapolation + ! DIIS extrapolation if(max_diis > 1) then @@ -149,22 +149,22 @@ subroutine RHF(dotest,maxSCF,thresh,max_diis,guess_type,level_shift,nNuc,ZNuc,rN end if -! Level-shifting + ! Level shift if(level_shift > 0d0 .and. Conv > thresh) call level_shifting(level_shift,nBas,nO,S,c,F) -! Diagonalize Fock matrix + ! Diagonalize Fock matrix Fp = matmul(transpose(X),matmul(F,X)) cp(:,:) = Fp(:,:) call diagonalize_matrix(nBas,cp,eHF) c = matmul(X,cp) -! Density matrix + ! Density matrix P(:,:) = 2d0*matmul(c(:,1:nO),transpose(c(:,1:nO))) -! Dump results + ! Dump results write(*,'(1X,A1,1X,I3,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,F16.10,1X,A1,1X,E10.2,1X,A1,1X)') & '|',nSCF,'|',ERHF + ENuc,'|',EJ,'|',EK,'|',Conv,'|'