diff --git a/input/dft b/input/dft index ed3a55d..b16601a 100644 --- a/input/dft +++ b/input/dft @@ -19,25 +19,25 @@ # Number of states in ensemble (nEns) 4 # occupation numbers -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 +0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 1.0 0. 0.0 + 1 0.0 0.0 # Ncentered ? F # Parameters for CC weight-dependent exchange functional 4 -0.60601 -0.0631565 -0.0289751 0.00244785 --1.28842 -0.173117 0.0900511 -0.0118975 +-0.718713,-0.133321,0.226288,-0.250718 +-0.525899,0.687216,-0.13866,-0.0226579 0.0 0.0 0.0 0.0 # choice of UCC exchange coefficient : 1 for Cx1, 2 for Cx2, 3 for Cx1*Cx2 1 diff --git a/input/methods b/input/methods index 8124aba..979b54d 100644 --- a/input/methods +++ b/input/methods @@ -1,5 +1,5 @@ # RHF UHF KS MOM - F F T F + F F T F # MP2* MP3 MP2-F12 F F F # CCD pCCD DCD CCSD CCSD(T) @@ -11,9 +11,9 @@ # RPA* RPAx* crRPA ppRPA F F F F # G0F2* evGF2* qsGF2* G0F3 evGF3 - F F F F F + F F F F F # G0W0* evGW* qsGW* ufG0W0 ufGW - F F F F F + F F F F F # G0T0 evGT qsGT F F F # MCMP2 diff --git a/input/options b/input/options index 2d0e99d..f4a2779 100644 --- a/input/options +++ b/input/options @@ -1,11 +1,11 @@ # HF: maxSCF thresh DIIS n_diis guess_type ortho_type mix_guess stability - 1024 0.00001 T 5 1 1 T F + 1024 0.00001 T 2 1 1 F F # MP: # CC: maxSCF thresh DIIS n_diis 64 0.00001 T 5 # spin: TDA singlet triplet spin_conserved spin_flip - F T F T T + F T T T T # GF: maxSCF thresh DIIS n_diis lin eta renorm reg 256 0.00001 T 5 T 0.0 3 F # GW: maxSCF thresh DIIS n_diis lin eta COHSEX SOSEX TDA_W G0W GW0 reg diff --git a/mol/h2.xyz b/mol/h2.xyz index a4e936a..6a4e902 100644 --- a/mol/h2.xyz +++ b/mol/h2.xyz @@ -1,4 +1,4 @@ 2 - -H 0. 0. 0. -H 0. 0. 2.000000 + +H 0.0 0.0 0.0 +H 0.0 0.0 0.7 diff --git a/src/CI/CID.f90 b/src/CI/CID.f90 index ca6a78d..b8f52d6 100644 --- a/src/CI/CID.f90 +++ b/src/CI/CID.f90 @@ -96,7 +96,7 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi H(ishift+1,jshift+1) = E0 - print*,'00 block done...' + print*,'00 block done...' ! 0D blocks @@ -193,7 +193,7 @@ subroutine CID(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERIi end do end do - print*,'DD block done...' + print*,'DD block done...' write(*,*) write(*,*) 'Diagonalizing CID matrix...' diff --git a/src/CI/CISD.f90 b/src/CI/CISD.f90 index 1037d11..079521c 100644 --- a/src/CI/CISD.f90 +++ b/src/CI/CISD.f90 @@ -86,7 +86,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI write(*,*) 'nH = ',nH write(*,*) - maxH = min(nH,21) + maxH = min(nH,41) ! Memory allocation @@ -99,7 +99,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI H(ishift+1,jshift+1) = E0 - print*,'00 block done...' + print*,'00 block done...' ! 0S blocks @@ -172,7 +172,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI end do end do - print*,'SS block done...' + print*,'SS block done...' ! SD blocks @@ -285,7 +285,7 @@ subroutine CISD(singlet_manifold,triplet_manifold,nBasin,nCin,nOin,nVin,nRin,ERI end do end do - print*,'DD block done...' + print*,'DD block done...' write(*,*) write(*,*) 'Diagonalizing CISD matrix...' diff --git a/src/GF/evGF2.f90 b/src/GF/evGF2.f90 index 9466f53..e0fbeef 100644 --- a/src/GF/evGF2.f90 +++ b/src/GF/evGF2.f90 @@ -70,6 +70,7 @@ subroutine evGF2(BSE,TDA,dBSE,dTDA,evDyn,maxSCF,thresh,max_diis,singlet,triplet, error_diis(:,:) = 0d0 eGF2(:) = eHF(:) eOld(:) = eHF(:) + rcond = 0d0 !------------------------------------------------------------------------ ! Main SCF loop diff --git a/src/GF/evUGF2.f90 b/src/GF/evUGF2.f90 index 3d7c844..7d8c635 100644 --- a/src/GF/evUGF2.f90 +++ b/src/GF/evUGF2.f90 @@ -102,6 +102,7 @@ subroutine evUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, eGF2(:,:) = eHF(:,:) eOld(:,:) = eHF(:,:) Z(:,:) = 1d0 + rcond(:) = 0d0 !------------------------------------------------------------------------ ! Main loop diff --git a/src/GF/qsGF2.f90 b/src/GF/qsGF2.f90 index 6a5a395..350a81c 100644 --- a/src/GF/qsGF2.f90 +++ b/src/GF/qsGF2.f90 @@ -94,6 +94,8 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, nBasSq = nBas*nBas + print*,maxSCF + ! TDA if(TDA) then @@ -119,7 +121,7 @@ subroutine qsGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,singlet,triplet, c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - rcond = 1d0 + rcond = 0d0 !------------------------------------------------------------------------ ! Main loop diff --git a/src/GF/qsUGF2.f90 b/src/GF/qsUGF2.f90 index e94a344..7a53c1d 100644 --- a/src/GF/qsUGF2.f90 +++ b/src/GF/qsUGF2.f90 @@ -135,7 +135,7 @@ subroutine qsUGF2(maxSCF,thresh,max_diis,BSE,TDA,dBSE,dTDA,evDyn,spin_conserved, c(:,:,:) = cHF(:,:,:) F_diis(:,:,:) = 0d0 error_diis(:,:,:) = 0d0 - rcond = 1d0 + rcond(:) = 0d0 !------------------------------------------------------------------------ ! Main loop diff --git a/src/GT/Bethe_Salpeter_Tmatrix.f90 b/src/GT/Bethe_Salpeter_Tmatrix.f90 index 05e019e..b038740 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix.f90 @@ -58,7 +58,8 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, integer :: iblock double precision :: EcRPA(nspin) - double precision,allocatable :: TA(:,:),TB(:,:) + double precision,allocatable :: TAs(:,:),TBs(:,:) + double precision,allocatable :: TAt(:,:),TBt(:,:) double precision,allocatable :: OmBSE(:,:) double precision,allocatable :: XpY_BSE(:,:,:) double precision,allocatable :: XmY_BSE(:,:,:) @@ -69,16 +70,12 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Memory allocation - allocate(TA(nS,nS),TB(nS,nS),OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) + allocate(TAs(nS,nS),TBs(nS,nS),TAt(nS,nS),TBt(nS,nS), & + OmBSE(nS,nspin),XpY_BSE(nS,nS,nspin),XmY_BSE(nS,nS,nspin)) -! Initialize T matrix - - TA(:,:) = 0d0 - TB(:,:) = 0d0 - -!---------------------------------------------- -! Compute T-matrix for alpha-beta block -!---------------------------------------------- +!---------------------------------------! +! Compute T-matrix for alpha-beta block ! +!---------------------------------------! ispin = 1 iblock = 3 @@ -88,17 +85,17 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TA) - if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,ERI,Omega1s,rho1s,Omega2s,rho2s,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Omega1s,rho1s,Omega2s,rho2s,TAs) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,Omega1s,rho1s,Omega2s,rho2s,TBs) -! print*,'aa block of TA' -! call matout(nS,nS,TA) -! print*,'aa block of TB' -! call matout(nS,nS,TB) +! print*,'ab block of TA' +! call matout(nS,nS,TAs) +! print*,'ab block of TB' +! call matout(nS,nS,TBs) -!---------------------------------------------- -! Compute T-matrix for alpha-alpha block -!---------------------------------------------- +!----------------------------------------! +! Compute T-matrix for alpha-alpha block ! +!----------------------------------------! ispin = 2 iblock = 4 @@ -108,17 +105,17 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TA) - if(.not.TDA) call static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,ERI,Omega1t,rho1t,Omega2t,rho2t,TB) + call static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Omega1t,rho1t,Omega2t,rho2t,TAt) + if(.not.TDA) call static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,Omega1t,rho1t,Omega2t,rho2t,TBt) -! print*,'aa+ab block of TA' -! call matout(nS,nS,TA) -! print*,'aa+ab block of TB' -! call matout(nS,nS,TB) +! print*,'aa block of TA' +! call matout(nS,nS,TAt) +! print*,'aa block of TB' +! call matout(nS,nS,TBt) -!------------------- -! Singlet manifold -!------------------- +!------------------! +! Singlet manifold ! +!------------------! if(singlet) then @@ -127,39 +124,37 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE singlet excitation energies - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & + call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt+TAs,TBt+TBs, & EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.true.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) - !------------------------------------------------- - ! Compute the dynamical screening at the BSE level - !------------------------------------------------- - if(dBSE) then - + ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - + if(evDyn) then - + print*, ' Iterative dynamical correction for BSE@GT NYI' ! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & ! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) else - - call Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,Omega1s,Omega2s,rho1s,rho2s, & - eT,eGT,dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + + call Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Omega1s,Omega2s,Omega1t,Omega2t,rho1s,rho2s,rho1t,rho2t,eT,eGT, & + dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin), & + TAs,TBs,TAt,TBt) end if - + end if end if -!------------------- -! Triplet manifold -!------------------- +!------------------! +! Triplet manifold ! +!------------------! if(triplet) then @@ -168,32 +163,29 @@ subroutine Bethe_Salpeter_Tmatrix(TDA_T,TDA,dBSE,dTDA,evDyn,singlet,triplet,eta, ! Compute BSE triplet excitation energies - call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TA,TB, & + call linear_response_Tmatrix(ispin,.false.,TDA,eta,nBas,nC,nO,nV,nR,nS,1d0,eGT,ERI,TAt-TAs,TBt-TBs, & EcBSE(ispin),OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) call print_excitation('BSE@GT ',ispin,nS,OmBSE(:,ispin)) call print_transition_vectors(.false.,nBas,nC,nO,nV,nR,nS,dipole_int, & OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) - !------------------------------------------------- - ! Compute the dynamical screening at the BSE level - !------------------------------------------------- - if(dBSE) then - + ! Compute dynamic correction for BSE via perturbation theory (iterative or renormalized) - + if(evDyn) then - - + print*, ' Iterative dynamical correction for BSE@GT NYI' ! call Bethe_Salpeter_dynamic_perturbation_iterative(dTDA,eta,nBas,nC,nO,nV,nR,nS,eGW,dipole_int,OmRPA,rho_RPA, & ! OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) else - - call Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,Omega1t,Omega2t,rho1t,rho2t, & - eT,eGT,dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin)) + + call Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Omega1s,Omega2s,Omega1t,Omega2t,rho1s,rho2s,rho1t,rho2t,eT,eGT, & + dipole_int,OmBSE(:,ispin),XpY_BSE(:,:,ispin),XmY_BSE(:,:,ispin), & + TAs,TBs,TAt,TBt) end if - + end if end if diff --git a/src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 b/src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 index 75eae8b..e1408b3 100644 --- a/src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 +++ b/src/GT/Bethe_Salpeter_Tmatrix_dynamic_perturbation.f90 @@ -1,5 +1,7 @@ -subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,Omega1,Omega2,rho1,rho2, & - eT,eGT,dipole_int,OmBSE,XpY,XmY) +subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(ispin,dTDA,eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,nOOt,nVVt, & + Omega1s,Omega2s,Omega1t,Omega2t,rho1s,rho2s,rho1t,rho2t,eT,eGT, & + dipole_int,OmBSE,XpY,XmY,TAs,TBs,TAt,TBt) + ! Compute dynamical effects via perturbation theory for BSE@GT implicit none @@ -7,6 +9,7 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR ! Input variables + integer,intent(in) :: ispin logical,intent(in) :: dTDA double precision,intent(in) :: eta integer,intent(in) :: nBas @@ -16,8 +19,10 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR integer,intent(in) :: nR integer,intent(in) :: nS - integer,intent(in) :: nOO - integer,intent(in) :: nVV + integer,intent(in) :: nOOs + integer,intent(in) :: nVVs + integer,intent(in) :: nOOt + integer,intent(in) :: nVVt double precision,intent(in) :: eT(nBas) double precision,intent(in) :: eGT(nBas) @@ -26,16 +31,25 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR double precision,intent(in) :: XpY(nS,nS) double precision,intent(in) :: XmY(nS,nS) - double precision,intent(in) :: Omega1(nVV) - double precision,intent(in) :: Omega2(nOO) - double precision,intent(in) :: rho1(nBas,nBas,nVV) - double precision,intent(in) :: rho2(nBas,nBas,nOO) + double precision,intent(in) :: Omega1s(nVVs) + double precision,intent(in) :: Omega2s(nOOs) + double precision,intent(in) :: rho1s(nBas,nBas,nVVs) + double precision,intent(in) :: rho2s(nBas,nBas,nOOs) + double precision,intent(in) :: Omega1t(nVVt) + double precision,intent(in) :: Omega2t(nOOt) + double precision,intent(in) :: rho1t(nBas,nBas,nVVt) + double precision,intent(in) :: rho2t(nBas,nBas,nOOt) + + double precision,intent(in) :: TAs(nS,nS) + double precision,intent(in) :: TBs(nS,nS) + double precision,intent(in) :: TAt(nS,nS) + double precision,intent(in) :: TBt(nS,nS) ! Local variables integer :: ia - integer,parameter :: maxS = 10 + integer :: maxS = 10 double precision :: gapGT double precision,allocatable :: OmDyn(:) @@ -43,30 +57,63 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR double precision,allocatable :: X(:) double precision,allocatable :: Y(:) - double precision,allocatable :: Ap_dyn(:,:) - double precision,allocatable :: ZAp_dyn(:,:) + double precision,allocatable :: dTAs(:,:) + double precision,allocatable :: ZAs(:,:) - double precision,allocatable :: Bp_dyn(:,:) - double precision,allocatable :: ZBp_dyn(:,:) - - double precision,allocatable :: Am_dyn(:,:) - double precision,allocatable :: ZAm_dyn(:,:) - - double precision,allocatable :: Bm_dyn(:,:) - double precision,allocatable :: ZBm_dyn(:,:) + double precision,allocatable :: dTAt(:,:) + double precision,allocatable :: ZAt(:,:) ! Memory allocation - allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) - - if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) + maxS = min(nS,maxS) + allocate(OmDyn(maxS),ZDyn(maxS),X(nS),Y(nS),dTAs(nS,nS),ZAs(nS,nS),dTAt(nS,nS),ZAt(nS,nS)) if(dTDA) then write(*,*) write(*,*) '*** dynamical TDA activated ***' write(*,*) + else + print*, ' Beyond-TDA dynamical correction for BSE@GT NYI' + return end if + OmDyn(:) = 0d0 + ZDyn(:) = 0d0 + + do ia=1,maxS + + ! Compute dynamical T-matrix for alpha-beta block + + call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOs,nVVs,1d0,eGT,Omega1s,Omega2s,rho1s,rho2s,OmBSE(ia),dTAs,ZAs) + + ! Compute dynamical T-matrix for alpha-beta block + + call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOOt,nVVt,1d0,eGT,Omega1t,Omega2t,rho1t,rho2t,OmBSE(ia),dTAt,ZAt) + + X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) + Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) + + ! First-order correction + + if(ispin == 1) then + ZDyn(ia) = - dot_product(X,matmul(ZAt+ZAs,X)) + OmDyn(ia) = - dot_product(X,matmul(dTAt+dTAs,X)) + dot_product(X,matmul(TAt+TAs,X)) + end if + + if(ispin == 2) then + ZDyn(ia) = - dot_product(X,matmul(ZAt-ZAs,X)) + OmDyn(ia) = - dot_product(X,matmul(dTAt-dTAs,X)) + dot_product(X,matmul(TAt-TAs,X)) + end if + + ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) + OmDyn(ia) = ZDyn(ia)*OmDyn(ia) + + end do + +!--------------! +! Dump results ! +!--------------! + gapGT = eGT(nO+1) - eGT(nO) write(*,*) '---------------------------------------------------------------------------------------------------' @@ -77,54 +124,11 @@ subroutine Bethe_Salpeter_Tmatrix_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' - do ia=1,min(nS,maxS) - - X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) - Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) - - ! First-order correction - - if(dTDA) then - - ! Resonant part of the BSE correction for dynamical TDA - - call dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia),Ap_dyn,Zap_dyn) - - ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) - OmDyn(ia) = dot_product(X,matmul( Ap_dyn,X)) - - else - - print*, ' Beyond-TDA dynamical correction for BSE@GT NYI' - ! Resonant and anti-resonant part of the BSE correction - -! call dynamic_Tmatrix_TAB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia), & -! Ap_dyn,Am_dyn,Bp_dyn,Bm_dyn) - - ! Renormalization factor of the resonant and anti-resonant parts - -! call dynamic_Tmatrix_ZAB(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,1d0,eGT,Omega1,Omega2,rho1,rho2,OmBSE(ia), & -! ZAp_dyn,ZAm_dyn,ZBp_dyn,ZBm_dyn) - - ZDyn(ia) = dot_product(X,matmul(ZAp_dyn,X)) & - - dot_product(Y,matmul(ZAm_dyn,Y)) & - + dot_product(X,matmul(ZBp_dyn,Y)) & - - dot_product(Y,matmul(ZBm_dyn,X)) - - OmDyn(ia) = dot_product(X,matmul(Ap_dyn,X)) & - - dot_product(Y,matmul(Am_dyn,Y)) & - + dot_product(X,matmul(Bp_dyn,Y)) & - - dot_product(Y,matmul(Bm_dyn,X)) - - end if - - ZDyn(ia) = 1d0/(1d0 - ZDyn(ia)) - OmDyn(ia) = ZDyn(ia)*OmDyn(ia) - + do ia=1,maxS write(*,'(2X,I5,5X,F15.6,5X,F15.6,5X,F15.6,5X,F15.6)') & ia,OmBSE(ia)*HaToeV,(OmBSE(ia)+OmDyn(ia))*HaToeV,OmDyn(ia)*HaToeV,ZDyn(ia) - end do + write(*,*) '---------------------------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/G0T0.f90 b/src/GT/G0T0.f90 index 94f886c..11f4415 100644 --- a/src/GT/G0T0.f90 +++ b/src/GT/G0T0.f90 @@ -48,6 +48,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcAC(nspin) + double precision :: EcGM double precision,allocatable :: Omega1s(:),Omega1t(:) double precision,allocatable :: X1s(:,:),X1t(:,:) double precision,allocatable :: Y1s(:,:),Y1t(:,:) @@ -60,11 +61,6 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing double precision,allocatable :: SigT(:) double precision,allocatable :: Z(:) - double precision,allocatable :: Omega(:,:) - double precision,allocatable :: XpY(:,:,:) - double precision,allocatable :: XmY(:,:,:) - double precision,allocatable :: rho(:,:,:,:) - ! Output variables double precision,intent(out) :: eG0T0(nBas) @@ -81,6 +77,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing nOOs = nO*nO nVVs = nV*nV +! nOOs = nO*(nO + 1)/2 +! nVVs = nV*(nV + 1)/2 nOOt = nO*(nO - 1)/2 nVVt = nV*(nV - 1)/2 @@ -101,6 +99,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ispin = 1 iblock = 3 +! iblock = 1 ! Compute linear response @@ -109,8 +108,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! EcRPA(ispin) = 1d0*EcRPA(ispin) -! call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:)) -! call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:)) + call print_excitation('pp-RPA (N+2)',iblock,nVVs,Omega1s(:)) + call print_excitation('pp-RPA (N-2)',iblock,nOOs,Omega2s(:)) !---------------------------------------------- ! alpha-alpha block @@ -127,35 +126,52 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ! EcRPA(ispin) = 2d0*EcRPA(ispin) ! EcRPA(ispin) = 3d0*EcRPA(ispin) -! call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:)) -! call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:)) + call print_excitation('pp-RPA (N+2)',iblock,nVVt,Omega1t(:)) + call print_excitation('pp-RPA (N-2)',iblock,nOOt,Omega2t(:)) !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- + EcGM = 0d0 SigT(:) = 0d0 Z(:) = 0d0 - iblock = 3 + iblock = 3 +! iblock = 1 call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO,X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,SigT) + if(regularize) then - call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z) + call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z) - iblock = 4 + else + + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eHF,Omega1s,rho1s,Omega2s,rho2s,Z) + + end if + + iblock = 4 call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO,X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,SigT) + if(regularize) then - call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,Z) + call regularized_self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,Z) + + else + + call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eHF,Omega1t,rho1t,Omega2t,rho2t,Z) + + end if Z(:) = 1d0/(1d0 - Z(:)) - !---------------------------------------------- ! Compute the exchange part of the self-energy !---------------------------------------------- @@ -184,6 +200,8 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing ispin = 1 iblock = 3 +! iblock = 1 + call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eG0T0,ERI_MO, & Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) ispin = 2 @@ -194,7 +212,7 @@ subroutine G0T0(doACFDT,exchange_kernel,doXBS,BSE,TDA_T,TDA,dBSE,dTDA,evDyn,sing EcRPA(1) = EcRPA(1) - EcRPA(2) EcRPA(2) = 3d0*EcRPA(2) - call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcRPA) + call print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eG0T0,EcGM,EcRPA) ! Perform BSE calculation diff --git a/src/GT/dynamic_Tmatrix_A.f90 b/src/GT/dynamic_Tmatrix_A.f90 index e6770a0..7476be2 100644 --- a/src/GT/dynamic_Tmatrix_A.f90 +++ b/src/GT/dynamic_Tmatrix_A.f90 @@ -1,4 +1,4 @@ -subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,A_dyn,ZA_dyn) +subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,TA,ZA) ! Compute the dynamic part of the Bethe-Salpeter equation matrices for GT @@ -36,13 +36,13 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O ! Output variables - double precision,intent(out) :: A_dyn(nS,nS) - double precision,intent(out) :: ZA_dyn(nS,nS) + double precision,intent(out) :: TA(nS,nS) + double precision,intent(out) :: ZA(nS,nS) ! Initialization - A_dyn(:,:) = 0d0 - ZA_dyn(:,:) = 0d0 + TA(:,:) = 0d0 + ZA(:,:) = 0d0 ! Build dynamic A matrix @@ -58,44 +58,30 @@ subroutine dynamic_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,O chi = 0d0 do cd=1,nVV - eps = - Omega1(cd) - chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*eps/(eps**2 + eta**2) + eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) + chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2) end do do kl=1,nOO - eps = + Omega2(kl) - chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2) + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) + chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2) end do - A_dyn(ia,jb) = A_dyn(ia,jb) - lambda*chi + TA(ia,jb) = TA(ia,jb) - lambda*chi chi = 0d0 do cd=1,nVV eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) - chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*eps/(eps**2 + eta**2) + chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do do kl=1,nOO eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) - chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*eps/(eps**2 + eta**2) + chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 end do - A_dyn(ia,jb) = A_dyn(ia,jb) + 1d0*lambda*chi - - chi = 0d0 - - do cd=1,nVV - eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(j)) - chi = chi + rho1(i,b,cd)*rho1(j,a,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do - - do kl=1,nOO - eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(b)) - chi = chi + rho2(i,b,kl)*rho2(j,a,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 - end do - - ZA_dyn(ia,jb) = ZA_dyn(ia,jb) - 1d0*lambda*chi + ZA(ia,jb) = ZA(ia,jb) + lambda*chi end do end do diff --git a/src/GT/dynamic_Tmatrix_B.f90 b/src/GT/dynamic_Tmatrix_B.f90 new file mode 100644 index 0000000..6d67b78 --- /dev/null +++ b/src/GT/dynamic_Tmatrix_B.f90 @@ -0,0 +1,91 @@ +subroutine dynamic_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,eGT,Omega1,Omega2,rho1,rho2,OmBSE,TB,ZB) + +! Compute the off-diagonal dynamic part of the Bethe-Salpeter equation matrices for GT + + implicit none + include 'parameters.h' + +! Input variables + + double precision,intent(in) :: eta + integer,intent(in) :: nBas + integer,intent(in) :: nC + integer,intent(in) :: nO + integer,intent(in) :: nV + integer,intent(in) :: nR + integer,intent(in) :: nS + + integer,intent(in) :: nOO + integer,intent(in) :: nVV + + double precision,intent(in) :: lambda + double precision,intent(in) :: eGT(nBas) + double precision,intent(in) :: OmBSE + + + double precision,intent(in) :: Omega1(nVV) + double precision,intent(in) :: Omega2(nOO) + double precision,intent(in) :: rho1(nBas,nBas,nVV) + double precision,intent(in) :: rho2(nBas,nBas,nOO) + +! Local variables + + double precision :: chi + double precision :: eps + integer :: i,j,a,b,ia,jb,cd,kl + +! Output variables + + double precision,intent(out) :: TB(nS,nS) + double precision,intent(out) :: ZB(nS,nS) + +! Initialization + + TB(:,:) = 0d0 + ZB(:,:) = 0d0 + +! Build dynamic A matrix + + ia = 0 + do i=nC+1,nO + do a=nO+1,nBas-nR + ia = ia + 1 + jb = 0 + do j=nC+1,nO + do b=nO+1,nBas-nR + jb = jb + 1 + + chi = 0d0 + + do cd=1,nVV + eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(b)) + chi = chi + rho1(i,j,cd)*rho1(b,a,cd)*eps/(eps**2 + eta**2) + end do + + do kl=1,nOO + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(j)) + chi = chi + rho2(i,j,kl)*rho2(b,a,kl)*eps/(eps**2 + eta**2) + end do + + TB(ia,jb) = TB(ia,jb) + lambda*chi + + chi = 0d0 + + do cd=1,nVV + eps = + OmBSE - Omega1(cd) + (eGT(i) + eGT(b)) + chi = chi + rho1(i,j,cd)*rho1(b,a,cd)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + + do kl=1,nOO + eps = + OmBSE + Omega2(kl) - (eGT(a) + eGT(j)) + chi = chi + rho2(i,a,kl)*rho2(b,a,kl)*(eps**2 - eta**2)/(eps**2 + eta**2)**2 + end do + + ZB(ia,jb) = ZB(ia,jb) - lambda*chi + + end do + end do + end do + end do + +end subroutine dynamic_Tmatrix_B diff --git a/src/GT/evGT.f90 b/src/GT/evGT.f90 index 65dcfa4..1c5e518 100644 --- a/src/GT/evGT.f90 +++ b/src/GT/evGT.f90 @@ -46,7 +46,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Local variables - logical :: linear_mixing integer :: nSCF integer :: n_diis double precision :: rcond @@ -55,6 +54,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & integer :: iblock integer :: nOOs,nOOt integer :: nVVs,nVVt + double precision :: EcGM double precision :: EcRPA(nspin) double precision :: EcBSE(nspin) double precision :: EcAC(nspin) @@ -117,6 +117,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & eGT(:) = eG0T0(:) eOld(:) = eGT(:) Z(:) = 1d0 + rcond = 0d0 !------------------------------------------------------------------------ ! Main loop @@ -148,10 +149,14 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) + EcRPA(1) = EcRPA(1) - EcRPA(2) + EcRPA(2) = 3d0*EcRPA(2) + !---------------------------------------------- ! Compute T-matrix version of the self-energy !---------------------------------------------- + EcGM = 0d0 SigT(:) = 0d0 Z(:) = 0d0 @@ -161,7 +166,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & X1s,Y1s,rho1s,X2s,Y2s,rho2s) call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Omega1s,rho1s,Omega2s,rho2s,SigT) + Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & Omega1s,rho1s,Omega2s,rho2s,Z) @@ -172,7 +177,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & X1t,Y1t,rho1t,X2t,Y2t,rho2t) call self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Omega1t,rho1t,Omega2t,rho2t,SigT) + Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & Omega1t,rho1t,Omega2t,rho2t,Z) @@ -195,7 +200,7 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! Dump results !---------------------------------------------- - call print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) + call print_evGT(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) ! DIIS extrapolation @@ -219,29 +224,6 @@ subroutine evGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS, & ! End main loop !------------------------------------------------------------------------ -! Compute the ppRPA correlation energy - - ispin = 1 - iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) - ispin = 2 - iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) - EcRPA(1) = EcRPA(1) - EcRPA(2) - EcRPA(2) = 3d0*EcRPA(2) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@evGT total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - - ! Perform BSE calculation if(BSE) then diff --git a/src/GT/excitation_density_Tmatrix.f90 b/src/GT/excitation_density_Tmatrix.f90 index 01ae36c..fa58d75 100644 --- a/src/GT/excitation_density_Tmatrix.f90 +++ b/src/GT/excitation_density_Tmatrix.f90 @@ -22,8 +22,8 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r ! Local variables - integer :: k,l - integer :: c,d + integer :: i,j,k,l + integer :: a,b,c,d integer :: p,q integer :: ab,cd,ij,kl double precision,external :: Kronecker_delta @@ -47,45 +47,71 @@ subroutine excitation_density_Tmatrix(ispin,nBas,nC,nO,nV,nR,nOO,nVV,ERI,X1,Y1,r do p=nC+1,nBas-nR do q=nC+1,nBas-nR - do ab=1,nVV +! do ab=1,nVV + ab = 0 + do a=nO+1,nBas-nR + do b=a,nBas-nR + ab = ab + 1 cd = 0 do c=nO+1,nBas-nR +! do d=nO+1,c do d=c,nBas-nR cd = cd + 1 - rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,c,d)*X1(cd,ab) + rho1(p,q,ab) = rho1(p,q,ab) & + + (1d0*ERI(p,q,c,d) + 0d0*ERI(p,q,d,c))*X1(cd,ab) +! + ERI(p,q,c,d)*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(c,d))) +! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X1(cd,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(c,d))) end do end do kl = 0 do k=nC+1,nO +! do l=nC+1,k do l=k,nO kl = kl + 1 - rho1(p,q,ab) = rho1(p,q,ab) + ERI(p,q,k,l)*Y1(kl,ab) + rho1(p,q,ab) = rho1(p,q,ab) & + + (1d0*ERI(p,q,k,l) + 0d0*ERI(p,q,l,k))*Y1(kl,ab) +! + ERI(p,q,k,l)*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(k,l))) +! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y1(kl,ab)/sqrt((1d0 + Kronecker_delta(a,b))*(1d0 + Kronecker_delta(k,l))) end do end do end do + end do - do ij=1,nOO +! do ij=1,nOO + ij = 0 + do i=nC+1,nO + do j=i,nO + ij = ij + 1 cd = 0 do c=nO+1,nBas-nR +! do d=nO+1,c do d=c,nBas-nR cd = cd + 1 - rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,c,d)*X2(cd,ij) + rho2(p,q,ij) = rho2(p,q,ij) & + + (1d0*ERI(p,q,c,d) + 0d0*ERI(p,q,d,c))*X2(cd,ij) +! + ERI(p,q,c,d)*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(c,d))) +! + (ERI(p,q,c,d) + ERI(p,q,d,c))*X2(cd,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(c,d))) end do end do kl = 0 do k=nC+1,nO +! do l=nC+1,k do l=k,nO kl = kl + 1 - rho2(p,q,ij) = rho2(p,q,ij) + ERI(p,q,k,l)*Y2(kl,ij) + rho2(p,q,ij) = rho2(p,q,ij) & + + (1d0*ERI(p,q,k,l) + 0d0*ERI(p,q,l,k))*Y2(kl,ij) +! + ERI(p,q,k,l)*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(k,l))) +! + (ERI(p,q,k,l) + ERI(p,q,l,k))*Y2(kl,ij)/sqrt((1d0 + Kronecker_delta(i,j))*(1d0 + Kronecker_delta(k,l))) end do end do end do + end do end do end do diff --git a/src/GT/print_G0T0.f90 b/src/GT/print_G0T0.f90 index 2150cd4..95b4949 100644 --- a/src/GT/print_G0T0.f90 +++ b/src/GT/print_G0T0.f90 @@ -1,13 +1,15 @@ -subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcRPA) +subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) ! Print one-electron energies and other stuff for G0T0 implicit none include 'parameters.h' - integer,intent(in) :: nBas,nO + integer,intent(in) :: nBas + integer,intent(in) :: nO double precision,intent(in) :: ENuc double precision,intent(in) :: ERHF + double precision,intent(in) :: EcGM double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: eHF(nBas) double precision,intent(in) :: SigT(nBas) @@ -50,10 +52,12 @@ subroutine print_G0T0(nBas,nO,eHF,ENuc,ERHF,SigT,Z,eGT,EcRPA) write(*,'(2X,A50,F15.6,A3)') 'G0T0 LUMO energy (eV) =',eGT(LUMO)*HaToeV,' eV' write(*,'(2X,A50,F15.6,A3)') 'G0T0 HOMO-LUMO gap (eV) =',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' - write(*,'(2X,A50,F20.10,A3)') 'Tr@RPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 correlation energy =',EcGM,' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 total energy =',ENuc + ERHF + EcGM,' au' write(*,*)'-------------------------------------------------------------------------------' write(*,*) diff --git a/src/GT/print_evGT.f90 b/src/GT/print_evGT.f90 index 9b15356..7a11544 100644 --- a/src/GT/print_evGT.f90 +++ b/src/GT/print_evGT.f90 @@ -1,4 +1,4 @@ -subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) +subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,ENuc,ERHF,SigT,Z,eGT,EcGM,EcRPA) ! Print one-electron energies and other stuff for evGT @@ -10,6 +10,10 @@ subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) integer,intent(in) :: nSCF double precision,intent(in) :: Conv double precision,intent(in) :: eHF(nBas) + double precision,intent(in) :: ENuc + double precision,intent(in) :: ERHF + double precision,intent(in) :: EcGM + double precision,intent(in) :: EcRPA(nspin) double precision,intent(in) :: SigT(nBas) double precision,intent(in) :: Z(nBas) double precision,intent(in) :: eGT(nBas) @@ -49,6 +53,13 @@ subroutine print_evGT(nBas,nO,nSCF,Conv,eHF,SigT,Z,eGT) write(*,'(2X,A30,F15.6)') 'evGT LUMO energy (eV):',eGT(LUMO)*HaToeV write(*,'(2X,A30,F15.6)') 'evGT HOMO-LUMO gap (eV):',Gap*HaToeV write(*,*)'-------------------------------------------------------------------------------' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (singlet) =',EcRPA(1),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy (triplet) =',EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 correlation energy =',EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' Tr@ppRPA@G0T0 total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2),' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 correlation energy =',EcGM,' au' + write(*,'(2X,A50,F20.10,A3)') ' GM@G0T0 total energy =',ENuc + ERHF + EcGM,' au' + write(*,*)'-------------------------------------------------------------------------------' write(*,*) end subroutine print_evGT diff --git a/src/GT/print_qsGT.f90 b/src/GT/print_qsGT.f90 index 5883b62..e383c6e 100644 --- a/src/GT/print_qsGT.f90 +++ b/src/GT/print_qsGT.f90 @@ -48,9 +48,9 @@ subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex write(*,*)'-------------------------------------------------------------------------------' if(nSCF < 10) then - write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I1,A1,I1,A12)')' Self-consistent qsG',nSCF,'T',nSCF,' calculation' else - write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'W',nSCF,' calculation' + write(*,'(1X,A21,I2,A1,I2,A12)')' Self-consistent qsG',nSCF,'T',nSCF,' calculation' endif write(*,*)'-------------------------------------------------------------------------------' write(*,'(1X,A1,1X,A3,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X,A15,1X,A1,1X)') & @@ -70,10 +70,10 @@ subroutine print_qsGT(nBas,nO,nSCF,Conv,thresh,eHF,eGT,c,SigC,Z,ENuc,ET,EV,EJ,Ex write(*,'(2X,A30,F15.6,A3)') 'qsGT LUMO energy:',eGT(LUMO)*HaToeV,' eV' write(*,'(2X,A30,F15.6,A3)') 'qsGT HOMO-LUMO gap :',Gap*HaToeV,' eV' write(*,*)'-------------------------------------------' - write(*,'(2X,A30,F15.6,A3)') ' qsGT total energy:',ENuc + EqsGT,' au' - write(*,'(2X,A30,F15.6,A3)') ' qsGT exchange energy:',Ex,' au' - write(*,'(2X,A30,F15.6,A3)') ' GM@qsGT correlation energy:',EcGM,' au' - write(*,'(2X,A30,F15.6,A3)') 'RPA@qsGT correlation energy:',sum(EcRPA(:)),' au' + write(*,'(2X,A30,F15.6,A3)') ' qsGT total energy:',ENuc + EqsGT,' au' + write(*,'(2X,A30,F15.6,A3)') ' qsGT exchange energy:',Ex,' au' + write(*,'(2X,A30,F15.6,A3)') ' GM@qsGT correlation energy:',EcGM,' au' + write(*,'(2X,A30,F15.6,A3)') 'ppRPA@qsGT correlation energy:',sum(EcRPA(:)),' au' write(*,*)'-------------------------------------------' write(*,*) diff --git a/src/GT/qsGT.f90 b/src/GT/qsGT.f90 index 6a5aa5a..1e4c349 100644 --- a/src/GT/qsGT.f90 +++ b/src/GT/qsGT.f90 @@ -70,7 +70,6 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T integer :: nOOs,nOOt integer :: nVVs,nVVt - logical :: print_W = .false. double precision,allocatable :: error_diis(:,:) double precision,allocatable :: F_diis(:,:) double precision,allocatable :: c(:,:) @@ -160,7 +159,7 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - rcond = 1d0 + rcond = 0d0 !------------------------------------------------------------------------ ! Main loop @@ -198,8 +197,12 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) + EcRPA(1) = EcRPA(1) - EcRPA(2) + EcRPA(2) = 3d0*EcRPA(2) + ! Compute correlation part of the self-energy + EcGM = 0d0 SigT(:,:) = 0d0 Z(:) = 0d0 @@ -208,22 +211,47 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOs,nVVs,ERI_MO, & X1s,Y1s,rho1s,X2s,Y2s,rho2s) - call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Omega1s,rho1s,Omega2s,rho2s,SigT) + if(regularize) then - call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & - Omega1s,rho1s,Omega2s,rho2s,Z) + call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + + call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,Z) + + else + + call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,EcGM,SigT) + + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOs,nVVs,eGT, & + Omega1s,rho1s,Omega2s,rho2s,Z) + + end if iblock = 4 call excitation_density_Tmatrix(iblock,nBas,nC,nO,nV,nR,nOOt,nVVt,ERI_MO, & X1t,Y1t,rho1t,X2t,Y2t,rho2t) - call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Omega1t,rho1t,Omega2t,rho2t,SigT) + if(regularize) then - call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & - Omega1t,rho1t,Omega2t,rho2t,Z) + call self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + + call renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,Z) + + else + + + call regularized_self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,EcGM,SigT) + + call regularized_renormalization_factor_Tmatrix(eta,nBas,nC,nO,nV,nR,nOOt,nVVt,eGT, & + Omega1t,rho1t,Omega2t,rho2t,Z) + + end if Z(:) = 1d0/(1d0 - Z(:)) @@ -302,28 +330,6 @@ subroutine qsGT(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,BSE,TDA_T,T ! End main loop !------------------------------------------------------------------------ -! Compute the ppRPA correlation energy - - ispin = 1 - iblock = 3 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOs,nVVs,1d0,eGT,ERI_MO, & - Omega1s,X1s,Y1s,Omega2s,X2s,Y2s,EcRPA(ispin)) - ispin = 2 - iblock = 4 - call linear_response_pp(iblock,TDA_T,nBas,nC,nO,nV,nR,nOOt,nVVt,1d0,eGT,ERI_MO, & - Omega1t,X1t,Y1t,Omega2t,X2t,Y2t,EcRPA(ispin)) - EcRPA(1) = EcRPA(1) - EcRPA(2) - EcRPA(2) = 3d0*EcRPA(2) - - write(*,*) - write(*,*)'-------------------------------------------------------------------------------' - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (singlet) =',EcRPA(1) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy (triplet) =',EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT correlation energy =',EcRPA(1) + EcRPA(2) - write(*,'(2X,A50,F20.10)') 'Tr@ppRPA@qsGT total energy =',ENuc + ERHF + EcRPA(1) + EcRPA(2) - write(*,*)'-------------------------------------------------------------------------------' - write(*,*) - ! Did it actually converge? if(nSCF == maxSCF+1) then diff --git a/src/GT/self_energy_Tmatrix.f90 b/src/GT/self_energy_Tmatrix.f90 index 3a6c515..13c211f 100644 --- a/src/GT/self_energy_Tmatrix.f90 +++ b/src/GT/self_energy_Tmatrix.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT) +subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) ! Compute the correlation part of the T-matrix self-energy @@ -23,12 +23,13 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2 ! Local variables - integer :: i,a,p,q,cd,kl + integer :: i,j,a,b,p,q,cd,kl double precision :: eps ! Output variables - double precision,intent(inout) :: SigT(nBas,nBas) + double precision,intent(inout):: EcGM + double precision,intent(inout):: SigT(nBas,nBas) !---------------------------------------------- ! Occupied part of the T-matrix self-energy @@ -46,7 +47,7 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2 enddo !---------------------------------------------- - ! Virtual part of the T-matrix self-energy +! Virtual part of the T-matrix self-energy !---------------------------------------------- do p=nC+1,nBas-nR @@ -60,4 +61,26 @@ subroutine self_energy_Tmatrix(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2 enddo enddo +!---------------------------------------------- +! Galitskii-Migdal correlation energy +!---------------------------------------------- + + do i=nC+1,nO + do j=nC+1,nO + do cd=1,nVV + eps = e(i) + e(j) - Omega1(cd) + EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + do kl=1,nOO + eps = e(a) + e(b) - Omega2(kl) + EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + end subroutine self_energy_Tmatrix diff --git a/src/GT/self_energy_Tmatrix_diag.f90 b/src/GT/self_energy_Tmatrix_diag.f90 index 76884d0..c52eb2d 100644 --- a/src/GT/self_energy_Tmatrix_diag.f90 +++ b/src/GT/self_energy_Tmatrix_diag.f90 @@ -1,4 +1,4 @@ -subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,SigT) +subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,Omega2,rho2,EcGM,SigT) ! Compute diagonal of the correlation part of the T-matrix self-energy @@ -23,11 +23,12 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O ! Local variables - integer :: i,a,p,cd,kl + integer :: i,j,a,b,p,cd,kl double precision :: eps ! Output variables + double precision,intent(inout) :: EcGM double precision,intent(inout) :: SigT(nBas) !---------------------------------------------- @@ -56,4 +57,26 @@ subroutine self_energy_Tmatrix_diag(eta,nBas,nC,nO,nV,nR,nOO,nVV,e,Omega1,rho1,O enddo enddo +!---------------------------------------------- +! Galitskii-Migdal correlation energy +!---------------------------------------------- + + do i=nC+1,nO + do j=nC+1,nO + do cd=1,nVV + eps = e(i) + e(j) - Omega1(cd) + EcGM = EcGM + rho1(i,j,cd)*rho1(i,j,cd)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + + do a=nO+1,nBas-nR + do b=nO+1,nBas-nR + do kl=1,nOO + eps = e(a) + e(b) - Omega2(kl) + EcGM = EcGM - rho2(a,b,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) + enddo + enddo + enddo + end subroutine self_energy_Tmatrix_diag diff --git a/src/GT/static_Tmatrix_A.f90 b/src/GT/static_Tmatrix_A.f90 index 9e3505e..ee8158d 100644 --- a/src/GT/static_Tmatrix_A.f90 +++ b/src/GT/static_Tmatrix_A.f90 @@ -1,4 +1,4 @@ -subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TA) +subroutine static_Tmatrix_A(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,TA) ! Compute the OOVV block of the static T-matrix for the resonant block @@ -7,7 +7,6 @@ subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome ! Input variables - integer,intent(in) :: ispin double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -18,7 +17,6 @@ subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: Omega1(nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Omega2(nOO) @@ -34,6 +32,8 @@ subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome double precision,intent(out) :: TA(nS,nS) + TA(:,:) = 0d0 + ia = 0 do i=nC+1,nO do a=nO+1,nBas-nR @@ -46,13 +46,13 @@ subroutine static_Tmatrix_A(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome chi = 0d0 do cd=1,nVV - eps = - Omega1(cd) + eps = + Omega1(cd) ! chi = chi + lambda*rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) chi = chi + rho1(i,b,cd)*rho1(a,j,cd)*eps/(eps**2 + eta**2) enddo do kl=1,nOO - eps = + Omega2(kl) + eps = - Omega2(kl) ! chi = chi - lambda*rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) chi = chi + rho2(i,b,kl)*rho2(a,j,kl)*eps/(eps**2 + eta**2) enddo diff --git a/src/GT/static_Tmatrix_B.f90 b/src/GT/static_Tmatrix_B.f90 index 7cf9bb0..2c6f7aa 100644 --- a/src/GT/static_Tmatrix_B.f90 +++ b/src/GT/static_Tmatrix_B.f90 @@ -1,4 +1,4 @@ -subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Omega1,rho1,Omega2,rho2,TB) +subroutine static_Tmatrix_B(eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,Omega1,rho1,Omega2,rho2,TB) ! Compute the OVVO block of the static T-matrix for the coupling block @@ -7,7 +7,6 @@ subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome ! Input variables - integer,intent(in) :: ispin double precision,intent(in) :: eta integer,intent(in) :: nBas integer,intent(in) :: nC @@ -18,7 +17,6 @@ subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome integer,intent(in) :: nOO integer,intent(in) :: nVV double precision,intent(in) :: lambda - double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) double precision,intent(in) :: Omega1(nVV) double precision,intent(in) :: rho1(nBas,nBas,nVV) double precision,intent(in) :: Omega2(nOO) @@ -34,6 +32,8 @@ subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome double precision,intent(out) :: TB(nS,nS) + TB(:,:) = 0d0 + ia = 0 do i=nC+1,nO do a=nO+1,nBas-nR @@ -46,13 +46,13 @@ subroutine static_Tmatrix_B(ispin,eta,nBas,nC,nO,nV,nR,nS,nOO,nVV,lambda,ERI,Ome chi = 0d0 do cd=1,nVV - eps = - Omega1(cd) + eps = + Omega1(cd) ! chi = chi + lambda*rho1(i,b,cd)*rho1(a,j,cd)*Omega1(cd)/Omega1(cd)**2 + eta**2 chi = chi + rho1(i,j,cd)*rho1(a,b,cd)*eps/(eps**2 + eta**2) enddo do kl=1,nOO - eps = + Omega2(kl) + eps = - Omega2(kl) ! chi = chi + lambda*rho2(i,b,kl)*rho2(a,j,kl)*Omega2(kl)/Omega2(kl)**2 + eta**2 chi = chi + rho2(i,j,kl)*rho2(a,b,kl)*eps/(eps**2 + eta**2) enddo diff --git a/src/GW/Bethe_Salpeter_dynamic_perturbation.f90 b/src/GW/Bethe_Salpeter_dynamic_perturbation.f90 index 0af6209..a8fddcf 100644 --- a/src/GW/Bethe_Salpeter_dynamic_perturbation.f90 +++ b/src/GW/Bethe_Salpeter_dynamic_perturbation.f90 @@ -29,7 +29,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,e integer :: ia - integer,parameter :: maxS = 10 + integer :: maxS = 10 double precision :: gapGW double precision,allocatable :: OmDyn(:) @@ -51,7 +51,8 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,e ! Memory allocation - allocate(OmDyn(nS),ZDyn(nS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) + maxS = min(nS,maxS) + allocate(OmDyn(maxS),ZDyn(maxS),X(nS),Y(nS),Ap_dyn(nS,nS),ZAp_dyn(nS,nS)) if(.not.dTDA) allocate(Am_dyn(nS,nS),ZAm_dyn(nS,nS),Bp_dyn(nS,nS),ZBp_dyn(nS,nS),Bm_dyn(nS,nS),ZBm_dyn(nS,nS)) @@ -71,7 +72,7 @@ subroutine Bethe_Salpeter_dynamic_perturbation(dTDA,eta,nBas,nC,nO,nV,nR,nS,eW,e write(*,'(2X,A5,1X,A20,1X,A20,1X,A20,1X,A20)') '#','Static (eV)','Dynamic (eV)','Correction (eV)','Renorm. (eV)' write(*,*) '---------------------------------------------------------------------------------------------------' - do ia=1,min(nS,maxS) + do ia=1,maxS X(:) = 0.5d0*(XpY(ia,:) + XmY(ia,:)) Y(:) = 0.5d0*(XpY(ia,:) - XmY(ia,:)) diff --git a/src/GW/UG0W0.f90 b/src/GW/UG0W0.f90 index 9c7364d..975e754 100644 --- a/src/GW/UG0W0.f90 +++ b/src/GW/UG0W0.f90 @@ -132,21 +132,25 @@ subroutine UG0W0(doACFDT,exchange_kernel,doXBS,COHSEX,BSE,TDA_W,TDA,dBSE,dTDA,ev call unrestricted_excitation_density(nBas,nC,nO,nR,nS_aa,nS_bb,nS_sc,ERI_aaaa,ERI_aabb,ERI_bbbb,XpY_RPA,rho_RPA) -!---------------------! -! Compute self-energy ! -!---------------------! +!------------------------------------------------! +! Compute self-energy and renormalization factor ! +!------------------------------------------------! do is=1,nspin call self_energy_exchange_diag(nBas,cHF(:,:,is),PHF(:,:,is),ERI,SigX(:,is)) end do - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) + if(regularize) then -!--------------------------------! -! Compute renormalization factor ! -!--------------------------------! + call unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + else + + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + + end if !-----------------------------------! ! Solve the quasi-particle equation ! diff --git a/src/GW/evGW.f90 b/src/GW/evGW.f90 index 4790cf4..f510975 100644 --- a/src/GW/evGW.f90 +++ b/src/GW/evGW.f90 @@ -139,6 +139,7 @@ subroutine evGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, eGW(:) = eG0W0(:) eOld(:) = eGW(:) Z(:) = 1d0 + rcond = 0d0 !------------------------------------------------------------------------ ! Main loop diff --git a/src/GW/evUGW.f90 b/src/GW/evUGW.f90 index 3e5faec..ceac1dd 100644 --- a/src/GW/evUGW.f90 +++ b/src/GW/evUGW.f90 @@ -153,6 +153,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE eGW(:,:) = eG0W0(:,:) eOld(:,:) = eGW(:,:) Z(:,:) = 1d0 + rcond(:) = 0d0 !------------------------------------------------------------------------ ! Main loop @@ -166,7 +167,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE call unrestricted_linear_response(ispin,.true.,TDA_W,.false.,eta,nBas,nC,nO,nV,nR,nS_aa,nS_bb,nS_sc,nS_sc,1d0, & eGW,ERI_aaaa,ERI_aabb,ERI_bbbb,OmRPA,rho_RPA,EcRPA,OmRPA,XpY_RPA,XmY_RPA) - endif + end if !----------------------! ! Excitation densities ! @@ -180,15 +181,33 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE if(G0W) then - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + if(regularize) then + + call unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + + else + + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + + end if else - call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + if(regularize) then - endif + call unrestricted_regularized_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + + else + + call unrestricted_self_energy_correlation_diag(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + + end if + + end if !-----------------------------------! ! Solve the quasi-particle equation ! @@ -222,7 +241,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE if(minval(rcond(:)) < 1d-15) n_diis = 0 - endif + end if ! Save quasiparticles energy for next cycle @@ -253,7 +272,7 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE stop - endif + end if ! Deallocate memory @@ -316,6 +335,6 @@ subroutine evUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE end if - endif + end if end subroutine evUGW diff --git a/src/GW/print_qsUGW.f90 b/src/GW/print_qsUGW.f90 index 983e148..9afcbde 100644 --- a/src/GW/print_qsUGW.f90 +++ b/src/GW/print_qsUGW.f90 @@ -1,4 +1,4 @@ -subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, & +subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,Ov, & ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigC,Z,dipole) ! Print one-electron energies and other stuff for qsUGW @@ -24,12 +24,7 @@ subroutine print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,cGW,PGW,Ov,T,V,J,K, & double precision,intent(in) :: eHF(nBas,nspin) double precision,intent(in) :: eGW(nBas,nspin) double precision,intent(in) :: cGW(nBas,nBas,nspin) - double precision,intent(in) :: PGW(nBas,nBas,nspin) double precision,intent(in) :: Ov(nBas,nBas) - double precision,intent(in) :: T(nBas,nBas) - double precision,intent(in) :: V(nBas,nBas) - double precision,intent(in) :: J(nBas,nBas,nspin) - double precision,intent(in) :: K(nBas,nBas,nspin) double precision,intent(in) :: SigC(nBas,nBas,nspin) double precision,intent(in) :: Z(nBas,nspin) double precision,intent(in) :: dipole(ncart) diff --git a/src/GW/qsGW.f90 b/src/GW/qsGW.f90 index 9a418aa..806ae6d 100644 --- a/src/GW/qsGW.f90 +++ b/src/GW/qsGW.f90 @@ -153,7 +153,7 @@ subroutine qsGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE, c(:,:) = cHF(:,:) F_diis(:,:) = 0d0 error_diis(:,:) = 0d0 - rcond = 1d0 + rcond = 0d0 !------------------------------------------------------------------------ ! Main loop diff --git a/src/GW/qsUGW.f90 b/src/GW/qsUGW.f90 index 22745fe..acfab9f 100644 --- a/src/GW/qsUGW.f90 +++ b/src/GW/qsUGW.f90 @@ -166,7 +166,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE c(:,:,:) = cHF(:,:,:) F_diis(:,:,:) = 0d0 error_diis(:,:,:) = 0d0 - rcond = 1d0 + rcond(:) = 0d0 !------------------------------------------------------------------------ ! Main loop @@ -227,13 +227,31 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE if(G0W) then - call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + if(regularize) then + + call unrestricted_regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + + else + + call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eHF,OmRPA,rho_RPA,Z) + + end if else - call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) - call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + if(regularize) then + + call unrestricted_regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_regularized_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + + else + + call unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,SigC,EcGM) + call unrestricted_renormalization_factor(eta,nBas,nC,nO,nV,nR,nS_sc,eGW,OmRPA,rho_RPA,Z) + + end if endif @@ -348,7 +366,7 @@ subroutine qsUGW(maxSCF,thresh,max_diis,doACFDT,exchange_kernel,doXBS,COHSEX,BSE !------------------------------------------------------------------------ call dipole_moment(nBas,P(:,:,1)+P(:,:,2),nNuc,ZNuc,rNuc,dipole_int_AO,dipole) - call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,P,S,T,V,J,K,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigCp,Z,dipole) + call print_qsUGW(nBas,nO,nSCF,Conv,thresh,eHF,eGW,c,S,ENuc,ET,EV,EJ,Ex,EcGM,EcRPA,EqsGW,SigCp,Z,dipole) enddo !------------------------------------------------------------------------ diff --git a/src/GW/self_energy_correlation.f90 b/src/GW/self_energy_correlation.f90 index 4d62b43..ebe0b67 100644 --- a/src/GW/self_energy_correlation.f90 +++ b/src/GW/self_energy_correlation.f90 @@ -28,8 +28,8 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec ! Output variables - double precision,intent(out) :: SigC(nBas,nBas) double precision,intent(out) :: EcGM + double precision,intent(out) :: SigC(nBas,nBas) ! Initialize @@ -102,7 +102,7 @@ subroutine self_energy_correlation(COHSEX,eta,nBas,nC,nO,nV,nR,nS,e,Omega,rho,Ec end do end do - ! GM correlation energy + ! Galitskii-Migdal correlation energy EcGM = 0d0 do i=nC+1,nO diff --git a/src/GW/ufG0W0.f90 b/src/GW/ufG0W0.f90 index 51ac629..ad9f079 100644 --- a/src/GW/ufG0W0.f90 +++ b/src/GW/ufG0W0.f90 @@ -33,6 +33,10 @@ subroutine ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) double precision,allocatable :: eGW(:) double precision,allocatable :: Z(:) + logical :: verbose = .true. + double precision,parameter :: cutoff1 = 0.01d0 + double precision,parameter :: cutoff2 = 0.01d0 + ! Output variables ! Hello world @@ -183,9 +187,8 @@ subroutine ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) ! Compute weights ! !-----------------! - Z(:) = 0d0 do s=1,nH - Z(s) = Z(s) + cGW(1,s)**2 + Z(s) = cGW(1,s)**2 end do !--------------! @@ -207,6 +210,64 @@ subroutine ufG0W0(nBas,nC,nO,nV,nR,nS,ENuc,ERHF,ERI,eHF) write(*,*)'-------------------------------------------' write(*,*) + if(verbose) then + + do s=1,nH + + if(Z(s) > cutoff1) then + + write(*,*)'*************************************************************' + write(*,'(1X,A20,I3,A6,I3)')'Vector for orbital ',p,' and #',s + write(*,'(1X,A7,F10.6,A13,F10.6,1X)')' e_QP = ',eGW(s)*HaToeV,' eV and Z = ',Z(s) + write(*,*)'*************************************************************' + write(*,'(1X,A20,1X,A20,1X,A15,1X)') & + ' Configuration ',' Coefficient ',' Weight ' + write(*,*)'*************************************************************' + + if(p <= nO) & + write(*,'(1X,A7,I3,A16,1X,F15.6,1X,F15.6)') & + ' (',p,') ',cGW(1,s),cGW(1,s)**2 + if(p > nO) & + write(*,'(1X,A16,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',p,') ',cGW(1,s),cGW(1,s)**2 + + klc = 0 + do k=nC+1,nO + do l=nC+1,nO + do c=nO+1,nBas-nR + + klc = klc + 1 +! if(abs(cGW(1+klc,s)) > cutoff2) & + write(*,'(1X,A3,I3,A1,I3,A6,I3,A7,1X,F15.6,1X,F15.6)') & + ' (',k,',',l,') -> (',c,') ',cGW(1+klc,s),cGW(1+klc,s)**2 + + end do + end do + end do + + kcd = 0 + do k=nC+1,nO + do c=nO+1,nBas-nR + do d=nO+1,nBas-nR + + kcd = kcd + 1 +! if(abs(cGW(1+n2h1p+kcd,s)) > cutoff2) & + write(*,'(1X,A7,I3,A6,I3,A1,I3,A3,1X,F15.6,1X,F15.6)') & + ' (',k,') -> (',c,',',d,') ',cGW(1+n2h1p+kcd,s),cGW(1+n2h1p+kcd,s)**2 + + end do + end do + end do + + write(*,*)'*************************************************************' + write(*,*) + + end if + + end do + + end if + end do end subroutine ufG0W0 diff --git a/src/GW/unrestricted_regularized_self_energy_correlation.f90 b/src/GW/unrestricted_regularized_self_energy_correlation.f90 index 9553e84..8c28dbc 100644 --- a/src/GW/unrestricted_regularized_self_energy_correlation.f90 +++ b/src/GW/unrestricted_regularized_self_energy_correlation.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM) +subroutine unrestricted_regularized_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega,rho,SigC,EcGM) ! Compute diagonal of the correlation part of the self-energy @@ -130,4 +130,4 @@ subroutine unrestricted_self_energy_correlation(eta,nBas,nC,nO,nV,nR,nSt,e,Omega end do end do -end subroutine unrestricted_self_energy_correlation +end subroutine unrestricted_regularized_self_energy_correlation diff --git a/src/HF/RHF.f90 b/src/HF/RHF.f90 index a4bf5fe..a66f2aa 100644 --- a/src/HF/RHF.f90 +++ b/src/HF/RHF.f90 @@ -94,6 +94,7 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T error_diis(:,:) = 0d0 Conv = 1d0 nSCF = 0 + rcond = 0d0 !------------------------------------------------------------------------ ! Main SCF loop @@ -204,6 +205,6 @@ subroutine RHF(maxSCF,thresh,max_diis,guess_type,nNuc,ZNuc,rNuc,ENuc,nBas,nO,S,T ! Compute Vx for post-HF calculations - call exchange_potential(nBas,c,K,Vx) + call mo_fock_exchange_potential(nBas,c,K,Vx) end subroutine RHF diff --git a/src/HF/UHF.f90 b/src/HF/UHF.f90 index 9443924..afa1fad 100644 --- a/src/HF/UHF.f90 +++ b/src/HF/UHF.f90 @@ -253,7 +253,7 @@ subroutine UHF(maxSCF,thresh,max_diis,guess_type,mix,nNuc,ZNuc,rNuc,ENuc,nBas,nO ! Compute Vx for post-HF calculations do ispin=1,nspin - call exchange_potential(nBas,c(:,:,ispin),K(:,:,ispin),Vx(:,ispin)) + call mo_fock_exchange_potential(nBas,c(:,:,ispin),K(:,:,ispin),Vx(:,ispin)) end do end subroutine UHF diff --git a/src/HF/exchange_potential.f90 b/src/HF/mo_fock_exchange_potential.f90 similarity index 85% rename from src/HF/exchange_potential.f90 rename to src/HF/mo_fock_exchange_potential.f90 index 9f72ce7..56203c3 100644 --- a/src/HF/exchange_potential.f90 +++ b/src/HF/mo_fock_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine exchange_potential(nBas,c,Fx,Vx) +subroutine mo_fock_exchange_potential(nBas,c,Fx,Vx) ! Compute the exchange potential in the MO basis @@ -31,4 +31,4 @@ subroutine exchange_potential(nBas,c,Fx,Vx) end do end do -end subroutine exchange_potential +end subroutine mo_fock_exchange_potential diff --git a/src/LR/linear_response_Tmatrix.f90 b/src/LR/linear_response_Tmatrix.f90 index c0a4acf..89b9207 100644 --- a/src/LR/linear_response_Tmatrix.f90 +++ b/src/LR/linear_response_Tmatrix.f90 @@ -42,14 +42,13 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda call linear_response_A_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,e,ERI,A) - if(ispin == 1) A(:,:) = A(:,:) + A_BSE(:,:) - if(ispin == 2) A(:,:) = A(:,:) - A_BSE(:,:) - ! print*,'A' ! call matout(nS,nS,A) ! print*,'TA' ! call matout(nS,nS,A_BSE) + A(:,:) = A(:,:) - A_BSE(:,:) + ! Tamm-Dancoff approximation if(TDA) then @@ -64,13 +63,12 @@ subroutine linear_response_Tmatrix(ispin,dRPA,TDA,eta,nBas,nC,nO,nV,nR,nS,lambda call linear_response_B_matrix(ispin,dRPA,nBas,nC,nO,nV,nR,nS,lambda,ERI,B) - if(ispin == 1) B(:,:) = B(:,:) + B_BSE(:,:) - if(ispin == 2) B(:,:) = B(:,:) - B_BSE(:,:) +! print*,'B' +! call matout(nS,nS,B) +! print*,'TB' +! call matout(nS,nS,B_BSE) -! print*,'B' -! call matout(nS,nS,B) -! print*,'TB' -! call matout(nS,nS,B_BSE) + B(:,:) = B(:,:) - B_BSE(:,:) ! Build A + B and A - B matrices diff --git a/src/QuAcK/QuAcK.f90 b/src/QuAcK/QuAcK.f90 index 48f2189..2959cec 100644 --- a/src/QuAcK/QuAcK.f90 +++ b/src/QuAcK/QuAcK.f90 @@ -451,8 +451,10 @@ program QuAcK ket1 = 1 ket2 = 1 call AOtoMO_integral_transform(bra1,bra2,ket1,ket2,nBas,cHF,ERI_AO,ERI_MO) + F_MO(:,:) = F_AO(:,:) call AOtoMO_transform(nBas,cHF,F_MO) + end if end if diff --git a/src/QuAcK/read_methods.f90 b/src/QuAcK/read_methods.f90 index 8ddd885..9b909e9 100644 --- a/src/QuAcK/read_methods.f90 +++ b/src/QuAcK/read_methods.f90 @@ -30,7 +30,7 @@ subroutine read_methods(doRHF,doUHF,doKS,doMOM, & ! Local variables - character(len=1) :: answer1,answer2,answer3,answer4,answer5,answer6,answer7 + character(len=1) :: answer1,answer2,answer3,answer4,answer5 ! Open file with method specification diff --git a/src/QuAcK/read_options.f90 b/src/QuAcK/read_options.f90 index 382ec95..383aaee 100644 --- a/src/QuAcK/read_options.f90 +++ b/src/QuAcK/read_options.f90 @@ -195,15 +195,15 @@ subroutine read_options(maxSCF_HF,thresh_HF,DIIS_HF,n_diis_HF,guess_type,ortho_t if(answer8 == 'T') regGW = .true. if(.not.DIIS_GW) n_diis_GW = 1 -! Read GF options +! Read GT options - maxSCF_GF = 64 - thresh_GF = 1d-5 - DIIS_GF = .false. - n_diis_GF = 5 - linGF = .false. - eta_GF = 0d0 - regGF = .false. + maxSCF_GT = 64 + thresh_GT = 1d-5 + DIIS_GT = .false. + n_diis_GT = 5 + linGT = .false. + eta_GT = 0d0 + regGT = .false. TDA_T = .false. read(1,*) diff --git a/src/eDFT/UB88_gga_exchange_energy.f90 b/src/eDFT/B88_gga_exchange_energy.f90 similarity index 89% rename from src/eDFT/UB88_gga_exchange_energy.f90 rename to src/eDFT/B88_gga_exchange_energy.f90 index 2ecf64d..ed7e221 100644 --- a/src/eDFT/UB88_gga_exchange_energy.f90 +++ b/src/eDFT/B88_gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) +subroutine B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Compute Becke's 88 GGA exchange energy @@ -45,4 +45,4 @@ subroutine UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) end do -end subroutine UB88_gga_exchange_energy +end subroutine B88_gga_exchange_energy diff --git a/src/eDFT/UB88_gga_exchange_potential.f90 b/src/eDFT/B88_gga_exchange_potential.f90 similarity index 94% rename from src/eDFT/UB88_gga_exchange_potential.f90 rename to src/eDFT/B88_gga_exchange_potential.f90 index fcec63c..02d2ca7 100644 --- a/src/eDFT/UB88_gga_exchange_potential.f90 +++ b/src/eDFT/B88_gga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Compute Becke's GGA exchange potential @@ -70,4 +70,4 @@ subroutine UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) end do end do -end subroutine UB88_gga_exchange_potential +end subroutine B88_gga_exchange_potential diff --git a/src/eDFT/UC16_lda_correlation_energy.f90 b/src/eDFT/C16_lda_correlation_energy.f90 similarity index 94% rename from src/eDFT/UC16_lda_correlation_energy.f90 rename to src/eDFT/C16_lda_correlation_energy.f90 index d069e4f..3b9df41 100644 --- a/src/eDFT/UC16_lda_correlation_energy.f90 +++ b/src/eDFT/C16_lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine UC16_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine C16_lda_correlation_energy(nGrid,weight,rho,Ec) ! Compute unrestricted Chachiyo's LDA correlation energy @@ -90,4 +90,4 @@ subroutine UC16_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = Ec(2) - Ec(1) - Ec(3) -end subroutine UC16_lda_correlation_energy +end subroutine C16_lda_correlation_energy diff --git a/src/eDFT/UC16_lda_correlation_potential.f90 b/src/eDFT/C16_lda_correlation_potential.f90 similarity index 96% rename from src/eDFT/UC16_lda_correlation_potential.f90 rename to src/eDFT/C16_lda_correlation_potential.f90 index 43fa599..aa58e0b 100644 --- a/src/eDFT/UC16_lda_correlation_potential.f90 +++ b/src/eDFT/C16_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine UC16_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine C16_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) ! Compute unrestricted LDA correlation potential @@ -128,4 +128,4 @@ include 'parameters.h' enddo enddo -end subroutine UC16_lda_correlation_potential +end subroutine C16_lda_correlation_potential diff --git a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/CC_lda_exchange_derivative_discontinuity.f90 similarity index 84% rename from src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 rename to src/eDFT/CC_lda_exchange_derivative_discontinuity.f90 index 78d6822..e7998d2 100644 --- a/src/eDFT/UCC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/CC_lda_exchange_derivative_discontinuity.f90 @@ -1,5 +1,6 @@ -subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,doNcentered,& - kappa,ExDD) +subroutine CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,Cx_choice,& + doNcentered,kappa,ExDD) + ! Compute the unrestricted version of the curvature-corrected exchange ensemble derivative @@ -152,26 +153,18 @@ subroutine UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,wei ExDD(:) = 0d0 - if (doNcentered) then + do iEns=1,nEns + do jEns=2,nEns - do iEns=1,nEns - do jEns=2,nEns + if(doNcentered) then - ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - kappa(iEns)* wEns(jEns))*dExdw(jEns) + ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - kappa(iEns)*wEns(jEns))*dExdw(jEns) + else + + ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) + end if - end do end do + end do - else - - do iEns=1,nEns - do jEns=2,nEns - - ExDD(iEns) = ExDD(iEns) + (Kronecker_delta(iEns,jEns) - wEns(jEns))*dExdw(jEns) - - end do - end do - - endif - -end subroutine UCC_lda_exchange_derivative_discontinuity +end subroutine CC_lda_exchange_derivative_discontinuity diff --git a/src/eDFT/UCC_lda_exchange_energy.f90 b/src/eDFT/CC_lda_exchange_energy.f90 similarity index 93% rename from src/eDFT/UCC_lda_exchange_energy.f90 rename to src/eDFT/CC_lda_exchange_energy.f90 index 0a5b615..d4a07d0 100644 --- a/src/eDFT/UCC_lda_exchange_energy.f90 +++ b/src/eDFT/CC_lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine UCC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) +subroutine CC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) ! Compute the unrestricted version of the curvature-corrected exchange functional @@ -107,4 +107,4 @@ subroutine UCC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice, enddo -end subroutine UCC_lda_exchange_energy +end subroutine CC_lda_exchange_energy diff --git a/src/eDFT/UCC_lda_exchange_individual_energy.f90 b/src/eDFT/CC_lda_exchange_individual_energy.f90 similarity index 93% rename from src/eDFT/UCC_lda_exchange_individual_energy.f90 rename to src/eDFT/CC_lda_exchange_individual_energy.f90 index 06835df..29b0778 100644 --- a/src/eDFT/UCC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/CC_lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,LZx,Ex) +subroutine CC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho,Cx_choice,doNcentered,LZx,Ex) ! Compute the unrestricted version of the curvature-corrected exchange functional @@ -128,4 +128,4 @@ subroutine UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho enddo -end subroutine UCC_lda_exchange_individual_energy +end subroutine CC_lda_exchange_individual_energy diff --git a/src/eDFT/UCC_lda_exchange_potential.f90 b/src/eDFT/CC_lda_exchange_potential.f90 similarity index 93% rename from src/eDFT/UCC_lda_exchange_potential.f90 rename to src/eDFT/CC_lda_exchange_potential.f90 index ea244a2..7cc753e 100644 --- a/src/eDFT/UCC_lda_exchange_potential.f90 +++ b/src/eDFT/CC_lda_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Cx_choice,doNcentered,Fx) +subroutine CC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Cx_choice,doNcentered,Fx) ! Compute the unrestricted version of the curvature-corrected exchange potential @@ -116,4 +116,4 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho enddo enddo -end subroutine UCC_lda_exchange_potential +end subroutine CC_lda_exchange_potential diff --git a/src/eDFT/UG96_gga_exchange_energy.f90 b/src/eDFT/G96_gga_exchange_energy.f90 similarity index 89% rename from src/eDFT/UG96_gga_exchange_energy.f90 rename to src/eDFT/G96_gga_exchange_energy.f90 index f09d2a7..93c3ece 100644 --- a/src/eDFT/UG96_gga_exchange_energy.f90 +++ b/src/eDFT/G96_gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) +subroutine G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Compute Gill's 96 GGA exchange energy @@ -45,4 +45,4 @@ subroutine UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) end do -end subroutine UG96_gga_exchange_energy +end subroutine G96_gga_exchange_energy diff --git a/src/eDFT/UG96_gga_exchange_potential.f90 b/src/eDFT/G96_gga_exchange_potential.f90 similarity index 93% rename from src/eDFT/UG96_gga_exchange_potential.f90 rename to src/eDFT/G96_gga_exchange_potential.f90 index 4b28a3c..029354c 100644 --- a/src/eDFT/UG96_gga_exchange_potential.f90 +++ b/src/eDFT/G96_gga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Compute Gill's GGA exchange poential @@ -61,4 +61,4 @@ subroutine UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) enddo enddo -end subroutine UG96_gga_exchange_potential +end subroutine G96_gga_exchange_potential diff --git a/src/eDFT/ULYP_gga_correlation_energy.f90 b/src/eDFT/LYP_gga_correlation_energy.f90 similarity index 95% rename from src/eDFT/ULYP_gga_correlation_energy.f90 rename to src/eDFT/LYP_gga_correlation_energy.f90 index e0788b8..a86295a 100644 --- a/src/eDFT/ULYP_gga_correlation_energy.f90 +++ b/src/eDFT/LYP_gga_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) +subroutine LYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) ! Compute unrestricted LYP GGA correlation energy @@ -70,4 +70,4 @@ subroutine ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) end do -end subroutine ULYP_gga_correlation_energy +end subroutine LYP_gga_correlation_energy diff --git a/src/eDFT/ULYP_gga_correlation_potential.f90 b/src/eDFT/LYP_gga_correlation_potential.f90 similarity index 98% rename from src/eDFT/ULYP_gga_correlation_potential.f90 rename to src/eDFT/LYP_gga_correlation_potential.f90 index 702aba3..d30507d 100644 --- a/src/eDFT/ULYP_gga_correlation_potential.f90 +++ b/src/eDFT/LYP_gga_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine LYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute LYP correlation potential @@ -153,4 +153,4 @@ subroutine ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) end do end do -end subroutine ULYP_gga_correlation_potential +end subroutine LYP_gga_correlation_potential diff --git a/src/eDFT/UPBE_gga_correlation_energy.f90 b/src/eDFT/PBE_gga_correlation_energy.f90 similarity index 97% rename from src/eDFT/UPBE_gga_correlation_energy.f90 rename to src/eDFT/PBE_gga_correlation_energy.f90 index 461e2bb..c93d812 100644 --- a/src/eDFT/UPBE_gga_correlation_energy.f90 +++ b/src/eDFT/PBE_gga_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine UPBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec) +subroutine PBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec) ! Compute unrestricted PBE GGA correlation energy @@ -169,4 +169,4 @@ subroutine UPBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec) Ec(2) = Ec(2) - Ec(1) - Ec(3) -end subroutine UPBE_gga_correlation_energy +end subroutine PBE_gga_correlation_energy diff --git a/src/eDFT/UPBE_gga_correlation_potential.f90 b/src/eDFT/PBE_gga_correlation_potential.f90 similarity index 93% rename from src/eDFT/UPBE_gga_correlation_potential.f90 rename to src/eDFT/PBE_gga_correlation_potential.f90 index af8794a..40e5f3d 100644 --- a/src/eDFT/UPBE_gga_correlation_potential.f90 +++ b/src/eDFT/PBE_gga_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine PBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute LYP correlation potential @@ -36,7 +36,7 @@ subroutine UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute matrix elements in the AO basis - call UPW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) + call PW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) do mu=1,nBas do nu=1,nBas @@ -85,4 +85,4 @@ subroutine UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) end do end do -end subroutine UPBE_gga_correlation_potential +end subroutine PBE_gga_correlation_potential diff --git a/src/eDFT/UPBE_gga_exchange_energy.f90 b/src/eDFT/PBE_gga_exchange_energy.f90 similarity index 90% rename from src/eDFT/UPBE_gga_exchange_energy.f90 rename to src/eDFT/PBE_gga_exchange_energy.f90 index 0cc35bd..5c76336 100644 --- a/src/eDFT/UPBE_gga_exchange_energy.f90 +++ b/src/eDFT/PBE_gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) +subroutine PBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) ! Compute PBE GGA exchange energy @@ -46,4 +46,4 @@ subroutine UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) end do -end subroutine UPBE_gga_exchange_energy +end subroutine PBE_gga_exchange_energy diff --git a/src/eDFT/UPBE_gga_exchange_potential.f90 b/src/eDFT/PBE_gga_exchange_potential.f90 similarity index 93% rename from src/eDFT/UPBE_gga_exchange_potential.f90 rename to src/eDFT/PBE_gga_exchange_potential.f90 index 866b90a..81a7687 100644 --- a/src/eDFT/UPBE_gga_exchange_potential.f90 +++ b/src/eDFT/PBE_gga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine PBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Compute PBE GGA exchange potential @@ -64,4 +64,4 @@ subroutine UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) end do end do -end subroutine UPBE_gga_exchange_potential +end subroutine PBE_gga_exchange_potential diff --git a/src/eDFT/UPW92_lda_correlation_energy.f90 b/src/eDFT/PW92_lda_correlation_energy.f90 similarity index 96% rename from src/eDFT/UPW92_lda_correlation_energy.f90 rename to src/eDFT/PW92_lda_correlation_energy.f90 index 9a70abc..7b13588 100644 --- a/src/eDFT/UPW92_lda_correlation_energy.f90 +++ b/src/eDFT/PW92_lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine UPW92_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine PW92_lda_correlation_energy(nGrid,weight,rho,Ec) ! Compute unrestricted PW92 LDA correlation energy @@ -117,4 +117,4 @@ subroutine UPW92_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = Ec(2) - Ec(1) - Ec(3) -end subroutine UPW92_lda_correlation_energy +end subroutine PW92_lda_correlation_energy diff --git a/src/eDFT/UPW92_lda_correlation_potential.f90 b/src/eDFT/PW92_lda_correlation_potential.f90 similarity index 98% rename from src/eDFT/UPW92_lda_correlation_potential.f90 rename to src/eDFT/PW92_lda_correlation_potential.f90 index d63d9e4..32dad64 100644 --- a/src/eDFT/UPW92_lda_correlation_potential.f90 +++ b/src/eDFT/PW92_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine UPW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine PW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) ! Compute unrestricted PW92 LDA correlation potential @@ -182,4 +182,4 @@ subroutine UPW92_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) end do end do -end subroutine UPW92_lda_correlation_potential +end subroutine PW92_lda_correlation_potential diff --git a/src/eDFT/US51_lda_exchange_energy.f90 b/src/eDFT/S51_lda_exchange_energy.f90 similarity index 85% rename from src/eDFT/US51_lda_exchange_energy.f90 rename to src/eDFT/S51_lda_exchange_energy.f90 index aad5346..1173c4c 100644 --- a/src/eDFT/US51_lda_exchange_energy.f90 +++ b/src/eDFT/S51_lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex) +subroutine S51_lda_exchange_energy(nGrid,weight,rho,Ex) ! Compute Slater's LDA exchange energy @@ -31,4 +31,4 @@ subroutine US51_lda_exchange_energy(nGrid,weight,rho,Ex) enddo -end subroutine US51_lda_exchange_energy +end subroutine S51_lda_exchange_energy diff --git a/src/eDFT/US51_lda_exchange_individual_energy.f90 b/src/eDFT/S51_lda_exchange_individual_energy.f90 similarity index 90% rename from src/eDFT/US51_lda_exchange_individual_energy.f90 rename to src/eDFT/S51_lda_exchange_individual_energy.f90 index c91d9fc..7bb4313 100644 --- a/src/eDFT/US51_lda_exchange_individual_energy.f90 +++ b/src/eDFT/S51_lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex) +subroutine S51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex) ! Compute the restricted version of Slater's LDA exchange individual energy @@ -58,4 +58,4 @@ subroutine US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex enddo -end subroutine US51_lda_exchange_individual_energy +end subroutine S51_lda_exchange_individual_energy diff --git a/src/eDFT/US51_lda_exchange_potential.f90 b/src/eDFT/S51_lda_exchange_potential.f90 similarity index 88% rename from src/eDFT/US51_lda_exchange_potential.f90 rename to src/eDFT/S51_lda_exchange_potential.f90 index 9c2bc0d..4c0978e 100644 --- a/src/eDFT/US51_lda_exchange_potential.f90 +++ b/src/eDFT/S51_lda_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) +subroutine S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) ! Compute Slater's LDA exchange potential @@ -42,4 +42,4 @@ subroutine US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) enddo enddo -end subroutine US51_lda_exchange_potential +end subroutine S51_lda_exchange_potential diff --git a/src/eDFT/UVWN3_lda_correlation_energy.f90 b/src/eDFT/VWN3_lda_correlation_energy.f90 similarity index 97% rename from src/eDFT/UVWN3_lda_correlation_energy.f90 rename to src/eDFT/VWN3_lda_correlation_energy.f90 index ff9d18e..b5b2aa5 100644 --- a/src/eDFT/UVWN3_lda_correlation_energy.f90 +++ b/src/eDFT/VWN3_lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine VWN3_lda_correlation_energy(nGrid,weight,rho,Ec) ! Compute unrestricted VWN3 LDA correlation energy @@ -134,4 +134,4 @@ subroutine UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = Ec(2) - Ec(1) - Ec(3) -end subroutine UVWN3_lda_correlation_energy +end subroutine VWN3_lda_correlation_energy diff --git a/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 b/src/eDFT/VWN3_lda_correlation_individual_energy.f90 similarity index 97% rename from src/eDFT/UVWN3_lda_correlation_individual_energy.f90 rename to src/eDFT/VWN3_lda_correlation_individual_energy.f90 index d7e4363..5c90110 100644 --- a/src/eDFT/UVWN3_lda_correlation_individual_energy.f90 +++ b/src/eDFT/VWN3_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine UVWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec) +subroutine VWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,doNcentered,LZc,Ec) ! Compute VWN3 LDA correlation potential @@ -178,4 +178,4 @@ subroutine UVWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,do end do -end subroutine UVWN3_lda_correlation_individual_energy +end subroutine VWN3_lda_correlation_individual_energy diff --git a/src/eDFT/UVWN3_lda_correlation_potential.f90 b/src/eDFT/VWN3_lda_correlation_potential.f90 similarity index 98% rename from src/eDFT/UVWN3_lda_correlation_potential.f90 rename to src/eDFT/VWN3_lda_correlation_potential.f90 index ce1c93a..2465cb7 100644 --- a/src/eDFT/UVWN3_lda_correlation_potential.f90 +++ b/src/eDFT/VWN3_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine UVWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine VWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) ! Compute unrestricted VWN3 LDA correlation potential @@ -193,4 +193,4 @@ subroutine UVWN3_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) end do end do -end subroutine UVWN3_lda_correlation_potential +end subroutine VWN3_lda_correlation_potential diff --git a/src/eDFT/UVWN5_lda_correlation_energy.f90 b/src/eDFT/VWN5_lda_correlation_energy.f90 similarity index 97% rename from src/eDFT/UVWN5_lda_correlation_energy.f90 rename to src/eDFT/VWN5_lda_correlation_energy.f90 index 85e50c7..4a6137e 100644 --- a/src/eDFT/UVWN5_lda_correlation_energy.f90 +++ b/src/eDFT/VWN5_lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine VWN5_lda_correlation_energy(nGrid,weight,rho,Ec) ! Compute unrestricted VWN5 LDA correlation energy @@ -134,4 +134,4 @@ subroutine UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = Ec(2) - Ec(1) - Ec(3) -end subroutine UVWN5_lda_correlation_energy +end subroutine VWN5_lda_correlation_energy diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/VWN5_lda_correlation_individual_energy.f90 similarity index 97% rename from src/eDFT/UVWN5_lda_correlation_individual_energy.f90 rename to src/eDFT/VWN5_lda_correlation_individual_energy.f90 index 99b3159..d1330c0 100644 --- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/VWN5_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) +subroutine VWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) ! Compute VWN5 LDA correlation potential @@ -181,4 +181,4 @@ subroutine UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZ end do -end subroutine UVWN5_lda_correlation_individual_energy +end subroutine VWN5_lda_correlation_individual_energy diff --git a/src/eDFT/UVWN5_lda_correlation_potential.f90 b/src/eDFT/VWN5_lda_correlation_potential.f90 similarity index 98% rename from src/eDFT/UVWN5_lda_correlation_potential.f90 rename to src/eDFT/VWN5_lda_correlation_potential.f90 index 0ffb8c0..48d6518 100644 --- a/src/eDFT/UVWN5_lda_correlation_potential.f90 +++ b/src/eDFT/VWN5_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine VWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) ! Compute unrestricted VWN5 LDA correlation potential @@ -190,4 +190,4 @@ subroutine UVWN5_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) end do end do -end subroutine UVWN5_lda_correlation_potential +end subroutine VWN5_lda_correlation_potential diff --git a/src/eDFT/UW38_lda_correlation_energy.f90 b/src/eDFT/W38_lda_correlation_energy.f90 similarity index 89% rename from src/eDFT/UW38_lda_correlation_energy.f90 rename to src/eDFT/W38_lda_correlation_energy.f90 index 391c8c7..4b97f23 100644 --- a/src/eDFT/UW38_lda_correlation_energy.f90 +++ b/src/eDFT/W38_lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine UW38_lda_correlation_energy(nGrid,weight,rho,Ec) +subroutine W38_lda_correlation_energy(nGrid,weight,rho,Ec) ! Compute the unrestricted version of the Wigner's LDA correlation energy @@ -49,4 +49,4 @@ subroutine UW38_lda_correlation_energy(nGrid,weight,rho,Ec) Ec(2) = -4d0*a*Ec(2) -end subroutine UW38_lda_correlation_energy +end subroutine W38_lda_correlation_energy diff --git a/src/eDFT/UW38_lda_correlation_individual_energy.f90 b/src/eDFT/W38_lda_correlation_individual_energy.f90 similarity index 90% rename from src/eDFT/UW38_lda_correlation_individual_energy.f90 rename to src/eDFT/W38_lda_correlation_individual_energy.f90 index 6c859c2..2973185 100644 --- a/src/eDFT/UW38_lda_correlation_individual_energy.f90 +++ b/src/eDFT/W38_lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine UW38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) +subroutine W38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) ! Compute the unrestricted version of the Wigner's LDA individual energy @@ -59,4 +59,4 @@ subroutine UW38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,Ec) Ec(2) = -4d0*a*Ec(2) -end subroutine UW38_lda_correlation_individual_energy +end subroutine W38_lda_correlation_individual_energy diff --git a/src/eDFT/UW38_lda_correlation_potential.f90 b/src/eDFT/W38_lda_correlation_potential.f90 similarity index 93% rename from src/eDFT/UW38_lda_correlation_potential.f90 rename to src/eDFT/W38_lda_correlation_potential.f90 index 8c1c200..5e7a865 100644 --- a/src/eDFT/UW38_lda_correlation_potential.f90 +++ b/src/eDFT/W38_lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine UW38_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) +subroutine W38_lda_correlation_potential(nGrid,weight,nBas,AO,rho,Fc) ! Compute the unrestricted version of the Wigner's LDA correlation potential @@ -73,4 +73,4 @@ include 'parameters.h' Fc(:,:,:) = -4d0*a*Fc(:,:,:) -end subroutine UW38_lda_correlation_potential +end subroutine W38_lda_correlation_potential diff --git a/src/eDFT/allocate_grid.f90 b/src/eDFT/allocate_grid.f90 index 22ab6cf..9bd50d4 100644 --- a/src/eDFT/allocate_grid.f90 +++ b/src/eDFT/allocate_grid.f90 @@ -1,5 +1,4 @@ -subroutine allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent, & - radial_precision,nAng,nGrid) +subroutine allocate_grid(nNuc,ZNuc,max_ang_mom,min_exponent,max_exponent,radial_precision,nAng,nGrid) ! Allocate quadrature grid with numgrid (Radovan Bast) diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/auxiliary_energy.f90 similarity index 90% rename from src/eDFT/unrestricted_auxiliary_energy.f90 rename to src/eDFT/auxiliary_energy.f90 index 00046e8..a5060f6 100644 --- a/src/eDFT/unrestricted_auxiliary_energy.f90 +++ b/src/eDFT/auxiliary_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,Eaux) +subroutine auxiliary_energy(nBas,nEns,eps,occnum,Eaux) ! Compute the auxiliary KS energies @@ -52,4 +52,4 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,Eaux) end do -end subroutine unrestricted_auxiliary_energy +end subroutine auxiliary_energy diff --git a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 b/src/eDFT/correlation_derivative_discontinuity.f90 similarity index 58% rename from src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 rename to src/eDFT/correlation_derivative_discontinuity.f90 index 38e8d1c..9bcd1bf 100644 --- a/src/eDFT/unrestricted_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) +subroutine correlation_derivative_discontinuity(rung,DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) ! Compute the correlation part of the derivative discontinuity @@ -34,26 +34,26 @@ subroutine unrestricted_correlation_derivative_discontinuity(rung,DFA,nEns,wEns, case(1) - call unrestricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + call lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! GGA functionals case(2) - call unrestricted_gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + call gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! MGGA functionals case(3) - call unrestricted_mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + call mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! Hybrid functionals case(4) - call unrestricted_hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) + call hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) end select -end subroutine unrestricted_correlation_derivative_discontinuity +end subroutine correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_correlation_energy.f90 b/src/eDFT/correlation_energy.f90 similarity index 62% rename from src/eDFT/unrestricted_correlation_energy.f90 rename to src/eDFT/correlation_energy.f90 index 61e8ee5..4763db3 100644 --- a/src/eDFT/unrestricted_correlation_energy.f90 +++ b/src/eDFT/correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) +subroutine correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) ! Compute the unrestricted version of the correlation energy @@ -34,26 +34,26 @@ subroutine unrestricted_correlation_energy(rung,DFA,nEns,wEns,nGrid,weight,rho,d case(1) - call unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) + call lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) ! GGA functionals case(2) - call unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) + call gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) ! MGGA functionals case(3) - call unrestricted_mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) + call mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) ! Hybrid functionals case(4) - call unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) + call hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) end select -end subroutine unrestricted_correlation_energy +end subroutine correlation_energy diff --git a/src/eDFT/unrestricted_correlation_individual_energy.f90 b/src/eDFT/correlation_individual_energy.f90 similarity index 75% rename from src/eDFT/unrestricted_correlation_individual_energy.f90 rename to src/eDFT/correlation_individual_energy.f90 index a7f96be..074870e 100644 --- a/src/eDFT/unrestricted_correlation_individual_energy.f90 +++ b/src/eDFT/correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight, & +subroutine correlation_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nGrid,weight, & rhow,drhow,rho,drho,LZc,Ec) ! Compute the correlation energy of individual states @@ -37,7 +37,7 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns case(1) - call unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec) + call lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec) ! GGA functionals @@ -55,8 +55,8 @@ subroutine unrestricted_correlation_individual_energy(rung,DFA,LDA_centered,nEns case(4) - call unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZc,Ec) + call hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZc,Ec) end select -end subroutine unrestricted_correlation_individual_energy +end subroutine correlation_individual_energy diff --git a/src/eDFT/unrestricted_correlation_potential.f90 b/src/eDFT/correlation_potential.f90 similarity index 66% rename from src/eDFT/unrestricted_correlation_potential.f90 rename to src/eDFT/correlation_potential.f90 index 9a488ed..5e1256b 100644 --- a/src/eDFT/unrestricted_correlation_potential.f90 +++ b/src/eDFT/correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute the correlation potential @@ -43,26 +43,26 @@ subroutine unrestricted_correlation_potential(rung,DFA,nEns,wEns,nGrid,weight,nB case(1) - call unrestricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) + call lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) ! GGA functionals case(2) - call unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + call gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! MGGA functionals case(3) - call unrestricted_mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + call mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Hybrid functionals case(4) - call unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + call hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) end select -end subroutine unrestricted_correlation_potential +end subroutine correlation_potential diff --git a/src/eDFT/unrestricted_density_matrix.f90 b/src/eDFT/density_matrix.f90 similarity index 90% rename from src/eDFT/unrestricted_density_matrix.f90 rename to src/eDFT/density_matrix.f90 index 02f1d22..eba8714 100644 --- a/src/eDFT/unrestricted_density_matrix.f90 +++ b/src/eDFT/density_matrix.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_density_matrix(nBas,nEns,c,P,occnum) +subroutine density_matrix(nBas,nEns,c,P,occnum) ! Calculate density matrices @@ -45,4 +45,4 @@ subroutine unrestricted_density_matrix(nBas,nEns,c,P,occnum) -end subroutine unrestricted_density_matrix +end subroutine density_matrix diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index f4a20e8..9ed88dc 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -163,6 +163,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max n_diis = 0 F_diis(:,:,:) = 0d0 err_diis(:,:,:) = 0d0 + rcond(:) = 1d0 !------------------------------------------------------------------------ ! Main SCF loop @@ -184,7 +185,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Compute density matrix !------------------------------------------------------------------------ - call unrestricted_density_matrix(nBas,nEns,c(:,:,:),P(:,:,:,:),occnum(:,:,:)) + call density_matrix(nBas,nEns,c(:,:,:),P(:,:,:,:),occnum(:,:,:)) ! Weight-dependent density matrix @@ -236,20 +237,20 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Compute Hartree potential do ispin=1,nspin - call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) + call hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) end do ! Compute exchange potential do ispin=1,nspin - call unrestricted_exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & + call exchange_potential(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & Pw(:,:,ispin),ERI,AO,dAO,rhow(:,ispin),drhow(:,:,ispin), & Cx_choice,doNcentered,Fx(:,:,ispin),FxHF(:,:,ispin)) end do ! Compute correlation potential - call unrestricted_correlation_potential(c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rhow,drhow,Fc) + call correlation_potential(c_rung,c_DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rhow,drhow,Fc) ! Build Fock operator @@ -268,7 +269,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! DIIS extrapolation n_diis = min(n_diis+1,max_diis) - if(minval(rcond(:)) > 1d-7) then + if(minval(rcond(:)) > 1d-15) then do ispin=1,nspin call DIIS_extrapolation(rcond(ispin),nBasSq,nBasSq,n_diis, & err_diis(:,:,ispin),F_diis(:,:,ispin),err(:,:,ispin),F(:,:,ispin)) @@ -314,19 +315,19 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Hartree energy - call unrestricted_hartree_energy(nBas,Pw,J,EH) + call hartree_energy(nBas,Pw,J,EH) ! Exchange energy do ispin=1,nspin - call unrestricted_exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & + call exchange_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & Pw(:,:,ispin),FxHF(:,:,ispin),rhow(:,ispin),drhow(:,:,ispin), & Cx_choice,doNcentered,Ex(ispin)) end do ! Correlation energy - call unrestricted_correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) + call correlation_energy(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,Ec) ! Total energy @@ -346,6 +347,9 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max end do write(*,*)'------------------------------------------------------------------------------------------' +! print*,'Ensemble energy:',Ew + ENuc,'au' + + !------------------------------------------------------------------------ ! End of SCF loop !------------------------------------------------------------------------ @@ -377,7 +381,7 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,nCC,aCC,nGrid,weight,max ! Compute individual energies from ensemble energy !------------------------------------------------------------------------ - call unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & + call individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & AO,dAO,T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew) end subroutine eDFT_UKS diff --git a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 b/src/eDFT/exchange_derivative_discontinuity.f90 similarity index 66% rename from src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 rename to src/eDFT/exchange_derivative_discontinuity.f90 index b5e3269..bc8485e 100644 --- a/src/eDFT/unrestricted_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,drhow,& +subroutine exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,drhow,& Cx_choice,doNcentered,kappa,ExDD) ! Compute the exchange part of the derivative discontinuity @@ -22,10 +22,10 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC logical,intent(in) :: doNcentered double precision,intent(in) :: kappa(nEns) -! Local variables +!Local variables -! Output variables +!Output variables double precision,intent(out) :: ExDD(nEns) @@ -41,27 +41,27 @@ subroutine unrestricted_exchange_derivative_discontinuity(rung,DFA,nEns,wEns,nCC case(1) - call unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),& + call lda_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),& rhow(:),Cx_choice,doNcentered,kappa,ExDD(:)) ! GGA functionals case(2) - call unrestricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) + call gga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) ! MGGA functionals case(3) - call unrestricted_mgga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) + call mgga_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nGrid,weight(:),rhow(:),drhow(:,:),ExDD(:)) ! Hybrid functionals case(4) - call unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),& + call hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns(:),nCC,aCC,nGrid,weight(:),& rhow(:),Cx_choice,doNcentered,ExDD(:)) end select -end subroutine unrestricted_exchange_derivative_discontinuity +end subroutine exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_exchange_energy.f90 b/src/eDFT/exchange_energy.f90 similarity index 69% rename from src/eDFT/unrestricted_exchange_energy.f90 rename to src/eDFT/exchange_energy.f90 index d61199b..a787910 100644 --- a/src/eDFT/unrestricted_exchange_energy.f90 +++ b/src/eDFT/exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, & +subroutine exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, & rho,drho,Cx_choice,doNcentered,Ex) ! Compute the exchange energy @@ -43,27 +43,27 @@ subroutine unrestricted_exchange_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC, case(1) - call unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) + call lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) ! GGA functionals case(2) - call unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) + call gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex) ! MGGA functionals case(3) - call unrestricted_mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) + call mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) ! Hybrid functionals case(4) - call unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, & + call hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, & rho,drho,Cx_choice,doNcentered,Ex) end select -end subroutine unrestricted_exchange_energy +end subroutine exchange_energy diff --git a/src/eDFT/unrestricted_exchange_individual_energy.f90 b/src/eDFT/exchange_individual_energy.f90 similarity index 72% rename from src/eDFT/unrestricted_exchange_individual_energy.f90 rename to src/eDFT/exchange_individual_energy.f90 index bb29a17..63a025f 100644 --- a/src/eDFT/unrestricted_exchange_individual_energy.f90 +++ b/src/eDFT/exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & +subroutine exchange_individual_energy(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas, & ERI,Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,LZx,Ex) ! Compute the exchange individual energy @@ -45,27 +45,27 @@ subroutine unrestricted_exchange_individual_energy(rung,DFA,LDA_centered,nEns,wE case(1) - call unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,& + call lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,& rhow,rho,Cx_choice,doNcentered,LZx,Ex) ! GGA functionals case(2) - call unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) + call gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) ! MGGA functionals case(3) - call unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) + call mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) ! Hybrid functionals case(4) - call unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow,P,rho,drho,LZx,Ex) + call hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow,P,rho,drho,LZx,Ex) end select -end subroutine unrestricted_exchange_individual_energy +end subroutine exchange_individual_energy diff --git a/src/eDFT/unrestricted_exchange_potential.f90 b/src/eDFT/exchange_potential.f90 similarity index 75% rename from src/eDFT/unrestricted_exchange_potential.f90 rename to src/eDFT/exchange_potential.f90 index 9012c43..3322ab0 100644 --- a/src/eDFT/unrestricted_exchange_potential.f90 +++ b/src/eDFT/exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, & +subroutine exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, & ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF) ! Compute the exchange potential @@ -52,28 +52,29 @@ subroutine unrestricted_exchange_potential(rung,DFA,LDA_centered,nEns,wEns,nCC,a case(1) - call unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,& + call lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,& Cx_choice,doNcentered,Fx) ! GGA functionals case(2) - call unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,& + Cx_choice,Fx) ! MGGA functionals case(3) - call unrestricted_mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Hybrid functionals case(4) - call unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, & + call hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, & ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF) end select -end subroutine unrestricted_exchange_potential +end subroutine exchange_potential diff --git a/src/eDFT/unrestricted_fock_exchange_energy.f90 b/src/eDFT/fock_exchange_energy.f90 similarity index 79% rename from src/eDFT/unrestricted_fock_exchange_energy.f90 rename to src/eDFT/fock_exchange_energy.f90 index fca45fa..17e7d49 100644 --- a/src/eDFT/unrestricted_fock_exchange_energy.f90 +++ b/src/eDFT/fock_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_fock_exchange_energy(nBas,P,Fx,Ex) +subroutine fock_exchange_energy(nBas,P,Fx,Ex) ! Compute the (exact) Fock exchange energy @@ -22,4 +22,4 @@ subroutine unrestricted_fock_exchange_energy(nBas,P,Fx,Ex) Ex = 0.5d0*trace_matrix(nBas,matmul(P,Fx)) -end subroutine unrestricted_fock_exchange_energy +end subroutine fock_exchange_energy diff --git a/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 b/src/eDFT/fock_exchange_individual_energy.f90 similarity index 81% rename from src/eDFT/unrestricted_fock_exchange_individual_energy.f90 rename to src/eDFT/fock_exchange_individual_energy.f90 index 06cf7c3..560576f 100644 --- a/src/eDFT/unrestricted_fock_exchange_individual_energy.f90 +++ b/src/eDFT/fock_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,Ex) +subroutine fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,Ex) ! Compute the HF individual energy in the unrestricted formalism @@ -32,7 +32,7 @@ subroutine unrestricted_fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,E do ispin=1,nspin - call unrestricted_fock_exchange_potential(nBas,Pw(:,:,ispin),ERI,Fx(:,:,ispin)) + call fock_exchange_potential(nBas,Pw(:,:,ispin),ERI,Fx(:,:,ispin)) LZx(ispin) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,ispin),Fx(:,:,ispin))) @@ -43,4 +43,4 @@ subroutine unrestricted_fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,E end do -end subroutine unrestricted_fock_exchange_individual_energy +end subroutine fock_exchange_individual_energy diff --git a/src/eDFT/unrestricted_fock_exchange_potential.f90 b/src/eDFT/fock_exchange_potential.f90 similarity index 83% rename from src/eDFT/unrestricted_fock_exchange_potential.f90 rename to src/eDFT/fock_exchange_potential.f90 index acb0ef7..f1b23f9 100644 --- a/src/eDFT/unrestricted_fock_exchange_potential.f90 +++ b/src/eDFT/fock_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_fock_exchange_potential(nBas,P,ERI,Fx) +subroutine fock_exchange_potential(nBas,P,ERI,Fx) ! Compute the Fock exchange potential @@ -31,4 +31,4 @@ subroutine unrestricted_fock_exchange_potential(nBas,P,ERI,Fx) enddo enddo -end subroutine unrestricted_fock_exchange_potential +end subroutine fock_exchange_potential diff --git a/src/eDFT/unrestricted_gga_correlation_derivative_discontinuity.f90 b/src/eDFT/gga_correlation_derivative_discontinuity.f90 similarity index 81% rename from src/eDFT/unrestricted_gga_correlation_derivative_discontinuity.f90 rename to src/eDFT/gga_correlation_derivative_discontinuity.f90 index de4fcc0..fbf8895 100644 --- a/src/eDFT/unrestricted_gga_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/gga_correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) +subroutine gga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! Compute the correlation GGA part of the derivative discontinuity @@ -41,4 +41,4 @@ subroutine unrestricted_gga_correlation_derivative_discontinuity(DFA,nEns,wEns,n end select -end subroutine unrestricted_gga_correlation_derivative_discontinuity +end subroutine gga_correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_gga_correlation_energy.f90 b/src/eDFT/gga_correlation_energy.f90 similarity index 74% rename from src/eDFT/unrestricted_gga_correlation_energy.f90 rename to src/eDFT/gga_correlation_energy.f90 index 04c0418..ef53764 100644 --- a/src/eDFT/unrestricted_gga_correlation_energy.f90 +++ b/src/eDFT/gga_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) +subroutine gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) ! Compute unrestricted GGA correlation energy @@ -28,11 +28,11 @@ subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,dr case (1) - call ULYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) + call LYP_gga_correlation_energy(nGrid,weight,rho,drho,Ec) case (2) - call UPBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec) + call PBE_gga_correlation_energy(nGrid,weight,rho,drho,Ec) case default @@ -41,4 +41,4 @@ subroutine unrestricted_gga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,dr end select -end subroutine unrestricted_gga_correlation_energy +end subroutine gga_correlation_energy diff --git a/src/eDFT/unrestricted_gga_correlation_potential.f90 b/src/eDFT/gga_correlation_potential.f90 similarity index 73% rename from src/eDFT/unrestricted_gga_correlation_potential.f90 rename to src/eDFT/gga_correlation_potential.f90 index 9d004dc..d7042dc 100644 --- a/src/eDFT/unrestricted_gga_correlation_potential.f90 +++ b/src/eDFT/gga_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute unrestricted GGA correlation potential @@ -30,11 +30,11 @@ subroutine unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBa case (1) - call ULYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + call LYP_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) case (2) - call UPBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + call PBE_gga_correlation_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fc) case default @@ -43,4 +43,4 @@ subroutine unrestricted_gga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBa end select -end subroutine unrestricted_gga_correlation_potential +end subroutine gga_correlation_potential diff --git a/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 b/src/eDFT/gga_exchange_derivative_discontinuity.f90 similarity index 82% rename from src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 rename to src/eDFT/gga_exchange_derivative_discontinuity.f90 index 5fc2c20..4373930 100644 --- a/src/eDFT/unrestricted_gga_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/gga_exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) +subroutine gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) ! Compute the exchange GGA part of the derivative discontinuity @@ -45,4 +45,4 @@ subroutine unrestricted_gga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGri end select -end subroutine unrestricted_gga_exchange_derivative_discontinuity +end subroutine gga_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_gga_exchange_energy.f90 b/src/eDFT/gga_exchange_energy.f90 similarity index 53% rename from src/eDFT/unrestricted_gga_exchange_energy.f90 rename to src/eDFT/gga_exchange_energy.f90 index cfbcf1c..fba9cd0 100644 --- a/src/eDFT/unrestricted_gga_exchange_energy.f90 +++ b/src/eDFT/gga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) +subroutine gga_exchange_energy(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,Cx_choice,Ex) ! Select GGA exchange functional for energy calculation @@ -11,11 +11,15 @@ subroutine unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho, integer,intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nCC + double precision,intent(in) :: aCC(nCC,nEns-1) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) double precision,intent(in) :: rho(nGrid) + integer,intent(in) :: Cx_choice double precision,intent(in) :: drho(ncart,nGrid) + ! Output variables double precision :: Ex @@ -24,15 +28,20 @@ subroutine unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho, case (1) - call UG96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + call G96_gga_exchange_energy(nGrid,weight,rho,drho,Ex) case (2) - call UB88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + call B88_gga_exchange_energy(nGrid,weight,rho,drho,Ex) case (3) - call UPBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + call PBE_gga_exchange_energy(nGrid,weight,rho,drho,Ex) + + case (4) + + call CC_B88_gga_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,drho,& + Cx_choice,Ex) case default @@ -41,4 +50,4 @@ subroutine unrestricted_gga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho, end select -end subroutine unrestricted_gga_exchange_energy +end subroutine gga_exchange_energy diff --git a/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 b/src/eDFT/gga_exchange_individual_energy.f90 similarity index 83% rename from src/eDFT/unrestricted_gga_exchange_individual_energy.f90 rename to src/eDFT/gga_exchange_individual_energy.f90 index f43e0ac..4b22aeb 100644 --- a/src/eDFT/unrestricted_gga_exchange_individual_energy.f90 +++ b/src/eDFT/gga_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) +subroutine gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) ! Compute GGA exchange energy for individual states @@ -33,4 +33,4 @@ subroutine unrestricted_gga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weigh end select -end subroutine unrestricted_gga_exchange_individual_energy +end subroutine gga_exchange_individual_energy diff --git a/src/eDFT/unrestricted_gga_exchange_potential.f90 b/src/eDFT/gga_exchange_potential.f90 similarity index 55% rename from src/eDFT/unrestricted_gga_exchange_potential.f90 rename to src/eDFT/gga_exchange_potential.f90 index 57de1e5..81f05a3 100644 --- a/src/eDFT/unrestricted_gga_exchange_potential.f90 +++ b/src/eDFT/gga_exchange_potential.f90 @@ -1,4 +1,5 @@ -subroutine unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine gga_exchange_potential(DFA,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,& + rho,drho,Cx_choice,Fx) ! Select GGA exchange functional for potential calculation @@ -10,6 +11,8 @@ subroutine unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,A integer,intent(in) :: DFA integer,intent(in) :: nEns double precision,intent(in) :: wEns(nEns) + integer,intent(in) :: nCC + double precision,intent(in) :: aCC(nCC,nEns-1) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: nBas @@ -17,6 +20,7 @@ subroutine unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,A double precision,intent(in) :: dAO(3,nBas,nGrid) double precision,intent(in) :: rho(nGrid) double precision,intent(in) :: drho(3,nGrid) + integer,intent(in) :: Cx_choice ! Output variables @@ -28,15 +32,20 @@ subroutine unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,A case (1) - call UG96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call G96_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) case (2) - call UB88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call B88_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) case (3) - call UPBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + call PBE_gga_exchange_potential(nGrid,weight,nBas,AO,dAO,rho,drho,Fx) + + case (4) + + call CC_B88_gga_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO,rho,drho,& + Cx_choice,Fx) case default @@ -45,4 +54,4 @@ subroutine unrestricted_gga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,A end select -end subroutine unrestricted_gga_exchange_potential +end subroutine gga_exchange_potential diff --git a/src/eDFT/unrestricted_hartree_energy.f90 b/src/eDFT/hartree_energy.f90 similarity index 87% rename from src/eDFT/unrestricted_hartree_energy.f90 rename to src/eDFT/hartree_energy.f90 index 9937809..d3f04a3 100644 --- a/src/eDFT/unrestricted_hartree_energy.f90 +++ b/src/eDFT/hartree_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hartree_energy(nBas,P,J,EH) +subroutine hartree_energy(nBas,P,J,EH) ! Compute the unrestricted version of the Hartree energy @@ -26,4 +26,4 @@ subroutine unrestricted_hartree_energy(nBas,P,J,EH) + 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,1))) EH(3) = 0.5d0*trace_matrix(nBas,matmul(P(:,:,2),J(:,:,2))) -end subroutine unrestricted_hartree_energy +end subroutine hartree_energy diff --git a/src/eDFT/unrestricted_hartree_individual_energy.f90 b/src/eDFT/hartree_individual_energy.f90 similarity index 86% rename from src/eDFT/unrestricted_hartree_individual_energy.f90 rename to src/eDFT/hartree_individual_energy.f90 index 04b7bcf..5c0b649 100644 --- a/src/eDFT/unrestricted_hartree_individual_energy.f90 +++ b/src/eDFT/hartree_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) +subroutine hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) ! Compute the hartree contribution to the individual energies @@ -35,7 +35,7 @@ subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) EH(:,:) = 0.d0 do ispin=1,nspin - call unrestricted_hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) + call hartree_potential(nBas,Pw(:,:,ispin),ERI,J(:,:,ispin)) end do LZH(1) = - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) @@ -52,4 +52,4 @@ subroutine unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) end do -end subroutine unrestricted_hartree_individual_energy +end subroutine hartree_individual_energy diff --git a/src/eDFT/unrestricted_hartree_potential.f90 b/src/eDFT/hartree_potential.f90 similarity index 84% rename from src/eDFT/unrestricted_hartree_potential.f90 rename to src/eDFT/hartree_potential.f90 index 5860e02..0aacdd2 100644 --- a/src/eDFT/unrestricted_hartree_potential.f90 +++ b/src/eDFT/hartree_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hartree_potential(nBas,P,ERI,J) +subroutine hartree_potential(nBas,P,ERI,J) ! Compute the unrestricted version of the Hartree potential @@ -30,4 +30,4 @@ subroutine unrestricted_hartree_potential(nBas,P,ERI,J) enddo -end subroutine unrestricted_hartree_potential +end subroutine hartree_potential diff --git a/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 b/src/eDFT/hybrid_correlation_derivative_discontinuity.f90 similarity index 81% rename from src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 rename to src/eDFT/hybrid_correlation_derivative_discontinuity.f90 index 09e608b..40a3bb3 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/hybrid_correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) +subroutine hybrid_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! Compute the correlation hybrid part of the derivative discontinuity @@ -43,4 +43,4 @@ subroutine unrestricted_hybrid_correlation_derivative_discontinuity(DFA,nEns,wEn end select -end subroutine unrestricted_hybrid_correlation_derivative_discontinuity +end subroutine hybrid_correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_hybrid_correlation_energy.f90 b/src/eDFT/hybrid_correlation_energy.f90 similarity index 66% rename from src/eDFT/unrestricted_hybrid_correlation_energy.f90 rename to src/eDFT/hybrid_correlation_energy.f90 index f54d7a3..4b272ca 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_energy.f90 +++ b/src/eDFT/hybrid_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) +subroutine hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) ! Compute the unrestricted version of the correlation energy for hybrid functionals @@ -35,18 +35,18 @@ subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho aC = 0.81d0 - call unrestricted_lda_correlation_energy(3,nEns,wEns,nGrid,weight,rho,EcLDA) - call unrestricted_gga_correlation_energy(1,nEns,wEns,nGrid,weight,rho,drho,EcGGA) + call lda_correlation_energy(3,nEns,wEns,nGrid,weight,rho,EcLDA) + call gga_correlation_energy(1,nEns,wEns,nGrid,weight,rho,drho,EcGGA) Ec(:) = EcLDA(:) + aC*(EcGGA(:) - EcLDA(:)) case(3) - call unrestricted_gga_correlation_energy(1,nEns,wEns,nGrid,weight,rho,drho,Ec) + call gga_correlation_energy(1,nEns,wEns,nGrid,weight,rho,drho,Ec) case(4) - call unrestricted_gga_correlation_energy(2,nEns,wEns,nGrid,weight,rho,drho,Ec) + call gga_correlation_energy(2,nEns,wEns,nGrid,weight,rho,drho,Ec) case default @@ -55,4 +55,4 @@ subroutine unrestricted_hybrid_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho end select -end subroutine unrestricted_hybrid_correlation_energy +end subroutine hybrid_correlation_energy diff --git a/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 b/src/eDFT/hybrid_correlation_individual_energy.f90 similarity index 86% rename from src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 rename to src/eDFT/hybrid_correlation_individual_energy.f90 index 54b73f2..6a2214d 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_individual_energy.f90 +++ b/src/eDFT/hybrid_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight, & +subroutine hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid,weight, & rhow,drhow,rho,drho,LZc,Ec) ! Compute the hybrid correlation energy for individual states @@ -39,4 +39,4 @@ subroutine unrestricted_hybrid_correlation_individual_energy(DFA,nEns,wEns,nGrid end select -end subroutine unrestricted_hybrid_correlation_individual_energy +end subroutine hybrid_correlation_individual_energy diff --git a/src/eDFT/unrestricted_hybrid_correlation_potential.f90 b/src/eDFT/hybrid_correlation_potential.f90 similarity index 69% rename from src/eDFT/unrestricted_hybrid_correlation_potential.f90 rename to src/eDFT/hybrid_correlation_potential.f90 index 815fbbd..9104a80 100644 --- a/src/eDFT/unrestricted_hybrid_correlation_potential.f90 +++ b/src/eDFT/hybrid_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute the correlation potential for hybrid functionals @@ -42,8 +42,8 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight, aC = 0.81d0 - call unrestricted_lda_correlation_potential(3,nEns,wEns,nGrid,weight,nBas,AO,rho,FcLDA) - call unrestricted_gga_correlation_potential(1,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA) + call lda_correlation_potential(3,nEns,wEns,nGrid,weight,nBas,AO,rho,FcLDA) + call gga_correlation_potential(1,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FcGGA) Fc(:,:,:) = FcLDA(:,:,:) + aC*(FcGGA(:,:,:) - FcLDA(:,:,:)) @@ -51,13 +51,13 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight, allocate(FcGGA(nBas,nBas,nspin)) - call unrestricted_gga_correlation_potential(1,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + call gga_correlation_potential(1,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) case(4) allocate(FcGGA(nBas,nBas,nspin)) - call unrestricted_gga_correlation_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) + call gga_correlation_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) case default @@ -66,4 +66,4 @@ subroutine unrestricted_hybrid_correlation_potential(DFA,nEns,wEns,nGrid,weight, end select -end subroutine unrestricted_hybrid_correlation_potential +end subroutine hybrid_correlation_potential diff --git a/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 b/src/eDFT/hybrid_exchange_derivative_discontinuity.f90 similarity index 85% rename from src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 rename to src/eDFT/hybrid_exchange_derivative_discontinuity.f90 index 8fe8e55..243d85c 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/hybrid_exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& +subroutine hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& Cx_choice,doNcentered,ExDD) ! Compute the exchange part of the derivative discontinuity for hybrid functionals @@ -50,4 +50,4 @@ subroutine unrestricted_hybrid_exchange_derivative_discontinuity(DFA,nEns,wEns,n end select -end subroutine unrestricted_hybrid_exchange_derivative_discontinuity +end subroutine hybrid_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_hybrid_exchange_energy.f90 b/src/eDFT/hybrid_exchange_energy.f90 similarity index 67% rename from src/eDFT/unrestricted_hybrid_exchange_energy.f90 rename to src/eDFT/hybrid_exchange_energy.f90 index 365722e..a584f4c 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_energy.f90 +++ b/src/eDFT/hybrid_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, & +subroutine hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P,FxHF, & rho,drho,Cx_choice,doNcentered,Ex) ! Compute the exchange energy for hybrid functionals @@ -37,17 +37,17 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aC case (1) - call unrestricted_fock_exchange_energy(nBas,P,FxHF,Ex) + call fock_exchange_energy(nBas,P,FxHF,Ex) case (2) a0 = 0.20d0 aX = 0.72d0 - call unrestricted_lda_exchange_energy(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,& + call lda_exchange_energy(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,& rho,Cx_choice,doNcentered,ExLDA) - call unrestricted_gga_exchange_energy(2,nEns,wEns,nGrid,weight,rho,drho,ExGGA) - call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + call gga_exchange_energy(2,nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call fock_exchange_energy(nBas,P,FxHF,ExHF) Ex = ExLDA & + a0*(ExHF - ExLDA) & @@ -55,15 +55,15 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aC case (3) - call unrestricted_gga_exchange_energy(2,nEns,wEns,nGrid,weight,rho,drho,ExGGA) - call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + call gga_exchange_energy(2,nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call fock_exchange_energy(nBas,P,FxHF,ExHF) Ex = 0.5d0*ExHF + 0.5d0*ExGGA case (4) - call unrestricted_gga_exchange_energy(3,nEns,wEns,nGrid,weight,rho,drho,ExGGA) - call unrestricted_fock_exchange_energy(nBas,P,FxHF,ExHF) + call gga_exchange_energy(3,nEns,wEns,nGrid,weight,rho,drho,ExGGA) + call fock_exchange_energy(nBas,P,FxHF,ExHF) Ex = 0.25d0*ExHF + 0.75d0*ExGGA @@ -74,4 +74,4 @@ subroutine unrestricted_hybrid_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aC end select -end subroutine unrestricted_hybrid_exchange_energy +end subroutine hybrid_exchange_energy diff --git a/src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90 b/src/eDFT/hybrid_exchange_individual_energy.f90 similarity index 81% rename from src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90 rename to src/eDFT/hybrid_exchange_individual_energy.f90 index bef5554..f3b4199 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_individual_energy.f90 +++ b/src/eDFT/hybrid_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow, & +subroutine hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,nBas,ERI,Pw,rhow,drhow, & P,rho,drho,LZx,Ex) ! Compute the hybrid exchange energy for individual states @@ -34,7 +34,7 @@ subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,we case (1) - call unrestricted_fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,Ex) + call fock_exchange_individual_energy(nBas,nEns,Pw,P,ERI,LZx,Ex) case default @@ -43,4 +43,4 @@ subroutine unrestricted_hybrid_exchange_individual_energy(DFA,nEns,wEns,nGrid,we end select -end subroutine unrestricted_hybrid_exchange_individual_energy +end subroutine hybrid_exchange_individual_energy diff --git a/src/eDFT/unrestricted_hybrid_exchange_potential.f90 b/src/eDFT/hybrid_exchange_potential.f90 similarity index 71% rename from src/eDFT/unrestricted_hybrid_exchange_potential.f90 rename to src/eDFT/hybrid_exchange_potential.f90 index c214eb5..5dfb81c 100644 --- a/src/eDFT/unrestricted_hybrid_exchange_potential.f90 +++ b/src/eDFT/hybrid_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, & +subroutine hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,P, & ERI,AO,dAO,rho,drho,Cx_choice,doNcentered,Fx,FxHF) ! Compute the exchange potential for hybrid functionals @@ -44,7 +44,7 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC case(1) - call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + call fock_exchange_potential(nBas,P,ERI,FxHF) Fx(:,:) = FxHF(:,:) case(2) @@ -54,10 +54,10 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC a0 = 0.20d0 aX = 0.72d0 - call unrestricted_lda_exchange_potential(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight, & + call lda_exchange_potential(1,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight, & nBas,AO,rho,Cx_choice,doNcentered,FxLDA) - call unrestricted_gga_exchange_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) - call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + call gga_exchange_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call fock_exchange_potential(nBas,P,ERI,FxHF) Fx(:,:) = FxLDA(:,:) & + a0*(FxHF(:,:) - FxLDA(:,:)) & @@ -67,8 +67,8 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC allocate(FxGGA(nBas,nBas)) - call unrestricted_gga_exchange_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) - call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + call gga_exchange_potential(2,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call fock_exchange_potential(nBas,P,ERI,FxHF) Fx(:,:) = 0.5d0*FxHF(:,:) + 0.5d0*FxGGA(:,:) @@ -76,8 +76,8 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC allocate(FxGGA(nBas,nBas)) - call unrestricted_gga_exchange_potential(3,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) - call unrestricted_fock_exchange_potential(nBas,P,ERI,FxHF) + call gga_exchange_potential(3,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,FxGGA) + call fock_exchange_potential(nBas,P,ERI,FxHF) Fx(:,:) = 0.25d0*FxHF(:,:) + 0.75d0*FxGGA(:,:) @@ -88,4 +88,4 @@ subroutine unrestricted_hybrid_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC end select -end subroutine unrestricted_hybrid_exchange_potential +end subroutine hybrid_exchange_potential diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/individual_energy.f90 similarity index 85% rename from src/eDFT/unrestricted_individual_energy.f90 rename to src/eDFT/individual_energy.f90 index 9af4ed7..b546bab 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO, & +subroutine individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,dAO, & T,V,ERI,ENuc,eKS,Pw,rhow,drhow,J,Fx,FxHF,Fc,P,rho,drho,occnum,Cx_choice,doNcentered,Ew) ! Compute unrestricted individual energies as well as excitation energies @@ -118,7 +118,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered LZH(:) = 0d0 EH(:,:) = 0d0 - call unrestricted_hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) + call hartree_individual_energy(nBas,nEns,Pw,P,ERI,LZH,EH) !------------------------------------------------------------------------ ! Individual exchange energy @@ -126,7 +126,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered LZx(:) = 0d0 Ex(:,:) = 0d0 - call unrestricted_exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, & + call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,ERI, & Pw,rhow,drhow,P,rho,drho,Cx_choice,doNcentered,LZx,Ex) !------------------------------------------------------------------------ @@ -135,26 +135,37 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered LZc(:) = 0d0 Ec(:,:) = 0d0 - call unrestricted_correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, & + call correlation_individual_energy(c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, & rhow,drhow,rho,drho,LZc,Ec) !------------------------------------------------------------------------ ! Compute auxiliary energies !------------------------------------------------------------------------ - call unrestricted_auxiliary_energy(nBas,nEns,eKS,occnum,Eaux) + call auxiliary_energy(nBas,nEns,eKS,occnum,Eaux) !------------------------------------------------------------------------ ! Compute derivative discontinuities !------------------------------------------------------------------------ do ispin=1,nspin - call unrestricted_exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,nCC,aCC,nGrid,weight, & - rhow(:,ispin),drhow(:,:,ispin),Cx_choice,doNcentered,kappa,ExDD(ispin,:)) + call exchange_derivative_discontinuity(x_rung,x_DFA,nEns,wEns,nCC,aCC,nGrid,weight, & + rhow(:,ispin),drhow(:,:,ispin),Cx_choice,& + doNcentered,kappa,ExDD(ispin,:)) end do - call unrestricted_correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD) + call correlation_derivative_discontinuity(c_rung,c_DFA,nEns,wEns,nGrid,weight,rhow,drhow,EcDD) +! Scaling derivative discontinuity for N-centered ensembles + +! if(doNcentered) then + +! do iEns=1,nEns +! ExDD(:,iEns) = (1d0 - kappa(iEns))*ExDD(:,iEns) +! EcDD(:,iEns) = (1d0 - kappa(iEns))*EcDD(:,iEns) +! end do + +! end if !------------------------------------------------------------------------ ! Total energy @@ -177,9 +188,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered end do end if - - print*,"LZ shift:",sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)),"au" - + +print*,'LZshift=',sum(LZH(:)) + sum(LZx(:)) + sum(LZc(:)) + ! do iEns=1,nEns ! E(iEns) = sum(ET(:,iEns)) + sum(EV(:,iEns)) & ! + sum(EH(:,iEns)) + sum(Ex(:,iEns)) + sum(Ec(:,iEns)) & @@ -224,7 +235,7 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered ! Dump results !------------------------------------------------------------------------ - call print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, & + call print_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, & Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD) -end subroutine unrestricted_individual_energy +end subroutine individual_energy diff --git a/src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/lda_correlation_derivative_discontinuity.f90 similarity index 83% rename from src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90 rename to src/eDFT/lda_correlation_derivative_discontinuity.f90 index e0ea35d..3229bf4 100644 --- a/src/eDFT/unrestricted_lda_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/lda_correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) +subroutine lda_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! Compute the correlation LDA part of the derivative discontinuity @@ -49,4 +49,4 @@ subroutine unrestricted_lda_correlation_derivative_discontinuity(DFA,nEns,wEns,n end select -end subroutine unrestricted_lda_correlation_derivative_discontinuity +end subroutine lda_correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_lda_correlation_energy.f90 b/src/eDFT/lda_correlation_energy.f90 similarity index 62% rename from src/eDFT/unrestricted_lda_correlation_energy.f90 rename to src/eDFT/lda_correlation_energy.f90 index 7b43cf2..d6bb7ee 100644 --- a/src/eDFT/unrestricted_lda_correlation_energy.f90 +++ b/src/eDFT/lda_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) +subroutine lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec) ! Select the unrestrited version of the LDA correlation functional @@ -24,23 +24,23 @@ subroutine unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec case (1) - call UW38_lda_correlation_energy(nGrid,weight,rho,Ec) + call W38_lda_correlation_energy(nGrid,weight,rho,Ec) case (2) - call UPW92_lda_correlation_energy(nGrid,weight,rho,Ec) + call PW92_lda_correlation_energy(nGrid,weight,rho,Ec) case (3) - call UVWN3_lda_correlation_energy(nGrid,weight,rho,Ec) + call VWN3_lda_correlation_energy(nGrid,weight,rho,Ec) case (4) - call UVWN5_lda_correlation_energy(nGrid,weight,rho,Ec) + call VWN5_lda_correlation_energy(nGrid,weight,rho,Ec) case (5) - call UC16_lda_correlation_energy(nGrid,weight,rho,Ec) + call C16_lda_correlation_energy(nGrid,weight,rho,Ec) case default @@ -49,4 +49,4 @@ subroutine unrestricted_lda_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,Ec end select -end subroutine unrestricted_lda_correlation_energy +end subroutine lda_correlation_energy diff --git a/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 b/src/eDFT/lda_correlation_individual_energy.f90 similarity index 63% rename from src/eDFT/unrestricted_lda_correlation_individual_energy.f90 rename to src/eDFT/lda_correlation_individual_energy.f90 index 68bec93..98b3b4f 100644 --- a/src/eDFT/unrestricted_lda_correlation_individual_energy.f90 +++ b/src/eDFT/lda_correlation_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec) +subroutine lda_correlation_individual_energy(DFA,LDA_centered,nEns,wEns,nGrid,weight,rhow,rho,LZc,Ec) ! Compute LDA correlation energy for individual states @@ -27,19 +27,19 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns, case (1) -! call UW38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec) +! call W38_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec) case (2) -! call UPW92_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec) +! call PW92_lda_correlation_individual_energy(nGrid,weight,rhow,rho,LZc,Ec) case (3) -! call UVWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) +! call VWN3_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) case (4) - call UVWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) + call VWN5_lda_correlation_individual_energy(nEns,nGrid,weight,rhow,rho,LZc,Ec) case default @@ -48,4 +48,4 @@ subroutine unrestricted_lda_correlation_individual_energy(DFA,LDA_centered,nEns, end select -end subroutine unrestricted_lda_correlation_individual_energy +end subroutine lda_correlation_individual_energy diff --git a/src/eDFT/unrestricted_lda_correlation_potential.f90 b/src/eDFT/lda_correlation_potential.f90 similarity index 60% rename from src/eDFT/unrestricted_lda_correlation_potential.f90 rename to src/eDFT/lda_correlation_potential.f90 index fa88e9e..e50b033 100644 --- a/src/eDFT/unrestricted_lda_correlation_potential.f90 +++ b/src/eDFT/lda_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) +subroutine lda_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,rho,Fc) ! Select LDA correlation potential @@ -28,23 +28,23 @@ include 'parameters.h' case (1) - call UW38_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call W38_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) case (2) - call UPW92_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call PW92_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) case (3) - call UVWN3_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call VWN3_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) case (4) - call UVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call VWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) case (5) - call UC16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) + call C16_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:,:),Fc(:,:,:)) case default @@ -53,4 +53,4 @@ include 'parameters.h' end select -end subroutine unrestricted_lda_correlation_potential +end subroutine lda_correlation_potential diff --git a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/lda_exchange_derivative_discontinuity.f90 similarity index 76% rename from src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 rename to src/eDFT/lda_exchange_derivative_discontinuity.f90 index d33af22..5065d69 100644 --- a/src/eDFT/unrestricted_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/lda_exchange_derivative_discontinuity.f90 @@ -1,5 +1,5 @@ -subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& - Cx_choice,doNcentered,kappa,ExDD) +subroutine lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& + Cx_choice,doNcentered,kappa,ExDD) ! Compute the exchange LDA part of the derivative discontinuity @@ -38,7 +38,7 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC, case (2) - call UCC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,& + call CC_lda_exchange_derivative_discontinuity(nEns,wEns,nCC,aCC,nGrid,weight,rhow,& Cx_choice,doNcentered,kappa,ExDD) case default @@ -48,4 +48,4 @@ subroutine unrestricted_lda_exchange_derivative_discontinuity(DFA,nEns,wEns,nCC, end select -end subroutine unrestricted_lda_exchange_derivative_discontinuity +end subroutine lda_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_lda_exchange_energy.f90 b/src/eDFT/lda_exchange_energy.f90 similarity index 72% rename from src/eDFT/unrestricted_lda_exchange_energy.f90 rename to src/eDFT/lda_exchange_energy.f90 index 9d4e8ec..e92b0cc 100644 --- a/src/eDFT/unrestricted_lda_exchange_energy.f90 +++ b/src/eDFT/lda_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) +subroutine lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) ! Select LDA exchange functional @@ -29,11 +29,11 @@ subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,n case (1) - call US51_lda_exchange_energy(nGrid,weight,rho,Ex) + call S51_lda_exchange_energy(nGrid,weight,rho,Ex) case (2) - call UCC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) + call CC_lda_exchange_energy(nEns,wEns,nCC,aCC,nGrid,weight,rho,Cx_choice,doNcentered,Ex) case default @@ -43,4 +43,4 @@ subroutine unrestricted_lda_exchange_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,n end select -end subroutine unrestricted_lda_exchange_energy +end subroutine lda_exchange_energy diff --git a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 b/src/eDFT/lda_exchange_individual_energy.f90 similarity index 77% rename from src/eDFT/unrestricted_lda_exchange_individual_energy.f90 rename to src/eDFT/lda_exchange_individual_energy.f90 index 67600b8..43cdf86 100644 --- a/src/eDFT/unrestricted_lda_exchange_individual_energy.f90 +++ b/src/eDFT/lda_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& +subroutine lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,rhow,& rho,Cx_choice,doNcentered,LZx,Ex) ! Compute LDA exchange energy for individual states @@ -32,11 +32,11 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn case (1) - call US51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex) + call S51_lda_exchange_individual_energy(nEns,nGrid,weight,rhow,rho,LZx,Ex) case (2) - call UCC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho, & + call CC_lda_exchange_individual_energy(nEns,wEns,nCC,aCC,nGrid,weight,rhow,rho, & Cx_choice,doNcentered,LZx,Ex) case default @@ -46,4 +46,4 @@ subroutine unrestricted_lda_exchange_individual_energy(DFA,LDA_centered,nEns,wEn end select -end subroutine unrestricted_lda_exchange_individual_energy +end subroutine lda_exchange_individual_energy diff --git a/src/eDFT/unrestricted_lda_exchange_potential.f90 b/src/eDFT/lda_exchange_potential.f90 similarity index 75% rename from src/eDFT/unrestricted_lda_exchange_potential.f90 rename to src/eDFT/lda_exchange_potential.f90 index 60c2c75..0abc0d1 100644 --- a/src/eDFT/unrestricted_lda_exchange_potential.f90 +++ b/src/eDFT/lda_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho, & +subroutine lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho, & Cx_choice,doNcentered,Fx) ! Select LDA correlation potential @@ -33,11 +33,11 @@ subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aC case (1) - call US51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) + call S51_lda_exchange_potential(nGrid,weight,nBas,AO,rho,Fx) case (2) - call UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Cx_choice,doNcentered,Fx) + call CC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho,Cx_choice,doNcentered,Fx) case default @@ -46,4 +46,4 @@ subroutine unrestricted_lda_exchange_potential(DFA,LDA_centered,nEns,wEns,nCC,aC end select -end subroutine unrestricted_lda_exchange_potential +end subroutine lda_exchange_potential diff --git a/src/eDFT/unrestricted_mgga_correlation_derivative_discontinuity.f90 b/src/eDFT/mgga_correlation_derivative_discontinuity.f90 similarity index 79% rename from src/eDFT/unrestricted_mgga_correlation_derivative_discontinuity.f90 rename to src/eDFT/mgga_correlation_derivative_discontinuity.f90 index 20e62c9..9050407 100644 --- a/src/eDFT/unrestricted_mgga_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/mgga_correlation_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) +subroutine mgga_correlation_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,Ec) ! Compute the correlation MGGA part of the derivative discontinuity @@ -31,4 +31,4 @@ subroutine unrestricted_mgga_correlation_derivative_discontinuity(DFA,nEns,wEns, end select -end subroutine unrestricted_mgga_correlation_derivative_discontinuity +end subroutine mgga_correlation_derivative_discontinuity diff --git a/src/eDFT/unrestricted_mgga_correlation_energy.f90 b/src/eDFT/mgga_correlation_energy.f90 similarity index 84% rename from src/eDFT/unrestricted_mgga_correlation_energy.f90 rename to src/eDFT/mgga_correlation_energy.f90 index 09a2901..d6b1f7a 100644 --- a/src/eDFT/unrestricted_mgga_correlation_energy.f90 +++ b/src/eDFT/mgga_correlation_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) +subroutine mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ec) ! Compute unrestricted MGGA correlation energy @@ -33,4 +33,4 @@ subroutine unrestricted_mgga_correlation_energy(DFA,nEns,wEns,nGrid,weight,rho,d end select -end subroutine unrestricted_mgga_correlation_energy +end subroutine mgga_correlation_energy diff --git a/src/eDFT/unrestricted_mgga_correlation_potential.f90 b/src/eDFT/mgga_correlation_potential.f90 similarity index 84% rename from src/eDFT/unrestricted_mgga_correlation_potential.f90 rename to src/eDFT/mgga_correlation_potential.f90 index 6cf959b..642fd8f 100644 --- a/src/eDFT/unrestricted_mgga_correlation_potential.f90 +++ b/src/eDFT/mgga_correlation_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) +subroutine mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fc) ! Compute unrestricted MGGA correlation potential @@ -35,4 +35,4 @@ subroutine unrestricted_mgga_correlation_potential(DFA,nEns,wEns,nGrid,weight,nB end select -end subroutine unrestricted_mgga_correlation_potential +end subroutine mgga_correlation_potential diff --git a/src/eDFT/unrestricted_mgga_exchange_derivative_discontinuity.f90 b/src/eDFT/mgga_exchange_derivative_discontinuity.f90 similarity index 80% rename from src/eDFT/unrestricted_mgga_exchange_derivative_discontinuity.f90 rename to src/eDFT/mgga_exchange_derivative_discontinuity.f90 index 4587ca1..5131d20 100644 --- a/src/eDFT/unrestricted_mgga_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/mgga_exchange_derivative_discontinuity.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_mgga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) +subroutine mgga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGrid,weight,rhow,drhow,ExDD) ! Compute the exchange MGGA part of the derivative discontinuity @@ -33,4 +33,4 @@ subroutine unrestricted_mgga_exchange_derivative_discontinuity(DFA,nEns,wEns,nGr end select -end subroutine unrestricted_mgga_exchange_derivative_discontinuity +end subroutine mgga_exchange_derivative_discontinuity diff --git a/src/eDFT/unrestricted_mgga_exchange_energy.f90 b/src/eDFT/mgga_exchange_energy.f90 similarity index 82% rename from src/eDFT/unrestricted_mgga_exchange_energy.f90 rename to src/eDFT/mgga_exchange_energy.f90 index 7347d67..c8f83af 100644 --- a/src/eDFT/unrestricted_mgga_exchange_energy.f90 +++ b/src/eDFT/mgga_exchange_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) +subroutine mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho,Ex) ! Select MGGA exchange functional for energy calculation @@ -29,4 +29,4 @@ subroutine unrestricted_mgga_exchange_energy(DFA,nEns,wEns,nGrid,weight,rho,drho end select -end subroutine unrestricted_mgga_exchange_energy +end subroutine mgga_exchange_energy diff --git a/src/eDFT/unrestricted_mgga_exchange_individual_energy.f90 b/src/eDFT/mgga_exchange_individual_energy.f90 similarity index 83% rename from src/eDFT/unrestricted_mgga_exchange_individual_energy.f90 rename to src/eDFT/mgga_exchange_individual_energy.f90 index edcdaff..08c837c 100644 --- a/src/eDFT/unrestricted_mgga_exchange_individual_energy.f90 +++ b/src/eDFT/mgga_exchange_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) +subroutine mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weight,rhow,drhow,rho,drho,LZx,Ex) ! Compute MGGA exchange energy for individual states @@ -33,4 +33,4 @@ subroutine unrestricted_mgga_exchange_individual_energy(DFA,nEns,wEns,nGrid,weig end select -end subroutine unrestricted_mgga_exchange_individual_energy +end subroutine mgga_exchange_individual_energy diff --git a/src/eDFT/unrestricted_mgga_exchange_potential.f90 b/src/eDFT/mgga_exchange_potential.f90 similarity index 84% rename from src/eDFT/unrestricted_mgga_exchange_potential.f90 rename to src/eDFT/mgga_exchange_potential.f90 index b6707b4..15a6e12 100644 --- a/src/eDFT/unrestricted_mgga_exchange_potential.f90 +++ b/src/eDFT/mgga_exchange_potential.f90 @@ -1,4 +1,4 @@ -subroutine unrestricted_mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) +subroutine mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas,AO,dAO,rho,drho,Fx) ! Select MGGA exchange functional for potential calculation @@ -33,4 +33,4 @@ subroutine unrestricted_mgga_exchange_potential(DFA,nEns,wEns,nGrid,weight,nBas, end select -end subroutine unrestricted_mgga_exchange_potential +end subroutine mgga_exchange_potential diff --git a/src/eDFT/print_UKS.f90 b/src/eDFT/print_UKS.f90 index cb56dda..0d41965 100644 --- a/src/eDFT/print_UKS.f90 +++ b/src/eDFT/print_UKS.f90 @@ -66,7 +66,7 @@ subroutine print_UKS(nBas,nEns,occnum,Ov,wEns,eps,c,ENuc,ET,EV,EH,Ex,Ec,Ew,dipol HOMOb = -huge(0d0) if(iHOMOb > 0) HOMOb = eps(iHOMOb,2) LUMOb = +huge(0d0) - if(iLUMOb <= nBas) LUMOb = eps(iLUMOb,1) + if(iLUMOb <= nBas) LUMOb = eps(iLUMOb,2) HOMO = max(HOMOa,HOMOb) LUMO = min(LUMOa,LUMOb) diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_individual_energy.f90 similarity index 96% rename from src/eDFT/print_unrestricted_individual_energy.f90 rename to src/eDFT/print_individual_energy.f90 index 257652b..88ff78b 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_individual_energy.f90 @@ -1,4 +1,4 @@ -subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, & +subroutine print_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux,LZH,LZx,LZc,ExDD,EcDD,E, & Om,OmH,Omx,Omc,Omaux,OmxDD,OmcDD) ! Print individual energies for eDFT calculation @@ -149,12 +149,12 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:)),' au' write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:)),' au' write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:)),' au' - write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZx(:)),' au' + write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',sum(LZH(:))+sum(LZx(:))+sum(LZc(:)),' au' write(*,*) write(*,'(A40,F16.10,A3)') ' H Levy-Zahariev shift: ',sum(LZH(:))*HaToeV,' eV' write(*,'(A40,F16.10,A3)') ' x Levy-Zahariev shift: ',sum(LZx(:))*HaToeV,' eV' write(*,'(A40,F16.10,A3)') ' c Levy-Zahariev shift: ',sum(LZc(:))*HaToeV,' eV' - write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZx(:)))*HaToeV,' eV' + write(*,'(A40,F16.10,A3)') ' Hxc Levy-Zahariev shift: ',(sum(LZH(:))+sum(LZx(:))+sum(LZc(:)))*HaToeV,' eV' write(*,'(A60)') '-------------------------------------------------' write(*,*) @@ -214,7 +214,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux write(*,'(A60)') '-------------------------------------------------' do iEns=2,nEns -! write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns), ' au' + write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns), ' au' write(*,*) write(*,'(A44, F16.10,A3)') ' H energy contribution : ',OmH(iEns), ' au' write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns), ' au' @@ -228,7 +228,7 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns)*HaToeV, ' eV' + write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Om(iEns)*HaToeV, ' eV' write(*,*) write(*,'(A44, F16.10,A3)') ' H energy contribution : ',OmH(iEns)*HaToeV, ' eV' write(*,'(A44, F16.10,A3)') ' x energy contribution : ',Omx(iEns)*HaToeV, ' eV' @@ -243,4 +243,4 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EH,Ex,Ec,Eaux write(*,'(A60)') '-------------------------------------------------' write(*,*) -end subroutine print_unrestricted_individual_energy +end subroutine print_individual_energy diff --git a/src/eDFT/read_options_dft.f90 b/src/eDFT/read_options_dft.f90 index d718477..ea972a5 100644 --- a/src/eDFT/read_options_dft.f90 +++ b/src/eDFT/read_options_dft.f90 @@ -117,7 +117,11 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, case ('PBE') x_DFA = 3 - + + case ('CC-B88') + + x_DFA = 4 + case default call print_warning('!!! GGA exchange functional not available !!!') diff --git a/src/eDFT/xc_potential_grid.f90 b/src/eDFT/xc_potential_grid.f90 new file mode 100644 index 0000000..3a29628 --- /dev/null +++ b/src/eDFT/xc_potential_grid.f90 @@ -0,0 +1,54 @@ +subroutine xc_potential_grid(nBas,nGrid,AO,rho,Fx,Vxgrid) + + +! Compute the exchange-correlation potential on the grid + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: nBas, nGrid + double precision,intent(in) :: rho(nGrid,nspin) + double precision,intent(in) :: Fx(nBas,nBas,nspin) + double precision,intent(in) :: AO(nBas,nGrid) + +! Local variables + + integer :: mu,nu + integer :: ispin,iG + double precision :: r + double precision :: Fxgrid(nGrid,nspin) + +! Output variables + + double precision,intent(out) :: Vxgrid(nGrid) + +! Compute Vx + + Vxgrid(:) = 0d0 + Fxgrid(:,:) = 0d0 + + do iG=1,nGrid + do ispin=1,nspin + do mu=1,nBas + do nu=1,nBas + r = max(0d0,rho(iG,ispin)) + if(r > threshold) then + Fxgrid(iG,ispin) = Fxgrid(iG,ispin) + AO(mu,iG)*AO(nu,iG)*4d0/3d0*CxLSDA*r**(1d0/3d0) + endif + enddo + enddo + enddo + enddo + + Vxgrid(:)=Fxgrid(:,1)+Fxgrid(:,2) + open(411, file = 'Vxgrid.dat', status = 'new') + do iG=1,nGrid + write(411,*) iG, Vxgrid(iG) + end do + close(411) + + +end subroutine xc_potential_grid +