LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 19:43:32 +01:00
Code Releases Activity

working on aos debug

Browse Source
This commit is contained in:
Abdallah Ammar 2024-05-07 01:56:14 +02:00
parent bd8d45b99b
commit 2a8b9e544b
2 changed files with 285 additions and 239 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

16
plugins/local/non_h_ints_mu/deb_aos.irp.f
View File

@ -31,12 +31,14 @@ subroutine print_aos()
integer :: i, ipoint integer :: i, ipoint
double precision :: r(3) double precision :: r(3)
double precision :: ao_val, ao_der(3), ao_lap double precision :: ao_val, ao_der(3), ao_lap
double precision :: mo_val, mo_der(3), mo_lap
PROVIDE final_grid_points aos_in_r_array aos_grad_in_r_array aos_lapl_in_r_array PROVIDE final_grid_points aos_in_r_array aos_grad_in_r_array aos_lapl_in_r_array
write(1000, *) n_points_final_grid
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
r(:) = final_grid_points(:,ipoint) r(:) = final_grid_points(:,ipoint)
print*, r write(1000, '(3(f15.7, 3X))') r
enddo enddo
do ipoint = 1, n_points_final_grid do ipoint = 1, n_points_final_grid
@ -45,7 +47,17 @@ subroutine print_aos()
ao_val = aos_in_r_array (i,ipoint) ao_val = aos_in_r_array (i,ipoint)
ao_der(:) = aos_grad_in_r_array(i,ipoint,:) ao_der(:) = aos_grad_in_r_array(i,ipoint,:)
ao_lap = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint) ao_lap = aos_lapl_in_r_array(1,i,ipoint) + aos_lapl_in_r_array(2,i,ipoint) + aos_lapl_in_r_array(3,i,ipoint)
write(*, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap write(1010, '(5(f15.7, 3X))') ao_val, ao_der, ao_lap
enddo
enddo
do ipoint = 1, n_points_final_grid
r(:) = final_grid_points(:,ipoint)
do i = 1, mo_num
mo_val = mos_in_r_array (i,ipoint)
mo_der(:) = mos_grad_in_r_array(i,ipoint,:)
mo_lap = mos_lapl_in_r_array(i,ipoint,1) + mos_lapl_in_r_array(i,ipoint,2) + mos_lapl_in_r_array(i,ipoint,3)
write(2010, '(5(f15.7, 3X))') mo_val, mo_der, mo_lap
enddo enddo
enddo enddo

134
src/ao_basis/aos_in_r.irp.f
View File

@ -1,22 +1,28 @@
! ---
double precision function ao_value(i, r) double precision function ao_value(i, r)
implicit none
BEGIN_DOC BEGIN_DOC
! Returns the value of the i-th ao at point $\textbf{r}$ ! Returns the value of the i-th ao at point $\textbf{r}$
END_DOC END_DOC
double precision, intent(in) :: r(3)
implicit none
integer, intent(in) :: i integer, intent(in) :: i
double precision, intent(in) :: r(3)
integer :: m, num_ao integer :: m, num_ao
integer :: power_ao(3)
double precision :: center_ao(3) double precision :: center_ao(3)
double precision :: beta double precision :: beta
integer :: power_ao(3)
double precision :: accu, dx, dy, dz, r2 double precision :: accu, dx, dy, dz, r2
num_ao = ao_nucl(i) num_ao = ao_nucl(i)
power_ao(1:3) = ao_power(i,1:3) power_ao(1:3) = ao_power(i,1:3)
center_ao(1:3) = nucl_coord(num_ao,1:3) center_ao(1:3) = nucl_coord(num_ao,1:3)
dx = (r(1) - center_ao(1)) dx = r(1) - center_ao(1)
dy = (r(2) - center_ao(2)) dy = r(2) - center_ao(2)
dz = (r(3) - center_ao(3)) dz = r(3) - center_ao(3)
r2 = dx*dx + dy*dy + dz*dz r2 = dx*dx + dy*dy + dz*dz
dx = dx**power_ao(1) dx = dx**power_ao(1)
dy = dy**power_ao(2) dy = dy**power_ao(2)
@ -33,25 +39,28 @@ end
double precision function primitive_value(i, j, r) double precision function primitive_value(i, j, r)
implicit none
BEGIN_DOC BEGIN_DOC
! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r} ! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r}
! **without the coefficient** ! **without the coefficient**
END_DOC END_DOC
double precision, intent(in) :: r(3)
implicit none
integer, intent(in) :: i, j integer, intent(in) :: i, j
double precision, intent(in) :: r(3)
integer :: m, num_ao integer :: m, num_ao
integer :: power_ao(3)
double precision :: center_ao(3) double precision :: center_ao(3)
double precision :: beta double precision :: beta
integer :: power_ao(3)
double precision :: accu, dx, dy, dz, r2 double precision :: accu, dx, dy, dz, r2
num_ao = ao_nucl(i) num_ao = ao_nucl(i)
power_ao(1:3)= ao_power(i,1:3) power_ao(1:3)= ao_power(i,1:3)
center_ao(1:3) = nucl_coord(num_ao,1:3) center_ao(1:3) = nucl_coord(num_ao,1:3)
dx = (r(1) - center_ao(1)) dx = r(1) - center_ao(1)
dy = (r(2) - center_ao(2)) dy = r(2) - center_ao(2)
dz = (r(3) - center_ao(3)) dz = r(3) - center_ao(3)
r2 = dx*dx + dy*dy + dz*dz r2 = dx*dx + dy*dy + dz*dz
dx = dx**power_ao(1) dx = dx**power_ao(1)
dy = dy**power_ao(2) dy = dy**power_ao(2)
@ -106,7 +115,7 @@ subroutine give_all_aos_at_r(r, tmp_array)
tmp_array(k) = 0.d0 tmp_array(k) = 0.d0
do l = 1, ao_prim_num(k) do l = 1, ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i) beta = ao_expo_ordered_transp_per_nucl(l,j,i)
if(dabs(beta*r2).gt.40.d0) cycle if(beta*r2.gt.50.d0) cycle
tmp_array(k) += ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) tmp_array(k) += ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
enddo enddo
@ -121,8 +130,9 @@ end
! --- ! ---
subroutine give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array) subroutine give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array)
implicit none
BEGIN_DOC BEGIN_DOC
!
! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z ! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
! !
! output : ! output :
@ -131,6 +141,8 @@ subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$ ! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$
! !
END_DOC END_DOC
implicit none
double precision, intent(in) :: r(3) double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num) double precision, intent(out) :: aos_array(ao_num)
double precision, intent(out) :: aos_grad_array(3,ao_num) double precision, intent(out) :: aos_grad_array(3,ao_num)
@ -138,70 +150,84 @@ subroutine give_all_aos_and_grad_at_r(r,aos_array,aos_grad_array)
integer :: power_ao(3) integer :: power_ao(3)
integer :: i, j, k, l, m integer :: i, j, k, l, m
double precision :: dx, dy, dz, r2 double precision :: dx, dy, dz, r2
double precision :: dx2,dy2,dz2
double precision :: dx1, dy1, dz1 double precision :: dx1, dy1, dz1
double precision :: dx2, dy2, dz2
double precision :: center_ao(3) double precision :: center_ao(3)
double precision :: beta, accu_1, accu_2, contrib double precision :: beta, accu_1, accu_2, contrib
do i = 1, nucl_num do i = 1, nucl_num
center_ao(1:3) = nucl_coord(i,1:3) center_ao(1:3) = nucl_coord(i,1:3)
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2)) dx = r(1) - center_ao(1)
dz = (r(3) - center_ao(3)) dy = r(2) - center_ao(2)
dz = r(3) - center_ao(3)
r2 = dx*dx + dy*dy + dz*dz r2 = dx*dx + dy*dy + dz*dz
do j = 1, Nucl_N_Aos(i) do j = 1, Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
aos_array(k) = 0.d0 aos_array(k) = 0.d0
aos_grad_array(1,k) = 0.d0 aos_grad_array(1,k) = 0.d0
aos_grad_array(2,k) = 0.d0 aos_grad_array(2,k) = 0.d0
aos_grad_array(3,k) = 0.d0 aos_grad_array(3,k) = 0.d0
power_ao(1:3) = ao_power_ordered_transp_per_nucl(1:3,j,i) power_ao(1:3) = ao_power_ordered_transp_per_nucl(1:3,j,i)
dx2 = dx**power_ao(1) dx2 = dx**power_ao(1)
dy2 = dy**power_ao(2) dy2 = dy**power_ao(2)
dz2 = dz**power_ao(3) dz2 = dz**power_ao(3)
dx1 = 0.d0
if(power_ao(1) .ne. 0) then if(power_ao(1) .ne. 0) then
dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1) dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
else
dx1 = 0.d0
endif endif
dy1 = 0.d0
if(power_ao(2) .ne. 0) then if(power_ao(2) .ne. 0) then
dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1) dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
else
dy1 = 0.d0
endif endif
dz1 = 0.d0
if(power_ao(3) .ne. 0) then if(power_ao(3) .ne. 0) then
dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1) dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
else
dz1 = 0.d0
endif endif
accu_1 = 0.d0 accu_1 = 0.d0
accu_2 = 0.d0 accu_2 = 0.d0
do l = 1, ao_prim_num(k) do l = 1, ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i) beta = ao_expo_ordered_transp_per_nucl(l,j,i)
contrib = 0.d0
if(beta*r2.gt.50.d0) cycle if(beta*r2.gt.50.d0) cycle
contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
accu_1 += contrib accu_1 += contrib
accu_2 += contrib * beta accu_2 += contrib * beta
enddo enddo
aos_array(k) = accu_1 * dx2 * dy2 * dz2 aos_array(k) = accu_1 * dx2 * dy2 * dz2
aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx * dy2 * dz2 * accu_2 aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx * dy2 * dz2 * accu_2
aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy * dz2 * accu_2 aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy * dz2 * accu_2
aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz * accu_2 aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz * accu_2
enddo enddo
enddo enddo
end end
! ---
subroutine give_all_aos_and_grad_and_lapl_at_r(r, aos_array, aos_grad_array, aos_lapl_array) subroutine give_all_aos_and_grad_and_lapl_at_r(r, aos_array, aos_grad_array, aos_lapl_array)
implicit none
BEGIN_DOC BEGIN_DOC
!
! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z ! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
! !
! output : ! output :
! !
! * aos_array(i) = ao(i) evaluated at $\textbf{r}$ ! * aos_array(i) = ao(i) evaluated at $\textbf{r}$
! * aos_grad_array(1,i) = $\nabla_x$ of the ao(i) evaluated at $\textbf{r}$ ! * aos_grad_array(1,i) = $\nabla_x$ of the ao(i) evaluated at $\textbf{r}$
!
END_DOC END_DOC
implicit none
double precision, intent(in) :: r(3) double precision, intent(in) :: r(3)
double precision, intent(out) :: aos_array(ao_num) double precision, intent(out) :: aos_array(ao_num)
double precision, intent(out) :: aos_grad_array(3,ao_num) double precision, intent(out) :: aos_grad_array(3,ao_num)
@ -210,26 +236,31 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
integer :: power_ao(3) integer :: power_ao(3)
integer :: i, j, k, l, m integer :: i, j, k, l, m
double precision :: dx, dy, dz, r2 double precision :: dx, dy, dz, r2
double precision :: dx2,dy2,dz2
double precision :: dx1, dy1, dz1 double precision :: dx1, dy1, dz1
double precision :: dx2, dy2, dz2
double precision :: dx3, dy3, dz3 double precision :: dx3, dy3, dz3
double precision :: dx4, dy4, dz4 double precision :: dx4, dy4, dz4
double precision :: dx5, dy5, dz5 double precision :: dx5, dy5, dz5
double precision :: center_ao(3) double precision :: center_ao(3)
double precision :: beta, accu_1, accu_2, accu_3, contrib double precision :: beta, accu_1, accu_2, accu_3, contrib
do i = 1, nucl_num do i = 1, nucl_num
center_ao(1:3) = nucl_coord(i,1:3) center_ao(1:3) = nucl_coord(i,1:3)
dx = (r(1) - center_ao(1))
dy = (r(2) - center_ao(2)) dx = r(1) - center_ao(1)
dz = (r(3) - center_ao(3)) dy = r(2) - center_ao(2)
dz = r(3) - center_ao(3)
r2 = dx*dx + dy*dy + dz*dz r2 = dx*dx + dy*dy + dz*dz
do j = 1, Nucl_N_Aos(i) do j = 1, Nucl_N_Aos(i)
k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
aos_array(k) = 0.d0 aos_array(k) = 0.d0
aos_grad_array(1,k) = 0.d0 aos_grad_array(1,k) = 0.d0
aos_grad_array(2,k) = 0.d0 aos_grad_array(2,k) = 0.d0
aos_grad_array(3,k) = 0.d0 aos_grad_array(3,k) = 0.d0
aos_lapl_array(1,k) = 0.d0 aos_lapl_array(1,k) = 0.d0
aos_lapl_array(2,k) = 0.d0 aos_lapl_array(2,k) = 0.d0
aos_lapl_array(3,k) = 0.d0 aos_lapl_array(3,k) = 0.d0
@ -238,17 +269,19 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
dx2 = dx**power_ao(1) dx2 = dx**power_ao(1)
dy2 = dy**power_ao(2) dy2 = dy**power_ao(2)
dz2 = dz**power_ao(3) dz2 = dz**power_ao(3)
! ---
dx1 = 0.d0
if(power_ao(1) .ne. 0) then if(power_ao(1) .ne. 0) then
dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1) dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
else
dx1 = 0.d0
endif endif
! For the Laplacian
dx3 = 0.d0
if(power_ao(1) .ge. 2) then if(power_ao(1) .ge. 2) then
dx3 = dble(power_ao(1)) * dble((power_ao(1)-1)) * dx**(power_ao(1)-2) dx3 = dble(power_ao(1)) * dble((power_ao(1)-1)) * dx**(power_ao(1)-2)
else
dx3 = 0.d0
endif endif
if(power_ao(1) .ge. 1) then if(power_ao(1) .ge. 1) then
dx4 = dble((2 * power_ao(1) + 1)) * dx**(power_ao(1)) dx4 = dble((2 * power_ao(1) + 1)) * dx**(power_ao(1))
else else
@ -257,16 +290,16 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
dx5 = dx**(power_ao(1)+2) dx5 = dx**(power_ao(1)+2)
! ---
dy1 = 0.d0
if(power_ao(2) .ne. 0) then if(power_ao(2) .ne. 0) then
dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1) dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
else
dy1 = 0.d0
endif endif
! For the Laplacian
dy3 = 0.d0
if(power_ao(2) .ge. 2) then if(power_ao(2) .ge. 2) then
dy3 = dble(power_ao(2)) * dble((power_ao(2)-1)) * dy**(power_ao(2)-2) dy3 = dble(power_ao(2)) * dble((power_ao(2)-1)) * dy**(power_ao(2)-2)
else
dy3 = 0.d0
endif endif
if(power_ao(2) .ge. 1) then if(power_ao(2) .ge. 1) then
@ -277,17 +310,16 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
dy5 = dy**(power_ao(2)+2) dy5 = dy**(power_ao(2)+2)
! ---
dz1 = 0.d0
if(power_ao(3) .ne. 0) then if(power_ao(3) .ne. 0) then
dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1) dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
else
dz1 = 0.d0
endif endif
! For the Laplacian
dz3 = 0.d0
if(power_ao(3) .ge. 2) then if(power_ao(3) .ge. 2) then
dz3 = dble(power_ao(3)) * dble((power_ao(3)-1)) * dz**(power_ao(3)-2) dz3 = dble(power_ao(3)) * dble((power_ao(3)-1)) * dz**(power_ao(3)-2)
else
dz3 = 0.d0
endif endif
if(power_ao(3) .ge. 1) then if(power_ao(3) .ge. 1) then
@ -298,29 +330,31 @@ subroutine give_all_aos_and_grad_and_lapl_at_r(r,aos_array,aos_grad_array,aos_la
dz5 = dz**(power_ao(3)+2) dz5 = dz**(power_ao(3)+2)
! ---
accu_1 = 0.d0 accu_1 = 0.d0
accu_2 = 0.d0 accu_2 = 0.d0
accu_3 = 0.d0 accu_3 = 0.d0
do l = 1,ao_prim_num(k) do l = 1,ao_prim_num(k)
beta = ao_expo_ordered_transp_per_nucl(l,j,i) beta = ao_expo_ordered_transp_per_nucl(l,j,i)
if(beta*r2.gt.50.d0) cycle
contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2) contrib = ao_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
accu_1 += contrib accu_1 += contrib
accu_2 += contrib * beta accu_2 += contrib * beta
accu_3 += contrib * beta**2 accu_3 += contrib * beta**2
enddo enddo
aos_array(k) = accu_1 * dx2 * dy2 * dz2
aos_array(k) = accu_1 * dx2 * dy2 * dz2
aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx * dy2 * dz2 * accu_2 aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx * dy2 * dz2 * accu_2
aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy * dz2 * accu_2 aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy * dz2 * accu_2
aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz * accu_2 aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz * accu_2
aos_lapl_array(1,k) = accu_1 * dx3 * dy2 * dz2 - 2.d0 * dx4 * dy2 * dz2 * accu_2 + 4.d0 * dx5 * dy2 * dz2 * accu_3 aos_lapl_array(1,k) = accu_1 * dx3 * dy2 * dz2 - 2.d0 * dx4 * dy2 * dz2 * accu_2 + 4.d0 * dx5 * dy2 * dz2 * accu_3
aos_lapl_array(2,k) = accu_1 * dx2 * dy3 * dz2 - 2.d0 * dx2 * dy4 * dz2 * accu_2 + 4.d0 * dx2 * dy5 * dz2 * accu_3 aos_lapl_array(2,k) = accu_1 * dx2 * dy3 * dz2 - 2.d0 * dx2 * dy4 * dz2 * accu_2 + 4.d0 * dx2 * dy5 * dz2 * accu_3
aos_lapl_array(3,k) = accu_1 * dx2 * dy2 * dz3 - 2.d0 * dx2 * dy2 * dz4 * accu_2 + 4.d0 * dx2 * dy2 * dz5 * accu_3 aos_lapl_array(3,k) = accu_1 * dx2 * dy2 * dz3 - 2.d0 * dx2 * dy2 * dz4 * accu_2 + 4.d0 * dx2 * dy2 * dz5 * accu_3
enddo
enddo
enddo
enddo
end end
! ---