From ed05965399c700a9b6f7d09c7ae8fa6d75efb128 Mon Sep 17 00:00:00 2001 From: Pierre-Francois Loos Date: Mon, 29 Nov 2021 13:11:01 +0100 Subject: [PATCH] commit b4 indiv energy --- input/dft | 20 ++++---- src/eDFT/UCC_lda_exchange_potential.f90 | 2 + ...VWN5_lda_correlation_individual_energy.f90 | 4 +- .../print_unrestricted_individual_energy.f90 | 49 +++++++++---------- src/eDFT/read_options_dft.f90 | 15 +++--- src/eDFT/unrestricted_auxiliary_energy.f90 | 5 ++ 6 files changed, 51 insertions(+), 44 deletions(-) diff --git a/input/dft b/input/dft index 7a4ae5a..7522ef8 100644 --- a/input/dft +++ b/input/dft @@ -13,17 +13,17 @@ # GGA = 2: LYP,PBE # MGGA = 3: # Hybrid = 4: HF,B3LYP,PBE -0 H +1 VWN5 # quadrature grid SG-n - 1 +0 # Number of states in ensemble (nEns) -4 +2 # occupation numbers -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 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 -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 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 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 +1 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 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 @@ -31,13 +31,13 @@ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 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 + 0.75 0.0 0.0 # Ncentered ? F # Parameters for CC weight-dependent exchange functional 4 -0.60601 -0.0631565 -0.0289751 0.00244785 --1.28839 -0.179261 0.105627 -0.0215269 +0.0 0.0 0.0 0.0 +0.0 0.0 0.0 0.0 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/src/eDFT/UCC_lda_exchange_potential.f90 b/src/eDFT/UCC_lda_exchange_potential.f90 index a7ca60e..98d3067 100644 --- a/src/eDFT/UCC_lda_exchange_potential.f90 +++ b/src/eDFT/UCC_lda_exchange_potential.f90 @@ -73,6 +73,8 @@ subroutine UCC_lda_exchange_potential(nEns,wEns,nCC,aCC,nGrid,weight,nBas,AO,rho w1 = wEns(2) Fx1 = 1d0 - w1*(1d0 - w1)*(a1 + b1*(w1 - 0.5d0) + c1*(w1 - 0.5d0)**2) + + w2 = wEns(3) Fx2 = 1d0 - w2*(1d0 - w2)*(a2 + b2*(w2 - 0.5d0) + c2*(w2 - 0.5d0)**2) endif diff --git a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 index f26936e..c8a0489 100644 --- a/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 +++ b/src/eDFT/UVWN5_lda_correlation_individual_energy.f90 @@ -175,9 +175,9 @@ subroutine UVWN5_lda_correlation_individual_energy(nGrid,weight,rhow,rho,doNcent decdrb = decdr_p + decdr_a*fz/d2fz*(1d0-z**4) + ec_a*dfzdrb/d2fz*(1d0-z**4) - 4d0*ec_a*fz/d2fz*dzdrb*z**3 & + (decdr_f - decdr_p)*fz*z**4 + (ec_f - ec_p)*dfzdrb*z**4 + 4d0*(ec_f - ec_p)*fz*dzdrb*z**3 - decdr = 0d0 + decdr = 0d0 if(ra > threshold) decdr = decdr + decdra - if(rb > threshold) decdr = decdr + decdrb +! if(rb > threshold) decdr = decdr + decdrb Ecrr(2) = Ecrr(2) - weight(iG)*decdr*r*r diff --git a/src/eDFT/print_unrestricted_individual_energy.f90 b/src/eDFT/print_unrestricted_individual_energy.f90 index f75eaa0..bd8e7d8 100644 --- a/src/eDFT/print_unrestricted_individual_energy.f90 +++ b/src/eDFT/print_unrestricted_individual_energy.f90 @@ -149,39 +149,38 @@ subroutine print_unrestricted_individual_energy(nEns,ENuc,Ew,ET,EV,EJ,Ex,Ec,Exc, ! Total Energy and IP and EA !------------------------------------------------------------------------ -! write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A60)') ' IP AND EA FROM AUXILIARY ENERGIES ' -! write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') '-------------------------------------------------' + write(*,'(A60)') ' ENERGY DIFFERENCES FROM AUXILIARY ENERGIES ' + write(*,'(A60)') '-------------------------------------------------' -! do iEns=2,nEns -! write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Omaux(iEns)+OmxcDD(iEns),' au' -! write(*,*) -! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au' -! 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(*,*) + do iEns=2,nEns + write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',Omaux(iEns)+OmxcDD(iEns),' au' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns), ' au' + 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(*,*) -! write(*,'(A60)') '-------------------------------------------------' -! write(*,*) + write(*,'(A60)') '-------------------------------------------------' -! write(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',(Omaux(iEns)+OmxcDD(iEns))*HaToeV,' eV' -! write(*,*) -! write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV' -! 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(*,'(A40,I2,A1,F16.10,A3)') ' Energy difference 1 -> ',iEns,':',(Omaux(iEns)+OmxcDD(iEns))*HaToeV,' eV' + write(*,*) + write(*,'(A44, F16.10,A3)') ' auxiliary energy contribution : ',Omaux(iEns)*HaToeV, ' eV' + 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(*,*) + write(*,'(A60)') '-------------------------------------------------' + write(*,*) write(*,'(A60)') '-------------------------------------------------' - write(*,'(A60)') ' IP and EA FROM INDIVIDUAL ENERGIES ' + write(*,'(A60)') ' ENERGY DIFFERENCES FROM INDIVIDUAL ENERGIES ' write(*,'(A60)') '-------------------------------------------------' do iEns=1,nEns -! write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' + write(*,'(A40,I2,A2,F16.10,A3)') ' Individual energy state ',iEns,': ',E(iEns) + ENuc,' au' end do write(*,'(A60)') '-------------------------------------------------' diff --git a/src/eDFT/read_options_dft.f90 b/src/eDFT/read_options_dft.f90 index a17b0ee..d718477 100644 --- a/src/eDFT/read_options_dft.f90 +++ b/src/eDFT/read_options_dft.f90 @@ -400,7 +400,7 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, ! Read parameters for weight-dependent functional read(1,*) read(1,*) nCC - do iEns=2,nEns + do iEns=2,maxEns read(1,*) (aCC(iCC,iEns-1),iCC=1,nCC) end do @@ -410,12 +410,13 @@ subroutine read_options_dft(nBas,method,x_rung,x_DFA,c_rung,c_DFA,SGn,nEns,wEns, write(*,*)'----------------------------------------------------------' write(*,*)' Parameters for weight-dependent exchange functional ' - write(*,*)'----------------------------------------------------------' - do iEns=2,nEns - write(*,'(A6,I2,A2)') 'State ',iEns,':' - do iCC=1,nCC - write(*,'(I2,F10.6)') iCC,aCC(iCC,iEns-1) - end do + do iEns=2,maxEns + write(*,*)'----------------------------------------------------------' + write(*,'(A10,I2,A2)') ' State ',iEns,':' + write(*,*)'----------------------------------------------------------' + write(*,*) + call matout(nCC,1,acc(:,iEns-1)) + write(*,*) end do write(*,*) diff --git a/src/eDFT/unrestricted_auxiliary_energy.f90 b/src/eDFT/unrestricted_auxiliary_energy.f90 index cb7c987..c609caf 100644 --- a/src/eDFT/unrestricted_auxiliary_energy.f90 +++ b/src/eDFT/unrestricted_auxiliary_energy.f90 @@ -26,7 +26,12 @@ subroutine unrestricted_auxiliary_energy(nBas,nEns,eps,occnum,doNcentered,Eaux) double precision,intent(out) :: Eaux(nspin,nEns) +! Memory allocation + allocate(nEl(nEns)) + +! Compute the number of electrons + nEl(:) = 0d0 do iEns=1,nEns do iBas=1,nBas