PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-10-21 07:18:19 +02:00
Code Releases Activity

Merge pull request #354 from AbdAmmar/dev-stable

Browse Source 
Dev stable
This commit is contained in:
AbdAmmar 2024-10-18 20:05:46 +02:00 committed by GitHub
commit 340ed7a240
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
12 changed files with 1216 additions and 749 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/ao_basis/EZFIO.cfg
View File

@ -73,7 +73,7 @@ doc: If true, use cgtos for AO integrals
interface: ezfio interface: ezfio
default: False default: False
[ao_expo_im_cgtos] [ao_expo_im]
type: double precision type: double precision
doc: imag part for Exponents for each primitive of each cGTOs |AO| doc: imag part for Exponents for each primitive of each cGTOs |AO|
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max) size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
@ -82,12 +82,12 @@ interface: ezfio, provider
[ao_expo_pw] [ao_expo_pw]
type: double precision type: double precision
doc: plane wave part for each primitive GTOs |AO| doc: plane wave part for each primitive GTOs |AO|
size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max) size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider interface: ezfio, provider
[ao_expo_phase] [ao_expo_phase]
type: double precision type: double precision
doc: phase shift for each primitive GTOs |AO| doc: phase shift for each primitive GTOs |AO|
size: (4,ao_basis.ao_num,ao_basis.ao_prim_num_max) size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider interface: ezfio, provider

2
src/ao_basis/README.rst
View File

@ -40,6 +40,6 @@ Complex Gaussian-Type Orbitals (cGTOs) are also supported:
\chi_i(\mathbf{r}) = x^a y^b z^c \sum_k c_{ki} \left( e^{-\alpha_{ki} \mathbf{r}^2 - \imath \mathbf{k}_{ki} \cdot \mathbf{r} - \imath \phi_{ki}} + \text{C.C.} \right) \chi_i(\mathbf{r}) = x^a y^b z^c \sum_k c_{ki} \left( e^{-\alpha_{ki} \mathbf{r}^2 - \imath \mathbf{k}_{ki} \cdot \mathbf{r} - \imath \phi_{ki}} + \text{C.C.} \right)
where: where:
- :math:`\alpha \in \mathbb{C}` and :math:`\Re(\alpha) > 0` (specified by ``ao_expo`` and ``ao_expo_im_cgtos``), - :math:`\alpha \in \mathbb{C}` and :math:`\Re(\alpha) > 0` (specified by ``ao_expo`` and ``ao_expo_im``),
- :math:`\mathbf{k} = (k_x, k_y, k_z) \in \mathbb{R}^3` (specified by ``ao_expo_pw``), - :math:`\mathbf{k} = (k_x, k_y, k_z) \in \mathbb{R}^3` (specified by ``ao_expo_pw``),
- :math:`\phi = \phi_x + \phi_y + \phi_z \in \mathbb{R}` (specified by ``ao_expo_phase``). - :math:`\phi = \phi_x + \phi_y + \phi_z \in \mathbb{R}` (specified by ``ao_expo_phase``).

4
src/ao_basis/aos.irp.f
View File

@ -75,10 +75,6 @@ END_PROVIDER
enddo enddo
endif endif
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
! Normalization of the contracted basis functions ! Normalization of the contracted basis functions
if (ao_normalized) then if (ao_normalized) then
norm = 0.d0 norm = 0.d0

76
src/ao_one_e_ints/aos_cgtos.irp.f
View File

@ -4,6 +4,7 @@
BEGIN_PROVIDER [double precision, ao_coef_cgtos_norm_ord_transp, (ao_prim_num_max, ao_num)] BEGIN_PROVIDER [double precision, ao_coef_cgtos_norm_ord_transp, (ao_prim_num_max, ao_num)]
implicit none implicit none
integer :: i, j integer :: i, j
do j = 1, ao_num do j = 1, ao_num
@ -17,25 +18,22 @@ END_PROVIDER
! --- ! ---
BEGIN_PROVIDER [complex*16, ao_expo_cgtos_ord_transp, (ao_prim_num_max, ao_num)] BEGIN_PROVIDER [complex*16, ao_expo_cgtos_ord_transp, (ao_prim_num_max, ao_num)]
&BEGIN_PROVIDER [complex*16, ao_expo_pw_ord_transp, (4, ao_prim_num_max, ao_num)] &BEGIN_PROVIDER [double precision, ao_expo_pw_ord_transp, (4, ao_prim_num_max, ao_num)]
&BEGIN_PROVIDER [complex*16, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)] &BEGIN_PROVIDER [double precision, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)]
implicit none implicit none
integer :: i, j, m integer :: i, j, m
do j = 1, ao_num do j = 1, ao_num
do i = 1, ao_prim_num_max do i = 1, ao_prim_num_max
ao_expo_cgtos_ord_transp(i,j) = ao_expo_cgtos_ord(j,i) ao_expo_cgtos_ord_transp(i,j) = ao_expo_cgtos_ord(j,i)
do m = 1, 3
do m = 1, 4
ao_expo_pw_ord_transp(m,i,j) = ao_expo_pw_ord(m,j,i) ao_expo_pw_ord_transp(m,i,j) = ao_expo_pw_ord(m,j,i)
ao_expo_phase_ord_transp(m,i,j) = ao_expo_phase_ord(m,j,i) ao_expo_phase_ord_transp(m,i,j) = ao_expo_phase_ord(m,j,i)
enddo enddo
ao_expo_pw_ord_transp(4,i,j) = ao_expo_pw_ord_transp(1,i,j) * ao_expo_pw_ord_transp(1,i,j) &
+ ao_expo_pw_ord_transp(2,i,j) * ao_expo_pw_ord_transp(2,i,j) &
+ ao_expo_pw_ord_transp(3,i,j) * ao_expo_pw_ord_transp(3,i,j)
ao_expo_phase_ord_transp(4,i,j) = ao_expo_phase_ord_transp(1,j,i) &
+ ao_expo_phase_ord_transp(2,j,i) &
+ ao_expo_phase_ord_transp(3,j,i)
enddo enddo
enddo enddo
@ -50,16 +48,12 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
integer :: i, j, ii, m, powA(3), nz integer :: i, j, ii, m, powA(3), nz
double precision :: norm double precision :: norm
double precision :: kA2, phiA double precision :: kA2, phiA
complex*16 :: expo, expo_inv, C_A(3) complex*16 :: expo, expo_inv, C_Ae(3), C_Ap(3)
complex*16 :: overlap_x, overlap_y, overlap_z complex*16 :: overlap_x, overlap_y, overlap_z
complex*16 :: integ1, integ2, C1, C2 complex*16 :: integ1, integ2, C1, C2
nz = 100 nz = 100
C_A(1) = (0.d0, 0.d0)
C_A(2) = (0.d0, 0.d0)
C_A(3) = (0.d0, 0.d0)
ao_coef_norm_cgtos = 0.d0 ao_coef_norm_cgtos = 0.d0
do i = 1, ao_num do i = 1, ao_num
@ -69,27 +63,34 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
powA(2) = ao_power(i,2) powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3) powA(3) = ao_power(i,3)
! Normalization of the primitives
if(primitives_normalized) then if(primitives_normalized) then
! Normalization of the primitives
do j = 1, ao_prim_num(i) do j = 1, ao_prim_num(i)
expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im_cgtos(i,j) expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im(i,j)
expo_inv = (1.d0, 0.d0) / expo expo_inv = (1.d0, 0.d0) / expo
do m = 1, 3 do m = 1, 3
C_A(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j) C_Ap(m) = nucl_coord(ii,m)
C_Ae(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j)
enddo enddo
phiA = ao_expo_phase(4,i,j) phiA = ao_expo_phase(1,i,j) + ao_expo_phase(2,i,j) + ao_expo_phase(3,i,j)
KA2 = ao_expo_pw(4,i,j) KA2 = ao_expo_pw(1,i,j) * ao_expo_pw(1,i,j) &
+ ao_expo_pw(2,i,j) * ao_expo_pw(2,i,j) &
+ ao_expo_pw(3,i,j) * ao_expo_pw(3,i,j)
C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2) C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2)
C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2) C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2)
call overlap_cgaussian_xyz(C_A, C_A, expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz) call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, &
call overlap_cgaussian_xyz(conjg(C_A), C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz) C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz)
call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, &
conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz)
norm = 2.d0 * real(C1 * integ1 + C2 * integ2) norm = 2.d0 * real(C1 * integ1 + C2 * integ2)
!ao_coef_norm_cgtos(i,j) = 1.d0 / dsqrt(norm)
ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm) ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
enddo enddo
@ -99,7 +100,7 @@ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
ao_coef_norm_cgtos(i,j) = ao_coef(i,j) ao_coef_norm_cgtos(i,j) = ao_coef(i,j)
enddo enddo
endif endif ! primitives_normalized
enddo enddo
@ -126,12 +127,17 @@ END_PROVIDER
iorder(j) = j iorder(j) = j
d(j,1) = ao_expo(i,j) d(j,1) = ao_expo(i,j)
d(j,2) = ao_coef_norm_cgtos(i,j) d(j,2) = ao_coef_norm_cgtos(i,j)
d(j,3) = ao_expo_im_cgtos(i,j) d(j,3) = ao_expo_im(i,j)
do m = 1, 4 do m = 1, 3
d(j,3+m) = ao_expo_pw(m,i,j) d(j,3+m) = ao_expo_pw(m,i,j)
enddo
d(j,7) = d(j,4) * d(j,4) + d(j,5) * d(j,5) + d(j,6) * d(j,6)
do m = 1, 3
d(j,7+m) = ao_expo_phase(m,i,j) d(j,7+m) = ao_expo_phase(m,i,j)
enddo enddo
d(j,11) = d(j,8) + d(j,9) + d(j,10)
enddo enddo
call dsort(d(1,1), iorder, ao_prim_num(i)) call dsort(d(1,1), iorder, ao_prim_num(i))
@ -165,8 +171,8 @@ END_PROVIDER
double precision :: c, overlap, overlap_x, overlap_y, overlap_z double precision :: c, overlap, overlap_x, overlap_y, overlap_z
double precision :: KA2(3), phiA(3) double precision :: KA2(3), phiA(3)
double precision :: KB2(3), phiB(3) double precision :: KB2(3), phiB(3)
complex*16 :: alpha, alpha_inv, A_center(3) complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3)
complex*16 :: beta, beta_inv, B_center(3) complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3)
complex*16 :: C1(1:4), C2(1:4) complex*16 :: C1(1:4), C2(1:4)
complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1 complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1
complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2 complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2
@ -181,8 +187,8 @@ END_PROVIDER
!$OMP PARALLEL DO SCHEDULE(GUIDED) & !$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) & !$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, &
!$OMP alpha, alpha_inv, A_center, power_A, KA2, phiA, & !$OMP alpha, alpha_inv, Ae_center, Ap_center, power_A, KA2, phiA, &
!$OMP beta, beta_inv, B_center, power_B, KB2, phiB, & !$OMP beta, beta_inv, Be_center, Bp_center, power_B, KB2, phiB, &
!$OMP overlap_x , overlap_y , overlap_z , overlap, & !$OMP overlap_x , overlap_y , overlap_z , overlap, &
!$OMP overlap_x1, overlap_y1, overlap_z1, overlap1, & !$OMP overlap_x1, overlap_y1, overlap_z1, overlap1, &
!$OMP overlap_x2, overlap_y2, overlap_z2, overlap2) & !$OMP overlap_x2, overlap_y2, overlap_z2, overlap2) &
@ -212,7 +218,8 @@ END_PROVIDER
alpha_inv = (1.d0, 0.d0) / alpha alpha_inv = (1.d0, 0.d0) / alpha
do m = 1, 3 do m = 1, 3
phiA(m) = ao_expo_phase_ord_transp(m,n,j) phiA(m) = ao_expo_phase_ord_transp(m,n,j)
A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_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_expo_pw_ord_transp(m,n,j)
KA2(m) = ao_expo_pw_ord_transp(m,n,j) * ao_expo_pw_ord_transp(m,n,j) KA2(m) = ao_expo_pw_ord_transp(m,n,j) * ao_expo_pw_ord_transp(m,n,j)
enddo enddo
@ -222,7 +229,8 @@ END_PROVIDER
beta_inv = (1.d0, 0.d0) / beta beta_inv = (1.d0, 0.d0) / beta
do m = 1, 3 do m = 1, 3
phiB(m) = ao_expo_phase_ord_transp(m,l,i) phiB(m) = ao_expo_phase_ord_transp(m,l,i)
B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_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_expo_pw_ord_transp(m,l,i)
KB2(m) = ao_expo_pw_ord_transp(m,l,i) * ao_expo_pw_ord_transp(m,l,i) KB2(m) = ao_expo_pw_ord_transp(m,l,i) * ao_expo_pw_ord_transp(m,l,i)
enddo enddo
@ -238,11 +246,11 @@ END_PROVIDER
C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3))) 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) C2(4) = C2(1) * C2(2) * C2(3)
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x1, overlap_y1, overlap_z1, overlap1, dim1) Ap_center, Bp_center, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1)
call overlap_cgaussian_xyz(conjg(A_center), B_center, conjg(alpha), beta, power_A, power_B, & call overlap_cgaussian_xyz(conjg(Ae_center), Be_center, conjg(alpha), beta, power_A, power_B, &
overlap_x2, overlap_y2, overlap_z2, overlap2, dim1) 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_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_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2)

50
src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f
View File

