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

Merge pull request #30 from AbdAmmar/dev-stable-tc-scf

Browse Source 
Dev stable tc scf
This commit is contained in:
AbdAmmar 2024-01-16 19:18:19 +01:00 committed by GitHub
commit a9f4f3324f
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
305 changed files with 5792 additions and 5840 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

1
bin/qp_convert_output_to_ezfio
View File

@ -256,6 +256,7 @@ def write_ezfio(res, filename):
MoTag = res.determinants_mo_type
ezfio.set_mo_basis_mo_label('Orthonormalized')
ezfio.set_determinants_mo_label('Orthonormalized')
MO_type = MoTag
allMOs = res.mo_sets[MO_type]

3
bin/zcat
View File

@ -16,7 +16,8 @@ with gzip.open("$1", "rt") as f:
EOF
fi
else
command=$(which -a zcat | grep -v 'qp2/bin/' | head -1)
SCRIPTPATH="$( cd -- "$(dirname "$0")" >/dev/null 2>&1 ; pwd -P )"
command=$(which -a zcat | grep -v "$SCRIPTPATH/" | head -1)
exec $command $@
fi

12
configure vendored
View File

@ -195,7 +195,7 @@ if [[ "${PACKAGES}.x" != ".x" ]] ; then
fi
if [[ ${PACKAGES} = all ]] ; then
PACKAGES="zlib ninja zeromq f77zmq gmp ocaml docopt resultsFile bats trexio qmckl"
PACKAGES="zlib ninja zeromq f77zmq gmp ocaml docopt resultsFile bats trexio"
fi
@ -402,11 +402,11 @@ if [[ ${TREXIO} = $(not_found) ]] ; then
fail
fi
QMCKL=$(find_lib -lqmckl)
if [[ ${QMCKL} = $(not_found) ]] ; then
error "QMCkl (qmckl | qmckl-intel) is not installed."
fail
fi
#QMCKL=$(find_lib -lqmckl)
#if [[ ${QMCKL} = $(not_found) ]] ; then
# error "QMCkl (qmckl | qmckl-intel) is not installed."
# fail
#fi
F77ZMQ=$(find_lib -lzmq -lf77zmq -lpthread)
if [[ ${F77ZMQ} = $(not_found) ]] ; then

25
ocaml/Input_determinants_by_hand.ml
View File

@ -13,6 +13,7 @@ module Determinants_by_hand : sig
psi_coef : Det_coef.t array;
psi_det : Determinant.t array;
state_average_weight : Positive_float.t array;
mo_label : MO_label.t;
} [@@deriving sexp]
val read : ?full:bool -> unit -> t option
val write : ?force:bool -> t -> unit
@ -34,11 +35,21 @@ end = struct
psi_coef : Det_coef.t array;
psi_det : Determinant.t array;
state_average_weight : Positive_float.t array;
mo_label : MO_label.t;
} [@@deriving sexp]
;;
let get_default = Qpackage.get_ezfio_default "determinants";;
let read_mo_label () =
if not (Ezfio.has_determinants_mo_label ()) then
if Ezfio.has_mo_basis_mo_label () then (
let label = Ezfio.get_mo_basis_mo_label () in
Ezfio.set_determinants_mo_label label) ;
Ezfio.get_determinants_mo_label ()
|> MO_label.of_string
;;
let read_n_int () =
if not (Ezfio.has_determinants_n_int()) then
Ezfio.get_mo_basis_mo_num ()
@ -290,6 +301,10 @@ end = struct
Ezfio.set_determinants_psi_det_qp_edit r
;;
let write_mo_label a =
MO_label.to_string a
|> Ezfio.set_determinants_mo_label
let read ?(full=true) () =
@ -311,6 +326,7 @@ end = struct
psi_det = read_psi_det ~read_only () ;
n_states = read_n_states () ;
state_average_weight = read_state_average_weight () ;
mo_label = read_mo_label () ;
}
with _ -> None
else
@ -328,6 +344,7 @@ end = struct
psi_det ;
n_states ;
state_average_weight ;
mo_label ;
} =
write_n_int n_int ;
write_bit_kind bit_kind;
@ -340,7 +357,9 @@ end = struct
write_psi_coef ~n_det:n_det ~n_states:n_states psi_coef ;
write_psi_det ~n_int:n_int ~n_det:n_det psi_det
end;
write_state_average_weight state_average_weight
write_state_average_weight state_average_weight ;
write_mo_label mo_label ;
()
;;
@ -545,6 +564,8 @@ psi_det = %s
let bitkind =
Printf.sprintf "(bit_kind %d)" (Lazy.force Qpackage.bit_kind
|> Bit_kind.to_int)
and mo_label =
Printf.sprintf "(mo_label %s)" (MO_label.to_string @@ read_mo_label ())
and n_int =
Printf.sprintf "(n_int %d)" (N_int_number.get_max ())
and n_states =
@ -553,7 +574,7 @@ psi_det = %s
Printf.sprintf "(n_det_qp_edit %d)" (Det_number.to_int @@ read_n_det_qp_edit ())
in
let s =
String.concat "" [ header ; bitkind ; n_int ; n_states ; psi_coef ; psi_det ; n_det_qp_edit ]
String.concat "" [ header ; mo_label ; bitkind ; n_int ; n_states ; psi_coef ; psi_det ; n_det_qp_edit ]
in

4
ocaml/qptypes_generator.ml
View File

@ -154,8 +154,8 @@ let input_ezfio = "
* N_int_number : int
determinants_n_int
1 : 30
N_int > 30
1 : 128
N_int > 128
* Det_number : int
determinants_n_det

1
plugins/.gitignore vendored
View File

@ -1,2 +1 @@
*

1
src/ao_many_one_e_ints/NEED → plugins/local/ao_many_one_e_ints/NEED
View File

@ -4,3 +4,4 @@ becke_numerical_grid
mo_one_e_ints
dft_utils_in_r
tc_keywords
hamiltonian

0
src/ao_many_one_e_ints/README.rst → plugins/local/ao_many_one_e_ints/README.rst
View File

38
src/ao_many_one_e_ints/ao_erf_gauss.irp.f → plugins/local/ao_many_one_e_ints/ao_erf_gauss.irp.f
View File

