From 2ce6b128547c74bbd405e4ebea9816b5611ef104 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Sun, 26 Apr 2020 23:17:27 +0200 Subject: [PATCH] Three-state extension and MOM --- examples/molecule.H2 | 2 +- input/basis | 78 ++-------- input/dft | 8 +- input/molecule | 5 +- input/molecule.xyz | 5 +- input/weight | 78 ++-------- src/eDFT/LIM_RKS.f90 | 66 +++++--- src/eDFT/MOM_RKS.f90 | 142 ++++++++++++++++++ ..._lda_exchange_derivative_discontinuity.f90 | 17 ++- src/eDFT/RGIC_lda_exchange_energy.f90 | 14 +- .../RGIC_lda_exchange_individual_energy.f90 | 14 +- src/eDFT/RGIC_lda_exchange_potential.f90 | 14 +- ...a_correlation_derivative_discontinuity.f90 | 10 +- src/eDFT/RMFL20_lda_correlation_energy.f90 | 11 +- ...FL20_lda_correlation_individual_energy.f90 | 11 +- src/eDFT/RMFL20_lda_correlation_potential.f90 | 11 +- src/eDFT/eDFT.f90 | 22 ++- .../print_restricted_individual_energy.f90 | 5 + src/eDFT/restricted_auxiliary_energy.f90 | 14 +- src/eDFT/restricted_density_matrix.f90 | 13 ++ 20 files changed, 344 insertions(+), 196 deletions(-) create mode 100644 src/eDFT/MOM_RKS.f90 diff --git a/examples/molecule.H2 b/examples/molecule.H2 index 02328b1..81c624a 100644 --- a/examples/molecule.H2 +++ b/examples/molecule.H2 @@ -2,4 +2,4 @@ 2 1 1 0 0 # Znuc x y z H 0. 0. 0. - H 0. 0. 2.5 + H 0. 0. 1.4 diff --git a/input/basis b/input/basis index 27669ff..fb05e68 100644 --- a/input/basis +++ b/input/basis @@ -1,66 +1,18 @@ -1 21 -S 10 - 1 54620.0000000 0.0000180 - 2 8180.0000000 0.0001380 - 3 1862.0000000 0.0007230 - 4 527.3000000 0.0030390 - 5 172.0000000 0.0109080 - 6 62.1000000 0.0340350 - 7 24.2100000 0.0911930 - 8 9.9930000 0.1992680 - 9 4.3050000 0.3293550 - 10 1.9210000 0.3404890 -S 10 - 1 54620.0000000 -0.0000030 - 2 8180.0000000 -0.0000250 - 3 1862.0000000 -0.0001310 - 4 527.3000000 -0.0005580 - 5 172.0000000 -0.0019880 - 6 62.1000000 -0.0063700 - 7 24.2100000 -0.0172170 - 8 9.9930000 -0.0408580 - 9 4.3050000 -0.0742370 - 10 1.9210000 -0.1192340 +1 3 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.8663000 1.0000000 + 1 0.1220000 1.0000000 +P 1 + 1 0.7270000 1.0000000 +2 3 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.2475000 1.0000000 -S 1 - 1 0.1009000 1.0000000 -S 1 - 1 0.0412900 1.0000000 -P 4 - 1 43.7500000 0.0006330 - 2 10.3300000 0.0048080 - 3 3.2260000 0.0205270 - 4 1.1270000 0.0678160 + 1 0.1220000 1.0000000 P 1 - 1 0.4334000 1.0000000 -P 1 - 1 0.1808000 1.0000000 -P 1 - 1 0.0782700 1.0000000 -P 1 - 1 0.0337200 1.0000000 -D 1 - 1 1.6350000 1.0000000 -D 1 - 1 0.7410000 1.0000000 -D 1 - 1 0.3350000 1.0000000 -D 1 - 1 0.1519000 1.0000000 -F 1 - 1 0.6860000 1.0000000 -F 1 - 1 0.4010000 1.0000000 -F 1 - 1 0.2350000 1.0000000 -G 1 - 1 0.6030000 1.0000000 -G 1 - 1 0.3240000 1.0000000 -H 1 - 1 0.5100000 1.0000000 - - + 1 0.7270000 1.0000000 diff --git a/input/dft b/input/dft index 1c8aaa6..3dda946 100644 --- a/input/dft +++ b/input/dft @@ -1,5 +1,5 @@ # Restricted or unrestricted KS calculation - GOK-RKS + MOM-RKS # exchange rung: # Hartree = 0 # LDA = 1: RS51,RMFL20 @@ -13,12 +13,12 @@ # GGA = 2: # Hybrid = 4: # Hartree-Fock = 666 - 1 RMFL20 + 1 RVWN5 # quadrature grid SG-n 1 # Number of states in ensemble (nEns) - 2 + 3 # Ensemble weights: wEns(1),...,wEns(nEns-1) - 0.0 + 0.0000000 0.0000000 # GOK-DFT: maxSCF thresh DIIS n_diis guess_type ortho_type 32 0.00001 T 5 1 1 diff --git a/input/molecule b/input/molecule index 6a6f6d1..81c624a 100644 --- a/input/molecule +++ b/input/molecule @@ -1,4 +1,5 @@ # nAt nEla nElb nCore nRyd - 1 2 2 0 0 + 2 1 1 0 0 # Znuc x y z - Be 0.0 0.0 0.0 + H 0. 0. 0. + H 0. 0. 1.4 diff --git a/input/molecule.xyz b/input/molecule.xyz index 8023e37..6edc99d 100644 --- a/input/molecule.xyz +++ b/input/molecule.xyz @@ -1,3 +1,4 @@ - 1 + 2 - Be 0.0000000000 0.0000000000 0.0000000000 + H 0.0000000000 0.0000000000 0.0000000000 + H 0.0000000000 0.0000000000 0.7408481486 diff --git a/input/weight b/input/weight index 27669ff..fb05e68 100644 --- a/input/weight +++ b/input/weight @@ -1,66 +1,18 @@ -1 21 -S 10 - 1 54620.0000000 0.0000180 - 2 8180.0000000 0.0001380 - 3 1862.0000000 0.0007230 - 4 527.3000000 0.0030390 - 5 172.0000000 0.0109080 - 6 62.1000000 0.0340350 - 7 24.2100000 0.0911930 - 8 9.9930000 0.1992680 - 9 4.3050000 0.3293550 - 10 1.9210000 0.3404890 -S 10 - 1 54620.0000000 -0.0000030 - 2 8180.0000000 -0.0000250 - 3 1862.0000000 -0.0001310 - 4 527.3000000 -0.0005580 - 5 172.0000000 -0.0019880 - 6 62.1000000 -0.0063700 - 7 24.2100000 -0.0172170 - 8 9.9930000 -0.0408580 - 9 4.3050000 -0.0742370 - 10 1.9210000 -0.1192340 +1 3 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.8663000 1.0000000 + 1 0.1220000 1.0000000 +P 1 + 1 0.7270000 1.0000000 +2 3 +S 3 + 1 13.0100000 0.0196850 + 2 1.9620000 0.1379770 + 3 0.4446000 0.4781480 S 1 - 1 0.2475000 1.0000000 -S 1 - 1 0.1009000 1.0000000 -S 1 - 1 0.0412900 1.0000000 -P 4 - 1 43.7500000 0.0006330 - 2 10.3300000 0.0048080 - 3 3.2260000 0.0205270 - 4 1.1270000 0.0678160 + 1 0.1220000 1.0000000 P 1 - 1 0.4334000 1.0000000 -P 1 - 1 0.1808000 1.0000000 -P 1 - 1 0.0782700 1.0000000 -P 1 - 1 0.0337200 1.0000000 -D 1 - 1 1.6350000 1.0000000 -D 1 - 1 0.7410000 1.0000000 -D 1 - 1 0.3350000 1.0000000 -D 1 - 1 0.1519000 1.0000000 -F 1 - 1 0.6860000 1.0000000 -F 1 - 1 0.4010000 1.0000000 -F 1 - 1 0.2350000 1.0000000 -G 1 - 1 0.6030000 1.0000000 -G 1 - 1 0.3240000 1.0000000 -H 1 - 1 0.5100000 1.0000000 - - + 1 0.7270000 1.0000000 diff --git a/src/eDFT/LIM_RKS.f90 b/src/eDFT/LIM_RKS.f90 index 7a9f80a..e81a09f 100644 --- a/src/eDFT/LIM_RKS.f90 +++ b/src/eDFT/LIM_RKS.f90 @@ -1,4 +1,4 @@ -subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight, & +subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c) ! Perform restricted Kohn-Sham calculation for ensembles @@ -12,7 +12,6 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight character(len=12),intent(in) :: x_DFA,c_DFA logical,intent(in) :: LDA_centered integer,intent(in) :: nEns - double precision,intent(in) :: wEns(nEns) integer,intent(in) :: nGrid double precision,intent(in) :: weight(nGrid) integer,intent(in) :: maxSCF,max_diis,guess_type @@ -35,8 +34,7 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight ! Local variables integer :: iEns - double precision :: EwZW - double precision :: EwEW + double precision :: Ew(nEnS) double precision :: wLIM(nEns) double precision :: Om(nEns) @@ -49,6 +47,11 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight write(*,*)'************************************************' write(*,*) +! Initializatio + + Ew(:) = 0d0 + Om(:) = 0d0 + !------------------------------------------------------------------------ ! Zero-weight calculation !------------------------------------------------------------------------ @@ -57,8 +60,9 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight write(*,'(A40)') ' ZERO-WEIGHT CALCULATION ' write(*,'(A40)') '*************************************************' - wLIM(1) = 1d0 - wLIM(2:nEns) = 0d0 + wLIM(1) = 1d0 + wLIM(2) = 0d0 + wLIM(3) = 0d0 do iEns=1,nEns write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns) @@ -67,43 +71,71 @@ subroutine LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight write(*,*) call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwZW,c) + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c) !------------------------------------------------------------------------ ! Equiensemble calculation !------------------------------------------------------------------------ write(*,'(A40)') '*************************************************' - write(*,'(A40)') ' NON-ZERO-WEIGHT CALCULATION ' + write(*,'(A40)') ' TWO-STATE EQUI-WEIGHT CALCULATION ' write(*,'(A40)') '*************************************************' + wLIM(1) = 0.5d0 + wLIM(2) = 0.5d0 + wLIM(3) = 0.0d0 + do iEns=1,nEns - write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wEns(iEns) + write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns) end do write(*,'(A40)') '*************************************************' write(*,*) - call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns,nGrid,weight,maxSCF,thresh, & - max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,EwEW,c) + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c) + +!------------------------------------------------------------------------ +! Equiensemble calculation +!------------------------------------------------------------------------ + + write(*,'(A40)') '*************************************************' + write(*,'(A40)') ' THREE-STATE EQUI-WEIGHT CALCULATION ' + write(*,'(A40)') '*************************************************' + + wLIM(1) = 1d0/3d0 + wLIM(2) = 1d0/3d0 + wLIM(3) = 1d0/3d0 + + do iEns=1,nEns + write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wLIM(iEns) + end do + write(*,'(A40)') '*************************************************' + write(*,*) + + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wLIM,nGrid,weight,maxSCF,thresh, & + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c) + !------------------------------------------------------------------------ ! LIM excitation energies !------------------------------------------------------------------------ - Om(:) = 0d0 - if(wEns(2) > 10d-3) then - Om(2) = (EwEW - EwZW)/wEns(2) - end if + Om(2) = (Ew(2) - Ew(1))/2d0 + Om(3) = (Ew(3) - Ew(1))/3d0 + 0.5d0*Om(2) write(*,'(A60)') '-------------------------------------------------' write(*,'(A60)') ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES ' write(*,'(A60)') '-------------------------------------------------' - write(*,'(A44,F16.10,A3)') ' Zero-weight ensemble energy',EwZW, ' au' - write(*,'(A44,F16.10,A3)') ' Equi-weight ensemble energy',EwEW, ' au' + write(*,'(A44,F16.10,A3)') ' Ensemble energy #1 ',Ew(1),' au' + write(*,'(A44,F16.10,A3)') ' Ensemble energy #2 ',Ew(2),' au' + write(*,'(A44,F16.10,A3)') ' Ensemble energy #3 ',Ew(3),' au' write(*,'(A60)') '-------------------------------------------------' do iEns=2,nEns write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au' + end do + write(*,*) + do iEns=2,nEns write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV' end do write(*,'(A60)') '-------------------------------------------------' diff --git a/src/eDFT/MOM_RKS.f90 b/src/eDFT/MOM_RKS.f90 new file mode 100644 index 0000000..5f41a01 --- /dev/null +++ b/src/eDFT/MOM_RKS.f90 @@ -0,0 +1,142 @@ +subroutine MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight, & + maxSCF,thresh,max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,c) + +! Perform restricted Kohn-Sham calculation for ensembles + + implicit none + include 'parameters.h' + +! Input variables + + integer,intent(in) :: x_rung,c_rung + character(len=12),intent(in) :: x_DFA,c_DFA + logical,intent(in) :: LDA_centered + integer,intent(in) :: nEns + integer,intent(in) :: nGrid + double precision,intent(in) :: weight(nGrid) + integer,intent(in) :: maxSCF,max_diis,guess_type + double precision,intent(in) :: thresh + integer,intent(in) :: nBas + double precision,intent(in) :: AO(nBas,nGrid) + double precision,intent(in) :: dAO(ncart,nBas,nGrid) + + integer,intent(in) :: nO,nV + double precision,intent(in) :: S(nBas,nBas) + double precision,intent(in) :: T(nBas,nBas) + double precision,intent(in) :: V(nBas,nBas) + double precision,intent(in) :: Hc(nBas,nBas) + double precision,intent(in) :: X(nBas,nBas) + double precision,intent(in) :: ERI(nBas,nBas,nBas,nBas) + double precision,intent(in) :: ENuc + + double precision,intent(out) :: c(nBas,nBas) + +! Local variables + + integer :: iEns + double precision :: Ew(nEns) + double precision :: wMOM(nEns) + double precision :: Om(nEns) + +! Hello world + + write(*,*) + write(*,*)'************************************************' + write(*,*)'* Maximum Overlao method *' + write(*,*)'* for excitation energies *' + write(*,*)'************************************************' + write(*,*) + +! Initialization + + Ew(:) = 0d0 + Om(:) = 0d0 + +!------------------------------------------------------------------------ +! Zero-weight calculation +!------------------------------------------------------------------------ + + write(*,'(A40)') '*************************************************' + write(*,'(A40)') ' STATE-SPECIFIC CALCULATION #1 ' + write(*,'(A40)') '*************************************************' + + wMOM(1) = 1d0 + wMOM(2) = 0d0 + wMOM(3) = 0d0 + + do iEns=1,nEns + write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wMOM(iEns) + end do + write(*,'(A40)') '*************************************************' + write(*,*) + + call GOK_RKS(.false.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(1),c) + +!------------------------------------------------------------------------ +! Equiensemble calculation +!------------------------------------------------------------------------ + + write(*,'(A40)') '*************************************************' + write(*,'(A40)') ' STATE-SPECIFIC CALCULATION #2 ' + write(*,'(A40)') '*************************************************' + + wMOM(1) = 0d0 + wMOM(2) = 1d0 + wMOM(3) = 0d0 + + do iEns=1,nEns + write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wMOM(iEns) + end do + write(*,'(A40)') '*************************************************' + write(*,*) + + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(2),c) + +!------------------------------------------------------------------------ +! Equiensemble calculation +!------------------------------------------------------------------------ + + write(*,'(A40)') '*************************************************' + write(*,'(A40)') ' STATE-SPECIFIC CALCULATION #3 ' + write(*,'(A40)') '*************************************************' + + wMOM(1) = 0d0 + wMOM(2) = 0d0 + wMOM(3) = 1.0d0 + + do iEns=1,nEns + write(*,'(A20,I2,A2,F16.10)') ' Weight of state ',iEns,': ',wMOM(iEns) + end do + write(*,'(A40)') '*************************************************' + write(*,*) + + call GOK_RKS(.true.,x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wMOM,nGrid,weight,maxSCF,thresh, & + max_diis,guess_type,nBas,AO,dAO,nO,nV,S,T,V,Hc,ERI,X,ENuc,Ew(3),c) + +!------------------------------------------------------------------------ +! MOM excitation energies +!------------------------------------------------------------------------ + + Om(:) = Ew(:) - Ew(1) + + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' LINEAR INTERPOLATION METHOD EXCITATION ENERGIES ' + write(*,'(A60)') '-------------------------------------------------' + + write(*,'(A44,F16.10,A3)') ' Ensemble energy #1 ',Ew(1),' au' + write(*,'(A44,F16.10,A3)') ' Ensemble energy #2 ',Ew(2),' au' + write(*,'(A44,F16.10,A3)') ' Ensemble energy #3 ',Ew(3),' au' + write(*,'(A60)') '-------------------------------------------------' + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns), ' au' + end do + write(*,*) + do iEns=2,nEns + write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',Om(iEns)*HaToeV,' eV' + end do + write(*,'(A60)') '-------------------------------------------------' + + +end subroutine MOM_RKS diff --git a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 index 479879a..11b2f4a 100644 --- a/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 +++ b/src/eDFT/RGIC_lda_exchange_derivative_discontinuity.f90 @@ -37,9 +37,9 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho ! Parameters for H2 at equilibrium -! a = + 0.5751782560799208d0 -! b = - 0.021108186591137282d0 -! c = - 0.36718902716347124d0 + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 ! Parameters for stretch H2 @@ -55,11 +55,11 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho ! Parameters for HNO - a = 0.0061158387543040335d0 - b = -0.00005968703047293955d0 - c = -0.00001692245714408755d0 +! a = 0.0061158387543040335d0 +! b = -0.00005968703047293955d0 +! c = -0.00001692245714408755d0 - w = wEns(2) + w = 0.5d0*wEns(2) + wEns(3) dCxGICdw = (0.5d0*b + (2d0*a + 0.5d0*c)*(w - 0.5d0) - (1d0 - w)*w*(3d0*b + 4d0*c*(w - 0.5d0))) dCxGICdw = CxLDA*dCxGICdw @@ -72,7 +72,8 @@ subroutine RGIC_lda_exchange_derivative_discontinuity(nEns,wEns,nGrid,weight,rho if(r > threshold) then dExdw(1) = 0d0 - dExdw(2) = dExdw(2) + weight(iG)*dCxGICdw*r**(4d0/3d0) + dExdw(2) = dExdw(2) + 0.5d0*weight(iG)*dCxGICdw*r**(4d0/3d0) + dExdw(3) = dExdw(3) + 1.0d0*weight(iG)*dCxGICdw*r**(4d0/3d0) end if diff --git a/src/eDFT/RGIC_lda_exchange_energy.f90 b/src/eDFT/RGIC_lda_exchange_energy.f90 index bee986b..901210a 100644 --- a/src/eDFT/RGIC_lda_exchange_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_energy.f90 @@ -29,9 +29,9 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Parameters for H2 at equilibrium -! a = + 0.5751782560799208d0 -! b = - 0.021108186591137282d0 -! c = - 0.36718902716347124d0 + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 ! Parameters for stretch H2 @@ -47,11 +47,11 @@ subroutine RGIC_lda_exchange_energy(nEns,wEns,nGrid,weight,rho,Ex) ! Parameters for HNO - a = 0.0061158387543040335d0 - b = -0.00005968703047293955d0 - c = -0.00001692245714408755d0 +! a = 0.0061158387543040335d0 +! b = -0.00005968703047293955d0 +! c = -0.00001692245714408755d0 - w = wEns(2) + w = 0.5d0*wEns(2) + wEns(3) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) CxGIC = CxLDA*CxGIC diff --git a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 index ecb5849..3788c6d 100644 --- a/src/eDFT/RGIC_lda_exchange_individual_energy.f90 +++ b/src/eDFT/RGIC_lda_exchange_individual_energy.f90 @@ -32,9 +32,9 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E ! Parameters for H2 at equilibrium -! a = + 0.5751782560799208d0 -! b = - 0.021108186591137282d0 -! c = - 0.36718902716347124d0 + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 ! Parameters for stretch H2 @@ -51,11 +51,11 @@ subroutine RGIC_lda_exchange_individual_energy(nEns,wEns,nGrid,weight,rhow,rho,E ! Parameters for HNO - a = 0.0061158387543040335d0 - b = -0.00005968703047293955d0 - c = -0.00001692245714408755d0 +! a = 0.0061158387543040335d0 +! b = -0.00005968703047293955d0 +! c = -0.00001692245714408755d0 - w = wEns(2) + w = 0.5d0*wEns(2) + wEns(3) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) CxGIC = CxLDA*CxGIC diff --git a/src/eDFT/RGIC_lda_exchange_potential.f90 b/src/eDFT/RGIC_lda_exchange_potential.f90 index b200c05..d27ec2d 100644 --- a/src/eDFT/RGIC_lda_exchange_potential.f90 +++ b/src/eDFT/RGIC_lda_exchange_potential.f90 @@ -32,9 +32,9 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Parameters for H2 at equilibrium -! a = + 0.5751782560799208d0 -! b = - 0.021108186591137282d0 -! c = - 0.36718902716347124d0 + a = + 0.5751782560799208d0 + b = - 0.021108186591137282d0 + c = - 0.36718902716347124d0 ! Parameters for stretch H2 @@ -51,11 +51,11 @@ subroutine RGIC_lda_exchange_potential(nEns,wEns,nGrid,weight,nBas,AO,rho,Fx) ! Parameters for HNO - a = 0.0061158387543040335d0 - b = -0.00005968703047293955d0 - c = -0.00001692245714408755d0 +! a = 0.0061158387543040335d0 +! b = -0.00005968703047293955d0 +! c = -0.00001692245714408755d0 - w = wEns(2) + w = 0.5d0*wEns(2) + wEns(3) CxGIC = 1d0 - w*(1d0 - w)*(a + b*(w - 0.5d0) + c*(w - 0.5d0)**2) CxGIC = CxLDA*CxGIC diff --git a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 index d6f7487..5f3e4a3 100644 --- a/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 +++ b/src/eDFT/RMFL20_lda_correlation_derivative_discontinuity.f90 @@ -34,9 +34,13 @@ subroutine RMFL20_lda_correlation_derivative_discontinuity(nEns,wEns,nGrid,weigh aMFL(2,1) = +0.00540994d0 aMFL(3,1) = +0.0830766d0 - aMFL(1,2) = -0.0144633d0 - aMFL(2,2) = -0.0506019d0 - aMFL(3,2) = +0.0331417d0 + aMFL(1,2) = -0.0282814d0 + aMFL(2,2) = +0.00273925d0 + aMFL(3,2) = +0.0664914d0 + + aMFL(1,3) = -0.0144633d0 + aMFL(2,3) = -0.0506019d0 + aMFL(3,3) = +0.0331417d0 ! Compute correlation energy for ground, singly-excited and doubly-excited states diff --git a/src/eDFT/RMFL20_lda_correlation_energy.f90 b/src/eDFT/RMFL20_lda_correlation_energy.f90 index 8826761..731378f 100644 --- a/src/eDFT/RMFL20_lda_correlation_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_energy.f90 @@ -36,9 +36,13 @@ subroutine RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns,nGrid,weight,rho aMFL(2,1) = +0.00540994d0 aMFL(3,1) = +0.0830766d0 - aMFL(1,2) = -0.0144633d0 - aMFL(2,2) = -0.0506019d0 - aMFL(3,2) = +0.0331417d0 + aMFL(1,2) = -0.0282814d0 + aMFL(2,2) = +0.00273925d0 + aMFL(3,2) = +0.0664914d0 + + aMFL(1,3) = -0.0144633d0 + aMFL(2,3) = -0.0506019d0 + aMFL(3,3) = +0.0331417d0 ! Compute correlation energy for ground and doubly-excited states @@ -52,6 +56,7 @@ subroutine RMFL20_lda_correlation_energy(LDA_centered,nEns,wEns,nGrid,weight,rho call RVWN5_lda_correlation_energy(nGrid,weight(:),rho(:),EcLDA) + EceLDA(3) = EcLDA + wEns(3)*(EceLDA(3) - EceLDA(1)) EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1)) EceLDA(1) = EcLDA diff --git a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 index f405063..93a6170 100644 --- a/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 +++ b/src/eDFT/RMFL20_lda_correlation_individual_energy.f90 @@ -36,9 +36,13 @@ subroutine RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid aMFL(2,1) = +0.00540994d0 aMFL(3,1) = +0.0830766d0 - aMFL(1,2) = -0.0144633d0 - aMFL(2,2) = -0.0506019d0 - aMFL(3,2) = +0.0331417d0 + aMFL(1,2) = -0.0282814d0 + aMFL(2,2) = +0.00273925d0 + aMFL(3,2) = +0.0664914d0 + + aMFL(1,3) = -0.0144633d0 + aMFL(2,3) = -0.0506019d0 + aMFL(3,3) = +0.0331417d0 ! Compute correlation energy for ground- and doubly-excited states @@ -52,6 +56,7 @@ subroutine RMFL20_lda_correlation_individual_energy(LDA_centered,nEns,wEns,nGrid call RVWN5_lda_correlation_individual_energy(nGrid,weight(:),rhow(:),rho(:),EcLDA) + EceLDA(3) = EcLDA + wEns(3)*(EceLDA(3) - EceLDA(1)) EceLDA(2) = EcLDA + wEns(2)*(EceLDA(2) - EceLDA(1)) EceLDA(1) = EcLDA diff --git a/src/eDFT/RMFL20_lda_correlation_potential.f90 b/src/eDFT/RMFL20_lda_correlation_potential.f90 index 0007e35..e1ffffc 100644 --- a/src/eDFT/RMFL20_lda_correlation_potential.f90 +++ b/src/eDFT/RMFL20_lda_correlation_potential.f90 @@ -37,9 +37,13 @@ subroutine RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns,nGrid,weight, aMFL(2,1) = +0.00540994d0 aMFL(3,1) = +0.0830766d0 - aMFL(1,2) = -0.0144633d0 - aMFL(2,2) = -0.0506019d0 - aMFL(3,2) = +0.0331417d0 + aMFL(1,2) = -0.0282814d0 + aMFL(2,2) = +0.00273925d0 + aMFL(3,2) = +0.0664914d0 + + aMFL(1,3) = -0.0144633d0 + aMFL(2,3) = -0.0506019d0 + aMFL(3,3) = +0.0331417d0 ! Compute correlation energy for ground- and doubly-excited states @@ -53,6 +57,7 @@ subroutine RMFL20_lda_correlation_potential(LDA_centered,nEns,wEns,nGrid,weight, call RVWN5_lda_correlation_potential(nGrid,weight(:),nBas,AO(:,:),rho(:),FcLDA(:,:)) + FceLDA(:,:,3) = FcLDA(:,:) + wEns(3)*(FceLDA(:,:,3) - FceLDA(:,:,1)) FceLDA(:,:,2) = FcLDA(:,:) + wEns(2)*(FceLDA(:,:,2) - FceLDA(:,:,1)) FceLDA(:,:,1) = FcLDA(:,:) diff --git a/src/eDFT/eDFT.f90 b/src/eDFT/eDFT.f90 index 4a6aec6..f8240a5 100644 --- a/src/eDFT/eDFT.f90 +++ b/src/eDFT/eDFT.f90 @@ -174,13 +174,13 @@ program eDFT end if !------------------------------------------------------------------------ -! Compute RKS energy +! Compute LIM excitation energies !------------------------------------------------------------------------ if(method == 'LIM-RKS') then call cpu_time(start_KS) - call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,wEns(:),nGrid,weight(:), & + call LIM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:)) call cpu_time(end_KS) @@ -191,6 +191,24 @@ program eDFT end if +!------------------------------------------------------------------------ +! Compute MOM excitation energies +!------------------------------------------------------------------------ + + if(method == 'MOM-RKS') then + + call cpu_time(start_KS) + call MOM_RKS(x_rung,x_DFA,c_rung,c_DFA,LDA_centered,nEns,nGrid,weight(:), & + maxSCF,thresh,max_diis,guess_type,nBas,AO(:,:),dAO(:,:,:),nO(1),nV(1), & + S(:,:),T(:,:),V(:,:),Hc(:,:),ERI(:,:,:,:),X(:,:),ENuc,c(:,:)) + call cpu_time(end_KS) + + t_KS = end_KS - start_KS + write(*,'(A65,1X,F9.3,A8)') 'Total CPU time for MOM-RKS = ',t_KS,' seconds' + write(*,*) + + end if + !------------------------------------------------------------------------ ! Compute GOK-UKS energy (BROKEN) !------------------------------------------------------------------------ diff --git a/src/eDFT/print_restricted_individual_energy.f90 b/src/eDFT/print_restricted_individual_energy.f90 index 295fde8..c67fea7 100644 --- a/src/eDFT/print_restricted_individual_energy.f90 +++ b/src/eDFT/print_restricted_individual_energy.f90 @@ -147,8 +147,10 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au' write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au' + write(*,*) end do write(*,'(A60)') '-------------------------------------------------' + write(*,*) do iEns=2,nEns write(*,'(A40,I2,A2,F16.10,A3)') ' Excitation energy 1 ->',iEns,': ',(Omaux(iEns)+OmxcDD(iEns))*HaToeV,' eV' @@ -157,6 +159,7 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV' write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' + write(*,*) end do write(*,'(A60)') '-------------------------------------------------' write(*,*) @@ -179,6 +182,7 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns), ' au' write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns), ' au' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns),' au' + write(*,*) end do write(*,'(A60)') '-------------------------------------------------' @@ -192,6 +196,7 @@ subroutine print_restricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc,Ea write(*,'(A44, F16.10,A3)') ' x ensemble derivative : ',OmxDD(iEns)*HaToeV, ' eV' write(*,'(A44, F16.10,A3)') ' c ensemble derivative : ',OmcDD(iEns)*HaToeV, ' eV' write(*,'(A44, F16.10,A3)') ' xc ensemble derivative : ',OmxcDD(iEns)*HaToeV,' eV' + write(*,*) end do write(*,'(A60)') '-------------------------------------------------' write(*,*) diff --git a/src/eDFT/restricted_auxiliary_energy.f90 b/src/eDFT/restricted_auxiliary_energy.f90 index 7e3f251..71b1c8f 100644 --- a/src/eDFT/restricted_auxiliary_energy.f90 +++ b/src/eDFT/restricted_auxiliary_energy.f90 @@ -27,7 +27,7 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux) Eaux(iEns) = 2d0*sum(eps(1:nO)) -! Doubly-excited state density matrix +! Singly-excited state density matrix iEns = 2 @@ -37,6 +37,18 @@ subroutine restricted_auxiliary_energy(nBas,nEns,nO,eps,Eaux) Eaux(iEns) = 0d0 end if + Eaux(iEns) = Eaux(iEns) + eps(nO) + eps(nO+1) + +! Doubly-excited state density matrix + + iEns = 3 + + if(nO > 1) then + Eaux(iEns) = 2d0*sum(eps(1:nO-1)) + else + Eaux(iEns) = 0d0 + end if + Eaux(iEns) = Eaux(iEns) + 2d0*eps(nO+1) end subroutine restricted_auxiliary_energy diff --git a/src/eDFT/restricted_density_matrix.f90 b/src/eDFT/restricted_density_matrix.f90 index 2fd36a8..907cc31 100644 --- a/src/eDFT/restricted_density_matrix.f90 +++ b/src/eDFT/restricted_density_matrix.f90 @@ -37,6 +37,19 @@ subroutine restricted_density_matrix(nBas,nEns,nO,c,P) P(:,:,iEns) = 0d0 end if + P(:,:,iEns) = P(:,:,iEns) + 1d0*matmul(c(:,nO :nO ),transpose(c(:,nO :nO ))) & + + 1d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1))) + +! Doubly-excited state density matrix + + iEns = 3 + + if(nO > 1) then + P(:,:,iEns) = 2d0*matmul(c(:,1:nO-1),transpose(c(:,1:nO-1))) + else + P(:,:,iEns) = 0d0 + end if + P(:,:,iEns) = P(:,:,iEns) + 2d0*matmul(c(:,nO+1:nO+1),transpose(c(:,nO+1:nO+1))) end subroutine restricted_density_matrix