@ -17,8 +17,8 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)]
double precision :: c, Z, C_center(3) double precision :: c, Z, C_center(3)
double precision :: phiA, KA2 double precision :: phiA, KA2
double precision :: phiB, KB2 double precision :: phiB, KB2
complex*16 :: alpha, alpha_inv, A_center(3) complex*16 :: alpha, alpha_inv, Ae_center(3), Ap_center(3)
complex*16 :: beta, beta_inv, B_center(3) complex*16 :: beta, beta_inv, Be_center(3), Bp_center(3)
complex*16 :: C1, C2, I1, I2 complex*16 :: C1, C2, I1, I2
complex*16 :: NAI_pol_mult_cgtos complex*16 :: NAI_pol_mult_cgtos
@ -27,9 +27,9 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)]
!$OMP PARALLEL & !$OMP PARALLEL &
!$OMP DEFAULT (NONE) & !$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, & !$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, C1, C2, I1, I2, &
!$OMP alpha, alpha_inv, A_center, phiA, KA2, power_A, C1, I1, & !$OMP alpha, alpha_inv, Ae_center, Ap_center, phiA, KA2, power_A, &
!$OMP beta, beta_inv, B_center, phiB, KB2, power_B, C2, I2) & !$OMP beta, beta_inv, Be_center, Bp_center, phiB, KB2, power_B) &
!$OMP SHARED (ao_num, ao_prim_num, ao_nucl, nucl_coord, & !$OMP SHARED (ao_num, ao_prim_num, ao_nucl, nucl_coord, &
!$OMP ao_power, nucl_num, nucl_charge, n_pt_max_integrals, & !$OMP ao_power, nucl_num, nucl_charge, n_pt_max_integrals, &
!$OMP ao_expo_cgtos_ord_transp, ao_coef_cgtos_norm_ord_transp, & !$OMP ao_expo_cgtos_ord_transp, ao_coef_cgtos_norm_ord_transp, &
@ -51,9 +51,9 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)]
alpha = ao_expo_cgtos_ord_transp(n,j) alpha = ao_expo_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha alpha_inv = (1.d0, 0.d0) / alpha
do m = 1, 3 do m = 1, 3
A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_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_expo_pw_ord_transp(m,n,j)
enddo enddo
phiA = ao_expo_phase_ord_transp(4,n,j) phiA = ao_expo_phase_ord_transp(4,n,j)
KA2 = ao_expo_pw_ord_transp(4,n,j) KA2 = ao_expo_pw_ord_transp(4,n,j)
@ -62,9 +62,9 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)]
beta = ao_expo_cgtos_ord_transp(l,i) beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta beta_inv = (1.d0, 0.d0) / beta
do m = 1, 3 do m = 1, 3
B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_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_expo_pw_ord_transp(m,l,i)
enddo enddo
phiB = ao_expo_phase_ord_transp(4,l,i) phiB = ao_expo_phase_ord_transp(4,l,i)
KB2 = ao_expo_pw_ord_transp(4,l,i) KB2 = ao_expo_pw_ord_transp(4,l,i)
@ -79,9 +79,11 @@ BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)]
C_center(1:3) = nucl_coord(k,1:3) C_center(1:3) = nucl_coord(k,1:3)
I1 = NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_max_integrals) 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(A_center), B_center, power_A, power_B, conjg(alpha), beta, 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) c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2)
enddo enddo
@ -100,7 +102,8 @@ END_PROVIDER
! --- ! ---
complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in) 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 BEGIN_DOC
! !
@ -115,7 +118,8 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp
integer, intent(in) :: n_pt_in, power_A(3), power_B(3) integer, intent(in) :: n_pt_in, power_A(3), power_B(3)
double precision, intent(in) :: C_center(3) double precision, intent(in) :: C_center(3)
complex*16, intent(in) :: alpha, beta, A_center(3), B_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 integer :: i, n_pt, n_pt_out
double precision :: dist_AB, dist_AC double precision :: dist_AB, dist_AC
@ -124,20 +128,20 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp
complex*16 :: d(0:n_pt_in) complex*16 :: d(0:n_pt_in)
complex*16, external :: V_n_e_cgtos complex*16, external :: V_n_e_cgtos
complex*16, external :: crint_2 complex*16, external :: crint_sum
complex*16, external :: crint_sum_2 complex*16, external :: crint_1
dist_AB = 0.d0 dist_AB = 0.d0
dist_AC = 0.d0 dist_AC = 0.d0
do i = 1, 3 do i = 1, 3
dist_AB += abs(A_center(i) - B_center(i)) dist_AB += abs(Ae_center(i) - Be_center(i))
dist_AC += abs(A_center(i) - C_center(i) * (1.d0, 0.d0)) dist_AC += abs(Ae_center(i) - C_center(i) * (1.d0, 0.d0))
enddo enddo
if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13)) then if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13) .or. use_pw) then
continue continue
@ -157,8 +161,8 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp
dist = (0.d0, 0.d0) dist = (0.d0, 0.d0)
dist_integral = (0.d0, 0.d0) dist_integral = (0.d0, 0.d0)
do i = 1, 3 do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv P_center(i) = (alpha * Ae_center(i) + beta * Be_center(i)) * p_inv
dist += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i)) dist += (Ae_center(i) - Be_center(i)) * (Ae_center(i) - Be_center(i))
dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i)) dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
enddo enddo
@ -175,12 +179,12 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp
n_pt = 2 * ((power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3))) 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 if(n_pt == 0) then
NAI_pol_mult_cgtos = coeff * crint_2(0, const) NAI_pol_mult_cgtos = coeff * crint_1(0, const)
return return
endif endif
d(0:n_pt_in) = (0.d0, 0.d0) d(0:n_pt_in) = (0.d0, 0.d0)
call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, & 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) power_A, power_B, C_center, n_pt_in, d, n_pt_out)
if(n_pt_out < 0) then if(n_pt_out < 0) then
@ -188,7 +192,7 @@ complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alp
return return
endif endif
NAI_pol_mult_cgtos = coeff * crint_sum_2(n_pt_out, const, d) NAI_pol_mult_cgtos = coeff * crint_sum(n_pt_out, const, d)
return return
end end

228
src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f
View File

@ -10,13 +10,15 @@
double precision :: c, deriv_tmp double precision :: c, deriv_tmp
double precision :: KA2, phiA double precision :: KA2, phiA
double precision :: KB2, phiB double precision :: KB2, phiB
complex*16 :: alpha, alpha_inv, A_center(3), C1 double precision :: aa
complex*16 :: beta, beta_inv, B_center(3), C2 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_x, overlap_y, overlap_z, overlap
complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1 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_x0_2, overlap_y0_2, overlap_z0_2
complex*16 :: overlap_m2_1, overlap_p2_1 complex*16 :: overlap_m2_1, overlap_m1_1, overlap_p1_1, overlap_p2_1
complex*16 :: overlap_m2_2, overlap_p2_2 complex*16 :: overlap_m2_2, overlap_m1_2, overlap_p1_2, overlap_p2_2
complex*16 :: deriv_tmp_1, deriv_tmp_2 complex*16 :: deriv_tmp_1, deriv_tmp_2
@ -24,32 +26,35 @@
! -- Dummy call to provide everything ! -- Dummy call to provide everything
A_center(:) = (0.0d0, 0.d0) Ae_center(:) = (0.0d0, 0.d0)
B_center(:) = (1.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) alpha = (1.0d0, 0.d0)
beta = (0.1d0, 0.d0) beta = (0.1d0, 0.d0)
power_A = 1 power_A = 1
power_B = 0 power_B = 0
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) Ap_center, Bp_center, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
! --- ! ---
!$OMP PARALLEL DO SCHEDULE(GUIDED) & !$OMP PARALLEL &
!$OMP DEFAULT(NONE) & !$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, & !$OMP PRIVATE(i, j, m, n, l, ii, jj, c, aa, xa, C1, C2, &
!$OMP A_center, power_A, alpha, alpha_inv, KA2, phiA, & !$OMP Ae_center, power_A, alpha, alpha_inv, KA2, phiA, Ap_center, &
!$OMP B_center, power_B, beta, beta_inv, KB2, phiB, & !$OMP Be_center, power_B, beta, beta_inv, KB2, phiB, Bp_center, &
!$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, & !$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, &
!$OMP overlap_x, overlap_y, overlap_z, overlap, & !$OMP overlap_x, overlap_y, overlap_z, overlap, &
!$OMP overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2, & !$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_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, &
!$OMP overlap_y0_2, overlap_z0_2) & !$OMP overlap_y0_2, overlap_z0_2) &
!$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, & !$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, &
!$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, & !$OMP ao_coef_cgtos_norm_ord_transp, ao_expo_cgtos_ord_transp, &
!$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, & !$OMP ao_expo_pw_ord_transp, ao_expo_phase_ord_transp, &
!$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z) !$OMP ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z)
!$OMP DO SCHEDULE(GUIDED)
do j = 1, ao_num do j = 1, ao_num
jj = ao_nucl(j) jj = ao_nucl(j)
@ -73,7 +78,8 @@
alpha = ao_expo_cgtos_ord_transp(n,j) alpha = ao_expo_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha alpha_inv = (1.d0, 0.d0) / alpha
do m = 1, 3 do m = 1, 3
A_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * alpha_inv * ao_expo_pw_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_expo_pw_ord_transp(m,n,j)
enddo enddo
phiA = ao_expo_phase_ord_transp(4,n,j) phiA = ao_expo_phase_ord_transp(4,n,j)
KA2 = ao_expo_pw_ord_transp(4,n,j) KA2 = ao_expo_pw_ord_transp(4,n,j)
@ -83,7 +89,8 @@
beta = ao_expo_cgtos_ord_transp(l,i) beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta beta_inv = (1.d0, 0.d0) / beta
do m = 1, 3 do m = 1, 3
B_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * beta_inv * ao_expo_pw_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_expo_pw_ord_transp(m,l,i)
enddo enddo
phiB = ao_expo_phase_ord_transp(4,l,i) phiB = ao_expo_phase_ord_transp(4,l,i)
KB2 = ao_expo_pw_ord_transp(4,l,i) KB2 = ao_expo_pw_ord_transp(4,l,i)
@ -93,42 +100,62 @@
C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2)) 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)) C2 = zexp((0.d0, 1.d0) * (-phiA + phiB) - 0.25d0 * (alpha_inv * KA2 + conjg(beta_inv) * KB2))
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1) 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, &
call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & Ap_center, conjg(Bp_center), overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1)
overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1)
! --- ! ---
power_A(1) = power_A(1) - 2 power_A(1) = power_A(1) - 1
if(power_A(1) > -1) then if(power_A(1) > -1) then
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_m2_1, overlap_y, overlap_z, overlap, dim1) 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, &
call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & Ap_center, conjg(Bp_center), overlap_m1_2, overlap_y, overlap_z, overlap, dim1)
overlap_m2_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 else
overlap_m2_1 = (0.d0, 0.d0) overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0) overlap_m2_2 = (0.d0, 0.d0)
endif 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) + 4 power_A(1) = power_A(1) + 1
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_p2_1, overlap_y, overlap_z, overlap, dim1) 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, &
call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & Ap_center, conjg(Bp_center), overlap_p1_2, overlap_y, overlap_z, overlap, dim1)
overlap_p2_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 power_A(1) = power_A(1) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(1)) + 1.d0) * overlap_x0_1 & aa = dble(power_A(1))
+ dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_1 & xa = Ap_center(1) - Ae_center(1)
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_y0_1 * overlap_z0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(1)) + 1.d0) * overlap_x0_2 & deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_2 & + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_x0_1 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2 + 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) deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
@ -136,34 +163,55 @@
! --- ! ---
power_A(2) = power_A(2) - 2 power_A(2) = power_A(2) - 1
if(power_A(2) > -1) then if(power_A(2) > -1) then
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_m2_1, overlap_y, overlap, dim1) 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, &
call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & Ap_center, conjg(Bp_center), overlap_x, overlap_m1_2, overlap_z, overlap, dim1)
overlap_x, overlap_m2_2, overlap_y, 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 else
overlap_m2_1 = (0.d0, 0.d0) overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0) overlap_m2_2 = (0.d0, 0.d0)
endif 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) + 4 power_A(2) = power_A(2) + 1
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_p2_1, overlap_y, overlap, dim1) 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, &
call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & Ap_center, conjg(Bp_center), overlap_x, overlap_p1_2, overlap_z, overlap, dim1)
overlap_x, overlap_p2_2, overlap_y, 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 power_A(2) = power_A(2) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(2)) + 1.d0) * overlap_y0_1 & aa = dble(power_A(2))
+ dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_1 & xa = Ap_center(2) - Ae_center(2)
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_z0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(2)) + 1.d0) * overlap_y0_2 & deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_2 & + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_y0_1 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2 + 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) deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
@ -171,34 +219,55 @@
! --- ! ---
power_A(3) = power_A(3) - 2 power_A(3) = power_A(3) - 1
if(power_A(3) > -1) then if(power_A(3) > -1) then
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_y, overlap_m2_1, overlap, dim1) 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, &
call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_m1_2, overlap, dim1)
overlap_x, overlap_y, overlap_m2_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 else
overlap_m2_1 = (0.d0, 0.d0) overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0) overlap_m2_2 = (0.d0, 0.d0)
endif 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) + 4 power_A(3) = power_A(3) + 1
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & call overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_y, overlap_p2_1, overlap, dim1) 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, &
call overlap_cgaussian_xyz(A_center, conjg(B_center), alpha, conjg(beta), power_A, power_B, & Ap_center, conjg(Bp_center), overlap_x, overlap_y, overlap_p1_2, overlap, dim1)
overlap_x, overlap_y, overlap_p2_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 power_A(3) = power_A(3) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(3)) + 1.d0) * overlap_z0_1 & aa = dble(power_A(3))
+ dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_1 & xa = Ap_center(3) - Ae_center(3)
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_y0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * dble(power_A(3)) + 1.d0) * overlap_z0_2 & deriv_tmp_1 = aa * (aa - 1.d0) * overlap_m2_1 - 4.d0 * alpha * aa * xa * overlap_m1_1 &
+ dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_2 & + 4.d0 * alpha * (alpha * xa * xa - aa - 0.5d0) * overlap_z0_1 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2 + 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) deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
@ -210,7 +279,8 @@
enddo enddo
enddo enddo
enddo enddo
!$OMP END PARALLEL DO !$OMP END DO
!$OMP END PARALLEL
END_PROVIDER END_PROVIDER