@ -98,7 +98,7 @@ double precision function phi_j_erf_mu_r_phi(i, j, mu_in, C_center)
enddo
enddo
end function phi_j_erf_mu_r_phi
end
! ---
@ -201,7 +201,7 @@ subroutine erf_mu_gauss_ij_ao(i, j, mu, C_center, delta, gauss_ints)
enddo
enddo
end subroutine erf_mu_gauss_ij_ao
end
! ---
@ -266,7 +266,7 @@ subroutine NAI_pol_x_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
enddo
enddo
end subroutine NAI_pol_x_mult_erf_ao
end
! ---
@ -340,7 +340,7 @@ subroutine NAI_pol_x_mult_erf_ao_v0(i_ao, j_ao, mu_in, C_center, LD_C, ints, LD_
deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_v0
end
! ---
@ -420,7 +420,7 @@ subroutine NAI_pol_x_mult_erf_ao_v(i_ao, j_ao, mu_in, C_center, LD_C, ints, LD_i
deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_v
end
! ---
@ -479,7 +479,7 @@ double precision function NAI_pol_x_mult_erf_ao_x(i_ao, j_ao, mu_in, C_center)
enddo
enddo
end function NAI_pol_x_mult_erf_ao_x
end
! ---
@ -538,7 +538,7 @@ double precision function NAI_pol_x_mult_erf_ao_y(i_ao, j_ao, mu_in, C_center)
enddo
enddo
end function NAI_pol_x_mult_erf_ao_y
end
! ---
@ -597,7 +597,7 @@ double precision function NAI_pol_x_mult_erf_ao_z(i_ao, j_ao, mu_in, C_center)
enddo
enddo
end function NAI_pol_x_mult_erf_ao_z
end
! ---
@ -667,7 +667,7 @@ double precision function NAI_pol_x_mult_erf_ao_with1s_x(i_ao, j_ao, beta, B_cen
enddo
enddo
end function NAI_pol_x_mult_erf_ao_with1s_x
end
! ---
@ -737,7 +737,7 @@ double precision function NAI_pol_x_mult_erf_ao_with1s_y(i_ao, j_ao, beta, B_cen
enddo
enddo
end function NAI_pol_x_mult_erf_ao_with1s_y
end
! ---
@ -807,7 +807,7 @@ double precision function NAI_pol_x_mult_erf_ao_with1s_z(i_ao, j_ao, beta, B_cen
enddo
enddo
end function NAI_pol_x_mult_erf_ao_with1s_z
end
! ---
@ -880,7 +880,7 @@ subroutine NAI_pol_x_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_cen
enddo
enddo
end subroutine NAI_pol_x_mult_erf_ao_with1s
end
! ---
@ -967,7 +967,7 @@ subroutine NAI_pol_x_mult_erf_ao_with1s_v0(i_ao, j_ao, beta, B_center, LD_B, mu_
deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_with1s_v0
end
! ---
@ -1057,7 +1057,7 @@ subroutine NAI_pol_x_mult_erf_ao_with1s_v(i_ao, j_ao, beta, B_center, LD_B, mu_i
deallocate(integral)
end subroutine NAI_pol_x_mult_erf_ao_with1s_v
end
! ---
@ -1175,7 +1175,7 @@ subroutine NAI_pol_x2_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_ce
enddo
enddo
end subroutine NAI_pol_x2_mult_erf_ao_with1s
end
! ---
@ -1241,7 +1241,7 @@ subroutine NAI_pol_x2_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
enddo
enddo
end subroutine NAI_pol_x2_mult_erf_ao
end
! ---
@ -1320,7 +1320,7 @@ subroutine NAI_pol_012_mult_erf_ao_with1s(i_ao, j_ao, beta, B_center, mu_in, C_c
enddo
enddo
end subroutine NAI_pol_012_mult_erf_ao_with1s
end
! ---
@ -1328,7 +1328,7 @@ subroutine NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
BEGIN_DOC
!
! Computes the following integral :
! Computes the following integrals :
!
! int(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) \frac{\erf(\mu | r - R_C | )}{ | r - R_C | }$.
!
@ -1395,7 +1395,7 @@ subroutine NAI_pol_012_mult_erf_ao(i_ao, j_ao, mu_in, C_center, ints)
enddo
enddo
end subroutine NAI_pol_012_mult_erf_ao
end
! ---

0
src/ao_many_one_e_ints/ao_erf_gauss_grad.irp.f → plugins/local/ao_many_one_e_ints/ao_erf_gauss_grad.irp.f
View File

89
src/ao_many_one_e_ints/ao_gaus_gauss.irp.f → plugins/local/ao_many_one_e_ints/ao_gaus_gauss.irp.f
View File

@ -152,7 +152,7 @@ double precision function overlap_gauss_r12_ao(D_center, delta, i, j)
enddo
enddo
end function overlap_gauss_r12_ao
end
! --
@ -199,7 +199,7 @@ double precision function overlap_abs_gauss_r12_ao(D_center, delta, i, j)
enddo
enddo
end function overlap_gauss_r12_ao
end
! --
@ -257,7 +257,7 @@ subroutine overlap_gauss_r12_ao_v(D_center, LD_D, delta, i, j, resv, LD_resv, n_
deallocate(analytical_j)
end subroutine overlap_gauss_r12_ao_v
end
! ---
@ -327,7 +327,7 @@ double precision function overlap_gauss_r12_ao_with1s(B_center, beta, D_center,
enddo
enddo
end function overlap_gauss_r12_ao_with1s
end
! ---
@ -420,7 +420,86 @@ subroutine overlap_gauss_r12_ao_with1s_v(B_center, beta, D_center, LD_D, delta,
deallocate(fact_g, G_center, analytical_j)
end subroutine overlap_gauss_r12_ao_with1s_v
end
! ---
subroutine overlap_gauss_r12_ao_012(D_center, delta, i, j, ints)
BEGIN_DOC
!
! Computes the following integrals :
!
! ints(1) = $\int_{-\infty}^{infty} dr x^0 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
! ints(2) = $\int_{-\infty}^{infty} dr x^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(3) = $\int_{-\infty}^{infty} dr y^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(4) = $\int_{-\infty}^{infty} dr z^1 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
! ints(5) = $\int_{-\infty}^{infty} dr x^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(6) = $\int_{-\infty}^{infty} dr y^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
! ints(7) = $\int_{-\infty}^{infty} dr z^2 * \chi_i(r) \chi_j(r) e^{-\delta (r - D_center)^2}
!
END_DOC
include 'utils/constants.include.F'
implicit none
integer, intent(in) :: i, j
double precision, intent(in) :: delta, D_center(3)
double precision, intent(out) :: ints(7)
integer :: k, l, m
integer :: power_A(3), power_B(3), power_A1(3), power_A2(3)
double precision :: A_center(3), B_center(3), alpha, beta, coef1, coef
double precision :: integral0, integral1, integral2
double precision, external :: overlap_gauss_r12
ints = 0.d0
if(ao_overlap_abs(j,i).lt.1.d-12) then
return
endif
power_A(1:3) = ao_power(i,1:3)
power_B(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(ao_nucl(i),1:3)
B_center(1:3) = nucl_coord(ao_nucl(j),1:3)
do l = 1, ao_prim_num(i)
alpha = ao_expo_ordered_transp (l,i)
coef1 = ao_coef_normalized_ordered_transp(l,i)
do k = 1, ao_prim_num(j)
beta = ao_expo_ordered_transp(k,j)
coef = coef1 * ao_coef_normalized_ordered_transp(k,j)
if(dabs(coef) .lt. 1d-12) cycle
integral0 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A, power_B, alpha, beta)
ints(1) += coef * integral0
do m = 1, 3
power_A1 = power_A
power_A1(m) += 1
integral1 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A1, power_B, alpha, beta)
ints(1+m) += coef * (integral1 + A_center(m)*integral0)
power_A2 = power_A
power_A2(m) += 2
integral2 = overlap_gauss_r12(D_center, delta, A_center, B_center, power_A2, power_B, alpha, beta)
ints(4+m) += coef * (integral2 + A_center(m) * (2.d0*integral1 + A_center(m)*integral0))
enddo
enddo ! k
enddo ! l
return
end
! ---

198
src/ao_many_one_e_ints/grad2_jmu_manu.irp.f → plugins/local/ao_many_one_e_ints/grad2_jmu_manu.irp.f
View File

@ -1,11 +1,11 @@
! ---
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, int2_grad1u2_grad2u2_env2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [1 - erf(mu r12)]^2
!
END_DOC
@ -15,30 +15,30 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
double precision :: coef, beta, B_center(3)
double precision :: tmp
double precision :: wall0, wall1
double precision :: int_gauss, dsqpi_3_2, int_j1b
double precision :: int_gauss, dsqpi_3_2, int_env
double precision :: factor_ij_1s, beta_ij, center_ij_1s(3), sq_pi_3_2
double precision, allocatable :: int_fit_v(:)
double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_grad1u2_grad2u2_j1b2_test ...'
print*, ' providing int2_grad1u2_grad2u2_env2_test ...'
sq_pi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points_transp j1b_pen List_comb_thr_b3_coef
provide mu_erf final_grid_points_transp List_comb_thr_b3_coef
call wall_time(wall0)
int2_grad1u2_grad2u2_j1b2_test(:,:,:) = 0.d0
int2_grad1u2_grad2u2_env2_test(:,:,:) = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit_v, tmp,int_gauss,int_j1b,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points,List_comb_thr_b3_size, &
!$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, &
!$OMP ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc)
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit_v, tmp,int_gauss,int_env,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points,List_comb_thr_b3_size, &
!$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_env2_test, ao_abs_comb_b3_env, &
!$OMP ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc)
!$OMP DO SCHEDULE(dynamic)
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
@ -54,13 +54,13 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
! i_1s = 1
! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i)
int_env = ao_abs_comb_b3_env(1,j,i)
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit)
! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
! if(dabs(coef_fit*int_env*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j)
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss
int2_grad1u2_grad2u2_env2_test(j,i,ipoint) += coef_fit * int_gauss
enddo
! --- --- ---
@ -71,7 +71,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
int_env = ao_abs_comb_b3_env(i_1s,j,i)
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -81,11 +81,11 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
!DIR$ FORCEINLINE
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef
! if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
! if(dabs(coef_fit*factor_ij_1s*int_env*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, &
! expo_fit, i, j, int_fit_v, n_points_final_grid)
int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss
int2_grad1u2_grad2u2_env2_test(j,i,ipoint) += coef_fit * int_gauss
enddo
enddo
@ -98,26 +98,26 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
do j = 1, i-1
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2_test(i,j,ipoint)
int2_grad1u2_grad2u2_env2_test(j,i,ipoint) = int2_grad1u2_grad2u2_env2_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test', wall1 - wall0
print*, ' wall time for int2_grad1u2_grad2u2_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao_num, n_points_final_grid)]
!
! BEGIN_DOC
! !
! ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2
! !
! END_DOC
!
BEGIN_PROVIDER [double precision, int2_grad1u2_grad2u2_env2_test_v, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [1 - erf(mu r12)]^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), expo_fit, coef_fit
@ -128,24 +128,24 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
double precision, allocatable :: int_fit_v(:),big_array(:,:,:)
double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_grad1u2_grad2u2_j1b2_test_v ...'
print*, ' providing int2_grad1u2_grad2u2_env2_test_v ...'
provide mu_erf final_grid_points_transp j1b_pen
provide mu_erf final_grid_points_transp
call wall_time(wall0)
double precision :: int_j1b
double precision :: int_env
big_array(:,:,:) = 0.d0
allocate(big_array(n_points_final_grid,ao_num, ao_num))
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,&
!$OMP coef_fit, expo_fit, int_fit_v, tmp,int_j1b) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size,&
!$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, big_array,&
!$OMP ao_abs_comb_b3_j1b,ao_overlap_abs,thrsh_cycle_tc)
!
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center,&
!$OMP coef_fit, expo_fit, int_fit_v, tmp,int_env) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size,&
!$OMP final_grid_points_transp, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, big_array,&
!$OMP ao_abs_comb_b3_env,ao_overlap_abs,thrsh_cycle_tc)
!
allocate(int_fit_v(n_points_final_grid))
!$OMP DO SCHEDULE(dynamic)
do i = 1, ao_num
@ -159,7 +159,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
int_env = ao_abs_comb_b3_env(i_1s,j,i)
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -187,7 +187,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
do i = 1, ao_num
do j = i, ao_num
do ipoint = 1, n_points_final_grid
int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,j,i)
int2_grad1u2_grad2u2_env2_test_v(j,i,ipoint) = big_array(ipoint,j,i)
enddo
enddo
enddo
@ -195,23 +195,23 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_grad1u2_grad2u2_j1b2_test_v(j,i,ipoint) = big_array(ipoint,i,j)
int2_grad1u2_grad2u2_env2_test_v(j,i,ipoint) = big_array(ipoint,i,j)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2_test_v', wall1 - wall0
print*, ' wall time for int2_grad1u2_grad2u2_env2_test_v (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, int2_u2_env2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [u_12^mu]^2
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [u_12^mu]^2
!
END_DOC
@ -219,29 +219,29 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), int_fit, expo_fit, coef_fit
double precision :: coef, beta, B_center(3), tmp
double precision :: wall0, wall1,int_j1b
double precision :: wall0, wall1,int_env
double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s
double precision :: factor_ij_1s,beta_ij,center_ij_1s(3),sq_pi_3_2
print*, ' providing int2_u2_j1b2_test ...'
print*, ' providing int2_u2_env2_test ...'
sq_pi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf final_grid_points
call wall_time(wall0)
int2_u2_j1b2_test = 0.d0
int2_u2_env2_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp, int_j1b,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP coef_fit, expo_fit, int_fit, tmp, int_env,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2, &
!$OMP List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b,thrsh_cycle_tc)
!$OMP List_comb_thr_b3_cent, int2_u2_env2_test,ao_abs_comb_b3_env,thrsh_cycle_tc)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
@ -257,12 +257,12 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
! i_1s = 1
! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i)
if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
int_env = ao_abs_comb_b3_env(1,j,i)
if(dabs(int_env).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_x_2(i_fit)
coef_fit = coef_gauss_j_mu_x_2(i_fit)
! if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
! if(dabs(coef_fit*int_env*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += coef_fit * int_fit
enddo
@ -275,8 +275,8 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
int_env = ao_abs_comb_b3_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -286,13 +286,13 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
coef_fit = coef_gauss_j_mu_x_2(i_fit)
!DIR$ FORCEINLINE
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
! if(dabs(coef_fit*coef*factor_ij_1s*int_env*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
tmp += coef * coef_fit * int_fit
enddo
enddo
int2_u2_j1b2_test(j,i,ipoint) = tmp
int2_u2_env2_test(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -302,23 +302,23 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_u2_j1b2_test(j,i,ipoint) = int2_u2_j1b2_test(i,j,ipoint)
int2_u2_env2_test(j,i,ipoint) = int2_u2_env2_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u2_j1b2_test', wall1 - wall0
print*, ' wall time for int2_u2_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_PROVIDER [double precision, int2_u_grad1u_x_env2_test, (ao_num,ao_num,n_points_final_grid,3)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
!
END_DOC
@ -327,27 +327,27 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
double precision :: r(3), int_fit(3), expo_fit, coef_fit
double precision :: coef, beta, B_center(3), dist
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, coef_tmp
double precision :: tmp_x, tmp_y, tmp_z, int_j1b
double precision :: tmp_x, tmp_y, tmp_z, int_env
double precision :: wall0, wall1, sq_pi_3_2,sq_alpha
print*, ' providing int2_u_grad1u_x_j1b2_test ...'
print*, ' providing int2_u_grad1u_x_env2_test ...'
sq_pi_3_2 = dacos(-1.D0)**(1.d0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf final_grid_points
call wall_time(wall0)
int2_u_grad1u_x_j1b2_test = 0.d0
int2_u_grad1u_x_env2_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, &
!$OMP tmp_x, tmp_y, tmp_z,int_j1b,sq_alpha) &
!$OMP tmp_x, tmp_y, tmp_z,int_env,sq_alpha) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2,thrsh_cycle_tc)
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_x_env2_test,ao_abs_comb_b3_env,sq_pi_3_2,thrsh_cycle_tc)
!$OMP DO
do ipoint = 1, n_points_final_grid
@ -365,8 +365,8 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
int_env = ao_abs_comb_b3_env(i_1s,j,i)
if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -389,7 +389,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv)
! if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle
! if(dabs(coef_tmp*int_env*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle
call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
@ -402,9 +402,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
enddo
int2_u_grad1u_x_j1b2_test(j,i,ipoint,1) = tmp_x
int2_u_grad1u_x_j1b2_test(j,i,ipoint,2) = tmp_y
int2_u_grad1u_x_j1b2_test(j,i,ipoint,3) = tmp_z
int2_u_grad1u_x_env2_test(j,i,ipoint,1) = tmp_x
int2_u_grad1u_x_env2_test(j,i,ipoint,2) = tmp_y
int2_u_grad1u_x_env2_test(j,i,ipoint,3) = tmp_z
enddo
enddo
enddo
@ -414,24 +414,25 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_u_grad1u_x_j1b2_test(j,i,ipoint,1) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,1)
int2_u_grad1u_x_j1b2_test(j,i,ipoint,2) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,2)
int2_u_grad1u_x_j1b2_test(j,i,ipoint,3) = int2_u_grad1u_x_j1b2_test(i,j,ipoint,3)
int2_u_grad1u_x_env2_test(j,i,ipoint,1) = int2_u_grad1u_x_env2_test(i,j,ipoint,1)
int2_u_grad1u_x_env2_test(j,i,ipoint,2) = int2_u_grad1u_x_env2_test(i,j,ipoint,2)
int2_u_grad1u_x_env2_test(j,i,ipoint,3) = int2_u_grad1u_x_env2_test(i,j,ipoint,3)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_x_j1b2_test', wall1 - wall0
print*, ' wall time for int2_u_grad1u_x_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, int2_u_grad1u_env2_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu]
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu]
!
END_DOC
@ -442,31 +443,31 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
double precision :: alpha_1s, alpha_1s_inv, centr_1s(3), expo_coef_1s, tmp
double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s
double precision :: j12_mu_r12,int_j1b
double precision :: j12_mu_r12,int_env
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
double precision :: beta_ij,center_ij_1s(3),factor_ij_1s
print*, ' providing int2_u_grad1u_j1b2_test ...'
print*, ' providing int2_u_grad1u_env2_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen ao_overlap_abs List_comb_thr_b3_cent
provide mu_erf final_grid_points ao_overlap_abs List_comb_thr_b3_cent
call wall_time(wall0)
int2_u_grad1u_j1b2_test = 0.d0
int2_u_grad1u_env2_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
!$OMP beta_ij,center_ij_1s,factor_ij_1s, &
!$OMP int_j1b,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP int_env,alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b3_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2, &
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_j1b, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test,thrsh_cycle_tc)
!$OMP List_comb_thr_b3_coef, List_comb_thr_b3_expo, ao_abs_comb_b3_env, &
!$OMP List_comb_thr_b3_cent, int2_u_grad1u_env2_test,thrsh_cycle_tc)
!$OMP DO
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
@ -484,11 +485,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
! i_1s = 1
! --- --- ---
int_j1b = ao_abs_comb_b3_j1b(1,j,i)
! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
int_env = ao_abs_comb_b3_env(1,j,i)
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
! if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.thrsh_cycle_tc) cycle
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r)
tmp += coef_fit * int_fit
@ -502,8 +501,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
coef = List_comb_thr_b3_coef (i_1s,j,i)
beta = List_comb_thr_b3_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
int_env = ao_abs_comb_b3_env(i_1s,j,i)
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@ -513,7 +511,6 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
! if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.thrsh_cycle_tc)cycle
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
alpha_1s = beta + expo_fit
@ -533,7 +530,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
enddo
enddo
int2_u_grad1u_j1b2_test(j,i,ipoint) = tmp
int2_u_grad1u_env2_test(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -543,14 +540,15 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_u_grad1u_j1b2_test(j,i,ipoint) = int2_u_grad1u_j1b2_test(i,j,ipoint)
int2_u_grad1u_env2_test(j,i,ipoint) = int2_u_grad1u_env2_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_j1b2_test', wall1 - wall0
print*, ' wall time for int2_u_grad1u_env2_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---

153
src/ao_many_one_e_ints/grad2_jmu_modif.irp.f → plugins/local/ao_many_one_e_ints/grad2_jmu_modif.irp.f
View File

@ -21,7 +21,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin
print*, ' providing int2_grad1u2_grad2u2 ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf
provide final_grid_points
int2_grad1u2_grad2u2 = 0.d0
@ -63,17 +64,17 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_poin
enddo
call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2 =', wall1 - wall0
print*, ' wall time for int2_grad1u2_grad2u2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, int2_grad1u2_grad2u2_env2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [1 - erf(mu r12)]^2
! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [1 - erf(mu r12)]^2
!
END_DOC
@ -87,21 +88,22 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_grad1u2_grad2u2_j1b2 ...'
print*, ' providing int2_grad1u2_grad2u2_env2 ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf
provide final_grid_points
int2_grad1u2_grad2u2_j1b2 = 0.d0
int2_grad1u2_grad2u2_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_grad1u2_grad2u2_j1b2)
!$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_env1s_square_cent, int2_grad1u2_grad2u2_env2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
@ -125,14 +127,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
! ---
do i_1s = 2, List_all_comb_b3_size
do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s)
coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
beta = List_env1s_square_expo (i_1s)
B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_env1s_square_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -143,7 +145,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
enddo
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = tmp
int2_grad1u2_grad2u2_env2(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -153,23 +155,23 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_grad1u2_grad2u2_j1b2(j,i,ipoint) = int2_grad1u2_grad2u2_j1b2(i,j,ipoint)
int2_grad1u2_grad2u2_env2(j,i,ipoint) = int2_grad1u2_grad2u2_env2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_grad1u2_grad2u2_j1b2 =', wall1 - wall0
print*, ' wall time for int2_grad1u2_grad2u2_env2 (min) =', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, int2_u2_env2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 [u_12^mu]^2
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 [u_12^mu]^2
!
END_DOC
@ -182,21 +184,22 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing int2_u2_j1b2 ...'
print*, ' providing int2_u2_env2 ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf
provide final_grid_points
int2_u2_j1b2 = 0.d0
int2_u2_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u2_j1b2)
!$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_env1s_square_cent, int2_u2_env2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
@ -220,14 +223,14 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
! ---
do i_1s = 2, List_all_comb_b3_size
do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s)
coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
beta = List_env1s_square_expo (i_1s)
B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_env1s_square_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -238,7 +241,7 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
enddo
int2_u2_j1b2(j,i,ipoint) = tmp
int2_u2_env2(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -248,23 +251,23 @@ BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_u2_j1b2(j,i,ipoint) = int2_u2_j1b2(i,j,ipoint)
int2_u2_env2(j,i,ipoint) = int2_u2_env2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u2_j1b2', wall1 - wall0
print*, ' wall time for int2_u2_env2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_PROVIDER [double precision, int2_u_grad1u_x_env2, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu] r2
!
END_DOC
@ -276,23 +279,24 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1
print*, ' providing int2_u_grad1u_x_j1b2 ...'
print*, ' providing int2_u_grad1u_x_env2 ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf
provide final_grid_points
int2_u_grad1u_x_j1b2 = 0.d0
int2_u_grad1u_x_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp, &
!$OMP tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_x_j1b2)
!$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_env1s_square_cent, int2_u_grad1u_x_env2)
!$OMP DO
do ipoint = 1, n_points_final_grid
@ -321,14 +325,14 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
! ---
do i_1s = 2, List_all_comb_b3_size
do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s)
coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
beta = List_env1s_square_expo (i_1s)
B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_env1s_square_cent(3,i_1s)
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3))
@ -355,9 +359,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
enddo
int2_u_grad1u_x_j1b2(j,i,ipoint,1) = tmp_x
int2_u_grad1u_x_j1b2(j,i,ipoint,2) = tmp_y
int2_u_grad1u_x_j1b2(j,i,ipoint,3) = tmp_z
int2_u_grad1u_x_env2(j,i,ipoint,1) = tmp_x
int2_u_grad1u_x_env2(j,i,ipoint,2) = tmp_y
int2_u_grad1u_x_env2(j,i,ipoint,3) = tmp_z
enddo
enddo
enddo
@ -367,25 +371,25 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2, (ao_num, ao_num, n_poin
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_u_grad1u_x_j1b2(j,i,ipoint,1) = int2_u_grad1u_x_j1b2(i,j,ipoint,1)
int2_u_grad1u_x_j1b2(j,i,ipoint,2) = int2_u_grad1u_x_j1b2(i,j,ipoint,2)
int2_u_grad1u_x_j1b2(j,i,ipoint,3) = int2_u_grad1u_x_j1b2(i,j,ipoint,3)
int2_u_grad1u_x_env2(j,i,ipoint,1) = int2_u_grad1u_x_env2(i,j,ipoint,1)
int2_u_grad1u_x_env2(j,i,ipoint,2) = int2_u_grad1u_x_env2(i,j,ipoint,2)
int2_u_grad1u_x_env2(j,i,ipoint,3) = int2_u_grad1u_x_env2(i,j,ipoint,3)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_x_j1b2 = ', wall1 - wall0
print*, ' wall time for int2_u_grad1u_x_env2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [ double precision, int2_u_grad1u_env2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2)^2 u_12^mu [\grad_1 u_12^mu]
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2)^2 u_12^mu [\grad_1 u_12^mu]
!
END_DOC
@ -397,22 +401,23 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing int2_u_grad1u_j1b2 ...'
print*, ' providing int2_u_grad1u_env2 ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf
provide final_grid_points
int2_u_grad1u_j1b2 = 0.d0
int2_u_grad1u_env2 = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp, alpha_1s, dist, &
!$OMP alpha_1s_inv, centr_1s, expo_coef_1s, coef_tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b3_size, &
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_square_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf, &
!$OMP List_all_comb_b3_coef, List_all_comb_b3_expo, &
!$OMP List_all_comb_b3_cent, int2_u_grad1u_j1b2)
!$OMP List_env1s_square_coef, List_env1s_square_expo, &
!$OMP List_env1s_square_cent, int2_u_grad1u_env2)
!$OMP DO
do ipoint = 1, n_points_final_grid
do i = 1, ao_num
@ -436,14 +441,14 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
! ---
do i_1s = 2, List_all_comb_b3_size
do i_1s = 2, List_env1s_square_size
coef = List_all_comb_b3_coef (i_1s)
coef = List_env1s_square_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b3_expo (i_1s)
B_center(1) = List_all_comb_b3_cent(1,i_1s)
B_center(2) = List_all_comb_b3_cent(2,i_1s)
B_center(3) = List_all_comb_b3_cent(3,i_1s)
beta = List_env1s_square_expo (i_1s)
B_center(1) = List_env1s_square_cent(1,i_1s)
B_center(2) = List_env1s_square_cent(2,i_1s)
B_center(3) = List_env1s_square_cent(3,i_1s)
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
+ (B_center(3) - r(3)) * (B_center(3) - r(3))
@ -468,7 +473,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
enddo
int2_u_grad1u_j1b2(j,i,ipoint) = tmp
int2_u_grad1u_env2(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -478,13 +483,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2, (ao_num, ao_num, n_points
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
int2_u_grad1u_j1b2(j,i,ipoint) = int2_u_grad1u_j1b2(i,j,ipoint)
int2_u_grad1u_env2(j,i,ipoint) = int2_u_grad1u_env2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for int2_u_grad1u_j1b2', wall1 - wall0
print*, ' wall time for int2_u_grad1u_env2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER

115
src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f → plugins/local/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f
View File

@ -1,11 +1,11 @@
! ---
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_erf_rk_cst_mu_env_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R|
! int dr phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R| - 1) / |r - R|
!
END_DOC
@ -13,24 +13,23 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
integer :: i, j, ipoint, i_1s
double precision :: r(3), int_mu, int_coulomb
double precision :: coef, beta, B_center(3)
double precision :: tmp,int_j1b
double precision :: tmp,int_env
double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2
print*, ' providing v_ij_erf_rk_cst_mu_j1b_test ...'
print*, ' providing v_ij_erf_rk_cst_mu_env_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
v_ij_erf_rk_cst_mu_j1b_test = 0.d0
v_ij_erf_rk_cst_mu_env_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_j1b)&
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp, int_env)&
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b, &
!$OMP v_ij_erf_rk_cst_mu_j1b_test, mu_erf, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_env, &
!$OMP v_ij_erf_rk_cst_mu_env_test, mu_erf, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO
!do ipoint = 1, 10
@ -48,8 +47,8 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -60,7 +59,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
tmp += coef * (int_mu - int_coulomb)
enddo
v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = tmp
v_ij_erf_rk_cst_mu_env_test(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -70,22 +69,22 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint)
v_ij_erf_rk_cst_mu_env_test(j,i,ipoint) = v_ij_erf_rk_cst_mu_env_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0
print*, ' wall time for v_ij_erf_rk_cst_mu_env_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_PROVIDER [double precision, x_v_ij_erf_rk_cst_mu_env_test, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_DOC
! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
! int dr x phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R|) - 1)/|r - R|
END_DOC
implicit none
@ -93,23 +92,23 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
double precision :: coef, beta, B_center(3), r(3), ints(3), ints_coulomb(3)
double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b,factor_ij_1s,beta_ij,center_ij_1s
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_env,factor_ij_1s,beta_ij,center_ij_1s
print*, ' providing x_v_ij_erf_rk_cst_mu_j1b_test ...'
print*, ' providing x_v_ij_erf_rk_cst_mu_env_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide expo_erfc_mu_gauss ao_prod_sigma ao_prod_center
call wall_time(wall0)
x_v_ij_erf_rk_cst_mu_j1b_test = 0.d0
x_v_ij_erf_rk_cst_mu_env_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, &
!$OMP int_j1b, tmp_x, tmp_y, tmp_z,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP int_env, tmp_x, tmp_y, tmp_z,factor_ij_1s,beta_ij,center_ij_1s) &
!$OMP SHARED (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,&
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent, &
!$OMP x_v_ij_erf_rk_cst_mu_j1b_test, mu_erf,ao_abs_comb_b2_j1b, &
!$OMP x_v_ij_erf_rk_cst_mu_env_test, mu_erf,ao_abs_comb_b2_env, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,thrsh_cycle_tc)
! !$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss)
!$OMP DO
@ -129,8 +128,8 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -143,9 +142,9 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
tmp_z += coef * (ints(3) - ints_coulomb(3))
enddo
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = tmp_x
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = tmp_y
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = tmp_z
x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,1) = tmp_x
x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,2) = tmp_y
x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,3) = tmp_z
enddo
enddo
enddo
@ -155,26 +154,26 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,1)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,2)
x_v_ij_erf_rk_cst_mu_j1b_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,3)
x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_env_test(i,j,ipoint,1)
x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_env_test(i,j,ipoint,2)
x_v_ij_erf_rk_cst_mu_env_test(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_env_test(i,j,ipoint,3)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b_test', wall1 - wall0
print*, ' wall time for x_v_ij_erf_rk_cst_mu_env_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
! TODO analytically
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
!
END_DOC
@ -185,29 +184,28 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
double precision :: tmp
double precision :: wall0, wall1
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_env
double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_test ...'
print*, ' providing v_ij_u_cst_mu_env_test ...'
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
v_ij_u_cst_mu_j1b_test = 0.d0
v_ij_u_cst_mu_env_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_env) &
!$OMP SHARED (n_points_final_grid, ao_num, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_env_test,ao_abs_comb_b2_env, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO
do ipoint = 1, n_points_final_grid
@ -225,8 +223,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
! i_1s = 1
! --- --- ---
int_j1b = ao_abs_comb_b2_j1b(1,j,i)
! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
int_env = ao_abs_comb_b2_env(1,j,i)
! if(dabs(int_env).lt.thrsh_cycle_tc) cycle
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_j_mu_x(i_fit)
coef_fit = coef_gauss_j_mu_x(i_fit)
@ -242,8 +240,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
do i_1s = 2, List_comb_thr_b2_size(j,i)
coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -259,7 +257,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
enddo
enddo
v_ij_u_cst_mu_j1b_test(j,i,ipoint) = tmp
v_ij_u_cst_mu_env_test(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -269,23 +267,23 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_u_cst_mu_j1b_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_test(i,j,ipoint)
v_ij_u_cst_mu_env_test(j,i,ipoint) = v_ij_u_cst_mu_env_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_test', wall1 - wall0
print*, ' wall time for v_ij_u_cst_mu_env_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_ng_1_test, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12) with u(mu,r12) \approx 1/2 mu e^{-2.5 * mu (r12)^2}
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12) with u(mu,r12) \approx 1/2 mu e^{-2.5 * mu (r12)^2}
!
END_DOC
@ -296,27 +294,26 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
double precision :: tmp
double precision :: wall0, wall1
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_env
double precision, external :: overlap_gauss_r12_ao
double precision, external :: overlap_gauss_r12_ao_with1s
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
provide mu_erf final_grid_points j1b_pen
call wall_time(wall0)
v_ij_u_cst_mu_j1b_ng_1_test = 0.d0
v_ij_u_cst_mu_env_ng_1_test = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_env) &
!$OMP SHARED (n_points_final_grid, ao_num, &
!$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, &
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_env_ng_1_test,ao_abs_comb_b2_env, &
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
!$OMP DO
do ipoint = 1, n_points_final_grid
@ -334,8 +331,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
! i_1s = 1
! --- --- ---
int_j1b = ao_abs_comb_b2_j1b(1,j,i)
! if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
int_env = ao_abs_comb_b2_env(1,j,i)
! if(dabs(int_env).lt.thrsh_cycle_tc) cycle
expo_fit = expo_good_j_mu_1gauss
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
tmp += int_fit
@ -347,8 +344,8 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
do i_1s = 2, List_comb_thr_b2_size(j,i)
coef = List_comb_thr_b2_coef (i_1s,j,i)
beta = List_comb_thr_b2_expo (i_1s,j,i)
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
! if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
int_env = ao_abs_comb_b2_env(i_1s,j,i)
! if(dabs(coef)*dabs(int_env).lt.thrsh_cycle_tc)cycle
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@ -364,7 +361,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
! enddo
enddo
v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint) = tmp
v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -374,13 +371,13 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_u_cst_mu_j1b_ng_1_test(j,i,ipoint) = v_ij_u_cst_mu_j1b_ng_1_test(i,j,ipoint)
v_ij_u_cst_mu_env_ng_1_test(j,i,ipoint) = v_ij_u_cst_mu_env_ng_1_test(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_ng_1_test', wall1 - wall0
print*, ' wall time for v_ij_u_cst_mu_env_ng_1_test (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER

237
src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f → plugins/local/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f
View File

@ -1,11 +1,11 @@
! ---
BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_erf_rk_cst_mu_env, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R| - 1) / |r - R|
! int dr phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R| - 1) / |r - R|
!
END_DOC
@ -17,18 +17,20 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
double precision :: wall0, wall1
double precision, external :: NAI_pol_mult_erf_ao_with1s
print *, ' providing v_ij_erf_rk_cst_mu_j1b ...'
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE env_expo
print *, ' providing v_ij_erf_rk_cst_mu_env ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
v_ij_erf_rk_cst_mu_j1b = 0.d0
v_ij_erf_rk_cst_mu_env = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, int_mu, int_coulomb, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, &
!$OMP v_ij_erf_rk_cst_mu_j1b, mu_erf)
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, final_grid_points, &
!$OMP List_env1s_coef, List_env1s_expo, List_env1s_cent, &
!$OMP v_ij_erf_rk_cst_mu_env, mu_erf)
!$OMP DO
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
@ -43,28 +45,27 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
! ---
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
coef = List_env1s_coef (1)
beta = List_env1s_expo (1)
B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_env1s_cent(3,1)
int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
! if(dabs(coef)*dabs(int_mu - int_coulomb) .lt. 1d-12) cycle
tmp += coef * (int_mu - int_coulomb)
! ---
do i_1s = 2, List_all_comb_b2_size
do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s)
coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
beta = List_env1s_expo (i_1s)
B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_env1s_cent(3,i_1s)
int_mu = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
int_coulomb = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
@ -74,7 +75,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
! ---
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = tmp
v_ij_erf_rk_cst_mu_env(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -84,22 +85,22 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_po
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_erf_rk_cst_mu_j1b(j,i,ipoint) = v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
v_ij_erf_rk_cst_mu_env(j,i,ipoint) = v_ij_erf_rk_cst_mu_env(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_erf_rk_cst_mu_j1b', wall1 - wall0
print*, ' wall time for v_ij_erf_rk_cst_mu_env (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_PROVIDER [double precision, x_v_ij_erf_rk_cst_mu_env, (ao_num, ao_num, n_points_final_grid, 3)]
BEGIN_DOC
! int dr x phi_i(r) phi_j(r) 1s_j1b(r) (erf(mu(R) |r - R|) - 1)/|r - R|
! int dr x phi_i(r) phi_j(r) 1s_env(r) (erf(mu(R) |r - R|) - 1)/|r - R|
END_DOC
implicit none
@ -108,17 +109,17 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
double precision :: tmp_x, tmp_y, tmp_z
double precision :: wall0, wall1
print*, ' providing x_v_ij_erf_rk_cst_mu_j1b ...'
print*, ' providing x_v_ij_erf_rk_cst_mu_env ...'
call wall_time(wall0)
x_v_ij_erf_rk_cst_mu_j1b = 0.d0
x_v_ij_erf_rk_cst_mu_env = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, ints, ints_coulomb, &
!$OMP tmp_x, tmp_y, tmp_z) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, final_grid_points,&
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, List_all_comb_b2_cent, &
!$OMP x_v_ij_erf_rk_cst_mu_j1b, mu_erf)
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, final_grid_points,&
!$OMP List_env1s_coef, List_env1s_expo, List_env1s_cent, &
!$OMP x_v_ij_erf_rk_cst_mu_env, mu_erf)
!$OMP DO
!do ipoint = 1, 10
do ipoint = 1, n_points_final_grid
@ -135,11 +136,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
! ---
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
coef = List_env1s_coef (1)
beta = List_env1s_expo (1)
B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_env1s_cent(3,1)
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
@ -152,14 +153,14 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
! ---
do i_1s = 2, List_all_comb_b2_size
do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s)
coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
beta = List_env1s_expo (i_1s)
B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_env1s_cent(3,i_1s)
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints )
call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, ints_coulomb)
@ -171,9 +172,9 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
! ---
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = tmp_x
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = tmp_y
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = tmp_z
x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,1) = tmp_x
x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,2) = tmp_y
x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,3) = tmp_z
enddo
enddo
enddo
@ -183,25 +184,25 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b, (ao_num, ao_num, n_
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1)
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2)
x_v_ij_erf_rk_cst_mu_j1b(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3)
x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,1) = x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,1)
x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,2) = x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,2)
x_v_ij_erf_rk_cst_mu_env(j,i,ipoint,3) = x_v_ij_erf_rk_cst_mu_env(i,j,ipoint,3)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for x_v_ij_erf_rk_cst_mu_j1b =', wall1 - wall0
print*, ' wall time for x_v_ij_erf_rk_cst_mu_env (min) =', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_fit, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
!
END_DOC
@ -214,23 +215,23 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
double precision, external :: overlap_gauss_r12_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_fit ...'
print*, ' providing v_ij_u_cst_mu_env_fit ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
provide mu_erf final_grid_points env_expo
PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
PROVIDE List_env1s_size List_env1s_coef List_env1s_expo List_env1s_cent
v_ij_u_cst_mu_j1b_fit = 0.d0
v_ij_u_cst_mu_env_fit = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
!$OMP coef_fit, expo_fit, int_fit, tmp) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, &
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, &
!$OMP final_grid_points, ng_fit_jast, &
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_fit)
!$OMP List_env1s_coef, List_env1s_expo, &
!$OMP List_env1s_cent, v_ij_u_cst_mu_env_fit)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
@ -247,11 +248,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
! ---
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
coef = List_env1s_coef (1)
beta = List_env1s_expo (1)
B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_env1s_cent(3,1)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -259,14 +260,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
! ---
do i_1s = 2, List_all_comb_b2_size
do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s)
coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
beta = List_env1s_expo (i_1s)
B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_env1s_cent(3,i_1s)
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@ -277,7 +278,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
enddo
v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = tmp
v_ij_u_cst_mu_env_fit(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -287,23 +288,23 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_poi
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = v_ij_u_cst_mu_j1b_fit(i,j,ipoint)
v_ij_u_cst_mu_env_fit(j,i,ipoint) = v_ij_u_cst_mu_env_fit(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_fit', wall1 - wall0
print*, ' wall time for v_ij_u_cst_mu_env_fit (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_an_old, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
!
END_DOC
@ -322,24 +323,24 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
double precision, external :: overlap_gauss_r12_ao_with1s
double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_an_old ...'
print*, ' providing v_ij_u_cst_mu_env_an_old ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
provide mu_erf final_grid_points env_expo
PROVIDE List_env1s_size List_env1s_coef List_env1s_expo List_env1s_cent
ct = inv_sq_pi_2 / mu_erf
v_ij_u_cst_mu_j1b_an_old = 0.d0
v_ij_u_cst_mu_env_an_old = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP r1_2, tmp, int_c1, int_e1, int_o, int_c2, &
!$OMP int_e2, int_c3, int_e3) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, &
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, &
!$OMP final_grid_points, mu_erf, ct, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an_old)
!$OMP List_env1s_coef, List_env1s_expo, &
!$OMP List_env1s_cent, v_ij_u_cst_mu_env_an_old)
!$OMP DO
do ipoint = 1, n_points_final_grid
@ -353,11 +354,11 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
! ---
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
coef = List_env1s_coef (1)
beta = List_env1s_expo (1)
B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_env1s_cent(3,1)
int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
@ -379,14 +380,14 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
! ---
do i_1s = 2, List_all_comb_b2_size
do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s)
coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
beta = List_env1s_expo (i_1s)
B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_env1s_cent(3,i_1s)
int_c1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r)
int_e1 = NAI_pol_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r)
@ -410,7 +411,7 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
! ---
v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = tmp
v_ij_u_cst_mu_env_an_old(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -420,23 +421,23 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an_old, (ao_num, ao_num, n_p
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_u_cst_mu_j1b_an_old(j,i,ipoint) = v_ij_u_cst_mu_j1b_an_old(i,j,ipoint)
v_ij_u_cst_mu_env_an_old(j,i,ipoint) = v_ij_u_cst_mu_env_an_old(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_an_old', wall1 - wall0
print*, ' wall time for v_ij_u_cst_mu_env_an_old (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_points_final_grid)]
BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_env_an, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! int dr2 phi_i(r2) phi_j(r2) 1s_j1b(r2) u(mu, r12)
! int dr2 phi_i(r2) phi_j(r2) 1s_env(r2) u(mu, r12)
!
END_DOC
@ -454,23 +455,23 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
double precision, external :: overlap_gauss_r12_ao_with1s
double precision, external :: NAI_pol_mult_erf_ao_with1s
print*, ' providing v_ij_u_cst_mu_j1b_an ...'
print*, ' providing v_ij_u_cst_mu_env_an ...'
call wall_time(wall0)
provide mu_erf final_grid_points j1b_pen
PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
provide mu_erf final_grid_points env_expo
PROVIDE List_env1s_size List_env1s_coef List_env1s_expo List_env1s_cent
ct = inv_sq_pi_2 / mu_erf
v_ij_u_cst_mu_j1b_an = 0.d0
v_ij_u_cst_mu_env_an = 0.d0
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
!$OMP r1_2, tmp, int_c, int_e, int_o) &
!$OMP SHARED (n_points_final_grid, ao_num, List_all_comb_b2_size, &
!$OMP SHARED (n_points_final_grid, ao_num, List_env1s_size, &
!$OMP final_grid_points, mu_erf, ct, &
!$OMP List_all_comb_b2_coef, List_all_comb_b2_expo, &
!$OMP List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an)
!$OMP List_env1s_coef, List_env1s_expo, &
!$OMP List_env1s_cent, v_ij_u_cst_mu_env_an)
!$OMP DO
do ipoint = 1, n_points_final_grid
@ -484,11 +485,11 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
! ---
coef = List_all_comb_b2_coef (1)
beta = List_all_comb_b2_expo (1)
B_center(1) = List_all_comb_b2_cent(1,1)
B_center(2) = List_all_comb_b2_cent(2,1)
B_center(3) = List_all_comb_b2_cent(3,1)
coef = List_env1s_coef (1)
beta = List_env1s_expo (1)
B_center(1) = List_env1s_cent(1,1)
B_center(2) = List_env1s_cent(2,1)
B_center(3) = List_env1s_cent(3,1)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e)
@ -504,14 +505,14 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
! ---
do i_1s = 2, List_all_comb_b2_size
do i_1s = 2, List_env1s_size
coef = List_all_comb_b2_coef (i_1s)
coef = List_env1s_coef (i_1s)
if(dabs(coef) .lt. 1d-15) cycle ! beta = 0.0
beta = List_all_comb_b2_expo (i_1s)
B_center(1) = List_all_comb_b2_cent(1,i_1s)
B_center(2) = List_all_comb_b2_cent(2,i_1s)
B_center(3) = List_all_comb_b2_cent(3,i_1s)
beta = List_env1s_expo (i_1s)
B_center(1) = List_env1s_cent(1,i_1s)
B_center(2) = List_env1s_cent(2,i_1s)
B_center(3) = List_env1s_cent(3,i_1s)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_c)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_e)
@ -529,7 +530,7 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
! ---
v_ij_u_cst_mu_j1b_an(j,i,ipoint) = tmp
v_ij_u_cst_mu_env_an(j,i,ipoint) = tmp
enddo
enddo
enddo
@ -539,13 +540,13 @@ BEGIN_PROVIDER [double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_point
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
v_ij_u_cst_mu_j1b_an(j,i,ipoint) = v_ij_u_cst_mu_j1b_an(i,j,ipoint)
v_ij_u_cst_mu_env_an(j,i,ipoint) = v_ij_u_cst_mu_env_an(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for v_ij_u_cst_mu_j1b_an', wall1 - wall0
print*, ' wall time for v_ij_u_cst_mu_env_an (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER

0
src/ao_many_one_e_ints/grad_related_ints.irp.f → plugins/local/ao_many_one_e_ints/grad_related_ints.irp.f
View File

574
plugins/local/ao_many_one_e_ints/lin_fc_rsdft.irp.f Normal file
View File

@ -0,0 +1,574 @@
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_0, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_x, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_y, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_z, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du_2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_long_Du_0 = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12]
!
! Ir2_Mu_long_Du_x = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * x2
! Ir2_Mu_long_Du_y = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * y2
! Ir2_Mu_long_Du_z = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * z2
!
! Ir2_Mu_long_Du_2 = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) [(1 - erf(mu r_12) / r_12] * r2^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3), int_clb(7), int_erf(7)
double precision :: c_1s, e_1s, R_1s(3)
double precision :: tmp_Du_0, tmp_Du_x, tmp_Du_y, tmp_Du_z, tmp_Du_2
double precision :: wall0, wall1
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_size List_env1s_expo List_env1s_coef List_env1s_cent
print *, ' providing Ir2_Mu_long_Du ...'
call wall_time(wall0)
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, c_1s, e_1s, R_1s, int_erf, int_clb, &
!$OMP tmp_Du_0, tmp_Du_x, tmp_Du_y, tmp_Du_z, tmp_Du_2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_erf, &
!$OMP List_env1s_size, List_env1s_expo, &
!$OMP List_env1s_coef, List_env1s_cent, &
!$OMP Ir2_Mu_long_Du_0, Ir2_Mu_long_Du_x, &
!$OMP Ir2_Mu_long_Du_y, Ir2_Mu_long_Du_z, &
!$OMP Ir2_Mu_long_Du_2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
do i = 1, ao_num
do j = i, ao_num
call NAI_pol_012_mult_erf_ao(i, j, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao(i, j, mu_erf, r, int_erf)
tmp_Du_0 = int_clb(1) - int_erf(1)
tmp_Du_x = int_clb(2) - int_erf(2)
tmp_Du_y = int_clb(3) - int_erf(3)
tmp_Du_z = int_clb(4) - int_erf(4)
tmp_Du_2 = int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7)
do i_1s = 2, List_env1s_size
e_1s = List_env1s_expo(i_1s)
c_1s = List_env1s_coef(i_1s)
R_1s(1) = List_env1s_cent(1,i_1s)
R_1s(2) = List_env1s_cent(2,i_1s)
R_1s(3) = List_env1s_cent(3,i_1s)
call NAI_pol_012_mult_erf_ao_with1s(i, j, e_1s, R_1s, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao_with1s(i, j, e_1s, R_1s, mu_erf, r, int_erf)
tmp_Du_0 = tmp_Du_0 + c_1s * (int_clb(1) - int_erf(1))
tmp_Du_x = tmp_Du_x + c_1s * (int_clb(2) - int_erf(2))
tmp_Du_y = tmp_Du_y + c_1s * (int_clb(3) - int_erf(3))
tmp_Du_z = tmp_Du_z + c_1s * (int_clb(4) - int_erf(4))
tmp_Du_2 = tmp_Du_2 + c_1s * (int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7))
enddo
Ir2_Mu_long_Du_0(j,i,ipoint) = tmp_Du_0
Ir2_Mu_long_Du_x(j,i,ipoint) = tmp_Du_x
Ir2_Mu_long_Du_y(j,i,ipoint) = tmp_Du_y
Ir2_Mu_long_Du_z(j,i,ipoint) = tmp_Du_z
Ir2_Mu_long_Du_2(j,i,ipoint) = tmp_Du_2
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_long_Du_0(j,i,ipoint) = Ir2_Mu_long_Du_0(i,j,ipoint)
Ir2_Mu_long_Du_x(j,i,ipoint) = Ir2_Mu_long_Du_x(i,j,ipoint)
Ir2_Mu_long_Du_y(j,i,ipoint) = Ir2_Mu_long_Du_y(i,j,ipoint)
Ir2_Mu_long_Du_z(j,i,ipoint) = Ir2_Mu_long_Du_z(i,j,ipoint)
Ir2_Mu_long_Du_2(j,i,ipoint) = Ir2_Mu_long_Du_2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_long_Du (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_gauss_Du, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_gauss_Du = int dr2 phi_i(r2) phi_j(r2) fc_env(r2) e^{-(mu r_12)^2}
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du
double precision :: mu_sq, dx, dy, dz, tmp_arg, rmu_sq(3)
double precision :: e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
double precision, external :: overlap_gauss_r12_ao
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_size List_env1s_expo List_env1s_coef List_env1s_cent
print *, ' providing Ir2_Mu_gauss_Du ...'
call wall_time(wall0)
mu_sq = mu_erf * mu_erf
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, dx, dy, dz, r, tmp_arg, coef, &
!$OMP rmu_sq, e_1s, c_1s, R_1s, beta, B_center, tmp_Du) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_sq, &
!$OMP List_env1s_size, List_env1s_expo, &
!$OMP List_env1s_coef, List_env1s_cent, &
!$OMP Ir2_Mu_gauss_Du)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
rmu_sq(1) = mu_sq * r(1)
rmu_sq(2) = mu_sq * r(2)
rmu_sq(3) = mu_sq * r(3)
do i = 1, ao_num
do j = i, ao_num
tmp_Du = overlap_gauss_r12_ao(r, mu_sq, j, i)
do i_1s = 2, List_env1s_size
e_1s = List_env1s_expo(i_1s)
c_1s = List_env1s_coef(i_1s)
R_1s(1) = List_env1s_cent(1,i_1s)
R_1s(2) = List_env1s_cent(2,i_1s)
R_1s(3) = List_env1s_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = mu_sq + e_1s
tmp_arg = mu_sq * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = c_1s * dexp(-tmp_arg)
B_center(1) = (rmu_sq(1) + e_1s * R_1s(1)) / beta
B_center(2) = (rmu_sq(2) + e_1s * R_1s(2)) / beta
B_center(3) = (rmu_sq(3) + e_1s * R_1s(3)) / beta
tmp_Du += coef * overlap_gauss_r12_ao(B_center, beta, j, i)
enddo
Ir2_Mu_gauss_Du(j,i,ipoint) = tmp_Du
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_gauss_Du(j,i,ipoint) = Ir2_Mu_gauss_Du(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_gauss_Du (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_0, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_x, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_y, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_z, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_long_Du2_2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_long_Du2_0 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12]
!
! Ir2_Mu_long_Du2_x = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * x2
! Ir2_Mu_long_Du2_y = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * y2
! Ir2_Mu_long_Du2_z = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * z2
!
! Ir2_Mu_long_Du2_2 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12) / r_12] * r2^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3), int_clb(7), int_erf(7)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2
double precision :: mu_sq, tmp_arg, dx, dy, dz, rmu_sq(3)
double precision :: e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_square_size List_env1s_square_expo List_env1s_square_coef List_env1s_square_cent
print *, ' providing Ir2_Mu_long_Du2 ...'
call wall_time(wall0)
mu_sq = mu_erf * mu_erf
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, r, rmu_sq, dx, dy, dz, &
!$OMP e_1s, c_1s, R_1s, tmp_arg, coef, beta, B_center, &
!$OMP int_erf, int_clb, &
!$OMP tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_sq, &
!$OMP mu_erf, List_env1s_square_size, List_env1s_square_expo, &
!$OMP List_env1s_square_coef, List_env1s_square_cent, &
!$OMP Ir2_Mu_long_Du2_0, Ir2_Mu_long_Du2_x, &
!$OMP Ir2_Mu_long_Du2_y, Ir2_Mu_long_Du2_z, &
!$OMP Ir2_Mu_long_Du2_2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
rmu_sq(1) = mu_sq * r(1)
rmu_sq(2) = mu_sq * r(2)
rmu_sq(3) = mu_sq * r(3)
do i = 1, ao_num
do j = i, ao_num
call NAI_pol_012_mult_erf_ao_with1s(i, j, mu_sq, r, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao_with1s(i, j, mu_sq, r, mu_erf, r, int_erf)
tmp_Du2_0 = int_clb(1) - int_erf(1)
tmp_Du2_x = int_clb(2) - int_erf(2)
tmp_Du2_y = int_clb(3) - int_erf(3)
tmp_Du2_z = int_clb(4) - int_erf(4)
tmp_Du2_2 = int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7)
do i_1s = 2, List_env1s_square_size
e_1s = List_env1s_square_expo(i_1s)
c_1s = List_env1s_square_coef(i_1s)
R_1s(1) = List_env1s_square_cent(1,i_1s)
R_1s(2) = List_env1s_square_cent(2,i_1s)
R_1s(3) = List_env1s_square_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = mu_sq + e_1s
tmp_arg = mu_sq * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = c_1s * dexp(-tmp_arg)
B_center(1) = (rmu_sq(1) + e_1s * R_1s(1)) / beta
B_center(2) = (rmu_sq(2) + e_1s * R_1s(2)) / beta
B_center(3) = (rmu_sq(3) + e_1s * R_1s(3)) / beta
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, 1.d+9, r, int_clb)
call NAI_pol_012_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, int_erf)
tmp_Du2_0 = tmp_Du2_0 + coef * (int_clb(1) - int_erf(1))
tmp_Du2_x = tmp_Du2_x + coef * (int_clb(2) - int_erf(2))
tmp_Du2_y = tmp_Du2_y + coef * (int_clb(3) - int_erf(3))
tmp_Du2_z = tmp_Du2_z + coef * (int_clb(4) - int_erf(4))
tmp_Du2_2 = tmp_Du2_2 + coef * (int_clb(5) + int_clb(6) + int_clb(7) - int_erf(5) - int_erf(6) - int_erf(7))
enddo
Ir2_Mu_long_Du2_0(j,i,ipoint) = tmp_Du2_0
Ir2_Mu_long_Du2_x(j,i,ipoint) = tmp_Du2_x
Ir2_Mu_long_Du2_y(j,i,ipoint) = tmp_Du2_y
Ir2_Mu_long_Du2_z(j,i,ipoint) = tmp_Du2_z
Ir2_Mu_long_Du2_2(j,i,ipoint) = tmp_Du2_2
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_long_Du2_0(j,i,ipoint) = Ir2_Mu_long_Du2_0(i,j,ipoint)
Ir2_Mu_long_Du2_x(j,i,ipoint) = Ir2_Mu_long_Du2_x(i,j,ipoint)
Ir2_Mu_long_Du2_y(j,i,ipoint) = Ir2_Mu_long_Du2_y(i,j,ipoint)
Ir2_Mu_long_Du2_z(j,i,ipoint) = Ir2_Mu_long_Du2_z(i,j,ipoint)
Ir2_Mu_long_Du2_2(j,i,ipoint) = Ir2_Mu_long_Du2_2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_long_Du2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_gauss_Du2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_gauss_Du2 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 e^{-(mu r_12)^2}
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s
double precision :: r(3)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du2
double precision :: mu_sq, dx, dy, dz, tmp_arg, rmu_sq(3)
double precision :: e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
double precision, external :: overlap_gauss_r12_ao
PROVIDE mu_erf
PROVIDE final_grid_points
PROVIDE List_env1s_square_size List_env1s_square_expo List_env1s_square_coef List_env1s_square_cent
print *, ' providing Ir2_Mu_gauss_Du2 ...'
call wall_time(wall0)
mu_sq = 2.d0 * mu_erf * mu_erf
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, dx, dy, dz, r, tmp_arg, coef, &
!$OMP rmu_sq, e_1s, c_1s, R_1s, beta, B_center, tmp_Du2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, mu_sq, &
!$OMP List_env1s_square_size, List_env1s_square_expo, &
!$OMP List_env1s_square_coef, List_env1s_square_cent, &
!$OMP Ir2_Mu_gauss_Du2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
rmu_sq(1) = mu_sq * r(1)
rmu_sq(2) = mu_sq * r(2)
rmu_sq(3) = mu_sq * r(3)
do i = 1, ao_num
do j = i, ao_num
tmp_Du2 = overlap_gauss_r12_ao(r, mu_sq, j, i)
do i_1s = 2, List_env1s_square_size
e_1s = List_env1s_square_expo(i_1s)
c_1s = List_env1s_square_coef(i_1s)
R_1s(1) = List_env1s_square_cent(1,i_1s)
R_1s(2) = List_env1s_square_cent(2,i_1s)
R_1s(3) = List_env1s_square_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = mu_sq + e_1s
tmp_arg = mu_sq * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = c_1s * dexp(-tmp_arg)
B_center(1) = (rmu_sq(1) + e_1s * R_1s(1)) / beta
B_center(2) = (rmu_sq(2) + e_1s * R_1s(2)) / beta
B_center(3) = (rmu_sq(3) + e_1s * R_1s(3)) / beta
tmp_Du2 += coef * overlap_gauss_r12_ao(B_center, beta, j, i)
enddo
Ir2_Mu_gauss_Du2(j,i,ipoint) = tmp_Du2
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_gauss_Du2(j,i,ipoint) = Ir2_Mu_gauss_Du2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_gauss_Du2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_0, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_x, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_y, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_z, (ao_num, ao_num, n_points_final_grid)]
&BEGIN_PROVIDER [double precision, Ir2_Mu_short_Du2_2, (ao_num, ao_num, n_points_final_grid)]
BEGIN_DOC
!
! Ir2_Mu_short_Du2_0 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2
!
! Ir2_Mu_short_Du2_x = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * x2
! Ir2_Mu_short_Du2_y = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * y2
! Ir2_Mu_short_Du2_z = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * z2
!
! Ir2_Mu_short_Du2_2 = int dr2 phi_i(r2) phi_j(r2) [fc_env(r2)]^2 [(1 - erf(mu r_12)]^2 * r2^2
!
END_DOC
implicit none
integer :: i, j, ipoint, i_1s, i_fit
double precision :: r(3), ints(7)
double precision :: coef, beta, B_center(3)
double precision :: tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2
double precision :: tmp_arg, dx, dy, dz
double precision :: expo_fit, coef_fit, e_1s, c_1s, R_1s(3)
double precision :: wall0, wall1
PROVIDE final_grid_points
PROVIDE List_env1s_square_size List_env1s_square_expo List_env1s_square_coef List_env1s_square_cent
PROVIDE ng_fit_jast expo_gauss_1_erf_x_2 coef_gauss_1_erf_x_2
print *, ' providing Ir2_Mu_short_Du2 ...'
call wall_time(wall0)
!$OMP PARALLEL DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, dx, dy, dz, &
!$OMP expo_fit, coef_fit, e_1s, c_1s, R_1s, &
!$OMP tmp_arg, coef, beta, B_center, ints, &
!$OMP tmp_Du2_0, tmp_Du2_x, tmp_Du2_y, tmp_Du2_z, tmp_Du2_2) &
!$OMP SHARED (n_points_final_grid, ao_num, final_grid_points, &
!$OMP ng_fit_jast, expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2, &
!$OMP List_env1s_square_size, List_env1s_square_expo, &
!$OMP List_env1s_square_coef, List_env1s_square_cent, &
!$OMP Ir2_Mu_short_Du2_0, Ir2_Mu_short_Du2_x, &
!$OMP Ir2_Mu_short_Du2_y, Ir2_Mu_short_Du2_z, &
!$OMP Ir2_Mu_short_Du2_2)
!$OMP DO
do ipoint = 1, n_points_final_grid
r(1) = final_grid_points(1,ipoint)
r(2) = final_grid_points(2,ipoint)
r(3) = final_grid_points(3,ipoint)
do i = 1, ao_num
do j = i, ao_num
tmp_Du2_0 = 0.d0
tmp_Du2_x = 0.d0
tmp_Du2_y = 0.d0
tmp_Du2_z = 0.d0
tmp_Du2_2 = 0.d0
do i_fit = 1, ng_fit_jast
expo_fit = expo_gauss_1_erf_x_2(i_fit)
coef_fit = coef_gauss_1_erf_x_2(i_fit)
call overlap_gauss_r12_ao_012(r, expo_fit, i, j, ints)
tmp_Du2_0 += coef_fit * ints(1)
tmp_Du2_x += coef_fit * ints(2)
tmp_Du2_y += coef_fit * ints(3)
tmp_Du2_z += coef_fit * ints(4)
tmp_Du2_2 += coef_fit * (ints(5) + ints(6) + ints(7))
do i_1s = 2, List_env1s_square_size
e_1s = List_env1s_square_expo(i_1s)
c_1s = List_env1s_square_coef(i_1s)
R_1s(1) = List_env1s_square_cent(1,i_1s)
R_1s(2) = List_env1s_square_cent(2,i_1s)
R_1s(3) = List_env1s_square_cent(3,i_1s)
dx = r(1) - R_1s(1)
dy = r(2) - R_1s(2)
dz = r(3) - R_1s(3)
beta = expo_fit + e_1s
tmp_arg = expo_fit * e_1s * (dx*dx + dy*dy + dz*dz) / beta
coef = coef_fit * c_1s * dexp(-tmp_arg)
B_center(1) = (expo_fit * r(1) + e_1s * R_1s(1)) / beta
B_center(2) = (expo_fit * r(2) + e_1s * R_1s(2)) / beta
B_center(3) = (expo_fit * r(3) + e_1s * R_1s(3)) / beta
call overlap_gauss_r12_ao_012(B_center, beta, i, j, ints)
tmp_Du2_0 += coef * ints(1)
tmp_Du2_x += coef * ints(2)
tmp_Du2_y += coef * ints(3)
tmp_Du2_z += coef * ints(4)
tmp_Du2_2 += coef * (ints(5) + ints(6) + ints(7))
enddo ! i_1s
enddo ! i_fit
Ir2_Mu_short_Du2_0(j,i,ipoint) = tmp_Du2_0
Ir2_Mu_short_Du2_x(j,i,ipoint) = tmp_Du2_x
Ir2_Mu_short_Du2_y(j,i,ipoint) = tmp_Du2_y
Ir2_Mu_short_Du2_z(j,i,ipoint) = tmp_Du2_z
Ir2_Mu_short_Du2_2(j,i,ipoint) = tmp_Du2_2
enddo ! j
enddo ! i
enddo ! ipoint
!$OMP END DO
!$OMP END PARALLEL
do ipoint = 1, n_points_final_grid
do i = 2, ao_num
do j = 1, i-1
Ir2_Mu_short_Du2_0(j,i,ipoint) = Ir2_Mu_short_Du2_0(i,j,ipoint)
Ir2_Mu_short_Du2_x(j,i,ipoint) = Ir2_Mu_short_Du2_x(i,j,ipoint)
Ir2_Mu_short_Du2_y(j,i,ipoint) = Ir2_Mu_short_Du2_y(i,j,ipoint)
Ir2_Mu_short_Du2_z(j,i,ipoint) = Ir2_Mu_short_Du2_z(i,j,ipoint)
Ir2_Mu_short_Du2_2(j,i,ipoint) = Ir2_Mu_short_Du2_2(i,j,ipoint)
enddo
enddo
enddo
call wall_time(wall1)
print*, ' wall time for Ir2_Mu_short_Du2 (min) = ', (wall1 - wall0) / 60.d0
END_PROVIDER
! ---

