mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-18 11:23:38 +01:00
numerical integrals: 1 shot --> blocks over r2
This commit is contained in:
AbdAmmar
parent
63f6404f8e
commit
38f8b96d41
@ -425,7 +425,6 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
|
|||||||
|
|
||||||
! an additional term is added here directly instead of
|
! an additional term is added here directly instead of
|
||||||
! being added in int2_grad1_u12_square_ao for performance
|
! being added in int2_grad1_u12_square_ao for performance
|
||||||
! note that the factor
|
|
||||||
|
|
||||||
PROVIDE int2_u2_j1b2
|
PROVIDE int2_u2_j1b2
|
||||||
|
|
||||||
@ -465,25 +464,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
|
|||||||
! ---
|
! ---
|
||||||
|
|
||||||
deallocate(b_mat)
|
deallocate(b_mat)
|
||||||
|
|
||||||
call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)
|
call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)
|
||||||
|
|
||||||
!!$OMP PARALLEL &
|
|
||||||
!!$OMP DEFAULT (NONE) &
|
|
||||||
!!$OMP PRIVATE (i, j, k, l) &
|
|
||||||
!!$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num)
|
|
||||||
!!$OMP DO SCHEDULE (static)
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! do l = 1, ao_num
|
|
||||||
! do i = 1, ao_num
|
|
||||||
! do k = 1, ao_num
|
|
||||||
! tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
!!$OMP END DO
|
|
||||||
!!$OMP END PARALLEL
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
if(write_tc_integ.and.mpi_master) then
|
if(write_tc_integ.and.mpi_master) then
|
||||||
|
|||||||
@ -67,72 +67,6 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
|
|||||||
deallocate(tmp, b_mat)
|
deallocate(tmp, b_mat)
|
||||||
|
|
||||||
call sum_A_At(tc_grad_square_ao_test(1,1,1,1), ao_num*ao_num)
|
call sum_A_At(tc_grad_square_ao_test(1,1,1,1), ao_num*ao_num)
|
||||||
!do i = 1, ao_num
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! do k = i, ao_num
|
|
||||||
|
|
||||||
! do l = max(j,k), ao_num
|
|
||||||
! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
|
|
||||||
! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
|
|
||||||
! end do
|
|
||||||
|
|
||||||
! !if (j.eq.k) then
|
|
||||||
! ! do l = j+1, ao_num
|
|
||||||
! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
|
|
||||||
! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
|
|
||||||
! ! end do
|
|
||||||
! !else
|
|
||||||
! ! do l = j, ao_num
|
|
||||||
! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
|
|
||||||
! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
|
|
||||||
! ! enddo
|
|
||||||
! !endif
|
|
||||||
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
!enddo
|
|
||||||
!tc_grad_square_ao_test = 2.d0 * tc_grad_square_ao_test
|
|
||||||
! !$OMP PARALLEL &
|
|
||||||
! !$OMP DEFAULT (NONE) &
|
|
||||||
! !$OMP PRIVATE (i, j, k, l) &
|
|
||||||
! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
|
|
||||||
! !$OMP DO SCHEDULE (static)
|
|
||||||
! integer :: ii
|
|
||||||
! ii = 0
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! do l = 1, ao_num
|
|
||||||
! do i = 1, ao_num
|
|
||||||
! do k = 1, ao_num
|
|
||||||
! if((i.lt.j) .and. (k.lt.l)) cycle
|
|
||||||
! ii = ii + 1
|
|
||||||
! tc_grad_square_ao_test(k,i,l,j) = tc_grad_square_ao_test(k,i,l,j) + tc_grad_square_ao_test(l,j,k,i)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! print *, ' ii =', ii
|
|
||||||
! !$OMP END DO
|
|
||||||
! !$OMP END PARALLEL
|
|
||||||
|
|
||||||
! !$OMP PARALLEL &
|
|
||||||
! !$OMP DEFAULT (NONE) &
|
|
||||||
! !$OMP PRIVATE (i, j, k, l) &
|
|
||||||
! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
|
|
||||||
! !$OMP DO SCHEDULE (static)
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! do l = 1, ao_num
|
|
||||||
! do i = 1, j-1
|
|
||||||
! do k = 1, l-1
|
|
||||||
! ii = ii + 1
|
|
||||||
! tc_grad_square_ao_test(k,i,l,j) = tc_grad_square_ao_test(l,j,k,i)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! print *, ' ii =', ii
|
|
||||||
! print *, ao_num * ao_num * ao_num * ao_num
|
|
||||||
! !$OMP END DO
|
|
||||||
! !$OMP END PARALLEL
|
|
||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
|||||||
@ -32,7 +32,8 @@
|
|||||||
grad1_u12_num = 0.d0
|
grad1_u12_num = 0.d0
|
||||||
grad1_u12_squared_num = 0.d0
|
grad1_u12_squared_num = 0.d0
|
||||||
|
|
||||||
if(j1b_type .eq. 100) then
|
if( (j1b_type .eq. 100) .or. &
|
||||||
|
(j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
@ -111,19 +112,120 @@
|
|||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then
|
else
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
|
||||||
!$OMP DEFAULT (NONE) &
|
stop
|
||||||
!$OMP PRIVATE (ipoint, jpoint, r1, r2, grad1_u2b, dx, dy, dz) &
|
|
||||||
!$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points, &
|
|
||||||
!$OMP final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num)
|
|
||||||
!$OMP DO SCHEDULE (static)
|
|
||||||
do ipoint = 1, n_points_final_grid ! r1
|
|
||||||
|
|
||||||
r1(1) = final_grid_points(1,ipoint)
|
endif
|
||||||
r1(2) = final_grid_points(2,ipoint)
|
|
||||||
r1(3) = final_grid_points(3,ipoint)
|
END_PROVIDER
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine get_grad1_u12_r1_seq(r1, n_grid2, resx, resy, resz)
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! grad_1 u(r1,r2)
|
||||||
|
!
|
||||||
|
! this will be integrated numerically over r2:
|
||||||
|
! we use grid for r1 and extra_grid for r2
|
||||||
|
!
|
||||||
|
! for 99 < j1b_type < 199
|
||||||
|
!
|
||||||
|
! u(r1,r2) = j12_mu(r12) x v(r1) x v(r2)
|
||||||
|
! grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2)
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: n_grid2
|
||||||
|
double precision, intent(in) :: r1(3)
|
||||||
|
double precision, intent(out) :: resx(n_grid2), resy(n_grid2), resz(n_grid2)
|
||||||
|
|
||||||
|
double precision :: v1b_r1
|
||||||
|
double precision :: grad1_v1b(3)
|
||||||
|
double precision, allocatable :: v1b_r2(:)
|
||||||
|
double precision, allocatable :: u2b_r12(:)
|
||||||
|
double precision, allocatable :: gradx1_u2b(:), grady1_u2b(:), gradz1_u2b(:)
|
||||||
|
double precision, external :: j1b_nucl
|
||||||
|
|
||||||
|
PROVIDE j1b_type
|
||||||
|
PROVIDE final_grid_points_extra
|
||||||
|
|
||||||
|
if( (j1b_type .eq. 100) .or. &
|
||||||
|
(j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then
|
||||||
|
|
||||||
|
call grad1_j12_mu_r1_seq(r1, n_grid2, resx, resy, resz)
|
||||||
|
|
||||||
|
elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then
|
||||||
|
|
||||||
|
allocate(v1b_r2(n_grid2))
|
||||||
|
allocate(u2b_r12(n_grid2))
|
||||||
|
allocate(gradx1_u2b(n_grid2))
|
||||||
|
allocate(grady1_u2b(n_grid2))
|
||||||
|
allocate(gradz1_u2b(n_grid2))
|
||||||
|
|
||||||
|
v1b_r1 = j1b_nucl(r1)
|
||||||
|
call grad1_j1b_nucl(r1, grad1_v1b)
|
||||||
|
|
||||||
|
call j1b_nucl_r1_seq(n_grid2, v1b_r2)
|
||||||
|
call j12_mu_r1_seq(r1, n_grid2, u2b_r12)
|
||||||
|
call grad1_j12_mu_r1_seq(r1, n_grid2, gradx1_u2b, grady1_u2b, gradz1_u2b)
|
||||||
|
|
||||||
|
resx(:) = (gradx1_u2b(:) * v1b_r1 + u2b_r12(:) * grad1_v1b(1)) * v1b_r2(:)
|
||||||
|
resy(:) = (grady1_u2b(:) * v1b_r1 + u2b_r12(:) * grad1_v1b(2)) * v1b_r2(:)
|
||||||
|
resz(:) = (gradz1_u2b(:) * v1b_r1 + u2b_r12(:) * grad1_v1b(3)) * v1b_r2(:)
|
||||||
|
|
||||||
|
deallocate(v1b_r2, u2b_r12, gradx1_u2b, grady1_u2b, gradz1_u2b)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine get_grad1_u12_r1_seq
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine get_grad1_u12_squared_r1_seq(r1, n_grid2, res)
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! grad_1 u(r1,r2)
|
||||||
|
!
|
||||||
|
! this will be integrated numerically over r2:
|
||||||
|
! we use grid for r1 and extra_grid for r2
|
||||||
|
!
|
||||||
|
! for 99 < j1b_type < 199
|
||||||
|
!
|
||||||
|
! u(r1,r2) = j12_mu(r12) x v(r1) x v(r2)
|
||||||
|
! grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2)
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: n_grid2
|
||||||
|
double precision, intent(in) :: r1(3)
|
||||||
|
double precision, intent(out) :: res(n_grid2)
|
||||||
|
|
||||||
|
|
||||||
|
integer :: jpoint
|
||||||
|
double precision :: r2(3)
|
||||||
|
double precision :: v1b_r1, v1b_r2, u2b_r12
|
||||||
|
double precision :: grad1_v1b(3), grad1_u2b(3)
|
||||||
|
double precision :: dx, dy, dz
|
||||||
|
double precision, external :: j12_mu, j1b_nucl
|
||||||
|
|
||||||
|
PROVIDE j1b_type
|
||||||
|
PROVIDE final_grid_points_extra
|
||||||
|
|
||||||
|
if( (j1b_type .eq. 100) .or. &
|
||||||
|
(j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then
|
||||||
|
|
||||||
do jpoint = 1, n_points_extra_final_grid ! r2
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
|
||||||
@ -137,15 +239,30 @@
|
|||||||
dy = grad1_u2b(2)
|
dy = grad1_u2b(2)
|
||||||
dz = grad1_u2b(3)
|
dz = grad1_u2b(3)
|
||||||
|
|
||||||
grad1_u12_num(jpoint,ipoint,1) = dx
|
res(jpoint) = dx*dx + dy*dy + dz*dz
|
||||||
grad1_u12_num(jpoint,ipoint,2) = dy
|
enddo
|
||||||
grad1_u12_num(jpoint,ipoint,3) = dz
|
|
||||||
|
|
||||||
grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
|
elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then
|
||||||
|
|
||||||
|
v1b_r1 = j1b_nucl(r1)
|
||||||
|
call grad1_j1b_nucl(r1, grad1_v1b)
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
|
||||||
|
r2(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r2(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r2(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
v1b_r2 = j1b_nucl(r2)
|
||||||
|
u2b_r12 = j12_mu(r1, r2)
|
||||||
|
call grad1_j12_mu(r1, r2, grad1_u2b)
|
||||||
|
|
||||||
|
dx = (grad1_u2b(1) * v1b_r1 + u2b_r12 * grad1_v1b(1)) * v1b_r2
|
||||||
|
dy = (grad1_u2b(2) * v1b_r1 + u2b_r12 * grad1_v1b(2)) * v1b_r2
|
||||||
|
dz = (grad1_u2b(3) * v1b_r1 + u2b_r12 * grad1_v1b(3)) * v1b_r2
|
||||||
|
|
||||||
|
res(jpoint) = dx*dx + dy*dy + dz*dz
|
||||||
enddo
|
enddo
|
||||||
enddo
|
|
||||||
!$OMP END DO
|
|
||||||
!$OMP END PARALLEL
|
|
||||||
|
|
||||||
else
|
else
|
||||||
|
|
||||||
@ -154,7 +271,88 @@
|
|||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
END_PROVIDER
|
return
|
||||||
|
end subroutine get_grad1_u12_squared_r1_seq
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine get_grad1_u12_withsq_r1_seq(r1, n_grid2, resx, resy, resz, res)
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! grad_1 u(r1,r2)
|
||||||
|
!
|
||||||
|
! this will be integrated numerically over r2:
|
||||||
|
! we use grid for r1 and extra_grid for r2
|
||||||
|
!
|
||||||
|
! for 99 < j1b_type < 199
|
||||||
|
!
|
||||||
|
! u(r1,r2) = j12_mu(r12) x v(r1) x v(r2)
|
||||||
|
! grad1 u(r1, r2) = [(grad1 j12_mu) v(r1) + j12_mu grad1 v(r1)] v(r2)
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: n_grid2
|
||||||
|
double precision, intent(in) :: r1(3)
|
||||||
|
double precision, intent(out) :: resx(n_grid2), resy(n_grid2), resz(n_grid2), res(n_grid2)
|
||||||
|
|
||||||
|
integer :: jpoint
|
||||||
|
double precision :: v1b_r1
|
||||||
|
double precision :: grad1_v1b(3)
|
||||||
|
double precision, allocatable :: v1b_r2(:)
|
||||||
|
double precision, allocatable :: u2b_r12(:)
|
||||||
|
double precision, allocatable :: gradx1_u2b(:), grady1_u2b(:), gradz1_u2b(:)
|
||||||
|
double precision, external :: j1b_nucl
|
||||||
|
|
||||||
|
PROVIDE j1b_type
|
||||||
|
PROVIDE final_grid_points_extra
|
||||||
|
|
||||||
|
if( (j1b_type .eq. 100) .or. &
|
||||||
|
(j1b_type .ge. 200) .and. (j1b_type .lt. 300) ) then
|
||||||
|
|
||||||
|
call grad1_j12_mu_r1_seq(r1, n_grid2, resx, resy, resz)
|
||||||
|
do jpoint = 1, n_points_extra_final_grid
|
||||||
|
res(jpoint) = resx(jpoint) * resx(jpoint) &
|
||||||
|
+ resy(jpoint) * resy(jpoint) &
|
||||||
|
+ resz(jpoint) * resz(jpoint)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then
|
||||||
|
|
||||||
|
allocate(v1b_r2(n_grid2))
|
||||||
|
allocate(u2b_r12(n_grid2))
|
||||||
|
allocate(gradx1_u2b(n_grid2))
|
||||||
|
allocate(grady1_u2b(n_grid2))
|
||||||
|
allocate(gradz1_u2b(n_grid2))
|
||||||
|
|
||||||
|
v1b_r1 = j1b_nucl(r1)
|
||||||
|
call grad1_j1b_nucl(r1, grad1_v1b)
|
||||||
|
|
||||||
|
call j1b_nucl_r1_seq(n_grid2, v1b_r2)
|
||||||
|
call j12_mu_r1_seq(r1, n_grid2, u2b_r12)
|
||||||
|
call grad1_j12_mu_r1_seq(r1, n_grid2, gradx1_u2b, grady1_u2b, gradz1_u2b)
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid
|
||||||
|
resx(jpoint) = (gradx1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(1)) * v1b_r2(jpoint)
|
||||||
|
resy(jpoint) = (grady1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(2)) * v1b_r2(jpoint)
|
||||||
|
resz(jpoint) = (gradz1_u2b(jpoint) * v1b_r1 + u2b_r12(jpoint) * grad1_v1b(3)) * v1b_r2(jpoint)
|
||||||
|
res (jpoint) = resx(jpoint) * resx(jpoint) &
|
||||||
|
+ resy(jpoint) * resy(jpoint) &
|
||||||
|
+ resz(jpoint) * resz(jpoint)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
deallocate(v1b_r2, u2b_r12, gradx1_u2b, grady1_u2b, gradz1_u2b)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine get_grad1_u12_withsq_r1_seq
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
@ -190,6 +388,7 @@ end function j12_mu
|
|||||||
subroutine grad1_j12_mu(r1, r2, grad)
|
subroutine grad1_j12_mu(r1, r2, grad)
|
||||||
|
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
|
!
|
||||||
! gradient of j(mu(r1,r2),r12) form of jastrow.
|
! gradient of j(mu(r1,r2),r12) form of jastrow.
|
||||||
!
|
!
|
||||||
! if mu(r1,r2) = cst ---> j1b_type < 200 and
|
! if mu(r1,r2) = cst ---> j1b_type < 200 and
|
||||||
@ -199,9 +398,10 @@ subroutine grad1_j12_mu(r1, r2, grad)
|
|||||||
! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and
|
! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and
|
||||||
!
|
!
|
||||||
! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2)
|
! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2)
|
||||||
!
|
|
||||||
! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
|
! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
|
||||||
|
!
|
||||||
END_DOC
|
END_DOC
|
||||||
|
|
||||||
include 'constants.include.F'
|
include 'constants.include.F'
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
@ -851,3 +1051,240 @@ subroutine f_mu_and_deriv_mu_simple(rho,alpha,mu0,beta,f_mu,d_drho_f_mu)
|
|||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine grad1_j12_mu_r1_seq(r1, n_grid2, gradx, grady, gradz)
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! gradient of j(mu(r1,r2),r12) form of jastrow.
|
||||||
|
!
|
||||||
|
! if mu(r1,r2) = cst ---> j1b_type < 200 and
|
||||||
|
!
|
||||||
|
! d/dx1 j(mu,r12) = 0.5 * (1 - erf(mu *r12))/r12 * (x1 - x2)
|
||||||
|
!
|
||||||
|
! if mu(r1,r2) /= cst ---> 200 < j1b_type < 300 and
|
||||||
|
!
|
||||||
|
! d/dx1 j(mu(r1,r2),r12) = exp(-(mu(r1,r2)*r12)**2) /(2 *sqrt(pi) * mu(r1,r2)**2 ) d/dx1 mu(r1,r2)
|
||||||
|
! + 0.5 * (1 - erf(mu(r1,r2) *r12))/r12 * (x1 - x2)
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
include 'constants.include.F'
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer , intent(in) :: n_grid2
|
||||||
|
double precision, intent(in) :: r1(3)
|
||||||
|
double precision, intent(out) :: gradx(n_grid2)
|
||||||
|
double precision, intent(out) :: grady(n_grid2)
|
||||||
|
double precision, intent(out) :: gradz(n_grid2)
|
||||||
|
|
||||||
|
integer :: jpoint
|
||||||
|
double precision :: r2(3)
|
||||||
|
double precision :: dx, dy, dz, r12, tmp
|
||||||
|
|
||||||
|
if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
|
||||||
|
r2(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r2(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r2(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
dx = r1(1) - r2(1)
|
||||||
|
dy = r1(2) - r2(2)
|
||||||
|
dz = r1(3) - r2(3)
|
||||||
|
|
||||||
|
r12 = dsqrt(dx * dx + dy * dy + dz * dz)
|
||||||
|
if(r12 .lt. 1d-10) return
|
||||||
|
|
||||||
|
tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
|
||||||
|
|
||||||
|
gradx(jpoint) = tmp * dx
|
||||||
|
grady(jpoint) = tmp * dy
|
||||||
|
gradz(jpoint) = tmp * dz
|
||||||
|
enddo
|
||||||
|
|
||||||
|
elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then
|
||||||
|
|
||||||
|
double precision :: mu_val, mu_tmp, mu_der(3)
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
|
||||||
|
r2(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r2(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r2(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
dx = r1(1) - r2(1)
|
||||||
|
dy = r1(2) - r2(2)
|
||||||
|
dz = r1(3) - r2(3)
|
||||||
|
r12 = dsqrt(dx * dx + dy * dy + dz * dz)
|
||||||
|
|
||||||
|
call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
|
||||||
|
mu_tmp = mu_val * r12
|
||||||
|
tmp = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
|
||||||
|
gradx(jpoint) = tmp * mu_der(1)
|
||||||
|
grady(jpoint) = tmp * mu_der(2)
|
||||||
|
gradz(jpoint) = tmp * mu_der(3)
|
||||||
|
|
||||||
|
if(r12 .lt. 1d-10) return
|
||||||
|
tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
|
||||||
|
|
||||||
|
gradx(jpoint) = gradx(jpoint) + tmp * dx
|
||||||
|
grady(jpoint) = grady(jpoint) + tmp * dy
|
||||||
|
gradz(jpoint) = gradz(jpoint) + tmp * dz
|
||||||
|
enddo
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine grad1_j12_mu_r1_seq
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine j12_mu_r1_seq(r1, n_grid2, res)
|
||||||
|
|
||||||
|
include 'constants.include.F'
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: n_grid2
|
||||||
|
double precision, intent(in) :: r1(3)
|
||||||
|
double precision, intent(out) :: res(n_grid2)
|
||||||
|
|
||||||
|
integer :: jpoint
|
||||||
|
double precision :: r2(3)
|
||||||
|
double precision :: mu_tmp, r12
|
||||||
|
|
||||||
|
PROVIDE final_grid_points_extra
|
||||||
|
|
||||||
|
if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
|
||||||
|
r2(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r2(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r2(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
|
||||||
|
+ (r1(2) - r2(2)) * (r1(2) - r2(2)) &
|
||||||
|
+ (r1(3) - r2(3)) * (r1(3) - r2(3)) )
|
||||||
|
mu_tmp = mu_erf * r12
|
||||||
|
|
||||||
|
res(jpoint) = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf
|
||||||
|
enddo
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu_r1_seq'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine j12_mu_r1_seq
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine j1b_nucl_r1_seq(n_grid2, res)
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer, intent(in) :: n_grid2
|
||||||
|
double precision, intent(out) :: res(n_grid2)
|
||||||
|
|
||||||
|
double precision :: r(3)
|
||||||
|
integer :: i, jpoint
|
||||||
|
double precision :: a, d, e, x, y, z
|
||||||
|
|
||||||
|
if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then
|
||||||
|
|
||||||
|
res = 1.d0
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
r(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
do i = 1, nucl_num
|
||||||
|
a = j1b_pen(i)
|
||||||
|
d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
|
||||||
|
+ (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
|
||||||
|
+ (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
|
||||||
|
|
||||||
|
res(jpoint) -= dexp(-a*dsqrt(d))
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then
|
||||||
|
|
||||||
|
res = 1.d0
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
r(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
do i = 1, nucl_num
|
||||||
|
a = j1b_pen(i)
|
||||||
|
d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
|
||||||
|
+ (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
|
||||||
|
+ (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
|
||||||
|
e = 1.d0 - dexp(-a*d)
|
||||||
|
|
||||||
|
res(jpoint) *= e
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then
|
||||||
|
|
||||||
|
res = 1.d0
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
r(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
do i = 1, nucl_num
|
||||||
|
a = j1b_pen(i)
|
||||||
|
d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
|
||||||
|
+ (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
|
||||||
|
+ (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
|
||||||
|
res(jpoint) -= j1b_pen_coef(i) * dexp(-a*d)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
|
||||||
|
|
||||||
|
res = 1.d0
|
||||||
|
|
||||||
|
do jpoint = 1, n_points_extra_final_grid ! r2
|
||||||
|
r(1) = final_grid_points_extra(1,jpoint)
|
||||||
|
r(2) = final_grid_points_extra(2,jpoint)
|
||||||
|
r(3) = final_grid_points_extra(3,jpoint)
|
||||||
|
|
||||||
|
do i = 1, nucl_num
|
||||||
|
a = j1b_pen(i)
|
||||||
|
x = r(1) - nucl_coord(i,1)
|
||||||
|
y = r(2) - nucl_coord(i,2)
|
||||||
|
z = r(3) - nucl_coord(i,3)
|
||||||
|
d = x*x + y*y + z*z
|
||||||
|
res(jpoint) -= dexp(-a*d*d)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_r1_seq'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine j1b_nucl_r1_seq
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|||||||
@ -149,22 +149,6 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
|
|||||||
deallocate(b_mat)
|
deallocate(b_mat)
|
||||||
|
|
||||||
call sum_A_At(tc_grad_and_lapl_ao(1,1,1,1), ao_num*ao_num)
|
call sum_A_At(tc_grad_and_lapl_ao(1,1,1,1), ao_num*ao_num)
|
||||||
! !$OMP PARALLEL &
|
|
||||||
! !$OMP DEFAULT (NONE) &
|
|
||||||
! !$OMP PRIVATE (i, j, k, l) &
|
|
||||||
! !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao, ao_num)
|
|
||||||
! !$OMP DO SCHEDULE (static)
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! do l = 1, ao_num
|
|
||||||
! do i = 1, ao_num
|
|
||||||
! do k = 1, ao_num
|
|
||||||
! tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! !$OMP END DO
|
|
||||||
! !$OMP END PARALLEL
|
|
||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
|
|||||||
@ -1,10 +1,11 @@
|
|||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)]
|
BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_final_grid, 3)]
|
||||||
|
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
!
|
!
|
||||||
|
! TODO
|
||||||
|
! combine with int2_grad1_u12_square_ao to avoid repeated calculation ?
|
||||||
|
!
|
||||||
! int2_grad1_u12_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2)
|
! int2_grad1_u12_ao(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2)
|
||||||
!
|
!
|
||||||
! where r1 = r(ipoint)
|
! where r1 = r(ipoint)
|
||||||
@ -106,7 +107,6 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
|
|||||||
elseif(j1b_type .ge. 100) then
|
elseif(j1b_type .ge. 100) then
|
||||||
|
|
||||||
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
|
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
|
||||||
PROVIDE grad1_u12_num
|
|
||||||
|
|
||||||
double precision, allocatable :: tmp(:,:,:)
|
double precision, allocatable :: tmp(:,:,:)
|
||||||
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
|
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
|
||||||
@ -126,39 +126,71 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
|
|||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
int2_grad1_u12_ao = 0.d0
|
integer :: n_blocks, n_rest, n_pass
|
||||||
do m = 1, 3
|
integer :: i_blocks, i_rest, i_pass, ii
|
||||||
!call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0 &
|
double precision, allocatable :: tmp_grad1_u12(:,:,:)
|
||||||
! this work also because of the symmetry in K(1,2) and sign compensation in L(1,2,3)
|
|
||||||
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -1.d0 &
|
! n_points_final_grid = n_blocks * n_pass + n_rest
|
||||||
, tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid &
|
n_blocks = 8
|
||||||
, 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num)
|
n_rest = int(mod(n_points_final_grid, n_blocks))
|
||||||
|
n_pass = int((n_points_final_grid - n_rest) / n_blocks)
|
||||||
|
|
||||||
|
if(n_pass .le. 1) then
|
||||||
|
print*, ' blocks are to large or grid is very small !'
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
|
||||||
|
allocate(tmp_grad1_u12(n_points_extra_final_grid,n_blocks,3))
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT (NONE) &
|
||||||
|
!$OMP PRIVATE (i_pass, i_blocks, ipoint, ii, m, tmp_grad1_u12) &
|
||||||
|
!$OMP SHARED (n_pass, n_blocks, n_points_extra_final_grid, ao_num, &
|
||||||
|
!$OMP final_grid_points, tmp, int2_grad1_u12_ao)
|
||||||
|
!$OMP DO
|
||||||
|
do i_pass = 1, n_pass
|
||||||
|
ii = (i_pass-1)*n_blocks + 1
|
||||||
|
|
||||||
|
do i_blocks = 1, n_blocks
|
||||||
|
ipoint = ii - 1 + i_blocks ! r1
|
||||||
|
call get_grad1_u12_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12(1,i_blocks,1) &
|
||||||
|
, tmp_grad1_u12(1,i_blocks,2) &
|
||||||
|
, tmp_grad1_u12(1,i_blocks,3) )
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
!! these dgemm are equivalent to
|
do m = 1, 3
|
||||||
!!$OMP PARALLEL &
|
call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, 1.d0 &
|
||||||
!!$OMP DEFAULT (NONE) &
|
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid &
|
||||||
!!$OMP PRIVATE (j, i, ipoint, jpoint, w) &
|
, 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num)
|
||||||
!!$OMP SHARED (int2_grad1_u12_ao, ao_num, n_points_final_grid, &
|
enddo
|
||||||
!!$OMP n_points_extra_final_grid, final_weight_at_r_vector_extra, &
|
enddo
|
||||||
!!$OMP aos_in_r_array_extra_transp, grad1_u12_num, tmp)
|
!$OMP END DO
|
||||||
!!$OMP DO SCHEDULE (static)
|
!$OMP END PARALLEL
|
||||||
!do ipoint = 1, n_points_final_grid
|
|
||||||
! do j = 1, ao_num
|
deallocate(tmp_grad1_u12)
|
||||||
! do i = 1, ao_num
|
|
||||||
! do jpoint = 1, n_points_extra_final_grid
|
! TODO
|
||||||
! w = -tmp(jpoint,i,j)
|
! OPENMP
|
||||||
! !w = tmp(jpoint,i,j) this work also because of the symmetry in K(1,2)
|
if(n_rest .ne. 0) then
|
||||||
! ! and sign compensation in L(1,2,3)
|
|
||||||
! int2_grad1_u12_ao(i,j,ipoint,1) += w * grad1_u12_num(jpoint,ipoint,1)
|
allocate(tmp_grad1_u12(n_points_extra_final_grid,n_rest,3))
|
||||||
! int2_grad1_u12_ao(i,j,ipoint,2) += w * grad1_u12_num(jpoint,ipoint,2)
|
|
||||||
! int2_grad1_u12_ao(i,j,ipoint,3) += w * grad1_u12_num(jpoint,ipoint,3)
|
ii = n_pass*n_blocks + 1
|
||||||
! enddo
|
do i_rest = 1, n_rest
|
||||||
! enddo
|
ipoint = ii - 1 + i_rest ! r1
|
||||||
! enddo
|
call get_grad1_u12_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12(1,i_rest,1) &
|
||||||
!enddo
|
, tmp_grad1_u12(1,i_rest,2) &
|
||||||
!!$OMP END DO
|
, tmp_grad1_u12(1,i_rest,3) )
|
||||||
!!$OMP END PARALLEL
|
enddo
|
||||||
|
|
||||||
|
do m = 1, 3
|
||||||
|
call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, 1.d0 &
|
||||||
|
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid &
|
||||||
|
, 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
deallocate(tmp_grad1_u12)
|
||||||
|
endif
|
||||||
|
|
||||||
deallocate(tmp)
|
deallocate(tmp)
|
||||||
else
|
else
|
||||||
@ -185,6 +217,72 @@ END_PROVIDER
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_num_1shot, (ao_num, ao_num, n_points_final_grid, 3)]
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! int2_grad1_u12_ao_num_1shot(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2)
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer :: ipoint, i, j, m, jpoint
|
||||||
|
double precision :: time0, time1
|
||||||
|
double precision :: x, y, z, w, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
|
||||||
|
|
||||||
|
print*, ' providing int2_grad1_u12_ao_num_1shot ...'
|
||||||
|
call wall_time(time0)
|
||||||
|
|
||||||
|
PROVIDE j1b_type
|
||||||
|
|
||||||
|
if(j1b_type .ge. 100) then
|
||||||
|
|
||||||
|
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
|
||||||
|
PROVIDE grad1_u12_num
|
||||||
|
|
||||||
|
double precision, allocatable :: tmp(:,:,:)
|
||||||
|
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
|
||||||
|
tmp = 0.d0
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT (NONE) &
|
||||||
|
!$OMP PRIVATE (j, i, jpoint) &
|
||||||
|
!$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
|
||||||
|
!$OMP DO SCHEDULE (static)
|
||||||
|
do j = 1, ao_num
|
||||||
|
do i = 1, ao_num
|
||||||
|
do jpoint = 1, n_points_extra_final_grid
|
||||||
|
tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
do m = 1, 3
|
||||||
|
!call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -1.d0 &
|
||||||
|
! this work also because of the symmetry in K(1,2) and sign compensation in L(1,2,3)
|
||||||
|
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0 &
|
||||||
|
, tmp(1,1,1), n_points_extra_final_grid, grad1_u12_num(1,1,m), n_points_extra_final_grid &
|
||||||
|
, 0.d0, int2_grad1_u12_ao_num_1shot(1,1,1,m), ao_num*ao_num)
|
||||||
|
enddo
|
||||||
|
|
||||||
|
deallocate(tmp)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
call wall_time(time1)
|
||||||
|
print*, ' wall time for int2_grad1_u12_ao_num_1shot =', time1-time0
|
||||||
|
call print_memory_usage()
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_points_final_grid)]
|
BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_points_final_grid)]
|
||||||
|
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
@ -275,16 +373,14 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
|
|||||||
elseif(j1b_type .ge. 100) then
|
elseif(j1b_type .ge. 100) then
|
||||||
|
|
||||||
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
|
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
|
||||||
PROVIDE grad1_u12_squared_num
|
|
||||||
|
|
||||||
double precision, allocatable :: tmp(:,:,:)
|
double precision, allocatable :: tmp(:,:,:)
|
||||||
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
|
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
|
||||||
tmp = 0.d0
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PRIVATE (j, i, jpoint) &
|
!$OMP PRIVATE (j, i, jpoint) &
|
||||||
!$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
|
!$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
|
||||||
!$OMP DO SCHEDULE (static)
|
!$OMP DO COLLAPSE(2)
|
||||||
do j = 1, ao_num
|
do j = 1, ao_num
|
||||||
do i = 1, ao_num
|
do i = 1, ao_num
|
||||||
do jpoint = 1, n_points_extra_final_grid
|
do jpoint = 1, n_points_extra_final_grid
|
||||||
@ -295,31 +391,63 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
|
|||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
int2_grad1_u12_square_ao = 0.d0
|
integer :: n_blocks, n_rest, n_pass
|
||||||
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 &
|
integer :: i_blocks, i_rest, i_pass, ii
|
||||||
, tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid &
|
double precision, allocatable :: tmp_grad1_u12_squared(:,:)
|
||||||
, 0.d0, int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num)
|
|
||||||
|
|
||||||
!! this dgemm is equivalen to
|
! n_points_final_grid = n_blocks * n_pass + n_rest
|
||||||
!!$OMP PARALLEL &
|
n_blocks = 16
|
||||||
!!$OMP DEFAULT (NONE) &
|
n_rest = int(mod(n_points_final_grid, n_blocks))
|
||||||
!!$OMP PRIVATE (i, j, ipoint, jpoint, w) &
|
n_pass = int((n_points_final_grid - n_rest) / n_blocks)
|
||||||
!!$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, &
|
|
||||||
!!$OMP n_points_extra_final_grid, final_weight_at_r_vector_extra, &
|
if(n_pass .le. 1) then
|
||||||
!!$OMP aos_in_r_array_extra_transp, grad1_u12_squared_num, tmp)
|
print*, ' blocks are to large or grid is very small !'
|
||||||
!!$OMP DO SCHEDULE (static)
|
stop
|
||||||
!do ipoint = 1, n_points_final_grid
|
endif
|
||||||
! do j = 1, ao_num
|
|
||||||
! do i = 1, ao_num
|
allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_blocks))
|
||||||
! do jpoint = 1, n_points_extra_final_grid
|
|
||||||
! w = -0.5d0 * tmp(jpoint,i,j)
|
!$OMP PARALLEL &
|
||||||
! int2_grad1_u12_square_ao(i,j,ipoint) += w * grad1_u12_squared_num(jpoint,ipoint)
|
!$OMP DEFAULT (NONE) &
|
||||||
! enddo
|
!$OMP PRIVATE (i_pass, i_blocks, ipoint, ii, tmp_grad1_u12_squared) &
|
||||||
! enddo
|
!$OMP SHARED (n_pass, n_blocks, n_points_extra_final_grid, ao_num, &
|
||||||
! enddo
|
!$OMP final_grid_points, tmp, int2_grad1_u12_square_ao)
|
||||||
!enddo
|
!$OMP DO
|
||||||
!!$OMP END DO
|
do i_pass = 1, n_pass
|
||||||
!!$OMP END PARALLEL
|
ii = (i_pass-1)*n_blocks + 1
|
||||||
|
|
||||||
|
do i_blocks = 1, n_blocks
|
||||||
|
ipoint = ii - 1 + i_blocks ! r1
|
||||||
|
call get_grad1_u12_squared_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12_squared(1,i_blocks))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, -0.5d0 &
|
||||||
|
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid &
|
||||||
|
, 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num)
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
deallocate(tmp_grad1_u12_squared)
|
||||||
|
|
||||||
|
! TODO
|
||||||
|
! OPENMP
|
||||||
|
if(n_rest .ne. 0) then
|
||||||
|
|
||||||
|
allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_rest))
|
||||||
|
|
||||||
|
ii = n_pass*n_blocks + 1
|
||||||
|
do i_rest = 1, n_rest
|
||||||
|
ipoint = ii - 1 + i_rest ! r1
|
||||||
|
call get_grad1_u12_squared_r1_seq(final_grid_points(1,ipoint), n_points_extra_final_grid, tmp_grad1_u12_squared(1,i_rest))
|
||||||
|
enddo
|
||||||
|
|
||||||
|
call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, -0.5d0 &
|
||||||
|
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid &
|
||||||
|
, 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num)
|
||||||
|
|
||||||
|
deallocate(tmp_grad1_u12_squared)
|
||||||
|
endif
|
||||||
|
|
||||||
deallocate(tmp)
|
deallocate(tmp)
|
||||||
|
|
||||||
@ -338,3 +466,65 @@ END_PROVIDER
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
|
BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao_num_1shot, (ao_num, ao_num, n_points_final_grid)]
|
||||||
|
|
||||||
|
BEGIN_DOC
|
||||||
|
!
|
||||||
|
! int2_grad1_u12_square_ao_num_1shot = -(1/2) x int dr2 chi_l(r2) chi_j(r2) [grad_1 u(r1,r2)]^2
|
||||||
|
!
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer :: i, j, jpoint
|
||||||
|
double precision :: time0, time1
|
||||||
|
|
||||||
|
print*, ' providing int2_grad1_u12_square_ao_num_1shot ...'
|
||||||
|
call wall_time(time0)
|
||||||
|
|
||||||
|
PROVIDE j1b_type
|
||||||
|
|
||||||
|
if(j1b_type .ge. 100) then
|
||||||
|
|
||||||
|
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
|
||||||
|
PROVIDE grad1_u12_squared_num
|
||||||
|
|
||||||
|
double precision, allocatable :: tmp(:,:,:)
|
||||||
|
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
|
||||||
|
|
||||||
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT (NONE) &
|
||||||
|
!$OMP PRIVATE (j, i, jpoint) &
|
||||||
|
!$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
|
||||||
|
!$OMP DO COLLAPSE(2)
|
||||||
|
do j = 1, ao_num
|
||||||
|
do i = 1, ao_num
|
||||||
|
do jpoint = 1, n_points_extra_final_grid
|
||||||
|
tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, -0.5d0 &
|
||||||
|
, tmp(1,1,1), n_points_extra_final_grid, grad1_u12_squared_num(1,1), n_points_extra_final_grid &
|
||||||
|
, 0.d0, int2_grad1_u12_square_ao_num_1shot(1,1,1), ao_num*ao_num)
|
||||||
|
|
||||||
|
FREE grad1_u12_squared_num
|
||||||
|
deallocate(tmp)
|
||||||
|
|
||||||
|
else
|
||||||
|
|
||||||
|
print *, ' j1b_type = ', j1b_type, 'not implemented yet'
|
||||||
|
stop
|
||||||
|
|
||||||
|
endif
|
||||||
|
|
||||||
|
call wall_time(time1)
|
||||||
|
print*, ' wall time for int2_grad1_u12_square_ao_num_1shot =', time1-time0
|
||||||
|
call print_memory_usage()
|
||||||
|
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|||||||
@ -11,12 +11,24 @@ program test_non_h
|
|||||||
my_n_pt_a_grid = tc_grid1_a
|
my_n_pt_a_grid = tc_grid1_a
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||||
|
|
||||||
|
if(j1b_type .ge. 100) then
|
||||||
|
my_extra_grid_becke = .True.
|
||||||
|
PROVIDE tc_grid2_a tc_grid2_r
|
||||||
|
my_n_pt_r_extra_grid = tc_grid2_r
|
||||||
|
my_n_pt_a_extra_grid = tc_grid2_a
|
||||||
|
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
|
||||||
|
endif
|
||||||
|
|
||||||
|
|
||||||
!call routine_grad_squared()
|
!call routine_grad_squared()
|
||||||
!call routine_fit()
|
!call routine_fit()
|
||||||
|
|
||||||
!call test_ipp()
|
!call test_ipp()
|
||||||
|
|
||||||
call test_v_ij_u_cst_mu_j1b_an()
|
!call test_v_ij_u_cst_mu_j1b_an()
|
||||||
|
|
||||||
|
call test_int2_grad1_u12_square_ao()
|
||||||
|
call test_int2_grad1_u12_ao()
|
||||||
end
|
end
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
@ -583,7 +595,91 @@ subroutine test_v_ij_u_cst_mu_j1b_an()
|
|||||||
print*, ' accuracy(%) = ', 100.d0 * accu / norm
|
print*, ' accuracy(%) = ', 100.d0 * accu / norm
|
||||||
|
|
||||||
return
|
return
|
||||||
end subroutine test_v_ij_u_cst_mu_j1b_an()
|
end subroutine test_v_ij_u_cst_mu_j1b_an
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine test_int2_grad1_u12_square_ao()
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer :: i, j, ipoint
|
||||||
|
double precision :: I_old, I_new
|
||||||
|
double precision :: norm, accu, thr, diff
|
||||||
|
|
||||||
|
PROVIDE int2_grad1_u12_square_ao
|
||||||
|
PROVIDE int2_grad1_u12_square_ao_num_1shot
|
||||||
|
|
||||||
|
thr = 1d-8
|
||||||
|
norm = 0.d0
|
||||||
|
accu = 0.d0
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
do i = 1, ao_num
|
||||||
|
do j = 1, ao_num
|
||||||
|
|
||||||
|
I_old = int2_grad1_u12_square_ao_num_1shot(j,i,ipoint)
|
||||||
|
I_new = int2_grad1_u12_square_ao (j,i,ipoint)
|
||||||
|
|
||||||
|
diff = dabs(I_new-I_old)
|
||||||
|
if(diff .gt. thr) then
|
||||||
|
print *, ' problem on:', j, i, ipoint
|
||||||
|
print *, ' old value :', I_old
|
||||||
|
print *, ' new value :', I_new
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
|
||||||
|
accu += diff
|
||||||
|
norm += dabs(I_old)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
print*, ' accuracy(%) = ', 100.d0 * accu / norm
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine test_int2_grad1_u12_square_ao
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
subroutine test_int2_grad1_u12_ao()
|
||||||
|
|
||||||
|
implicit none
|
||||||
|
integer :: i, j, ipoint, m
|
||||||
|
double precision :: I_old, I_new
|
||||||
|
double precision :: norm, accu, thr, diff
|
||||||
|
|
||||||
|
PROVIDE int2_grad1_u12_ao
|
||||||
|
PROVIDE int2_grad1_u12_ao_num_1shot
|
||||||
|
|
||||||
|
thr = 1d-8
|
||||||
|
norm = 0.d0
|
||||||
|
accu = 0.d0
|
||||||
|
do ipoint = 1, n_points_final_grid
|
||||||
|
do i = 1, ao_num
|
||||||
|
do j = 1, ao_num
|
||||||
|
|
||||||
|
do m = 1, 3
|
||||||
|
I_old = int2_grad1_u12_ao_num_1shot(j,i,ipoint,m)
|
||||||
|
I_new = int2_grad1_u12_ao (j,i,ipoint,m)
|
||||||
|
|
||||||
|
diff = dabs(I_new-I_old)
|
||||||
|
if(diff .gt. thr) then
|
||||||
|
print *, ' problem on:', j, i, ipoint, m
|
||||||
|
print *, ' old value :', I_old
|
||||||
|
print *, ' new value :', I_new
|
||||||
|
stop
|
||||||
|
endif
|
||||||
|
|
||||||
|
accu += diff
|
||||||
|
norm += dabs(I_old)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
|
||||||
|
print*, ' accuracy(%) = ', 100.d0 * accu / norm
|
||||||
|
|
||||||
|
return
|
||||||
|
end subroutine test_int2_grad1_u12_ao
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|||||||
@ -1,7 +1,4 @@
|
|||||||
|
|
||||||
! TODO
|
|
||||||
! remove ao_two_e_coul and use map directly
|
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, ao_num)]
|
BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, ao_num)]
|
||||||
@ -160,8 +157,6 @@ BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num
|
|||||||
END_DOC
|
END_DOC
|
||||||
|
|
||||||
integer :: i, j, k, l
|
integer :: i, j, k, l
|
||||||
double precision :: integral
|
|
||||||
double precision, allocatable :: tmp(:)
|
|
||||||
double precision, external :: get_ao_two_e_integral
|
double precision, external :: get_ao_two_e_integral
|
||||||
|
|
||||||
PROVIDE ao_integrals_map
|
PROVIDE ao_integrals_map
|
||||||
@ -183,25 +178,6 @@ BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num
|
|||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
|
|
||||||
! TODO
|
|
||||||
! allocate(tmp(ao_num))
|
|
||||||
!
|
|
||||||
! !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,l,j,k,tmp)
|
|
||||||
! do j = 1, ao_num
|
|
||||||
! do l = 1, ao_num
|
|
||||||
! do i = 1, ao_num
|
|
||||||
! call get_ao_two_e_integrals(i, l, l, ao_num, tmp(1))
|
|
||||||
! do k = 1, ao_num
|
|
||||||
! ao_two_e_coul(k,i,l,j) = tmp(k)
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! enddo
|
|
||||||
! !$OMP END PARALLEL DO
|
|
||||||
!
|
|
||||||
! deallocate(tmp)
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|||||||
@ -13,11 +13,9 @@ program tc_scf
|
|||||||
print *, ' starting ...'
|
print *, ' starting ...'
|
||||||
|
|
||||||
my_grid_becke = .True.
|
my_grid_becke = .True.
|
||||||
|
|
||||||
PROVIDE tc_grid1_a tc_grid1_r
|
PROVIDE tc_grid1_a tc_grid1_r
|
||||||
my_n_pt_r_grid = tc_grid1_r
|
my_n_pt_r_grid = tc_grid1_r
|
||||||
my_n_pt_a_grid = tc_grid1_a
|
my_n_pt_a_grid = tc_grid1_a
|
||||||
|
|
||||||
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
|
||||||
|
|
||||||
PROVIDE mu_erf
|
PROVIDE mu_erf
|
||||||
@ -26,6 +24,14 @@ program tc_scf
|
|||||||
print *, ' j1b_type = ', j1b_type
|
print *, ' j1b_type = ', j1b_type
|
||||||
print *, j1b_pen
|
print *, j1b_pen
|
||||||
|
|
||||||
|
if(j1b_type .ge. 100) then
|
||||||
|
my_extra_grid_becke = .True.
|
||||||
|
PROVIDE tc_grid2_a tc_grid2_r
|
||||||
|
my_n_pt_r_extra_grid = tc_grid2_r
|
||||||
|
my_n_pt_a_extra_grid = tc_grid2_a
|
||||||
|
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
|
||||||
|
endif
|
||||||
|
|
||||||
!call create_guess()
|
!call create_guess()
|
||||||
!call orthonormalize_mos()
|
!call orthonormalize_mos()
|
||||||
|
|
||||||
|
|||||||
Loading…
Reference in New Issue
Block a user