233
src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f
View File

@ -21,10 +21,10 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
double precision :: KJ2, phiJ double precision :: KJ2, phiJ
double precision :: KK2, phiK double precision :: KK2, phiK
double precision :: KL2, phiL double precision :: KL2, phiL
complex*16 :: expo1, expo1_inv, I_center(3) complex*16 :: expo1, expo1_inv, Ie_center(3), Ip_center(3)
complex*16 :: expo2, expo2_inv, J_center(3) complex*16 :: expo2, expo2_inv, Je_center(3), Jp_center(3)
complex*16 :: expo3, expo3_inv, K_center(3) complex*16 :: expo3, expo3_inv, Ke_center(3), Kp_center(3)
complex*16 :: expo4, expo4_inv, L_center(3) complex*16 :: expo4, expo4_inv, Le_center(3), Lp_center(3)
complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv
complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv
complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv
@ -44,7 +44,8 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
ao_two_e_integral_cgtos = ao_2e_cgtos_schwartz_accel(i, j, k, l) ao_two_e_integral_cgtos = ao_2e_cgtos_schwartz_accel(i, j, k, l)
else return
endif
dim1 = n_pt_max_integrals dim1 = n_pt_max_integrals
@ -63,7 +64,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
ao_two_e_integral_cgtos = 0.d0 ao_two_e_integral_cgtos = 0.d0
if(ii /= jj .or. kk /= ll .or. jj /= kk) then if(use_pw .or. ii /= jj .or. kk /= ll .or. jj /= kk) then
do p = 1, ao_prim_num(i) do p = 1, ao_prim_num(i)
@ -71,7 +72,8 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
expo1 = ao_expo_cgtos_ord_transp(p,i) expo1 = ao_expo_cgtos_ord_transp(p,i)
expo1_inv = (1.d0, 0.d0) / expo1 expo1_inv = (1.d0, 0.d0) / expo1
do m = 1, 3 do m = 1, 3
I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) Ip_center(m) = nucl_coord(ii,m)
Ie_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i)
enddo enddo
phiI = ao_expo_phase_ord_transp(4,p,i) phiI = ao_expo_phase_ord_transp(4,p,i)
KI2 = ao_expo_pw_ord_transp(4,p,i) KI2 = ao_expo_pw_ord_transp(4,p,i)
@ -82,17 +84,20 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
expo2 = ao_expo_cgtos_ord_transp(q,j) expo2 = ao_expo_cgtos_ord_transp(q,j)
expo2_inv = (1.d0, 0.d0) / expo2 expo2_inv = (1.d0, 0.d0) / expo2
do m = 1, 3 do m = 1, 3
J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) Jp_center(m) = nucl_coord(jj,m)
Je_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j)
enddo enddo
phiJ = ao_expo_phase_ord_transp(4,q,j) phiJ = ao_expo_phase_ord_transp(4,q,j)
KJ2 = ao_expo_pw_ord_transp(4,q,j) KJ2 = ao_expo_pw_ord_transp(4,q,j)
call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, &
expo1, expo2, I_power, J_power, I_center, J_center, dim1) expo1, expo2, I_power, J_power, Ie_center, Je_center, Ip_center, Jp_center, dim1)
p1_inv = (1.d0, 0.d0) / pp1 p1_inv = (1.d0, 0.d0) / pp1
call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, &
conjg(expo1), expo2, I_power, J_power, conjg(I_center), J_center, dim1) conjg(expo1), expo2, I_power, J_power, conjg(Ie_center), Je_center, conjg(Ip_center), Jp_center, dim1)
p2_inv = (1.d0, 0.d0) / pp2 p2_inv = (1.d0, 0.d0) / pp2
do r = 1, ao_prim_num(k) do r = 1, ao_prim_num(k)
@ -101,7 +106,8 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
expo3 = ao_expo_cgtos_ord_transp(r,k) expo3 = ao_expo_cgtos_ord_transp(r,k)
expo3_inv = (1.d0, 0.d0) / expo3 expo3_inv = (1.d0, 0.d0) / expo3
do m = 1, 3 do m = 1, 3
K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) Kp_center(m) = nucl_coord(kk,m)
Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k)
enddo enddo
phiK = ao_expo_phase_ord_transp(4,r,k) phiK = ao_expo_phase_ord_transp(4,r,k)
KK2 = ao_expo_pw_ord_transp(4,r,k) KK2 = ao_expo_pw_ord_transp(4,r,k)
@ -112,17 +118,20 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
expo4 = ao_expo_cgtos_ord_transp(s,l) expo4 = ao_expo_cgtos_ord_transp(s,l)
expo4_inv = (1.d0, 0.d0) / expo4 expo4_inv = (1.d0, 0.d0) / expo4
do m = 1, 3 do m = 1, 3
L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) Lp_center(m) = nucl_coord(ll,m)
Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l)
enddo enddo
phiL = ao_expo_phase_ord_transp(4,s,l) phiL = ao_expo_phase_ord_transp(4,s,l)
KL2 = ao_expo_pw_ord_transp(4,s,l) KL2 = ao_expo_pw_ord_transp(4,s,l)
call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, &
expo3, expo4, K_power, L_power, K_center, L_center, dim1) expo3, expo4, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1)
q1_inv = (1.d0, 0.d0) / qq1 q1_inv = (1.d0, 0.d0) / qq1
call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, &
conjg(expo3), expo4, K_power, L_power, conjg(K_center), L_center, dim1) conjg(expo3), expo4, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1)
q2_inv = (1.d0, 0.d0) / qq2 q2_inv = (1.d0, 0.d0) / qq2
C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) &
@ -188,62 +197,34 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
coef1 = ao_coef_cgtos_norm_ord_transp(p,i) coef1 = ao_coef_cgtos_norm_ord_transp(p,i)
expo1 = ao_expo_cgtos_ord_transp(p,i) expo1 = ao_expo_cgtos_ord_transp(p,i)
expo1_inv = (1.d0, 0.d0) / expo1
do m = 1, 3
I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i)
enddo
phiI = ao_expo_phase_ord_transp(4,p,i) phiI = ao_expo_phase_ord_transp(4,p,i)
KI2 = ao_expo_pw_ord_transp(4,p,i)
do q = 1, ao_prim_num(j) do q = 1, ao_prim_num(j)
coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j)
expo2 = ao_expo_cgtos_ord_transp(q,j) expo2 = ao_expo_cgtos_ord_transp(q,j)
expo2_inv = (1.d0, 0.d0) / expo2
do m = 1, 3
J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j)
enddo
phiJ = ao_expo_phase_ord_transp(4,q,j) phiJ = ao_expo_phase_ord_transp(4,q,j)
KJ2 = ao_expo_pw_ord_transp(4,q,j)
do r = 1, ao_prim_num(k) do r = 1, ao_prim_num(k)
coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k)
expo3 = ao_expo_cgtos_ord_transp(r,k) expo3 = ao_expo_cgtos_ord_transp(r,k)
expo3_inv = (1.d0, 0.d0) / expo3
do m = 1, 3
K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k)
enddo
phiK = ao_expo_phase_ord_transp(4,r,k) phiK = ao_expo_phase_ord_transp(4,r,k)
KK2 = ao_expo_pw_ord_transp(4,r,k)
do s = 1, ao_prim_num(l) do s = 1, ao_prim_num(l)
coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l)
expo4 = ao_expo_cgtos_ord_transp(s,l) expo4 = ao_expo_cgtos_ord_transp(s,l)
expo4_inv = (1.d0, 0.d0) / expo4
do m = 1, 3
L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l)
enddo
phiL = ao_expo_phase_ord_transp(4,s,l) phiL = ao_expo_phase_ord_transp(4,s,l)
KL2 = ao_expo_pw_ord_transp(4,s,l)
C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL))
- 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL))
C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL))
- 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL))
C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL))
- 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL))
C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL))
- 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL))
C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) &
- 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2))
C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) &
- 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2))
C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) &
- 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2))
C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) &
- 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2))
int1 = ERI_cgtos(expo1, expo2, expo3, expo4, & int1 = ERI_cgtos(expo1, expo2, expo3, expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), & I_power(1), J_power(1), K_power(1), L_power(1), &
@ -294,8 +275,7 @@ double precision function ao_two_e_integral_cgtos(i, j, k, l)
enddo ! q enddo ! q
enddo ! p enddo ! p
endif ! same centers endif
endif ! do schwartz
end end
@ -321,10 +301,10 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
double precision :: KJ2, phiJ double precision :: KJ2, phiJ
double precision :: KK2, phiK double precision :: KK2, phiK
double precision :: KL2, phiL double precision :: KL2, phiL
complex*16 :: expo1, expo1_inv, I_center(3) complex*16 :: expo1, expo1_inv, Ie_center(3), Ip_center(3)
complex*16 :: expo2, expo2_inv, J_center(3) complex*16 :: expo2, expo2_inv, Je_center(3), Jp_center(3)
complex*16 :: expo3, expo3_inv, K_center(3) complex*16 :: expo3, expo3_inv, Ke_center(3), Kp_center(3)
complex*16 :: expo4, expo4_inv, L_center(3) complex*16 :: expo4, expo4_inv, Le_center(3), Lp_center(3)
complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv complex*16 :: P1_new(0:max_dim,3), P1_center(3), fact_p1, pp1, p1_inv
complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv complex*16 :: P2_new(0:max_dim,3), P2_center(3), fact_p2, pp2, p2_inv
complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv complex*16 :: Q1_new(0:max_dim,3), Q1_center(3), fact_q1, qq1, q1_inv
@ -362,7 +342,7 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k))) allocate(schwartz_kl(0:ao_prim_num(l),0:ao_prim_num(k)))
if(ii /= jj .or. kk /= ll .or. jj /= kk) then if(use_pw .or. ii /= jj .or. kk /= ll .or. jj /= kk) then
schwartz_kl(0,0) = 0.d0 schwartz_kl(0,0) = 0.d0
do r = 1, ao_prim_num(k) do r = 1, ao_prim_num(k)
@ -371,7 +351,8 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
expo1 = ao_expo_cgtos_ord_transp(r,k) expo1 = ao_expo_cgtos_ord_transp(r,k)
expo1_inv = (1.d0, 0.d0) / expo1 expo1_inv = (1.d0, 0.d0) / expo1
do m = 1, 3 do m = 1, 3
K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k) Kp_center(m) = nucl_coord(kk,m)
Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k)
enddo enddo
phiK = ao_expo_phase_ord_transp(4,r,k) phiK = ao_expo_phase_ord_transp(4,r,k)
KK2 = ao_expo_pw_ord_transp(4,r,k) KK2 = ao_expo_pw_ord_transp(4,r,k)
@ -383,17 +364,20 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
expo2 = ao_expo_cgtos_ord_transp(s,l) expo2 = ao_expo_cgtos_ord_transp(s,l)
expo2_inv = (1.d0, 0.d0) / expo2 expo2_inv = (1.d0, 0.d0) / expo2
do m = 1, 3 do m = 1, 3
L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l) Lp_center(m) = nucl_coord(ll,m)
Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l)
enddo enddo
phiL = ao_expo_phase_ord_transp(4,s,l) phiL = ao_expo_phase_ord_transp(4,s,l)
KL2 = ao_expo_pw_ord_transp(4,s,l) KL2 = ao_expo_pw_ord_transp(4,s,l)
call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, &
expo1, expo2, K_power, L_power, K_center, L_center, dim1) expo1, expo2, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1)
p1_inv = (1.d0, 0.d0) / pp1 p1_inv = (1.d0, 0.d0) / pp1
call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, &
conjg(expo1), expo2, K_power, L_power, conjg(K_center), L_center, dim1) conjg(expo1), expo2, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1)
p2_inv = (1.d0, 0.d0) / pp2 p2_inv = (1.d0, 0.d0) / pp2
C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL) - 0.5d0 * (expo1_inv * KK2 + expo2_inv * KL2)) C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL) - 0.5d0 * (expo1_inv * KK2 + expo2_inv * KL2))
@ -457,7 +441,8 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
expo1 = ao_expo_cgtos_ord_transp(p,i) expo1 = ao_expo_cgtos_ord_transp(p,i)
expo1_inv = (1.d0, 0.d0) / expo1 expo1_inv = (1.d0, 0.d0) / expo1
do m = 1, 3 do m = 1, 3
I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i) Ip_center(m) = nucl_coord(ii,m)
Ie_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i)
enddo enddo
phiI = ao_expo_phase_ord_transp(4,p,i) phiI = ao_expo_phase_ord_transp(4,p,i)
KI2 = ao_expo_pw_ord_transp(4,p,i) KI2 = ao_expo_pw_ord_transp(4,p,i)
@ -468,17 +453,19 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
expo2 = ao_expo_cgtos_ord_transp(q,j) expo2 = ao_expo_cgtos_ord_transp(q,j)
expo2_inv = (1.d0, 0.d0) / expo2 expo2_inv = (1.d0, 0.d0) / expo2
do m = 1, 3 do m = 1, 3
J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j) Jp_center(m) = nucl_coord(jj,m)
Je_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j)
enddo enddo
phiJ = ao_expo_phase_ord_transp(4,q,j) phiJ = ao_expo_phase_ord_transp(4,q,j)
KJ2 = ao_expo_pw_ord_transp(4,q,j) KJ2 = ao_expo_pw_ord_transp(4,q,j)
call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, & call give_explicit_cpoly_and_cgaussian(P1_new, P1_center, pp1, fact_p1, iorder_p1, &
expo1, expo2, I_power, J_power, I_center, J_center, dim1) expo1, expo2, I_power, J_power, Ie_center, Je_center, Ip_center, Jp_center, dim1)
p1_inv = (1.d0, 0.d0) / pp1 p1_inv = (1.d0, 0.d0) / pp1
call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, & call give_explicit_cpoly_and_cgaussian(P2_new, P2_center, pp2, fact_p2, iorder_p2, &
conjg(expo1), expo2, I_power, J_power, conjg(I_center), J_center, dim1) conjg(expo1), expo2, I_power, J_power, conjg(Ie_center), Je_center, conjg(Ip_center), Jp_center, dim1)
p2_inv = (1.d0, 0.d0) / pp2 p2_inv = (1.d0, 0.d0) / pp2
C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ) - 0.5d0 * (expo1_inv * KI2 + expo2_inv * KJ2)) C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ) - 0.5d0 * (expo1_inv * KI2 + expo2_inv * KJ2))
@ -538,7 +525,8 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
expo3 = ao_expo_cgtos_ord_transp(r,k) expo3 = ao_expo_cgtos_ord_transp(r,k)
expo3_inv = (1.d0, 0.d0) / expo3 expo3_inv = (1.d0, 0.d0) / expo3
do m = 1, 3 do m = 1, 3
K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k) Kp_center(m) = nucl_coord(kk,m)
Ke_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k)
enddo enddo
phiK = ao_expo_phase_ord_transp(4,r,k) phiK = ao_expo_phase_ord_transp(4,r,k)
KK2 = ao_expo_pw_ord_transp(4,r,k) KK2 = ao_expo_pw_ord_transp(4,r,k)
@ -550,17 +538,20 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
expo4 = ao_expo_cgtos_ord_transp(s,l) expo4 = ao_expo_cgtos_ord_transp(s,l)
expo4_inv = (1.d0, 0.d0) / expo4 expo4_inv = (1.d0, 0.d0) / expo4
do m = 1, 3 do m = 1, 3
L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l) Lp_center(m) = nucl_coord(ll,m)
Le_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l)
enddo enddo
phiL = ao_expo_phase_ord_transp(4,s,l) phiL = ao_expo_phase_ord_transp(4,s,l)
KL2 = ao_expo_pw_ord_transp(4,s,l) KL2 = ao_expo_pw_ord_transp(4,s,l)
call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, & call give_explicit_cpoly_and_cgaussian(Q1_new, Q1_center, qq1, fact_q1, iorder_q1, &
expo3, expo4, K_power, L_power, K_center, L_center, dim1) expo3, expo4, K_power, L_power, Ke_center, Le_center, Kp_center, Lp_center, dim1)
q1_inv = (1.d0, 0.d0) / qq1 q1_inv = (1.d0, 0.d0) / qq1
call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, & call give_explicit_cpoly_and_cgaussian(Q2_new, Q2_center, qq2, fact_q2, iorder_q2, &
conjg(expo3), expo4, K_power, L_power, conjg(K_center), L_center, dim1) conjg(expo3), expo4, K_power, L_power, conjg(Ke_center), Le_center, conjg(Kp_center), Lp_center, dim1)
q2_inv = (1.d0, 0.d0) / qq2 q2_inv = (1.d0, 0.d0) / qq2
C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) &
@ -627,31 +618,21 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
coef1 = ao_coef_cgtos_norm_ord_transp(r,k) * ao_coef_cgtos_norm_ord_transp(r,k) coef1 = ao_coef_cgtos_norm_ord_transp(r,k) * ao_coef_cgtos_norm_ord_transp(r,k)
expo1 = ao_expo_cgtos_ord_transp(r,k) expo1 = ao_expo_cgtos_ord_transp(r,k)
expo1_inv = (1.d0, 0.d0) / expo1
do m = 1, 3
K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,r,k)
enddo
phiK = ao_expo_phase_ord_transp(4,r,k) phiK = ao_expo_phase_ord_transp(4,r,k)
KK2 = ao_expo_pw_ord_transp(4,r,k)
schwartz_kl(0,r) = 0.d0 schwartz_kl(0,r) = 0.d0
do s = 1, ao_prim_num(l) do s = 1, ao_prim_num(l)
coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(s,l) * ao_coef_cgtos_norm_ord_transp(s,l) coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(s,l) * ao_coef_cgtos_norm_ord_transp(s,l)
expo2 = ao_expo_cgtos_ord_transp(s,l) expo2 = ao_expo_cgtos_ord_transp(s,l)
expo2_inv = (1.d0, 0.d0) / expo2
do m = 1, 3
L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,s,l)
enddo
phiL = ao_expo_phase_ord_transp(4,s,l) phiL = ao_expo_phase_ord_transp(4,s,l)
KL2 = ao_expo_pw_ord_transp(4,s,l)
C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL) - 0.5d0 * (expo1_inv * KK2 + expo2_inv * KL2)) C1 = zexp(-(0.d0, 2.d0) * (phiK + phiL))
C2 = zexp(-(0.d0, 2.d0) * phiL - 0.5d0 * (real(expo1_inv) * KK2 + expo2_inv * KL2)) C2 = zexp(-(0.d0, 2.d0) * phiL)
!C3 = C2 !C3 = C2
C4 = zexp((0.d0, 2.d0) * (phiK - phiL) - 0.5d0 * (conjg(expo1_inv) * KK2 + expo2_inv * KL2)) C4 = zexp((0.d0, 2.d0) * (phiK - phiL))
C5 = zexp(-(0.d0, 2.d0) * phiK - 0.5d0 * (expo1_inv * KK2 + real(expo2_inv) * KL2)) C5 = zexp(-(0.d0, 2.d0) * phiK)
C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KK2 + real(expo2_inv) * KL2)) C6 = (1.d0, 0.d0)
!C7 = C6 !C7 = C6
!C8 = conjg(C5) !C8 = conjg(C5)
@ -711,30 +692,20 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
coef1 = ao_coef_cgtos_norm_ord_transp(p,i) coef1 = ao_coef_cgtos_norm_ord_transp(p,i)
expo1 = ao_expo_cgtos_ord_transp(p,i) expo1 = ao_expo_cgtos_ord_transp(p,i)
expo1_inv = (1.d0, 0.d0) / expo1
do m = 1, 3
I_center(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo1_inv * ao_expo_pw_ord_transp(m,p,i)
enddo
phiI = ao_expo_phase_ord_transp(4,p,i) phiI = ao_expo_phase_ord_transp(4,p,i)
KI2 = ao_expo_pw_ord_transp(4,p,i)
do q = 1, ao_prim_num(j) do q = 1, ao_prim_num(j)
coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j) coef2 = coef1 * ao_coef_cgtos_norm_ord_transp(q,j)
expo2 = ao_expo_cgtos_ord_transp(q,j) expo2 = ao_expo_cgtos_ord_transp(q,j)
expo2_inv = (1.d0, 0.d0) / expo2
do m = 1, 3
J_center(m) = nucl_coord(jj,m) - (0.d0, 0.5d0) * expo2_inv * ao_expo_pw_ord_transp(m,q,j)
enddo
phiJ = ao_expo_phase_ord_transp(4,q,j) phiJ = ao_expo_phase_ord_transp(4,q,j)
KJ2 = ao_expo_pw_ord_transp(4,q,j)
C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ) - 0.5d0 * (expo1_inv * KI2 + expo2_inv * KJ2)) C1 = zexp(-(0.d0, 2.d0) * (phiI + phiJ))
C2 = zexp(-(0.d0, 2.d0) * phiJ - 0.5d0 * (real(expo1_inv) * KI2 + expo2_inv * KJ2)) C2 = zexp(-(0.d0, 2.d0) * phiJ)
!C3 = C2 !C3 = C2
C4 = zexp((0.d0, 2.d0) * (phiI - phiJ) - 0.5d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2)) C4 = zexp((0.d0, 2.d0) * (phiI - phiJ))
C5 = zexp(-(0.d0, 2.d0) * phiI - 0.5d0 * (expo1_inv * KI2 + real(expo2_inv) * KJ2)) C5 = zexp(-(0.d0, 2.d0) * phiI)
C6 = zexp(-(0.5d0, 0.d0) * (real(expo1_inv) * KI2 + real(expo2_inv) * KJ2)) C6 = (1.d0, 0.d0)
!C7 = C6 !C7 = C6
!C8 = conjg(C5) !C8 = conjg(C5)
@ -791,41 +762,23 @@ double precision function ao_2e_cgtos_schwartz_accel(i, j, k, l)
coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k) coef3 = coef2 * ao_coef_cgtos_norm_ord_transp(r,k)
expo3 = ao_expo_cgtos_ord_transp(r,k) expo3 = ao_expo_cgtos_ord_transp(r,k)
expo3_inv = (1.d0, 0.d0) / expo3
do m = 1, 3
K_center(m) = nucl_coord(kk,m) - (0.d0, 0.5d0) * expo3_inv * ao_expo_pw_ord_transp(m,r,k)
enddo
phiK = ao_expo_phase_ord_transp(4,r,k) phiK = ao_expo_phase_ord_transp(4,r,k)
KK2 = ao_expo_pw_ord_transp(4,r,k)
do s = 1, ao_prim_num(l) do s = 1, ao_prim_num(l)
if(schwartz_kl(s,r)*schwartz_ij < thr) cycle if(schwartz_kl(s,r)*schwartz_ij < thr) cycle
coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l) coef4 = coef3 * ao_coef_cgtos_norm_ord_transp(s,l)
expo4 = ao_expo_cgtos_ord_transp(s,l) expo4 = ao_expo_cgtos_ord_transp(s,l)
expo4_inv = (1.d0, 0.d0) / expo4
do m = 1, 3
L_center(m) = nucl_coord(ll,m) - (0.d0, 0.5d0) * expo4_inv * ao_expo_pw_ord_transp(m,s,l)
enddo
phiL = ao_expo_phase_ord_transp(4,s,l) phiL = ao_expo_phase_ord_transp(4,s,l)
KL2 = ao_expo_pw_ord_transp(4,s,l)
C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL) & C1 = zexp((0.d0, 1.d0) * (-phiI - phiJ - phiK - phiL))
- 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL))
C2 = zexp((0.d0, 1.d0) * (-phiI - phiJ + phiK - phiL) & C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL))
- 0.25d0 * (expo1_inv * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL))
C3 = zexp((0.d0, 1.d0) * ( phiI - phiJ - phiK - phiL) & C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL))
- 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + expo3_inv * KK2 + expo4_inv * KL2)) C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL))
C4 = zexp((0.d0, 1.d0) * ( phiI - phiJ + phiK - phiL) & C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL))
- 0.25d0 * (conjg(expo1_inv) * KI2 + expo2_inv * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2)) C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL))
C5 = zexp((0.d0, 1.d0) * (-phiI + phiJ - phiK - phiL) &
- 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2))
C6 = zexp((0.d0, 1.d0) * (-phiI + phiJ + phiK - phiL) &
- 0.25d0 * (expo1_inv * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2))
C7 = zexp((0.d0, 1.d0) * ( phiI + phiJ - phiK - phiL) &
- 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + expo3_inv * KK2 + expo4_inv * KL2))
C8 = zexp((0.d0, 1.d0) * ( phiI + phiJ + phiK - phiL) &
- 0.25d0 * (conjg(expo1_inv) * KI2 + conjg(expo2_inv) * KJ2 + conjg(expo3_inv) * KK2 + expo4_inv * KL2))
int1 = ERI_cgtos(expo1, expo2, expo3, expo4, & int1 = ERI_cgtos(expo1, expo2, expo3, expo4, &
I_power(1), J_power(1), K_power(1), L_power(1), & I_power(1), J_power(1), K_power(1), L_power(1), &
@ -937,7 +890,7 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_
complex*16 :: dx(0:max_dim), Ix_pol(0:max_dim), dy(0:max_dim), Iy_pol(0:max_dim), dz(0:max_dim), Iz_pol(0:max_dim) complex*16 :: dx(0:max_dim), Ix_pol(0:max_dim), dy(0:max_dim), Iy_pol(0:max_dim), dz(0:max_dim), Iz_pol(0:max_dim)
complex*16 :: d1(0:max_dim), d_poly(0:max_dim) complex*16 :: d1(0:max_dim), d_poly(0:max_dim)
complex*16 :: crint_sum_2 complex*16 :: crint_sum
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: dx, Ix_pol, dy, Iy_pol, dz, Iz_pol
@ -1076,7 +1029,7 @@ complex*16 function general_primitive_integral_cgtos(dim, P_new, P_center, fact_
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out) call multiply_cpoly(d_poly, n_pt_tmp, Iz_pol, n_Iz, d1, n_pt_out)
accu = crint_sum_2(n_pt_out, const, d1) accu = crint_sum(n_pt_out, const, d1)
general_primitive_integral_cgtos = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / sq_ppq general_primitive_integral_cgtos = fact_p * fact_q * accu * pi_5_2 * p_inv * q_inv / sq_ppq
@ -1108,10 +1061,10 @@ complex*16 function ERI_cgtos(alpha, beta, delta, gama, a_x, b_x, c_x, d_x, a_y,
ERI_cgtos = (0.d0, 0.d0) ERI_cgtos = (0.d0, 0.d0)
ASSERT (REAL(alpha) >= 0.d0) ASSERT (real(alpha) >= 0.d0)
ASSERT (REAL(beta ) >= 0.d0) ASSERT (real(beta ) >= 0.d0)
ASSERT (REAL(delta) >= 0.d0) ASSERT (real(delta) >= 0.d0)
ASSERT (REAL(gama ) >= 0.d0) ASSERT (real(gama ) >= 0.d0)
nx = a_x + b_x + c_x + d_x nx = a_x + b_x + c_x + d_x
if(iand(nx, 1) == 1) then if(iand(nx, 1) == 1) then
@ -1727,3 +1680,21 @@ end
! --- ! ---
BEGIN_PROVIDER [logical, use_pw]
implicit none
logical :: exist
use_pw = .false.
call ezfio_has_ao_basis_ao_expo_pw(exist)
if(exist) then
PROVIDE ao_expo_pw_ord_transp
if(maxval(dabs(ao_expo_pw_ord_transp(4,:,:))) .gt. 1d-15) use_pw = .true.
endif
END_PROVIDER
! ---

293
src/hartree_fock/deb_ao_2e_int.irp.f
View File

@ -3,10 +3,13 @@ program deb_ao_2e_int
implicit none implicit none
call check_ao_two_e_integral_cgtos() !call check_ao_two_e_integral_cgtos()
!call check_crint1() !call check_crint1()
!call check_crint2() !call check_crint2()
!call check_crint3() !call check_crint3()
!call check_crint4()
call check_crint5()
!call check_crint6()
end end
@ -344,3 +347,291 @@ end
! --- ! ---
subroutine check_crint4()
implicit none
integer :: i_test, n_test
integer :: i, seed_size, clock_time
double precision :: xr(1), x, shift
double precision :: yr(1), y
double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
double precision :: t1, t2, t_int1, t_int2
complex*16 :: rho
complex*16 :: int1, int2, int3
integer, allocatable :: seed(:)
call random_seed(size=seed_size)
allocate(seed(seed_size))
call system_clock(count=clock_time)
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
!seed = 123456789
call random_seed(put=seed)
t_int1 = 0.d0
t_int2 = 0.d0
n_test = 100
shift = 15.d0
acc_re = 0.d0
nrm_re = 0.d0
acc_im = 0.d0
nrm_im = 0.d0
do i_test = 1, n_test
call random_number(xr)
call random_number(yr)
x = 1.d0 * (2.d0 * shift * xr(1) - shift)
y = 1.d0 * (2.d0 * shift * yr(1) - shift)
rho = x + (0.d0, 1.d0) * y
call wall_time(t1)
call zboysfun00_1(rho, int1)
call wall_time(t2)
t_int1 = t_int1 + t2 - t1
call wall_time(t1)
call zboysfun00_2(rho, int2)
call wall_time(t2)
t_int2 = t_int2 + t2 - t1
dif_re = dabs(real(int1) - real(int2))
dif_im = dabs(aimag(int1) - aimag(int2))
if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then
print*, ' important error found: '
print*, " rho = ", x, y
print*, real(int1), real(int2), dif_re
print*, aimag(int1), aimag(int2), dif_im
call crint_quad_12(0, rho, 10000000, int3)
if(zabs(int1 - int3) .lt. zabs(int2 - int3)) then
print*, ' implementation 2 seems to be wrong'
else
print*, ' implementation 1 seems to be wrong'
print*, ' quad 10000000:', real(int3), aimag(int3)
call crint_quad_12(0, rho, 100000000, int3)
print*, ' quad 100000000:', real(int3), aimag(int3)
endif
!print*, ' quad:', real(int3), aimag(int3)
!stop
endif
if((real(int1) /= real(int1)) .or. (aimag(int1) /= aimag(int1)) .or. &
(real(int2) /= real(int2)) .or. (aimag(int2) /= aimag(int2)) ) then
cycle
else
acc_re += dif_re
acc_im += dif_im
nrm_re += dabs(real(int1))
nrm_im += dabs(aimag(int1))
endif
enddo
print*, "accuracy on real part (%):", 100.d0 * acc_re / (nrm_re + 1d-15)
print*, "accuracy on imag part (%):", 100.d0 * acc_im / (nrm_im + 1d-15)
print*, "zerf_1 wall time (sec) = ", t_int1
print*, "zerf_2 wall time (sec) = ", t_int2
deallocate(seed)
end
! ---
subroutine check_crint5()
implicit none
integer :: i_test, n_test
integer :: i, seed_size, clock_time
integer :: n
double precision :: xr(1), yr(1), nr(1), x, shift, y
double precision :: dif1_re, dif1_im, acc1_re, acc1_im
double precision :: dif2_re, dif2_im, acc2_re, acc2_im
double precision :: nrm_re, nrm_im
double precision :: t1, t2, t_int1, t_int2
complex*16 :: rho
complex*16 :: int1, int2, int_ref
integer, allocatable :: seed(:)
complex*16, external :: crint_1, crint_2
call random_seed(size=seed_size)
allocate(seed(seed_size))
call system_clock(count=clock_time)
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
!seed = 123456789
call random_seed(put=seed)
t_int1 = 0.d0
t_int2 = 0.d0
n_test = 100
acc1_re = 0.d0
acc1_im = 0.d0
acc2_re = 0.d0
acc2_im = 0.d0
nrm_re = 0.d0
nrm_im = 0.d0
do i_test = 1, n_test
call random_number(xr)
call random_number(yr)
call random_number(nr)
x = 1.d+1 * (30.d0 * xr(1) - 15.d0)
y = 1.d+1 * (30.d0 * yr(1) - 15.d0)
n = int(16.d0 * nr(1))
rho = x + (0.d0, 1.d0) * y
call wall_time(t1)
int1 = crint_1(n, rho)
call wall_time(t2)
t_int1 = t_int1 + t2 - t1
call wall_time(t1)
int2 = crint_2(n, rho)
call wall_time(t2)
t_int2 = t_int2 + t2 - t1
call crint_quad_12(n, rho, 10000000, int_ref)
dif1_re = dabs(real(int1) - real(int_ref))
dif1_im = dabs(aimag(int1) - aimag(int_ref))
dif2_re = dabs(real(int2) - real(int_ref))
dif2_im = dabs(aimag(int2) - aimag(int_ref))
if((dif2_re .gt. 1d-7) .or. (dif2_im .gt. 1d-7)) then
print*, ' important error found: '
print*, " n, rho = ", n, x, y
print*, real(int1), real(int2), real(int_ref)
print*, aimag(int1), aimag(int2), aimag(int_ref)
!stop
endif
acc1_re += dif1_re
acc1_im += dif1_im
acc2_re += dif2_re
acc2_im += dif2_im
nrm_re += dabs(real(int_ref))
nrm_im += dabs(aimag(int_ref))
enddo
print*, "accuracy on boys_1 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc1_im / (nrm_im + 1d-15)
print*, "accuracy on boys_2 (%):", 100.d0 * acc1_re / (nrm_re + 1d-15), 100.d0 * acc2_im / (nrm_im + 1d-15)
print*, "boys_1 wall time (sec) = ", t_int1
print*, "boys_2 wall time (sec) = ", t_int2
deallocate(seed)
end
! ---
subroutine check_crint6()
implicit none
integer :: i_test, n_test
integer :: i, seed_size, clock_time
integer :: n
double precision :: xr(1), yr(1), nr(1), x, shift, y
double precision :: dif_re, dif_im, acc_re, acc_im
double precision :: nrm_re, nrm_im
double precision :: t1, t2, t_int1, t_int2
complex*16 :: rho
complex*16 :: int1, int2, int3
integer, allocatable :: seed(:)
complex*16, external :: crint_1, crint_2
call random_seed(size=seed_size)
allocate(seed(seed_size))
call system_clock(count=clock_time)
seed = clock_time + 37 * (/ (i, i=0, seed_size-1) /)
!seed = 123456789
call random_seed(put=seed)
t_int1 = 0.d0
t_int2 = 0.d0
n_test = 100
acc_re = 0.d0
acc_im = 0.d0
nrm_re = 0.d0
nrm_im = 0.d0
do i_test = 1, n_test
call random_number(xr)
call random_number(yr)
call random_number(nr)
x = 1.d0 * (30.d0 * xr(1) - 15.d0)
y = 1.d0 * (30.d0 * yr(1) - 15.d0)
n = int(16.d0 * nr(1))
rho = x + (0.d0, 1.d0) * y
call wall_time(t1)
int1 = crint_1(n, rho)
call wall_time(t2)
t_int1 = t_int1 + t2 - t1
call wall_time(t1)
int2 = crint_2(n, rho)
call wall_time(t2)
t_int2 = t_int2 + t2 - t1
dif_re = dabs(real(int1) - real(int2))
dif_im = dabs(aimag(int1) - aimag(int2))
if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then
print*, ' important error found: '
print*, " n, rho = ", n, x, y
print*, real(int1), real(int2), dif_re
print*, aimag(int1), aimag(int2), dif_im
call crint_quad_12(n, rho, 100000000, int3)
print*, ' quad 100000000:', real(int3), aimag(int3)
!print*, ' quad 100000000:', dabs(real(int1) - real(int3)), dabs(aimag(int1) - aimag(int3))
!stop
endif
acc_re += dif_re
acc_im += dif_im
nrm_re += dabs(real(int1))
nrm_im += dabs(aimag(int1))
enddo
print*, "diff (%):", 100.d0 * acc_re / (nrm_re + 1d-15), 100.d0 * acc_im / (nrm_im + 1d-15)
print*, "boys_1 wall time (sec) = ", t_int1
print*, "boys_2 wall time (sec) = ", t_int2
deallocate(seed)
end
! ---

86
src/utils/cgtos_one_e.irp.f
View File

@ -1,11 +1,12 @@
! --- ! ---
complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A, power_B, dim) complex*16 function overlap_cgaussian_x(Ae_center, Be_center, alpha, beta, power_A, power_B, Ap_center, Bp_center, dim)
BEGIN_DOC BEGIN_DOC
! !
! \int_{-infty}^{+infty} (x-A_x)^ax (x-B_x)^bx exp(-alpha (x-A_x)^2) exp(- beta(x-B_X)^2) dx ! \int_{-infty}^{+infty} (x - Ap_x)^ax (x - Bp_x)^bx exp(-alpha (x - Ae_x)^2) exp(-beta (x - Be_X)^2) dx
!
! with complex arguments ! with complex arguments
! !
END_DOC END_DOC
@ -14,20 +15,19 @@ complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A
include 'constants.include.F' include 'constants.include.F'
integer, intent(in) :: dim, power_A, power_B integer, intent(in) :: dim, power_A, power_B
complex*16, intent(in) :: A_center, B_center, alpha, beta complex*16, intent(in) :: Ae_center, alpha, Ap_center
complex*16, intent(in) :: Be_center, beta, Bp_center
integer :: i, iorder_p integer :: i, iorder_p
double precision :: fact_p_mod
complex*16 :: P_new(0:max_dim), P_center, fact_p, p, inv_sq_p complex*16 :: P_new(0:max_dim), P_center, fact_p, p, inv_sq_p
complex*16 :: Fc_integral complex*16 :: Fc_integral
call give_explicit_cpoly_and_cgaussian_x( P_new, P_center, p, fact_p, iorder_p & call give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_p, iorder_p, &
, alpha, beta, power_A, power_B, A_center, B_center, dim) alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim)
fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p)) if(zabs(fact_p) .lt. 1.d-14) then
if(fact_p_mod .lt. 1.d-14) then
overlap_cgaussian_x = (0.d0, 0.d0) overlap_cgaussian_x = (0.d0, 0.d0)
return return
endif endif
@ -37,22 +37,25 @@ complex*16 function overlap_cgaussian_x(A_center, B_center, alpha, beta, power_A
overlap_cgaussian_x = (0.d0, 0.d0) overlap_cgaussian_x = (0.d0, 0.d0)
do i = 0, iorder_p do i = 0, iorder_p
overlap_cgaussian_x += P_new(i) * Fc_integral(i, inv_sq_p) overlap_cgaussian_x = overlap_cgaussian_x + P_new(i) * Fc_integral(i, inv_sq_p)
enddo enddo
overlap_cgaussian_x *= fact_p overlap_cgaussian_x = overlap_cgaussian_x * fact_p
end end
! --- ! ---
subroutine overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, & subroutine overlap_cgaussian_xyz(Ae_center, Be_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_y, overlap_z, overlap, dim ) Ap_center, Bp_center, overlap_x, overlap_y, overlap_z, overlap, dim)
BEGIN_DOC BEGIN_DOC
! !
! S_x = \int (x-A_x)^{a_x} exp(-\alpha(x-A_x)^2) (x-B_x)^{b_x} exp(-beta(x-B_x)^2) dx ! S_x = \int (x - Ap_x)^{a_x} exp(-\alpha (x - Ae_x)^2)
! (x - Bp_x)^{b_x} exp(-\beta (x - Be_x)^2) dx
!
! S = S_x S_y S_z ! S = S_x S_y S_z
!
! for complex arguments ! for complex arguments
! !
END_DOC END_DOC
@ -61,58 +64,75 @@ subroutine overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power
include 'constants.include.F' include 'constants.include.F'
integer, intent(in) :: dim, power_A(3), power_B(3) integer, intent(in) :: dim, power_A(3), power_B(3)
complex*16, intent(in) :: A_center(3), B_center(3), alpha, beta complex*16, intent(in) :: Ae_center(3), alpha, Ap_center(3)
complex*16, intent(in) :: Be_center(3), beta, Bp_center(3)
complex*16, intent(out) :: overlap_x, overlap_y, overlap_z, overlap complex*16, intent(out) :: overlap_x, overlap_y, overlap_z, overlap
integer :: i, nmax, iorder_p(3) integer :: i, nmax, iorder_p(3)
double precision :: fact_p_mod
complex*16 :: P_new(0:max_dim,3), P_center(3), fact_p, p, inv_sq_p complex*16 :: P_new(0:max_dim,3), P_center(3), fact_p, p, inv_sq_p
complex*16 :: F_integral_tab(0:max_dim) complex*16 :: F_integral_tab(0:max_dim)
complex*16 :: ab, arg
complex*16 :: Fc_integral complex*16, external :: Fc_integral
call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, alpha, beta, power_A, power_B, A_center, B_center, dim) call give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_p, iorder_p, &
alpha, beta, power_A, power_B, Ae_center, Be_center, Ap_center, Bp_center, dim)
fact_p_mod = dsqrt(real(fact_p)*real(fact_p) + aimag(fact_p)*aimag(fact_p)) if(zabs(fact_p) .lt. 1.d-14) then
if(fact_p_mod .lt. 1.d-14) then overlap_x = (0.d0, 0.d0)
overlap_x = (1.d-10, 0.d0) overlap_y = (0.d0, 0.d0)
overlap_y = (1.d-10, 0.d0) overlap_z = (0.d0, 0.d0)
overlap_z = (1.d-10, 0.d0) overlap = (0.d0, 0.d0)
overlap = (1.d-10, 0.d0)
return return
endif endif
nmax = maxval(iorder_p) nmax = maxval(iorder_p)
inv_sq_p = (1.d0, 0.d0) / zsqrt(p) inv_sq_p = (1.d0, 0.d0) / zsqrt(p)
do i = 0, nmax do i = 0, nmax
F_integral_tab(i) = Fc_integral(i, inv_sq_p) F_integral_tab(i) = Fc_integral(i, inv_sq_p)
enddo enddo
overlap_x = P_new(0,1) * F_integral_tab(0) ab = alpha * beta * inv_sq_p
overlap_y = P_new(0,2) * F_integral_tab(0)
overlap_z = P_new(0,3) * F_integral_tab(0)
arg = ab * (Ae_center(1) - Be_center(1)) &
* (Ae_center(1) - Be_center(1))
if(real(arg) > 40.d0) then
overlap_x = (0.d0, 0.d0)
else
overlap_x = P_new(0,1) * F_integral_tab(0)
do i = 1, iorder_p(1) do i = 1, iorder_p(1)
overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i) overlap_x = overlap_x + P_new(i,1) * F_integral_tab(i)
enddo enddo
call cgaussian_product_x(alpha, A_center(1), beta, B_center(1), fact_p, p, P_center(1)) overlap_x = overlap_x * zexp(-arg)
overlap_x *= fact_p endif
arg = ab * (Ae_center(2) - Be_center(2)) &
* (Ae_center(2) - Be_center(2))
if(real(arg) > 40.d0) then
overlap_y = (0.d0, 0.d0)
else
overlap_y = P_new(0,2) * F_integral_tab(0)
do i = 1, iorder_p(2) do i = 1, iorder_p(2)
overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i) overlap_y = overlap_y + P_new(i,2) * F_integral_tab(i)
enddo enddo
call cgaussian_product_x(alpha, A_center(2), beta, B_center(2), fact_p, p, P_center(2)) overlap_y = overlap_y * zexp(-arg)
overlap_y *= fact_p endif
arg = ab * (Ae_center(3) - Be_center(3)) &
* (Ae_center(3) - Be_center(3))
if(real(arg) > 40.d0) then
overlap_z = (0.d0, 0.d0)
else
overlap_z = P_new(0,3) * F_integral_tab(0)
do i = 1, iorder_p(3) do i = 1, iorder_p(3)
overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i) overlap_z = overlap_z + P_new(i,3) * F_integral_tab(i)
enddo enddo
call cgaussian_product_x(alpha, A_center(3), beta, B_center(3), fact_p, p, P_center(3)) overlap_z = overlap_z * zexp(-arg)
overlap_z *= fact_p endif
overlap = overlap_x * overlap_y * overlap_z overlap = overlap_x * overlap_y * overlap_z
return
end end
! --- ! ---