0
src/ao_many_one_e_ints/list_grid.irp.f → plugins/local/ao_many_one_e_ints/list_grid.irp.f
View File

351
plugins/local/ao_many_one_e_ints/listj1b.irp.f Normal file
View File

@ -0,0 +1,351 @@
! ---
BEGIN_PROVIDER [integer, List_env1s_size]
implicit none
PROVIDE env_type
if(env_type .eq. "Prod_Gauss") then
List_env1s_size = 2**nucl_num
elseif(env_type .eq. "Sum_Gauss") then
List_env1s_size = nucl_num + 1
else
print *, ' Error in List_env1s_size: Unknown env_type = ', env_type
stop
endif
print *, ' nb of 1s-Gaussian in the envelope = ', List_env1s_size
END_PROVIDER
! ---
BEGIN_PROVIDER [integer, List_env1s, (nucl_num, List_env1s_size)]
implicit none
integer :: i, j
if(nucl_num .gt. 32) then
print *, ' nucl_num = ', nucl_num, '> 32'
stop
endif
List_env1s = 0
do i = 0, List_env1s_size-1
do j = 0, nucl_num-1
if (btest(i,j)) then
List_env1s(j+1,i+1) = 1
endif
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_env1s_coef, ( List_env1s_size)]
&BEGIN_PROVIDER [ double precision, List_env1s_expo, ( List_env1s_size)]
&BEGIN_PROVIDER [ double precision, List_env1s_cent, (3, List_env1s_size)]
implicit none
integer :: i, j, k, phase
double precision :: tmp_alphaj, tmp_alphak
double precision :: tmp_cent_x, tmp_cent_y, tmp_cent_z
provide env_type env_expo env_coef
List_env1s_coef = 0.d0
List_env1s_expo = 0.d0
List_env1s_cent = 0.d0
if(env_type .eq. "Prod_Gauss") then
do i = 1, List_env1s_size
tmp_cent_x = 0.d0
tmp_cent_y = 0.d0
tmp_cent_z = 0.d0
do j = 1, nucl_num
tmp_alphaj = dble(List_env1s(j,i)) * env_expo(j)
List_env1s_expo(i) += tmp_alphaj
tmp_cent_x += tmp_alphaj * nucl_coord(j,1)
tmp_cent_y += tmp_alphaj * nucl_coord(j,2)
tmp_cent_z += tmp_alphaj * nucl_coord(j,3)
enddo
if(List_env1s_expo(i) .lt. 1d-10) cycle
List_env1s_cent(1,i) = tmp_cent_x / List_env1s_expo(i)
List_env1s_cent(2,i) = tmp_cent_y / List_env1s_expo(i)
List_env1s_cent(3,i) = tmp_cent_z / List_env1s_expo(i)
enddo
! ---
do i = 1, List_env1s_size
do j = 2, nucl_num, 1
tmp_alphaj = dble(List_env1s(j,i)) * env_expo(j)
do k = 1, j-1, 1
tmp_alphak = dble(List_env1s(k,i)) * env_expo(k)
List_env1s_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+ (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+ (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
enddo
enddo
if(List_env1s_expo(i) .lt. 1d-10) cycle
List_env1s_coef(i) = List_env1s_coef(i) / List_env1s_expo(i)
enddo
! ---
do i = 1, List_env1s_size
phase = 0
do j = 1, nucl_num
phase += List_env1s(j,i)
enddo
List_env1s_coef(i) = (-1.d0)**dble(phase) * dexp(-List_env1s_coef(i))
enddo
elseif(env_type .eq. "Sum_Gauss") then
List_env1s_coef( 1) = 1.d0
List_env1s_expo( 1) = 0.d0
List_env1s_cent(1:3,1) = 0.d0
do i = 1, nucl_num
List_env1s_coef( i+1) = -1.d0 * env_coef(i)
List_env1s_expo( i+1) = env_expo(i)
List_env1s_cent(1,i+1) = nucl_coord(i,1)
List_env1s_cent(2,i+1) = nucl_coord(i,2)
List_env1s_cent(3,i+1) = nucl_coord(i,3)
enddo
else
print *, ' Error in List_env1s: Unknown env_type = ', env_type
stop
endif
END_PROVIDER
! ---
BEGIN_PROVIDER [integer, List_env1s_square_size]
implicit none
double precision :: tmp
if(env_type .eq. "Prod_Gauss") then
List_env1s_square_size = 3**nucl_num
elseif(env_type .eq. "Sum_Gauss") then
tmp = 0.5d0 * dble(nucl_num) * (dble(nucl_num) + 3.d0)
List_env1s_square_size = int(tmp) + 1
else
print *, ' Error in List_env1s_square_size: Unknown env_type = ', env_type
stop
endif
print *, ' nb of 1s-Gaussian in the square of envelope = ', List_env1s_square_size
END_PROVIDER
! ---
BEGIN_PROVIDER [integer, List_env1s_square, (nucl_num, List_env1s_square_size)]
implicit none
integer :: i, j, ii, jj
integer, allocatable :: M(:,:), p(:)
if(nucl_num .gt. 32) then
print *, ' nucl_num = ', nucl_num, '> 32'
stop
endif
List_env1s_square(:,:) = 0
List_env1s_square(:,List_env1s_square_size) = 2
allocate(p(nucl_num))
p = 0
do i = 2, List_env1s_square_size-1
do j = 1, nucl_num
ii = 0
do jj = 1, j-1, 1
ii = ii + p(jj) * 3**(jj-1)
enddo
p(j) = modulo(i-1-ii, 3**j) / 3**(j-1)
List_env1s_square(j,i) = p(j)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_env1s_square_coef, ( List_env1s_square_size)]
&BEGIN_PROVIDER [ double precision, List_env1s_square_expo, ( List_env1s_square_size)]
&BEGIN_PROVIDER [ double precision, List_env1s_square_cent, (3, List_env1s_square_size)]
implicit none
integer :: i, j, k, phase
integer :: ii
double precision :: tmp_alphaj, tmp_alphak, facto
double precision :: tmp1, tmp2, tmp3, tmp4
double precision :: xi, yi, zi, xj, yj, zj
double precision :: dx, dy, dz, r2
provide env_type env_expo env_coef
List_env1s_square_coef = 0.d0
List_env1s_square_expo = 0.d0
List_env1s_square_cent = 0.d0
if(env_type .eq. "Prod_Gauss") then
do i = 1, List_env1s_square_size
do j = 1, nucl_num
tmp_alphaj = dble(List_env1s_square(j,i)) * env_expo(j)
List_env1s_square_expo(i) += tmp_alphaj
List_env1s_square_cent(1,i) += tmp_alphaj * nucl_coord(j,1)
List_env1s_square_cent(2,i) += tmp_alphaj * nucl_coord(j,2)
List_env1s_square_cent(3,i) += tmp_alphaj * nucl_coord(j,3)
enddo
if(List_env1s_square_expo(i) .lt. 1d-10) cycle
List_env1s_square_cent(1,i) = List_env1s_square_cent(1,i) / List_env1s_square_expo(i)
List_env1s_square_cent(2,i) = List_env1s_square_cent(2,i) / List_env1s_square_expo(i)
List_env1s_square_cent(3,i) = List_env1s_square_cent(3,i) / List_env1s_square_expo(i)
enddo
! ---
do i = 1, List_env1s_square_size
do j = 2, nucl_num, 1
tmp_alphaj = dble(List_env1s_square(j,i)) * env_expo(j)
do k = 1, j-1, 1
tmp_alphak = dble(List_env1s_square(k,i)) * env_expo(k)
List_env1s_square_coef(i) += tmp_alphaj * tmp_alphak * ( (nucl_coord(j,1) - nucl_coord(k,1)) * (nucl_coord(j,1) - nucl_coord(k,1)) &
+ (nucl_coord(j,2) - nucl_coord(k,2)) * (nucl_coord(j,2) - nucl_coord(k,2)) &
+ (nucl_coord(j,3) - nucl_coord(k,3)) * (nucl_coord(j,3) - nucl_coord(k,3)) )
enddo
enddo
if(List_env1s_square_expo(i) .lt. 1d-10) cycle
List_env1s_square_coef(i) = List_env1s_square_coef(i) / List_env1s_square_expo(i)
enddo
! ---
do i = 1, List_env1s_square_size
facto = 1.d0
phase = 0
do j = 1, nucl_num
tmp_alphaj = dble(List_env1s_square(j,i))
facto *= 2.d0 / (gamma(tmp_alphaj+1.d0) * gamma(3.d0-tmp_alphaj))
phase += List_env1s_square(j,i)
enddo
List_env1s_square_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_env1s_square_coef(i))
enddo
elseif(env_type .eq. "Sum_Gauss") then
ii = 1
List_env1s_square_coef( ii) = 1.d0
List_env1s_square_expo( ii) = 0.d0
List_env1s_square_cent(1:3,ii) = 0.d0
do i = 1, nucl_num
ii = ii + 1
List_env1s_square_coef( ii) = -2.d0 * env_coef(i)
List_env1s_square_expo( ii) = env_expo(i)
List_env1s_square_cent(1,ii) = nucl_coord(i,1)
List_env1s_square_cent(2,ii) = nucl_coord(i,2)
List_env1s_square_cent(3,ii) = nucl_coord(i,3)
enddo
do i = 1, nucl_num
ii = ii + 1
List_env1s_square_coef( ii) = 1.d0 * env_coef(i) * env_coef(i)
List_env1s_square_expo( ii) = 2.d0 * env_expo(i)
List_env1s_square_cent(1,ii) = nucl_coord(i,1)
List_env1s_square_cent(2,ii) = nucl_coord(i,2)
List_env1s_square_cent(3,ii) = nucl_coord(i,3)
enddo
do i = 1, nucl_num-1
tmp1 = env_expo(i)
xi = nucl_coord(i,1)
yi = nucl_coord(i,2)
zi = nucl_coord(i,3)
do j = i+1, nucl_num
tmp2 = env_expo(j)
tmp3 = tmp1 + tmp2
tmp4 = 1.d0 / tmp3
xj = nucl_coord(j,1)
yj = nucl_coord(j,2)
zj = nucl_coord(j,3)
dx = xi - xj
dy = yi - yj
dz = zi - zj
r2 = dx*dx + dy*dy + dz*dz
ii = ii + 1
! x 2 to avoid doing integrals twice
List_env1s_square_coef( ii) = 2.d0 * dexp(-tmp1*tmp2*tmp4*r2) * env_coef(i) * env_coef(j)
List_env1s_square_expo( ii) = tmp3
List_env1s_square_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj)
List_env1s_square_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj)
List_env1s_square_cent(3,ii) = tmp4 * (tmp1 * zi + tmp2 * zj)
enddo
enddo
else
print *, ' Error in List_env1s_square: Unknown env_type = ', env_type
stop
endif
END_PROVIDER
! ---

