mirror of https://gitlab.com/scemama/eplf
Multideterminant corrected. Bug in MPI?
This commit is contained in:
Anthony Scemama
parent
f19db06333
commit
0d69dec304
|
|
@ -38,15 +38,25 @@ BEGIN_PROVIDER [ real, density_alpha_value_p ]
|
|||
enddo
|
||||
|
||||
! TODO vectorization
|
||||
integer :: k,j,l
|
||||
integer :: k,j,l, ik, il
|
||||
real :: buffer
|
||||
PROVIDE det
|
||||
PROVIDE elec_alpha_num
|
||||
do k=1,det_num
|
||||
do l=1,det_num
|
||||
buffer = 0.
|
||||
do i=1,elec_alpha_num-mo_closed_num
|
||||
buffer += mo_value_p(det(i,k,1))*mo_value_p(det(i,l,1))
|
||||
enddo
|
||||
density_alpha_value_p += det_coef(k)*det_coef(l)*buffer
|
||||
|
||||
if (det_exc(k,l,3) == 0) then
|
||||
buffer = 0.
|
||||
do i=1,elec_alpha_num-mo_closed_num
|
||||
buffer += mo_value_p(det(i,k,1))*mo_value_p(det(i,l,1))
|
||||
enddo
|
||||
density_alpha_value_p += det_coef(k)*det_coef(l)*buffer
|
||||
else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,1) == 1) ) then
|
||||
call get_single_excitation(k,l,ik,il,1)
|
||||
buffer = mo_value_p(ik)*mo_value_p(il)
|
||||
density_alpha_value_p += det_coef(k)*det_coef(l)*buffer
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
|
@ -65,15 +75,23 @@ BEGIN_PROVIDER [ real, density_beta_value_p ]
|
|||
enddo
|
||||
|
||||
! TODO vectorization
|
||||
integer :: k,j,l
|
||||
integer :: k,j,l, ik, il
|
||||
real :: buffer
|
||||
PROVIDE det
|
||||
PROVIDE elec_beta_num
|
||||
do k=1,det_num
|
||||
do l=1,det_num
|
||||
buffer = 0.
|
||||
do i=1,elec_beta_num-mo_closed_num
|
||||
buffer += mo_value_p(det(i,k,2))*mo_value_p(det(i,l,2))
|
||||
enddo
|
||||
density_beta_value_p += det_coef(k)*det_coef(l)*buffer
|
||||
if (det_exc(k,l,3) == 0) then
|
||||
buffer = 0.
|
||||
do i=1,elec_beta_num-mo_closed_num
|
||||
buffer += mo_value_p(det(i,k,2))*mo_value_p(det(i,l,2))
|
||||
enddo
|
||||
density_beta_value_p += det_coef(k)*det_coef(l)*buffer
|
||||
else if ( (det_exc(k,l,3) == 1).and.(det_exc(k,l,2) == 1) ) then
|
||||
call get_single_excitation(k,l,ik,il,2)
|
||||
buffer = mo_value_p(ik)*mo_value_p(il)
|
||||
density_beta_value_p += det_coef(k)*det_coef(l)*buffer
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
|
||||
|
|
@ -98,9 +116,12 @@ BEGIN_PROVIDER [ double precision, density_alpha_grad_p, (3) ]
|
|||
density_alpha_grad_p(3) += 2.*mo_grad_p(i,3)*mo_value_p(i)
|
||||
enddo
|
||||
|
||||
! TODO vectorization
|
||||
! TODO Faux !! => Regles de Slater
|
||||
integer :: k,j,l
|
||||
real :: buffer(3)
|
||||
if (det_num > 1) then
|
||||
call abrt(irp_here,'Need to implement Slater rules for gradient')
|
||||
endif
|
||||
do k=1,det_num
|
||||
do l=1,det_num
|
||||
buffer(1) = 0.
|
||||
|
|
@ -142,6 +163,9 @@ BEGIN_PROVIDER [ double precision, density_beta_grad_p, (3) ]
|
|||
enddo
|
||||
|
||||
! TODO vectorization
|
||||
if (det_num > 1) then
|
||||
call abrt(irp_here,'Need to implement Slater rules for gradient')
|
||||
endif
|
||||
integer :: k,j,l
|
||||
real :: buffer(3)
|
||||
do k=1,det_num
|
||||
|
|
@ -189,6 +213,9 @@ BEGIN_PROVIDER [ double precision, density_alpha_lapl_p ]
|
|||
! TODO vectorization
|
||||
integer :: k,j,l
|
||||
real :: buffer
|
||||
if (det_num > 1) then
|
||||
call abrt(irp_here,'Need to implement Slater rules for gradient')
|
||||
endif
|
||||
do k=1,det_num
|
||||
do l=1,det_num
|
||||
buffer = 0.
|
||||
|
|
@ -225,6 +252,9 @@ BEGIN_PROVIDER [ double precision, density_beta_lapl_p ]
|
|||
density_beta_lapl_p *= 2.
|
||||
|
||||
! TODO vectorization
|
||||
if (det_num > 1) then
|
||||
call abrt(irp_here,'Need to implement Slater rules for laplacian')
|
||||
endif
|
||||
integer :: k,j,l
|
||||
real :: buffer
|
||||
do k=1,det_num
|
||||
|
|
|
|||
164
src/det.irp.f
164
src/det.irp.f
|
|
@ -33,4 +33,168 @@ BEGIN_PROVIDER [ integer, det, (elec_alpha_num-mo_closed_num,det_num,2) ]
|
|||
END_PROVIDER
|
||||
|
||||
|
||||
BEGIN_PROVIDER [ integer, det_exc, (det_num, det_num, 3) ]
|
||||
implicit none
|
||||
BEGIN_DOC
|
||||
! Degree of excitation between two determinants. The sign is the phase.
|
||||
END_DOC
|
||||
|
||||
integer :: p
|
||||
do p=1,2
|
||||
|
||||
integer :: k, l
|
||||
do l=1,det_num
|
||||
det_exc(l,l,p) = 0
|
||||
do k=l+1,det_num
|
||||
det_exc(k,l,p) = 0
|
||||
|
||||
! Excitation degree
|
||||
integer :: i, j
|
||||
do i=1,elec_num_2(p)-mo_closed_num
|
||||
|
||||
logical :: found
|
||||
found = .False.
|
||||
do j=1,elec_num_2(p)-mo_closed_num
|
||||
if (det(j,l,p) == det(i,k,p)) then
|
||||
found = .True.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (.not.found) then
|
||||
det_exc(k,l,p) += 1
|
||||
endif
|
||||
|
||||
enddo
|
||||
|
||||
det_exc(l,k,p) = det_exc(k,l,p)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
do l=1,det_num
|
||||
do k=l+1,det_num
|
||||
det_exc(k,l,3) = det_exc(k,l,1) + det_exc(k,l,2)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Phase
|
||||
do p=1,2
|
||||
do i=mo_closed_num,mo_num
|
||||
integer :: det_pos(det_num)
|
||||
do k=1,det_num
|
||||
det_pos(k) = 0
|
||||
do j=1,elec_num_2(p)-mo_closed_num
|
||||
if (det(j,k,p) == i) then
|
||||
det_pos(k) = j
|
||||
endif
|
||||
enddo
|
||||
enddo
|
||||
do k=1,det_num
|
||||
do l=k+1,det_num
|
||||
det_exc(k,l,3) *= -2*mod( (det_pos(k)+det_pos(l)), 2 )+1
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
|
||||
do l=1,det_num
|
||||
do k=l+1,det_num
|
||||
det_exc(l,k,3) = det_exc(k,l,3)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
subroutine get_single_excitation(k,l,m,n,p)
|
||||
implicit none
|
||||
integer, intent(in) :: k, l ! determinant indices
|
||||
integer, intent(out) :: m, n ! m->n excitation
|
||||
integer, intent(in) :: p ! spin
|
||||
|
||||
logical :: found
|
||||
integer :: i,j
|
||||
m=0
|
||||
n=0
|
||||
do j=1,elec_num_2(p)-mo_closed_num
|
||||
found = .False.
|
||||
do i=1,elec_num_2(p)-mo_closed_num
|
||||
if (det(j,k,p) == det(i,l,p)) then
|
||||
found = .True.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (.not.found) then
|
||||
m = det(j,k,p)
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
do i=1,elec_num_2(p)-mo_closed_num
|
||||
found = .False.
|
||||
do j=1,elec_num_2(p)-mo_closed_num
|
||||
if (det(i,k,p) == det(i,l,p)) then
|
||||
found = .True.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (.not.found) then
|
||||
n = det(i,l,p)
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
subroutine get_double_excitation(k,l,m,n,r,s,p)
|
||||
implicit none
|
||||
integer, intent(in) :: k, l ! determinant indices
|
||||
integer, intent(out) :: m, n ! m->n excitation
|
||||
integer, intent(out) :: r, s ! r->s excitation
|
||||
integer, intent(in) :: p ! spin
|
||||
|
||||
logical :: found
|
||||
integer :: i,j
|
||||
m=0
|
||||
n=0
|
||||
r=0
|
||||
s=0
|
||||
do j=1,elec_num_2(p)-mo_closed_num
|
||||
found = .False.
|
||||
do i=1,elec_num_2(p)-mo_closed_num
|
||||
if (det(j,k,p) == det(i,l,p)) then
|
||||
found = .True.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (.not.found) then
|
||||
if (m == 0) then
|
||||
m = det(j,k,p)
|
||||
else
|
||||
r = det(j,k,p)
|
||||
exit
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
|
||||
do i=1,elec_num_2(p)-mo_closed_num
|
||||
found = .False.
|
||||
do j=1,elec_num_2(p)-mo_closed_num
|
||||
if (det(i,k,p) == det(i,l,p)) then
|
||||
found = .True.
|
||||
exit
|
||||
endif
|
||||
enddo
|
||||
if (.not.found) then
|
||||
if (n == 0) then
|
||||
n = det(i,l,p)
|
||||
else
|
||||
s = det(i,l,p)
|
||||
exit
|
||||
endif
|
||||
endif
|
||||
enddo
|
||||
|
||||
end
|
||||
|
||||
|
|
|
|||
|
|
@ -11,6 +11,7 @@ BEGIN_PROVIDER [ real, eplf_gamma ]
|
|||
!eplf_gamma = (4./(3.*N)*pi*density_value_p)**(2./3.) * eps
|
||||
eplf_gamma = density_value_p * eps
|
||||
!eplf_gamma = 1.e10
|
||||
!eplf_gamma = 1.e5
|
||||
END_PROVIDER
|
||||
|
||||
BEGIN_PROVIDER [ double precision, ao_eplf_integral_matrix, (ao_num,ao_num) ]
|
||||
|
|
@ -79,41 +80,116 @@ END_PROVIDER
|
|||
|
||||
do j=1,mo_closed_num
|
||||
do i=1,mo_closed_num
|
||||
eplf_up_up += 2.d0*mo_value_p(i)* ( &
|
||||
eplf_up_up += 2.d0* mo_value_p(i)* ( &
|
||||
mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
|
||||
mo_value_p(j)*mo_eplf_integral_matrix(i,j) )
|
||||
eplf_up_dn += 2.d0*mo_value_p(j)**2 * &
|
||||
mo_eplf_integral_matrix(i,i)
|
||||
eplf_up_dn += 2.d0* mo_value_p(i)*mo_value_p(i)* &
|
||||
mo_eplf_integral_matrix(j,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
integer :: k,l,m,n,p
|
||||
integer :: ik,il,jk,jl
|
||||
double precision :: ckl
|
||||
double precision :: phase
|
||||
integer :: exc
|
||||
|
||||
PROVIDE det
|
||||
PROVIDE elec_num_2
|
||||
|
||||
do k=1,det_num
|
||||
do l=1,det_num
|
||||
|
||||
ckl = det_coef(k)*det_coef(l)
|
||||
|
||||
do p=1,2
|
||||
do m=1,elec_num_2(p)-mo_closed_num
|
||||
j = det(m,k,p)
|
||||
do n=1,elec_num_2(p)-mo_closed_num
|
||||
i = det(n,l,p)
|
||||
eplf_up_up += ckl*mo_value_p(i)* ( &
|
||||
mo_value_p(i)*mo_eplf_integral_matrix(j,j) - &
|
||||
mo_value_p(j)*mo_eplf_integral_matrix(i,j) )
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
exc = det_exc(k,l,3)
|
||||
|
||||
do m=1,elec_beta_num-mo_closed_num
|
||||
j = det(m,k,2)
|
||||
do n=1,elec_alpha_num-mo_closed_num
|
||||
i = det(n,k,1)
|
||||
eplf_up_dn += ckl * ( mo_value_p(i)**2 * mo_eplf_integral_matrix(j,j) &
|
||||
+ mo_value_p(j)**2 * mo_eplf_integral_matrix(i,i) )
|
||||
enddo
|
||||
enddo
|
||||
if ( exc < 0 ) then
|
||||
phase = -0.5d0
|
||||
exc = -exc
|
||||
else
|
||||
phase = 0.5d0
|
||||
endif
|
||||
|
||||
if ( exc == 0 ) then
|
||||
! Sum over all alpha-alpha and beta-beta interactions
|
||||
do p=1,2
|
||||
do m=1,elec_num_2(p)-mo_closed_num
|
||||
jk = det(m,k,p)
|
||||
jl = det(m,l,p)
|
||||
do n=1,elec_num_2(p)-mo_closed_num
|
||||
ik = det(n,k,p)
|
||||
il = det(n,l,p)
|
||||
eplf_up_up += phase*ckl*mo_value_p(ik)* ( &
|
||||
mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
|
||||
mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Sum over all alpha-beta interactions
|
||||
do m=1,elec_beta_num-mo_closed_num
|
||||
jk = det(m,k,2)
|
||||
jl = det(m,l,2)
|
||||
do n=1,elec_alpha_num-mo_closed_num
|
||||
ik = det(n,k,1)
|
||||
il = det(n,l,1)
|
||||
eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
|
||||
+ mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
else if ( exc == 1 ) then
|
||||
|
||||
do p=1,2
|
||||
if ( det_exc(k,l,p) == 1 ) then
|
||||
! Sum over only the sigma-sigma interactions involving the excitation
|
||||
call get_single_excitation(k,l,ik,il,p)
|
||||
do m=1,elec_num_2(p)-mo_closed_num
|
||||
jk = det(m,k,p)
|
||||
jl = det(m,l,p)
|
||||
eplf_up_up += phase*ckl*mo_value_p(ik)* ( &
|
||||
mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
|
||||
mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
|
||||
enddo
|
||||
|
||||
! Sum over only the sigma-(sigma_bar) interactions involving the excitation
|
||||
integer :: p2
|
||||
p2 = 1+mod(p,2)
|
||||
do m=1,elec_num_2(p2)-mo_closed_num
|
||||
jk = det(m,k,p2)
|
||||
jl = det(m,l,p2)
|
||||
eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
|
||||
+ mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
|
||||
enddo
|
||||
endif
|
||||
enddo
|
||||
|
||||
else if (exc == 2) then
|
||||
|
||||
if ( ( det_exc(k,l,1) == 2 ).or.( det_exc(k,l,2) == 2 ) ) then
|
||||
|
||||
! Consider only the double excitations of same-spin electrons
|
||||
if ( det_exc(k,l,1) == 2 ) then
|
||||
call get_double_excitation(k,l,ik,jk,il,jl,1)
|
||||
else if ( det_exc(k,l,2) == 2 ) then
|
||||
call get_double_excitation(k,l,ik,jk,il,jl,2)
|
||||
endif
|
||||
|
||||
eplf_up_up += phase*ckl*mo_value_p(ik)* ( &
|
||||
mo_value_p(il)*mo_eplf_integral_matrix(jk,jl) - &
|
||||
mo_value_p(jl)*mo_eplf_integral_matrix(jk,il) )
|
||||
|
||||
else if ( det_exc(k,l,1) == 1 ) then
|
||||
! Consider only the double excitations of opposite-spin electrons
|
||||
call get_single_excitation(k,l,ik,jk,1)
|
||||
call get_single_excitation(k,l,il,jl,2)
|
||||
|
||||
eplf_up_dn += phase*ckl * ( mo_value_p(ik)*mo_value_p(il) * mo_eplf_integral_matrix(jk,jl) &
|
||||
+ mo_value_p(jk)*mo_value_p(jl) * mo_eplf_integral_matrix(ik,il) )
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
|
|
|||
|
|
@ -100,7 +100,7 @@ BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
|
|||
do iy=1,grid_y_num
|
||||
point(2) = grid_origin(2)+(iy-1)*grid_step(2)
|
||||
do ix=1,grid_x_num
|
||||
icount = icount-1
|
||||
icount -= 1
|
||||
if (icount == mpi_rank) then
|
||||
point(1) = grid_origin(1)+(ix-1)*grid_step(1)
|
||||
TOUCH point
|
||||
|
|
@ -116,7 +116,7 @@ BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
|
|||
IRP_IF MPI
|
||||
integer :: dim, ierr
|
||||
do iz=1,grid_z_num
|
||||
real :: buffer(grid_x_num*grid_y_num)
|
||||
real :: buffer(grid_x_num*grid_y_num+1)
|
||||
icount = 0
|
||||
do iy=1,grid_y_num
|
||||
do ix=1,grid_x_num
|
||||
|
|
@ -127,6 +127,10 @@ BEGIN_PROVIDER [ real, grid_$X, (grid_x_num,grid_y_num,grid_z_num) ]
|
|||
dim = grid_x_num * grid_y_num
|
||||
call MPI_REDUCE(buffer,grid_$X(1,1,iz),dim,mpi_real, &
|
||||
mpi_sum,0,MPI_COMM_WORLD,ierr)
|
||||
if (ierr /= MPI_SUCCESS) then
|
||||
call abrt(irp_here,'Unable to fetch buffer')
|
||||
endif
|
||||
call barrier
|
||||
enddo
|
||||
IRP_ENDIF
|
||||
|
||||
|
|
|
|||
Loading…
Reference in New Issue