PTEROSOR/QuAcK
10
1
Fork 0
mirror of https://github.com/pfloos/quack synced 2024-11-05 05:33:50 +01:00
Code Releases Activity

still looking for some bug

Browse Source
This commit is contained in:
Clotilde Marut 2020-07-05 19:04:08 +02:00
parent f3e7222e49
commit 372ccd2e4e
5 changed files with 103 additions and 19 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
input/basis
View File

@ -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

30
input/weight
View File

@ -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

3
src/eDFT/eDFT_UKS.f90
View File

@ -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

5
src/eDFT/unrestricted_density_matrix.f90
View File

@ -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)))
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

52
src/eDFT/unrestricted_individual_energy.f90
View File

@ -65,6 +65,8 @@ subroutine unrestricted_individual_energy(x_rung,x_DFA,c_rung,c_DFA,LDA_centered
integer :: ispin,iEns
double precision,external :: electron_number
! Output variables
double precision,intent(out) :: E(nEns)
@ -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,10 +140,39 @@ 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
!------------------------------------------------------------------------