From 9b2ba694d9e7f71801c5dac7c8073c06d8605b47 Mon Sep 17 00:00:00 2001
From: AbdAmmar
Date: Wed, 24 Jan 2024 19:25:17 +0100
Subject: [PATCH] Improved AosxAos representations of 1e-Jastrow
---
plugins/local/jastrow/README.md | 5 +-
plugins/local/non_h_ints_mu/jast_1e.irp.f | 195 +++++++++---------
.../local/non_h_ints_mu/jast_1e_utils.irp.f | 93 ++++-----
3 files changed, 149 insertions(+), 144 deletions(-)
diff --git a/plugins/local/jastrow/README.md b/plugins/local/jastrow/README.md
index 22486edd..67898e23 100644
--- a/plugins/local/jastrow/README.md
+++ b/plugins/local/jastrow/README.md
@@ -65,5 +65,8 @@ are defined by the tables `j1e_coef` and `j1e_expo`, respectively.
-- if `j1e_type` is **Charge_Harmonizer_AO**: The one-electron Jastrow factor **Charge_Harmonizer** is fitted by the atomic orbitals
+- if `j1e_type` is **Charge_Harmonizer_AO**: The one-electron Jastrow factor **Charge_Harmonizer** is fitted by the product of atomic orbitals:
+
+
+
diff --git a/plugins/local/non_h_ints_mu/jast_1e.irp.f b/plugins/local/non_h_ints_mu/jast_1e.irp.f
index 37ac0092..fbd032ed 100644
--- a/plugins/local/non_h_ints_mu/jast_1e.irp.f
+++ b/plugins/local/non_h_ints_mu/jast_1e.irp.f
@@ -231,96 +231,22 @@ END_PROVIDER
! !$OMP END PARALLEL
!
! deallocate(coef_fit)
-!
-! elseif(j1e_type .eq. "Charge_Harmonizer_AO2") then
-!
-! ! \grad_1 \sum_{\eta,\beta} C_{\eta,\beta} \chi_{\eta} \chi_{\beta}
-! ! where
-! ! \chi_{\eta} are the AOs
-! ! C_{\eta,\beta} are fitted to mimic (j1e_type .eq. "Charge_Harmonizer")
-! !
-! ! The - sign is in the parameters C_{\eta,\beta}
-!
-! PROVIDE aos_grad_in_r_array
-!
-! allocate(coef_fit2(ao_num*ao_num))
-!
-! if(mpi_master) then
-! call ezfio_has_jastrow_j1e_coef_ao2(exists)
-! endif
-! IRP_IF MPI_DEBUG
-! print *, irp_here, mpi_rank
-! call MPI_BARRIER(MPI_COMM_WORLD, ierr)
-! IRP_ENDIF
-! IRP_IF MPI
-! call MPI_BCAST(coef_fit2, ao_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
-! if (ierr /= MPI_SUCCESS) then
-! stop 'Unable to read j1e_coef_ao2 with MPI'
-! endif
-! IRP_ENDIF
-! if(exists) then
-! if(mpi_master) then
-! write(6,'(A)') '.. >>>>> [ IO READ: j1e_coef_ao2 ] <<<<< ..'
-! call ezfio_get_jastrow_j1e_coef_ao2(coef_fit2)
-! IRP_IF MPI
-! call MPI_BCAST(coef_fit2, ao_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
-! if (ierr /= MPI_SUCCESS) then
-! stop 'Unable to read j1e_coef_ao2 with MPI'
-! endif
-! IRP_ENDIF
-! endif
-! else
-!
-! call get_j1e_coef_fit_ao2(ao_num*ao_num, coef_fit2)
-! call ezfio_set_jastrow_j1e_coef_ao2(coef_fit2)
-!
-! endif
-!
-! !$OMP PARALLEL &
-! !$OMP DEFAULT (NONE) &
-! !$OMP PRIVATE (i, j, ij, ipoint, c) &
-! !$OMP SHARED (n_points_final_grid, ao_num, &
-! !$OMP aos_grad_in_r_array, coef_fit2, &
-! !$OMP aos_in_r_array, j1e_gradx, j1e_grady, j1e_gradz)
-! !$OMP DO SCHEDULE (static)
-! do ipoint = 1, n_points_final_grid
-!
-! j1e_gradx(ipoint) = 0.d0
-! j1e_grady(ipoint) = 0.d0
-! j1e_gradz(ipoint) = 0.d0
-!
-! do i = 1, ao_num
-! do j = 1, ao_num
-! ij = (i-1)*ao_num + j
-!
-! c = coef_fit2(ij)
-!
-! j1e_gradx(ipoint) += c * (aos_in_r_array(i,ipoint) * aos_grad_in_r_array(j,ipoint,1) + aos_grad_in_r_array(i,ipoint,1) * aos_in_r_array(j,ipoint))
-! j1e_grady(ipoint) += c * (aos_in_r_array(i,ipoint) * aos_grad_in_r_array(j,ipoint,2) + aos_grad_in_r_array(i,ipoint,2) * aos_in_r_array(j,ipoint))
-! j1e_gradz(ipoint) += c * (aos_in_r_array(i,ipoint) * aos_grad_in_r_array(j,ipoint,3) + aos_grad_in_r_array(i,ipoint,3) * aos_in_r_array(j,ipoint))
-! enddo
-! enddo
-! enddo
-! !$OMP END DO
-! !$OMP END PARALLEL
-!
-! deallocate(coef_fit2)
elseif(j1e_type .eq. "Charge_Harmonizer_AO") then
- ! \sum_{\eta} \vec{C}_{\eta} \chi_{\eta}
+ ! \grad_1 \sum_{\eta,\beta} C_{\eta,\beta} \chi_{\eta} \chi_{\beta}
! where
! \chi_{\eta} are the AOs
- ! \vec{C}_{\eta} are fitted to mimic (j1e_type .eq. "Charge_Harmonizer")
+ ! C_{\eta,\beta} are fitted to mimic (j1e_type .eq. "Charge_Harmonizer")
!
- ! The - sign is in the parameters \vec{C}_{\eta}
+ ! The - sign is in the parameters C_{\eta,\beta}
PROVIDE aos_grad_in_r_array
- allocate(coef_fit3(ao_num,3))
+ allocate(coef_fit2(ao_num*ao_num))
if(mpi_master) then
- call ezfio_has_jastrow_j1e_coef_ao3(exists)
+ call ezfio_has_jastrow_j1e_coef_ao2(exists)
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
@@ -328,34 +254,34 @@ END_PROVIDER
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
- call MPI_BCAST(coef_fit3, (ao_num*3), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+ call MPI_BCAST(coef_fit2, ao_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
- stop 'Unable to read j1e_coef_ao3 with MPI'
+ stop 'Unable to read j1e_coef_ao2 with MPI'
endif
IRP_ENDIF
if(exists) then
if(mpi_master) then
- write(6,'(A)') '.. >>>>> [ IO READ: j1e_coef_ao3 ] <<<<< ..'
- call ezfio_get_jastrow_j1e_coef_ao3(coef_fit3)
+ write(6,'(A)') '.. >>>>> [ IO READ: j1e_coef_ao2 ] <<<<< ..'
+ call ezfio_get_jastrow_j1e_coef_ao2(coef_fit2)
IRP_IF MPI
- call MPI_BCAST(coef_fit3, (ao_num*3), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+ call MPI_BCAST(coef_fit2, ao_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
- stop 'Unable to read j1e_coef_ao3 with MPI'
+ stop 'Unable to read j1e_coef_ao2 with MPI'
endif
IRP_ENDIF
endif
else
- call get_j1e_coef_fit_ao3(ao_num, coef_fit3)
- call ezfio_set_jastrow_j1e_coef_ao3(coef_fit3)
+ call get_j1e_coef_fit_ao2(ao_num*ao_num, coef_fit2)
+ call ezfio_set_jastrow_j1e_coef_ao2(coef_fit2)
endif
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
- !$OMP PRIVATE (i, ipoint, cx, cy, cz) &
+ !$OMP PRIVATE (i, j, ij, ipoint, c) &
!$OMP SHARED (n_points_final_grid, ao_num, &
- !$OMP aos_grad_in_r_array, coef_fit3, &
+ !$OMP aos_grad_in_r_array, coef_fit2, &
!$OMP aos_in_r_array, j1e_gradx, j1e_grady, j1e_gradz)
!$OMP DO SCHEDULE (static)
do ipoint = 1, n_points_final_grid
@@ -363,20 +289,95 @@ END_PROVIDER
j1e_gradx(ipoint) = 0.d0
j1e_grady(ipoint) = 0.d0
j1e_gradz(ipoint) = 0.d0
- do i = 1, ao_num
- cx = coef_fit3(i,1)
- cy = coef_fit3(i,2)
- cz = coef_fit3(i,3)
- j1e_gradx(ipoint) += cx * aos_in_r_array(i,ipoint)
- j1e_grady(ipoint) += cy * aos_in_r_array(i,ipoint)
- j1e_gradz(ipoint) += cz * aos_in_r_array(i,ipoint)
+ do i = 1, ao_num
+ do j = 1, ao_num
+ ij = (i-1)*ao_num + j
+
+ c = coef_fit2(ij)
+
+ j1e_gradx(ipoint) += c * (aos_in_r_array(i,ipoint) * aos_grad_in_r_array(j,ipoint,1) + aos_grad_in_r_array(i,ipoint,1) * aos_in_r_array(j,ipoint))
+ j1e_grady(ipoint) += c * (aos_in_r_array(i,ipoint) * aos_grad_in_r_array(j,ipoint,2) + aos_grad_in_r_array(i,ipoint,2) * aos_in_r_array(j,ipoint))
+ j1e_gradz(ipoint) += c * (aos_in_r_array(i,ipoint) * aos_grad_in_r_array(j,ipoint,3) + aos_grad_in_r_array(i,ipoint,3) * aos_in_r_array(j,ipoint))
+ enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
- deallocate(coef_fit3)
+ deallocate(coef_fit2)
+
+! elseif(j1e_type .eq. "Charge_Harmonizer_AO3") then
+!
+! ! \sum_{\eta} \vec{C}_{\eta} \chi_{\eta}
+! ! where
+! ! \chi_{\eta} are the AOs
+! ! \vec{C}_{\eta} are fitted to mimic (j1e_type .eq. "Charge_Harmonizer")
+! !
+! ! The - sign is in the parameters \vec{C}_{\eta}
+!
+! PROVIDE aos_grad_in_r_array
+!
+! allocate(coef_fit3(ao_num,3))
+!
+! if(mpi_master) then
+! call ezfio_has_jastrow_j1e_coef_ao3(exists)
+! endif
+! IRP_IF MPI_DEBUG
+! print *, irp_here, mpi_rank
+! call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+! IRP_ENDIF
+! IRP_IF MPI
+! !include 'mpif.h'
+! call MPI_BCAST(coef_fit3, (ao_num*3), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+! if (ierr /= MPI_SUCCESS) then
+! stop 'Unable to read j1e_coef_ao3 with MPI'
+! endif
+! IRP_ENDIF
+! if(exists) then
+! if(mpi_master) then
+! write(6,'(A)') '.. >>>>> [ IO READ: j1e_coef_ao3 ] <<<<< ..'
+! call ezfio_get_jastrow_j1e_coef_ao3(coef_fit3)
+! IRP_IF MPI
+! call MPI_BCAST(coef_fit3, (ao_num*3), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+! if (ierr /= MPI_SUCCESS) then
+! stop 'Unable to read j1e_coef_ao3 with MPI'
+! endif
+! IRP_ENDIF
+! endif
+! else
+!
+! call get_j1e_coef_fit_ao3(ao_num, coef_fit3)
+! call ezfio_set_jastrow_j1e_coef_ao3(coef_fit3)
+!
+! endif
+!
+! !$OMP PARALLEL &
+! !$OMP DEFAULT (NONE) &
+! !$OMP PRIVATE (i, ipoint, cx, cy, cz) &
+! !$OMP SHARED (n_points_final_grid, ao_num, &
+! !$OMP aos_grad_in_r_array, coef_fit3, &
+! !$OMP aos_in_r_array, j1e_gradx, j1e_grady, j1e_gradz)
+! !$OMP DO SCHEDULE (static)
+! do ipoint = 1, n_points_final_grid
+!
+! j1e_gradx(ipoint) = 0.d0
+! j1e_grady(ipoint) = 0.d0
+! j1e_gradz(ipoint) = 0.d0
+! do i = 1, ao_num
+! cx = coef_fit3(i,1)
+! cy = coef_fit3(i,2)
+! cz = coef_fit3(i,3)
+!
+! j1e_gradx(ipoint) += cx * aos_in_r_array(i,ipoint)
+! j1e_grady(ipoint) += cy * aos_in_r_array(i,ipoint)
+! j1e_gradz(ipoint) += cz * aos_in_r_array(i,ipoint)
+! enddo
+! enddo
+! !$OMP END DO
+! !$OMP END PARALLEL
+!
+! deallocate(coef_fit3)
else
diff --git a/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f b/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f
index 9dc0d5b0..842908a7 100644
--- a/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f
+++ b/plugins/local/non_h_ints_mu/jast_1e_utils.irp.f
@@ -80,24 +80,33 @@ subroutine get_j1e_coef_fit_ao(dim_fit, coef_fit)
allocate(A_inv(ao_num,ao_num))
call get_inverse(A, ao_num, ao_num, A_inv, ao_num)
- deallocate(A)
! coef_fit = A_inv x b
call dgemv("N", ao_num, ao_num, 1.d0, A_inv, ao_num, b, 1, 0.d0, coef_fit, 1)
- !integer :: j, k
- !double precision :: tmp
- !print *, ' check A_inv'
- !do i = 1, ao_num
- ! tmp = 0.d0
- ! do j = 1, ao_num
- ! tmp += ao_overlap(i,j) * coef_fit(j)
- ! enddo
- ! tmp = tmp - b(i)
- ! print*, i, tmp
- !enddo
+ integer :: j
+ double precision :: tmp, acc, nrm
- deallocate(A_inv, b)
+ acc = 0.d0
+ nrm = 0.d0
+ print *, ' check A_inv'
+ do i = 1, ao_num
+ tmp = 0.d0
+ do j = 1, ao_num
+ tmp += ao_overlap(i,j) * coef_fit(j)
+ enddo
+ tmp = tmp - b(i)
+ if(dabs(tmp) .gt. 1d-8) then
+ print*, ' problem found in fitting 1e-Jastrow'
+ print*, i, tmp
+ endif
+
+ acc += dabs(tmp)
+ nrm += dabs(b(i))
+ enddo
+ print *, ' Relative Error (%) =', 100.d0*acc/nrm
+
+ deallocate(A, A_inv, b)
call wall_time(t1)
print*, ' END after (min) ', (t1-t0)/60.d0
@@ -128,7 +137,7 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit)
PROVIDE mo_coef
call wall_time(t0)
- print*, ' PROVIDING the representation of 1e-Jastrow in AOs x AOx ... '
+ print*, ' PROVIDING the representation of 1e-Jastrow in AOs x AOs ... '
! --- --- ---
! get u1e(r)
@@ -188,10 +197,10 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit)
!$OMP END DO
!$OMP END PARALLEL
- print *, ' A'
- do ij = 1, ao_num*ao_num
- write(*, '(100000(f15.7))') (A(ij,kl), kl = 1, ao_num*ao_num)
- enddo
+! print *, ' A'
+! do ij = 1, ao_num*ao_num
+! write(*, '(100000(f15.7))') (A(ij,kl), kl = 1, ao_num*ao_num)
+! enddo
! --- --- ---
! get b
@@ -223,44 +232,35 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit)
! solve Ax = b
allocate(A_inv(ao_num*ao_num,ao_num*ao_num))
- call get_inverse(A, ao_num*ao_num, ao_num*ao_num, A_inv, ao_num*ao_num)
-
- integer :: mn
- print *, ' check A_inv'
- do ij = 1, ao_num*ao_num
- do kl = 1, ao_num*ao_num
-
- tmp = 0.d0
- do mn = 1, ao_num*ao_num
- tmp += A(ij,mn) * A_inv(mn,kl)
- enddo
-
- print*, ij, kl, tmp
- enddo
- enddo
+ !call get_inverse(A, ao_num*ao_num, ao_num*ao_num, A_inv, ao_num*ao_num)
+ call get_pseudo_inverse(A, ao_num*ao_num, ao_num*ao_num, ao_num*ao_num, A_inv, ao_num*ao_num, 5d-8)
! coef_fit = A_inv x b
- !call dgemv("N", ao_num*ao_num, ao_num*ao_num, 1.d0, A_inv, ao_num*ao_num, b, 1, 0.d0, coef_fit(1,1), 1)
- do ij = 1, ao_num*ao_num
- coef_fit(ij) = 0.d0
- do kl = 1, ao_num*ao_num
- coef_fit(ij) += A_inv(ij,kl) * b(kl)
- enddo
- enddo
+ call dgemv("N", ao_num*ao_num, ao_num*ao_num, 1.d0, A_inv, ao_num*ao_num, b, 1, 0.d0, coef_fit, 1)
- double precision :: tmp
- print *, ' check A_inv'
+ integer :: mn
+ double precision :: tmp, acc, nrm
+
+ acc = 0.d0
+ nrm = 0.d0
do ij = 1, ao_num*ao_num
tmp = 0.d0
do kl = 1, ao_num*ao_num
tmp += A(ij,kl) * coef_fit(kl)
enddo
tmp = tmp - b(ij)
- print*, ij, tmp
- enddo
+ if(dabs(tmp) .gt. 1d-7) then
+ print*, ' problem found in fitting 1e-Jastrow'
+ print*, ij, tmp
+ endif
- deallocate(A)
- deallocate(A_inv, b)
+ acc += dabs(tmp)
+ nrm += dabs(b(ij))
+ enddo
+ print *, ' Relative Error (%) =', 100.d0*acc/nrm
+
+
+ deallocate(A, A_inv, b)
call wall_time(t1)
print*, ' END after (min) ', (t1-t0)/60.d0
@@ -373,6 +373,7 @@ subroutine get_j1e_coef_fit_ao3(dim_fit, coef_fit)
enddo
tmp = tmp - b(i,d)
if(dabs(tmp) .gt. 1d-8) then
+ print*, ' problem found in fitting 1e-Jastrow'
print*, d, i, tmp
endif