197
plugins/local/ao_many_one_e_ints/listj1b_sorted.irp.f Normal file
View File

@ -0,0 +1,197 @@
! ---
BEGIN_PROVIDER [integer, List_comb_thr_b2_size, (ao_num, ao_num)]
&BEGIN_PROVIDER [integer, max_List_comb_thr_b2_size]
implicit none
integer :: i_1s, i, j, ipoint
integer :: list(ao_num)
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
List_comb_thr_b2_size = 0
print*,'List_env1s_size = ',List_env1s_size
do i = 1, ao_num
do j = i, ao_num
do i_1s = 1, List_env1s_size
coef = List_env1s_coef(i_1s)
if(dabs(coef).lt.thrsh_cycle_tc) cycle
beta = List_env1s_expo(i_1s)
beta = max(beta,1.d-12)
center(1:3) = List_env1s_cent(1:3,i_1s)
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc)then
List_comb_thr_b2_size(j,i) += 1
endif
enddo
enddo
enddo
do i = 1, ao_num
do j = 1, i-1
List_comb_thr_b2_size(j,i) = List_comb_thr_b2_size(i,j)
enddo
enddo
do i = 1, ao_num
list(i) = maxval(List_comb_thr_b2_size(:,i))
enddo
max_List_comb_thr_b2_size = maxval(list)
print*, ' max_List_comb_thr_b2_size = ',max_List_comb_thr_b2_size
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, List_comb_thr_b2_coef, ( max_List_comb_thr_b2_size,ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_expo, ( max_List_comb_thr_b2_size,ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, List_comb_thr_b2_cent, (3,max_List_comb_thr_b2_size,ao_num,ao_num)]
&BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_env , ( max_List_comb_thr_b2_size,ao_num,ao_num)]
implicit none
integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
ao_abs_comb_b2_env = 10000000.d0
do i = 1, ao_num
do j = i, ao_num
icount = 0
do i_1s = 1, List_env1s_size
coef = List_env1s_coef (i_1s)
if(dabs(coef).lt.thrsh_cycle_tc)cycle
beta = List_env1s_expo (i_1s)
center(1:3) = List_env1s_cent(1:3,i_1s)
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc)then
icount += 1
List_comb_thr_b2_coef(icount,j,i) = coef
List_comb_thr_b2_expo(icount,j,i) = beta
List_comb_thr_b2_cent(1:3,icount,j,i) = center(1:3)
ao_abs_comb_b2_env(icount,j,i) = int_env
endif
enddo
enddo
enddo
do i = 1, ao_num
do j = 1, i-1
do icount = 1, List_comb_thr_b2_size(j,i)
List_comb_thr_b2_coef(icount,j,i) = List_comb_thr_b2_coef(icount,i,j)
List_comb_thr_b2_expo(icount,j,i) = List_comb_thr_b2_expo(icount,i,j)
List_comb_thr_b2_cent(1:3,icount,j,i) = List_comb_thr_b2_cent(1:3,icount,i,j)
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [integer, List_comb_thr_b3_size, (ao_num,ao_num)]
&BEGIN_PROVIDER [integer, max_List_comb_thr_b3_size]
implicit none
integer :: i_1s,i,j,ipoint
integer :: list(ao_num)
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
List_comb_thr_b3_size = 0
print*,'List_env1s_square_size = ',List_env1s_square_size
do i = 1, ao_num
do j = 1, ao_num
do i_1s = 1, List_env1s_square_size
coef = List_env1s_square_coef (i_1s)
beta = List_env1s_square_expo (i_1s)
center(1:3) = List_env1s_square_cent(1:3,i_1s)
if(dabs(coef).lt.thrsh_cycle_tc)cycle
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc) then
List_comb_thr_b3_size(j,i) += 1
endif
enddo
enddo
enddo
do i = 1, ao_num
list(i) = maxval(List_comb_thr_b3_size(:,i))
enddo
max_List_comb_thr_b3_size = maxval(list)
print*, ' max_List_comb_thr_b3_size = ',max_List_comb_thr_b3_size
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, List_comb_thr_b3_coef, ( max_List_comb_thr_b3_size,ao_num,ao_num)]
&BEGIN_PROVIDER [double precision, List_comb_thr_b3_expo, ( max_List_comb_thr_b3_size,ao_num,ao_num)]
&BEGIN_PROVIDER [double precision, List_comb_thr_b3_cent, (3, max_List_comb_thr_b3_size,ao_num,ao_num)]
&BEGIN_PROVIDER [double precision, ao_abs_comb_b3_env , ( max_List_comb_thr_b3_size,ao_num,ao_num)]
implicit none
integer :: i_1s,i,j,ipoint,icount
double precision :: coef,beta,center(3),int_env
double precision :: r(3),weight,dist
ao_abs_comb_b3_env = 10000000.d0
do i = 1, ao_num
do j = 1, ao_num
icount = 0
do i_1s = 1, List_env1s_square_size
coef = List_env1s_square_coef (i_1s)
beta = List_env1s_square_expo (i_1s)
beta = max(beta,1.d-12)
center(1:3) = List_env1s_square_cent(1:3,i_1s)
if(dabs(coef).lt.thrsh_cycle_tc)cycle
int_env = 0.d0
do ipoint = 1, n_points_extra_final_grid
r(1:3) = final_grid_points_extra(1:3,ipoint)
weight = final_weight_at_r_vector_extra(ipoint)
dist = ( center(1) - r(1) )*( center(1) - r(1) )
dist += ( center(2) - r(2) )*( center(2) - r(2) )
dist += ( center(3) - r(3) )*( center(3) - r(3) )
int_env += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
enddo
if(dabs(coef)*dabs(int_env).gt.thrsh_cycle_tc)then
icount += 1
List_comb_thr_b3_coef(icount,j,i) = coef
List_comb_thr_b3_expo(icount,j,i) = beta
List_comb_thr_b3_cent(1:3,icount,j,i) = center(1:3)
ao_abs_comb_b3_env(icount,j,i) = int_env
endif
enddo
enddo
enddo
END_PROVIDER
! ---