76
src/utils/cgtos_utils.irp.f
View File

@ -1,13 +1,19 @@
! --- ! ---
subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim) subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorder, &
alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim)
BEGIN_DOC BEGIN_DOC
!
! Transform the product of ! Transform the product of
! (x-x_A)^a (x-x_B)^b exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) !
! (x - x_Ap)^a (x - x_Bp)^b exp(-alpha (r - Ae)^2) exp(-beta (r - Be)^2)
!
! into ! into
!
! fact_k \sum_{i=0}^{iorder} (x - x_P)^i exp(-p (r - P)^2) ! fact_k \sum_{i=0}^{iorder} (x - x_P)^i exp(-p (r - P)^2)
!
END_DOC END_DOC
implicit none implicit none
@ -15,13 +21,13 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde
integer, intent(in) :: dim integer, intent(in) :: dim
integer, intent(in) :: a, b integer, intent(in) :: a, b
complex*16, intent(in) :: alpha, beta, A_center, B_center complex*16, intent(in) :: alpha, Ae_center, Ap_center
complex*16, intent(in) :: beta, Be_center, Bp_center
integer, intent(out) :: iorder integer, intent(out) :: iorder
complex*16, intent(out) :: p, P_center, fact_k complex*16, intent(out) :: p, P_center, fact_k
complex*16, intent(out) :: P_new(0:max_dim) complex*16, intent(out) :: P_new(0:max_dim)
integer :: n_new, i, j integer :: n_new, i, j
double precision :: tmp_mod
complex*16 :: P_a(0:max_dim), P_b(0:max_dim) complex*16 :: P_a(0:max_dim), P_b(0:max_dim)
complex*16 :: p_inv, ab, d_AB, tmp complex*16 :: p_inv, ab, d_AB, tmp
@ -35,13 +41,12 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde
! new center ! new center
p_inv = (1.d0, 0.d0) / p p_inv = (1.d0, 0.d0) / p
ab = alpha * beta ab = alpha * beta
P_center = (alpha * A_center + beta * B_center) * p_inv P_center = (alpha * Ae_center + beta * Be_center) * p_inv
! get the factor ! get the factor
d_AB = (A_center - B_center) * (A_center - B_center) d_AB = (Ae_center - Be_center) * (Ae_center - Be_center)
tmp = ab * p_inv * d_AB tmp = ab * p_inv * d_AB
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) if(zabs(tmp) .lt. 50.d0) then
if(tmp_mod .lt. 50.d0) then
fact_k = zexp(-tmp) fact_k = zexp(-tmp)
else else
fact_k = (0.d0, 0.d0) fact_k = (0.d0, 0.d0)
@ -49,7 +54,7 @@ subroutine give_explicit_cpoly_and_cgaussian_x(P_new, P_center, p, fact_k, iorde
! Recenter the polynomials P_a and P_b on P_center ! Recenter the polynomials P_a and P_b on P_center
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0), A_center, P_center, a, P_b(0), B_center, P_center, b) call recentered_cpoly2(P_a(0), Ap_center, P_center, a, P_b(0), Bp_center, P_center, b)
n_new = 0 n_new = 0
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
@ -60,11 +65,17 @@ end
! --- ! ---
subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, alpha, beta, a, b, A_center, B_center, dim) subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder, &
alpha, beta, a, b, Ae_center, Be_center, Ap_center, Bp_center, dim)
BEGIN_DOC BEGIN_DOC
!
! Transforms the product of ! Transforms the product of
! (x-x_A)^a(1) (x-x_B)^b(1) (y-y_A)^a(2) (y-y_B)^b(2) (z-z_A)^a(3) (z-z_B)^b(3) exp(-(r-A)^2 alpha) exp(-(r-B)^2 beta) !
! (x - x_Ap)^a(1) (x - x_Bp)^b(1) exp(-alpha (x - x_Ae)^2) exp(-beta (x - x_Be)^2) x
! (y - y_Ap)^a(2) (y - y_Bp)^b(2) exp(-alpha (y - y_Ae)^2) exp(-beta (y - y_Be)^2) x
! (z - z_Ap)^a(3) (z - z_Bp)^b(3) exp(-alpha (z - z_Ae)^2) exp(-beta (z - z_Be)^2)
!
! into ! into
! fact_k * [sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x] exp (-p (x-P_center(1))^2) ! fact_k * [sum (l_x = 0,i_order(1)) P_new(l_x,1) * (x-P_center(1))^l_x] exp (-p (x-P_center(1))^2)
! * [sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y] exp (-p (y-P_center(2))^2) ! * [sum (l_y = 0,i_order(2)) P_new(l_y,2) * (y-P_center(2))^l_y] exp (-p (y-P_center(2))^2)
@ -73,18 +84,19 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
! WARNING ::: IF fact_k is too smal then: ! WARNING ::: IF fact_k is too smal then:
! returns a "s" function centered in zero ! returns a "s" function centered in zero
! with an inifinite exponent and a zero polynom coef ! with an inifinite exponent and a zero polynom coef
!
END_DOC END_DOC
implicit none implicit none
include 'constants.include.F' include 'constants.include.F'
integer, intent(in) :: dim, a(3), b(3) integer, intent(in) :: dim, a(3), b(3)
complex*16, intent(in) :: alpha, beta, A_center(3), B_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, intent(out) :: iorder(3) integer, intent(out) :: iorder(3)
complex*16, intent(out) :: p, P_center(3), fact_k, P_new(0:max_dim,3) complex*16, intent(out) :: p, P_center(3), fact_k, P_new(0:max_dim,3)
integer :: n_new, i, j integer :: n_new, i, j
double precision :: tmp_mod
complex*16 :: P_a(0:max_dim,3), P_b(0:max_dim,3) complex*16 :: P_a(0:max_dim,3), P_b(0:max_dim,3)
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: P_a, P_b
@ -97,12 +109,11 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
P_new(0,3) = (0.d0, 0.d0) P_new(0,3) = (0.d0, 0.d0)
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call cgaussian_product(alpha, A_center, beta, B_center, fact_k, p, P_center) call cgaussian_product(alpha, Ae_center, beta, Be_center, fact_k, p, P_center)
! IF fact_k is too smal then: returns a "s" function centered in zero ! IF fact_k is too smal then: returns a "s" function centered in zero
! with an inifinite exponent and a zero polynom coef ! with an inifinite exponent and a zero polynom coef
tmp_mod = dsqrt(real(fact_k)*real(fact_k) + aimag(fact_k)*aimag(fact_k)) if(zabs(fact_k) < 1d-14) then
if(tmp_mod < 1d-14) then
iorder = 0 iorder = 0
p = (1.d+14, 0.d0) p = (1.d+14, 0.d0)
fact_k = (0.d0 , 0.d0) fact_k = (0.d0 , 0.d0)
@ -112,7 +123,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
endif endif
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0,1), A_center(1), P_center(1), a(1), P_b(0,1), B_center(1), P_center(1), b(1)) call recentered_cpoly2(P_a(0,1), Ap_center(1), P_center(1), a(1), P_b(0,1), Bp_center(1), P_center(1), b(1))
iorder(1) = a(1) + b(1) iorder(1) = a(1) + b(1)
do i = 0, iorder(1) do i = 0, iorder(1)
P_new(i,1) = 0.d0 P_new(i,1) = 0.d0
@ -122,7 +133,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
call multiply_cpoly(P_a(0,1), a(1), P_b(0,1), b(1), P_new(0,1), n_new) call multiply_cpoly(P_a(0,1), a(1), P_b(0,1), b(1), P_new(0,1), n_new)
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0,2), A_center(2), P_center(2), a(2), P_b(0,2), B_center(2), P_center(2), b(2)) call recentered_cpoly2(P_a(0,2), Ap_center(2), P_center(2), a(2), P_b(0,2), Bp_center(2), P_center(2), b(2))
iorder(2) = a(2) + b(2) iorder(2) = a(2) + b(2)
do i = 0, iorder(2) do i = 0, iorder(2)
P_new(i,2) = 0.d0 P_new(i,2) = 0.d0
@ -132,7 +143,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
call multiply_cpoly(P_a(0,2), a(2), P_b(0,2), b(2), P_new(0,2), n_new) call multiply_cpoly(P_a(0,2), a(2), P_b(0,2), b(2), P_new(0,2), n_new)
!DIR$ FORCEINLINE !DIR$ FORCEINLINE
call recentered_cpoly2(P_a(0,3), A_center(3), P_center(3), a(3), P_b(0,3), B_center(3), P_center(3), b(3)) call recentered_cpoly2(P_a(0,3), Ap_center(3), P_center(3), a(3), P_b(0,3), Bp_center(3), P_center(3), b(3))
iorder(3) = a(3) + b(3) iorder(3) = a(3) + b(3)
do i = 0, iorder(3) do i = 0, iorder(3)
P_new(i,3) = 0.d0 P_new(i,3) = 0.d0
@ -156,13 +167,12 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
complex*16, intent(in) :: a, b, xa(3), xb(3) complex*16, intent(in) :: a, b, xa(3), xb(3)
complex*16, intent(out) :: p, k, xp(3) complex*16, intent(out) :: p, k, xp(3)
double precision :: tmp_mod
complex*16 :: p_inv, xab(3), ab complex*16 :: p_inv, xab(3), ab
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: xab
ASSERT (REAL(a) > 0.) ASSERT (real(a) > 0.)
ASSERT (REAL(b) > 0.) ASSERT (real(b) > 0.)
! new exponent ! new exponent
p = a + b p = a + b
@ -175,8 +185,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
ab = a * b * p_inv ab = a * b * p_inv
k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3)) k = ab * (xab(1)*xab(1) + xab(2)*xab(2) + xab(3)*xab(3))
tmp_mod = dsqrt(real(k)*real(k) + aimag(k)*aimag(k)) if(real(k) .gt. 40.d0) then
if(tmp_mod .gt. 40.d0) then
k = (0.d0, 0.d0) k = (0.d0, 0.d0)
xp(1:3) = (0.d0, 0.d0) xp(1:3) = (0.d0, 0.d0)
return return
@ -199,15 +208,15 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp)
END_DOC END_DOC
implicit none implicit none
complex*16, intent(in) :: a, b, xa, xb complex*16, intent(in) :: a, b, xa, xb
complex*16, intent(out) :: p, k, xp complex*16, intent(out) :: p, k, xp
double precision :: tmp_mod
complex*16 :: p_inv complex*16 :: p_inv
complex*16 :: xab, ab complex*16 :: xab, ab
ASSERT (REAL(a) > 0.) ASSERT (real(a) > 0.)
ASSERT (REAL(b) > 0.) ASSERT (real(b) > 0.)
! new center ! new center
p = a + b p = a + b
@ -218,8 +227,7 @@ subroutine cgaussian_product_x(a, xa, b, xb, k, p, xp)
ab = a * b * p_inv ab = a * b * p_inv
k = ab * xab * xab k = ab * xab * xab
tmp_mod = dsqrt(REAL(k)*REAL(k) + AIMAG(k)*AIMAG(k)) if(real(k) > 40.d0) then
if(tmp_mod > 40.d0) then
k = (0.d0, 0.d0) k = (0.d0, 0.d0)
xp = (0.d0, 0.d0) xp = (0.d0, 0.d0)
return return
@ -290,7 +298,6 @@ subroutine add_cpoly(b, nb, c, nc, d, nd)
complex*16, intent(out) :: d(0:nb+nc) complex*16, intent(out) :: d(0:nb+nc)
integer :: ib integer :: ib
double precision :: tmp_mod
complex*16 :: tmp complex*16 :: tmp
nd = nb + nc nd = nb + nc
@ -299,11 +306,9 @@ subroutine add_cpoly(b, nb, c, nc, d, nd)
enddo enddo
tmp = d(nd) tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) do while( (zabs(tmp) .lt. 1.d-15) .and. (nd >= 0) )
do while( (tmp_mod .lt. 1.d-15) .and. (nd >= 0) )
nd -= 1 nd -= 1
tmp = d(nd) tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
if(nd < 0) exit if(nd < 0) exit
enddo enddo
@ -326,7 +331,6 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd)
complex*16, intent(inout) :: d(0:max(nb, nd)) complex*16, intent(inout) :: d(0:max(nb, nd))
integer :: ib integer :: ib
double precision :: tmp_mod
complex*16 :: tmp complex*16 :: tmp
nd = max(nd, nb) nd = max(nd, nb)
@ -337,12 +341,10 @@ subroutine add_cpoly_multiply(b, nb, cst, d, nd)
enddo enddo
tmp = d(nd) tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp)) do while(zabs(tmp) .lt. 1.d-15)
do while(tmp_mod .lt. 1.d-15)
nd -= 1 nd -= 1
if(nd < 0) exit if(nd < 0) exit
tmp = d(nd) tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
enddo enddo
endif endif

