mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-06-20 04:02:22 +02:00
553 lines
14 KiB
Fortran
553 lines
14 KiB
Fortran
|
|
! ---
|
|
|
|
double precision function num_v_ij_u_cst_mu_env(i, j, ipoint)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! \int dr2 u12 \phi_i(r2) \phi_j(r2) x v_env(r2)
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
|
|
integer :: jpoint
|
|
double precision :: r1(3), r2(3)
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: j12_mu, env_nucl, j12_mu_gauss
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
num_v_ij_u_cst_mu_env = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(3,jpoint)
|
|
|
|
num_v_ij_u_cst_mu_env += ao_value(i, r2) * ao_value(j, r2) * j12_mu_gauss(r1, r2) * env_nucl(r2) * final_weight_at_r_vector(jpoint)
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
double precision function num_int2_u2_env2(i, j, ipoint)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! \int dr2 u12^2 \phi_i(r2) \phi_j(r2) x v_env(r2)^2
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
|
|
integer :: jpoint, i_fit
|
|
double precision :: r1(3), r2(3)
|
|
double precision :: dx, dy, dz, r12, x2, tmp1, tmp2, tmp3, coef, expo
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
double precision, external :: j12_mu
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
num_int2_u2_env2 = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(3,jpoint)
|
|
dx = r1(1) - r2(1)
|
|
dy = r1(2) - r2(2)
|
|
dz = r1(3) - r2(3)
|
|
x2 = dx * dx + dy * dy + dz * dz
|
|
r12 = dsqrt(x2)
|
|
|
|
tmp1 = env_nucl(r2)
|
|
tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
|
|
|
|
!tmp3 = 0.d0
|
|
!do i_fit = 1, n_max_fit_slat
|
|
! expo = expo_gauss_j_mu_x_2(i_fit)
|
|
! coef = coef_gauss_j_mu_x_2(i_fit)
|
|
! tmp3 += coef * dexp(-expo*x2)
|
|
!enddo
|
|
tmp3 = j12_mu(r1, r2)
|
|
tmp3 = tmp3 * tmp3
|
|
|
|
num_int2_u2_env2 += tmp2 * tmp3
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
double precision function num_int2_grad1u2_grad2u2_env2(i, j, ipoint)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! \int dr2 \frac{-[erf(mu r12) -1]^2}{4} \phi_i(r2) \phi_j(r2) x v_env(r2)^2
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
|
|
integer :: jpoint, i_fit
|
|
double precision :: r1(3), r2(3)
|
|
double precision :: dx, dy, dz, r12, x2, tmp1, tmp2, tmp3, coef, expo
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
num_int2_grad1u2_grad2u2_env2 = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(3,jpoint)
|
|
dx = r1(1) - r2(1)
|
|
dy = r1(2) - r2(2)
|
|
dz = r1(3) - r2(3)
|
|
x2 = dx * dx + dy * dy + dz * dz
|
|
r12 = dsqrt(x2)
|
|
|
|
tmp1 = env_nucl(r2)
|
|
tmp2 = tmp1 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
|
|
|
|
!tmp3 = 0.d0
|
|
!do i_fit = 1, n_max_fit_slat
|
|
! expo = expo_gauss_1_erf_x_2(i_fit)
|
|
! coef = coef_gauss_1_erf_x_2(i_fit)
|
|
! tmp3 += coef * dexp(-expo*x2)
|
|
!enddo
|
|
tmp3 = derf(mu_erf*r12) - 1.d0
|
|
tmp3 = tmp3 * tmp3
|
|
|
|
tmp3 = -0.25d0 * tmp3
|
|
|
|
num_int2_grad1u2_grad2u2_env2 += tmp2 * tmp3
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
double precision function num_v_ij_erf_rk_cst_mu_env(i, j, ipoint)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! \int dr2 [erf(mu r12) -1]/r12 \phi_i(r2) \phi_j(r2) x v_env(r2)
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
|
|
integer :: jpoint
|
|
double precision :: r1(3), r2(3)
|
|
double precision :: dx, dy, dz, r12, tmp1, tmp2
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
num_v_ij_erf_rk_cst_mu_env = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(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) cycle
|
|
|
|
tmp1 = (derf(mu_erf * r12) - 1.d0) / r12
|
|
tmp2 = tmp1 * ao_value(i, r2) * ao_value(j, r2) * env_nucl(r2) * final_weight_at_r_vector(jpoint)
|
|
|
|
num_v_ij_erf_rk_cst_mu_env += tmp2
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
subroutine num_x_v_ij_erf_rk_cst_mu_env(i, j, ipoint, integ)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! \int dr2 [erf(mu r12) -1]/r12 \phi_i(r2) \phi_j(r2) x v_env(r2) x r2
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
double precision, intent(out) :: integ(3)
|
|
|
|
integer :: jpoint
|
|
double precision :: r1(3), r2(3), grad(3)
|
|
double precision :: dx, dy, dz, r12, tmp1, tmp2
|
|
double precision :: tmp_x, tmp_y, tmp_z
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
tmp_x = 0.d0
|
|
tmp_y = 0.d0
|
|
tmp_z = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(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) cycle
|
|
|
|
tmp1 = (derf(mu_erf * r12) - 1.d0) / r12
|
|
tmp2 = tmp1 * ao_value(i, r2) * ao_value(j, r2) * env_nucl(r2) * final_weight_at_r_vector(jpoint)
|
|
|
|
tmp_x += tmp2 * r2(1)
|
|
tmp_y += tmp2 * r2(2)
|
|
tmp_z += tmp2 * r2(3)
|
|
enddo
|
|
|
|
integ(1) = tmp_x
|
|
integ(2) = tmp_y
|
|
integ(3) = tmp_z
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
subroutine num_int2_grad1_u12_ao(i, j, ipoint, integ)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! \int dr2 [-grad_1 u12] \phi_i(r2) \phi_j(r2) x v12_env(r1, r2)
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
double precision, intent(out) :: integ(3)
|
|
|
|
integer :: jpoint
|
|
double precision :: tmp, r1(3), r2(3), grad(3)
|
|
double precision :: tmp_x, tmp_y, tmp_z
|
|
|
|
double precision, external :: ao_value
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
tmp_x = 0.d0
|
|
tmp_y = 0.d0
|
|
tmp_z = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(3,jpoint)
|
|
tmp = ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
|
|
|
|
call grad1_jmu_modif_num(r1, r2, grad)
|
|
|
|
tmp_x += tmp * (-1.d0 * grad(1))
|
|
tmp_y += tmp * (-1.d0 * grad(2))
|
|
tmp_z += tmp * (-1.d0 * grad(3))
|
|
enddo
|
|
|
|
integ(1) = tmp_x
|
|
integ(2) = tmp_y
|
|
integ(3) = tmp_z
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
double precision function num_grad12_j12(i, j, ipoint)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [v1^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) ]
|
|
!
|
|
END_DOC
|
|
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
|
|
integer :: jpoint
|
|
double precision :: r1(3), r2(3)
|
|
double precision :: tmp_x, tmp_y, tmp_z, r12
|
|
double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3)
|
|
double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp
|
|
double precision :: fst_term, scd_term, thd_term, tmp
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
double precision, external :: j12_mu
|
|
double precision, external :: grad_x_env_nucl_num
|
|
double precision, external :: grad_y_env_nucl_num
|
|
double precision, external :: grad_z_env_nucl_num
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
num_grad12_j12 = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(3,jpoint)
|
|
|
|
tmp_x = r1(1) - r2(1)
|
|
tmp_y = r1(2) - r2(2)
|
|
tmp_z = r1(3) - r2(3)
|
|
r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
|
|
|
|
dx1_v1 = grad_x_env_nucl_num(r1)
|
|
dy1_v1 = grad_y_env_nucl_num(r1)
|
|
dz1_v1 = grad_z_env_nucl_num(r1)
|
|
|
|
call grad1_j12_mu(r1, r2, grad_u12)
|
|
|
|
tmp1 = 1.d0 - derf(mu_erf * r12)
|
|
v1_tmp = env_nucl(r1)
|
|
v2_tmp = env_nucl(r2)
|
|
u12_tmp = j12_mu(r1, r2)
|
|
|
|
fst_term = 0.5d0 * tmp1 * tmp1 * v1_tmp * v1_tmp
|
|
|
|
tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * fst_term * v2_tmp * v2_tmp
|
|
|
|
num_grad12_j12 += tmp
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
double precision function num_u12sq_envsq(i, j, ipoint)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [ u12^2 (grad_1 v1)^2 ]
|
|
!
|
|
END_DOC
|
|
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
|
|
integer :: jpoint
|
|
double precision :: r1(3), r2(3)
|
|
double precision :: tmp_x, tmp_y, tmp_z, r12
|
|
double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3)
|
|
double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp
|
|
double precision :: fst_term, scd_term, thd_term, tmp
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
double precision, external :: j12_mu
|
|
double precision, external :: grad_x_env_nucl_num
|
|
double precision, external :: grad_y_env_nucl_num
|
|
double precision, external :: grad_z_env_nucl_num
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
num_u12sq_envsq = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(3,jpoint)
|
|
|
|
tmp_x = r1(1) - r2(1)
|
|
tmp_y = r1(2) - r2(2)
|
|
tmp_z = r1(3) - r2(3)
|
|
r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
|
|
|
|
dx1_v1 = grad_x_env_nucl_num(r1)
|
|
dy1_v1 = grad_y_env_nucl_num(r1)
|
|
dz1_v1 = grad_z_env_nucl_num(r1)
|
|
|
|
call grad1_j12_mu(r1, r2, grad_u12)
|
|
|
|
tmp1 = 1.d0 - derf(mu_erf * r12)
|
|
v1_tmp = env_nucl(r1)
|
|
v2_tmp = env_nucl(r2)
|
|
u12_tmp = j12_mu(r1, r2)
|
|
|
|
scd_term = u12_tmp * u12_tmp * (dx1_v1*dx1_v1 + dy1_v1*dy1_v1 + dz1_v1*dz1_v1)
|
|
|
|
tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * scd_term * v2_tmp * v2_tmp
|
|
|
|
num_u12sq_envsq += tmp
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
double precision function num_u12_grad1_u12_env_grad1_env(i, j, ipoint)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! -0.50 x \int r2 \phi_i(2) \phi_j(2) x v2^2 [ 2 u12 v1 (grad_1 u12) . (grad_1 v1) ]
|
|
!
|
|
END_DOC
|
|
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
|
|
integer :: jpoint
|
|
double precision :: r1(3), r2(3)
|
|
double precision :: tmp_x, tmp_y, tmp_z, r12
|
|
double precision :: dx1_v1, dy1_v1, dz1_v1, grad_u12(3)
|
|
double precision :: tmp1, v1_tmp, v2_tmp, u12_tmp
|
|
double precision :: fst_term, scd_term, thd_term, tmp
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
double precision, external :: j12_mu
|
|
double precision, external :: grad_x_env_nucl_num
|
|
double precision, external :: grad_y_env_nucl_num
|
|
double precision, external :: grad_z_env_nucl_num
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
num_u12_grad1_u12_env_grad1_env = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(3,jpoint)
|
|
|
|
tmp_x = r1(1) - r2(1)
|
|
tmp_y = r1(2) - r2(2)
|
|
tmp_z = r1(3) - r2(3)
|
|
r12 = dsqrt(tmp_x*tmp_x + tmp_y*tmp_y + tmp_z*tmp_z)
|
|
|
|
dx1_v1 = grad_x_env_nucl_num(r1)
|
|
dy1_v1 = grad_y_env_nucl_num(r1)
|
|
dz1_v1 = grad_z_env_nucl_num(r1)
|
|
|
|
call grad1_j12_mu(r1, r2, grad_u12)
|
|
|
|
tmp1 = 1.d0 - derf(mu_erf * r12)
|
|
v1_tmp = env_nucl(r1)
|
|
v2_tmp = env_nucl(r2)
|
|
u12_tmp = j12_mu(r1, r2)
|
|
|
|
thd_term = 2.d0 * v1_tmp * u12_tmp * (dx1_v1*grad_u12(1) + dy1_v1*grad_u12(2) + dz1_v1*grad_u12(3))
|
|
|
|
tmp = -0.5d0 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint) * thd_term * v2_tmp * v2_tmp
|
|
|
|
num_u12_grad1_u12_env_grad1_env += tmp
|
|
enddo
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|
|
|
|
subroutine num_int2_u_grad1u_total_env2(i, j, ipoint, integ)
|
|
|
|
BEGIN_DOC
|
|
!
|
|
! \int dr2 u12 (grad_1 u12) \phi_i(r2) \phi_j(r2) x v_env(r2)^2
|
|
!
|
|
END_DOC
|
|
|
|
implicit none
|
|
|
|
integer, intent(in) :: i, j, ipoint
|
|
double precision, intent(out) :: integ(3)
|
|
|
|
integer :: jpoint
|
|
double precision :: r1(3), r2(3), grad(3)
|
|
double precision :: dx, dy, dz, r12, tmp0, tmp1, tmp2
|
|
double precision :: tmp_x, tmp_y, tmp_z
|
|
|
|
double precision, external :: ao_value
|
|
double precision, external :: env_nucl
|
|
double precision, external :: j12_mu
|
|
|
|
r1(1) = final_grid_points(1,ipoint)
|
|
r1(2) = final_grid_points(2,ipoint)
|
|
r1(3) = final_grid_points(3,ipoint)
|
|
|
|
tmp_x = 0.d0
|
|
tmp_y = 0.d0
|
|
tmp_z = 0.d0
|
|
do jpoint = 1, n_points_final_grid
|
|
r2(1) = final_grid_points(1,jpoint)
|
|
r2(2) = final_grid_points(2,jpoint)
|
|
r2(3) = final_grid_points(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) cycle
|
|
|
|
tmp0 = env_nucl(r2)
|
|
tmp1 = 0.5d0 * j12_mu(r1, r2) * (1.d0 - derf(mu_erf * r12)) / r12
|
|
tmp2 = tmp0 * tmp0 * tmp1 * ao_value(i, r2) * ao_value(j, r2) * final_weight_at_r_vector(jpoint)
|
|
|
|
tmp_x += tmp2 * dx
|
|
tmp_y += tmp2 * dy
|
|
tmp_z += tmp2 * dz
|
|
enddo
|
|
|
|
integ(1) = tmp_x
|
|
integ(2) = tmp_y
|
|
integ(3) = tmp_z
|
|
|
|
return
|
|
end
|
|
|
|
! ---
|