0
src/ao_many_one_e_ints/prim_int_erf_gauss.irp.f → plugins/local/ao_many_one_e_ints/prim_int_erf_gauss.irp.f
View File

2
src/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f → plugins/local/ao_many_one_e_ints/prim_int_gauss_gauss.irp.f
View File

@ -200,7 +200,7 @@ subroutine overlap_gauss_r12_v(D_center, LD_D, delta, A_center, B_center, power_
deallocate(A_new, A_center_new, fact_a_new, iorder_a_new, overlap)
end subroutine overlap_gauss_r12_v
end
!---

0
src/ao_many_one_e_ints/stg_gauss_int.irp.f → plugins/local/ao_many_one_e_ints/stg_gauss_int.irp.f
View File

0
src/ao_many_one_e_ints/taylor_exp.irp.f → plugins/local/ao_many_one_e_ints/taylor_exp.irp.f
View File

0
src/ao_many_one_e_ints/xyz_grad_xyz_ao_pol.irp.f → plugins/local/ao_many_one_e_ints/xyz_grad_xyz_ao_pol.irp.f
View File

2
src/ao_tc_eff_map/NEED → plugins/local/ao_tc_eff_map/NEED
View File

@ -3,3 +3,5 @@ mo_one_e_ints
ao_many_one_e_ints
dft_utils_in_r
tc_keywords
hamiltonian
jastrow

0
src/ao_tc_eff_map/README.rst → plugins/local/ao_tc_eff_map/README.rst
View File

11
src/ao_tc_eff_map/compute_ints_eff_pot.irp.f → plugins/local/ao_tc_eff_map/compute_ints_eff_pot.irp.f
View File

@ -23,10 +23,9 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va
logical, external :: ao_two_e_integral_zero
double precision :: ao_tc_sym_two_e_pot, ao_two_e_integral_erf
double precision :: j1b_gauss_2e_j1, j1b_gauss_2e_j2
double precision :: env_gauss_2e_j1, env_gauss_2e_j2
PROVIDE j1b_type
thr = ao_integrals_threshold
@ -53,14 +52,6 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va
integral_erf = ao_two_e_integral_erf(i, k, j, l)
integral = integral_erf + integral_pot
!if( j1b_type .eq. 1 ) then
! !print *, ' j1b type 1 is added'
! integral = integral + j1b_gauss_2e_j1(i, k, j, l)
!elseif( j1b_type .eq. 2 ) then
! !print *, ' j1b type 2 is added'
! integral = integral + j1b_gauss_2e_j2(i, k, j, l)
!endif
if(abs(integral) < thr) then
cycle
endif

0
src/ao_tc_eff_map/integrals_eff_pot_in_map_slave.irp.f → plugins/local/ao_tc_eff_map/integrals_eff_pot_in_map_slave.irp.f
View File

0
src/ao_tc_eff_map/map_integrals_eff_pot.irp.f → plugins/local/ao_tc_eff_map/map_integrals_eff_pot.irp.f
View File

131
src/ao_tc_eff_map/one_e_1bgauss_grad2.irp.f → plugins/local/ao_tc_eff_map/one_e_1bgauss_grad2.irp.f
View File

@ -1,10 +1,10 @@
! ---
BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)]
BEGIN_PROVIDER [double precision, env_gauss_hermII, (ao_num,ao_num)]
BEGIN_DOC
!
! :math:`\langle \chi_A | -0.5 \grad \tau_{1b} \cdot \grad \tau_{1b} | \chi_B \rangle`
! :math:`\langle \chi_A | -0.5 \grad \tau_{env} \cdot \grad \tau_{env} | \chi_B \rangle`
!
END_DOC
@ -22,8 +22,6 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)]
double precision :: int_gauss_4G
PROVIDE j1b_type j1b_pen j1b_coeff
! --------------------------------------------------------------------------------
! -- Dummy call to provide everything
dim1 = 100
@ -38,10 +36,7 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)]
! --------------------------------------------------------------------------------
j1b_gauss_hermII(1:ao_num,1:ao_num) = 0.d0
if(j1b_type .eq. 1) then
! \tau_1b = \sum_iA -[1 - exp(-alpha_A r_iA^2)]
env_gauss_hermII(1:ao_num,1:ao_num) = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -51,113 +46,51 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermII, (ao_num,ao_num)]
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermII)
!$OMP nucl_num, env_expo, env_gauss_hermII)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k1 = 1, nucl_num
gama1 = j1b_pen(k1)
C_center1(1:3) = nucl_coord(k1,1:3)
c = 0.d0
do k1 = 1, nucl_num
gama1 = env_expo(k1)
C_center1(1:3) = nucl_coord(k1,1:3)
do k2 = 1, nucl_num
gama2 = j1b_pen(k2)
C_center2(1:3) = nucl_coord(k2,1:3)
do k2 = 1, nucl_num
gama2 = env_expo(k2)
C_center2(1:3) = nucl_coord(k2,1:3)
! < XA | exp[-gama1 r_C1^2 -gama2 r_C2^2] r_C1 \cdot r_C2 | XB >
c1 = int_gauss_4G( A_center, B_center, C_center1, C_center2 &
, power_A, power_B, alpha, beta, gama1, gama2 )
! < XA | exp[-gama1 r_C1^2 -gama2 r_C2^2] r_C1 \cdot r_C2 | XB >
c1 = int_gauss_4G( A_center, B_center, C_center1, C_center2 &
, power_A, power_B, alpha, beta, gama1, gama2 )
c = c - 2.d0 * gama1 * gama2 * c1
enddo
c = c - 2.d0 * gama1 * gama2 * c1
enddo
j1b_gauss_hermII(i,j) = j1b_gauss_hermII(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
env_gauss_hermII(i,j) = env_gauss_hermII(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
elseif(j1b_type .eq. 2) then
! \tau_1b = \sum_iA [c_A exp(-alpha_A r_iA^2)]
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k1, k2, l, m, alpha, beta, gama1, gama2, &
!$OMP A_center, B_center, C_center1, C_center2, &
!$OMP power_A, power_B, num_A, num_B, c1, c, &
!$OMP coef1, coef2) &
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermII, &
!$OMP j1b_coeff)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k1 = 1, nucl_num
gama1 = j1b_pen (k1)
coef1 = j1b_coeff(k1)
C_center1(1:3) = nucl_coord(k1,1:3)
do k2 = 1, nucl_num
gama2 = j1b_pen (k2)
coef2 = j1b_coeff(k2)
C_center2(1:3) = nucl_coord(k2,1:3)
! < XA | exp[-gama1 r_C1^2 -gama2 r_C2^2] r_C1 \cdot r_C2 | XB >
c1 = int_gauss_4G( A_center, B_center, C_center1, C_center2 &
, power_A, power_B, alpha, beta, gama1, gama2 )
c = c - 2.d0 * gama1 * gama2 * coef1 * coef2 * c1
enddo
enddo
j1b_gauss_hermII(i,j) = j1b_gauss_hermII(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
END_PROVIDER

128
src/ao_tc_eff_map/one_e_1bgauss_lap.irp.f → plugins/local/ao_tc_eff_map/one_e_1bgauss_lap.irp.f
View File

@ -1,10 +1,10 @@
! ---
BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)]
BEGIN_PROVIDER [double precision, env_gauss_hermI, (ao_num,ao_num)]
BEGIN_DOC
!
! :math:`\langle \chi_A | -0.5 \Delta \tau_{1b} | \chi_B \rangle`
! :math:`\langle \chi_A | -0.5 \Delta \tau_{env} | \chi_B \rangle`
!
END_DOC
@ -22,8 +22,6 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)]
double precision :: int_gauss_r0, int_gauss_r2
PROVIDE j1b_type j1b_pen j1b_coeff
! --------------------------------------------------------------------------------
! -- Dummy call to provide everything
dim1 = 100
@ -37,10 +35,7 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)]
, overlap_y, d_a_2, overlap_z, overlap, dim1 )
! --------------------------------------------------------------------------------
j1b_gauss_hermI(1:ao_num,1:ao_num) = 0.d0
if(j1b_type .eq. 1) then
! \tau_1b = \sum_iA -[1 - exp(-alpha_A r_iA^2)]
env_gauss_hermI(1:ao_num,1:ao_num) = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -50,109 +45,50 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_hermI, (ao_num,ao_num)]
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermI)
!$OMP nucl_num, env_expo, env_gauss_hermI)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
gama = j1b_pen(k)
C_center(1:3) = nucl_coord(k,1:3)
c = 0.d0
do k = 1, nucl_num
gama = env_expo(k)
C_center(1:3) = nucl_coord(k,1:3)
! < XA | exp[-gama r_C^2] | XB >
c1 = int_gauss_r0( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
! < XA | exp[-gama r_C^2] | XB >
c1 = int_gauss_r0( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
! < XA | r_A^2 exp[-gama r_C^2] | XB >
c2 = int_gauss_r2( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
! < XA | r_A^2 exp[-gama r_C^2] | XB >
c2 = int_gauss_r2( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
c = c + 3.d0 * gama * c1 - 2.d0 * gama * gama * c2
enddo
j1b_gauss_hermI(i,j) = j1b_gauss_hermI(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
c = c + 3.d0 * gama * c1 - 2.d0 * gama * gama * c2
enddo
env_gauss_hermI(i,j) = env_gauss_hermI(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
elseif(j1b_type .eq. 2) then
! \tau_1b = \sum_iA [c_A exp(-alpha_A r_iA^2)]
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, m, alpha, beta, gama, coef, &
!$OMP A_center, B_center, C_center, power_A, power_B, &
!$OMP num_A, num_B, c1, c2, c) &
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_hermI, &
!$OMP j1b_coeff)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
gama = j1b_pen (k)
coef = j1b_coeff(k)
C_center(1:3) = nucl_coord(k,1:3)
! < XA | exp[-gama r_C^2] | XB >
c1 = int_gauss_r0( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
! < XA | r_A^2 exp[-gama r_C^2] | XB >
c2 = int_gauss_r2( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
c = c + 3.d0 * gama * coef * c1 - 2.d0 * gama * gama * coef * c2
enddo
j1b_gauss_hermI(i,j) = j1b_gauss_hermI(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
END_PROVIDER

118
src/ao_tc_eff_map/one_e_1bgauss_nonherm.irp.f → plugins/local/ao_tc_eff_map/one_e_1bgauss_nonherm.irp.f
View File

@ -1,10 +1,11 @@
! ---
BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)]
BEGIN_PROVIDER [double precision, env_gauss_nonherm, (ao_num,ao_num)]
BEGIN_DOC
!
! j1b_gauss_nonherm(i,j) = \langle \chi_j | - grad \tau_{1b} \cdot grad | \chi_i \rangle
! env_gauss_nonherm(i,j) = \langle \chi_j | - grad \tau_{env} \cdot grad | \chi_i \rangle
!
END_DOC
@ -22,8 +23,6 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)]
double precision :: int_gauss_deriv
PROVIDE j1b_type j1b_pen j1b_coeff
! --------------------------------------------------------------------------------
! -- Dummy call to provide everything
dim1 = 100
@ -38,10 +37,8 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)]
! --------------------------------------------------------------------------------
j1b_gauss_nonherm(1:ao_num,1:ao_num) = 0.d0
env_gauss_nonherm(1:ao_num,1:ao_num) = 0.d0
if(j1b_type .eq. 1) then
! \tau_1b = \sum_iA -[1 - exp(-alpha_A r_iA^2)]
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -51,101 +48,46 @@ BEGIN_PROVIDER [ double precision, j1b_gauss_nonherm, (ao_num,ao_num)]
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_nonherm)
!$OMP nucl_num, env_expo, env_gauss_nonherm)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
gama = j1b_pen(k)
C_center(1:3) = nucl_coord(k,1:3)
c = 0.d0
do k = 1, nucl_num
gama = env_expo(k)
C_center(1:3) = nucl_coord(k,1:3)
! \langle \chi_A | exp[-gama r_C^2] r_C \cdot grad | \chi_B \rangle
c1 = int_gauss_deriv( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
! \langle \chi_A | exp[-gama r_C^2] r_C \cdot grad | \chi_B \rangle
c1 = int_gauss_deriv( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
c = c + 2.d0 * gama * c1
enddo
j1b_gauss_nonherm(i,j) = j1b_gauss_nonherm(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
c = c + 2.d0 * gama * c1
enddo
env_gauss_nonherm(i,j) = env_gauss_nonherm(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
elseif(j1b_type .eq. 2) then
! \tau_1b = \sum_iA [c_A exp(-alpha_A r_iA^2)]
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, m, alpha, beta, gama, coef, &
!$OMP A_center, B_center, C_center, power_A, power_B, &
!$OMP num_A, num_B, c1, c) &
!$OMP SHARED (ao_num, ao_prim_num, ao_expo_ordered_transp, &
!$OMP ao_power, ao_nucl, nucl_coord, &
!$OMP ao_coef_normalized_ordered_transp, &
!$OMP nucl_num, j1b_pen, j1b_gauss_nonherm, &
!$OMP j1b_coeff)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ordered_transp(l,j)
do m = 1, ao_prim_num(i)
beta = ao_expo_ordered_transp(m,i)
c = 0.d0
do k = 1, nucl_num
gama = j1b_pen (k)
coef = j1b_coeff(k)
C_center(1:3) = nucl_coord(k,1:3)
! \langle \chi_A | exp[-gama r_C^2] r_C \cdot grad | \chi_B \rangle
c1 = int_gauss_deriv( A_center, B_center, C_center &
, power_A, power_B, alpha, beta, gama )
c = c + 2.d0 * gama * coef * c1
enddo
j1b_gauss_nonherm(i,j) = j1b_gauss_nonherm(i,j) &
+ ao_coef_normalized_ordered_transp(l,j) &
* ao_coef_normalized_ordered_transp(m,i) * c
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
endif
END_PROVIDER

3
src/ao_tc_eff_map/providers_ao_eff_pot.irp.f → plugins/local/ao_tc_eff_map/providers_ao_eff_pot.irp.f
View File

@ -22,9 +22,6 @@ BEGIN_PROVIDER [ logical, ao_tc_sym_two_e_pot_in_map ]
integer :: kk, m, j1, i1, lmax
character*(64) :: fmt
!double precision :: j1b_gauss_coul_debug
!integral = j1b_gauss_coul_debug(1,1,1,1)
integral = ao_tc_sym_two_e_pot(1,1,1,1)
double precision :: map_mb

27
src/ao_tc_eff_map/two_e_1bgauss_j1.irp.f → plugins/local/ao_tc_eff_map/two_e_1bgauss_j1.irp.f
View File

@ -1,6 +1,6 @@
! ---
double precision function j1b_gauss_2e_j1(i, j, k, l)
double precision function env_gauss_2e_j1(i, j, k, l)
BEGIN_DOC
!
@ -36,10 +36,10 @@ double precision function j1b_gauss_2e_j1(i, j, k, l)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: ff, gg, cx, cy, cz
double precision :: j1b_gauss_2e_j1_schwartz
double precision :: env_gauss_2e_j1_schwartz
if( ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
j1b_gauss_2e_j1 = j1b_gauss_2e_j1_schwartz(i, j, k, l)
env_gauss_2e_j1 = env_gauss_2e_j1_schwartz(i, j, k, l)
return
endif
@ -59,7 +59,7 @@ double precision function j1b_gauss_2e_j1(i, j, k, l)
L_center(p) = nucl_coord(num_l,p)
enddo
j1b_gauss_2e_j1 = 0.d0
env_gauss_2e_j1 = 0.d0
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p, i)
@ -89,18 +89,18 @@ double precision function j1b_gauss_2e_j1(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j1 = j1b_gauss_2e_j1 + coef4 * ( cx + cy + cz )
env_gauss_2e_j1 = env_gauss_2e_j1 + coef4 * ( cx + cy + cz )
enddo ! s
enddo ! r
enddo ! q
enddo ! p
return
end function j1b_gauss_2e_j1
end
! ---
double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
double precision function env_gauss_2e_j1_schwartz(i, j, k, l)
BEGIN_DOC
!
@ -137,8 +137,6 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
double precision :: schwartz_ij, thr
double precision, allocatable :: schwartz_kl(:,:)
PROVIDE j1b_pen
dim1 = n_pt_max_integrals
thr = ao_integrals_threshold * ao_integrals_threshold
@ -186,8 +184,7 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
schwartz_kl(0,0) = max( schwartz_kl(0,r) , schwartz_kl(0,0) )
enddo
j1b_gauss_2e_j1_schwartz = 0.d0
env_gauss_2e_j1_schwartz = 0.d0
do p = 1, ao_prim_num(i)
expo1 = ao_expo_ordered_transp(p, i)
@ -226,7 +223,7 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j1_schwartz = j1b_gauss_2e_j1_schwartz + coef4 * ( cx + cy + cz )
env_gauss_2e_j1_schwartz = env_gauss_2e_j1_schwartz + coef4 * ( cx + cy + cz )
enddo ! s
enddo ! r
enddo ! q
@ -235,7 +232,7 @@ double precision function j1b_gauss_2e_j1_schwartz(i, j, k, l)
deallocate( schwartz_kl )
return
end function j1b_gauss_2e_j1_schwartz
end
! ---
@ -263,14 +260,12 @@ subroutine get_cxcycz_j1( dim1, cx, cy, cz &
double precision :: general_primitive_integral_erf_shifted
double precision :: general_primitive_integral_coul_shifted
PROVIDE j1b_pen
cx = 0.d0
cy = 0.d0
cz = 0.d0
do ii = 1, nucl_num
expoii = j1b_pen(ii)
expoii = env_expo(ii)
Centerii(1:3) = nucl_coord(ii, 1:3)
call gaussian_product(pp1, P1_center, expoii, Centerii, factii, pp2, P2_center)

26
src/ao_tc_eff_map/two_e_1bgauss_j2.irp.f → plugins/local/ao_tc_eff_map/two_e_1bgauss_j2.irp.f
View File

@ -1,6 +1,6 @@
! ---
double precision function j1b_gauss_2e_j2(i, j, k, l)
double precision function env_gauss_2e_j2(i, j, k, l)
BEGIN_DOC
!
@ -36,12 +36,12 @@ double precision function j1b_gauss_2e_j2(i, j, k, l)
double precision :: I_center(3), J_center(3), K_center(3), L_center(3)
double precision :: ff, gg, cx, cy, cz
double precision :: j1b_gauss_2e_j2_schwartz
double precision :: env_gauss_2e_j2_schwartz
dim1 = n_pt_max_integrals
if( ao_prim_num(i) * ao_prim_num(j) * ao_prim_num(k) * ao_prim_num(l) > 1024 ) then
j1b_gauss_2e_j2 = j1b_gauss_2e_j2_schwartz(i, j, k, l)
env_gauss_2e_j2 = env_gauss_2e_j2_schwartz(i, j, k, l)
return
endif
@ -61,7 +61,7 @@ double precision function j1b_gauss_2e_j2(i, j, k, l)
L_center(p) = nucl_coord(num_l,p)
enddo
j1b_gauss_2e_j2 = 0.d0
env_gauss_2e_j2 = 0.d0
do p = 1, ao_prim_num(i)
coef1 = ao_coef_normalized_ordered_transp(p, i)
@ -91,18 +91,18 @@ double precision function j1b_gauss_2e_j2(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j2 = j1b_gauss_2e_j2 + coef4 * ( cx + cy + cz )
env_gauss_2e_j2 = env_gauss_2e_j2 + coef4 * ( cx + cy + cz )
enddo ! s
enddo ! r
enddo ! q
enddo ! p
return
end function j1b_gauss_2e_j2
end
! ---
double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l)
double precision function env_gauss_2e_j2_schwartz(i, j, k, l)
BEGIN_DOC
!
@ -187,7 +187,7 @@ double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l)
enddo
j1b_gauss_2e_j2_schwartz = 0.d0
env_gauss_2e_j2_schwartz = 0.d0
do p = 1, ao_prim_num(i)
expo1 = ao_expo_ordered_transp(p, i)
@ -226,7 +226,7 @@ double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l)
, P1_center, P1_new, pp1, fact_p1, p1_inv, iorder_p &
, Q1_center, Q1_new, qq1, fact_q1, q1_inv, iorder_q )
j1b_gauss_2e_j2_schwartz = j1b_gauss_2e_j2_schwartz + coef4 * ( cx + cy + cz )
env_gauss_2e_j2_schwartz = env_gauss_2e_j2_schwartz + coef4 * ( cx + cy + cz )
enddo ! s
enddo ! r
enddo ! q
@ -235,7 +235,7 @@ double precision function j1b_gauss_2e_j2_schwartz(i, j, k, l)
deallocate( schwartz_kl )
return
end function j1b_gauss_2e_j2_schwartz
end
! ---
@ -263,15 +263,13 @@ subroutine get_cxcycz_j2( dim1, cx, cy, cz &
double precision :: general_primitive_integral_erf_shifted
double precision :: general_primitive_integral_coul_shifted
PROVIDE j1b_pen j1b_coeff
cx = 0.d0
cy = 0.d0
cz = 0.d0
do ii = 1, nucl_num
expoii = j1b_pen (ii)
coefii = j1b_coeff(ii)
expoii = env_expo(ii)
coefii = env_coef(ii)
Centerii(1:3) = nucl_coord(ii, 1:3)
call gaussian_product(pp1, P1_center, expoii, Centerii, factii, pp2, P2_center)

0
src/ao_tc_eff_map/two_e_ints_gauss.irp.f → plugins/local/ao_tc_eff_map/two_e_ints_gauss.irp.f
View File

49
src/ao_tc_eff_map/useful_sub.irp.f → plugins/local/ao_tc_eff_map/useful_sub.irp.f
View File

@ -174,7 +174,7 @@ double precision function general_primitive_integral_coul_shifted( dim
general_primitive_integral_coul_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q)
return
end function general_primitive_integral_coul_shifted
end
!______________________________________________________________________________________________________________________
!______________________________________________________________________________________________________________________
@ -354,7 +354,7 @@ double precision function general_primitive_integral_erf_shifted( dim
general_primitive_integral_erf_shifted = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / dsqrt(p_plus_q)
return
end function general_primitive_integral_erf_shifted
end
!______________________________________________________________________________________________________________________
!______________________________________________________________________________________________________________________
@ -362,3 +362,48 @@ end function general_primitive_integral_erf_shifted
! ---
subroutine inv_r_times_poly(r, dist_r, dist_vec, poly)
BEGIN_DOC
!
! returns
!
! poly(1) = x / sqrt(x^2+y^2+z^2), poly(2) = y / sqrt(x^2+y^2+z^2), poly(3) = z / sqrt(x^2+y^2+z^2)
!
! with the arguments
!
! r(1) = x, r(2) = y, r(3) = z, dist_r = sqrt(x^2+y^2+z^2)
!
! dist_vec(1) = sqrt(y^2+z^2), dist_vec(2) = sqrt(x^2+z^2), dist_vec(3) = sqrt(x^2+y^2)
!
END_DOC
implicit none
double precision, intent(in) :: r(3), dist_r, dist_vec(3)
double precision, intent(out) :: poly(3)
integer :: i
double precision :: inv_dist
if (dist_r .gt. 1.d-8)then
inv_dist = 1.d0/dist_r
do i = 1, 3
poly(i) = r(i) * inv_dist
enddo
else
do i = 1, 3
if(dabs(r(i)).lt.dist_vec(i)) then
inv_dist = 1.d0/dist_r
poly(i) = r(i) * inv_dist
else
poly(i) = 1.d0
endif
enddo
endif
end
! ---

0
src/basis_correction/51.basis_c.bats → plugins/local/basis_correction/51.basis_c.bats
View File

0
src/basis_correction/NEED → plugins/local/basis_correction/NEED
View File

2
src/basis_correction/README.rst → plugins/local/basis_correction/README.rst
View File

@ -12,7 +12,7 @@ This basis set correction relies mainy on :
When HF is a qualitative representation of the electron pairs (i.e. weakly correlated systems), such an approach for \mu(r) is OK.
See for instance JPCL, 10, 2931-2937 (2019) for typical flavours of the results.
Thanks to the trivial nature of such a two-body rdm, the equation (22) of J. Chem. Phys. 149, 194301 (2018) can be rewritten in a very efficient way, and therefore the limiting factor of such an approach is the AO->MO four-index transformation of the two-electron integrals.
b) "mu_of_r_potential = cas_ful" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder).
b) "mu_of_r_potential = cas_full" uses the two-body rdm of CAS-like wave function (i.e. linear combination of Slater determinants developped in an active space with the MOs stored in the EZFIO folder).
If the CAS is properly chosen (i.e. the CAS-like wave function qualitatively represents the wave function of the systems), then such an approach is OK for \mu(r) even in the case of strong correlation.
+) The use of DFT correlation functionals with multi-determinant reference (Ecmd). These functionals are originally defined in the RS-DFT framework (see for instance Theor. Chem. Acc.114, 305(2005)) and design to capture short-range correlation effects. A important quantity arising in the Ecmd is the exact on-top pair density of the system, and the main differences of approximated Ecmd relies on different approximations for the exact on-top pair density.