308
src/utils/cpx_boys.irp.f
View File

@ -10,161 +10,219 @@ complex*16 function crint_1(n, rho)
complex*16, intent(in) :: rho complex*16, intent(in) :: rho
integer :: i, mmax integer :: i, mmax
double precision :: rho_mod, rho_re, rho_im double precision :: rho_mod
double precision :: sq_rho_re, sq_rho_im double precision :: tmp
double precision :: n_tmp complex*16 :: rho_inv, rho_exp
complex*16 :: sq_rho, rho_inv, rho_exp
complex*16 :: crint_smallz, cpx_erf_1 complex*16 :: crint_smallz
rho_re = real (rho) rho_mod = zabs(rho)
rho_im = aimag(rho)
rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im) if(rho_mod < 3.5d0) then
if(rho_mod .lt. 0.35d0) then
select case(n)
case(0)
crint_1 = (((((((((1.3122532963802805073d-08 * rho &
- 1.450385222315046877d-07) * rho &
+ 1.458916900093370682d-06) * rho &
- 0.132275132275132275d-04) * rho &
+ 0.106837606837606838d-03) * rho &
- 0.757575757575757576d-03) * rho &
+ 0.462962962962962963d-02) * rho &
- 0.238095238095238095d-01) * rho &
+ 0.10000000000000000000d0) * rho &
- 0.33333333333333333333d0) * rho &
+ 1.0d0
case(1)
crint_1 = (((((((((1.198144314086343d-08 * rho &
- 1.312253296380281d-07) * rho &
+ 1.305346700083542d-06) * rho &
- 1.167133520074696d-05) * rho &
+ 9.259259259259259d-05) * rho &
- 6.410256410256410d-04) * rho &
+ 3.787878787878788d-03) * rho &
- 1.851851851851852d-02) * rho &
+ 7.142857142857142d-02) * rho &
- 2.000000000000000d-01) * rho &
+ 3.333333333333333d-01
case(2)
crint_1 = (((((((((1.102292768959436d-08 * rho &
- 1.198144314086343d-07) * rho &
+ 1.181027966742252d-06) * rho &
- 1.044277360066834d-05) * rho &
+ 8.169934640522875d-05) * rho &
- 5.555555555555556d-04) * rho &
+ 3.205128205128205d-03) * rho &
- 1.515151515151515d-02) * rho &
+ 5.555555555555555d-02) * rho &
- 1.428571428571428d-01) * rho &
+ 2.000000000000000d-01
case(3)
crint_1 = (((((((((1.020641452740218d-08 * rho &
- 1.102292768959436d-07) * rho &
+ 1.078329882677709d-06) * rho &
- 9.448223733938020d-06) * rho &
+ 7.309941520467836d-05) * rho &
- 4.901960784313725d-04) * rho &
+ 2.777777777777778d-03) * rho &
- 1.282051282051282d-02) * rho &
+ 4.545454545454546d-02) * rho &
- 1.111111111111111d-01) * rho &
+ 1.428571428571428d-01
case default
tmp = dble(n + n + 1)
crint_1 = (((((((((2.755731922398589d-07 * rho / (tmp + 20.d0) &
- 2.755731922398589d-06 / (tmp + 18.d0)) * rho &
+ 2.480158730158730d-05 / (tmp + 16.d0)) * rho &
- 1.984126984126984d-04 / (tmp + 14.d0)) * rho &
+ 1.388888888888889d-03 / (tmp + 12.d0)) * rho &
- 8.333333333333333d-03 / (tmp + 10.d0)) * rho &
+ 4.166666666666666d-02 / (tmp + 8.d0)) * rho &
- 1.666666666666667d-01 / (tmp + 6.d0)) * rho &
+ 5.000000000000000d-01 / (tmp + 4.d0)) * rho &
- 1.000000000000000d+00 / (tmp + 2.d0)) * rho &
+ 1.0d0 / tmp
end select
if(rho_mod < 10.d0) then
! small z
if(rho_mod .lt. 1.d-15) then
crint_1 = 1.d0 / dble(n + n + 1)
else else
crint_1 = crint_smallz(n, rho) crint_1 = crint_smallz(n, rho)
endif endif
else else
! large z
if(rho_mod .gt. 40.d0) then
n_tmp = dble(n) + 0.5d0
crint_1 = 0.5d0 * gamma(n_tmp) / (rho**n_tmp)
else
! get \sqrt(rho)
sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
sq_rho_im = 0.5d0 * rho_im / sq_rho_re
sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
rho_exp = 0.5d0 * zexp(-rho) rho_exp = 0.5d0 * zexp(-rho)
rho_inv = (1.d0, 0.d0) / rho rho_inv = (1.d0, 0.d0) / rho
crint_1 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho call zboysfun00_1(rho, crint_1)
mmax = n do i = 1, n
if(mmax .gt. 0) then crint_1 = ((dble(i) - 0.5d0) * crint_1 - rho_exp) * rho_inv
do i = 0, mmax-1
crint_1 = ((dble(i) + 0.5d0) * crint_1 - rho_exp) * rho_inv
enddo enddo
endif
endif
endif endif
return
end end
! --- ! ---
complex*16 function crint_sum_1(n_pt_out, rho, d1) subroutine crint_1_vec(n_max, rho, vals)
implicit none implicit none
include 'constants.include.F' include 'constants.include.F'
integer, intent(in) :: n_pt_out integer, intent(in) :: n_max
complex*16, intent(in) :: rho, d1(0:n_pt_out) complex*16, intent(in) :: rho
complex*16, intent(out) :: vals(0:n_max)
integer :: n, i, mmax integer :: n
double precision :: rho_mod, rho_re, rho_im double precision :: rho_mod
double precision :: sq_rho_re, sq_rho_im double precision :: tmp
complex*16 :: sq_rho, F0 complex*16 :: rho_inv, rho_exp
complex*16 :: rho_tmp, rho_inv, rho_exp
complex*16, allocatable :: Fm(:)
complex*16 :: crint_smallz, cpx_erf_1 complex*16 :: crint_smallz
rho_mod = zabs(rho)
rho_re = real (rho) if(rho_mod < 3.5d0) then
rho_im = aimag(rho)
rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
! ! debug if(rho_mod .lt. 0.35d0) then
! double precision :: d1_real(0:n_pt_out)
! double precision :: rint_sum
! do i = 0, n_pt_out
! d1_real(i) = real(d1(i))
! enddo
! crint_sum_1 = rint_sum(n_pt_out, rho_re, d1_real)
! return
if(rho_mod < 10.d0) then vals(0) = (((((((((1.3122532963802805073d-08 * rho &
! small z - 1.450385222315046877d-07) * rho &
+ 1.458916900093370682d-06) * rho &
- 0.132275132275132275d-04) * rho &
+ 0.106837606837606838d-03) * rho &
- 0.757575757575757576d-03) * rho &
+ 0.462962962962962963d-02) * rho &
- 0.238095238095238095d-01) * rho &
+ 0.10000000000000000000d0) * rho &
- 0.33333333333333333333d0) * rho &
+ 1.0d0
if(rho_mod .lt. 1.d-15) then if(n_max > 0) then
crint_sum_1 = d1(0) vals(1) = (((((((((1.198144314086343d-08 * rho &
do i = 2, n_pt_out, 2 - 1.312253296380281d-07) * rho &
n = shiftr(i, 1) + 1.305346700083542d-06) * rho &
crint_sum_1 = crint_sum_1 + d1(i) / dble(n+n+1) - 1.167133520074696d-05) * rho &
+ 9.259259259259259d-05) * rho &
- 6.410256410256410d-04) * rho &
+ 3.787878787878788d-03) * rho &
- 1.851851851851852d-02) * rho &
+ 7.142857142857142d-02) * rho &
- 2.000000000000000d-01) * rho &
+ 3.333333333333333d-01
if(n_max > 1) then
vals(2) = (((((((((1.102292768959436d-08 * rho &
- 1.198144314086343d-07) * rho &
+ 1.181027966742252d-06) * rho &
- 1.044277360066834d-05) * rho &
+ 8.169934640522875d-05) * rho &
- 5.555555555555556d-04) * rho &
+ 3.205128205128205d-03) * rho &
- 1.515151515151515d-02) * rho &
+ 5.555555555555555d-02) * rho &
- 1.428571428571428d-01) * rho &
+ 2.000000000000000d-01
if(n_max > 2) then
vals(3) = (((((((((1.020641452740218d-08 * rho &
- 1.102292768959436d-07) * rho &
+ 1.078329882677709d-06) * rho &
- 9.448223733938020d-06) * rho &
+ 7.309941520467836d-05) * rho &
- 4.901960784313725d-04) * rho &
+ 2.777777777777778d-03) * rho &
- 1.282051282051282d-02) * rho &
+ 4.545454545454546d-02) * rho &
- 1.111111111111111d-01) * rho &
+ 1.428571428571428d-01
do n = 4, n_max
tmp = dble(n + n + 1)
vals(n) = (((((((((2.755731922398589d-07 * rho / (tmp + 20.d0) &
- 2.755731922398589d-06 / (tmp + 18.d0)) * rho &
+ 2.480158730158730d-05 / (tmp + 16.d0)) * rho &
- 1.984126984126984d-04 / (tmp + 14.d0)) * rho &
+ 1.388888888888889d-03 / (tmp + 12.d0)) * rho &
- 8.333333333333333d-03 / (tmp + 10.d0)) * rho &
+ 4.166666666666666d-02 / (tmp + 8.d0)) * rho &
- 1.666666666666667d-01 / (tmp + 6.d0)) * rho &
+ 5.000000000000000d-01 / (tmp + 4.d0)) * rho &
- 1.000000000000000d+00 / (tmp + 2.d0)) * rho &
+ 1.0d0 / tmp
enddo enddo
endif ! n_max > 2
endif ! n_max > 1
endif ! n_max > 0
else else
crint_sum_1 = d1(0) * crint_smallz(0, rho) call crint_smallz_vec(n_max, rho, vals)
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
crint_sum_1 = crint_sum_1 + d1(i) * crint_smallz(n, rho)
enddo
endif endif
else else
! large z
if(rho_mod .gt. 40.d0) then
rho_inv = (1.d0, 0.d0) / rho
rho_tmp = 0.5d0 * sqpi * zsqrt(rho_inv)
crint_sum_1 = rho_tmp * d1(0)
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
rho_tmp = rho_tmp * (dble(n) + 0.5d0) * rho_inv
crint_sum_1 = crint_sum_1 + rho_tmp * d1(i)
enddo
else
! get \sqrt(rho)
sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
sq_rho_im = 0.5d0 * rho_im / sq_rho_re
sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
F0 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho
crint_sum_1 = F0 * d1(0)
rho_exp = 0.5d0 * zexp(-rho) rho_exp = 0.5d0 * zexp(-rho)
rho_inv = (1.d0, 0.d0) / rho rho_inv = (1.d0, 0.d0) / rho
mmax = shiftr(n_pt_out, 1) call zboysfun00_1(rho, vals(0))
if(mmax .gt. 0) then do n = 1, n_max
vals(n) = ((dble(n) - 0.5d0) * vals(n-1) - rho_exp) * rho_inv
allocate(Fm(mmax))
Fm(1:mmax) = (0.d0, 0.d0)
do n = 0, mmax-1
F0 = ((dble(n) + 0.5d0) * F0 - rho_exp) * rho_inv
Fm(n+1) = F0
enddo enddo
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
crint_sum_1 = crint_sum_1 + Fm(n) * d1(i)
enddo
deallocate(Fm)
endif endif
endif ! rho_mod return
endif ! rho_mod
end end
! --- ! ---
@ -219,15 +277,16 @@ complex*16 function crint_2(n, rho)
integer, intent(in) :: n integer, intent(in) :: n
complex*16, intent(in) :: rho complex*16, intent(in) :: rho
double precision :: tmp double precision :: tmp, abs_rho
complex*16 :: rho2
complex*16 :: vals(0:n) complex*16 :: vals(0:n)
complex*16, external :: crint_smallz complex*16, external :: crint_smallz
if(abs(rho) < 3.5d0) then abs_rho = zabs(rho)
if(abs(rho) .lt. 0.35d0) then if(abs_rho < 3.5d0) then
if(abs_rho .lt. 0.35d0) then
select case(n) select case(n)
case(0) case(0)
@ -343,7 +402,7 @@ subroutine zboysfun(n_max, x, vals)
integer :: n integer :: n
complex*16 :: yy, x_inv complex*16 :: yy, x_inv
call zboysfun00(x, vals(0)) call zboysfun00_2(x, vals(0))
yy = 0.5d0 * zexp(-x) yy = 0.5d0 * zexp(-x)
x_inv = (1.d0, 0.d0) / x x_inv = (1.d0, 0.d0) / x
@ -393,7 +452,7 @@ end
! --- ! ---
complex*16 function crint_sum_2(n_pt_out, rho, d1) complex*16 function crint_sum(n_pt_out, rho, d1)
implicit none implicit none
@ -408,13 +467,14 @@ complex*16 function crint_sum_2(n_pt_out, rho, d1)
n_max = shiftr(n_pt_out, 1) n_max = shiftr(n_pt_out, 1)
allocate(vals(0:n_max)) allocate(vals(0:n_max))
call crint_1_vec(n_max, rho, vals)
!call crint_2_vec(n_max, rho, vals) !call crint_2_vec(n_max, rho, vals)
! FOR DEBUG ! FOR DEBUG
call crint_quad_12_vec(n_max, rho, vals) !call crint_quad_12_vec(n_max, rho, vals)
crint_sum_2 = d1(0) * vals(0) crint_sum = d1(0) * vals(0)
do i = 2, n_pt_out, 2 do i = 2, n_pt_out, 2
crint_sum_2 += d1(i) * vals(shiftr(i, 1)) crint_sum += d1(i) * vals(shiftr(i, 1))
enddo enddo
deallocate(vals) deallocate(vals)
@ -434,7 +494,7 @@ subroutine crint_2_vec(n_max, rho, vals)
integer :: n integer :: n
double precision :: tmp, abs_rho double precision :: tmp, abs_rho
complex*16 :: rho2, rho3, erho complex*16 :: erho
abs_rho = abs(rho) abs_rho = abs(rho)
@ -455,7 +515,7 @@ subroutine crint_2_vec(n_max, rho, vals)
- 0.33333333333333333333d0) * rho & - 0.33333333333333333333d0) * rho &
+ 1.0d0 + 1.0d0
if(n > 0) then if(n_max > 0) then
vals(1) = (((((((((1.198144314086343d-08 * rho & vals(1) = (((((((((1.198144314086343d-08 * rho &
- 1.312253296380281d-07) * rho & - 1.312253296380281d-07) * rho &
@ -469,7 +529,7 @@ subroutine crint_2_vec(n_max, rho, vals)
- 2.000000000000000d-01) * rho & - 2.000000000000000d-01) * rho &
+ 3.333333333333333d-01 + 3.333333333333333d-01
if(n > 1) then if(n_max > 1) then
vals(2) = (((((((((1.102292768959436d-08 * rho & vals(2) = (((((((((1.102292768959436d-08 * rho &
- 1.198144314086343d-07) * rho & - 1.198144314086343d-07) * rho &
@ -483,7 +543,7 @@ subroutine crint_2_vec(n_max, rho, vals)
- 1.428571428571428d-01) * rho & - 1.428571428571428d-01) * rho &
+ 2.000000000000000d-01 + 2.000000000000000d-01
if(n > 2) then if(n_max > 2) then
vals(3) = (((((((((1.020641452740218d-08 * rho & vals(3) = (((((((((1.020641452740218d-08 * rho &
- 1.102292768959436d-07) * rho & - 1.102292768959436d-07) * rho &
@ -512,9 +572,9 @@ subroutine crint_2_vec(n_max, rho, vals)
+ 1.0d0 / tmp + 1.0d0 / tmp
enddo enddo
endif ! n > 2 endif ! n_max > 2
endif ! n > 1 endif ! n_max > 1
endif ! n > 0 endif ! n_max > 0
else else
@ -855,8 +915,6 @@ subroutine crint_quad_12_vec(n_max, rho, vals)
complex*16, intent(out) :: vals(0:n_max) complex*16, intent(out) :: vals(0:n_max)
integer :: n integer :: n
double precision :: tmp, abs_rho
complex*16 :: rho2, rho3, erho
do n = 0, n_max do n = 0, n_max
call crint_quad_12(n, rho, 10000000, vals(n)) call crint_quad_12(n, rho, 10000000, vals(n))

95
src/utils/cpx_erf.irp.f
View File

@ -1,6 +1,37 @@
! --- ! ---
subroutine zboysfun00_1(rho, F0)
implicit none
include 'constants.include.F'
complex*16, intent(in) :: rho
complex*16, intent(out) :: F0
double precision :: rho_re, rho_im, rho_mod
double precision :: sq_rho_re, sq_rho_im
complex*16 :: sq_rho
complex*16, external :: cpx_erf_1
rho_re = real(rho)
rho_im = aimag(rho)
rho_mod = dsqrt(rho_re*rho_re + rho_im*rho_im)
sq_rho_re = sq_op5 * dsqrt(rho_re + rho_mod)
sq_rho_im = 0.5d0 * rho_im / sq_rho_re
sq_rho = sq_rho_re + (0.d0, 1.d0) * sq_rho_im
F0 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho
return
end
! ---
complex*16 function cpx_erf_1(x, y) complex*16 function cpx_erf_1(x, y)
BEGIN_DOC BEGIN_DOC
@ -43,6 +74,7 @@ complex*16 function cpx_erf_1(x, y)
cpx_erf_1 = conjg(erf_tot) cpx_erf_1 = conjg(erf_tot)
endif endif
return
end end
! --- ! ---
@ -54,22 +86,42 @@ complex*16 function erf_E(x, yabs)
double precision, intent(in) :: x, yabs double precision, intent(in) :: x, yabs
if((dabs(x).gt.6.d0) .or. (x==0.d0)) then double precision :: xy
if(dabs(x) .gt. 6.d0) then
erf_E = (0.d0, 0.d0) erf_E = (0.d0, 0.d0)
return return
endif endif
if(dabs(x) .lt. 1.d-7) then xy = x * yabs
erf_E = -inv_pi * (0.d0, 1.d0) * yabs if(dabs(xy) .lt. 0.1d0) then
erf_E = ((((((((((((-(0.d0, 9.3968347936601903d-8) * xy + (6.5777843555621328d-7, 0.d0) &
) * xy + (0.d0, 4.2755598311153869d-6) &
) * xy - (2.5653358986692322d-5, 0.d0) &
) * xy - (0.d0, 0.00014109347442680775d0) &
) * xy + (0.00070546737213403885d0, 0.d0) &
) * xy + (0.d0, 0.0031746031746031746d0) &
) * xy - (0.012698412698412698d0, 0.d0) &
) * xy - (0.d0, 0.044444444444444446d0) &
) * xy + (0.13333333333333333d0, 0.d0) &
) * xy + (0.d0, 0.33333333333333331d0) &
) * xy - (0.66666666666666663d0, 0.d0) &
) * xy - (0.d0, 1.d0) &
) * xy + (1.d0, 0.d0)
erf_E = erf_E * (0.d0, 1.d0) * yabs * inv_pi * dexp(-x*x)
else else
erf_E = 0.5d0 * inv_pi * dexp(-x*x) & erf_E = 0.5d0 * inv_pi * dexp(-x*x) &
* ((1.d0, 0.d0) - zexp(-(2.d0, 0.d0) * x * yabs)) / x * ((1.d0, 0.d0) - zexp(-(0.d0, 2.d0) * xy)) / x
endif endif
return
end end
! --- ! ---
@ -84,31 +136,30 @@ double precision function erf_F(x, yabs)
integer, parameter :: Nmax = 13 integer, parameter :: Nmax = 13
integer :: i integer :: i
double precision :: tmp1, tmp2, x2, ct double precision :: tmp1, tmp2, x2
if(dabs(x) .gt. 5.8d0) then if(dabs(x) .gt. 5.8d0) then
erf_F = 0.d0 erf_F = 0.d0
return
endif
else
x2 = x * x x2 = x * x
ct = x * inv_pi
erf_F = 0.d0 erf_F = 0.d0
do i = 1, Nmax do i = 1, Nmax
tmp1 = 0.25d0 * dble(i) * dble(i) + x2 tmp1 = 0.25d0 * dble(i*i) + x2
tmp2 = dexp(-tmp1) / tmp1 tmp2 = dexp(-tmp1) / tmp1
erf_F = erf_F + tmp2 erf_F = erf_F + tmp2
if(dabs(tmp2) .lt. 1d-15) exit if(tmp2 .lt. 1d-15) exit
enddo enddo
erf_F = ct * erf_F
endif erf_F = x * inv_pi * erf_F
return
end end
! --- ! ---
@ -149,6 +200,7 @@ complex*16 function erf_G(x, yabs)
enddo enddo
return
end end
! --- ! ---
@ -163,19 +215,14 @@ complex*16 function erf_H(x, yabs)
integer, parameter :: Nmax = 13 integer, parameter :: Nmax = 13
integer :: i integer :: i
double precision :: tmp0, tmp1, tmp_mod, x2, ct, idble double precision :: tmp0, tmp1, x2, idble
complex*16 :: tmp2 complex*16 :: tmp2
if(x .eq. 0.d0) then
erf_H = (0.d0, 0.d0)
return
endif
if((dabs(x) .lt. 107d0) .and. (yabs .lt. 6.1d0)) then
if((dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0)) then !if((dabs(x) .lt. 10d0) .and. (yabs .lt. 6.1d0)) then
x2 = x * x x2 = x * x
ct = 0.5d0 * inv_pi
erf_H = 0.d0 erf_H = 0.d0
do i = 1, Nmax do i = 1, Nmax
@ -186,10 +233,10 @@ complex*16 function erf_H(x, yabs)
tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1 tmp2 = dexp(-tmp1-idble*yabs) * (x + (0.d0, 1.d0)*tmp0) / tmp1
erf_H = erf_H + tmp2 erf_H = erf_H + tmp2
tmp_mod = dsqrt(real(tmp2)*real(tmp2) + aimag(tmp2)*aimag(tmp2)) if(zabs(tmp2) .lt. 1d-15) exit
if(tmp_mod .lt. 1d-15) exit
enddo enddo
erf_H = ct * erf_H
erf_H = 0.5d0 * inv_pi * erf_H
else else
@ -201,7 +248,7 @@ end
! --- ! ---
subroutine zboysfun00(z, val) subroutine zboysfun00_2(z, val)
BEGIN_DOC BEGIN_DOC
! !