From da2ee2072305b36a9dc2b555566483d67c611a74 Mon Sep 17 00:00:00 2001 From: AbdAmmar Date: Sun, 4 Feb 2024 19:56:23 +0100 Subject: [PATCH] added 1e-term to Mu_Nu --- .../local/non_h_ints_mu/jast_1e_utils.irp.f | 34 +++------ .../non_h_ints_mu/jast_deriv_utils_vect.irp.f | 73 +++++++++++++++++++ .../local/non_h_ints_mu/tc_integ_num.irp.f | 1 - 3 files changed, 85 insertions(+), 23 deletions(-) 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 7aa85148..9cfabf58 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 @@ -126,9 +126,9 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit) integer :: ij, kl, mn integer :: info, n_svd, LWORK double precision :: g - double precision :: t0, t1 + double precision :: t0, t1, svd_t0, svd_t1 double precision :: cutoff_svd, D1_inv - double precision :: accu, norm, diff + double precision, allocatable :: diff(:) double precision, allocatable :: A(:,:,:,:), b(:), A_tmp(:,:,:,:) double precision, allocatable :: Pa(:,:), Pb(:,:), Pt(:,:) double precision, allocatable :: u1e_tmp(:), tmp(:,:,:) @@ -211,9 +211,6 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit) enddo call dgemv("T", n_points_final_grid, ao_num*ao_num, 1.d0, tmp(1,1,1), n_points_final_grid, u1e_tmp(1), 1, 0.d0, b(1), 1) - !call dgemm( "T", "N", ao_num*ao_num, 1, n_points_final_grid, 1.d0 & - ! , tmp(1,1,1), n_points_final_grid, u1e_tmp(1), n_points_final_grid & - ! , 0.d0, b(1), ao_num*ao_num) deallocate(u1e_tmp) deallocate(tmp) @@ -223,6 +220,8 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit) allocate(D(ao_num*ao_num), U(ao_num*ao_num,ao_num*ao_num), Vt(ao_num*ao_num,ao_num*ao_num)) + call wall_time(svd_t0) + allocate(work(1)) lwork = -1 call dgesvd( 'S', 'A', ao_num*ao_num, ao_num*ao_num, A(1,1,1,1), ao_num*ao_num & @@ -244,6 +243,9 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit) deallocate(work) + call wall_time(svd_t1) + print*, ' SVD time (min) ', (svd_t1-svd_t0)/60.d0 + if(D(1) .lt. 1d-14) then print*, ' largest singular value is very small:', D(1) n_svd = 1 @@ -289,24 +291,12 @@ subroutine get_j1e_coef_fit_ao2(dim_fit, coef_fit) ! --- - accu = 0.d0 - norm = 0.d0 - do k = 1, ao_num - do l = 1, ao_num - kl = (k-1)*ao_num + l - diff = 0.d0 - do i = 1, ao_num - do j = 1, ao_num - diff += A_tmp(k,l,i,j) * coef_fit(j,i) - enddo - enddo + allocate(diff(ao_num*ao_num)) - !print*, kl, b(kl) - accu += dabs(diff - b(kl)) - norm += dabs(b(kl)) - enddo - enddo - print*, ' accu total on Ax = b (%) = ', 100.d0*accu/norm + call dgemv("N", ao_num*ao_num, ao_num*ao_num, 1.d0, A_tmp(1,1,1,1), ao_num*ao_num, coef_fit(1,1), 1, 0.d0, diff(1), 1) + print*, ' accu total on Ax = b (%) = ', 100.d0*sum(dabs(diff-b))/sum(dabs(b)) + + deallocate(diff) deallocate(A_tmp) ! --- diff --git a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f index ffb7b513..5777a44a 100644 --- a/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f +++ b/plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f @@ -558,6 +558,8 @@ subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz) double precision, allocatable :: env_r2(:) double precision, allocatable :: u2b_r12(:) double precision, allocatable :: gradx1_u2b(:), grady1_u2b(:), gradz1_u2b(:) + double precision, allocatable :: u2b_mu(:), gradx1_mu(:), grady1_mu(:), gradz1_mu(:) + double precision, allocatable :: u2b_nu(:), gradx1_nu(:), grady1_nu(:), gradz1_nu(:) double precision, external :: env_nucl PROVIDE j1e_type j2e_type env_type @@ -604,6 +606,46 @@ subroutine get_grad1_u12_2e_r1_seq(ipoint, n_grid2, resx, resy, resz) endif ! env_type + elseif(j2e_type .eq. "Mu_Nu") then + + if(env_type .eq. "None") then + + call grad1_jmu_r1_seq(mu_erf, r1, n_grid2, resx, resy, resz) + + else + + ! u(r1,r2) = jmu(r12) x v(r1) x v(r2) + jnu(r12) x [1 - v(r1) x v(r2)] + + allocate(env_r2(n_grid2)) + allocate(u2b_mu(n_grid2)) + allocate(u2b_nu(n_grid2)) + allocate(gradx1_mu(n_grid2), grady1_mu(n_grid2), gradz1_mu(n_grid2)) + allocate(gradx1_nu(n_grid2), grady1_nu(n_grid2), gradz1_nu(n_grid2)) + + env_r1 = env_nucl(r1) + call grad1_env_nucl(r1, grad1_env) + call env_nucl_r1_seq(n_grid2, env_r2) + + call jmu_r1_seq(mu_erf, r1, n_grid2, u2b_mu) + call jmu_r1_seq(nu_erf, r1, n_grid2, u2b_nu) + + call grad1_jmu_r1_seq(mu_erf, r1, n_grid2, gradx1_mu, grady1_mu, gradz1_mu) + call grad1_jmu_r1_seq(nu_erf, r1, n_grid2, gradx1_nu, grady1_nu, gradz1_nu) + + do jpoint = 1, n_points_extra_final_grid + resx(jpoint) = gradx1_nu(jpoint) + ((gradx1_mu(jpoint) - gradx1_nu(jpoint)) * env_r1 + (u2b_mu(jpoint) - u2b_nu(jpoint)) * grad1_env(1)) * env_r2(jpoint) + resy(jpoint) = grady1_nu(jpoint) + ((grady1_mu(jpoint) - grady1_nu(jpoint)) * env_r1 + (u2b_mu(jpoint) - u2b_nu(jpoint)) * grad1_env(2)) * env_r2(jpoint) + resz(jpoint) = gradz1_nu(jpoint) + ((gradz1_mu(jpoint) - gradz1_nu(jpoint)) * env_r1 + (u2b_mu(jpoint) - u2b_nu(jpoint)) * grad1_env(3)) * env_r2(jpoint) + enddo + + deallocate(env_r2) + deallocate(u2b_mu) + deallocate(u2b_nu) + deallocate(gradx1_mu, grady1_mu, gradz1_mu) + deallocate(gradx1_nu, grady1_nu, gradz1_nu) + + endif ! env_type + else print *, ' Error in get_grad1_u12_withsq_r1_seq: Unknown Jastrow' @@ -635,6 +677,7 @@ subroutine get_u12_2e_r1_seq(ipoint, n_grid2, res) double precision :: grad1_env(3), r1(3) double precision, allocatable :: env_r2(:) double precision, allocatable :: u2b_r12(:) + double precision, allocatable :: u2b_mu(:), u2b_nu(:) double precision, external :: env_nucl PROVIDE j1e_type j2e_type env_type @@ -672,6 +715,36 @@ subroutine get_u12_2e_r1_seq(ipoint, n_grid2, res) endif ! env_type + elseif(j2e_type .eq. "Mu_Nu") then + + if(env_type .eq. "None") then + + call jmu_r1_seq(mu_erf, r1, n_grid2, res) + + else + + ! u(r1,r2) = jmu(r12) x v(r1) x v(r2) + jnu(r12) x [1 - v(r1) x v(r2)] + + allocate(env_r2(n_grid2)) + allocate(u2b_mu(n_grid2)) + allocate(u2b_nu(n_grid2)) + + env_r1 = env_nucl(r1) + call env_nucl_r1_seq(n_grid2, env_r2) + + call jmu_r1_seq(mu_erf, r1, n_grid2, u2b_mu) + call jmu_r1_seq(nu_erf, r1, n_grid2, u2b_nu) + + do jpoint = 1, n_points_extra_final_grid + res(jpoint) = u2b_nu(jpoint) + (u2b_mu(jpoint) - u2b_nu(jpoint)) * env_r1 * env_r2(jpoint) + enddo + + deallocate(env_r2) + deallocate(u2b_mu) + deallocate(u2b_nu) + + endif ! env_type + else print *, ' Error in get_u12_withsq_r1_seq: Unknown Jastrow' diff --git a/plugins/local/non_h_ints_mu/tc_integ_num.irp.f b/plugins/local/non_h_ints_mu/tc_integ_num.irp.f index 6b6e755d..e5d75c3d 100644 --- a/plugins/local/non_h_ints_mu/tc_integ_num.irp.f +++ b/plugins/local/non_h_ints_mu/tc_integ_num.irp.f @@ -45,7 +45,6 @@ !$OMP END DO !$OMP END PARALLEL - ! n_points_final_grid = n_blocks * n_pass + n_rest call total_memory(mem) mem = max(1.d0, qp_max_mem - mem) n_double = mem * 1.d8