From b4b6ef4266bbb5f40685395a29937201b09148c7 Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 12 Apr 2025 16:28:42 +0200 Subject: [PATCH] added ~/qp2/src/ao_cart_one_e_ints --- src/ao_cart_one_e_ints/EZFIO.cfg | 92 + src/ao_cart_one_e_ints/README.rst | 13 + src/ao_cart_one_e_ints/ao_one_e_ints.irp.f | 47 + src/ao_cart_one_e_ints/ao_overlap.irp.f | 233 ++ src/ao_cart_one_e_ints/aos_cgtos.irp.f | 124 + src/ao_cart_one_e_ints/kin_ao_ints.irp.f | 191 ++ .../one_e_coul_integrals_cgtos.irp.f | 558 +++++ .../one_e_kin_integrals_cgtos.irp.f | 318 +++ src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f | 794 ++++++ src/ao_cart_one_e_ints/pot_ao_ints.irp.f | 610 +++++ .../pot_ao_pseudo_ints.irp.f | 311 +++ src/ao_cart_one_e_ints/pot_pt_charges.irp.f | 108 + src/ao_cart_one_e_ints/pseudopot.f90 | 2127 +++++++++++++++++ src/ao_cart_one_e_ints/spread_dipole_ao.irp.f | 376 +++ 14 files changed, 5902 insertions(+) create mode 100644 src/ao_cart_one_e_ints/EZFIO.cfg create mode 100644 src/ao_cart_one_e_ints/README.rst create mode 100644 src/ao_cart_one_e_ints/ao_one_e_ints.irp.f create mode 100644 src/ao_cart_one_e_ints/ao_overlap.irp.f create mode 100644 src/ao_cart_one_e_ints/aos_cgtos.irp.f create mode 100644 src/ao_cart_one_e_ints/kin_ao_ints.irp.f create mode 100644 src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f create mode 100644 src/ao_cart_one_e_ints/one_e_kin_integrals_cgtos.irp.f create mode 100644 src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f create mode 100644 src/ao_cart_one_e_ints/pot_ao_ints.irp.f create mode 100644 src/ao_cart_one_e_ints/pot_ao_pseudo_ints.irp.f create mode 100644 src/ao_cart_one_e_ints/pot_pt_charges.irp.f create mode 100644 src/ao_cart_one_e_ints/pseudopot.f90 create mode 100644 src/ao_cart_one_e_ints/spread_dipole_ao.irp.f diff --git a/src/ao_cart_one_e_ints/EZFIO.cfg b/src/ao_cart_one_e_ints/EZFIO.cfg new file mode 100644 index 00000000..620afd80 --- /dev/null +++ b/src/ao_cart_one_e_ints/EZFIO.cfg @@ -0,0 +1,92 @@ +[ao_cart_integrals_n_e] +type: double precision +doc: Nucleus-electron integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[ao_cart_integrals_n_e_imag] +type: double precision +doc: Imaginary part of the nucleus-electron integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[io_ao_cart_integrals_n_e] +type: Disk_access +doc: Read/Write |AO| nucleus-electron attraction integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + +[ao_cart_integrals_kinetic] +type: double precision +doc: Kinetic energy integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[ao_cart_integrals_kinetic_imag] +type: double precision +doc: Imaginary part of the kinetic energy integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[io_ao_cart_integrals_kinetic] +type: Disk_access +doc: Read/Write |AO| kinetic integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + + +[ao_cart_integrals_pseudo] +type: double precision +doc: Pseudopotential integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[ao_cart_integrals_pseudo_imag] +type: double precision +doc: Imaginary part of the pseudopotential integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[io_ao_cart_integrals_pseudo] +type: Disk_access +doc: Read/Write |AO| pseudopotential integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[ao_cart_integrals_overlap] +type: double precision +doc: Overlap integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[ao_cart_integrals_overlap_imag] +type: double precision +doc: Imaginary part of the overlap integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[io_ao_cart_integrals_overlap] +type: Disk_access +doc: Read/Write |AO| overlap integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + +[ao_cart_one_e_integrals] +type: double precision +doc: Combined integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[ao_cart_one_e_integrals_imag] +type: double precision +doc: Imaginary part of the combined integrals in |AO| basis set +size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_num) +interface: ezfio + +[io_ao_cart_one_e_integrals] +type: Disk_access +doc: Read/Write |AO| one-electron integrals from/to disk [ Write | Read | None ] +interface: ezfio,provider,ocaml +default: None + diff --git a/src/ao_cart_one_e_ints/README.rst b/src/ao_cart_one_e_ints/README.rst new file mode 100644 index 00000000..3715e5d7 --- /dev/null +++ b/src/ao_cart_one_e_ints/README.rst @@ -0,0 +1,13 @@ +======================= +ao_cart_one_e_integrals +======================= + +All the one-electron integrals in the cartesian |AO| basis are here. + +The most important providers for usual quantum-chemistry calculation are: + +* `ao_cart_kinetic_integrals` which are the kinetic operator integrals on the |AO| basis +* `ao_cart_integrals_n_e` which are the nuclear-elctron operator integrals on the |AO| basis +* `ao_cart_one_e_integrals` which are the the h_core operator integrals on the |AO| basis + + diff --git a/src/ao_cart_one_e_ints/ao_one_e_ints.irp.f b/src/ao_cart_one_e_ints/ao_one_e_ints.irp.f new file mode 100644 index 00000000..be1a4c1c --- /dev/null +++ b/src/ao_cart_one_e_ints/ao_one_e_ints.irp.f @@ -0,0 +1,47 @@ + BEGIN_PROVIDER [ double precision, ao_cart_one_e_integrals,(ao_cart_num,ao_cart_num)] +&BEGIN_PROVIDER [ double precision, ao_cart_one_e_integrals_diag,(ao_cart_num)] + implicit none + integer :: i,j,n,l + BEGIN_DOC + ! One-electron Hamiltonian in the |AO| basis. + END_DOC + + IF (read_ao_cart_one_e_integrals) THEN + call ezfio_get_ao_cart_one_e_ints_ao_cart_one_e_integrals(ao_cart_one_e_integrals) + ELSE + ao_cart_one_e_integrals = ao_cart_integrals_n_e + ao_cart_kinetic_integrals + + ENDIF + + DO j = 1, ao_cart_num + ao_cart_one_e_integrals_diag(j) = ao_cart_one_e_integrals(j,j) + ENDDO + + IF (write_ao_cart_one_e_integrals) THEN + call ezfio_set_ao_cart_one_e_ints_ao_cart_one_e_integrals(ao_cart_one_e_integrals) + print *, 'AO one-e integrals written to disk' + ENDIF + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_one_e_integrals_imag,(ao_cart_num,ao_cart_num)] + implicit none + integer :: i,j,n,l + BEGIN_DOC + ! One-electron Hamiltonian in the |AO| basis. + END_DOC + + IF (read_ao_cart_one_e_integrals) THEN + call ezfio_get_ao_cart_one_e_ints_ao_cart_one_e_integrals(ao_cart_one_e_integrals_imag) + ELSE + print *, irp_here, ': Not yet implemented' + stop -1 + ENDIF + + IF (write_ao_cart_one_e_integrals) THEN + call ezfio_set_ao_cart_one_e_ints_ao_cart_one_e_integrals(ao_cart_one_e_integrals_imag) + print *, 'AO one-e integrals written to disk' + ENDIF + +END_PROVIDER + diff --git a/src/ao_cart_one_e_ints/ao_overlap.irp.f b/src/ao_cart_one_e_ints/ao_overlap.irp.f new file mode 100644 index 00000000..57dc890a --- /dev/null +++ b/src/ao_cart_one_e_ints/ao_overlap.irp.f @@ -0,0 +1,233 @@ + +! --- + + BEGIN_PROVIDER [double precision, ao_cart_overlap , (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_overlap_x, (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_overlap_y, (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_overlap_z, (ao_cart_num, ao_cart_num)] + + BEGIN_DOC + ! Overlap between atomic basis functions: + ! + ! :math:`\int \chi_i(r) \chi_j(r) dr` + END_DOC + + implicit none + integer :: i, j, n, l, dim1, power_A(3), power_B(3) + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta, c + double precision :: A_center(3), B_center(3) + + ao_cart_overlap = 0.d0 + ao_cart_overlap_x = 0.d0 + ao_cart_overlap_y = 0.d0 + ao_cart_overlap_z = 0.d0 + + if(read_ao_cart_integrals_overlap) then + + call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_overlap(ao_cart_overlap(1:ao_cart_num, 1:ao_cart_num)) + print *, 'AO overlap integrals read from disk' + + else + + if(use_cgtos) then + + do j = 1, ao_cart_num + do i = 1, ao_cart_num + ao_cart_overlap (i,j) = ao_cart_overlap_cgtos (i,j) + ao_cart_overlap_x(i,j) = ao_cart_overlap_cgtos_x(i,j) + ao_cart_overlap_y(i,j) = ao_cart_overlap_cgtos_y(i,j) + ao_cart_overlap_z(i,j) = ao_cart_overlap_cgtos_z(i,j) + enddo + enddo + + else + + dim1=100 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,n,l,c) & + !$OMP SHARED(nucl_coord,ao_cart_power,ao_cart_prim_num, & + !$OMP ao_cart_overlap_x,ao_cart_overlap_y,ao_cart_overlap_z,ao_cart_overlap,ao_cart_num,ao_cart_coef_normalized_ordered_transp,ao_cart_nucl, & + !$OMP ao_cart_expo_ordered_transp,dim1) + do j=1,ao_cart_num + A_center(1) = nucl_coord( ao_cart_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_cart_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_cart_nucl(j), 3 ) + power_A(1) = ao_cart_power( j, 1 ) + power_A(2) = ao_cart_power( j, 2 ) + power_A(3) = ao_cart_power( j, 3 ) + do i= 1,ao_cart_num + B_center(1) = nucl_coord( ao_cart_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_cart_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_cart_nucl(i), 3 ) + power_B(1) = ao_cart_power( i, 1 ) + power_B(2) = ao_cart_power( i, 2 ) + power_B(3) = ao_cart_power( i, 3 ) + do n = 1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(n,j) + do l = 1, ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + c = ao_cart_coef_normalized_ordered_transp(n,j) * ao_cart_coef_normalized_ordered_transp(l,i) + ao_cart_overlap(i,j) += c * overlap + if(isnan(ao_cart_overlap(i,j)))then + print*,'i,j',i,j + print*,'l,n',l,n + print*,'c,overlap',c,overlap + print*,overlap_x,overlap_y,overlap_z + stop + endif + ao_cart_overlap_x(i,j) += c * overlap_x + ao_cart_overlap_y(i,j) += c * overlap_y + ao_cart_overlap_z(i,j) += c * overlap_z + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + + endif + + endif + + if (write_ao_cart_integrals_overlap) then + call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_overlap(ao_cart_overlap(1:ao_cart_num, 1:ao_cart_num)) + print *, 'AO overlap integrals written to disk' + endif + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, ao_cart_overlap_imag, (ao_cart_num, ao_cart_num) ] + implicit none + BEGIN_DOC + ! Imaginary part of the overlap + END_DOC + ao_cart_overlap_imag = 0.d0 +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ complex*16, ao_cart_overlap_complex, (ao_cart_num, ao_cart_num) ] + implicit none + BEGIN_DOC + ! Overlap for complex AOs + END_DOC + integer :: i,j + do j=1,ao_cart_num + do i=1,ao_cart_num + ao_cart_overlap_complex(i,j) = dcmplx( ao_cart_overlap(i,j), ao_cart_overlap_imag(i,j) ) + enddo + enddo +END_PROVIDER + +! --- + +BEGIN_PROVIDER [ double precision, ao_cart_overlap_abs, (ao_cart_num, ao_cart_num) ] + + BEGIN_DOC + ! Overlap between absolute values of atomic basis functions: + ! + ! :math:`\int |\chi_i(r)| |\chi_j(r)| dr` + END_DOC + + implicit none + integer :: i, j, n, l, dim1, power_A(3), power_B(3) + double precision :: overlap_x, overlap_y, overlap_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + double precision :: lower_exp_val, dx + + if(is_periodic) then + + do j = 1, ao_cart_num + do i = 1, ao_cart_num + ao_cart_overlap_abs(i,j) = cdabs(ao_cart_overlap_complex(i,j)) + enddo + enddo + + else + + dim1=100 + lower_exp_val = 40.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B, & + !$OMP overlap_x,overlap_y, overlap_z, & + !$OMP alpha, beta,i,j,dx) & + !$OMP SHARED(nucl_coord,ao_cart_power,ao_cart_prim_num, & + !$OMP ao_cart_overlap_abs,ao_cart_num,ao_cart_coef_normalized_ordered_transp,ao_cart_nucl,& + !$OMP ao_cart_expo_ordered_transp,dim1,lower_exp_val) + do j=1,ao_cart_num + A_center(1) = nucl_coord( ao_cart_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_cart_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_cart_nucl(j), 3 ) + power_A(1) = ao_cart_power( j, 1 ) + power_A(2) = ao_cart_power( j, 2 ) + power_A(3) = ao_cart_power( j, 3 ) + do i= 1,ao_cart_num + ao_cart_overlap_abs(i,j)= 0.d0 + B_center(1) = nucl_coord( ao_cart_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_cart_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_cart_nucl(i), 3 ) + power_B(1) = ao_cart_power( i, 1 ) + power_B(2) = ao_cart_power( i, 2 ) + power_B(3) = ao_cart_power( i, 3 ) + do n = 1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(n,j) + do l = 1, ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(l,i) + call overlap_x_abs(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),overlap_x,lower_exp_val,dx,dim1) + call overlap_x_abs(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),overlap_y,lower_exp_val,dx,dim1) + call overlap_x_abs(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),overlap_z,lower_exp_val,dx,dim1) + ao_cart_overlap_abs(i,j) += abs(ao_cart_coef_normalized_ordered_transp(n,j) * ao_cart_coef_normalized_ordered_transp(l,i)) * overlap_x * overlap_y * overlap_z + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + + endif + +END_PROVIDER + +! --- + + +BEGIN_PROVIDER [ double precision, S_cart_half, (ao_cart_num,ao_cart_num) ] + implicit none + BEGIN_DOC + ! :math:`S^{1/2}` + END_DOC + + integer :: i,j,k + double precision, allocatable :: U(:,:) + double precision, allocatable :: Vt(:,:) + double precision, allocatable :: D(:) + + allocate(U(ao_cart_num,ao_cart_num),Vt(ao_cart_num,ao_cart_num),D(ao_cart_num)) + + call svd(ao_cart_overlap,size(ao_cart_overlap,1),U,size(U,1),D,Vt,size(Vt,1),ao_cart_num,ao_cart_num) + + do i=1,ao_cart_num + D(i) = dsqrt(D(i)) + do j=1,ao_cart_num + S_cart_half(j,i) = 0.d0 + enddo + enddo + + do k=1,ao_cart_num + do j=1,ao_cart_num + do i=1,ao_cart_num + S_cart_half(i,j) = S_cart_half(i,j) + U(i,k)*D(k)*Vt(k,j) + enddo + enddo + enddo + + deallocate(U,Vt,D) + +END_PROVIDER diff --git a/src/ao_cart_one_e_ints/aos_cgtos.irp.f b/src/ao_cart_one_e_ints/aos_cgtos.irp.f new file mode 100644 index 00000000..b9a0ca6f --- /dev/null +++ b/src/ao_cart_one_e_ints/aos_cgtos.irp.f @@ -0,0 +1,124 @@ +! --- + + BEGIN_PROVIDER [double precision, ao_cart_overlap_cgtos, (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_overlap_cgtos_x, (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_overlap_cgtos_y, (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_overlap_cgtos_z, (ao_cart_num, ao_cart_num)] + + implicit none + + integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3) + double precision :: c, overlap, overlap_x, overlap_y, overlap_z + double precision :: KA2(3), phiA(3) + double precision :: KB2(3), phiB(3) + complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) + complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) + complex*16 :: C1(1:4), C2(1:4) + complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1 + complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2 + + ao_cart_overlap_cgtos = 0.d0 + ao_cart_overlap_cgtos_x = 0.d0 + ao_cart_overlap_cgtos_y = 0.d0 + ao_cart_overlap_cgtos_z = 0.d0 + + dim1 = 100 + + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & + !$OMP alpha, alpha_inv, Ae_center, Ap_center, power_A, KA2, phiA, & + !$OMP beta, beta_inv, Be_center, Bp_center, power_B, KB2, phiB, & + !$OMP overlap_x , overlap_y , overlap_z , overlap, & + !$OMP overlap_x1, overlap_y1, overlap_z1, overlap1, & + !$OMP overlap_x2, overlap_y2, overlap_z2, overlap2) & + !$OMP SHARED(nucl_coord, ao_cart_power, ao_cart_prim_num, ao_cart_num, ao_cart_nucl, dim1, & + !$OMP ao_cart_coef_cgtos_norm_ord_transp, ao_cart_expo_cgtos_ord_transp, & + !$OMP ao_cart_expo_pw_ord_transp, ao_cart_expo_phase_ord_transp, & + !$OMP ao_cart_overlap_cgtos_x, ao_cart_overlap_cgtos_y, ao_cart_overlap_cgtos_z, & + !$OMP ao_cart_overlap_cgtos) + + do j = 1, ao_cart_num + + jj = ao_cart_nucl(j) + power_A(1) = ao_cart_power(j,1) + power_A(2) = ao_cart_power(j,2) + power_A(3) = ao_cart_power(j,3) + + do i = 1, ao_cart_num + + ii = ao_cart_nucl(i) + power_B(1) = ao_cart_power(i,1) + power_B(2) = ao_cart_power(i,2) + power_B(3) = ao_cart_power(i,3) + + do n = 1, ao_cart_prim_num(j) + + alpha = ao_cart_expo_cgtos_ord_transp(n,j) + alpha_inv = (1.d0, 0.d0) / alpha + do m = 1, 3 + phiA(m) = ao_cart_expo_phase_ord_transp(m,n,j) + Ap_center(m) = nucl_coord(jj,m) + Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_cart_expo_pw_ord_transp(m,n,j) + KA2(m) = ao_cart_expo_pw_ord_transp(m,n,j) * ao_cart_expo_pw_ord_transp(m,n,j) + enddo + + do l = 1, ao_cart_prim_num(i) + + beta = ao_cart_expo_cgtos_ord_transp(l,i) + beta_inv = (1.d0, 0.d0) / beta + do m = 1, 3 + phiB(m) = ao_cart_expo_phase_ord_transp(m,l,i) + Bp_center(m) = nucl_coord(ii,m) + Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_cart_expo_pw_ord_transp(m,l,i) + KB2(m) = ao_cart_expo_pw_ord_transp(m,l,i) * ao_cart_expo_pw_ord_transp(m,l,i) + enddo + + c = ao_cart_coef_cgtos_norm_ord_transp(n,j) * ao_cart_coef_cgtos_norm_ord_transp(l,i) + + C1(1) = zexp((0.d0, 1.d0) * (-phiA(1) - phiB(1)) - 0.25d0 * (alpha_inv * KA2(1) + beta_inv * KB2(1))) + C1(2) = zexp((0.d0, 1.d0) * (-phiA(2) - phiB(2)) - 0.25d0 * (alpha_inv * KA2(2) + beta_inv * KB2(2))) + C1(3) = zexp((0.d0, 1.d0) * (-phiA(3) - phiB(3)) - 0.25d0 * (alpha_inv * KA2(3) + beta_inv * KB2(3))) + C1(4) = C1(1) * C1(2) * C1(3) + + C2(1) = zexp((0.d0, 1.d0) * (phiA(1) - phiB(1)) - 0.25d0 * (conjg(alpha_inv) * KA2(1) + beta_inv * KB2(1))) + C2(2) = zexp((0.d0, 1.d0) * (phiA(2) - phiB(2)) - 0.25d0 * (conjg(alpha_inv) * KA2(2) + beta_inv * KB2(2))) + C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3))) + C2(4) = C2(1) * C2(2) * C2(3) + + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1) + + call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, & + conjg(Ap_center), Bp_center, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1) + + overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2) + overlap_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2) + overlap_z = 2.d0 * real(C1(3) * overlap_z1 + C2(3) * overlap_z2) + overlap = 2.d0 * real(C1(4) * overlap1 + C2(4) * overlap2 ) + + ao_cart_overlap_cgtos(i,j) = ao_cart_overlap_cgtos(i,j) + c * overlap + + if(isnan(ao_cart_overlap_cgtos(i,j))) then + print*,'i, j', i, j + print*,'l, n', l, n + print*,'c, overlap', c, overlap + print*, overlap_x, overlap_y, overlap_z + stop + endif + + ao_cart_overlap_cgtos_x(i,j) = ao_cart_overlap_cgtos_x(i,j) + c * overlap_x + ao_cart_overlap_cgtos_y(i,j) = ao_cart_overlap_cgtos_y(i,j) + c * overlap_y + ao_cart_overlap_cgtos_z(i,j) = ao_cart_overlap_cgtos_z(i,j) + c * overlap_z + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +! --- + + + diff --git a/src/ao_cart_one_e_ints/kin_ao_ints.irp.f b/src/ao_cart_one_e_ints/kin_ao_ints.irp.f new file mode 100644 index 00000000..4cc60d15 --- /dev/null +++ b/src/ao_cart_one_e_ints/kin_ao_ints.irp.f @@ -0,0 +1,191 @@ + +! --- + + BEGIN_PROVIDER [ double precision, ao_cart_deriv2_x, (ao_cart_num, ao_cart_num) ] +&BEGIN_PROVIDER [ double precision, ao_cart_deriv2_y, (ao_cart_num, ao_cart_num) ] +&BEGIN_PROVIDER [ double precision, ao_cart_deriv2_z, (ao_cart_num, ao_cart_num) ] + + BEGIN_DOC + ! Second derivative matrix elements in the |AO| basis. + ! + ! .. math:: + ! + ! {\tt ao\_deriv2\_x} = + ! \langle \chi_i(x,y,z) | \frac{\partial^2}{\partial x^2} |\chi_j (x,y,z) \rangle + ! + END_DOC + + implicit none + integer :: i, j, n, l, dim1, power_A(3), power_B(3) + double precision :: overlap, overlap_y, overlap_z + double precision :: overlap_x0, overlap_y0, overlap_z0 + double precision :: alpha, beta, c + double precision :: A_center(3), B_center(3) + double precision :: d_a_2,d_2 + + if(use_cgtos) then + + do j = 1, ao_cart_num + do i = 1, ao_cart_num + ao_cart_deriv2_x(i,j) = ao_cart_deriv2_cgtos_x(i,j) + ao_cart_deriv2_y(i,j) = ao_cart_deriv2_cgtos_y(i,j) + ao_cart_deriv2_z(i,j) = ao_cart_deriv2_cgtos_z(i,j) + enddo + enddo + + else + + dim1=100 + + ! -- Dummy call to provide everything + A_center(:) = 0.d0 + B_center(:) = 1.d0 + alpha = 1.d0 + beta = .1d0 + power_A = 1 + power_B = 0 + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1) + ! -- + + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_y, overlap_z, overlap, & + !$OMP alpha, beta, n, l, i,j,c,d_a_2,d_2,deriv_tmp, & + !$OMP overlap_x0,overlap_y0,overlap_z0) & + !$OMP SHARED(nucl_coord,ao_cart_power,ao_cart_prim_num, & + !$OMP ao_cart_deriv2_x,ao_cart_deriv2_y,ao_cart_deriv2_z,ao_cart_num,ao_cart_coef_normalized_ordered_transp,ao_cart_nucl, & + !$OMP ao_cart_expo_ordered_transp,dim1) + do j=1,ao_cart_num + A_center(1) = nucl_coord( ao_cart_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_cart_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_cart_nucl(j), 3 ) + power_A(1) = ao_cart_power( j, 1 ) + power_A(2) = ao_cart_power( j, 2 ) + power_A(3) = ao_cart_power( j, 3 ) + do i= 1,ao_cart_num + ao_cart_deriv2_x(i,j)= 0.d0 + ao_cart_deriv2_y(i,j)= 0.d0 + ao_cart_deriv2_z(i,j)= 0.d0 + B_center(1) = nucl_coord( ao_cart_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_cart_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_cart_nucl(i), 3 ) + power_B(1) = ao_cart_power( i, 1 ) + power_B(2) = ao_cart_power( i, 2 ) + power_B(3) = ao_cart_power( i, 3 ) + do n = 1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(n,j) + do l = 1, ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x0,overlap_y0,overlap_z0,overlap,dim1) + c = ao_cart_coef_normalized_ordered_transp(n,j) * ao_cart_coef_normalized_ordered_transp(l,i) + + power_A(1) = power_A(1)-2 + if (power_A(1)>-1) then + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_a_2,overlap_y,overlap_z,overlap,dim1) + else + d_a_2 = 0.d0 + endif + power_A(1) = power_A(1)+4 + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,d_2,overlap_y,overlap_z,overlap,dim1) + power_A(1) = power_A(1)-2 + + double precision :: deriv_tmp + deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(1)) +1.d0) * overlap_x0 & + +dble(power_A(1)) * (dble(power_A(1))-1.d0) * d_a_2 & + +4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0 + + ao_cart_deriv2_x(i,j) += c*deriv_tmp + power_A(2) = power_A(2)-2 + if (power_A(2)>-1) then + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_a_2,overlap_z,overlap,dim1) + else + d_a_2 = 0.d0 + endif + power_A(2) = power_A(2)+4 + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1) + power_A(2) = power_A(2)-2 + + deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(2)) +1.d0 ) * overlap_y0 & + +dble(power_A(2)) * (dble(power_A(2))-1.d0) * d_a_2 & + +4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0 + ao_cart_deriv2_y(i,j) += c*deriv_tmp + + power_A(3) = power_A(3)-2 + if (power_A(3)>-1) then + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_a_2,overlap,dim1) + else + d_a_2 = 0.d0 + endif + power_A(3) = power_A(3)+4 + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1) + power_A(3) = power_A(3)-2 + + deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(3)) +1.d0 ) * overlap_z0 & + +dble(power_A(3)) * (dble(power_A(3))-1.d0) * d_a_2 & + +4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0 + ao_cart_deriv2_z(i,j) += c*deriv_tmp + + enddo + enddo + enddo + enddo + !$OMP END PARALLEL DO + + endif + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, ao_cart_kinetic_integrals, (ao_cart_num,ao_cart_num)] + implicit none + BEGIN_DOC + ! Kinetic energy integrals in the |AO| basis. + ! + ! $\langle \chi_i |\hat{T}| \chi_j \rangle$ + ! + END_DOC + integer :: i,j,k,l + + if (read_ao_cart_integrals_kinetic) then + call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_kinetic(ao_cart_kinetic_integrals) + print *, 'AO kinetic integrals read from disk' + else + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP PRIVATE(i,j) & + !$OMP SHARED(ao_cart_num, ao_cart_kinetic_integrals,ao_cart_deriv2_x,ao_cart_deriv2_y,ao_cart_deriv2_z) + do j = 1, ao_cart_num + do i = 1, ao_cart_num + ao_cart_kinetic_integrals(i,j) = -0.5d0 * (ao_cart_deriv2_x(i,j) + ao_cart_deriv2_y(i,j) + ao_cart_deriv2_z(i,j) ) + enddo + enddo + !$OMP END PARALLEL DO + endif + if (write_ao_cart_integrals_kinetic) then + call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_kinetic(ao_cart_kinetic_integrals) + print *, 'AO kinetic integrals written to disk' + endif +END_PROVIDER + +BEGIN_PROVIDER [double precision, ao_cart_kinetic_integrals_imag, (ao_cart_num,ao_cart_num)] + implicit none + BEGIN_DOC + ! Kinetic energy integrals in the |AO| basis. + ! + ! $\langle \chi_i |\hat{T}| \chi_j \rangle$ + ! + END_DOC + integer :: i,j,k,l + + if (read_ao_cart_integrals_kinetic) then + call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_kinetic(ao_cart_kinetic_integrals_imag) + print *, 'AO kinetic integrals read from disk' + else + print *, irp_here, ': Not yet implemented' + endif + if (write_ao_cart_integrals_kinetic) then + call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_kinetic(ao_cart_kinetic_integrals_imag) + print *, 'AO kinetic integrals written to disk' + endif +END_PROVIDER diff --git a/src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f b/src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f new file mode 100644 index 00000000..6919b51f --- /dev/null +++ b/src/ao_cart_one_e_ints/one_e_coul_integrals_cgtos.irp.f @@ -0,0 +1,558 @@ + +! --- + +BEGIN_PROVIDER [double precision, ao_cart_integrals_n_e_cgtos, (ao_cart_num, ao_cart_num)] + + BEGIN_DOC + ! + ! Nucleus-electron interaction, in the cgtos |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` + ! + END_DOC + + implicit none + + integer :: power_A(3), power_B(3) + integer :: i, j, k, l, m, n, ii, jj + double precision :: c, Z, C_center(3) + double precision :: phiA, KA2 + double precision :: phiB, KB2 + complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3) + complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3) + complex*16 :: C1, C2, I1, I2 + + complex*16, external :: NAI_pol_mult_cgtos + + + + ao_cart_integrals_n_e_cgtos = 0.d0 + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, C1, C2, I1, I2, & + !$OMP alpha, alpha_inv, Ae_center, Ap_center, phiA, KA2, power_A, & + !$OMP beta, beta_inv, Be_center, Bp_center, phiB, KB2, power_B) & + !$OMP SHARED (ao_cart_num, ao_cart_prim_num, ao_cart_nucl, nucl_coord, & + !$OMP ao_cart_power, nucl_num, nucl_charge, n_pt_max_integrals, & + !$OMP ao_cart_expo_cgtos_ord_transp, ao_cart_coef_cgtos_norm_ord_transp, & + !$OMP ao_cart_expo_pw_ord_transp, ao_cart_expo_phase_ord_transp, & + !$OMP ao_cart_integrals_n_e_cgtos) + !$OMP DO SCHEDULE (dynamic) + + do j = 1, ao_cart_num + + jj = ao_cart_nucl(j) + power_A(1:3) = ao_cart_power(j,1:3) + + do i = 1, ao_cart_num + + ii = ao_cart_nucl(i) + power_B(1:3) = ao_cart_power(i,1:3) + + do n = 1, ao_cart_prim_num(j) + + alpha = ao_cart_expo_cgtos_ord_transp(n,j) + alpha_inv = (1.d0, 0.d0) / alpha + do m = 1, 3 + Ap_center(m) = nucl_coord(jj,m) + Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_cart_expo_pw_ord_transp(m,n,j) + enddo + phiA = ao_cart_expo_phase_ord_transp(4,n,j) + KA2 = ao_cart_expo_pw_ord_transp(4,n,j) + + do l = 1, ao_cart_prim_num(i) + + beta = ao_cart_expo_cgtos_ord_transp(l,i) + beta_inv = (1.d0, 0.d0) / beta + do m = 1, 3 + Bp_center(m) = nucl_coord(ii,m) + Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_cart_expo_pw_ord_transp(m,l,i) + enddo + phiB = ao_cart_expo_phase_ord_transp(4,l,i) + KB2 = ao_cart_expo_pw_ord_transp(4,l,i) + + C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2)) + C2 = zexp((0.d0, 1.d0) * ( phiA - phiB) - 0.25d0 * (conjg(alpha_inv) * KA2 + beta_inv * KB2)) + + c = 0.d0 + do k = 1, nucl_num + + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + I1 = NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, & + Ap_center, Bp_center, C_center, n_pt_max_integrals) + + I2 = NAI_pol_mult_cgtos(conjg(Ae_center), Be_center, power_A, power_B, conjg(alpha), beta, & + conjg(Ap_center), Bp_center, C_center, n_pt_max_integrals) + + c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2) + enddo + + ao_cart_integrals_n_e_cgtos(i,j) += c * ao_cart_coef_cgtos_norm_ord_transp(n,j) & + * ao_cart_coef_cgtos_norm_ord_transp(l,i) + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + +END_PROVIDER + +! --- + +complex*16 function NAI_pol_mult_cgtos(Ae_center, Be_center, power_A, power_B, alpha, beta, & + Ap_center, Bp_center, C_center, n_pt_in) + + BEGIN_DOC + ! + ! Computes the electron-nucleus attraction with two primitves cgtos. + ! + ! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle` + ! + END_DOC + + implicit none + include 'utils/constants.include.F' + + integer, intent(in) :: n_pt_in, power_A(3), power_B(3) + double precision, intent(in) :: C_center(3) + complex*16, intent(in) :: alpha, Ae_center(3), Ap_center(3) + complex*16, intent(in) :: beta, Be_center(3), Bp_center(3) + + integer :: i, n_pt, n_pt_out + double precision :: dist_AB, dist_AC + complex*16 :: p, p_inv, rho, dist, dist_integral, const, const_factor, coeff, factor + complex*16 :: P_center(3) + complex*16 :: d(0:n_pt_in) + + complex*16, external :: V_n_e_cgtos + complex*16, external :: crint_sum + complex*16, external :: crint_1 + + + + dist_AB = 0.d0 + dist_AC = 0.d0 + do i = 1, 3 + dist_AB += abs(Ae_center(i) - Be_center(i)) + dist_AC += abs(Ae_center(i) - C_center(i) * (1.d0, 0.d0)) + enddo + + if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13) .or. use_pw) then + + continue + + else + + NAI_pol_mult_cgtos = V_n_e_cgtos(power_A(1), power_A(2), power_A(3), & + power_B(1), power_B(2), power_B(3), & + alpha, beta) + return + + endif + + p = alpha + beta + p_inv = (1.d0, 0.d0) / p + rho = alpha * beta * p_inv + + dist = (Ae_center(1) - Be_center(1)) * (Ae_center(1) - Be_center(1)) & + + (Ae_center(2) - Be_center(2)) * (Ae_center(2) - Be_center(2)) & + + (Ae_center(3) - Be_center(3)) * (Ae_center(3) - Be_center(3)) + + const_factor = dist * rho + if(real(const_factor) > 80.d0) then + NAI_pol_mult_cgtos = (0.d0, 0.d0) + return + endif + + P_center(1) = (alpha * Ae_center(1) + beta * Be_center(1)) * p_inv + P_center(2) = (alpha * Ae_center(2) + beta * Be_center(2)) * p_inv + P_center(3) = (alpha * Ae_center(3) + beta * Be_center(3)) * p_inv + + dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) & + + (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) & + + (P_center(3) - C_center(3)) * (P_center(3) - C_center(3)) + + const = p * dist_integral + factor = zexp(-const_factor) + coeff = dtwo_pi * factor * p_inv + + n_pt = 2 * ((power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3))) + if(n_pt == 0) then + NAI_pol_mult_cgtos = coeff * crint_1(0, const) + return + endif + + d(0:n_pt_in) = (0.d0, 0.d0) + call give_cpolynomial_mult_center_one_e(Ap_center, Bp_center, alpha, beta, & + power_A, power_B, C_center, n_pt_in, d, n_pt_out) + + if(n_pt_out < 0) then + NAI_pol_mult_cgtos = (0.d0, 0.d0) + return + endif + + NAI_pol_mult_cgtos = coeff * crint_sum(n_pt_out, const, d) + + return +end + +! --- + +subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & + power_A, power_B, C_center, n_pt_in, d, n_pt_out) + + BEGIN_DOC + ! Returns the explicit polynomial in terms of the "t" variable of the following + ! + ! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$. + END_DOC + + implicit none + + integer, intent(in) :: n_pt_in, power_A(3), power_B(3) + double precision, intent(in) :: C_center(3) + complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3) + integer, intent(out) :: n_pt_out + complex*16, intent(inout) :: d(0:n_pt_in) + + integer :: a_x, b_x, a_y, b_y, a_z, b_z + integer :: n_pt1, n_pt2, n_pt3, dim, i, n_pt_tmp + complex*16 :: p, P_center(3), rho, p_inv, p_inv_2 + complex*16 :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2) + complex*16 :: d1(0:n_pt_in), d2(0:n_pt_in), d3(0:n_pt_in) + + ASSERT (n_pt_in > 1) + + p = alpha + beta + p_inv = (1.d0, 0.d0) / p + p_inv_2 = 0.5d0 * p_inv + + P_center(1) = (alpha * A_center(1) + beta * B_center(1)) * p_inv + P_center(2) = (alpha * A_center(2) + beta * B_center(2)) * p_inv + P_center(3) = (alpha * A_center(3) + beta * B_center(3)) * p_inv + + do i = 0, n_pt_in + d(i) = (0.d0, 0.d0) + d1(i) = (0.d0, 0.d0) + d2(i) = (0.d0, 0.d0) + d3(i) = (0.d0, 0.d0) + enddo + + ! --- + + n_pt1 = n_pt_in + + R1x(0) = (P_center(1) - A_center(1)) + R1x(1) = (0.d0, 0.d0) + R1x(2) = -(P_center(1) - C_center(1)) + + R1xp(0) = (P_center(1) - B_center(1)) + R1xp(1) = (0.d0, 0.d0) + R1xp(2) = -(P_center(1) - C_center(1)) + + R2x(0) = p_inv_2 + R2x(1) = (0.d0, 0.d0) + R2x(2) = -p_inv_2 + + a_x = power_A(1) + b_x = power_B(1) + + call I_x1_pol_mult_one_e_cgtos(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in) + + if(n_pt1 < 0) then + n_pt_out = -1 + do i = 0, n_pt_in + d(i) = (0.d0, 0.d0) + enddo + return + endif + + ! --- + + n_pt2 = n_pt_in + + R1x(0) = (P_center(2) - A_center(2)) + R1x(1) = (0.d0, 0.d0) + R1x(2) = -(P_center(2) - C_center(2)) + + R1xp(0) = (P_center(2) - B_center(2)) + R1xp(1) = (0.d0, 0.d0) + R1xp(2) = -(P_center(2) - C_center(2)) + + a_y = power_A(2) + b_y = power_B(2) + + call I_x1_pol_mult_one_e_cgtos(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in) + + if(n_pt2 < 0) then + n_pt_out = -1 + do i = 0, n_pt_in + d(i) = (0.d0, 0.d0) + enddo + return + endif + + ! --- + + n_pt3 = n_pt_in + + R1x(0) = (P_center(3) - A_center(3)) + R1x(1) = (0.d0, 0.d0) + R1x(2) = -(P_center(3) - C_center(3)) + + R1xp(0) = (P_center(3) - B_center(3)) + R1xp(1) = (0.d0, 0.d0) + R1xp(2) = -(P_center(3) - C_center(3)) + + a_z = power_A(3) + b_z = power_B(3) + + call I_x1_pol_mult_one_e_cgtos(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in) + + if(n_pt3 < 0) then + n_pt_out = -1 + do i = 0, n_pt_in + d(i) = (0.d0, 0.d0) + enddo + return + endif + + ! --- + + n_pt_tmp = 0 + call multiply_cpoly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp) + + n_pt_out = 0 + d1(0:n_pt_tmp) = (0.d0, 0.d0) + call multiply_cpoly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out) + do i = 0, n_pt_out + d(i) = d1(i) + enddo + +end + +! --- + +recursive subroutine I_x1_pol_mult_one_e_cgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt_in) + + BEGIN_DOC + ! Recursive routine involved in the electron-nucleus potential + END_DOC + + implicit none + include 'utils/constants.include.F' + + integer, intent(in) :: a, c, n_pt_in + complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2) + integer, intent(inout) :: nd + complex*16, intent(inout) :: d(0:n_pt_in) + + integer :: nx, ix, dim, iy, ny + complex*16 :: X(0:max_dim) + complex*16 :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y + + dim = n_pt_in + + if((a == 0) .and. (c == 0)) then + + nd = 0 + d(0) = (1.d0, 0.d0) + return + + elseif((c < 0) .or. (nd < 0)) then + + nd = -1 + return + + elseif((a == 0) .and. (c .ne. 0)) then + + call I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, d, nd, n_pt_in) + + elseif(a == 1) then + + nx = nd + do ix = 0, n_pt_in + X(ix) = (0.d0, 0.d0) + Y(ix) = (0.d0, 0.d0) + enddo + + call I_x2_pol_mult_one_e_cgtos(c-1, R1x, R1xp, R2x, X, nx, n_pt_in) + + do ix = 0, nx + X(ix) *= dble(c) + enddo + + call multiply_cpoly(X, nx, R2x, 2, d, nd) + + ny = 0 + call I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, Y, ny, n_pt_in) + call multiply_cpoly(Y, ny, R1x, 2, d, nd) + + else + + nx = 0 + do ix = 0, n_pt_in + X(ix) = (0.d0, 0.d0) + Y(ix) = (0.d0, 0.d0) + enddo + + call I_x1_pol_mult_one_e_cgtos(a-2, c, R1x, R1xp, R2x, X, nx, n_pt_in) + + do ix = 0, nx + X(ix) *= dble(a-1) + enddo + call multiply_cpoly(X, nx, R2x, 2, d, nd) + + nx = nd + do ix = 0, n_pt_in + X(ix) = (0.d0, 0.d0) + enddo + + call I_x1_pol_mult_one_e_cgtos(a-1, c-1, R1x, R1xp, R2x, X, nx, n_pt_in) + do ix = 0, nx + X(ix) *= dble(c) + enddo + + call multiply_cpoly(X, nx, R2x, 2, d, nd) + + ny = 0 + call I_x1_pol_mult_one_e_cgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in) + call multiply_cpoly(Y, ny, R1x, 2, d, nd) + + endif + +end + +! --- + +recursive subroutine I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, d, nd, dim) + + BEGIN_DOC + ! Recursive routine involved in the electron-nucleus potential + END_DOC + + implicit none + include 'utils/constants.include.F' + + integer, intent(in) :: dim, c + complex*16, intent(in) :: R1x(0:2), R1xp(0:2), R2x(0:2) + integer, intent(inout) :: nd + complex*16, intent(out) :: d(0:max_dim) + + integer :: i, nx, ix, ny + complex*16 :: X(0:max_dim), Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y + + if(c == 0) then + + nd = 0 + d(0) = (1.d0, 0.d0) + return + + elseif((nd < 0) .or. (c < 0)) then + + nd = -1 + return + + else + + nx = 0 + do ix = 0, dim + X(ix) = (0.d0, 0.d0) + Y(ix) = (0.d0, 0.d0) + enddo + + call I_x1_pol_mult_one_e_cgtos(0, c-2, R1x, R1xp, R2x, X, nx, dim) + + do ix = 0, nx + X(ix) *= dble(c-1) + enddo + + call multiply_cpoly(X, nx, R2x, 2, d, nd) + + ny = 0 + do ix = 0, dim + Y(ix) = (0.d0, 0.d0) + enddo + + call I_x1_pol_mult_one_e_cgtos(0, c-1, R1x, R1xp, R2x, Y, ny, dim) + + if(ny .ge. 0) then + call multiply_cpoly(Y, ny, R1xp, 2, d, nd) + endif + + endif + +end + +! --- + +complex*16 function V_n_e_cgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta) + + BEGIN_DOC + ! Primitve nuclear attraction between the two primitves centered on the same atom. + ! + ! $p_1 = x^{a_x} y^{a_y} z^{a_z} \exp(-\alpha r^2)$ + ! + ! $p_2 = x^{b_x} y^{b_y} z^{b_z} \exp(-\beta r^2)$ + END_DOC + + implicit none + + integer, intent(in) :: a_x, a_y, a_z, b_x, b_y, b_z + complex*16, intent(in) :: alpha, beta + + double precision :: V_phi, V_theta + complex*16 :: V_r_cgtos + + if( (iand(a_x + b_x, 1) == 1) .or. & + (iand(a_y + b_y, 1) == 1) .or. & + (iand(a_z + b_z, 1) == 1) ) then + + V_n_e_cgtos = (0.d0, 0.d0) + + else + + V_n_e_cgtos = V_r_cgtos(a_x + b_x + a_y + b_y + a_z + b_z + 1, alpha + beta) & + * V_phi(a_x + b_x, a_y + b_y) & + * V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1) + endif + +end + +! --- + +complex*16 function V_r_cgtos(n, alpha) + + BEGIN_DOC + ! Computes the radial part of the nuclear attraction integral: + ! + ! $\int_{0}^{\infty} r^n \exp(-\alpha r^2) dr$ + ! + END_DOC + + implicit none + include 'utils/constants.include.F' + + integer , intent(in) :: n + complex*16, intent(in) :: alpha + + double precision :: fact + + if(iand(n, 1) .eq. 1) then + V_r_cgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1)) + else + V_r_cgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1) + endif + +end + +! --- + diff --git a/src/ao_cart_one_e_ints/one_e_kin_integrals_cgtos.irp.f b/src/ao_cart_one_e_ints/one_e_kin_integrals_cgtos.irp.f new file mode 100644 index 00000000..58104210 --- /dev/null +++ b/src/ao_cart_one_e_ints/one_e_kin_integrals_cgtos.irp.f @@ -0,0 +1,318 @@ + +! --- + + BEGIN_PROVIDER [double precision, ao_cart_deriv2_cgtos_x, (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_deriv2_cgtos_y, (ao_cart_num, ao_cart_num)] +&BEGIN_PROVIDER [double precision, ao_cart_deriv2_cgtos_z, (ao_cart_num, ao_cart_num)] + + implicit none + integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3) + double precision :: c, deriv_tmp + double precision :: KA2, phiA + double precision :: KB2, phiB + double precision :: aa + complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3), C1 + complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3), C2 + complex*16 :: xa + complex*16 :: overlap_x, overlap_y, overlap_z, overlap + complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1 + complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2 + complex*16 :: overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1 + complex*16 :: overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2 + complex*16 :: deriv_tmp_1, deriv_tmp_2 + + + dim1 = 100 + + ! -- Dummy call to provide everything + + Ae_center(:) = (0.0d0, 0.d0) + Be_center(:) = (1.0d0, 0.d0) + Ap_center(:) = (0.0d0, 0.d0) + Bp_center(:) = (1.0d0, 0.d0) + alpha = (1.0d0, 0.d0) + beta = (0.1d0, 0.d0) + power_A = 1 + power_B = 0 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) + + ! --- + + !$OMP PARALLEL & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, aa, xa, C1, C2, & + !$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, & + !$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, & + !$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, & + !$OMP overlap_x, overlap_y, overlap_z, overlap, & + !$OMP overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1, & + !$OMP overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2, & + !$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, & + !$OMP overlap_y0_2, overlap_z0_2) & + !$OMP SHARED(nucl_coord, ao_cart_power, ao_cart_prim_num, ao_cart_num, ao_cart_nucl, dim1, & + !$OMP ao_cart_coef_cgtos_norm_ord_transp, ao_cart_expo_cgtos_ord_transp, & + !$OMP ao_cart_expo_pw_ord_transp, ao_cart_expo_phase_ord_transp, & + !$OMP ao_cart_deriv2_cgtos_x, ao_cart_deriv2_cgtos_y, ao_cart_deriv2_cgtos_z) + !$OMP DO SCHEDULE(GUIDED) + do j = 1, ao_cart_num + + jj = ao_cart_nucl(j) + power_A(1) = ao_cart_power(j,1) + power_A(2) = ao_cart_power(j,2) + power_A(3) = ao_cart_power(j,3) + + do i = 1, ao_cart_num + + ii = ao_cart_nucl(i) + power_B(1) = ao_cart_power(i,1) + power_B(2) = ao_cart_power(i,2) + power_B(3) = ao_cart_power(i,3) + + ao_cart_deriv2_cgtos_x(i,j) = 0.d0 + ao_cart_deriv2_cgtos_y(i,j) = 0.d0 + ao_cart_deriv2_cgtos_z(i,j) = 0.d0 + + do n = 1, ao_cart_prim_num(j) + + alpha = ao_cart_expo_cgtos_ord_transp(n,j) + alpha_inv = (1.d0, 0.d0) / alpha + do m = 1, 3 + Ap_center(m) = nucl_coord(jj,m) + Ae_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_cart_expo_pw_ord_transp(m,n,j) + enddo + phiA = ao_cart_expo_phase_ord_transp(4,n,j) + KA2 = ao_cart_expo_pw_ord_transp(4,n,j) + + do l = 1, ao_cart_prim_num(i) + + beta = ao_cart_expo_cgtos_ord_transp(l,i) + beta_inv = (1.d0, 0.d0) / beta + do m = 1, 3 + Bp_center(m) = nucl_coord(ii,m) + Be_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_cart_expo_pw_ord_transp(m,l,i) + enddo + phiB = ao_cart_expo_phase_ord_transp(4,l,i) + KB2 = ao_cart_expo_pw_ord_transp(4,l,i) + + c = ao_cart_coef_cgtos_norm_ord_transp(n,j) * ao_cart_coef_cgtos_norm_ord_transp(l,i) + + C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2)) + C2 = zexp((0.d0, 1.d0) * (-phiA + phiB) - 0.25d0 * (alpha_inv * KA2 + conjg(beta_inv) * KB2)) + + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1) + + ! --- + + power_A(1) = power_A(1) - 1 + if(power_A(1) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_m1_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_m1_2, overlap_y, overlap_z, overlap, dim1) + power_A(1) = power_A(1) - 1 + if(power_A(1) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_m2_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_m2_2, overlap_y, overlap_z, overlap, dim1) + else + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(1) = power_A(1) + 1 + else + overlap_m1_1 = (0.d0, 0.d0) + overlap_m1_2 = (0.d0, 0.d0) + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(1) = power_A(1) + 1 + + power_A(1) = power_A(1) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_p1_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_p1_2, overlap_y, overlap_z, overlap, dim1) + power_A(1) = power_A(1) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_p2_1, overlap_y, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_p2_2, overlap_y, overlap_z, overlap, dim1) + power_A(1) = power_A(1) - 2 + + aa = dble(power_A(1)) + xa = Ap_center(1) - Ae_center(1) + + deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_1 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) + deriv_tmp_1 = deriv_tmp_1 * overlap_y0_1 * overlap_z0_1 + + deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_2 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) + deriv_tmp_2 = deriv_tmp_2 * overlap_y0_2 * overlap_z0_2 + + deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) + + ao_cart_deriv2_cgtos_x(i,j) += c * deriv_tmp + + ! --- + + power_A(2) = power_A(2) - 1 + if(power_A(2) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_m1_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_m1_2, overlap_z, overlap, dim1) + power_A(2) = power_A(2) - 1 + if(power_A(2) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_m2_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_m2_2, overlap_z, overlap, dim1) + else + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(2) = power_A(2) + 1 + else + overlap_m1_1 = (0.d0, 0.d0) + overlap_m1_2 = (0.d0, 0.d0) + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(2) = power_A(2) + 1 + + power_A(2) = power_A(2) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_p1_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_p1_2, overlap_z, overlap, dim1) + power_A(2) = power_A(2) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_p2_1, overlap_z, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_p2_2, overlap_z, overlap, dim1) + power_A(2) = power_A(2) - 2 + + aa = dble(power_A(2)) + xa = Ap_center(2) - Ae_center(2) + + deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_1 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) + deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_z0_1 + + deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_2 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) + deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_z0_2 + + deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) + + ao_cart_deriv2_cgtos_y(i,j) += c * deriv_tmp + + ! --- + + power_A(3) = power_A(3) - 1 + if(power_A(3) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_m1_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m1_2, overlap, dim1) + power_A(3) = power_A(3) - 1 + if(power_A(3) > -1) then + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_m2_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m2_2, overlap, dim1) + else + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(3) = power_A(3) + 1 + else + overlap_m1_1 = (0.d0, 0.d0) + overlap_m1_2 = (0.d0, 0.d0) + overlap_m2_1 = (0.d0, 0.d0) + overlap_m2_2 = (0.d0, 0.d0) + endif + power_A(3) = power_A(3) + 1 + + power_A(3) = power_A(3) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_p1_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p1_2, overlap, dim1) + power_A(3) = power_A(3) + 1 + call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, & + Ap_center, Bp_center, overlap_x, overlap_y, overlap_p2_1, overlap, dim1) + call overlap_cgaussian_xyz(Ae_center, conjg(Be_center), alpha, conjg(beta), power_A, power_B, & + Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p2_2, overlap, dim1) + power_A(3) = power_A(3) - 2 + + aa = dble(power_A(3)) + xa = Ap_center(3) - Ae_center(3) + + deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_1 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_1 + 0.5d0 * overlap_p2_1) + deriv_tmp_1 = deriv_tmp_1 * overlap_x0_1 * overlap_y0_1 + + deriv_tmp_2 = aa * (aa - 1.d0) * overlap_m2_2 - 4.d0 * alpha * aa * xa * overlap_m1_2 & + + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_2 & + + 8.d0 * alpha * alpha * (xa * overlap_p1_2 + 0.5d0 * overlap_p2_2) + deriv_tmp_2 = deriv_tmp_2 * overlap_x0_2 * overlap_y0_2 + + deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2) + + ao_cart_deriv2_cgtos_z(i,j) += c * deriv_tmp + + ! --- + + enddo + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + +END_PROVIDER + +! --- + +BEGIN_PROVIDER [double precision, ao_cart_kinetic_integrals_cgtos, (ao_cart_num, ao_cart_num)] + + BEGIN_DOC + ! + ! Kinetic energy integrals in the cgtos |AO| basis. + ! + ! $\langle \chi_i |\hat{T}| \chi_j \rangle$ + ! + END_DOC + + implicit none + + integer :: i, j + + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP PRIVATE(i, j) & + !$OMP SHARED(ao_cart_num, ao_cart_kinetic_integrals_cgtos, ao_cart_deriv2_cgtos_x, ao_cart_deriv2_cgtos_y, ao_cart_deriv2_cgtos_z) + do j = 1, ao_cart_num + do i = 1, ao_cart_num + ao_cart_kinetic_integrals_cgtos(i,j) = -0.5d0 * (ao_cart_deriv2_cgtos_x(i,j) + & + ao_cart_deriv2_cgtos_y(i,j) + & + ao_cart_deriv2_cgtos_z(i,j)) + enddo + enddo + !$OMP END PARALLEL DO + +END_PROVIDER + +! --- + diff --git a/src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f b/src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f new file mode 100644 index 00000000..a81206c1 --- /dev/null +++ b/src/ao_cart_one_e_ints/pot_ao_erf_ints.irp.f @@ -0,0 +1,794 @@ + +! --- + +subroutine give_all_erf_kl_ao(integrals_ao,mu_in,C_center) + implicit none + BEGIN_DOC + ! Subroutine that returns all integrals over $r$ of type + ! $\frac{ \erf(\mu * | r - R_C | ) }{ | r - R_C | }$ + END_DOC + double precision, intent(in) :: mu_in,C_center(3) + double precision, intent(out) :: integrals_ao(ao_cart_num,ao_cart_num) + double precision :: NAI_pol_mult_erf_ao + integer :: i,j,l,k,m + do k = 1, ao_cart_num + do m = 1, ao_cart_num + integrals_ao(m,k) = NAI_pol_mult_erf_ao(m,k,mu_in,C_center) + enddo + enddo +end + +! --- + +double precision function NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center) + + BEGIN_DOC + ! + ! Computes the following integral : + ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + END_DOC + + implicit none + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: mu_in, C_center(3) + + integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in + double precision :: A_center(3), B_center(3), integral, alpha, beta + + double precision :: NAI_pol_mult_erf + + num_A = ao_cart_nucl(i_ao) + power_A(1:3) = ao_cart_power(i_ao,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_cart_nucl(j_ao) + power_B(1:3) = ao_cart_power(j_ao,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + n_pt_in = n_pt_max_integrals + + NAI_pol_mult_erf_ao = 0.d0 + do i = 1, ao_cart_prim_num(i_ao) + alpha = ao_cart_expo_ordered_transp(i,i_ao) + do j = 1, ao_cart_prim_num(j_ao) + beta = ao_cart_expo_ordered_transp(j,j_ao) + + integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in,mu_in) + + NAI_pol_mult_erf_ao += integral * ao_cart_coef_normalized_ordered_transp(j,j_ao) * ao_cart_coef_normalized_ordered_transp(i,i_ao) + enddo + enddo + +end function NAI_pol_mult_erf_ao + +! --- + +double precision function NAI_pol_mult_erf_ao_cart_with1s(i_ao, j_ao, beta, B_center, mu_in, C_center) + + BEGIN_DOC + ! + ! Computes the following integral : + ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) e^{-\beta (r - B_center)^2} \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + END_DOC + + implicit none + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: beta, B_center(3) + double precision, intent(in) :: mu_in, C_center(3) + + integer :: i, j, power_A1(3), power_A2(3), n_pt_in + double precision :: A1_center(3), A2_center(3), alpha1, alpha2, coef12, coef1, integral + + double precision, external :: NAI_pol_mult_erf_with1s, NAI_pol_mult_erf_ao + + ASSERT(beta .ge. 0.d0) + if(beta .lt. 1d-10) then + NAI_pol_mult_erf_ao_cart_with1s = NAI_pol_mult_erf_ao(i_ao, j_ao, mu_in, C_center) + return + endif + + power_A1(1:3) = ao_cart_power(i_ao,1:3) + power_A2(1:3) = ao_cart_power(j_ao,1:3) + + A1_center(1:3) = nucl_coord(ao_cart_nucl(i_ao),1:3) + A2_center(1:3) = nucl_coord(ao_cart_nucl(j_ao),1:3) + + n_pt_in = n_pt_max_integrals + + NAI_pol_mult_erf_ao_cart_with1s = 0.d0 + do i = 1, ao_cart_prim_num(i_ao) + alpha1 = ao_cart_expo_ordered_transp (i,i_ao) + coef1 = ao_cart_coef_normalized_ordered_transp(i,i_ao) + + do j = 1, ao_cart_prim_num(j_ao) + alpha2 = ao_cart_expo_ordered_transp(j,j_ao) + coef12 = coef1 * ao_cart_coef_normalized_ordered_transp(j,j_ao) + if(dabs(coef12) .lt. 1d-14) cycle + + integral = NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 & + , beta, B_center, C_center, n_pt_in, mu_in ) + + NAI_pol_mult_erf_ao_cart_with1s += integral * coef12 + enddo + enddo + +end function NAI_pol_mult_erf_ao_cart_with1s + +! --- + +double precision function NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in, mu_in) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! .. math:: + ! + ! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 ) + ! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + integer, intent(in) :: n_pt_in + integer, intent(in) :: power_A(3), power_B(3) + double precision, intent(in) :: C_center(3), A_center(3), B_center(3), alpha, beta, mu_in + + integer :: i, n_pt, n_pt_out + double precision :: P_center(3) + double precision :: d(0:n_pt_in), coeff, dist, const, factor + double precision :: const_factor, dist_integral + double precision :: accu, p_inv, p, rho, p_inv_2 + double precision :: p_new + + double precision :: rint + + p = alpha + beta + p_inv = 1.d0 / p + p_inv_2 = 0.5d0 * p_inv + rho = alpha * beta * p_inv + + dist = 0.d0 + dist_integral = 0.d0 + do i = 1, 3 + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i)) + dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) + enddo + const_factor = dist * rho + if(const_factor > 80.d0) then + NAI_pol_mult_erf = 0.d0 + return + endif + + p_new = mu_in / dsqrt(p + mu_in * mu_in) + factor = dexp(-const_factor) + coeff = dtwo_pi * factor * p_inv * p_new + + n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) ) + const = p * dist_integral * p_new * p_new + if(n_pt == 0) then + NAI_pol_mult_erf = coeff * rint(0, const) + return + endif + + do i = 0, n_pt_in + d(i) = 0.d0 + enddo + ! call give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in) + p_new = p_new * p_new + call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center) + + if(n_pt_out < 0) then + NAI_pol_mult_erf = 0.d0 + return + endif + + ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i + accu = 0.d0 + do i = 0, n_pt_out, 2 + accu += d(i) * rint(i/2, const) + enddo + NAI_pol_mult_erf = accu * coeff + +end function NAI_pol_mult_erf + +! --- + +subroutine NAI_pol_mult_erf_v(A_center, B_center, power_A, power_B, alpha, beta, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! .. math:: + ! + ! \int dr (x-A_x)^a (x-B_x)^b \exp(-\alpha (x-A_x)^2 - \beta (x-B_x)^2 ) + ! \frac{\erf(\mu |r - R_C |)}{| r - R_C |}$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + + integer, intent(in) :: n_pt_in, n_points, LD_C, LD_resv + integer, intent(in) :: power_A(3), power_B(3) + double precision, intent(in) :: A_center(3), B_center(3), alpha, beta, mu_in + double precision, intent(in) :: C_center(LD_C,3) + double precision, intent(out) :: res_v(LD_resv) + + integer :: i, n_pt, n_pt_out, ipoint + double precision :: P_center(3) + double precision :: d(0:n_pt_in), coeff, dist, const, factor + double precision :: const_factor, dist_integral + double precision :: accu, p_inv, p, rho, p_inv_2 + double precision :: p_new, p_new2, coef_tmp + + double precision :: rint + + res_V(1:LD_resv) = 0.d0 + + p = alpha + beta + p_inv = 1.d0 / p + p_inv_2 = 0.5d0 * p_inv + rho = alpha * beta * p_inv + p_new = mu_in / dsqrt(p + mu_in * mu_in) + p_new2 = p_new * p_new + coef_tmp = p * p_new2 + + dist = 0.d0 + do i = 1, 3 + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i)) + enddo + + const_factor = dist * rho + if(const_factor > 80.d0) then + return + endif + factor = dexp(-const_factor) + coeff = dtwo_pi * factor * p_inv * p_new + + n_pt = 2 * ( power_A(1) + power_B(1) + power_A(2) + power_B(2) + power_A(3) + power_B(3) ) + + if(n_pt == 0) then + + do ipoint = 1, n_points + dist_integral = 0.d0 + do i = 1, 3 + dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i)) + enddo + const = coef_tmp * dist_integral + + res_v(ipoint) = coeff * rint(0, const) + enddo + + else + + do ipoint = 1, n_points + dist_integral = 0.d0 + do i = 1, 3 + dist_integral += (P_center(i) - C_center(ipoint,i)) * (P_center(i) - C_center(ipoint,i)) + enddo + const = coef_tmp * dist_integral + + do i = 0, n_pt_in + d(i) = 0.d0 + enddo + call give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center) + + if(n_pt_out < 0) then + res_v(ipoint) = 0.d0 + cycle + endif + + ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i + accu = 0.d0 + do i = 0, n_pt_out, 2 + accu += d(i) * rint(i/2, const) + enddo + + res_v(ipoint) = accu * coeff + enddo + + endif + +end subroutine NAI_pol_mult_erf_v + +! --- + +double precision function NAI_pol_mult_erf_with1s( A1_center, A2_center, power_A1, power_A2, alpha1, alpha2 & + , beta, B_center, C_center, n_pt_in, mu_in ) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! .. math:: + ! + ! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2) + ! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2) + ! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2) + ! \exp(-\beta (r - B)^2) + ! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + integer, intent(in) :: n_pt_in + integer, intent(in) :: power_A1(3), power_A2(3) + double precision, intent(in) :: C_center(3), A1_center(3), A2_center(3), B_center(3) + double precision, intent(in) :: alpha1, alpha2, beta, mu_in + + integer :: i, n_pt, n_pt_out + double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12 + double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor + double precision :: dist_integral + double precision :: d(0:n_pt_in), coeff, const, factor + double precision :: accu + double precision :: p_new + + double precision :: rint + + + ! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2} + alpha12 = alpha1 + alpha2 + alpha12_inv = 1.d0 / alpha12 + alpha12_inv_2 = 0.5d0 * alpha12_inv + rho12 = alpha1 * alpha2 * alpha12_inv + A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv + A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv + A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv + dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1)) & + + (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2)) & + + (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) + + const_factor12 = dist12 * rho12 + if(const_factor12 > 80.d0) then + NAI_pol_mult_erf_with1s = 0.d0 + return + endif + + ! --- + + ! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2} + p = alpha12 + beta + p_inv = 1.d0 / p + p_inv_2 = 0.5d0 * p_inv + rho = alpha12 * beta * p_inv + P_center(1) = (alpha12 * A12_center(1) + beta * B_center(1)) * p_inv + P_center(2) = (alpha12 * A12_center(2) + beta * B_center(2)) * p_inv + P_center(3) = (alpha12 * A12_center(3) + beta * B_center(3)) * p_inv + dist = (A12_center(1) - B_center(1)) * (A12_center(1) - B_center(1)) & + + (A12_center(2) - B_center(2)) * (A12_center(2) - B_center(2)) & + + (A12_center(3) - B_center(3)) * (A12_center(3) - B_center(3)) + + const_factor = const_factor12 + dist * rho + if(const_factor > 80.d0) then + NAI_pol_mult_erf_with1s = 0.d0 + return + endif + + dist_integral = (P_center(1) - C_center(1)) * (P_center(1) - C_center(1)) & + + (P_center(2) - C_center(2)) * (P_center(2) - C_center(2)) & + + (P_center(3) - C_center(3)) * (P_center(3) - C_center(3)) + + ! --- + + p_new = mu_in / dsqrt(p + mu_in * mu_in) + factor = dexp(-const_factor) + coeff = dtwo_pi * factor * p_inv * p_new + + n_pt = 2 * ( (power_A1(1) + power_A2(1)) + (power_A1(2) + power_A2(2)) + (power_A1(3) + power_A2(3)) ) + const = p * dist_integral * p_new * p_new + if(n_pt == 0) then + NAI_pol_mult_erf_with1s = coeff * rint(0, const) + return + endif + + do i = 0, n_pt_in + d(i) = 0.d0 + enddo + p_new = p_new * p_new + call give_polynomial_mult_center_one_e_erf_opt( A1_center, A2_center, power_A1, power_A2, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center) + + if(n_pt_out < 0) then + NAI_pol_mult_erf_with1s = 0.d0 + return + endif + + ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i + accu = 0.d0 + do i = 0, n_pt_out, 2 + accu += d(i) * rint(i/2, const) + enddo + NAI_pol_mult_erf_with1s = accu * coeff + +end function NAI_pol_mult_erf_with1s + +! --- + +subroutine NAI_pol_mult_erf_with1s_v(A1_center, A2_center, power_A1, power_A2, alpha1, alpha2, beta, B_center, LD_B, C_center, LD_C, n_pt_in, mu_in, res_v, LD_resv, n_points) + + BEGIN_DOC + ! + ! Computes the following integral : + ! + ! .. math :: + ! + ! \int dx (x - A1_x)^a_1 (x - B1_x)^a_2 \exp(-\alpha_1 (x - A1_x)^2 - \alpha_2 (x - A2_x)^2) + ! \int dy (y - A1_y)^b_1 (y - B1_y)^b_2 \exp(-\alpha_1 (y - A1_y)^2 - \alpha_2 (y - A2_y)^2) + ! \int dz (x - A1_z)^c_1 (z - B1_z)^c_2 \exp(-\alpha_1 (z - A1_z)^2 - \alpha_2 (z - A2_z)^2) + ! \exp(-\beta (r - B)^2) + ! \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + END_DOC + + include 'utils/constants.include.F' + + implicit none + integer, intent(in) :: n_pt_in, LD_B, LD_C, LD_resv, n_points + integer, intent(in) :: power_A1(3), power_A2(3) + double precision, intent(in) :: A1_center(3), A2_center(3) + double precision, intent(in) :: C_center(LD_C,3), B_center(LD_B,3) + double precision, intent(in) :: alpha1, alpha2, beta, mu_in + double precision, intent(out) :: res_v(LD_resv) + + integer :: i, n_pt, n_pt_out, ipoint + double precision :: alpha12, alpha12_inv, alpha12_inv_2, rho12, A12_center(3), dist12, const_factor12 + double precision :: p, p_inv, p_inv_2, rho, P_center(3), dist, const_factor + double precision :: dist_integral + double precision :: d(0:n_pt_in), coeff, const, factor + double precision :: accu + double precision :: p_new, p_new2, coef_tmp, cons_tmp + + double precision :: rint + + + res_V(1:LD_resv) = 0.d0 + + ! e^{-alpha1 (r - A1)^2} e^{-alpha2 (r - A2)^2} = e^{-K12} e^{-alpha12 (r - A12)^2} + alpha12 = alpha1 + alpha2 + alpha12_inv = 1.d0 / alpha12 + alpha12_inv_2 = 0.5d0 * alpha12_inv + rho12 = alpha1 * alpha2 * alpha12_inv + A12_center(1) = (alpha1 * A1_center(1) + alpha2 * A2_center(1)) * alpha12_inv + A12_center(2) = (alpha1 * A1_center(2) + alpha2 * A2_center(2)) * alpha12_inv + A12_center(3) = (alpha1 * A1_center(3) + alpha2 * A2_center(3)) * alpha12_inv + dist12 = (A1_center(1) - A2_center(1)) * (A1_center(1) - A2_center(1))& + + (A1_center(2) - A2_center(2)) * (A1_center(2) - A2_center(2))& + + (A1_center(3) - A2_center(3)) * (A1_center(3) - A2_center(3)) + + const_factor12 = dist12 * rho12 + if(const_factor12 > 80.d0) then + return + endif + + ! e^{-K12} e^{-alpha12 (r - A12)^2} e^{-beta (r - B)^2} = e^{-K} e^{-p (r - P)^2} + p = alpha12 + beta + p_inv = 1.d0 / p + p_inv_2 = 0.5d0 * p_inv + rho = alpha12 * beta * p_inv + p_new = mu_in / dsqrt(p + mu_in * mu_in) + p_new2 = p_new * p_new + coef_tmp = dtwo_pi * p_inv * p_new + cons_tmp = p * p_new2 + n_pt = 2 * (power_A1(1) + power_A2(1) + power_A1(2) + power_A2(2) + power_A1(3) + power_A2(3) ) + + if(n_pt == 0) then + + do ipoint = 1, n_points + + dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))& + + (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))& + + (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3)) + const_factor = const_factor12 + dist * rho + if(const_factor > 80.d0) cycle + coeff = coef_tmp * dexp(-const_factor) + + P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv + P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv + P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv + dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))& + + (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))& + + (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3)) + const = cons_tmp * dist_integral + + res_v(ipoint) = coeff * rint(0, const) + enddo + + else + + do ipoint = 1, n_points + + dist = (A12_center(1) - B_center(ipoint,1)) * (A12_center(1) - B_center(ipoint,1))& + + (A12_center(2) - B_center(ipoint,2)) * (A12_center(2) - B_center(ipoint,2))& + + (A12_center(3) - B_center(ipoint,3)) * (A12_center(3) - B_center(ipoint,3)) + const_factor = const_factor12 + dist * rho + if(const_factor > 80.d0) cycle + coeff = coef_tmp * dexp(-const_factor) + + P_center(1) = (alpha12 * A12_center(1) + beta * B_center(ipoint,1)) * p_inv + P_center(2) = (alpha12 * A12_center(2) + beta * B_center(ipoint,2)) * p_inv + P_center(3) = (alpha12 * A12_center(3) + beta * B_center(ipoint,3)) * p_inv + dist_integral = (P_center(1) - C_center(ipoint,1)) * (P_center(1) - C_center(ipoint,1))& + + (P_center(2) - C_center(ipoint,2)) * (P_center(2) - C_center(ipoint,2))& + + (P_center(3) - C_center(ipoint,3)) * (P_center(3) - C_center(ipoint,3)) + const = cons_tmp * dist_integral + + do i = 0, n_pt_in + d(i) = 0.d0 + enddo + !TODO: VECTORIZE HERE + call give_polynomial_mult_center_one_e_erf_opt(A1_center, A2_center, power_A1, power_A2, C_center(ipoint,1:3), n_pt_in, d, n_pt_out, p_inv_2, p_new2, P_center) + + if(n_pt_out < 0) then + cycle + endif + + ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i + accu = 0.d0 + do i = 0, n_pt_out, 2 + accu += d(i) * rint(i/2, const) + enddo + + res_v(ipoint) = accu * coeff + enddo + + endif + +end subroutine NAI_pol_mult_erf_with1s_v + +! --- + +subroutine give_polynomial_mult_center_one_e_erf_opt(A_center, B_center, power_A, power_B, C_center, n_pt_in, d, n_pt_out, p_inv_2, p_new, P_center) + + BEGIN_DOC + ! Returns the explicit polynomial in terms of the $t$ variable of the + ! following polynomial: + ! + ! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$. + END_DOC + + implicit none + integer, intent(in) :: n_pt_in + integer, intent(in) :: power_A(3), power_B(3) + double precision, intent(in) :: A_center(3), B_center(3), C_center(3), p_inv_2, p_new, P_center(3) + integer, intent(out) :: n_pt_out + double precision, intent(out) :: d(0:n_pt_in) + + integer :: a_x, b_x, a_y, b_y, a_z, b_z + integer :: n_pt1, n_pt2, n_pt3, dim, i + integer :: n_pt_tmp + double precision :: d1(0:n_pt_in) + double precision :: d2(0:n_pt_in) + double precision :: d3(0:n_pt_in) + double precision :: accu + double precision :: R1x(0:2), B01(0:2), R1xp(0:2), R2x(0:2) + + accu = 0.d0 + ASSERT (n_pt_in > 1) + + R1x(0) = (P_center(1) - A_center(1)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(1) - C_center(1))* p_new + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(1) - B_center(1)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(1) - C_center(1))* p_new + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R2x(0) = p_inv_2 + R2x(1) = 0.d0 + R2x(2) = -p_inv_2 * p_new + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + + do i = 0, n_pt_in + d (i) = 0.d0 + d1(i) = 0.d0 + d2(i) = 0.d0 + d3(i) = 0.d0 + enddo + + n_pt1 = n_pt_in + n_pt2 = n_pt_in + n_pt3 = n_pt_in + a_x = power_A(1) + b_x = power_B(1) + call I_x1_pol_mult_one_e(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in) + if(n_pt1<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + R1x(0) = (P_center(2) - A_center(2)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(2) - C_center(2))* p_new + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(2) - B_center(2)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(2) - C_center(2))* p_new + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + a_y = power_A(2) + b_y = power_B(2) + call I_x1_pol_mult_one_e(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in) + if(n_pt2<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + R1x(0) = (P_center(3) - A_center(3)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(3) - C_center(3)) * p_new + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(3) - B_center(3)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(3) - C_center(3)) * p_new + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + a_z = power_A(3) + b_z = power_B(3) + + call I_x1_pol_mult_one_e(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in) + if(n_pt3 < 0) then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + n_pt_tmp = 0 + call multiply_poly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp) + do i = 0, n_pt_tmp + d1(i) = 0.d0 + enddo + n_pt_out = 0 + call multiply_poly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out) + do i = 0, n_pt_out + d(i) = d1(i) + enddo + +end subroutine give_polynomial_mult_center_one_e_erf_opt + +! --- + +subroutine give_polynomial_mult_center_one_e_erf(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out,mu_in) + + BEGIN_DOC + ! Returns the explicit polynomial in terms of the $t$ variable of the + ! following polynomial: + ! + ! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$. + END_DOC + + implicit none + integer, intent(in) :: n_pt_in + integer,intent(out) :: n_pt_out + double precision, intent(in) :: A_center(3), B_center(3),C_center(3) + double precision, intent(in) :: alpha,beta,mu_in + integer, intent(in) :: power_A(3), power_B(3) + integer :: a_x,b_x,a_y,b_y,a_z,b_z + double precision :: d(0:n_pt_in) + double precision :: d1(0:n_pt_in) + double precision :: d2(0:n_pt_in) + double precision :: d3(0:n_pt_in) + double precision :: accu, pq_inv, p10_1, p10_2, p01_1, p01_2 + double precision :: p,P_center(3),rho,p_inv,p_inv_2 + accu = 0.d0 + ASSERT (n_pt_in > 1) + p = alpha+beta + p_inv = 1.d0/p + p_inv_2 = 0.5d0/p + do i =1, 3 + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + enddo + + double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2) + R1x(0) = (P_center(1) - A_center(1)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in) + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(1) - B_center(1)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(1) - C_center(1))* mu_in**2 / (p+mu_in*mu_in) + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R2x(0) = p_inv_2 + R2x(1) = 0.d0 + R2x(2) = -p_inv_2* mu_in**2 / (p+mu_in*mu_in) + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + do i = 0,n_pt_in + d1(i) = 0.d0 + enddo + do i = 0,n_pt_in + d2(i) = 0.d0 + enddo + do i = 0,n_pt_in + d3(i) = 0.d0 + enddo + integer :: n_pt1,n_pt2,n_pt3,dim,i + n_pt1 = n_pt_in + n_pt2 = n_pt_in + n_pt3 = n_pt_in + a_x = power_A(1) + b_x = power_B(1) + call I_x1_pol_mult_one_e(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in) + ! print*,'passed the first I_x1' + if(n_pt1<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + R1x(0) = (P_center(2) - A_center(2)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in) + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(2) - B_center(2)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(2) - C_center(2))* mu_in**2 / (p+mu_in*mu_in) + !R1xp = (P_x - B_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + a_y = power_A(2) + b_y = power_B(2) + call I_x1_pol_mult_one_e(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in) + ! print*,'passed the second I_x1' + if(n_pt2<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + + R1x(0) = (P_center(3) - A_center(3)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in) + ! R1x = (P_x - A_x) - (P_x - C_x) ( t * mu/sqrt(p+mu^2) )^2 + R1xp(0) = (P_center(3) - B_center(3)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(3) - C_center(3))* mu_in**2 / (p+mu_in*mu_in) + !R2x = 0.5 / p - 0.5/p ( t * mu/sqrt(p+mu^2) )^2 + a_z = power_A(3) + b_z = power_B(3) + + ! print*,'a_z = ',a_z + ! print*,'b_z = ',b_z + call I_x1_pol_mult_one_e(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in) + ! print*,'passed the third I_x1' + if(n_pt3<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + integer :: n_pt_tmp + n_pt_tmp = 0 + call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp) + do i = 0,n_pt_tmp + d1(i) = 0.d0 + enddo + n_pt_out = 0 + call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out) + do i = 0, n_pt_out + d(i) = d1(i) + enddo + +end + diff --git a/src/ao_cart_one_e_ints/pot_ao_ints.irp.f b/src/ao_cart_one_e_ints/pot_ao_ints.irp.f new file mode 100644 index 00000000..4c59afb7 --- /dev/null +++ b/src/ao_cart_one_e_ints/pot_ao_ints.irp.f @@ -0,0 +1,610 @@ + +! --- + +BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e, (ao_cart_num,ao_cart_num)] + + BEGIN_DOC + ! Nucleus-electron interaction, in the |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` + ! + ! These integrals also contain the pseudopotential integrals. + END_DOC + + implicit none + integer :: num_A, num_B, power_A(3), power_B(3) + integer :: i, j, k, l, n_pt_in, m + double precision :: alpha, beta + double precision :: A_center(3),B_center(3),C_center(3) + double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + + ao_cart_integrals_n_e = 0.d0 + + if (read_ao_cart_integrals_n_e) then + + call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_n_e(ao_cart_integrals_n_e) + print *, 'AO N-e integrals read from disk' + + else + + if(use_cgtos) then + + do j = 1, ao_cart_num + do i = 1, ao_cart_num + ao_cart_integrals_n_e(i,j) = ao_cart_integrals_n_e_cgtos(i,j) + enddo + enddo + + else + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c,c1,n_pt_in) & + !$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,& + !$OMP n_pt_max_integrals,ao_cart_integrals_n_e,nucl_num,nucl_charge) + + n_pt_in = n_pt_max_integrals + + !$OMP DO SCHEDULE (dynamic) + + do j = 1, ao_cart_num + num_A = ao_cart_nucl(j) + power_A(1:3)= ao_cart_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_cart_num + + num_B = ao_cart_nucl(i) + power_B(1:3)= ao_cart_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(l,j) + + do m=1,ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(m,i) + + double precision :: c, c1 + c = 0.d0 + + do k = 1, nucl_num + double precision :: Z + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + !print *, ' ' + !print *, A_center, B_center, C_center, power_A, power_B + !print *, alpha, beta + + c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & + , alpha, beta, C_center, n_pt_in ) + + !print *, ' c1 = ', c1 + + c = c - Z * c1 + + enddo + ao_cart_integrals_n_e(i,j) = ao_cart_integrals_n_e(i,j) & + + ao_cart_coef_normalized_ordered_transp(l,j) & + * ao_cart_coef_normalized_ordered_transp(m,i) * c + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + + endif + + + IF(do_pseudo) THEN + ao_cart_integrals_n_e += ao_cart_pseudo_integrals + ENDIF + IF(point_charges) THEN + ao_cart_integrals_n_e += ao_cart_integrals_pt_chrg + ENDIF + + endif + + + if (write_ao_cart_integrals_n_e) then + call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_n_e(ao_cart_integrals_n_e) + print *, 'AO N-e integrals written to disk' + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_imag, (ao_cart_num,ao_cart_num)] + BEGIN_DOC + ! Nucleus-electron interaction, in the |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` + END_DOC + implicit none + double precision :: alpha, beta + integer :: num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt_in,m + double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + + if (read_ao_cart_integrals_n_e) then + call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_n_e_imag(ao_cart_integrals_n_e_imag) + print *, 'AO N-e integrals read from disk' + else + print *, irp_here, ': Not yet implemented' + endif +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, ao_cart_integrals_n_e_per_atom, (ao_cart_num,ao_cart_num,nucl_num)] + BEGIN_DOC +! Nucleus-electron interaction in the |AO| basis set, per atom A. +! +! :math:`\langle \chi_i | -\frac{1}{|r-R_A|} | \chi_j \rangle` + END_DOC + implicit none + double precision :: alpha, beta + integer :: i_c,num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt_in,m + double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + + ao_cart_integrals_n_e_per_atom = 0.d0 + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,power_A,power_B,& + !$OMP num_A,num_B,c,n_pt_in,C_center) & + !$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,& + !$OMP n_pt_max_integrals,ao_cart_integrals_n_e_per_atom,nucl_num) + n_pt_in = n_pt_max_integrals + !$OMP DO SCHEDULE (dynamic) + + double precision :: c + do j = 1, ao_cart_num + power_A(1)= ao_cart_power(j,1) + power_A(2)= ao_cart_power(j,2) + power_A(3)= ao_cart_power(j,3) + num_A = ao_cart_nucl(j) + A_center(1) = nucl_coord(num_A,1) + A_center(2) = nucl_coord(num_A,2) + A_center(3) = nucl_coord(num_A,3) + do k = 1, nucl_num + C_center(1) = nucl_coord(k,1) + C_center(2) = nucl_coord(k,2) + C_center(3) = nucl_coord(k,3) + do i = 1, ao_cart_num + power_B(1)= ao_cart_power(i,1) + power_B(2)= ao_cart_power(i,2) + power_B(3)= ao_cart_power(i,3) + num_B = ao_cart_nucl(i) + B_center(1) = nucl_coord(num_B,1) + B_center(2) = nucl_coord(num_B,2) + B_center(3) = nucl_coord(num_B,3) + c = 0.d0 + do l=1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(l,j) + do m=1,ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(m,i) + c = c + NAI_pol_mult(A_center,B_center,power_A,power_B, & + alpha,beta,C_center,n_pt_in) & + * ao_cart_coef_normalized_ordered_transp(l,j) & + * ao_cart_coef_normalized_ordered_transp(m,i) + enddo + enddo + ao_cart_integrals_n_e_per_atom(i,j,k) = -c + enddo + enddo + enddo + !$OMP END DO + !$OMP END PARALLEL + +END_PROVIDER + + + +double precision function NAI_pol_mult(A_center,B_center,power_A,power_B,alpha,beta,C_center,n_pt_in) + BEGIN_DOC +! Computes the electron-nucleus attraction with two primitves. +! +! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle` + END_DOC + + implicit none + integer, intent(in) :: n_pt_in + double precision,intent(in) :: C_center(3),A_center(3),B_center(3),alpha,beta + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,n_pt + double precision :: P_center(3) + double precision :: d(0:n_pt_in),pouet,coeff,rho,dist,const,pouet_2,p,p_inv,factor + double precision :: I_n_special_exact,integrate_bourrin,I_n_bibi + double precision :: V_n_e,const_factor,dist_integral,tmp + double precision :: accu,epsilo,rint + integer :: n_pt_out,lmax + include 'utils/constants.include.F' + if ( (A_center(1)/=B_center(1)).or. & + (A_center(2)/=B_center(2)).or. & + (A_center(3)/=B_center(3)).or. & + (A_center(1)/=C_center(1)).or. & + (A_center(2)/=C_center(2)).or. & + (A_center(3)/=C_center(3))) then + continue + else + NAI_pol_mult = V_n_e(power_A(1),power_A(2),power_A(3), & + power_B(1),power_B(2),power_B(3),alpha,beta) + return + endif + p = alpha + beta + p_inv = 1.d0/p + rho = alpha * beta * p_inv + dist = 0.d0 + dist_integral = 0.d0 + do i = 1, 3 + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + dist += (A_center(i) - B_center(i))*(A_center(i) - B_center(i)) + dist_integral += (P_center(i) - C_center(i))*(P_center(i) - C_center(i)) + enddo + const_factor = dist*rho + const = p * dist_integral + if(const_factor > 80.d0)then + NAI_pol_mult = 0.d0 + return + endif + factor = dexp(-const_factor) + coeff = dtwo_pi * factor * p_inv + lmax = 20 + + ! print*, "b" + do i = 0, n_pt_in + d(i) = 0.d0 + enddo + n_pt = 2 * ( (power_A(1) + power_B(1)) +(power_A(2) + power_B(2)) +(power_A(3) + power_B(3)) ) + if (n_pt == 0) then + epsilo = 1.d0 + pouet = rint(0,const) + NAI_pol_mult = coeff * pouet + return + endif + + call give_polynomial_mult_center_one_e(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out) + + + if(n_pt_out<0)then + NAI_pol_mult = 0.d0 + return + endif + accu = 0.d0 + + ! 1/r1 standard attraction integral + epsilo = 1.d0 + ! sum of integrals of type : int {t,[0,1]} exp-(rho.(P-Q)^2 * t^2) * t^i + do i =0 ,n_pt_out,2 + accu += d(i) * rint(i/2,const) + enddo + NAI_pol_mult = accu * coeff + +end + +! --- + +subroutine give_polynomial_mult_center_one_e(A_center,B_center,alpha,beta,power_A,power_B,C_center,n_pt_in,d,n_pt_out) + implicit none + BEGIN_DOC + ! Returns the explicit polynomial in terms of the "t" variable of the following + ! + ! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$. + END_DOC + integer, intent(in) :: n_pt_in + integer,intent(out) :: n_pt_out + double precision, intent(in) :: A_center(3), B_center(3),C_center(3) + double precision, intent(in) :: alpha,beta + integer, intent(in) :: power_A(3), power_B(3) + integer :: a_x,b_x,a_y,b_y,a_z,b_z + double precision :: d(0:n_pt_in) + double precision :: d1(0:n_pt_in) + double precision :: d2(0:n_pt_in) + double precision :: d3(0:n_pt_in) + double precision :: accu, pq_inv, p10_1, p10_2, p01_1, p01_2 + double precision :: p,P_center(3),rho,p_inv,p_inv_2 + + accu = 0.d0 + + ASSERT (n_pt_in > 1) + p = alpha+beta + p_inv = 1.d0/p + p_inv_2 = 0.5d0/p + do i =1, 3 + P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv + enddo + + double precision :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2) + R1x(0) = (P_center(1) - A_center(1)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(1) - C_center(1)) + + R1xp(0) = (P_center(1) - B_center(1)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(1) - C_center(1)) + + R2x(0) = p_inv_2 + R2x(1) = 0.d0 + R2x(2) = -p_inv_2 + + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + do i = 0,n_pt_in + d1(i) = 0.d0 + enddo + do i = 0,n_pt_in + d2(i) = 0.d0 + enddo + do i = 0,n_pt_in + d3(i) = 0.d0 + enddo + integer :: n_pt1,n_pt2,n_pt3,dim,i + n_pt1 = n_pt_in + n_pt2 = n_pt_in + n_pt3 = n_pt_in + a_x = power_A(1) + b_x = power_B(1) + call I_x1_pol_mult_one_e(a_x,b_x,R1x,R1xp,R2x,d1,n_pt1,n_pt_in) + + if(n_pt1<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + R1x(0) = (P_center(2) - A_center(2)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(2) - C_center(2)) + + R1xp(0) = (P_center(2) - B_center(2)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(2) - C_center(2)) + + a_y = power_A(2) + b_y = power_B(2) + call I_x1_pol_mult_one_e(a_y,b_y,R1x,R1xp,R2x,d2,n_pt2,n_pt_in) + + if(n_pt2<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + + + R1x(0) = (P_center(3) - A_center(3)) + R1x(1) = 0.d0 + R1x(2) = -(P_center(3) - C_center(3)) + + R1xp(0) = (P_center(3) - B_center(3)) + R1xp(1) = 0.d0 + R1xp(2) =-(P_center(3) - C_center(3)) + + a_z = power_A(3) + b_z = power_B(3) + + call I_x1_pol_mult_one_e(a_z,b_z,R1x,R1xp,R2x,d3,n_pt3,n_pt_in) + + if(n_pt3<0)then + n_pt_out = -1 + do i = 0,n_pt_in + d(i) = 0.d0 + enddo + return + endif + integer :: n_pt_tmp + n_pt_tmp = 0 + call multiply_poly(d1,n_pt1,d2,n_pt2,d,n_pt_tmp) + do i = 0,n_pt_tmp + d1(i) = 0.d0 + enddo + n_pt_out = 0 + call multiply_poly(d ,n_pt_tmp ,d3,n_pt3,d1,n_pt_out) + do i = 0, n_pt_out + d(i) = d1(i) + enddo + +end + + +recursive subroutine I_x1_pol_mult_one_e(a,c,R1x,R1xp,R2x,d,nd,n_pt_in) + implicit none + BEGIN_DOC +! Recursive routine involved in the electron-nucleus potential + END_DOC + integer , intent(in) :: n_pt_in + double precision,intent(inout) :: d(0:n_pt_in) + integer,intent(inout) :: nd + integer, intent(in) :: a,c + double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2) + include 'utils/constants.include.F' + double precision :: X(0:max_dim) + double precision :: Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y + integer :: nx, ix,dim,iy,ny + dim = n_pt_in + ! print*,'a,c = ',a,c + ! print*,'nd_in = ',nd + + if( (a==0) .and. (c==0))then + nd = 0 + d(0) = 1.d0 + return + elseif( (c<0).or.(nd<0) )then + nd = -1 + return + else if ((a==0).and.(c.ne.0)) then + call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,n_pt_in) + else if (a==1) then + nx = nd + do ix=0,n_pt_in + X(ix) = 0.d0 + Y(ix) = 0.d0 + enddo + call I_x2_pol_mult_one_e(c-1,R1x,R1xp,R2x,X,nx,n_pt_in) + do ix=0,nx + X(ix) *= dble(c) + enddo +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) + ny=0 + call I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,Y,ny,n_pt_in) +! call multiply_poly(Y,ny,R1x,2,d,nd) + call multiply_poly_c2(Y,ny,R1x,d,nd) + else + do ix=0,n_pt_in + X(ix) = 0.d0 + Y(ix) = 0.d0 + enddo + nx = 0 + call I_x1_pol_mult_one_e(a-2,c,R1x,R1xp,R2x,X,nx,n_pt_in) + do ix=0,nx + X(ix) *= dble(a-1) + enddo +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) + + nx = nd + do ix=0,n_pt_in + X(ix) = 0.d0 + enddo + call I_x1_pol_mult_one_e(a-1,c-1,R1x,R1xp,R2x,X,nx,n_pt_in) + do ix=0,nx + X(ix) *= dble(c) + enddo +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) + ny=0 + call I_x1_pol_mult_one_e(a-1,c,R1x,R1xp,R2x,Y,ny,n_pt_in) +! call multiply_poly(Y,ny,R1x,2,d,nd) + call multiply_poly_c2(Y,ny,R1x,d,nd) + endif +end + +recursive subroutine I_x2_pol_mult_one_e(c,R1x,R1xp,R2x,d,nd,dim) + implicit none + BEGIN_DOC +! Recursive routine involved in the electron-nucleus potential + END_DOC + integer , intent(in) :: dim + include 'utils/constants.include.F' + double precision :: d(0:max_dim) + integer,intent(inout) :: nd + integer, intent(in) :: c + double precision, intent(in) :: R1x(0:2),R1xp(0:2),R2x(0:2) + integer :: i + + if(c==0) then + nd = 0 + d(0) = 1.d0 + return + elseif ((nd<0).or.(c<0))then + nd = -1 + return + else + integer :: nx, ix,ny + double precision :: X(0:max_dim),Y(0:max_dim) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y + do ix=0,dim + X(ix) = 0.d0 + Y(ix) = 0.d0 + enddo + nx = 0 + call I_x1_pol_mult_one_e(0,c-2,R1x,R1xp,R2x,X,nx,dim) + do ix=0,nx + X(ix) *= dble(c-1) + enddo +! call multiply_poly(X,nx,R2x,2,d,nd) + call multiply_poly_c2(X,nx,R2x,d,nd) + ny = 0 + do ix=0,dim + Y(ix) = 0.d0 + enddo + + call I_x1_pol_mult_one_e(0,c-1,R1x,R1xp,R2x,Y,ny,dim) + if(ny.ge.0)then +! call multiply_poly(Y,ny,R1xp,2,d,nd) + call multiply_poly_c2(Y,ny,R1xp,d,nd) + endif + endif +end + +double precision function V_n_e(a_x,a_y,a_z,b_x,b_y,b_z,alpha,beta) + implicit none + BEGIN_DOC +! Primitve nuclear attraction between the two primitves centered on the same atom. +! +! $p_1 = x^{a_x} y^{a_y} z^{a_z} \exp(-\alpha r^2)$ +! +! $p_2 = x^{b_x} y^{b_y} z^{b_z} \exp(-\beta r^2)$ + END_DOC + integer :: a_x,a_y,a_z,b_x,b_y,b_z + double precision :: alpha,beta + double precision :: V_r, V_phi, V_theta + if(iand((a_x+b_x),1)==1.or.iand(a_y+b_y,1)==1.or.iand((a_z+b_z),1)==1)then + V_n_e = 0.d0 + else + V_n_e = V_r(a_x+b_x+a_y+b_y+a_z+b_z+1,alpha+beta) & + * V_phi(a_x+b_x,a_y+b_y) & + * V_theta(a_z+b_z,a_x+b_x+a_y+b_y+1) + endif + +end + + +double precision function int_gaus_pol(alpha,n) + implicit none + BEGIN_DOC +! Computes the integral: +! +! $\int_{-\infty}^{\infty} x^n \exp(-\alpha x^2) dx$. + END_DOC + double precision :: alpha + integer :: n + double precision :: dble_fact + include 'utils/constants.include.F' + + int_gaus_pol = 0.d0 + if(iand(n,1).eq.0)then + int_gaus_pol = dsqrt(alpha/pi) + double precision :: two_alpha + two_alpha = alpha+alpha + integer :: i + do i=1,n,2 + int_gaus_pol = int_gaus_pol * two_alpha + enddo + int_gaus_pol = dble_fact(n -1) / int_gaus_pol + endif + +end + +double precision function V_r(n,alpha) + implicit none + BEGIN_DOC + ! Computes the radial part of the nuclear attraction integral: + ! + ! $\int_{0}^{\infty} r^n \exp(-\alpha r^2) dr$ + ! + END_DOC + double precision :: alpha, fact + integer :: n + include 'utils/constants.include.F' + if(iand(n,1).eq.1)then + V_r = 0.5d0 * fact(shiftr(n,1)) / (alpha ** (shiftr(n,1) + 1)) + else + V_r = sqpi * fact(n) / fact(shiftr(n,1)) * (0.5d0/sqrt(alpha)) ** (n+1) + endif +end + diff --git a/src/ao_cart_one_e_ints/pot_ao_pseudo_ints.irp.f b/src/ao_cart_one_e_ints/pot_ao_pseudo_ints.irp.f new file mode 100644 index 00000000..b075faa9 --- /dev/null +++ b/src/ao_cart_one_e_ints/pot_ao_pseudo_ints.irp.f @@ -0,0 +1,311 @@ +BEGIN_PROVIDER [ double precision, ao_cart_pseudo_integrals, (ao_cart_num,ao_cart_num)] + implicit none + BEGIN_DOC + ! Pseudo-potential integrals in the |AO| basis set. + END_DOC + + if (read_ao_cart_integrals_pseudo) then + call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_pseudo(ao_cart_pseudo_integrals) + print *, 'AO pseudopotential integrals read from disk' + else + + ao_cart_pseudo_integrals = 0.d0 + if (do_pseudo) then + if (pseudo_klocmax > 0) then + ao_cart_pseudo_integrals += ao_cart_pseudo_integrals_local + endif + if (pseudo_kmax > 0) then + ao_cart_pseudo_integrals += ao_cart_pseudo_integrals_non_local + endif + endif + endif + + if (write_ao_cart_integrals_pseudo) then + call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_pseudo(ao_cart_pseudo_integrals) + print *, 'AO pseudopotential integrals written to disk' + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, ao_cart_pseudo_integrals_local, (ao_cart_num,ao_cart_num)] + use omp_lib + implicit none + BEGIN_DOC + ! Local pseudo-potential + END_DOC + include 'utils/constants.include.F' + double precision :: alpha, beta, gama, delta + integer :: num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,m + double precision :: Vloc, Vpseudo + + double precision :: wall_1, wall_2, wall_0 + integer :: thread_num + double precision :: c + double precision :: Z + + PROVIDE ao_cart_coef_normalized_ordered_transp + PROVIDE pseudo_v_k_transp pseudo_n_k_transp pseudo_klocmax pseudo_dz_k_transp + + ao_cart_pseudo_integrals_local = 0.d0 + + print*, 'Providing the nuclear electron pseudo integrals (local)' + + ! Dummy iteration for OpenMP + j=1 + i=1 + l=1 + m=1 + num_A = ao_cart_nucl(j) + power_A(1:3)= ao_cart_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + num_B = ao_cart_nucl(i) + power_B(1:3)= ao_cart_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + alpha = ao_cart_expo_ordered_transp(l,j) + beta = ao_cart_expo_ordered_transp(m,i) + c = 0.d0 + do k = 1, nucl_num + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + c = c + Vloc(pseudo_klocmax, & + pseudo_v_k_transp (1,k), & + pseudo_n_k_transp (1,k), & + pseudo_dz_k_transp(1,k), & + A_center,power_A,alpha,B_center,power_B,beta,C_center) + + enddo + + + ao_cart_pseudo_integrals_local = 0.d0 + call wall_time(wall_1) + + thread_num = 0 + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c, & + !$OMP wall_0,wall_2,thread_num) & + !$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,& + !$OMP ao_cart_pseudo_integrals_local,nucl_num,nucl_charge, & + !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_v_k_transp,pseudo_n_k_transp, pseudo_dz_k_transp,& + !$OMP wall_1) + + !$ thread_num = omp_get_thread_num() + + wall_0 = wall_1 + !$OMP DO + + do j = 1, ao_cart_num + + num_A = ao_cart_nucl(j) + power_A(1:3)= ao_cart_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_cart_num + + num_B = ao_cart_nucl(i) + power_B(1:3)= ao_cart_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(l,j) + + do m=1,ao_cart_prim_num(i) + if (dabs(ao_cart_coef_normalized_ordered_transp(l,j)*ao_cart_coef_normalized_ordered_transp(m,i))& + < thresh) then + cycle + endif + + beta = ao_cart_expo_ordered_transp(m,i) + c = 0.d0 + + do k = 1, nucl_num + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + c = c + Vloc(pseudo_klocmax, & + pseudo_v_k_transp (1,k), & + pseudo_n_k_transp (1,k), & + pseudo_dz_k_transp(1,k), & + A_center,power_A,alpha,B_center,power_B,beta,C_center) + + enddo + ao_cart_pseudo_integrals_local(i,j) = ao_cart_pseudo_integrals_local(i,j) +& + ao_cart_coef_normalized_ordered_transp(l,j)*ao_cart_coef_normalized_ordered_transp(m,i)*c + enddo + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(j)/float(ao_cart_num), '% in ', & + wall_2-wall_1, 's' + endif + endif + enddo + + !$OMP END DO + !$OMP END PARALLEL + + do i=1,ao_cart_num + do j=1,i + ao_cart_pseudo_integrals_local(j,i) = 0.5d0*(ao_cart_pseudo_integrals_local(i,j) + ao_cart_pseudo_integrals_local(i,j)) + ao_cart_pseudo_integrals_local(i,j) = ao_cart_pseudo_integrals_local(i,j) + enddo + enddo + END_PROVIDER + + + BEGIN_PROVIDER [ double precision, ao_cart_pseudo_integrals_non_local, (ao_cart_num,ao_cart_num)] + use omp_lib + implicit none + BEGIN_DOC + ! Non-local pseudo-potential + END_DOC + include 'utils/constants.include.F' + double precision :: alpha, beta, gama, delta + integer :: num_A,num_B + double precision :: A_center(3),B_center(3),C_center(3) + integer :: power_A(3),power_B(3) + integer :: i,j,k,l,m + double precision :: Vloc, Vpseudo + + double precision :: wall_1, wall_2, wall_0 + integer :: thread_num + double precision :: c + double precision :: Z + + PROVIDE ao_cart_coef_normalized_ordered_transp + PROVIDE pseudo_lmax pseudo_kmax pseudo_v_kl_transp pseudo_n_kl_transp pseudo_dz_kl_transp + ao_cart_pseudo_integrals_non_local = 0.d0 + + print*, 'Providing the nuclear electron pseudo integrals (non-local)' + + call wall_time(wall_1) + thread_num = 0 + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c, & + !$OMP wall_0,wall_2,thread_num) & + !$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,nucl_coord,ao_cart_coef_normalized_ordered_transp,& + !$OMP ao_cart_pseudo_integrals_non_local,nucl_num,nucl_charge,& + !$OMP pseudo_klocmax,pseudo_lmax,pseudo_kmax,pseudo_n_kl_transp, pseudo_v_kl_transp, pseudo_dz_kl_transp,& + !$OMP wall_1) + + !$ thread_num = omp_get_thread_num() + + wall_0 = wall_1 + !$OMP DO SCHEDULE (guided) +! + do j = 1, ao_cart_num + + num_A = ao_cart_nucl(j) + power_A(1:3)= ao_cart_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_cart_num + + num_B = ao_cart_nucl(i) + power_B(1:3)= ao_cart_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(l,j) + + do m=1,ao_cart_prim_num(i) + if (dabs(ao_cart_coef_normalized_ordered_transp(l,j)*ao_cart_coef_normalized_ordered_transp(m,i))& + < thresh) then + cycle + endif + + beta = ao_cart_expo_ordered_transp(m,i) + c = 0.d0 + + do k = 1, nucl_num + Z = nucl_charge(k) + + C_center(1:3) = nucl_coord(k,1:3) + + c = c + Vpseudo(pseudo_lmax,pseudo_kmax, & + pseudo_v_kl_transp(1,0,k), & + pseudo_n_kl_transp(1,0,k), & + pseudo_dz_kl_transp(1,0,k), & + A_center,power_A,alpha,B_center,power_B,beta,C_center) + enddo + ao_cart_pseudo_integrals_non_local(i,j) = ao_cart_pseudo_integrals_non_local(i,j) +& + ao_cart_coef_normalized_ordered_transp(l,j)*ao_cart_coef_normalized_ordered_transp(m,i)*c + enddo + enddo + enddo + + call wall_time(wall_2) + if (thread_num == 0) then + if (wall_2 - wall_0 > 1.d0) then + wall_0 = wall_2 + print*, 100.*float(j)/float(ao_cart_num), '% in ', & + wall_2-wall_1, 's' + endif + endif + enddo + + !$OMP END DO + + !$OMP END PARALLEL + + + do i=1,ao_cart_num + do j=1,i + ao_cart_pseudo_integrals_non_local(j,i) = 0.5d0*(ao_cart_pseudo_integrals_non_local(i,j) + ao_cart_pseudo_integrals_non_local(i,j)) + ao_cart_pseudo_integrals_non_local(i,j) = ao_cart_pseudo_integrals_non_local(i,j) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, pseudo_v_k_transp, (pseudo_klocmax,nucl_num) ] +&BEGIN_PROVIDER [ integer , pseudo_n_k_transp, (pseudo_klocmax,nucl_num) ] +&BEGIN_PROVIDER [ double precision, pseudo_dz_k_transp, (pseudo_klocmax,nucl_num)] + implicit none + BEGIN_DOC + ! Transposed arrays for pseudopotentials + END_DOC + + integer :: i,j + do j=1,nucl_num + do i=1,pseudo_klocmax + pseudo_v_k_transp (i,j) = pseudo_v_k (j,i) + pseudo_n_k_transp (i,j) = pseudo_n_k (j,i) + pseudo_dz_k_transp(i,j) = pseudo_dz_k(j,i) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ double precision, pseudo_v_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num) ] +&BEGIN_PROVIDER [ integer , pseudo_n_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num) ] +&BEGIN_PROVIDER [ double precision, pseudo_dz_kl_transp, (pseudo_kmax,0:pseudo_lmax,nucl_num)] + implicit none + BEGIN_DOC + ! Transposed arrays for pseudopotentials + END_DOC + + integer :: i,j,l + do j=1,nucl_num + do l=0,pseudo_lmax + do i=1,pseudo_kmax + pseudo_v_kl_transp (i,l,j) = pseudo_v_kl (j,i,l) + pseudo_n_kl_transp (i,l,j) = pseudo_n_kl (j,i,l) + pseudo_dz_kl_transp(i,l,j) = pseudo_dz_kl(j,i,l) + enddo + enddo + enddo +END_PROVIDER + diff --git a/src/ao_cart_one_e_ints/pot_pt_charges.irp.f b/src/ao_cart_one_e_ints/pot_pt_charges.irp.f new file mode 100644 index 00000000..940bbd89 --- /dev/null +++ b/src/ao_cart_one_e_ints/pot_pt_charges.irp.f @@ -0,0 +1,108 @@ + +BEGIN_PROVIDER [ double precision, ao_cart_integrals_pt_chrg, (ao_cart_num,ao_cart_num)] + + BEGIN_DOC + ! Point charge-electron interaction, in the |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_charge charge * \frac{1}{|r-R_charge|} | \chi_j \rangle` + ! + ! Notice the minus sign convention as it is supposed to be for electrons. + END_DOC + + implicit none + integer :: num_A, num_B, power_A(3), power_B(3) + integer :: i, j, k, l, n_pt_in, m + double precision :: alpha, beta + double precision :: A_center(3),B_center(3),C_center(3) + double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + + ao_cart_integrals_pt_chrg = 0.d0 + +! if (read_ao_cart_integrals_pt_chrg) then +! +! call ezfio_get_ao_cart_one_e_ints_ao_cart_integrals_pt_chrg(ao_cart_integrals_pt_chrg) +! print *, 'AO N-e integrals read from disk' +! +! else + +! if(use_cgtos) then +! !print *, " use_cgtos for ao_cart_integrals_pt_chrg ?", use_cgtos +! +! do j = 1, ao_cart_num +! do i = 1, ao_cart_num +! ao_cart_integrals_pt_chrg(i,j) = ao_cart_integrals_pt_chrg_cgtos(i,j) +! enddo +! enddo +! +! else + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c,c1,n_pt_in) & + !$OMP SHARED (ao_cart_num,ao_cart_prim_num,ao_cart_expo_ordered_transp,ao_cart_power,ao_cart_nucl,pts_charge_coord,ao_cart_coef_normalized_ordered_transp,nucl_coord,& + !$OMP n_pt_max_integrals,ao_cart_integrals_pt_chrg,n_pts_charge,pts_charge_z) + + n_pt_in = n_pt_max_integrals + + !$OMP DO SCHEDULE (dynamic) + + do j = 1, ao_cart_num + num_A = ao_cart_nucl(j) + power_A(1:3)= ao_cart_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_cart_num + + num_B = ao_cart_nucl(i) + power_B(1:3)= ao_cart_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(l,j) + + do m=1,ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(m,i) + + double precision :: c, c1 + c = 0.d0 + + do k = 1, n_pts_charge + double precision :: Z + Z = pts_charge_z(k) + + C_center(1:3) = pts_charge_coord(k,1:3) + + c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & + , alpha, beta, C_center, n_pt_in ) + + c = c - Z * c1 + + enddo + ao_cart_integrals_pt_chrg(i,j) = ao_cart_integrals_pt_chrg(i,j) & + + ao_cart_coef_normalized_ordered_transp(l,j) & + * ao_cart_coef_normalized_ordered_transp(m,i) * c + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + +! endif + + +! IF(do_pseudo) THEN +! ao_cart_integrals_pt_chrg += ao_cart_pseudo_integrals +! ENDIF + +! endif + + +! if (write_ao_cart_integrals_pt_chrg) then +! call ezfio_set_ao_cart_one_e_ints_ao_cart_integrals_pt_chrg(ao_cart_integrals_pt_chrg) +! print *, 'AO N-e integrals written to disk' +! endif + +END_PROVIDER diff --git a/src/ao_cart_one_e_ints/pseudopot.f90 b/src/ao_cart_one_e_ints/pseudopot.f90 new file mode 100644 index 00000000..141292d8 --- /dev/null +++ b/src/ao_cart_one_e_ints/pseudopot.f90 @@ -0,0 +1,2127 @@ +!! INFO : You can display equations using : http://www.codecogs.com/latex/eqneditor.php + +!! +!! {\tt Vps}(C) = \langle \Phi_A|{\tt Vloc}(C)+{\tt Vpp}(C)| \Phi_B \rangle +!! +!! with : +!! +!! {\tt Vloc}(C)=\sum_{k=1}^{\tt klocmax} v_k r_C^{n_k} \exp(-dz_k r_C^2) \\ +!! +!! {\tt Vpp}(C)=\sum_{l=0}^{\tt lmax}\left( \sum_{k=1}^{\tt kmax} v_{kl} +!! r_C^{n_{kl}} \exp(-dz_{kl} r_C)^2 \right) |l\rangle \langle l| +!! +double precision function Vps & +(a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) +implicit none +integer n_a(3),n_b(3) +double precision g_a,g_b,a(3),b(3),c(3) +integer lmax,kmax,n_kl(kmax,0:lmax) +double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) +integer klocmax,n_k(klocmax) +double precision v_k(klocmax),dz_k(klocmax) +double precision Vloc,Vpseudo + +Vps=Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) & + +Vpseudo(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +end +!! +!! Vps_num: brute force numerical evaluation of the same matrix element Vps +!! +double precision function Vps_num & +(npts,rmax,a,n_a,g_a,b,n_b,g_b,c,klocmax,v_k,n_k,dz_k,lmax,kmax,v_kl,n_kl,dz_kl) +implicit none +integer n_a(3),n_b(3) +double precision g_a,g_b,a(3),b(3),c(3),rmax +integer lmax,kmax,n_kl(kmax,0:lmax) +double precision v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) +integer klocmax,n_k(klocmax) +double precision v_k(klocmax),dz_k(klocmax) +double precision Vloc_num,Vpseudo_num,v1,v2 +integer npts,nptsgrid +nptsgrid=50 +call initpseudos(nptsgrid) +v1=Vloc_num(npts,rmax,klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) +v2=Vpseudo_num(nptsgrid,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +Vps_num=v1+v2 +end + +double precision function Vloc_num(npts_over,xmax,klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) +implicit none +integer klocmax +double precision v_k(klocmax),dz_k(klocmax) +integer n_k(klocmax) +integer npts_over,ix,iy,iz +double precision xmax,dx,x,y,z +double precision a(3),b(3),c(3),term,r,orb_phi,g_a,g_b,ac(3),bc(3) +integer n_a(3),n_b(3),k,l +do l=1,3 + ac(l)=a(l)-c(l) + bc(l)=b(l)-c(l) +enddo +dx=2.d0*xmax/npts_over +Vloc_num=0.d0 +do ix=1,npts_over +do iy=1,npts_over +do iz=1,npts_over + x=-xmax+dx*ix+dx/2.d0 + y=-xmax+dx*iy+dx/2.d0 + z=-xmax+dx*iz+dx/2.d0 + term=orb_phi(x,y,z,n_a,ac,g_a)*orb_phi(x,y,z,n_b,bc,g_b) + r=dsqrt(x**2+y**2+z**2) + do k=1,klocmax + Vloc_num=Vloc_num+dx**3*v_k(k)*r**n_k(k)*dexp(-dz_k(k)*r**2)*term + enddo +enddo +enddo +enddo +end + +double precision function orb_phi(x,y,z,npower,center,gamma) +implicit none +integer npower(3) +double precision x,y,z,r2,gamma,center(3) +r2=(x-center(1))**2+(y-center(2))**2+(z-center(3))**2 +orb_phi=(x-center(1))**npower(1)*(y-center(2))**npower(2)*(z-center(3))**npower(3) +orb_phi=orb_phi*dexp(-gamma*r2) +end + +!! Real spherical harmonics Ylm + +! factor = ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 +! Y_lm(theta,phi) = +! m > 0 factor* P_l^|m|(cos(theta)) cos (|m| phi) +! m = 0 1/sqrt(2) *factor* P_l^0(cos(theta)) +! m < 0 factor* P_l^|m|(cos(theta)) sin (|m| phi) +! +! x=cos(theta) + + double precision function ylm_real(l,m,x,phi) + implicit none + integer :: MM, iabs_m, m, l + double precision :: pi, fourpi, factor, x, phi, coef + double precision :: xchap, ychap, zchap + double precision, external :: fact + double precision :: PM(0:100,0:100), plm + MM=100 + pi=dacos(-1.d0) + fourpi=4.d0*pi + iabs_m=iabs(m) + if(iabs_m.gt.l)stop 'm must be between -l and l' + factor= dsqrt( ((l+l+1)*fact(l-iabs_m))/(fourpi*fact(l+iabs_m)) ) + if(dabs(x).gt.1.d0)then + print*,'pb. in ylm_no' + print*,'x=',x + stop + endif + call LPMN(MM,l,l,X,PM) + plm=PM(iabs_m,l) + coef=factor*plm + if(m.gt.0)ylm_real=dsqrt(2.d0)*coef*dcos(iabs_m*phi) + if(m.eq.0)ylm_real=coef + if(m.lt.0)ylm_real=dsqrt(2.d0)*coef*dsin(iabs_m*phi) + + if(l.eq.0)ylm_real=dsqrt(1.d0/fourpi) + + xchap=dsqrt(1.d0-x**2)*dcos(phi) + ychap=dsqrt(1.d0-x**2)*dsin(phi) + zchap=x + if(l.eq.1.and.m.eq.1)ylm_real=dsqrt(3.d0/fourpi)*xchap + if(l.eq.1.and.m.eq.0)ylm_real=dsqrt(3.d0/fourpi)*zchap + if(l.eq.1.and.m.eq.-1)ylm_real=dsqrt(3.d0/fourpi)*ychap + + if(l.eq.2.and.m.eq.2)ylm_real=dsqrt(15.d0/16.d0/pi)*(xchap*xchap-ychap*ychap) + if(l.eq.2.and.m.eq.1)ylm_real=dsqrt(15.d0/fourpi)*xchap*zchap + if(l.eq.2.and.m.eq.0)ylm_real=dsqrt(5.d0/16.d0/pi)*(2.d0*zchap*zchap-xchap*xchap-ychap*ychap) + if(l.eq.2.and.m.eq.-1)ylm_real=dsqrt(15.d0/fourpi)*ychap*zchap + if(l.eq.2.and.m.eq.-2)ylm_real=dsqrt(15.d0/fourpi)*xchap*ychap + + if(l.gt.2)stop 'l > 2 not coded!' + + end +! _ +! | | +! __ __ _ __ ___ ___ _ _ __| | ___ +! \ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ +! \ V / | |_) \__ \ __/ |_| | (_| | (_) | +! \_/ | .__/|___/\___|\__,_|\____|\___/ +! | | +! |_| + +!! Routine Vpseudo is based on formumla (66) +!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976): +!! +!! Vpseudo= (4pi)**2* \sum_{l=0}^lmax \sum_{m=-l}^{l} +!! \sum{lambda=0}^{l+nA} \sum_{mu=-lambda}^{lambda} +!! \sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz} +!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3)* Y_{lambda mu}(AC_unit) +!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)* +!! bigI(lambda,mu,l,m,k1,k2,k3) +!! \sum{lambdap=0}^{l+nB} \sum_{mup=-lambdap}^{lambdap} +!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz} +!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p)* Y_{lambdap mup}(BC_unit) +!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)* +!! bigI(lambdap,mup,l,m,k1p,k2p,k3p)* +!! \sum_{k=1}^{kmax} v_kl(k,l)* +!! bigR(lambda,lambdap,k1+k2+k3+k1p+k2p+k3p+n_kl(k,l),g_a,g_b,AC,BC,dz_kl(k,l)) +!! +!! nA=nAx+nAy+nAz +!! nB=nBx+nBy+nBz +!! AC=|A-C| +!! AC_x= A_x-C_x, etc. +!! BC=|B-C| +!! AC_unit= vect(AC)/AC +!! BC_unit= vect(BC)/BC +!! bigI(lambda,mu,l,m,k1,k2,k3)= +!! \int dOmega Y_{lambda mu}(Omega) xchap^k1 ychap^k2 zchap^k3 Y_{l m}(Omega) +!! +!! bigR(lambda,lambdap,N,g_a,g_b,gamm_k,AC,BC) +!! = exp(-g_a* AC**2 -g_b* BC**2) * int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) +!! /int dx x^{ktot} exp(-g_k)*x**2) M_lambda(2 g_k D x) + +double precision function Vpseudo & +(lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +implicit none + +! ___ +! | ._ ._ _|_ +! _|_ | | |_) |_| |_ +! | +double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) +integer, intent(in) :: lmax,kmax,n_kl(kmax,0:lmax) +integer, intent(in) :: n_a(3),n_b(3) +double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) + +! +! | _ _ _. | +! |_ (_) (_ (_| | +! + +double precision :: fourpi,f,prod,prodp,binom_func,accu,bigR,bigI,ylm +double precision :: theta_AC0,phi_AC0,theta_BC0,phi_BC0,ac,bc,big +double precision :: areal,freal,breal,t1,t2,int_prod_bessel +double precision :: arg + +integer :: ntot,ntotA,m,mu,mup,k1,k2,k3,ntotB,k1p,k2p,k3p,lambda,lambdap,ktot +integer :: l,k, nkl_max + +! _ +! |_) o _ _. ._ ._ _. +! |_) | (_| (_| | | (_| \/ +! _| / + +double precision, allocatable :: array_coefs_A(:,:) +double precision, allocatable :: array_coefs_B(:,:) + +double precision, allocatable :: array_R(:,:,:,:,:) +double precision, allocatable :: array_I_A(:,:,:,:,:) +double precision, allocatable :: array_I_B(:,:,:,:,:) + +double precision :: f1, f2, f3 + +if (kmax.eq.1.and.lmax.eq.0.and.v_kl(1,0).eq.0.d0) then + Vpseudo=0.d0 + return +end if + +fourpi=4.d0*dacos(-1.d0) +ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) +bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) +arg= g_a*ac*ac + g_b*bc*bc + +if(arg.gt.-dlog(1.d-20))then + Vpseudo=0.d0 + return +endif + +freal=dexp(-arg) + +areal=2.d0*g_a*ac +breal=2.d0*g_b*bc +ntotA=n_a(1)+n_a(2)+n_a(3) +ntotB=n_b(1)+n_b(2)+n_b(3) +ntot=ntotA+ntotB + +nkl_max=4 + +allocate (array_coefs_A(0:ntot,3)) +allocate (array_coefs_B(0:ntot,3)) + +allocate (array_R(kmax,0:ntot+nkl_max,0:lmax,0:lmax+ntot,0:lmax+ntot)) + +allocate (array_I_A(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot)) + +allocate (array_I_B(-(lmax+ntot):lmax+ntot,0:lmax+ntot,0:ntot,0:ntot,0:ntot)) + +if(ac.eq.0.d0.and.bc.eq.0.d0)then + + + accu=0.d0 + + do k=1,kmax + do l=0,lmax + ktot=ntot+n_kl(k,l) + if (v_kl(k,l) == 0.d0) cycle + do m=-l,l + prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) + if (prod == 0.d0) cycle + prodp=bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) + if (prodp == 0.d0) cycle + accu=accu+prod*prodp*v_kl(k,l)*int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,0,areal,breal,arg) + enddo + enddo + enddo + + Vpseudo=accu*fourpi + +else if(ac.ne.0.d0.and.bc.ne.0.d0)then + + f=fourpi*fourpi + + theta_AC0=dacos( (a(3)-c(3))/ac ) + phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) + theta_BC0=dacos( (b(3)-c(3))/bc ) + phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc) + + + do lambdap=0,lmax+ntotB + do lambda=0,lmax+ntotA + do l=0,lmax + do ktot=0,ntotA+ntotB+nkl_max + do k=1,kmax + array_R(k,ktot,l,lambda,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,lambdap,areal,breal,arg) + enddo + enddo + enddo + enddo + enddo + + do k1=0,n_a(1) + array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1) + enddo + do k2=0,n_a(2) + array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2) + enddo + do k3=0,n_a(3) + array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3) + enddo + + do k1p=0,n_b(1) + array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p) + enddo + do k2p=0,n_b(2) + array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p) + enddo + do k3p=0,n_b(3) + array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p) + enddo + + accu=0.d0 + do l=0,lmax + do m=-l,l + + do k3=0,n_a(3) + do k2=0,n_a(2) + do k1=0,n_a(1) + do lambda=0,l+ntotA + do mu=-lambda,lambda + array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) + enddo + enddo + enddo + enddo + enddo + + do k3p=0,n_b(3) + do k2p=0,n_b(2) + do k1p=0,n_b(1) + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) + enddo + enddo + enddo + enddo + enddo + + do k3=0,n_a(3) + if (array_coefs_A(k3,3) == 0.d0) cycle + do k2=0,n_a(2) + if (array_coefs_A(k2,2) == 0.d0) cycle + do k1=0,n_a(1) + if (array_coefs_A(k1,1) == 0.d0) cycle + + do lambda=0,l+ntotA + do mu=-lambda,lambda + + prod=ylm(lambda,mu,theta_AC0,phi_AC0)*array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*array_I_A(mu,lambda,k1,k2,k3) + if (prod == 0.d0) cycle + + do k3p=0,n_b(3) + do k2p=0,n_b(2) + do k1p=0,n_b(1) + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + + prodp=prod*ylm(lambdap,mup,theta_BC0,phi_BC0)* & + array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)* & + array_I_B(mup,lambdap,k1p,k2p,k3p) + + if (prodp == 0.d0) cycle + do k=1,kmax + ktot=k1+k2+k3+k1p+k2p+k3p+n_kl(k,l) + accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,lambdap) + enddo + + enddo + enddo + enddo + + enddo + enddo + + enddo + enddo + enddo + + enddo + enddo + + enddo + enddo + + Vpseudo=f*accu + +else if(ac.eq.0.d0.and.bc.ne.0.d0)then + + f=fourpi**1.5d0 + theta_BC0=dacos( (b(3)-c(3))/bc ) + phi_BC0=datan2((b(2)-c(2))/bc,(b(1)-c(1))/bc) + + areal=2.d0*g_a*ac + breal=2.d0*g_b*bc + freal=dexp(-g_a*ac**2-g_b*bc**2) + + do lambdap=0,lmax+ntotB + do l=0,lmax + do ktot=0,ntotA+ntotB+nkl_max + do k=1,kmax + array_R(k,ktot,l,0,lambdap)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),0,lambdap,areal,breal,arg) + enddo + enddo + enddo + enddo + + do k1p=0,n_b(1) + array_coefs_B(k1p,1) = binom_func(n_b(1),k1p)*(c(1)-b(1))**(n_b(1)-k1p) + enddo + do k2p=0,n_b(2) + array_coefs_B(k2p,2) = binom_func(n_b(2),k2p)*(c(2)-b(2))**(n_b(2)-k2p) + enddo + do k3p=0,n_b(3) + array_coefs_B(k3p,3) = binom_func(n_b(3),k3p)*(c(3)-b(3))**(n_b(3)-k3p) + enddo + + accu=0.d0 + do l=0,lmax + do m=-l,l + + do k3p=0,n_b(3) + do k2p=0,n_b(2) + do k1p=0,n_b(1) + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + array_I_B(mup,lambdap,k1p,k2p,k3p)=bigI(lambdap,mup,l,m,k1p,k2p,k3p) + enddo + enddo + enddo + enddo + enddo + + prod=bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) + + do k3p=0,n_b(3) + if (array_coefs_B(k3p,3) == 0.d0) cycle + do k2p=0,n_b(2) + if (array_coefs_B(k2p,2) == 0.d0) cycle + do k1p=0,n_b(1) + if (array_coefs_B(k1p,1) == 0.d0) cycle + do lambdap=0,l+ntotB + do mup=-lambdap,lambdap + + prodp=prod*array_coefs_B(k1p,1)*array_coefs_B(k2p,2)*array_coefs_B(k3p,3)*ylm(lambdap,mup,theta_BC0,phi_BC0)*array_I_B(mup,lambdap,k1p,k2p,k3p) + + if (prodp == 0.d0) cycle + + do k=1,kmax + + ktot=ntotA+k1p+k2p+k3p+n_kl(k,l) + accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,0,lambdap) + + enddo + + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + Vpseudo=f*accu + +else if(ac.ne.0.d0.and.bc.eq.0.d0)then + + f=fourpi**1.5d0 + theta_AC0=dacos( (a(3)-c(3))/ac ) + phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) + + areal=2.d0*g_a*ac + breal=2.d0*g_b*bc + freal=dexp(-g_a*ac**2-g_b*bc**2) + + do lambda=0,lmax+ntotA + do l=0,lmax + do ktot=0,ntotA+ntotB+nkl_max + do k=1,kmax + array_R(k,ktot,l,lambda,0)= int_prod_bessel(ktot+2,g_a+g_b+dz_kl(k,l),lambda,0,areal,breal,arg) + enddo + enddo + enddo + enddo + + do k1=0,n_a(1) + array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1) + enddo + do k2=0,n_a(2) + array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2) + enddo + do k3=0,n_a(3) + array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3) + enddo + + accu=0.d0 + do l=0,lmax + do m=-l,l + + do k3=0,n_a(3) + do k2=0,n_a(2) + do k1=0,n_a(1) + do lambda=0,l+ntotA + do mu=-lambda,lambda + array_I_A(mu,lambda,k1,k2,k3)=bigI(lambda,mu,l,m,k1,k2,k3) + enddo + enddo + enddo + enddo + enddo + + do k3=0,n_a(3) + if (array_coefs_A(k3,3) == 0.d0) cycle + do k2=0,n_a(2) + if (array_coefs_A(k2,2) == 0.d0) cycle + do k1=0,n_a(1) + if (array_coefs_A(k1,1) == 0.d0) cycle + do lambda=0,l+ntotA + do mu=-lambda,lambda + + prod=array_coefs_A(k1,1)*array_coefs_A(k2,2)*array_coefs_A(k3,3)*ylm(lambda,mu,theta_AC0,phi_AC0)*array_I_A(mu,lambda,k1,k2,k3) + if (prod == 0.d0) cycle + prodp=prod*bigI(0,0,l,m,n_b(1),n_b(2),n_b(3)) + + if (prodp == 0.d0) cycle + + do k=1,kmax + ktot=k1+k2+k3+ntotB+n_kl(k,l) + accu=accu+prodp*v_kl(k,l)*array_R(k,ktot,l,lambda,0) + enddo + + enddo + enddo + enddo + enddo + enddo + + enddo + enddo + + Vpseudo=f*accu +endif + +! _ +! |_ o ._ _. | o _ _ +! | | | | (_| | | _> (/_ +! + deallocate (array_R, array_I_A, array_I_B) + deallocate (array_coefs_A, array_coefs_B) + return +end + +! _ +! | | +!__ __ _ __ ___ ___ _ _ __| | ___ _ __ _ _ _ __ ___ +!\ \ / / | '_ \/ __|/ _ \ | | |/ _` |/ _ \ | '_ \| | | | '_ ` _ \ +! \ V / | |_) \__ \ __/ |_| | (_| | (_) | | | | | |_| | | | | | | +! \_/ | .__/|___/\___|\__,_|\__,_|\___/ |_| |_|\__,_|_| |_| |_| +! | | +! |_| + +double precision function Vpseudo_num(npts,rmax,lmax,kmax,v_kl,n_kl,dz_kl,a,n_a,g_a,b,n_b,g_b,c) +implicit none + + +! ___ +! | ._ ._ _|_ +! _|_ | | |_) |_| |_ +! | +double precision, intent(in) :: a(3),g_a,b(3),g_b,c(3) +integer, intent(in) :: lmax,kmax,npts +integer, intent(in) :: n_a(3),n_b(3), n_kl(kmax,0:lmax) +double precision, intent(in) :: v_kl(kmax,0:lmax),dz_kl(kmax,0:lmax) +double precision, intent(in) :: rmax + +! +! | _ _ _. | +! |_ (_) (_ (_| | +! + +integer :: l,m,k,kk +double precision ac(3),bc(3) +double precision dr,sum,rC +double precision overlap_orb_ylm_brute + +! _ +! / _. | _ | +! \_ (_| | (_ |_| | +! + +do l=1,3 + ac(l)=a(l)-c(l) + bc(l)=b(l)-c(l) +enddo + +dr=rmax/npts +sum=0.d0 +do l=0,lmax + do m=-l,l + do k=1,npts + rC=(k-1)*dr+dr/2.d0 + do kk=1,kmax + sum=sum+dr*v_kl(kk,l)*rC**(n_kl(kk,l)+2)*dexp(-dz_kl(kk,l)*rC**2) & + *overlap_orb_ylm_brute(npts,rC,n_a,ac,g_a,l,m) & + *overlap_orb_ylm_brute(npts,rC,n_b,bc,g_b,l,m) + enddo + enddo + enddo +enddo +Vpseudo_num=sum +return +end +!! Routine Vloc is a variation of formumla (66) +!! of Kahn Baybutt TRuhlar J.Chem.Phys. vol.65 3826 (1976) +!! without the projection operator +!! +!! Vloc= (4pi)**3/2* \sum_{k=1}^{klocmax} \sum_{l=0}^lmax \sum_{m=-l}^{l} +!!\sum{k1=0}^{nAx} \sum{k2=0}^{nAy} \sum{k3=0}^{nAz} +!! binom(nAx,k1)*binom(nAy,k2)*binom(nAz,k3) +!! *CAx**(nAx-k1)*CAy**(nAy-k2)*CAz**(nAz-k3)* +!! \sum{k1p=0}^{nBx} \sum{k2p=0}^{nBy} \sum{k3p=0}^{nBz} +!! binom(nBx,k1p)*binom(nBy,k2p)*binom(nBz,k3p) +!! *CBx**(nBx-k1p)*CBy**(nBy-k2p)*CBz**(nBz-k3p)* +!!\sum_{l=0}^lmax \sum_{m=-l}^{l} + +!! bigI(0,0,l,m,k1+k1p,k2+k2p,k3+k3p)*Y_{l m}(D_unit) +!! *v_k(k)* bigR(lambda,k1+k2+k3+k1p+k2p+k3p+n_k(k),g_a,g_b,AC,BC,dz_k(k)) +!! +!! nA=nAx+nAy+nAz +!! nB=nBx+nBy+nBz +!! D=(g_a AC+g_b BC) +!! D_unit= vect(D)/D +!! AC_x= A_x-C_x, etc. +!! BC=|B-C| +!! AC_unit= vect(AC)/AC +!! BC_unit= vect(BC)/BCA +!! +!! bigR(lambda,g_a,g_b,g_k,AC,BC) +!! = exp(-g_a* AC**2 -g_b* BC**2)* +!! I_loc= \int dx x**l *exp(-gam*x**2) M_n(ax) l=ktot+2 gam=g_a+g_b+dz_k(k) a=dreal n=l +!! M_n(x) modified spherical bessel function + + +double precision function Vloc(klocmax,v_k,n_k,dz_k,a,n_a,g_a,b,n_b,g_b,c) +implicit none +integer klocmax +double precision v_k(klocmax),dz_k(klocmax),crochet,bigA +integer n_k(klocmax) +double precision a(3),g_a,b(3),g_b,c(3),d(3) +integer n_a(3),n_b(3),ntotA,ntotB,ntot,m +integer i,l,k,ktot,k1,k2,k3,k1p,k2p,k3p +double precision f,fourpi,ac,bc,freal,d2,dreal,theta_DC0,phi_DC0,coef +double precision,allocatable :: array_R_loc(:,:,:) +double precision,allocatable :: array_coefs(:,:,:,:,:,:) +double precision int_prod_bessel_loc,binom_func,accu,prod,ylm,bigI,arg + + + fourpi=4.d0*dacos(-1.d0) + f=fourpi**1.5d0 + ac=dsqrt((a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2) + bc=dsqrt((b(1)-c(1))**2+(b(2)-c(2))**2+(b(3)-c(3))**2) + arg=g_a*ac**2+g_b*bc**2 + if(arg.gt.-dlog(1.d-20))then + Vloc=0.d0 + return + endif + + ntotA=n_a(1)+n_a(2)+n_a(3) + ntotB=n_b(1)+n_b(2)+n_b(3) + ntot=ntotA+ntotB + + if(ac.eq.0.d0.and.bc.eq.0.d0)then + accu=0.d0 + + do k=1,klocmax + accu=accu+v_k(k)*crochet(n_k(k)+2+ntot,g_a+g_b+dz_k(k)) + enddo + Vloc=accu*fourpi*bigI(0,0,0,0,n_a(1)+n_b(1),n_a(2)+n_b(2),n_a(3)+n_b(3)) + !bigI frequently is null + return + endif + + freal=dexp(-g_a*ac**2-g_b*bc**2) + + d2 = 0.d0 + do i=1,3 + d(i)=g_a*(a(i)-c(i))+g_b*(b(i)-c(i)) + d2=d2+d(i)*d(i) + enddo + d2=dsqrt(d2) + dreal=2.d0*d2 + + + allocate (array_R_loc(-2:ntot+klocmax,klocmax,0:ntot)) + allocate (array_coefs(0:ntot,0:ntot,0:ntot,0:ntot,0:ntot,0:ntot)) + + do ktot=-2,ntotA+ntotB+klocmax + do l=0,ntot + do k=1,klocmax + array_R_loc(ktot,k,l)=freal*int_prod_bessel_loc(ktot+2,g_a+g_b+dz_k(k),l,dreal) + enddo + enddo + enddo + + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + array_coefs(k1,k2,k3,k1p,k2p,k3p)=binom_func(n_a(1),k1)*binom_func(n_a(2),k2)*binom_func(n_a(3),k3)& + *(c(1)-a(1))**(n_a(1)-k1)*(c(2)-a(2))**(n_a(2)-k2)*(c(3)-a(3))**(n_a(3)-k3)& + *binom_func(n_b(1),k1p)*binom_func(n_b(2),k2p)*binom_func(n_b(3),k3p)& + *(c(1)-b(1))**(n_b(1)-k1p)*(c(2)-b(2))**(n_b(2)-k2p)*(c(3)-b(3))**(n_b(3)-k3p) + enddo + enddo + enddo + enddo + enddo + enddo + + + accu=0.d0 + if(d2 == 0.d0)then + l=0 + m=0 + coef=1.d0/dsqrt(4.d0*dacos(-1.d0)) + do k=1,klocmax + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) & + *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) + ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) + accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) + enddo + enddo + enddo + enddo + enddo + enddo + enddo + + else + theta_DC0=dacos(d(3)/d2) + phi_DC0=datan2(d(2)/d2,d(1)/d2) + + do k=1,klocmax + if (v_k(k) == 0.d0) cycle + do k1=0,n_a(1) + do k2=0,n_a(2) + do k3=0,n_a(3) + do k1p=0,n_b(1) + do k2p=0,n_b(2) + do k3p=0,n_b(3) + if (array_coefs(k1,k2,k3,k1p,k2p,k3p) == 0.d0) cycle + do l=0,ntot + do m=-l,l + coef=ylm(l,m,theta_DC0,phi_DC0) + if (coef == 0.d0) cycle + ktot=k1+k2+k3+k1p+k2p+k3p+n_k(k) + if (array_R_loc(ktot,k,l) == 0.d0) cycle + prod=coef*array_coefs(k1,k2,k3,k1p,k2p,k3p) & + *bigI(l,m,0,0,k1+k1p,k2+k2p,k3+k3p) + accu=accu+prod*v_k(k)*array_R_loc(ktot,k,l) + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + enddo + endif + Vloc=f*accu + + deallocate (array_R_loc) + deallocate (array_coefs) +end + +double precision function bigA(i,j,k) +implicit none +integer i,j,k +double precision fourpi,dble_fact +fourpi=4.d0*dacos(-1.d0) +bigA=0.d0 +if(mod(i,2).eq.1)return +if(mod(j,2).eq.1)return +if(mod(k,2).eq.1)return +bigA=fourpi*dble_fact(i-1)*dble_fact(j-1)*dble_fact(k-1)/dble_fact(i+j+k+1) +end +!! +!! I_{lambda,mu,l,m}^{k1,k2,k3} = /int dOmega Y_{lambda mu} xchap^k1 ychap^k2 zchap^k3 Y_{lm} +!! + +double precision function bigI(lambda,mu,l,m,k1,k2,k3) +implicit none +integer lambda,mu,l,m,k1,k2,k3 +integer k,i,kp,ip +double precision pi,sum,factor1,factor2,cylm,cylmp,bigA,binom_func,fact,coef_pm +double precision sgn, sgnp +pi=dacos(-1.d0) + +bigI=0.d0 +if(mu.gt.0.and.m.gt.0)then +sum=0.d0 +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) +if (factor2== 0.d0) return +sgn = 1.d0 +do k=0,mu/2 + do i=0,lambda-mu + if (coef_pm(lambda,i+mu) == 0.d0) cycle + sgnp = 1.d0 + do kp=0,m/2 + do ip=0,l-m + cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + if (cylm == 0.d0) cycle + cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(mu-2*k+m-2*kp+k1,2*k+2*kp+k2,i+ip+k3) + enddo + sgnp = -sgnp + enddo + enddo + sgn = -sgn +enddo +bigI=sum +return +endif + +if(mu.eq.0.and.m.eq.0)then +factor1=dsqrt((2*lambda+1)/(4.d0*pi)) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)/(4.d0*pi)) +if (factor2== 0.d0) return +sum=0.d0 +do i=0,lambda + do ip=0,l + cylm=factor1*coef_pm(lambda,i) + if (cylm == 0.d0) cycle + cylmp=factor2*coef_pm(l,ip) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(k1,k2,i+ip+k3) + enddo +enddo +bigI=sum +return +endif + +if(mu.eq.0.and.m.gt.0)then +factor1=dsqrt((2*lambda+1)/(4.d0*pi)) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) +if (factor2== 0.d0) return +sum=0.d0 +do i=0,lambda + sgnp = 1.d0 + do kp=0,m/2 + do ip=0,l-m + cylm=factor1*coef_pm(lambda,i) + if (cylm == 0.d0) cycle + cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(m-2*kp+k1,2*kp+k2,i+ip+k3) + enddo + sgnp = -sgnp + enddo +enddo +bigI=sum +return +endif + +if(mu.gt.0.and.m.eq.0)then +sum=0.d0 +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)/(4.d0*pi)) +if (factor2== 0.d0) return +sgn = 1.d0 +do k=0,mu/2 + do i=0,lambda-mu + if (coef_pm(lambda,i+mu) == 0.d0) cycle + do ip=0,l + cylm=sgn*factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + if (cylm == 0.d0) cycle + cylmp=factor2*coef_pm(l,ip) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(mu-2*k +k1,2*k +k2,i+ip +k3) + enddo + enddo + sgn = -sgn +enddo +bigI=sum +return +endif + +if(mu.lt.0.and.m.lt.0)then +mu=-mu +m=-m +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) +if (factor2== 0.d0) return +sum=0.d0 +sgn = 1.d0 +do k=0,(mu-1)/2 + do i=0,lambda-mu + if (coef_pm(lambda,i+mu) == 0.d0) cycle + sgnp = 1.d0 + do kp=0,(m-1)/2 + do ip=0,l-m + if (coef_pm(l,ip+m) == 0.d0) cycle + cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + if (cylm == 0.d0) cycle + cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-(2*kp+1)+k1,(2*k+1)+(2*kp+1)+k2,i+ip+k3) + enddo + sgnp = -sgnp + enddo + enddo + sgn = -sgn +enddo +mu=-mu +m=-m +bigI=sum +return +endif + +if(mu.eq.0.and.m.lt.0)then +m=-m +factor1=dsqrt((2*lambda+1)/(4.d0*pi)) +if (factor1 == 0.d0) return +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) +if (factor2 == 0.d0) return +sum=0.d0 +do i=0,lambda + sgnp = 1.d0 + do kp=0,(m-1)/2 + do ip=0,l-m + cylm=factor1*coef_pm(lambda,i) + if (cylm == 0.d0) cycle + cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(m-(2*kp+1)+k1,2*kp+1+k2,i+ip+k3) + enddo + sgnp = -sgnp + enddo +enddo +m=-m +bigI=sum +return +endif + +if(mu.lt.0.and.m.eq.0)then +sum=0.d0 +mu=-mu +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)/(4.d0*pi)) +if (factor2== 0.d0) return +sgn = 1.d0 +do k=0,(mu-1)/2 + do i=0,lambda-mu + do ip=0,l + cylm=sgn*factor1*binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + if (cylm == 0.d0) cycle + cylmp=factor2*coef_pm(l,ip) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+k1,2*k+1+k2,i+ip+k3) + enddo + enddo + sgn = -sgn +enddo +mu=-mu +bigI=sum +return +endif + +if(mu.gt.0.and.m.lt.0)then +sum=0.d0 +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) +if (factor2== 0.d0) return +m=-m +sgn=1.d0 +do k=0,mu/2 + do i=0,lambda-mu + if (coef_pm(lambda,i+mu) == 0.d0) cycle + sgnp=1.d0 + do kp=0,(m-1)/2 + do ip=0,l-m + if (coef_pm(l,ip+m) == 0.d0) cycle + cylm =sgn *factor1*binom_func(mu,2*k)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + if (cylm == 0.d0) cycle + cylmp=sgnp*factor2*binom_func(m,2*kp+1)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(mu-2*k+m-(2*kp+1)+k1,2*k+2*kp+1+k2,i+ip+k3) + enddo + sgnp = -sgnp + enddo + enddo + sgn = -sgn +enddo +m=-m +bigI=sum +return +endif + +if(mu.lt.0.and.m.gt.0)then +mu=-mu +factor1=dsqrt((2*lambda+1)*fact(lambda-iabs(mu))/(2.d0*pi*fact(lambda+iabs(mu)))) +if (factor1== 0.d0) return +factor2=dsqrt((2*l+1)*fact(l-iabs(m))/(2.d0*pi*fact(l+iabs(m)))) +if (factor2== 0.d0) return +sum=0.d0 +sgn = 1.d0 +do k=0,(mu-1)/2 + do i=0,lambda-mu + if (coef_pm(lambda,i+mu) == 0.d0) cycle + sgnp = 1.d0 + do kp=0,m/2 + do ip=0,l-m + if (coef_pm(l,ip+m) == 0.d0) cycle + cylm=sgn*factor1 *binom_func(mu,2*k+1)*fact(mu+i)/fact(i)*coef_pm(lambda,i+mu) + if (cylm == 0.d0) cycle + cylmp=sgnp*factor2*binom_func(m,2*kp)*fact(m+ip)/fact(ip)*coef_pm(l,ip+m) + if (cylmp == 0.d0) cycle + sum=sum+cylm*cylmp*bigA(mu-(2*k+1)+m-2*kp+k1,2*k+1+2*kp+k2,i+ip+k3) + enddo + sgnp = -sgnp + enddo + enddo + sgn = -sgn +enddo +bigI=sum +mu=-mu +return +endif + +stop 'pb in bigI!' +end + +double precision function crochet(n,g) +implicit none +integer n +double precision g,dble_fact,expo +double precision, parameter :: sq_pi_ov_2=dsqrt(dacos(-1.d0)*0.5d0) +expo=0.5d0*dfloat(n+1) +crochet=dble_fact(n-1)/(g+g)**expo +if(mod(n,2).eq.0)crochet=crochet*sq_pi_ov_2 +end + +!! +!! overlap= = /int dOmega Ylm (x-center_x)**nx*(y-center_y)**nx*(z-center)**nx +!! *exp(-g*(r-center)**2) +!! +double precision function overlap_orb_ylm_brute(npts,r,npower_orb,center_orb,g_orb,l,m) +implicit none +integer npower_orb(3),l,m,i,j,npts +double precision u,g_orb,du,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3) +double precision x_orb,y_orb,z_orb,twopi,r +twopi=2.d0*dacos(-1.d0) +du=2.d0/npts +dphi=twopi/npts +overlap_orb_ylm_brute=0.d0 +do i=1,npts + u=-1.d0+du*(i-1)+du/2.d0 + sintheta=dsqrt(1.d0-u**2) + do j=1,npts + phi=dphi*(j-1)+dphi/2.d0 + x_orb=r*dcos(phi)*sintheta + y_orb=r*dsin(phi)*sintheta + z_orb=r*u + term=orb_phi(x_orb,y_orb,z_orb,npower_orb,center_orb,g_orb)*ylm_real(l,m,u,phi) + overlap_orb_ylm_brute= overlap_orb_ylm_brute+term*du*dphi + enddo +enddo +end + +double precision function overlap_orb_ylm_grid(nptsgrid,r_orb,npower_orb,center_orb,g_orb,l,m) +implicit none +!! PSEUDOS +integer nptsgridmax,nptsgrid +double precision coefs_pseudo,ptsgrid +parameter(nptsgridmax=50) +common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) +!!!!! +integer npower_orb(3),l,m,i +double precision x,g_orb,two_pi,dx,dphi,term,orb_phi,ylm_real,sintheta,r_orb,phi,center_orb(3) +double precision x_orb,y_orb,z_orb,twopi,pi,cosphi,sinphi,xbid +pi=dacos(-1.d0) +twopi=2.d0*pi +overlap_orb_ylm_grid=0.d0 +do i=1,nptsgrid + x_orb=r_orb*ptsgrid(i,1) + y_orb=r_orb*ptsgrid(i,2) + z_orb=r_orb*ptsgrid(i,3) + x=ptsgrid(i,3) + phi=datan2(ptsgrid(i,2),ptsgrid(i,1)) + term=orb_phi(x_orb,y_orb,z_orb,npower_orb,center_orb,g_orb)*ylm_real(l,m,x,phi) + overlap_orb_ylm_grid= overlap_orb_ylm_grid+coefs_pseudo(i)*term +enddo +overlap_orb_ylm_grid=2.d0*twopi*overlap_orb_ylm_grid +end + +! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) exp(i m phi) +! l=0,1,2,.... +! m=0,1,...,l +! Here: +! n=l (n=0,1,...) +! m=0,1,...,n +! x=cos(theta) 0 < x < 1 +! +! +! This routine computes: PM(m,n) for n=0,...,N (number N in input) and m=0,..,n + +! Exemples (see 'Associated Legendre Polynomilas wikipedia') +! P_{0}^{0}(x)=1 +! P_{1}^{-1}(x)=-1/2 P_{1}^{1}(x) +! P_{1}^{0}(x)=x +! P_{1}^{1}(x)=-(1-x^2)^{1/2} +! P_{2}^{-2}(x)=1/24 P_{2}^{2}(x) +! P_{2}^{-1}(x)=-1/6 P_{2}^{1}(x) +! P_{2}^{0}(x)=1/2 (3x^{2}-1) +! P_{2}^{1}(x)=-3x(1-x^2)^{1/2} +! P_{2}^{2}(x)=3(1-x^2) + + + SUBROUTINE LPMN(MM,M,N,X,PM) +! +! Here N = LMAX +! Here M= MMAX (we take M=LMAX in input) +! +! ===================================================== +! Purpose: Compute the associated Legendre functions Pmn(x) +! Input : x --- Argument of Pmn(x) +! m --- Order of Pmn(x), m = 0,1,2,...,n +! n --- Degree of Pmn(x), n = 0,1,2,...,N +! mm --- Physical dimension of PM +! Output: PM(m,n) --- Pmn(x) +! ===================================================== +! + implicit none +! IMPLICIT DOUBLE PRECISION (P,X) + integer :: MM, N, I, J, M + double precision :: PM(0:MM,0:(N+1)), X, XQ, XS + DOUBLE PRECISION, SAVE :: INVERSE(100) = 0.D0 + DOUBLE PRECISION :: LS, II, JJ + IF (INVERSE(1) == 0.d0) THEN + DO I=1,100 + INVERSE(I) = 1.D0/DBLE(I) + ENDDO + ENDIF + DO I=0,N + DO J=0,M + PM(J,I)=0.0D0 + ENDDO + ENDDO + PM(0,0)=1.0D0 + IF (DABS(X).EQ.1.0D0) THEN + DO I=1,N + PM(0,I)=X**I + ENDDO + RETURN + ENDIF + LS=1.D0 + IF (DABS(X).GT.1.0D0) LS=-1.D0 + XQ=DSQRT(LS*(1.0D0-X*X)) + XS=LS*(1.0D0-X*X) + II = 1.D0 + DO I=1,M + PM(I,I)=-LS*II*XQ*PM(I-1,I-1) + II = II+2.D0 + ENDDO + II = 1.D0 + DO I=0,M + PM(I,I+1)=II*X*PM(I,I) + II = II+2.D0 + ENDDO + + II = 0.D0 + DO I=0,M + JJ = II+2.D0 + DO J=I+2,N + PM(I,J)=((2.0D0*JJ-1.0D0)*X*PM(I,J-1)- (II+JJ-1.0D0)*PM(I,J-2))*INVERSE(J-I) + JJ = JJ+1.D0 + ENDDO + II = II+1.D0 + ENDDO + END + + +! Y_l^m(theta,phi) = i^(m+|m|) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 +! P_l^|m|(cos(theta)) exp(i m phi) + + subroutine erreur(x,n,rmoy,error) + implicit none + integer :: i, n + double precision :: x(n), rn, rn1, error, rmoy +! calcul de la moyenne + rmoy=0.d0 + do i=1,n + rmoy=rmoy+x(i) + enddo + rmoy=rmoy/dfloat(n) +! calcul de l'erreur + error=0.d0 + do i=1,n + error=error+(x(i)-rmoy)**2 + enddo + if(n.gt.1)then + rn=dfloat(n) + rn1=dfloat(n-1) + error=dsqrt(error)/dsqrt(rn*rn1) + else + write(2,*)'Seulement un block Erreur nondefinie' + error=0.d0 + endif + end + + subroutine initpseudos(nptsgrid) + implicit none + integer nptsgridmax,nptsgrid,ik + double precision coefs_pseudo,ptsgrid + double precision p,q,r,s + parameter(nptsgridmax=50) + common/pseudos/coefs_pseudo(nptsgridmax),ptsgrid(nptsgridmax,3) + + p=1.d0/dsqrt(2.d0) + q=1.d0/dsqrt(3.d0) + r=1.d0/dsqrt(11.d0) + s=3.d0/dsqrt(11.d0) + + if(nptsgrid.eq.4)then + + ptsgrid(1,1)=q + ptsgrid(1,2)=q + ptsgrid(1,3)=q + + ptsgrid(2,1)=q + ptsgrid(2,2)=-q + ptsgrid(2,3)=-q + + ptsgrid(3,1)=-q + ptsgrid(3,2)=q + ptsgrid(3,3)=-q + + ptsgrid(4,1)=-q + ptsgrid(4,2)=-q + ptsgrid(4,3)=q + + do ik=1,4 + coefs_pseudo(ik)=1.d0/4.d0 + enddo + return + endif + + ptsgrid(1,1)=1.d0 + ptsgrid(1,2)=0.d0 + ptsgrid(1,3)=0.d0 + + ptsgrid(2,1)=-1.d0 + ptsgrid(2,2)=0.d0 + ptsgrid(2,3)=0.d0 + + ptsgrid(3,1)=0.d0 + ptsgrid(3,2)=1.d0 + ptsgrid(3,3)=0.d0 + + ptsgrid(4,1)=0.d0 + ptsgrid(4,2)=-1.d0 + ptsgrid(4,3)=0.d0 + + ptsgrid(5,1)=0.d0 + ptsgrid(5,2)=0.d0 + ptsgrid(5,3)=1.d0 + + ptsgrid(6,1)=0.d0 + ptsgrid(6,2)=0.d0 + ptsgrid(6,3)=-1.d0 + + do ik=1,6 + coefs_pseudo(ik)=1.d0/6.d0 + enddo + + if(nptsgrid.eq.6)return + + ptsgrid(7,1)=p + ptsgrid(7,2)=p + ptsgrid(7,3)=0.d0 + + ptsgrid(8,1)=p + ptsgrid(8,2)=-p + ptsgrid(8,3)=0.d0 + + ptsgrid(9,1)=-p + ptsgrid(9,2)=p + ptsgrid(9,3)=0.d0 + + ptsgrid(10,1)=-p + ptsgrid(10,2)=-p + ptsgrid(10,3)=0.d0 + + ptsgrid(11,1)=p + ptsgrid(11,2)=0.d0 + ptsgrid(11,3)=p + + ptsgrid(12,1)=p + ptsgrid(12,2)=0.d0 + ptsgrid(12,3)=-p + + ptsgrid(13,1)=-p + ptsgrid(13,2)=0.d0 + ptsgrid(13,3)=p + + ptsgrid(14,1)=-p + ptsgrid(14,2)=0.d0 + ptsgrid(14,3)=-p + + ptsgrid(15,1)=0.d0 + ptsgrid(15,2)=p + ptsgrid(15,3)=p + + ptsgrid(16,1)=0.d0 + ptsgrid(16,2)=p + ptsgrid(16,3)=-p + + ptsgrid(17,1)=0.d0 + ptsgrid(17,2)=-p + ptsgrid(17,3)=p + + ptsgrid(18,1)=0.d0 + ptsgrid(18,2)=-p + ptsgrid(18,3)=-p + + do ik=1,6 + coefs_pseudo(ik)=1.d0/30.d0 + enddo + do ik=7,18 + coefs_pseudo(ik)=1.d0/15.d0 + enddo + + if(nptsgrid.eq.18)return + + ptsgrid(19,1)=q + ptsgrid(19,2)=q + ptsgrid(19,3)=q + + ptsgrid(20,1)=-q + ptsgrid(20,2)=q + ptsgrid(20,3)=q + + ptsgrid(21,1)=q + ptsgrid(21,2)=-q + ptsgrid(21,3)=q + + ptsgrid(22,1)=q + ptsgrid(22,2)=q + ptsgrid(22,3)=-q + + ptsgrid(23,1)=-q + ptsgrid(23,2)=-q + ptsgrid(23,3)=q + + ptsgrid(24,1)=-q + ptsgrid(24,2)=q + ptsgrid(24,3)=-q + + ptsgrid(25,1)=q + ptsgrid(25,2)=-q + ptsgrid(25,3)=-q + + ptsgrid(26,1)=-q + ptsgrid(26,2)=-q + ptsgrid(26,3)=-q + + do ik=1,6 + coefs_pseudo(ik)=1.d0/21.d0 + enddo + do ik=7,18 + coefs_pseudo(ik)=4.d0/105.d0 + enddo + do ik=19,26 + coefs_pseudo(ik)=27.d0/840.d0 + enddo + + if(nptsgrid.eq.26)return + + ptsgrid(27,1)=r + ptsgrid(27,2)=r + ptsgrid(27,3)=s + + ptsgrid(28,1)=r + ptsgrid(28,2)=-r + ptsgrid(28,3)=s + + ptsgrid(29,1)=-r + ptsgrid(29,2)=r + ptsgrid(29,3)=s + + ptsgrid(30,1)=-r + ptsgrid(30,2)=-r + ptsgrid(30,3)=s + + ptsgrid(31,1)=r + ptsgrid(31,2)=r + ptsgrid(31,3)=-s + + ptsgrid(32,1)=r + ptsgrid(32,2)=-r + ptsgrid(32,3)=-s + + ptsgrid(33,1)=-r + ptsgrid(33,2)=r + ptsgrid(33,3)=-s + + ptsgrid(34,1)=-r + ptsgrid(34,2)=-r + ptsgrid(34,3)=-s + + ptsgrid(35,1)=r + ptsgrid(35,2)=s + ptsgrid(35,3)=r + + ptsgrid(36,1)=-r + ptsgrid(36,2)=s + ptsgrid(36,3)=r + + ptsgrid(37,1)=r + ptsgrid(37,2)=s + ptsgrid(37,3)=-r + + ptsgrid(38,1)=-r + ptsgrid(38,2)=s + ptsgrid(38,3)=-r + + ptsgrid(39,1)=r + ptsgrid(39,2)=-s + ptsgrid(39,3)=r + + ptsgrid(40,1)=r + ptsgrid(40,2)=-s + ptsgrid(40,3)=-r + + ptsgrid(41,1)=-r + ptsgrid(41,2)=-s + ptsgrid(41,3)=r + + ptsgrid(42,1)=-r + ptsgrid(42,2)=-s + ptsgrid(42,3)=-r + + ptsgrid(43,1)=s + ptsgrid(43,2)=r + ptsgrid(43,3)=r + + ptsgrid(44,1)=s + ptsgrid(44,2)=r + ptsgrid(44,3)=-r + + ptsgrid(45,1)=s + ptsgrid(45,2)=-r + ptsgrid(45,3)=r + + ptsgrid(46,1)=s + ptsgrid(46,2)=-r + ptsgrid(46,3)=-r + + ptsgrid(47,1)=-s + ptsgrid(47,2)=r + ptsgrid(47,3)=r + + ptsgrid(48,1)=-s + ptsgrid(48,2)=r + ptsgrid(48,3)=-r + + ptsgrid(49,1)=-s + ptsgrid(49,2)=-r + ptsgrid(49,3)=r + + ptsgrid(50,1)=-s + ptsgrid(50,2)=-r + ptsgrid(50,3)=-r + + do ik=1,6 + coefs_pseudo(ik)=4.d0/315.d0 + enddo + do ik=7,18 + coefs_pseudo(ik)=64.d0/2835.d0 + enddo + do ik=19,26 + coefs_pseudo(ik)=27.d0/1280.d0 + enddo + do ik=27,50 + coefs_pseudo(ik)=14641.d0/725760.d0 + enddo + + if(nptsgrid.eq.50)return + + write(*,*)'Grid for pseudos not available!' + write(*,*)'N=4-6-18-26-50 only!' + stop + end + +!! +!! R_{lambda,lamba',N}= exp(-ga_a AC**2 -g_b BC**2) /int_{0}{+infty} r**(2+n) exp(-(g_a+g_b+g_k)r**2) +!! * M_{lambda}( 2g_a ac r) M_{lambda'}(2g_b bc r) +!! + double precision function bigR(lambda,lambdap,n,g_a,g_b,ac,bc,g_k) + implicit none + integer lambda,lambdap,n,npts,i + double precision g_a,g_b,ac,bc,g_k,arg,factor,delta1,delta2,cc,rmax,dr,sum,x1,x2,r + double precision bessel_mod + arg=g_a*ac**2+g_b*bc**2 + factor=dexp(-arg) + delta1=2.d0*g_a*ac + delta2=2.d0*g_b*bc + cc=g_a+g_b+g_k + if(cc.eq.0.d0)stop 'pb. in bigR' + rmax=dsqrt(-dlog(10.d-20)/cc) + npts=500 + dr=rmax/npts + sum=0.d0 + do i=1,npts + r=(i-1)*dr + x1=delta1*r + x2=delta2*r + sum=sum+dr*r**(n+2)*dexp(-cc*r*r)*bessel_mod(x1,lambda)*bessel_mod(x2,lambdap) + enddo + bigR=sum*factor + end + + double precision function bessel_mod(x,n) + implicit none + integer n + double precision x,bessel_mod_exp,bessel_mod_recur + if(x.le.0.8d0)then + bessel_mod=bessel_mod_exp(n,x) + else + bessel_mod=bessel_mod_recur(n,x) + endif + end + + recursive function bessel_mod_recur(n,x) result(a) + implicit none + integer n + double precision x,a,bessel_mod_exp + if(x.le.0.8d0)then + a=bessel_mod_exp(n,x) + return + endif + if(n.eq.0)a=dsinh(x)/x + if(n.eq.1)a=(x*dcosh(x)-dsinh(x))/(x*x) + if(n.ge.2)a=bessel_mod_recur(n-2,x)-(n+n-1)/x*bessel_mod_recur(n-1,x) + end + + double precision function bessel_mod_exp(n,x) + implicit none + integer n,k + double precision x,coef,accu,fact,dble_fact + accu=0.d0 + do k=0,10 + coef=1.d0/(fact(k)*dble_fact(2*(n+k)+1)) + accu=accu+(0.5d0*x*x)**k*coef + enddo + bessel_mod_exp=x**n*accu + end + + + + +!c Computation of real spherical harmonics Ylm(theta,phi) +!c +!c l=0,1,.... +!c m=-l,l +!c +!c m>0: Y_lm = sqrt(2) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) cos(m phi) +!c m=0: Y_l0 = ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^0 (cos(theta)) +!c m<0: Y_lm = sqrt(2) ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 P_l^|m|(cos(theta)) sin(|m|phi) + +!Examples(wikipedia http://en.wikipedia.org/wiki/Table_of_spherical_harmonics#Real_spherical_harmonics) + +! l = 0 + +! Y_00 = \sqrt{1/(4pi)} + +! l = 1 + +! Y_1-1= \sqrt{3/(4pi)} y/r +! Y_10 = \sqrt{3/(4pi)} z/r +! Y_11 = \sqrt{3/(4pi)} x/r +! +! l = 2 +! +! Y_2,-2= 1/2 \sqrt{15/pi} xy/r^2 +! Y_2,-1= 1/2 \sqrt{15/pi} yz/r^2 +! Y_20 = 1/4 \sqrt{15/pi} (-x^2-y^2 +2z^2)/r^2 +! Y_21 = 1/2 \sqrt{15/pi} zx/r^2 +! Y_22 = 1/4 \sqrt{15/pi} (x^2-y^2)/r^2 +! +!c +double precision function ylm(l,m,theta,phi) +implicit none +integer l,m,i +double precision theta,phi,pm,factor,twopi,x,fact,sign +DIMENSION PM(0:100,0:100) +twopi=2.d0*dacos(-1.d0) +x=dcos(theta) +if (iand(m,1) == 1) then + sign=-1.d0 +else + sign=1.d0 +endif +CALL LPMN(100,l,l,X,PM) +if (m > 0) then + factor=dsqrt((l+l+1)*fact(l-m) /(twopi*fact(l+m)) ) +! factor = dble(l+m) +! do i=-m,m-1 +! factor = factor * (factor - 1.d0) +! enddo +! factor=dsqrt(dble(l+l+1)/(twopi*factor) ) + ylm=sign*factor*pm(m,l)*dcos(dfloat(m)*phi) +else if (m == 0) then + factor=dsqrt( 0.5d0*(l+l+1) /twopi ) + ylm=factor*pm(m,l) +else if (m < 0) then + factor=dsqrt( (l+l+1)*fact(l+m) /(twopi*fact(l-m)) ) +! factor = dble(l-m) +! do i=m,-m-1 +! factor = factor * (factor - 1.d0) +! enddo +! factor=dsqrt(dble(l+l+1)/(twopi*factor) ) + ylm=sign*factor*pm(-m,l)*dsin(dfloat(-m)*phi) +endif +end + +!c Explicit representation of Legendre polynomials P_n(x) +!! +!! P_n0(x) = P_n(x)= \sum_{k=0}^n a_k x^k +!! +!! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n ) +!! coef_pm(n,k) is the k_th coefficient of P_n(x) +double precision function coef_pm(n,k) +implicit none +integer n,k +double precision arg,binom_func,binom_gen +if(n.eq.0.and.k.ne.0)stop 'coef_pm not defined' +if(n.eq.0.and.k.eq.0)then +coef_pm=1.d0 +return +endif +arg=0.5d0*dfloat(n+k-1) +coef_pm=2.d0**n*binom_func(n,k)*binom_gen(arg,n) +end + +!! Ylm_bis uses the series expansion of Ylm in xchap^i ychap^j zchap^k +!! xchap=x/r etc. +!c m>0: Y_lm = sqrt(2)*factor* P_l^|m|(cos(theta)) cos(m phi) +!c m=0: Y_l0 = factor* P_l^0 (cos(theta)) +!c m<0: Y_lm = sqrt(2) factor* P_l^|m|(cos(theta)) sin(|m|phi) +!c factor= ([(2l+1)*(l-|m|)!]/[4pi*(l+|m|)!])^1/2 + +!! P_l^m (x) = (-1)**m (1-x**2)^m/2 d^m/dx^m P_l(x) m >0 or 0 +!! the series expansion of P_m (x) is known +!! +!! sin(theta)**m cos(mphi) = \sum_0^[m/2] binom(m,2k) x^(m-2k) y^2k (-1)**k (easy to proove with +!! Moivre formula) +!! (here x = xchap...) +!! +!! Ylm m> 0 = \sum_{k=0}^[m/2] \sum_{i=0}^(l-m) c_ki x^(m-2k) y^2k z^i +!! +!! c_ki= (-1)^k sqrt(2)*factor*binom(m,2k)*(m+i)!/i!*coef_pm(l,i+m) +!! +!! Ylm m< 0 = \sum_{k=0}^[(m-1)/2] \sum_{i=0}^(l-m) c_ki x^(m-(2k+1)) y^(2k+1) z^i +!! +!! c_ki= (-1)^k sqrt(2)*factor*binom(m,2k+1)*(m+i)!/i!*coef_pm(l,i+m) + + +double precision function ylm_bis(l,m,theta,phi) +implicit none +integer l,m,k,i +double precision x,y,z,theta,phi,sum,factor,pi,binom_func,fact,coef_pm,cylm +pi=dacos(-1.d0) +x=dsin(theta)*dcos(phi) +y=dsin(theta)*dsin(phi) +z=dcos(theta) +factor=dsqrt((2*l+1)*fact(l-iabs(m))/(4.d0*pi*fact(l+iabs(m)))) +if(m.gt.0)then +sum=0.d0 +do k=0,m/2 + do i=0,l-m + cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k)*fact(m+i)/fact(i)*coef_pm(l,i+m) + sum=sum+cylm*x**(m-2*k)*y**(2*k)*z**i + enddo +enddo +ylm_bis=sum +return +endif +if(m.eq.0)then +sum=0.d0 +do i=0,l + sum=sum+factor*coef_pm(l,i)*z**i +enddo +ylm_bis=sum +return +endif +if(m.lt.0)then +m=-m +sum=0.d0 +do k=0,(m-1)/2 + do i=0,l-m + cylm=(-1.d0)**k*factor*dsqrt(2.d0)*binom_func(m,2*k+1)*fact(m+i)/fact(i)*coef_pm(l,i+m) + sum=sum+cylm*x**(m-(2*k+1))*y**(2*k+1)*z**i + enddo +enddo +ylm_bis=sum +m=-m +return +endif +end + +!c +!c Computation of associated Legendre Polynomials PM(m,n) for n=0,...,N +!c Here: +!c n=l (n=0,1,...) +!c m=0,1,...,n +!c x=cos(theta) 0 < x < 1 +!c +!c This routine computes: PM(m,n) for n=0,...,N (number N in input) and m=0,..,n +!c Exemples (see 'Associated Legendre Polynomilas wikipedia') +!c P_{0}^{0}(x)=1 +!c P_{1}^{-1}(x)=-1/2 P_{1}^{1}(x) +!c P_{1}^{0}(x)=x +!c P_{1}^{1}(x)=-(1-x^2)^{1/2} +!c P_{2}^{-2}(x)=1/24 P_{2}^{2}(x) +!c P_{2}^{-1}(x)=-1/6 P_{2}^{1}(x) +!c P_{2}^{0}(x)=1/2 (3x^{2}-1) +!c P_{2}^{1}(x)=-3x(1-x^2)^{1/2} +!c P_{2}^{2}(x)=3(1-x^2) +!c +!c Explicit representation: +!! +!! P_n0(x) = P_n(x)= \sum_{k=0}^n a_k x^k +!! +!! with a_k= 2^n binom(n,k) binom( (n+k-1)/2, n ) + +double precision function binom_gen(alpha,n) + implicit none + integer :: n,k + double precision :: fact,alpha,prod, factn_inv + prod=1.d0 + factn_inv = 1.d0/fact(n) + do k=0,n-1 + prod=prod*(alpha-k) + binom_gen = prod*factn_inv + enddo +end + + +double precision function coef_nk(n,k) + implicit none + integer n,k + + double precision gam,dble_fact,fact + + if (k<0) stop 'pseudopot.f90 : coef_nk' + if (k>63) stop 'pseudopot.f90 : coef_nk' + gam=dble_fact(n+n+k+k+1) +! coef_nk=1.d0/(2.d0**k*fact(k)*gam) + coef_nk=1.d0/(dble(ibset(0_8,k))*fact(k)*gam) + + return + +end + +!! Calculation of +!! +!! I= \int dx x**l *exp(-gam*x**2) M_n(ax) M_m(bx) +!! +!! M_n(x) modified spherical bessel function +!! + +double precision function int_prod_bessel(l,gam,n,m,a,b,arg) + + implicit none + integer n,k,m,q,l,kcp + double precision gam,dble_fact,fact,pi,a,b + double precision int,intold,sum,coef_nk,crochet,u,int_prod_bessel_large,freal,arg + + integer :: n_1,m_1,nlm + double precision :: term_A, term_B, term_rap, expo + double precision :: s_q_0, s_q_k, s_0_0, a_over_b_square + double precision :: int_prod_bessel_loc + double precision :: inverses(0:300) + double precision :: two_qkmp1, qk, mk, nk + + logical done + + u=(a+b)*0.5d0/dsqrt(gam) + freal=dexp(-arg) + + if(a.eq.0.d0.and.b.eq.0.d0)then + if(n.ne.0.or.m.ne.0)then + int_prod_bessel=0.d0 + return + endif + + int_prod_bessel=crochet(l,gam)*freal + return + endif + + if(u.gt.6.d0)then + int_prod_bessel=int_prod_bessel_large(l,gam,n,m,a,b,arg) + return + endif + + if(a.ne.0.d0.and.b.ne.0.d0)then + + q=0 + intold=-1.d0 + int=0.d0 + done=.false. + + n_1 = n+n+1 + m_1 = m+m+1 + nlm = n+m+l + pi=dacos(-1.d0) + a_over_b_square = (a*a)/(b*b) + + ! First term of the sequence + + term_a =dble_fact(nlm-1) / (dble_fact(n_1)*dble_fact(m_1)) + expo=0.5d0*dfloat(nlm+1) + term_rap = term_a / (2.d0*gam)**expo + + s_0_0=term_rap*a**(n)*b**(m) + if(mod(nlm,2).eq.0)s_0_0=s_0_0*dsqrt(pi*.5d0) + + ! Initialize the first recurrence term for the q loop + s_q_0 = s_0_0 + + + mk = dble(m) + ! Loop over q for the convergence of the sequence + do while (.not.done) + + ! Init + s_q_k=s_q_0 + sum=s_q_0 + + if (q>300) then + stop 'pseudopot.f90 : q > 300' + endif + + qk = dble(q) + two_qkmp1 = 2.d0*(qk+mk)+1.d0 + do k=0,q-1 + if (s_q_k < 1.d-32) then + s_q_k = 0.d0 + exit + endif + s_q_k = two_qkmp1*qk*inverses(k)*s_q_k + sum=sum+s_q_k + two_qkmp1 = two_qkmp1-2.d0 + qk = qk-1.d0 + enddo + inverses(q) = a_over_b_square/(dble(q+n+q+n+3) * dble(q+1)) + + int=int+sum + + if(dabs(int-intold).lt.1d-15)then + done=.true. + else + + !Compute the s_q+1_0 + s_q_0=s_q_0*(q+q+nlm+1)*b*b/(dble(8*(m+q)+12)*(q+1)*gam) + + if(mod(n+m+l,2).eq.1)s_q_0=s_q_0*dsqrt(pi*.5d0) + ! Increment q + q=q+1 + intold=int + endif + + enddo + + int_prod_bessel=int*freal + + return + endif + + if(a.eq.0.d0.and.b.ne.0.d0)then + + int = int_prod_bessel_loc(l,gam,m,b) + int_prod_bessel=int*freal + return + endif + + if(a.ne.0.d0.and.b.eq.0.d0)then + + int = int_prod_bessel_loc(l,gam,n,a) + int_prod_bessel=int*freal + return + + endif + + stop 'pb in int_prod_bessel!!' +end + +double precision function int_prod_bessel_large(l,gam,n,m,a,b,arg) + implicit none + integer n,m,i,npts,l + double precision gam,a,b + double precision sum,x,bessel_mod,u,factor,arg + double precision xq(100),wq(100) + + u=(a+b)/(2.d0*dsqrt(gam)) + factor=dexp(u*u-arg)/dsqrt(gam) + +xq(1)= 5.38748089001123 +xq(2)= 4.60368244955074 +xq(3)= 3.94476404011563 +xq(4)= 3.34785456738322 +xq(5)= 2.78880605842813 +xq(6)= 2.25497400208928 +xq(7)= 1.73853771211659 +xq(8)= 1.23407621539532 +xq(9)= 0.737473728545394 +xq(10)= 0.245340708300901 +xq(11)=-0.245340708300901 +xq(12)=-0.737473728545394 +xq(13)=-1.23407621539532 +xq(14)=-1.73853771211659 +xq(15)=-2.25497400208928 +xq(16)=-2.78880605842813 +xq(17)=-3.34785456738322 +xq(18)=-3.94476404011563 +xq(19)=-4.60368244955074 +xq(20)=-5.38748089001123 +wq(1)= 2.229393645534151E-013 +wq(2)= 4.399340992273176E-010 +wq(3)= 1.086069370769280E-007 +wq(4)= 7.802556478532063E-006 +wq(5)= 2.283386360163528E-004 +wq(6)= 3.243773342237853E-003 +wq(7)= 2.481052088746362E-002 +wq(8)= 0.109017206020022 +wq(9)= 0.286675505362834 +wq(10)= 0.462243669600610 +wq(11)= 0.462243669600610 +wq(12)= 0.286675505362834 +wq(13)= 0.109017206020022 +wq(14)= 2.481052088746362E-002 +wq(15)= 3.243773342237853E-003 +wq(16)= 2.283386360163528E-004 +wq(17)= 7.802556478532063E-006 +wq(18)= 1.086069370769280E-007 +wq(19)= 4.399340992273176E-010 +wq(20)= 2.229393645534151E-013 + + npts=20 +! call gauher(xq,wq,npts) + + sum=0.d0 + do i=1,npts + x=(xq(i)+u)/dsqrt(gam) + sum=sum+wq(i)*x**l*bessel_mod(a*x,n)*bessel_mod(b*x,m)*dexp(-(a+b)*x) + enddo + int_prod_bessel_large=sum*factor +end + +!! Calculation of +!! +!! I= \int dx x**l *exp(-gam*x**2) M_n(ax) +!! +!! M_n(x) modified spherical bessel function +!! +double precision function int_prod_bessel_loc(l,gam,n,a) + implicit none + integer n,k,l,kcp + double precision gam,a + double precision int,intold,coef_nk,crochet,dble_fact, fact, pi, expo + double precision :: f_0, f_k + logical done + + pi=dacos(-1.d0) + intold=-1.d0 + done=.false. + int=0 + + ! Int f_0 + coef_nk=1.d0/dble_fact( n+n+1 ) + expo=0.5d0*dfloat(n+l+1) + crochet=dble_fact(n+l-1)/(gam+gam)**expo + if(mod(n+l,2).eq.0)crochet=crochet*dsqrt(0.5d0*pi) + + f_0 = coef_nk*a**n*crochet + + k=0 + + f_k=f_0 + do while (.not.done) + + int=int+f_k + +! f_k = f_k*(a**2*(2*(k+1)+n+l-1)) / (2*(k+1)*(2*(n+k+1)+1)*2*gam) + f_k = f_k*(a*a*dble(k+k+1+n+l)) / (dble((k+k+2)*(4*(n+k+1)+2))*gam) + + if(dabs(int-intold).lt.1d-15)then + done=.true. + else + k=k+1 + intold=int + endif + + enddo + + int_prod_bessel_loc=int +end + + double precision function int_prod_bessel_num(l,gam,n,m,a,b) + implicit none + integer n,m,l,i,npoints + double precision gam,a,b + double precision sum,dx,x,bessel_mod + sum=0.d0 + npoints=20000 + dx=30.d0/npoints + do i=1,npoints + x=(i-1)*dx+0.5d0*dx + sum=sum+dx*x**l*dexp(-gam*x**2)*bessel_mod(a*x,n)*bessel_mod(b*x,m) + enddo + int_prod_bessel_num=sum + end + + +! l,m : Y(l,m) parameters +! c(3) : pseudopotential center +! a(3) : Atomic Orbital center +! n_a(3) : Powers of x,y,z in the Atomic Orbital +! g_a : Atomic Orbital exponent +! r : Distance between the Atomic Orbital center and the considered point +double precision function ylm_orb(l,m,c,a,n_a,g_a,r) +implicit none +integer lmax_max +integer l,m +double precision a(3),g_a,c(3) +double precision prod,binom_func,accu,bigI,ylm,bessel_mod +double precision theta_AC0,phi_AC0,ac,ac2,factor,fourpi,arg,r,areal +integer ntotA,mu,k1,k2,k3,lambda +integer n_a(3) +double precision y, f1, f2 +double precision, allocatable :: array_coefs_A(:,:) + +ac2=(a(1)-c(1))**2+(a(2)-c(2))**2+(a(3)-c(3))**2 +ac=dsqrt(ac2) +arg=g_a*(ac2+r*r) +fourpi=4.d0*dacos(-1.d0) +factor=fourpi*dexp(-arg) +areal=2.d0*g_a*ac +ntotA=n_a(1)+n_a(2)+n_a(3) + + +if(ac.eq.0.d0)then + ylm_orb=dsqrt(fourpi)*r**ntotA*dexp(-g_a*r**2)*bigI(0,0,l,m,n_a(1),n_a(2),n_a(3)) + return +else + + theta_AC0=dacos( (a(3)-c(3))/ac ) + phi_AC0=datan2((a(2)-c(2))/ac,(a(1)-c(1))/ac) + + allocate (array_coefs_A(0:ntotA,3)) + do k1=0,n_a(1) + array_coefs_A(k1,1) = binom_func(n_a(1),k1)*(c(1)-a(1))**(n_a(1)-k1)*r**(k1) + enddo + do k2=0,n_a(2) + array_coefs_A(k2,2) = binom_func(n_a(2),k2)*(c(2)-a(2))**(n_a(2)-k2)*r**(k2) + enddo + do k3=0,n_a(3) + array_coefs_A(k3,3) = binom_func(n_a(3),k3)*(c(3)-a(3))**(n_a(3)-k3)*r**(k3) + enddo + + accu=0.d0 + do lambda=0,l+ntotA + do mu=-lambda,lambda + y = ylm(lambda,mu,theta_AC0,phi_AC0) + if (y == 0.d0) then + cycle + endif + do k3=0,n_a(3) + f1 = y*array_coefs_A(k3,3) + if (f1 == 0.d0) cycle + do k2=0,n_a(2) + f2 = f1*array_coefs_A(k2,2) + if (f2 == 0.d0) cycle + do k1=0,n_a(1) + prod=f2*array_coefs_A(k1,1)*bigI(lambda,mu,l,m,k1,k2,k3) + if (prod == 0.d0) then + cycle + endif + if (areal*r < 100.d0) then ! overflow! + accu=accu+prod*bessel_mod(areal*r,lambda) + endif + enddo + enddo + enddo + enddo + enddo + ylm_orb=factor*accu + deallocate (array_coefs_A) + return +endif + +end diff --git a/src/ao_cart_one_e_ints/spread_dipole_ao.irp.f b/src/ao_cart_one_e_ints/spread_dipole_ao.irp.f new file mode 100644 index 00000000..81e77668 --- /dev/null +++ b/src/ao_cart_one_e_ints/spread_dipole_ao.irp.f @@ -0,0 +1,376 @@ + BEGIN_PROVIDER [ double precision, ao_cart_spread_x, (ao_cart_num,ao_cart_num)] + &BEGIN_PROVIDER [ double precision, ao_cart_spread_y, (ao_cart_num,ao_cart_num)] + &BEGIN_PROVIDER [ double precision, ao_cart_spread_z, (ao_cart_num,ao_cart_num)] + BEGIN_DOC + ! * array of the integrals of ao_cart_i * x^2 ao_cart_j + ! + ! * array of the integrals of ao_cart_i * y^2 ao_cart_j + ! + ! * array of the integrals of ao_cart_i * z^2 ao_cart_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z + dim1=500 + lower_exp_val = 40.d0 + ao_cart_spread_x= 0.d0 + ao_cart_spread_y= 0.d0 + ao_cart_spread_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_cart_power,ao_cart_prim_num, & + !$OMP ao_cart_spread_x,ao_cart_spread_y,ao_cart_spread_z,ao_cart_num,ao_cart_coef_normalized_ordered_transp,ao_cart_nucl, & + !$OMP ao_cart_expo_ordered_transp,dim1,lower_exp_val) + do j=1,ao_cart_num + A_center(1) = nucl_coord( ao_cart_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_cart_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_cart_nucl(j), 3 ) + power_A(1) = ao_cart_power( j, 1 ) + power_A(2) = ao_cart_power( j, 2 ) + power_A(3) = ao_cart_power( j, 3 ) + do i= 1,ao_cart_num + B_center(1) = nucl_coord( ao_cart_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_cart_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_cart_nucl(i), 3 ) + power_B(1) = ao_cart_power( i, 1 ) + power_B(2) = ao_cart_power( i, 2 ) + power_B(3) = ao_cart_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(n,j) + do l = 1, ao_cart_prim_num(i) + c = ao_cart_coef_normalized_ordered_transp(n,j)*ao_cart_coef_normalized_ordered_transp(l,i) + beta = ao_cart_expo_ordered_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + call overlap_bourrin_spread(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) + accu_x += c*tmp*overlap_y*overlap_z + call overlap_bourrin_spread(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) + accu_y += c*tmp*overlap_x*overlap_z + call overlap_bourrin_spread(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) + accu_z += c*tmp*overlap_y*overlap_x + enddo + enddo + ao_cart_spread_x(i,j) = accu_x + ao_cart_spread_y(i,j) = accu_y + ao_cart_spread_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + + + BEGIN_PROVIDER [ double precision, ao_cart_dipole_x, (ao_cart_num,ao_cart_num)] + &BEGIN_PROVIDER [ double precision, ao_cart_dipole_y, (ao_cart_num,ao_cart_num)] + &BEGIN_PROVIDER [ double precision, ao_cart_dipole_z, (ao_cart_num,ao_cart_num)] + BEGIN_DOC + ! * array of the integrals of ao_cart_i * x ao_cart_j + ! + ! * array of the integrals of ao_cart_i * y ao_cart_j + ! + ! * array of the integrals of ao_cart_i * z ao_cart_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z,accu_x,accu_y,accu_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c + dim1=500 + lower_exp_val = 40.d0 + ao_cart_dipole_x= 0.d0 + ao_cart_dipole_y= 0.d0 + ao_cart_dipole_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_cart_power,ao_cart_prim_num, & + !$OMP ao_cart_dipole_x,ao_cart_dipole_y,ao_cart_dipole_z,ao_cart_num,ao_cart_coef_normalized_ordered_transp,ao_cart_nucl, & + !$OMP ao_cart_expo_ordered_transp,dim1,lower_exp_val) + do j=1,ao_cart_num + A_center(1) = nucl_coord( ao_cart_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_cart_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_cart_nucl(j), 3 ) + power_A(1) = ao_cart_power( j, 1 ) + power_A(2) = ao_cart_power( j, 2 ) + power_A(3) = ao_cart_power( j, 3 ) + do i= 1,ao_cart_num + B_center(1) = nucl_coord( ao_cart_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_cart_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_cart_nucl(i), 3 ) + power_B(1) = ao_cart_power( i, 1 ) + power_B(2) = ao_cart_power( i, 2 ) + power_B(3) = ao_cart_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(n,j) + do l = 1, ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(l,i) + c = ao_cart_coef_normalized_ordered_transp(l,i)*ao_cart_coef_normalized_ordered_transp(n,j) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + + call overlap_bourrin_dipole(A_center(1),B_center(1),alpha,beta,power_A(1),power_B(1),tmp,lower_exp_val,dx,dim1) + accu_x = accu_x + c*tmp*overlap_y*overlap_z + call overlap_bourrin_dipole(A_center(2),B_center(2),alpha,beta,power_A(2),power_B(2),tmp,lower_exp_val,dx,dim1) + accu_y = accu_y + c*tmp*overlap_x*overlap_z + call overlap_bourrin_dipole(A_center(3),B_center(3),alpha,beta,power_A(3),power_B(3),tmp,lower_exp_val,dx,dim1) + accu_z = accu_z + c*tmp*overlap_y*overlap_x + enddo + enddo + ao_cart_dipole_x(i,j) = accu_x + ao_cart_dipole_y(i,j) = accu_y + ao_cart_dipole_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + BEGIN_PROVIDER [ double precision, ao_cart_deriv_1_x, (ao_cart_num,ao_cart_num)] + &BEGIN_PROVIDER [ double precision, ao_cart_deriv_1_y, (ao_cart_num,ao_cart_num)] + &BEGIN_PROVIDER [ double precision, ao_cart_deriv_1_z, (ao_cart_num,ao_cart_num)] + BEGIN_DOC + ! * array of the integrals of ao_cart_i * d/dx ao_cart_j + ! + ! * array of the integrals of ao_cart_i * d/dy ao_cart_j + ! + ! * array of the integrals of ao_cart_i * d/dz ao_cart_j + END_DOC + implicit none + integer :: i,j,n,l + double precision :: f, tmp + integer :: dim1 + double precision :: overlap, overlap_x, overlap_y, overlap_z + double precision :: alpha, beta + double precision :: A_center(3), B_center(3) + integer :: power_A(3), power_B(3) + double precision :: lower_exp_val, dx, c,accu_x,accu_y,accu_z + integer :: i_component + dim1=500 + lower_exp_val = 40.d0 + ao_cart_deriv_1_x= 0.d0 + ao_cart_deriv_1_y= 0.d0 + ao_cart_deriv_1_z= 0.d0 + !$OMP PARALLEL DO SCHEDULE(GUIDED) & + !$OMP DEFAULT(NONE) & + !$OMP PRIVATE(A_center,B_center,power_A,power_B,& + !$OMP overlap_x,overlap_y, overlap_z, overlap, & + !$OMP alpha, beta,i,j,dx,tmp,c,i_component,accu_x,accu_y,accu_z) & + !$OMP SHARED(nucl_coord,ao_cart_power,ao_cart_prim_num, & + !$OMP ao_cart_deriv_1_x,ao_cart_deriv_1_y,ao_cart_deriv_1_z,ao_cart_num,ao_cart_coef_normalized_ordered_transp,ao_cart_nucl, & + !$OMP ao_cart_expo_ordered_transp,dim1,lower_exp_val) + do j=1,ao_cart_num + A_center(1) = nucl_coord( ao_cart_nucl(j), 1 ) + A_center(2) = nucl_coord( ao_cart_nucl(j), 2 ) + A_center(3) = nucl_coord( ao_cart_nucl(j), 3 ) + power_A(1) = ao_cart_power( j, 1 ) + power_A(2) = ao_cart_power( j, 2 ) + power_A(3) = ao_cart_power( j, 3 ) + do i= 1,ao_cart_num + B_center(1) = nucl_coord( ao_cart_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_cart_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_cart_nucl(i), 3 ) + power_B(1) = ao_cart_power( i, 1 ) + power_B(2) = ao_cart_power( i, 2 ) + power_B(3) = ao_cart_power( i, 3 ) + accu_x = 0.d0 + accu_y = 0.d0 + accu_z = 0.d0 + do n = 1,ao_cart_prim_num(j) + alpha = ao_cart_expo_ordered_transp(n,j) + do l = 1, ao_cart_prim_num(i) + beta = ao_cart_expo_ordered_transp(l,i) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + c = ao_cart_coef_normalized_ordered_transp(l,i) * ao_cart_coef_normalized_ordered_transp(n,j) + i_component = 1 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_x += c*(tmp*overlap_y*overlap_z) + i_component = 2 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_y += c*(tmp*overlap_x*overlap_z) + i_component = 3 + call overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,tmp,dim1) + accu_z += c*(tmp*overlap_y*overlap_x) + enddo + enddo + ao_cart_deriv_1_x(i,j) = accu_x + ao_cart_deriv_1_y(i,j) = accu_y + ao_cart_deriv_1_z(i,j) = accu_z + enddo + enddo + !$OMP END PARALLEL DO + END_PROVIDER + + + + + subroutine overlap_bourrin_spread(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) + BEGIN_DOC +! Computes the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x^2 ] +! needed for the dipole and those things + END_DOC + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + if(factor.lt.0.000001d0)then +! print*,'factor = ',factor + dx = 0.d0 + overlap_x = 0.d0 + return + endif + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) * x * x + enddo + overlap_x *= factor * dx + + end + + + subroutine overlap_bourrin_dipole(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) * x ] +! needed for the dipole and those things + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho + double precision :: P_center + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + if(power_B == 0 .and. power_A ==0)then + double precision :: F_integral + overlap_x = P_center * F_integral(0,p) * factor + dx = 0.d0 + return + endif + double precision :: pouet_timy + + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) * x + enddo + overlap_x *= factor * dx + + end + + subroutine overlap_bourrin_deriv_x(i_component,A_center,B_center,alpha,beta,power_A,power_B,dx,lower_exp_val,overlap_x,nx) + implicit none + integer :: i,j,k,l + integer,intent(in) :: power_A(3),power_B(3),i_component + double precision,intent(in) :: A_center(3), B_center(3),alpha,beta,lower_exp_val + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: overlap_first, overlap_second +! computes : = (a_x_i - 2 alpha ) + + call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)-1,power_B(i_component),overlap_first,lower_exp_val,dx,nx) + call overlap_bourrin_x(A_center(i_component),B_center(i_component),alpha,beta,power_A(i_component)+1,power_B(i_component),overlap_second,lower_exp_val,dx,nx) + overlap_x = (power_A(i_component) * overlap_first - 2.d0 * alpha * overlap_second) + end + + subroutine overlap_bourrin_x(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,lower_exp_val,dx,nx) + implicit none +! compute the following integral : +! int [-infty ; +infty] of [(x-A_center)^(power_A) * (x-B_center)^power_B * exp(-alpha(x-A_center)^2) * exp(-beta(x-B_center)^2) ] + integer :: i,j,k,l + integer,intent(in) :: power_A,power_B + double precision, intent(in) :: lower_exp_val + double precision,intent(in) :: A_center, B_center,alpha,beta + double precision, intent(out) :: overlap_x,dx + integer, intent(in) :: nx + double precision :: x_min,x_max,domain,x,factor,dist,p,p_inv,rho + double precision :: P_center,pouet_timy + if(power_A.lt.0.or.power_B.lt.0)then + overlap_x = 0.d0 + dx = 0.d0 + return + endif + p = alpha + beta + p_inv= 1.d0/p + rho = alpha * beta * p_inv + dist = (A_center - B_center)*(A_center - B_center) + P_center = (alpha * A_center + beta * B_center) * p_inv + factor = dexp(-rho * dist) + if(factor.lt.0.000001d0)then + dx = 0.d0 + overlap_x = 0.d0 + return + endif + + pouet_timy = dsqrt(lower_exp_val/p) + x_min = P_center - pouet_timy + x_max = P_center + pouet_timy + domain = x_max-x_min + dx = domain/dble(nx) + overlap_x = 0.d0 + x = x_min + do i = 1, nx + x += dx + overlap_x += (x-A_center)**(power_A) * (x-B_center)**(power_B) * dexp(-p * (x-P_center)*(x-P_center)) + enddo + overlap_x *= factor * dx + end +