0
src/basis_correction/TODO → plugins/local/basis_correction/TODO
View File

0
src/basis_correction/basis_correction.irp.f → plugins/local/basis_correction/basis_correction.irp.f
View File

0
src/basis_correction/eff_xi_based_func.irp.f → plugins/local/basis_correction/eff_xi_based_func.irp.f
View File

6
src/basis_correction/pbe_on_top.irp.f → plugins/local/basis_correction/pbe_on_top.irp.f
View File

@ -39,7 +39,7 @@
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
if(mu_of_r_potential == "cas_full")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else
@ -101,7 +101,7 @@
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
if(mu_of_r_potential == "cas_full")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else
@ -163,7 +163,7 @@
grad_rho_a(1:3) = one_e_dm_and_grad_alpha_in_r(1:3,ipoint,istate)
grad_rho_b(1:3) = one_e_dm_and_grad_beta_in_r(1:3,ipoint,istate)
if(mu_of_r_potential == "cas_ful")then
if(mu_of_r_potential == "cas_full")then
! You take the on-top of the CAS wave function which is computed with mu(r)
on_top = on_top_cas_mu_r(ipoint,istate)
else

4
src/basis_correction/print_routine.irp.f → plugins/local/basis_correction/print_routine.irp.f
View File

@ -4,7 +4,7 @@ subroutine print_basis_correction
provide mu_average_prov
if(mu_of_r_potential.EQ."hf")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated")then
else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated")then
provide ecmd_lda_mu_of_r ecmd_pbe_ueg_mu_of_r
provide ecmd_pbe_on_top_mu_of_r ecmd_pbe_on_top_su_mu_of_r
endif
@ -38,7 +38,7 @@ subroutine print_basis_correction
write(*, '(A29,X,I3,X,A3,X,F16.10)') ' ECMD PBE-UEG , state ',istate,' = ',ecmd_pbe_ueg_mu_of_r(istate)
enddo
else if(mu_of_r_potential.EQ."cas_ful".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
else if(mu_of_r_potential.EQ."cas_full".or.mu_of_r_potential.EQ."cas_truncated".or.mu_of_r_potential.EQ."pure_act")then
print*, ''
print*,'Using a CAS-like two-body density to define mu(r)'
print*,'This assumes that the CAS is a qualitative representation of the wave function '

0
src/basis_correction/print_su_pbe_ot.irp.f → plugins/local/basis_correction/print_su_pbe_ot.irp.f
View File

0
src/basis_correction/weak_corr_func.irp.f → plugins/local/basis_correction/weak_corr_func.irp.f
View File

0
src/bi_ort_ints/NEED → plugins/local/bi_ort_ints/NEED
View File

0
src/bi_ort_ints/README.rst → plugins/local/bi_ort_ints/README.rst
View File

0
src/bi_ort_ints/bi_ort_ints.irp.f → plugins/local/bi_ort_ints/bi_ort_ints.irp.f
View File

37
src/bi_ort_ints/biorthog_mo_for_h.irp.f → plugins/local/bi_ort_ints/biorthog_mo_for_h.irp.f
View File

@ -1,4 +1,39 @@
! ---
BEGIN_PROVIDER [double precision, ao_two_e_coul, (ao_num, ao_num, ao_num, ao_num) ]
BEGIN_DOC
!
! ao_two_e_coul(k,i,l,j) = ( k i | 1/r12 | l j ) = < l k | 1/r12 | j i >
!
END_DOC
integer :: i, j, k, l
double precision, external :: get_ao_two_e_integral
PROVIDE ao_integrals_map
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(ao_num, ao_two_e_coul, ao_integrals_map) &
!$OMP PRIVATE(i, j, k, l)
!$OMP DO
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! < 1:k, 2:l | 1:i, 2:j >
ao_two_e_coul(k,i,l,j) = get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
END_PROVIDER
! ---
double precision function bi_ortho_mo_coul_ints(l, k, j, i)
@ -25,7 +60,7 @@ double precision function bi_ortho_mo_coul_ints(l, k, j, i)
enddo
enddo
end function bi_ortho_mo_coul_ints
end
! ---

0
src/bi_ort_ints/no_dressing.irp.f → plugins/local/bi_ort_ints/no_dressing.irp.f
View File

0
src/bi_ort_ints/no_dressing_energy.irp.f → plugins/local/bi_ort_ints/no_dressing_energy.irp.f
View File

0
src/bi_ort_ints/no_dressing_naive.irp.f → plugins/local/bi_ort_ints/no_dressing_naive.irp.f
View File

17
src/bi_ort_ints/one_e_bi_ort.irp.f → plugins/local/bi_ort_ints/one_e_bi_ort.irp.f
View File

@ -8,23 +8,6 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
ao_one_e_integrals_tc_tot = ao_one_e_integrals
!provide j1b_type
!if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
!
! print *, ' do things properly !'
! stop
! !do i = 1, ao_num
! ! do j = 1, ao_num
! ! ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI (j,i) &
! ! + j1b_gauss_hermII (j,i) &
! ! + j1b_gauss_nonherm(j,i) )
! ! enddo
! !enddo
!endif
END_PROVIDER
! ---

