PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-06-17 10:45:33 +02:00
Code Releases Activity

added 1e-term to Mu_Nu

Browse Source
This commit is contained in:
AbdAmmar 2024-02-04 19:56:23 +01:00
parent acd26fdeb0
commit da2ee20723
3 changed files with 85 additions and 23 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

34
plugins/local/non_h_ints_mu/jast_1e_utils.irp.f
View File

@ -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)
! ---

73
plugins/local/non_h_ints_mu/jast_deriv_utils_vect.irp.f
View File

@ -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'

1
plugins/local/non_h_ints_mu/tc_integ_num.irp.f
View File

@ -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