From 372ccd2e4ea31f8cee4fa09e2c88ab58f98acbbd Mon Sep 17 00:00:00 2001 From: Clotilde Marut Date: Sun, 5 Jul 2020 19:04:08 +0200 Subject: [PATCH] still looking for some bug --- input/basis | 30 +++++++++--- input/weight | 30 +++++++++--- src/eDFT/eDFT_UKS.f90 | 3 +- src/eDFT/unrestricted_density_matrix.f90 | 7 ++- src/eDFT/unrestricted_individual_energy.f90 | 52 ++++++++++++++++++++- 5 files changed, 103 insertions(+), 19 deletions(-) diff --git a/input/basis b/input/basis index 6796e3b..cd6cb22 100644 --- a/input/basis +++ b/input/basis @@ -1,9 +1,25 @@ -1 3 -S 3 - 1 38.3600000 0.0238090 - 2 5.7700000 0.1548910 - 3 1.2400000 0.4699870 +1 10 +S 4 + 1 528.5000000 0.0009400 + 2 79.3100000 0.0072140 + 3 18.0500000 0.0359750 + 4 5.0850000 0.1277820 S 1 - 1 0.2976000 1.0000000 + 1 1.6090000 1.0000000 +S 1 + 1 0.5363000 1.0000000 +S 1 + 1 0.1833000 1.0000000 P 1 - 1 1.2750000 1.0000000 + 1 5.9940000 1.0000000 +P 1 + 1 1.7450000 1.0000000 +P 1 + 1 0.5600000 1.0000000 +D 1 + 1 4.2990000 1.0000000 +D 1 + 1 1.2230000 1.0000000 +F 1 + 1 2.6800000 1.0000000 + diff --git a/input/weight b/input/weight index 6796e3b..cd6cb22 100644 --- a/input/weight +++ b/input/weight @@ -1,9 +1,25 @@ -1 3 -S 3 - 1 38.3600000 0.0238090 - 2 5.7700000 0.1548910 - 3 1.2400000 0.4699870 +1 10 +S 4 + 1 528.5000000 0.0009400 + 2 79.3100000 0.0072140 + 3 18.0500000 0.0359750 + 4 5.0850000 0.1277820 S 1 - 1 0.2976000 1.0000000 + 1 1.6090000 1.0000000 +S 1 + 1 0.5363000 1.0000000 +S 1 + 1 0.1833000 1.0000000 P 1 - 1 1.2750000 1.0000000 + 1 5.9940000 1.0000000 +P 1 + 1 1.7450000 1.0000000 +P 1 + 1 0.5600000 1.0000000 +D 1 + 1 4.2990000 1.0000000 +D 1 + 1 1.2230000 1.0000000 +F 1 + 1 2.6800000 1.0000000 + diff --git a/src/eDFT/eDFT_UKS.f90 b/src/eDFT/eDFT_UKS.f90 index fdec76d..315359e 100644 --- a/src/eDFT/eDFT_UKS.f90 +++ b/src/eDFT/eDFT_UKS.f90 @@ -310,7 +310,8 @@ subroutine eDFT_UKS(x_rung,x_DFA,c_rung,c_DFA,nEns,wEns,aCC_w1,aCC_w2,nGrid,weig ! Coulomb energy EJ(1) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - EJ(2) = trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) + EJ(2) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & + + 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) !!!!!! EJ(3) = 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,2))) ! Exchange energy diff --git a/src/eDFT/unrestricted_density_matrix.f90 b/src/eDFT/unrestricted_density_matrix.f90 index 922b06d..53ec50d 100644 --- a/src/eDFT/unrestricted_density_matrix.f90 +++ b/src/eDFT/unrestricted_density_matrix.f90 @@ -33,8 +33,11 @@ subroutine unrestricted_density_matrix(nBas,nEns,nO,c,P) iEns = 2 P(:,:,1,iEns) = matmul(c(:,1:nO(1) ,1),transpose(c(:,1:nO(1) ,1))) - P(:,:,2,iEns) = matmul(c(:,1:nO(2)-1,2),transpose(c(:,1:nO(2)-1,2))) - + if (nO(2) > 1) then + P(:,:,2,iEns) = matmul(c(:,1:nO(2)-1,2),transpose(c(:,1:nO(2)-1,2))) + else + P(:,:,2,iEns) = 0.d0 + end if ! (N+1)-electron state: remove spin-up electrons iEns = 3 diff --git a/src/eDFT/unrestricted_individual_energy.f90 b/src/eDFT/unrestricted_individual_energy.f90 index 0872900..b349cdb 100644 --- a/src/eDFT/unrestricted_individual_energy.f90 +++ b/src/eDFT/unrestricted_individual_energy.f90 @@ -64,6 +64,8 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered double precision,external :: trace_matrix integer :: ispin,iEns + + double precision,external :: electron_number ! Output variables @@ -103,7 +105,9 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered EJ(1,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,1))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,1))) - EJ(2,iEns) = 2.0d0*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & + EJ(2,iEns) = trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & + + trace_matrix(nBas,matmul(P(:,:,2,iEns),J(:,:,1))) & + ! 2.0d0*trace_matrix(nBas,matmul(P(:,:,1,iEns),J(:,:,2))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,1),J(:,:,2))) & - 0.5d0*trace_matrix(nBas,matmul(Pw(:,:,2),J(:,:,1))) @@ -112,6 +116,21 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered end do +!------------------------------------------------------------------------ +! Checking Hartree contributions for each individual states +!------------------------------------------------------------------------ + + print*,'Hartree contributions for each individual states' + print*,'' + print*,'' + print*,'EJ(aa,1)=',EJ(1,1),'EJ(ab,1)=',EJ(2,1),'EJ(bb,1)=',EJ(3,1) + print*,'' + print*,'EJ(aa,2)=',EJ(1,2),'EJ(ab,2)=',EJ(2,2),'EJ(bb,2)=',EJ(3,2) + print*,'' + print*,'EJ(aa,3)=',EJ(1,3),'EJ(ab,3)=',EJ(2,3),'EJ(bb,3)=',EJ(3,3) + print*,'' + + !------------------------------------------------------------------------ ! Individual exchange energy !------------------------------------------------------------------------ @@ -121,9 +140,38 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered call exchange_individual_energy(x_rung,x_DFA,LDA_centered,nEns,wEns,aCC_w1,aCC_w2,nGrid,weight,nBas,ERI, & Pw(:,:,ispin),P(:,:,ispin,iEns),rhow(:,ispin),drhow(:,:,ispin), & rho(:,ispin,iEns),drho(:,:,ispin,iEns),Ex(ispin,iEns)) - print*,'Ex = ',Ex(ispin,iEns) end do end do + +!------------------------------------------------------------------------ +! Checking exchange contributions for each individual states +!------------------------------------------------------------------------ + print*,'' + print*,'' + print*,'Exchange contributions for each individual states' + print*,'' + print*,'' + print*,'Ex(aa,1) =' ,Ex(1,1),'Ex(bb,1) =' ,Ex(2,1) + print*,'' + print*,'Ex(aa,2) =' ,Ex(1,2),'Ex(bb,2) =' ,Ex(2,2) + print*,'' + print*,'Ex(aa,3) =' ,Ex(1,3),'Ex(bb,3) =' ,Ex(2,3) + +!------------------------------------------------------------------------ +! Checking number of alpha and beta electrons for each individual states +!------------------------------------------------------------------------ + print*,'' + print*,'' + print*,'Checking number of alpha and beta electrons for each individual states' + print*,'' + print*,'' + print*,'nEl(a,1) = ',electron_number(nGrid,weight,rho(:,1,1)),'nEl(b,1) = ',electron_number(nGrid,weight,rho(:,2,1)) + print*,'' + print*,'nEl(a,2) = ',electron_number(nGrid,weight,rho(:,1,2)),'nEl(b,2) = ',electron_number(nGrid,weight,rho(:,2,2)) + print*,'' + print*,'nEl(a,3) = ',electron_number(nGrid,weight,rho(:,1,3)),'nEl(b,3) = ',electron_number(nGrid,weight,rho(:,2,3)) + + !------------------------------------------------------------------------ ! Individual correlation energy