0
src/bi_ort_ints/semi_num_ints_mo.irp.f → plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f
View File

0
src/bi_ort_ints/three_body_ijm.irp.f → plugins/local/bi_ort_ints/three_body_ijm.irp.f
View File

0
src/bi_ort_ints/three_body_ijmk.irp.f → plugins/local/bi_ort_ints/three_body_ijmk.irp.f
View File

0
src/bi_ort_ints/three_body_ijmk_n4.irp.f → plugins/local/bi_ort_ints/three_body_ijmk_n4.irp.f
View File

0
src/bi_ort_ints/three_body_ijmk_old.irp.f → plugins/local/bi_ort_ints/three_body_ijmk_old.irp.f
View File

0
src/bi_ort_ints/three_body_ijmkl.irp.f → plugins/local/bi_ort_ints/three_body_ijmkl.irp.f
View File

0
src/bi_ort_ints/three_body_ijmkl_old.irp.f → plugins/local/bi_ort_ints/three_body_ijmkl_old.irp.f
View File

0
src/bi_ort_ints/three_body_ints_bi_ort.irp.f → plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f
View File

90
src/bi_ort_ints/total_twoe_pot.irp.f → plugins/local/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -1,91 +1,4 @@
! ---
BEGIN_PROVIDER [double precision, ao_two_e_vartc_tot, (ao_num, ao_num, ao_num, ao_num) ]
integer :: i, j, k, l
provide j1b_type
provide mo_r_coef mo_l_coef
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_two_e_vartc_tot(k,i,l,j) = ao_vartc_int_chemist(k,i,l,j)
enddo
enddo
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_num) ]
BEGIN_DOC
!
! ao_two_e_tc_tot(k,i,l,j) = (ki|V^TC(r_12)|lj) = <lk| V^TC(r_12) |ji> where V^TC(r_12) is the total TC operator
!
! including both hermitian and non hermitian parts. THIS IS IN CHEMIST NOTATION.
!
! WARNING :: non hermitian ! acts on "the right functions" (i,j)
!
END_DOC
integer :: i, j, k, l
double precision :: integral_sym, integral_nsym
double precision, external :: get_ao_tc_sym_two_e_pot
provide j1b_type
if(j1b_type .eq. 0) then
PROVIDE ao_tc_sym_two_e_pot_in_map
!!! TODO :: OPENMP
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
integral_sym = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map)
! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
integral_nsym = ao_non_hermit_term_chemist(k,i,l,j)
!print *, ' sym integ = ', integral_sym
!print *, ' non-sym integ = ', integral_nsym
ao_two_e_tc_tot(k,i,l,j) = integral_sym + integral_nsym
!write(111,*) ao_two_e_tc_tot(k,i,l,j)
enddo
enddo
enddo
enddo
else
PROVIDE ao_tc_int_chemist
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
ao_two_e_tc_tot(k,i,l,j) = ao_tc_int_chemist(k,i,l,j)
!write(222,*) ao_two_e_tc_tot(k,i,l,j)
enddo
enddo
enddo
enddo
FREE ao_tc_int_chemist
endif
END_PROVIDER
! ---
double precision function bi_ortho_mo_ints(l, k, j, i)
@ -118,8 +31,6 @@ end function bi_ortho_mo_ints
! ---
! TODO :: transform into DEGEMM
BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num, mo_num, mo_num)]
BEGIN_DOC
@ -267,7 +178,6 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num)]
&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num)]
&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num)]

0
src/bi_ortho_mos/EZFIO.cfg → plugins/local/bi_ortho_mos/EZFIO.cfg
View File

0
src/bi_ortho_mos/NEED → plugins/local/bi_ortho_mos/NEED
View File

0
src/bi_ortho_mos/bi_density.irp.f → plugins/local/bi_ortho_mos/bi_density.irp.f
View File

0
src/bi_ortho_mos/bi_ort_mos_in_r.irp.f → plugins/local/bi_ortho_mos/bi_ort_mos_in_r.irp.f
View File

0
src/bi_ortho_mos/grad_bi_ort_mos_in_r.irp.f → plugins/local/bi_ortho_mos/grad_bi_ort_mos_in_r.irp.f
View File

0
src/bi_ortho_mos/mos_rl.irp.f → plugins/local/bi_ortho_mos/mos_rl.irp.f
View File

0
src/bi_ortho_mos/overlap.irp.f → plugins/local/bi_ortho_mos/overlap.irp.f
View File

0
src/casscf_tc_bi/NEED → plugins/local/casscf_tc_bi/NEED
View File

0
src/casscf_tc_bi/det_manip.irp.f → plugins/local/casscf_tc_bi/det_manip.irp.f
View File

0
src/casscf_tc_bi/grad_dm.irp.f → plugins/local/casscf_tc_bi/grad_dm.irp.f
View File

0
src/casscf_tc_bi/grad_old.irp.f → plugins/local/casscf_tc_bi/grad_old.irp.f
View File

0
src/casscf_tc_bi/gradient.irp.f → plugins/local/casscf_tc_bi/gradient.irp.f
View File

0
src/casscf_tc_bi/test_tc_casscf.irp.f → plugins/local/casscf_tc_bi/test_tc_casscf.irp.f
View File

0
src/cipsi_tc_bi_ortho/EZFIO.cfg → plugins/local/cipsi_tc_bi_ortho/EZFIO.cfg
View File

0
src/cipsi_tc_bi_ortho/NEED → plugins/local/cipsi_tc_bi_ortho/NEED
View File

0
src/cipsi_tc_bi_ortho/cipsi.irp.f → plugins/local/cipsi_tc_bi_ortho/cipsi.irp.f
View File

0
src/cipsi_tc_bi_ortho/energy.irp.f → plugins/local/cipsi_tc_bi_ortho/energy.irp.f
View File

0
src/cipsi_tc_bi_ortho/environment.irp.f → plugins/local/cipsi_tc_bi_ortho/environment.irp.f
View File

0
src/cipsi_tc_bi_ortho/fock_diag.irp.f → plugins/local/cipsi_tc_bi_ortho/fock_diag.irp.f
View File

0
src/cipsi_tc_bi_ortho/get_d.irp.f → plugins/local/cipsi_tc_bi_ortho/get_d.irp.f
View File

0
src/cipsi_tc_bi_ortho/get_d0_good.irp.f → plugins/local/cipsi_tc_bi_ortho/get_d0_good.irp.f
View File

0
src/cipsi_tc_bi_ortho/get_d1_good.irp.f → plugins/local/cipsi_tc_bi_ortho/get_d1_good.irp.f
View File

0
src/cipsi_tc_bi_ortho/get_d2_good.irp.f → plugins/local/cipsi_tc_bi_ortho/get_d2_good.irp.f
View File

0
src/cipsi_tc_bi_ortho/lock_2rdm.irp.f → plugins/local/cipsi_tc_bi_ortho/lock_2rdm.irp.f
View File

0
src/cipsi_tc_bi_ortho/pouet → plugins/local/cipsi_tc_bi_ortho/pouet
View File

0
src/cipsi_tc_bi_ortho/pt2.irp.f → plugins/local/cipsi_tc_bi_ortho/pt2.irp.f
View File

0
src/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f → plugins/local/cipsi_tc_bi_ortho/pt2_stoch_routines.irp.f
View File

0
src/cipsi_tc_bi_ortho/pt2_type.irp.f → plugins/local/cipsi_tc_bi_ortho/pt2_type.irp.f
View File

0
src/cipsi_tc_bi_ortho/run_pt2_slave.irp.f → plugins/local/cipsi_tc_bi_ortho/run_pt2_slave.irp.f
View File

0
src/cipsi_tc_bi_ortho/run_selection_slave.irp.f → plugins/local/cipsi_tc_bi_ortho/run_selection_slave.irp.f
View File

0
src/cipsi_tc_bi_ortho/selection.irp.f → plugins/local/cipsi_tc_bi_ortho/selection.irp.f
View File

0
src/cipsi_tc_bi_ortho/selection_buffer.irp.f → plugins/local/cipsi_tc_bi_ortho/selection_buffer.irp.f
View File

0
src/cipsi_tc_bi_ortho/selection_types.f90 → plugins/local/cipsi_tc_bi_ortho/selection_types.f90
View File

0
src/cipsi_tc_bi_ortho/selection_weight.irp.f → plugins/local/cipsi_tc_bi_ortho/selection_weight.irp.f
View File

0
src/cipsi_tc_bi_ortho/slave_cipsi.irp.f → plugins/local/cipsi_tc_bi_ortho/slave_cipsi.irp.f
View File

Some files were not shown because too many files have changed in this diff Show More