PTEROSOR/QuantumPackage
10
0
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-12-22 04:13:55 +01:00
Code Releases Activity

Merge branch 'dev-stable' into dev-stable

Browse Source
This commit is contained in:
AbdAmmar 2024-10-15 00:30:08 +02:00 committed by GitHub
commit 600326e975
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
24 changed files with 3378 additions and 2513 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

28
src/ao_basis/EZFIO.cfg
View File

@ -67,15 +67,27 @@ doc: Use normalized primitive functions
interface: ezfio, provider
default: true
[ao_expoim_cosgtos]
type: double precision
doc: imag part for Exponents for each primitive of each cosGTOs |AO|
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider
[use_cosgtos]
[use_cgtos]
type: logical
doc: If true, use cosgtos for AO integrals
doc: If true, use cgtos for AO integrals
interface: ezfio
default: False
[ao_expo_im_cgtos]
type: double precision
doc: imag part for Exponents for each primitive of each cGTOs |AO|
size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider
[ao_expo_pw]
type: double precision
doc: plane wave part for each primitive GTOs |AO|
size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider
[ao_expo_phase]
type: double precision
doc: phase shift for each primitive GTOs |AO|
size: (3,ao_basis.ao_num,ao_basis.ao_prim_num_max)
interface: ezfio, provider

37
src/ao_basis/cgtos.irp.f Normal file
View File

@ -0,0 +1,37 @@
BEGIN_PROVIDER [logical, use_cgtos]
implicit none
BEGIN_DOC
! If true, use cgtos for AO integrals
END_DOC
logical :: has
PROVIDE ezfio_filename
use_cgtos = .False.
if (mpi_master) then
call ezfio_has_ao_basis_use_cgtos(has)
if (has) then
! write(6,'(A)') '.. >>>>> [ IO READ: use_cgtos ] <<<<< ..'
call ezfio_get_ao_basis_use_cgtos(use_cgtos)
else
call ezfio_set_ao_basis_use_cgtos(use_cgtos)
endif
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( use_cgtos, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read use_cgtos with MPI'
endif
IRP_ENDIF
! call write_time(6)
END_PROVIDER

34
src/ao_basis/cosgtos.irp.f
View File

@ -1,34 +0,0 @@
BEGIN_PROVIDER [ logical, use_cosgtos ]
implicit none
BEGIN_DOC
! If true, use cosgtos for AO integrals
END_DOC
logical :: has
PROVIDE ezfio_filename
use_cosgtos = .False.
if (mpi_master) then
call ezfio_has_ao_basis_use_cosgtos(has)
if (has) then
! write(6,'(A)') '.. >>>>> [ IO READ: use_cosgtos ] <<<<< ..'
call ezfio_get_ao_basis_use_cosgtos(use_cosgtos)
else
call ezfio_set_ao_basis_use_cosgtos(use_cosgtos)
endif
endif
IRP_IF MPI_DEBUG
print *, irp_here, mpi_rank
call MPI_BARRIER(MPI_COMM_WORLD, ierr)
IRP_ENDIF
IRP_IF MPI
include 'mpif.h'
integer :: ierr
call MPI_BCAST( use_cosgtos, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
if (ierr /= MPI_SUCCESS) then
stop 'Unable to read use_cosgtos with MPI'
endif
IRP_ENDIF
! call write_time(6)
END_PROVIDER

21
src/ao_one_e_ints/ao_overlap.irp.f
View File

@ -1,10 +1,10 @@
! ---
BEGIN_PROVIDER [ double precision, ao_overlap , (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_x, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_y, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_z, (ao_num, ao_num) ]
BEGIN_PROVIDER [double precision, ao_overlap , (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_x, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_y, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_z, (ao_num, ao_num)]
BEGIN_DOC
! Overlap between atomic basis functions:
@ -30,15 +30,14 @@
else
if(use_cosgtos) then
!print*, ' use_cosgtos for ao_overlap ?', use_cosgtos
if(use_cgtos) then
do j = 1, ao_num
do i = 1, ao_num
ao_overlap (i,j) = ao_overlap_cosgtos (i,j)
ao_overlap_x(i,j) = ao_overlap_cosgtos_x(i,j)
ao_overlap_y(i,j) = ao_overlap_cosgtos_y(i,j)
ao_overlap_z(i,j) = ao_overlap_cosgtos_z(i,j)
ao_overlap (i,j) = ao_overlap_cgtos (i,j)
ao_overlap_x(i,j) = ao_overlap_cgtos_x(i,j)
ao_overlap_y(i,j) = ao_overlap_cgtos_y(i,j)
ao_overlap_z(i,j) = ao_overlap_cgtos_z(i,j)
enddo
enddo
@ -49,7 +48,7 @@
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(A_center,B_center,power_A,power_B,&
!$OMP overlap_x,overlap_y, overlap_z, overlap, &
!$OMP alpha, beta,i,j,c) &
!$OMP alpha, beta,i,j,n,l,c) &
!$OMP SHARED(nucl_coord,ao_power,ao_prim_num, &
!$OMP ao_overlap_x,ao_overlap_y,ao_overlap_z,ao_overlap,ao_num,ao_coef_normalized_ordered_transp,ao_nucl, &
!$OMP ao_expo_ordered_transp,dim1)

267
src/ao_one_e_ints/aos_cgtos.irp.f Normal file
View File

@ -0,0 +1,267 @@
! ---
BEGIN_PROVIDER [double precision, ao_coef_cgtos_norm_ord_transp, (ao_prim_num_max, ao_num)]
implicit none
integer :: i, j
do j = 1, ao_num
do i = 1, ao_prim_num_max
ao_coef_cgtos_norm_ord_transp(i,j) = ao_coef_norm_cgtos_ord(j,i)
enddo
enddo
END_PROVIDER
! ---
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 [complex*16, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)]
implicit none
integer :: i, j, m
do j = 1, ao_num
do i = 1, ao_prim_num_max
ao_expo_cgtos_ord_transp(i,j) = ao_expo_cgtos_ord(j,i)
do m = 1, 3
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)
enddo
ao_expo_pw_ord_transp(4,i,j) = ao_expo_pw_ord_transp(1,i,j) &
+ ao_expo_pw_ord_transp(2,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
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
implicit none
integer :: i, j, powA(3), nz
double precision :: norm
complex*16 :: overlap_x, overlap_y, overlap_z, C_A(3)
complex*16 :: integ1, integ2, expo
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
do i = 1, ao_num
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
! TODO
! Normalization of the primitives
if(primitives_normalized) then
do j = 1, ao_prim_num(i)
expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im_cgtos(i,j)
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_A, C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz)
norm = 2.d0 * real(integ1 + integ2)
ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
enddo
else
do j = 1, ao_prim_num(i)
ao_coef_norm_cgtos(i,j) = ao_coef(i,j)
enddo
endif
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos_ord, (ao_num, ao_prim_num_max)]
&BEGIN_PROVIDER [complex*16 , ao_expo_cgtos_ord, (ao_num, ao_prim_num_max)]
&BEGIN_PROVIDER [double precision, ao_expo_pw_ord, (3, ao_num, ao_prim_num_max)]
&BEGIN_PROVIDER [double precision, ao_expo_phase_ord, (3, ao_num, ao_prim_num_max)]
implicit none
integer :: i, j
integer :: iorder(ao_prim_num_max)
double precision :: d(ao_prim_num_max,9)
d = 0.d0
do i = 1, ao_num
do j = 1, ao_prim_num(i)
iorder(j) = j
d(j,1) = ao_expo(i,j)
d(j,2) = ao_coef_norm_cgtos(i,j)
d(j,3) = ao_expo_im_cgtos(i,j)
d(j,4) = ao_expo_pw(1,i,j)
d(j,5) = ao_expo_pw(2,i,j)
d(j,6) = ao_expo_pw(3,i,j)
d(j,7) = ao_expo_phase(1,i,j)
d(j,8) = ao_expo_phase(2,i,j)
d(j,9) = ao_expo_phase(3,i,j)
enddo
call dsort(d(1,1), iorder, ao_prim_num(i))
do j = 2, 9
call dset_order(d(1,j), iorder, ao_prim_num(i))
enddo
do j = 1, ao_prim_num(i)
ao_expo_cgtos_ord (i,j) = d(j,1) + (0.d0, 1.d0) * d(j,3)
ao_coef_norm_cgtos_ord(i,j) = d(j,2)
ao_expo_pw_ord(i,j,1) = d(j,4)
ao_expo_pw_ord(i,j,2) = d(j,5)
ao_expo_pw_ord(i,j,3) = d(j,6)
ao_expo_phase_ord(i,j,1) = d(j,7)
ao_expo_phase_ord(i,j,2) = d(j,8)
ao_expo_phase_ord(i,j,3) = d(j,9)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_overlap_cgtos, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_x, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_y, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_overlap_cgtos_z, (ao_num, ao_num)]
implicit none
integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3)
double precision :: c, overlap, overlap_x, overlap_y, overlap_z
complex*16 :: alpha, alpha_inv, A_center(3), KA2(3), phiA(3)
complex*16 :: beta, beta_inv, B_center(3), KB2(3), phiB(3)
complex*16 :: C1(1:4), C2(1:4)
complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1
complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2
ao_overlap_cgtos = 0.d0
ao_overlap_cgtos_x = 0.d0
ao_overlap_cgtos_y = 0.d0
ao_overlap_cgtos_z = 0.d0
dim1 = 100
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, &
!$OMP alpha, alpha_inv, A_center, power_A, KA2, phiA, &
!$OMP beta, beta_inv, B_center, power_B, KB2, phiB, &
!$OMP overlap_x , overlap_y , overlap_z , overlap, &
!$OMP overlap_x1, overlap_y1, overlap_z1, overlap1, &
!$OMP overlap_x2, overlap_y2, overlap_z2, overlap2) &
!$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, &
!$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_overlap_cgtos_x, ao_overlap_cgtos_y, ao_overlap_cgtos_z, &
!$OMP ao_overlap_cgtos)
do j = 1, ao_num
jj = ao_nucl(j)
power_A(1) = ao_power(j,1)
power_A(2) = ao_power(j,2)
power_A(3) = ao_power(j,3)
do i = 1, ao_num
ii = ao_nucl(i)
power_B(1) = ao_power(i,1)
power_B(2) = ao_power(i,2)
power_B(3) = ao_power(i,3)
do n = 1, ao_prim_num(j)
alpha = ao_expo_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha
do m = 1, 3
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)
KA2(m) = ao_expo_pw_ord_transp(m,n,j) * ao_expo_pw_ord_transp(m,n,j)
enddo
do l = 1, ao_prim_num(i)
beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta
do m = 1, 3
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)
KB2(m) = ao_expo_pw_ord_transp(m,l,i) * ao_expo_pw_ord_transp(m,l,i)
enddo
c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i)
C1(1) = zexp((0.d0, 1.d0) * (-phiA(1) - phiB(1)) - 0.25d0 * (alpha_inv * KA2(1) + beta_inv * KB2(1)))
C1(2) = zexp((0.d0, 1.d0) * (-phiA(2) - phiB(2)) - 0.25d0 * (alpha_inv * KA2(2) + beta_inv * KB2(2)))
C1(3) = zexp((0.d0, 1.d0) * (-phiA(3) - phiB(3)) - 0.25d0 * (alpha_inv * KA2(3) + beta_inv * KB2(3)))
C1(4) = C1(1) * C1(2) * C1(3)
C2(1) = zexp((0.d0, 1.d0) * (phiA(1) - phiB(1)) - 0.25d0 * (conjg(alpha_inv) * KA2(1) + beta_inv * KB2(1)))
C2(2) = zexp((0.d0, 1.d0) * (phiA(2) - phiB(2)) - 0.25d0 * (conjg(alpha_inv) * KA2(2) + beta_inv * KB2(2)))
C2(3) = zexp((0.d0, 1.d0) * (phiA(3) - phiB(3)) - 0.25d0 * (conjg(alpha_inv) * KA2(3) + beta_inv * KB2(3)))
C2(4) = C2(1) * C2(2) * C2(3)
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_x1, overlap_y1, overlap_z1, overlap1, dim1)
call overlap_cgaussian_xyz(A_center, B_center, conjg(alpha), beta, power_A, power_B, &
overlap_x2, overlap_y2, overlap_z2, overlap2, dim1)
overlap_x = 2.d0 * real(C1(1) * overlap_x1 + C2(1) * overlap_x2)
overlap_y = 2.d0 * real(C1(2) * overlap_y1 + C2(2) * overlap_y2)
overlap_z = 2.d0 * real(C1(3) * overlap_z1 + C2(3) * overlap_z2)
overlap = 2.d0 * real(C1(4) * overlap1 + C2(4) * overlap2 )
ao_overlap_cgtos(i,j) = ao_overlap_cgtos(i,j) + c * overlap
if(isnan(ao_overlap_cgtos(i,j))) then
print*,'i, j', i, j
print*,'l, n', l, n
print*,'c, overlap', c, overlap
print*, overlap_x, overlap_y, overlap_z
stop
endif
ao_overlap_cgtos_x(i,j) = ao_overlap_cgtos_x(i,j) + c * overlap_x
ao_overlap_cgtos_y(i,j) = ao_overlap_cgtos_y(i,j) + c * overlap_y
ao_overlap_cgtos_z(i,j) = ao_overlap_cgtos_z(i,j) + c * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

210
src/ao_one_e_ints/aos_cosgtos.irp.f
View File

@ -1,210 +0,0 @@
! ---
BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
implicit none
integer :: i, j
do j = 1, ao_num
do i = 1, ao_prim_num_max
ao_coef_norm_ord_transp_cosgtos(i,j) = ao_coef_norm_ord_cosgtos(j,i)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ complex*16, ao_expo_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
implicit none
integer :: i, j
do j = 1, ao_num
do i = 1, ao_prim_num_max
ao_expo_ord_transp_cosgtos(i,j) = ao_expo_ord_cosgtos(j,i)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_coef_norm_cosgtos, (ao_num, ao_prim_num_max) ]
implicit none
integer :: i, j, powA(3), nz
double precision :: norm
complex*16 :: overlap_x, overlap_y, overlap_z, C_A(3)
complex*16 :: integ1, integ2, expo
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_cosgtos = 0.d0
do i = 1, ao_num
powA(1) = ao_power(i,1)
powA(2) = ao_power(i,2)
powA(3) = ao_power(i,3)
! Normalization of the primitives
if(primitives_normalized) then
do j = 1, ao_prim_num(i)
expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expoim_cosgtos(i,j)
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_A, C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz)
norm = 2.d0 * real( integ1 + integ2 )
ao_coef_norm_cosgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
enddo
else
do j = 1, ao_prim_num(i)
ao_coef_norm_cosgtos(i,j) = ao_coef(i,j)
enddo
endif
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_cosgtos, (ao_num, ao_prim_num_max) ]
&BEGIN_PROVIDER [ complex*16 , ao_expo_ord_cosgtos, (ao_num, ao_prim_num_max) ]
implicit none
integer :: i, j
integer :: iorder(ao_prim_num_max)
double precision :: d(ao_prim_num_max,3)
d = 0.d0
do i = 1, ao_num
do j = 1, ao_prim_num(i)
iorder(j) = j
d(j,1) = ao_expo(i,j)
d(j,2) = ao_coef_norm_cosgtos(i,j)
d(j,3) = ao_expoim_cosgtos(i,j)
enddo
call dsort (d(1,1), iorder, ao_prim_num(i))
call dset_order(d(1,2), iorder, ao_prim_num(i))
call dset_order(d(1,3), iorder, ao_prim_num(i))
do j = 1, ao_prim_num(i)
ao_expo_ord_cosgtos (i,j) = d(j,1) + (0.d0, 1.d0) * d(j,3)
ao_coef_norm_ord_cosgtos(i,j) = d(j,2)
enddo
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_x, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_y, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_z, (ao_num, ao_num) ]
implicit none
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: c, overlap, overlap_x, overlap_y, overlap_z
complex*16 :: alpha, beta, A_center(3), B_center(3)
complex*16 :: overlap1, overlap_x1, overlap_y1, overlap_z1
complex*16 :: overlap2, overlap_x2, overlap_y2, overlap_z2
ao_overlap_cosgtos = 0.d0
ao_overlap_cosgtos_x = 0.d0
ao_overlap_cosgtos_y = 0.d0
ao_overlap_cosgtos_z = 0.d0
dim1 = 100
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, n, l, c &
!$OMP , overlap_x , overlap_y , overlap_z , overlap &
!$OMP , overlap_x1, overlap_y1, overlap_z1, overlap1 &
!$OMP , overlap_x2, overlap_y2, overlap_z2, overlap2 ) &
!$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1 &
!$OMP , ao_overlap_cosgtos_x, ao_overlap_cosgtos_y, ao_overlap_cosgtos_z, ao_overlap_cosgtos &
!$OMP , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos )
do j = 1, ao_num
A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0)
A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
power_A(1) = ao_power(j,1)
power_A(2) = ao_power(j,2)
power_A(3) = ao_power(j,3)
do i = 1, ao_num
B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0)
B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
power_B(1) = ao_power(i,1)
power_B(2) = ao_power(i,2)
power_B(3) = ao_power(i,3)
do n = 1, ao_prim_num(j)
alpha = ao_expo_ord_transp_cosgtos(n,j)
do l = 1, ao_prim_num(i)
c = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
beta = ao_expo_ord_transp_cosgtos(l,i)
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x1, overlap_y1, overlap_z1, overlap1, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, conjg(alpha), beta, power_A, power_B &
, overlap_x2, overlap_y2, overlap_z2, overlap2, dim1 )
overlap_x = 2.d0 * real( overlap_x1 + overlap_x2 )
overlap_y = 2.d0 * real( overlap_y1 + overlap_y2 )
overlap_z = 2.d0 * real( overlap_z1 + overlap_z2 )
overlap = 2.d0 * real( overlap1 + overlap2 )
ao_overlap_cosgtos(i,j) = ao_overlap_cosgtos(i,j) + c * overlap
if( isnan(ao_overlap_cosgtos(i,j)) ) then
print*,'i, j', i, j
print*,'l, n', l, n
print*,'c, overlap', c, overlap
print*, overlap_x, overlap_y, overlap_z
stop
endif
ao_overlap_cosgtos_x(i,j) = ao_overlap_cosgtos_x(i,j) + c * overlap_x
ao_overlap_cosgtos_y(i,j) = ao_overlap_cosgtos_y(i,j) + c * overlap_y
ao_overlap_cosgtos_z(i,j) = ao_overlap_cosgtos_z(i,j) + c * overlap_z
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

21
src/ao_one_e_ints/kin_ao_ints.irp.f
View File

@ -23,14 +23,13 @@
double precision :: A_center(3), B_center(3)
double precision :: d_a_2,d_2
if(use_cosgtos) then
!print*, 'use_cosgtos for ao_kinetic_integrals ?', use_cosgtos
if(use_cgtos) then
do j = 1, ao_num
do i = 1, ao_num
ao_deriv2_x(i,j) = ao_deriv2_cosgtos_x(i,j)
ao_deriv2_y(i,j) = ao_deriv2_cosgtos_y(i,j)
ao_deriv2_z(i,j) = ao_deriv2_cosgtos_z(i,j)
ao_deriv2_x(i,j) = ao_deriv2_cgtos_x(i,j)
ao_deriv2_y(i,j) = ao_deriv2_cgtos_y(i,j)
ao_deriv2_z(i,j) = ao_deriv2_cgtos_z(i,j)
enddo
enddo
@ -92,8 +91,8 @@
power_A(1) = power_A(1)-2
double precision :: deriv_tmp
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(1) +1.d0) * overlap_x0 &
+power_A(1) * (power_A(1)-1.d0) * d_a_2 &
deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(1)) +1.d0) * overlap_x0 &
+dble(power_A(1)) * (dble(power_A(1))-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_y0*overlap_z0
ao_deriv2_x(i,j) += c*deriv_tmp
@ -107,8 +106,8 @@
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,d_2,overlap_z,overlap,dim1)
power_A(2) = power_A(2)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(2) +1.d0 ) * overlap_y0 &
+power_A(2) * (power_A(2)-1.d0) * d_a_2 &
deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(2)) +1.d0 ) * overlap_y0 &
+dble(power_A(2)) * (dble(power_A(2))-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_z0
ao_deriv2_y(i,j) += c*deriv_tmp
@ -122,8 +121,8 @@
call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_y,overlap_z,d_2,overlap,dim1)
power_A(3) = power_A(3)-2
deriv_tmp = (-2.d0 * alpha * (2.d0 * power_A(3) +1.d0 ) * overlap_z0 &
+power_A(3) * (power_A(3)-1.d0) * d_a_2 &
deriv_tmp = (-2.d0 * alpha * (2.d0 * dble(power_A(3)) +1.d0 ) * overlap_z0 &
+dble(power_A(3)) * (dble(power_A(3))-1.d0) * d_a_2 &
+4.d0 * alpha * alpha * d_2 )*overlap_x0*overlap_y0
ao_deriv2_z(i,j) += c*deriv_tmp

212
src/ao_one_e_ints/one_e_Coul_integrals_cosgtos.irp.f → src/ao_one_e_ints/one_e_coul_integrals_cgtos.irp.f
View File

@ -1,53 +1,74 @@
! ---
BEGIN_PROVIDER [ double precision, ao_integrals_n_e_cosgtos, (ao_num, ao_num)]
BEGIN_PROVIDER [double precision, ao_integrals_n_e_cgtos, (ao_num, ao_num)]
BEGIN_DOC
!
! Nucleus-electron interaction, in the cosgtos |AO| basis set.
! Nucleus-electron interaction, in the cgtos |AO| basis set.
!
! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
!
END_DOC
implicit none
integer :: num_A, num_B, power_A(3), power_B(3)
integer :: i, j, k, l, n_pt_in, m
double precision :: c, Z, A_center(3), B_center(3), C_center(3)
complex*16 :: alpha, beta, c1, c2
integer :: power_A(3), power_B(3)
integer :: i, j, k, l, m, n, ii, jj
double precision :: c, Z, C_center(3)
complex*16 :: alpha, alpha_inv, A_center(3), phiA, KA2
complex*16 :: beta, beta_inv, B_center(3), phiB, KB2
complex*16 :: C1, C2, I1, I2
complex*16 :: NAI_pol_mult_cosgtos
complex*16 :: NAI_pol_mult_cgtos
ao_integrals_n_e_cosgtos = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE ( i, j, k, l, m, alpha, beta, A_center, B_center, C_center &
!$OMP , power_A, power_B, num_A, num_B, Z, c, c1, c2, n_pt_in ) &
!$OMP SHARED ( ao_num, ao_prim_num, ao_nucl, nucl_coord, ao_power, nucl_num, nucl_charge &
!$OMP , ao_expo_ord_transp_cosgtos, ao_coef_norm_ord_transp_cosgtos &
!$OMP , n_pt_max_integrals, ao_integrals_n_e_cosgtos )
n_pt_in = n_pt_max_integrals
ao_integrals_n_e_cgtos = 0.d0
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, m, n, ii, jj, C_center, Z, c, &
!$OMP alpha, alpha_inv, A_center, phiA, KA2, power_A, C1, I1, &
!$OMP beta, beta_inv, B_center, phiB, KB2, power_B, C2, I2) &
!$OMP SHARED (ao_num, ao_prim_num, ao_nucl, nucl_coord, &
!$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_pw_ord_transp, ao_expo_phase_ord_transp, &
!$OMP ao_integrals_n_e_cgtos)
!$OMP DO SCHEDULE (dynamic)
do j = 1, ao_num
num_A = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
A_center(1:3) = nucl_coord(num_A,1:3)
jj = ao_nucl(j)
power_A(1:3) = ao_power(j,1:3)
do i = 1, ao_num
num_B = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
B_center(1:3) = nucl_coord(num_B,1:3)
do l = 1, ao_prim_num(j)
alpha = ao_expo_ord_transp_cosgtos(l,j)
ii = ao_nucl(i)
power_B(1:3) = ao_power(i,1:3)
do m = 1, ao_prim_num(i)
beta = ao_expo_ord_transp_cosgtos(m,i)
do n = 1, ao_prim_num(j)
alpha = ao_expo_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha
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)
enddo
phiA = ao_expo_phase_ord_transp(4,n,j)
KA2 = ao_expo_pw_ord_transp(4,n,j) * ao_expo_pw_ord_transp(4,n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta
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)
enddo
phiB = ao_expo_phase_ord_transp(4,l,i)
KB2 = ao_expo_pw_ord_transp(4,l,i) * ao_expo_pw_ord_transp(4,l,i)
C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2))
C2 = zexp((0.d0, 1.d0) * ( phiA - phiB) - 0.25d0 * (conjg(alpha_inv) * KA2 + beta_inv * KB2))
c = 0.d0
do k = 1, nucl_num
@ -56,26 +77,15 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e_cosgtos, (ao_num, ao_num)]
C_center(1:3) = nucl_coord(k,1:3)
!print *, ' '
!print *, A_center, B_center, C_center, power_A, power_B
!print *, real(alpha), real(beta)
I1 = NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_max_integrals)
c1 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B &
, alpha, beta, C_center, n_pt_in )
!c2 = c1
c2 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B &
, conjg(alpha), beta, C_center, n_pt_in )
!print *, ' c1 = ', real(c1)
!print *, ' c2 = ', real(c2)
c = c - Z * 2.d0 * real(c1 + c2)
I2 = NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, conjg(alpha), beta, C_center, n_pt_max_integrals)
c = c - Z * 2.d0 * real(C1 * I1 + C2 * I2)
enddo
ao_integrals_n_e_cosgtos(i,j) = ao_integrals_n_e_cosgtos(i,j) &
+ ao_coef_norm_ord_transp_cosgtos(l,j) &
* ao_coef_norm_ord_transp_cosgtos(m,i) * c
ao_integrals_n_e_cgtos(i,j) += c * ao_coef_cgtos_norm_ord_transp(n,j) &
* ao_coef_cgtos_norm_ord_transp(l,i)
enddo
enddo
enddo
@ -88,11 +98,11 @@ END_PROVIDER
! ---
complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in)
complex*16 function NAI_pol_mult_cgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in)
BEGIN_DOC
!
! Computes the electron-nucleus attraction with two primitves cosgtos.
! Computes the electron-nucleus attraction with two primitves cgtos.
!
! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle`
!
@ -102,29 +112,38 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a
include 'utils/constants.include.F'
integer, intent(in) :: n_pt_in, power_A(3), power_B(3)
double precision, intent(in) :: C_center(3), A_center(3), B_center(3)
complex*16, intent(in) :: alpha, beta
double precision, intent(in) :: C_center(3)
complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3)
integer :: i, n_pt, n_pt_out
double precision :: dist, const_mod
complex*16 :: p, p_inv, rho, dist_integral, const, const_factor, coeff, factor
complex*16 :: accu, P_center(3)
double precision :: dist_AB, dist_AC
complex*16 :: p, p_inv, rho, dist, dist_integral, const, const_factor, coeff, factor
complex*16 :: P_center(3)
complex*16 :: d(0:n_pt_in)
complex*16, external :: V_n_e_cosgtos
complex*16, external :: V_n_e_cgtos
complex*16, external :: crint_2
complex*16, external :: crint_sum_2
if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. &
(A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then
dist_AB = 0.d0
dist_AC = 0.d0
do i = 1, 3
dist_AB += abs(A_center(i) - B_center(i))
dist_AC += abs(A_center(i) - C_center(i) * (1.d0, 0.d0))
enddo
if((dist_AB .gt. 1d-13) .or. (dist_AC .gt. 1d-13)) then
continue
else
NAI_pol_mult_cosgtos = V_n_e_cosgtos( power_A(1), power_A(2), power_A(3) &
, power_B(1), power_B(2), power_B(3) &
, alpha, beta )
NAI_pol_mult_cgtos = V_n_e_cgtos(power_A(1), power_A(2), power_A(3), &
power_B(1), power_B(2), power_B(3), &
alpha, beta)
return
endif
@ -133,7 +152,7 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a
p_inv = (1.d0, 0.d0) / p
rho = alpha * beta * p_inv
dist = 0.d0
dist = (0.d0, 0.d0)
dist_integral = (0.d0, 0.d0)
do i = 1, 3
P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
@ -144,47 +163,38 @@ complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, a
const_factor = dist * rho
const = p * dist_integral
const_mod = dsqrt(real(const_factor)*real(const_factor) + aimag(const_factor)*aimag(const_factor))
if(const_mod > 80.d0) then
NAI_pol_mult_cosgtos = (0.d0, 0.d0)
if(abs(const_factor) > 80.d0) then
NAI_pol_mult_cgtos = (0.d0, 0.d0)
return
endif
factor = zexp(-const_factor)
coeff = dtwo_pi * factor * p_inv
do i = 0, n_pt_in
d(i) = (0.d0, 0.d0)
enddo
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
NAI_pol_mult_cosgtos = coeff * crint_2(0, const)
NAI_pol_mult_cgtos = coeff * crint_2(0, const)
return
endif
d(0:n_pt_in) = (0.d0, 0.d0)
call give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, &
power_A, power_B, C_center, n_pt_in, d, n_pt_out)
if(n_pt_out < 0) then
NAI_pol_mult_cosgtos = (0.d0, 0.d0)
NAI_pol_mult_cgtos = (0.d0, 0.d0)
return
endif
!accu = (0.d0, 0.d0)
!do i = 0, n_pt_out, 2
! accu += crint_2(shiftr(i, 1), const) * d(i)
!enddo
accu = crint_sum_2(n_pt_out, const, d)
NAI_pol_mult_cosgtos = accu * coeff
NAI_pol_mult_cgtos = coeff * crint_sum_2(n_pt_out, const, d)
return
end
! ---
subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
, power_A, power_B, C_center, n_pt_in, d, n_pt_out)
subroutine give_cpolynomial_mult_center_one_e(A_center, B_center, alpha, beta, &
power_A, power_B, C_center, n_pt_in, d, n_pt_out)
BEGIN_DOC
! Returns the explicit polynomial in terms of the "t" variable of the following
@ -195,8 +205,8 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
implicit none
integer, intent(in) :: n_pt_in, power_A(3), power_B(3)
double precision, intent(in) :: A_center(3), B_center(3), C_center(3)
complex*16, intent(in) :: alpha, beta
double precision, intent(in) :: C_center(3)
complex*16, intent(in) :: alpha, beta, A_center(3), B_center(3)
integer, intent(out) :: n_pt_out
complex*16, intent(out) :: d(0:n_pt_in)
@ -241,7 +251,7 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
a_x = power_A(1)
b_x = power_B(1)
call I_x1_pol_mult_one_e_cosgtos(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in)
call I_x1_pol_mult_one_e_cgtos(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in)
if(n_pt1 < 0) then
n_pt_out = -1
@ -265,7 +275,7 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
a_y = power_A(2)
b_y = power_B(2)
call I_x1_pol_mult_one_e_cosgtos(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in)
call I_x1_pol_mult_one_e_cgtos(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in)
if(n_pt2 < 0) then
n_pt_out = -1
@ -289,7 +299,7 @@ subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
a_z = power_A(3)
b_z = power_B(3)
call I_x1_pol_mult_one_e_cosgtos(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in)
call I_x1_pol_mult_one_e_cgtos(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in)
if(n_pt3 < 0) then
n_pt_out = -1
@ -317,7 +327,7 @@ end
! ---
recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt_in)
recursive subroutine I_x1_pol_mult_one_e_cgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt_in)
BEGIN_DOC
! Recursive routine involved in the electron-nucleus potential
@ -351,7 +361,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
elseif((a == 0) .and. (c .ne. 0)) then
call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, n_pt_in)
call I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, d, nd, n_pt_in)
elseif(a == 1) then
@ -361,7 +371,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
Y(ix) = (0.d0, 0.d0)
enddo
call I_x2_pol_mult_one_e_cosgtos(c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
call I_x2_pol_mult_one_e_cgtos(c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
do ix = 0, nx
X(ix) *= dble(c)
@ -370,7 +380,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
call multiply_cpoly(X, nx, R2x, 2, d, nd)
ny = 0
call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, Y, ny, n_pt_in)
call I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, Y, ny, n_pt_in)
call multiply_cpoly(Y, ny, R1x, 2, d, nd)
else
@ -381,7 +391,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
Y(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(a-2, c, R1x, R1xp, R2x, X, nx, n_pt_in)
call I_x1_pol_mult_one_e_cgtos(a-2, c, R1x, R1xp, R2x, X, nx, n_pt_in)
do ix = 0, nx
X(ix) *= dble(a-1)
@ -393,7 +403,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
X(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(a-1, c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
call I_x1_pol_mult_one_e_cgtos(a-1, c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
do ix = 0, nx
X(ix) *= dble(c)
enddo
@ -401,7 +411,7 @@ recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_
call multiply_cpoly(X, nx, R2x, 2, d, nd)
ny = 0
call I_x1_pol_mult_one_e_cosgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in)
call I_x1_pol_mult_one_e_cgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in)
call multiply_cpoly(Y, ny, R1x, 2, d, nd)
endif
@ -410,7 +420,7 @@ end
! ---
recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
recursive subroutine I_x2_pol_mult_one_e_cgtos(c, R1x, R1xp, R2x, d, nd, dim)
BEGIN_DOC
! Recursive routine involved in the electron-nucleus potential
@ -447,7 +457,7 @@ recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
Y(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(0, c-2, R1x, R1xp, R2x, X, nx, dim)
call I_x1_pol_mult_one_e_cgtos(0, c-2, R1x, R1xp, R2x, X, nx, dim)
do ix = 0, nx
X(ix) *= dble(c-1)
@ -460,7 +470,7 @@ recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
Y(ix) = (0.d0, 0.d0)
enddo
call I_x1_pol_mult_one_e_cosgtos(0, c-1, R1x, R1xp, R2x, Y, ny, dim)
call I_x1_pol_mult_one_e_cgtos(0, c-1, R1x, R1xp, R2x, Y, ny, dim)
if(ny .ge. 0) then
call multiply_cpoly(Y, ny, R1xp, 2, d, nd)
@ -472,7 +482,7 @@ end
! ---
complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta)
complex*16 function V_n_e_cgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta)
BEGIN_DOC
! Primitve nuclear attraction between the two primitves centered on the same atom.
@ -488,26 +498,26 @@ complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta)
complex*16, intent(in) :: alpha, beta
double precision :: V_phi, V_theta
complex*16 :: V_r_cosgtos
complex*16 :: V_r_cgtos
if( (iand(a_x + b_x, 1) == 1) .or. &
(iand(a_y + b_y, 1) == 1) .or. &
(iand(a_z + b_z, 1) == 1) ) then
V_n_e_cosgtos = (0.d0, 0.d0)
V_n_e_cgtos = (0.d0, 0.d0)
else
V_n_e_cosgtos = V_r_cosgtos(a_x + b_x + a_y + b_y + a_z + b_z + 1, alpha + beta) &
* V_phi(a_x + b_x, a_y + b_y) &
* V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1)
V_n_e_cgtos = V_r_cgtos(a_x + b_x + a_y + b_y + a_z + b_z + 1, alpha + beta) &
* V_phi(a_x + b_x, a_y + b_y) &
* V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1)
endif
end
! ---
complex*16 function V_r_cosgtos(n, alpha)
complex*16 function V_r_cgtos(n, alpha)
BEGIN_DOC
! Computes the radial part of the nuclear attraction integral:
@ -525,9 +535,9 @@ complex*16 function V_r_cosgtos(n, alpha)
double precision :: fact
if(iand(n, 1) .eq. 1) then
V_r_cosgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1))
V_r_cgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1))
else
V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
V_r_cgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
endif
end

246
src/ao_one_e_ints/one_e_kin_integrals_cgtos.irp.f Normal file
View File

@ -0,0 +1,246 @@
! ---
BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_x, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_y, (ao_num, ao_num)]
&BEGIN_PROVIDER [double precision, ao_deriv2_cgtos_z, (ao_num, ao_num)]
implicit none
integer :: i, j, m, n, l, ii, jj, dim1, power_A(3), power_B(3)
double precision :: c, deriv_tmp
complex*16 :: alpha, alpha_inv, A_center(3), KA2, phiA, C1
complex*16 :: beta, beta_inv, B_center(3), KB2, phiB, C2
complex*16 :: overlap_x, overlap_y, overlap_z, overlap
complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1
complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2
complex*16 :: overlap_m2_1, overlap_p2_1
complex*16 :: overlap_m2_2, overlap_p2_2
complex*16 :: deriv_tmp_1, deriv_tmp_2
dim1 = 100
! -- Dummy call to provide everything
A_center(:) = (0.0d0, 0.d0)
B_center(:) = (1.0d0, 0.d0)
alpha = (1.0d0, 0.d0)
beta = (0.1d0, 0.d0)
power_A = 1
power_B = 0
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
! ---
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE(i, j, m, n, l, ii, jj, c, C1, C2, &
!$OMP A_center, power_A, alpha, alpha_inv, KA2, phiA, &
!$OMP B_center, power_B, beta, beta_inv, KB2, phiB, &
!$OMP deriv_tmp, deriv_tmp_1, deriv_tmp_2, &
!$OMP overlap_x, overlap_y, overlap_z, overlap, &
!$OMP overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2, &
!$OMP overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, &
!$OMP overlap_y0_2, overlap_z0_2) &
!$OMP SHARED(nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1, &
!$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_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z)
do j = 1, ao_num
jj = ao_nucl(j)
power_A(1) = ao_power(j,1)
power_A(2) = ao_power(j,2)
power_A(3) = ao_power(j,3)
do i = 1, ao_num
ii = ao_nucl(i)
power_B(1) = ao_power(i,1)
power_B(2) = ao_power(i,2)
power_B(3) = ao_power(i,3)
ao_deriv2_cgtos_x(i,j) = 0.d0
ao_deriv2_cgtos_y(i,j) = 0.d0
ao_deriv2_cgtos_z(i,j) = 0.d0
do n = 1, ao_prim_num(j)
alpha = ao_expo_cgtos_ord_transp(n,j)
alpha_inv = (1.d0, 0.d0) / alpha
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)
enddo
phiA = ao_expo_phase_ord_transp(4,n,j)
KA2 = ao_expo_pw_ord_transp(4,n,j) * ao_expo_pw_ord_transp(4,n,j)
do l = 1, ao_prim_num(i)
beta = ao_expo_cgtos_ord_transp(l,i)
beta_inv = (1.d0, 0.d0) / beta
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)
enddo
phiB = ao_expo_phase_ord_transp(4,l,i)
KB2 = ao_expo_pw_ord_transp(4,l,i) * ao_expo_pw_ord_transp(4,l,i)
c = ao_coef_cgtos_norm_ord_transp(n,j) * ao_coef_cgtos_norm_ord_transp(l,i)
C1 = zexp((0.d0, 1.d0) * (-phiA - phiB) - 0.25d0 * (alpha_inv * KA2 + beta_inv * KB2))
C2 = zexp((0.d0, 1.d0) * ( phiA - phiB) - 0.25d0 * (conjg(alpha_inv) * KA2 + beta_inv * KB2))
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1)
call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, &
overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1)
! ---
power_A(1) = power_A(1) - 2
if(power_A(1) > -1) then
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_m2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, &
overlap_m2_2, overlap_y, overlap_z, overlap, dim1)
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(1) = power_A(1) + 4
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_p2_1, overlap_y, overlap_z, overlap, dim1)
call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, &
overlap_p2_2, overlap_y, overlap_z, overlap, dim1)
power_A(1) = power_A(1) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(1)) + 1.d0) * overlap_x0_1 &
+ dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_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 &
+ dble(power_A(1)) * (dble(power_A(1)) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
ao_deriv2_cgtos_x(i,j) += c * deriv_tmp
! ---
power_A(2) = power_A(2) - 2
if(power_A(2) > -1) then
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_m2_1, overlap_y, overlap, dim1)
call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, &
overlap_x, overlap_m2_2, overlap_y, overlap, dim1)
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(2) = power_A(2) + 4
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_p2_1, overlap_y, overlap, dim1)
call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, &
overlap_x, overlap_p2_2, overlap_y, overlap, dim1)
power_A(2) = power_A(2) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(2)) + 1.d0) * overlap_y0_1 &
+ dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_1 &
+ 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 &
+ dble(power_A(2)) * (dble(power_A(2)) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
ao_deriv2_cgtos_y(i,j) += c * deriv_tmp
! ---
power_A(3) = power_A(3) - 2
if(power_A(3) > -1) then
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_y, overlap_m2_1, overlap, dim1)
call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, &
overlap_x, overlap_y, overlap_m2_2, overlap, dim1)
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(3) = power_A(3) + 4
call overlap_cgaussian_xyz(A_center, B_center, alpha, beta, power_A, power_B, &
overlap_x, overlap_y, overlap_p2_1, overlap, dim1)
call overlap_cgaussian_xyz(A_center, B_center, alpha, conjg(beta), power_A, power_B, &
overlap_x, overlap_y, overlap_p2_2, overlap, dim1)
power_A(3) = power_A(3) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * dble(power_A(3)) + 1.d0) * overlap_z0_1 &
+ dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_1 &
+ 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 &
+ dble(power_A(3)) * (dble(power_A(3)) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2
deriv_tmp = 2.d0 * real(C1 * deriv_tmp_1 + C2 * deriv_tmp_2)
ao_deriv2_cgtos_z(i,j) += c * deriv_tmp
! ---
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cgtos, (ao_num, ao_num)]
BEGIN_DOC
!
! Kinetic energy integrals in the cgtos |AO| basis.
!
! $\langle \chi_i |\hat{T}| \chi_j \rangle$
!
END_DOC
implicit none
integer :: i, j
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i, j) &
!$OMP SHARED(ao_num, ao_kinetic_integrals_cgtos, ao_deriv2_cgtos_x, ao_deriv2_cgtos_y, ao_deriv2_cgtos_z)
do j = 1, ao_num
do i = 1, ao_num
ao_kinetic_integrals_cgtos(i,j) = -0.5d0 * (ao_deriv2_cgtos_x(i,j) + &
ao_deriv2_cgtos_y(i,j) + &
ao_deriv2_cgtos_z(i,j))
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

223
src/ao_one_e_ints/one_e_kin_integrals_cosgtos.irp.f
View File

@ -1,223 +0,0 @@
! ---
BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_x, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_y, (ao_num, ao_num) ]
&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_z, (ao_num, ao_num) ]
implicit none
integer :: i, j, n, l, dim1, power_A(3), power_B(3)
double precision :: c, deriv_tmp
complex*16 :: alpha, beta, A_center(3), B_center(3)
complex*16 :: overlap_x, overlap_y, overlap_z, overlap
complex*16 :: overlap_x0_1, overlap_y0_1, overlap_z0_1
complex*16 :: overlap_x0_2, overlap_y0_2, overlap_z0_2
complex*16 :: overlap_m2_1, overlap_p2_1
complex*16 :: overlap_m2_2, overlap_p2_2
complex*16 :: deriv_tmp_1, deriv_tmp_2
dim1 = 100
! -- Dummy call to provide everything
A_center(:) = (0.0d0, 0.d0)
B_center(:) = (1.0d0, 0.d0)
alpha = (1.0d0, 0.d0)
beta = (0.1d0, 0.d0)
power_A = 1
power_B = 0
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
! ---
!$OMP PARALLEL DO SCHEDULE(GUIDED) &
!$OMP DEFAULT(NONE) &
!$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, l, n, c &
!$OMP , deriv_tmp, deriv_tmp_1, deriv_tmp_2 &
!$OMP , overlap_x, overlap_y, overlap_z, overlap &
!$OMP , overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2 &
!$OMP , overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, overlap_y0_2, overlap_z0_2 ) &
!$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1 &
!$OMP , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos &
!$OMP , ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z )
do j = 1, ao_num
A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0)
A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
power_A(1) = ao_power(j,1)
power_A(2) = ao_power(j,2)
power_A(3) = ao_power(j,3)
do i = 1, ao_num
B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0)
B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
power_B(1) = ao_power(i,1)
power_B(2) = ao_power(i,2)
power_B(3) = ao_power(i,3)
ao_deriv2_cosgtos_x(i,j) = 0.d0
ao_deriv2_cosgtos_y(i,j) = 0.d0
ao_deriv2_cosgtos_z(i,j) = 0.d0
do n = 1, ao_prim_num(j)
alpha = ao_expo_ord_transp_cosgtos(n,j)
do l = 1, ao_prim_num(i)
c = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
beta = ao_expo_ord_transp_cosgtos(l,i)
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1 )
! ---
power_A(1) = power_A(1) - 2
if(power_A(1) > -1) then
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_m2_1, overlap_y, overlap_z, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_m2_2, overlap_y, overlap_z, overlap, dim1 )
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(1) = power_A(1) + 4
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_p2_1, overlap_y, overlap_z, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_p2_2, overlap_y, overlap_z, overlap, dim1 )
power_A(1) = power_A(1) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_1 &
+ power_A(1) * (power_A(1) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_y0_1 * overlap_z0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_2 &
+ power_A(1) * (power_A(1) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
ao_deriv2_cosgtos_x(i,j) += c * deriv_tmp
! ---
power_A(2) = power_A(2) - 2
if(power_A(2) > -1) then
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_m2_1, overlap_y, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_m2_2, overlap_y, overlap, dim1 )
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(2) = power_A(2) + 4
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_p2_1, overlap_y, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_p2_2, overlap_y, overlap, dim1 )
power_A(2) = power_A(2) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_1 &
+ power_A(2) * (power_A(2) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_z0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_2 &
+ power_A(2) * (power_A(2) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
ao_deriv2_cosgtos_y(i,j) += c * deriv_tmp
! ---
power_A(3) = power_A(3) - 2
if(power_A(3) > -1) then
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_y, overlap_m2_1, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_y, overlap_m2_2, overlap, dim1 )
else
overlap_m2_1 = (0.d0, 0.d0)
overlap_m2_2 = (0.d0, 0.d0)
endif
power_A(3) = power_A(3) + 4
call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
, overlap_x, overlap_y, overlap_p2_1, overlap, dim1 )
call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
, overlap_x, overlap_y, overlap_p2_2, overlap, dim1 )
power_A(3) = power_A(3) - 2
deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_1 &
+ power_A(3) * (power_A(3) - 1.d0) * overlap_m2_1 &
+ 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_y0_1
deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_2 &
+ power_A(3) * (power_A(3) - 1.d0) * overlap_m2_2 &
+ 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2
deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
ao_deriv2_cosgtos_z(i,j) += c * deriv_tmp
! ---
enddo
enddo
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cosgtos, (ao_num, ao_num)]
BEGIN_DOC
!
! Kinetic energy integrals in the cosgtos |AO| basis.
!
! $\langle \chi_i |\hat{T}| \chi_j \rangle$
!
END_DOC
implicit none
integer :: i, j
!$OMP PARALLEL DO DEFAULT(NONE) &
!$OMP PRIVATE(i, j) &
!$OMP SHARED(ao_num, ao_kinetic_integrals_cosgtos, ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z)
do j = 1, ao_num
do i = 1, ao_num
ao_kinetic_integrals_cosgtos(i,j) = -0.5d0 * ( ao_deriv2_cosgtos_x(i,j) &
+ ao_deriv2_cosgtos_y(i,j) &
+ ao_deriv2_cosgtos_z(i,j) )
enddo
enddo
!$OMP END PARALLEL DO
END_PROVIDER
! ---

5
src/ao_one_e_ints/pot_ao_ints.irp.f
View File

@ -27,12 +27,11 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
else
if(use_cosgtos) then
!print *, " use_cosgtos for ao_integrals_n_e ?", use_cosgtos
if(use_cgtos) then
do j = 1, ao_num
do i = 1, ao_num
ao_integrals_n_e(i,j) = ao_integrals_n_e_cosgtos(i,j)
ao_integrals_n_e(i,j) = ao_integrals_n_e_cgtos(i,j)
enddo
enddo

6
src/ao_one_e_ints/pot_pt_charges.irp.f
View File

@ -25,12 +25,12 @@ BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)]
!
! else
! if(use_cosgtos) then
! !print *, " use_cosgtos for ao_integrals_pt_chrg ?", use_cosgtos
! if(use_cgtos) then
! !print *, " use_cgtos for ao_integrals_pt_chrg ?", use_cgtos
!
! do j = 1, ao_num
! do i = 1, ao_num
! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cosgtos(i,j)
! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cgtos(i,j)
! enddo
! enddo
!

1
src/ao_two_e_ints/EZFIO.cfg
View File

@ -23,6 +23,7 @@ doc: If | (ii|jj) | < `ao_cholesky_threshold` then (ii|jj) is zero
interface: ezfio,provider,ocaml
default: 1.e-12
[do_ao_cholesky]
type: logical
doc: Perform Cholesky decomposition of AO integrals

153
src/ao_two_e_ints/deb_2eint_cgtos.irp.f Normal file
View File

@ -0,0 +1,153 @@
! ---
subroutine deb_ao_2eint_cgtos(i, j, k, l)
BEGIN_DOC
! integral of the AO basis <ik|jl> or (ij|kl)
! i(r1) j(r1) 1/r12 k(r2) l(r2)
END_DOC
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: i, j, k, l
integer :: p, q, r, s
integer :: num_i, num_j, num_k, num_l, dim1, I_power(3), J_power(3), K_power(3), L_power(3)
integer :: iorder_p1(3), iorder_p2(3), iorder_q1(3), iorder_q2(3)
complex*16 :: I_center(3), J_center(3), K_center(3), L_center(3)
complex*16 :: expo1, expo2, expo3, expo4
complex*16 :: P1_center(3), pp1
complex*16 :: P2_center(3), pp2
complex*16 :: Q1_center(3), qq1
complex*16 :: Q2_center(3), qq2
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
if(num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k) then
!print*, ao_prim_num(i), ao_prim_num(j), ao_prim_num(k), ao_prim_num(l)
do p = 1, 3
I_power(p) = ao_power(i,p)
J_power(p) = ao_power(j,p)
K_power(p) = ao_power(k,p)
L_power(p) = ao_power(l,p)
I_center(p) = nucl_coord(num_i,p) * (1.d0, 0.d0)
J_center(p) = nucl_coord(num_j,p) * (1.d0, 0.d0)
K_center(p) = nucl_coord(num_k,p) * (1.d0, 0.d0)
L_center(p) = nucl_coord(num_l,p) * (1.d0, 0.d0)
enddo
do p = 1, ao_prim_num(i)
expo1 = ao_expo_cgtos_ord_transp(p,i)
!print*, "expo1 = ", expo1
!print*, "center1 = ", I_center
do q = 1, ao_prim_num(j)
expo2 = ao_expo_cgtos_ord_transp(q,j)
!print*, "expo2 = ", expo2
!print*, "center2 = ", J_center
pp1 = expo1 + expo2
P1_center(1:3) = (expo1 * I_center(1:3) + expo2 * J_center(1:3)) / pp1
iorder_p1(1:3) = I_power(1:3) + J_power(1:3)
pp2 = conjg(expo1) + expo2
P2_center(1:3) = (conjg(expo1) * I_center(1:3) + expo2 * J_center(1:3)) / pp2
iorder_p2(1:3) = I_power(1:3) + J_power(1:3)
do r = 1, ao_prim_num(k)
expo3 = ao_expo_cgtos_ord_transp(r,k)
!print*, "expo3 = ", expo3
!print*, "center3 = ", K_center
do s = 1, ao_prim_num(l)
expo4 = ao_expo_cgtos_ord_transp(s,l)
!print*, "expo4 = ", expo4
!print*, "center4 = ", L_center
qq1 = expo3 + expo4
Q1_center(1:3) = (expo3 * K_center(1:3) + expo4 * L_center(1:3)) / qq1
iorder_q1(1:3) = K_power(1:3) + L_power(1:3)
qq2 = conjg(expo3) + expo4
Q2_center(1:3) = (conjg(expo3) * K_center(1:3) + expo4 * L_center(1:3)) / qq2
iorder_q2(1:3) = K_power(1:3) + L_power(1:3)
call deb_cboys(P1_center, pp1, iorder_p1, Q1_center, qq1, iorder_q1)
call deb_cboys(P1_center, pp1, iorder_p1, Q2_center, qq2, iorder_q2)
call deb_cboys(P2_center, pp2, iorder_p2, Q1_center, qq1, iorder_q1)
call deb_cboys(P2_center, pp2, iorder_p2, Q2_center, qq2, iorder_q2)
call deb_cboys(conjg(P2_center), conjg(pp2), iorder_p2, Q1_center, qq1, iorder_q1)
call deb_cboys(conjg(P2_center), conjg(pp2), iorder_p2, Q2_center, qq2, iorder_q2)
call deb_cboys(conjg(P1_center), conjg(pp1), iorder_p1, Q1_center, qq1, iorder_q1)
call deb_cboys(conjg(P1_center), conjg(pp1), iorder_p1, Q2_center, qq2, iorder_q2)
enddo ! s
enddo ! r
enddo ! q
enddo ! p
endif ! same centers
return
end
! ---
subroutine deb_cboys(P_center, p, iorder_p, Q_center, q, iorder_q)
implicit none
include 'utils/constants.include.F'
integer, intent(in) :: iorder_p(3), iorder_q(3)
complex*16, intent(in) :: P_center(3), p
complex*16, intent(in) :: Q_center(3), q
integer :: iorder, n
complex*16 :: dist, rho
complex*16 :: int1, int2
complex*16, external :: crint_2
dist = (P_center(1) - Q_center(1)) * (P_center(1) - Q_center(1)) &
+ (P_center(2) - Q_center(2)) * (P_center(2) - Q_center(2)) &
+ (P_center(3) - Q_center(3)) * (P_center(3) - Q_center(3))
rho = dist * p * q / (p + q)
if(abs(rho) .lt. 1d-15) return
iorder = 2*iorder_p(1)+2*iorder_q(1) + 2*iorder_p(2)+2*iorder_q(2) + 2*iorder_p(3)+2*iorder_q(3)
n = shiftr(iorder, 1)
!write(33,*) n, real(rho), aimag(rho)
!print*, n, real(rho), aimag(rho)
int1 = crint_2(n, rho)
call crint_quad_12(n, rho, 1000000, int2)
if(abs(int1 - int2) .gt. 1d-5) then
print*, ' important error found: '
print*, p!, P_center
print*, q!, Q_center
print*, dist
print*, " n, tho = ", n, real(rho), aimag(rho)
print*, real(int1), real(int2), dabs(real(int1-int2))
print*, aimag(int1), aimag(int2), dabs(aimag(int1-int2))
stop
endif
end
! ---

2
src/ao_two_e_ints/screening.irp.f
View File

@ -3,7 +3,7 @@ logical function ao_two_e_integral_zero(i,j,k,l)
integer, intent(in) :: i,j,k,l
ao_two_e_integral_zero = .False.
if (.not.(read_ao_two_e_integrals.or.is_periodic)) then
if (.not.(read_ao_two_e_integrals.or.is_periodic.or.use_cgtos)) then
if (ao_overlap_abs(j,l)*ao_overlap_abs(i,k) < ao_integrals_threshold) then
ao_two_e_integral_zero = .True.
return

1504
src/ao_two_e_ints/two_e_Coul_integrals_cosgtos.irp.f
View File

File diff suppressed because it is too large Load Diff

1721
src/ao_two_e_ints/two_e_coul_integrals_cgtos.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

21
src/ao_two_e_ints/two_e_integrals.irp.f
View File

@ -38,15 +38,14 @@ double precision function ao_two_e_integral(i, j, k, l)
double precision :: Q_new(0:max_dim,3),Q_center(3),fact_q,qq
double precision, external :: ao_two_e_integral_erf
double precision, external :: ao_two_e_integral_cosgtos
double precision, external :: ao_two_e_integral_cgtos
double precision, external :: ao_two_e_integral_schwartz_accel
logical, external :: do_schwartz_accel
if(use_cosgtos) then
!print *, ' use_cosgtos for ao_two_e_integral ?', use_cosgtos
if(use_cgtos) then
ao_two_e_integral = ao_two_e_integral_cosgtos(i, j, k, l)
ao_two_e_integral = ao_two_e_integral_cgtos(i, j, k, l)
else if (use_only_lr) then
@ -54,17 +53,17 @@ double precision function ao_two_e_integral(i, j, k, l)
else if (do_schwartz_accel(i,j,k,l)) then
ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
ao_two_e_integral = ao_two_e_integral_schwartz_accel(i,j,k,l)
else
dim1 = n_pt_max_integrals
dim1 = n_pt_max_integrals
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_two_e_integral = 0.d0
num_i = ao_nucl(i)
num_j = ao_nucl(j)
num_k = ao_nucl(k)
num_l = ao_nucl(l)
ao_two_e_integral = 0.d0
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
do p = 1, 3

28
src/determinants/tr_density_matrix.irp.f
View File

@ -6,28 +6,28 @@ BEGIN_PROVIDER [double precision, one_e_tr_dm_mo, (mo_num, mo_num, N_states, N_s
! One body transition density matrix for all pairs of states n and m, < Psi^n | a_i^\dagger a_a | Psi^m >
END_DOC
integer :: j,k,l,m,k_a,k_b,n
integer :: j,l,m,k_a,k_b,n
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck, cl, ckl
double precision :: ck, ckl
double precision :: phase
integer :: h1,h2,p1,p2,s1,s2, degree
integer :: h1,h2,p1,p2,degree
integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int)
integer :: exc(0:2,2),n_occ(2)
double precision, allocatable :: tmp_a(:,:,:,:), tmp_b(:,:,:,:)
integer :: krow, kcol, lrow, lcol
PROVIDE psi_det
PROVIDE psi_det_alpha_unique psi_det_beta_unique
one_e_tr_dm_mo = 0d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,k_a,k_b,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,&
!$OMP PRIVATE(j,k_a,k_b,l,m,n,occ,ck, ckl,phase,h1,h2,p1,p2,degree,exc,&
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num, &
!$OMP SHARED(N_int,N_states,elec_alpha_num, &
!$OMP elec_beta_num,one_e_tr_dm_mo,N_det,&
!$OMP mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,&
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,&
!$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values,&
!$OMP N_det_alpha_unique,N_det_beta_unique,irp_here)
allocate(tmp_a(mo_num,mo_num,N_states,N_states), tmp_b(mo_num,mo_num,N_states,N_states) )
@ -163,28 +163,28 @@ END_PROVIDER
! $\alpha$ and $\beta$ one-body transition density matrices for all pairs of states
END_DOC
integer :: j,k,l,m,n,k_a,k_b
integer :: j,l,m,n,k_a,k_b
integer :: occ(N_int*bit_kind_size,2)
double precision :: ck, cl, ckl
double precision :: ck, ckl
double precision :: phase
integer :: h1,h2,p1,p2,s1,s2, degree
integer :: h1,h2,p1,p2,degree
integer(bit_kind) :: tmp_det(N_int,2), tmp_det2(N_int)
integer :: exc(0:2,2),n_occ(2)
double precision, allocatable :: tmp_a(:,:,:,:), tmp_b(:,:,:,:)
integer :: krow, kcol, lrow, lcol
PROVIDE psi_det
PROVIDE psi_det_alpha_unique psi_det_beta_unique
one_e_tr_dm_mo_alpha = 0.d0
one_e_tr_dm_mo_beta = 0.d0
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(j,k,k_a,k_b,l,m,n,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc,&
!$OMP PRIVATE(j,k_a,k_b,l,m,n,occ,ck, ckl,phase,h1,h2,p1,p2,degree,exc,&
!$OMP tmp_a, tmp_b, n_occ, krow, kcol, lrow, lcol, tmp_det, tmp_det2)&
!$OMP SHARED(psi_det,psi_coef,N_int,N_states,elec_alpha_num, &
!$OMP SHARED(N_int,N_states,elec_alpha_num, &
!$OMP elec_beta_num,one_e_tr_dm_mo_alpha,one_e_tr_dm_mo_beta,N_det,&
!$OMP mo_num,psi_bilinear_matrix_rows,psi_bilinear_matrix_columns,&
!$OMP psi_bilinear_matrix_transp_rows, psi_bilinear_matrix_transp_columns,&
!$OMP psi_bilinear_matrix_order_reverse, psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP psi_det_alpha_unique, psi_det_beta_unique,&
!$OMP psi_bilinear_matrix_values, psi_bilinear_matrix_transp_values,&
!$OMP N_det_alpha_unique,N_det_beta_unique,irp_here)
allocate(tmp_a(mo_num,mo_num,N_states,N_states), tmp_b(mo_num,mo_num,N_states,N_states) )

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

@ -1,55 +1,67 @@
program deb_ao_2e_int
!call check_ao_two_e_integral_cosgtos()
call check_crint1()
implicit none
call check_ao_two_e_integral_cgtos()
!call check_crint1()
!call check_crint2()
!call check_crint3()
end
! ---
subroutine check_ao_two_e_integral_cosgtos()
subroutine check_ao_two_e_integral_cgtos()
implicit none
integer :: i, j, k, l
double precision :: tmp1, tmp2
double precision :: acc, nrm, dif
double precision :: tmp1, tmp2
double precision :: pw, pw0
double precision :: t1, t2, tt
double precision, external :: ao_two_e_integral
double precision, external :: ao_two_e_integral_cosgtos
double precision, external :: ao_two_e_integral_cgtos
acc = 0.d0
nrm = 0.d0
i = 1
j = 6
k = 1
l = 16
! do i = 1, ao_num
! do k = 1, ao_num
! do j = 1, ao_num
! do l = 1, ao_num
pw0 = dble(ao_num**3)
pw = 0.d0
tt = 0.d0
do i = 1, ao_num
call wall_time(t1)
do j = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
tmp1 = ao_two_e_integral (i, j, k, l)
tmp2 = ao_two_e_integral_cosgtos(i, j, k, l)
call deb_ao_2eint_cgtos(i, j, k, l)
dif = dabs(tmp1 - tmp2)
if(dif .gt. 1d-12) then
print*, ' error on:', i, j, k, l
print*, tmp1, tmp2, dif
stop
endif
!tmp1 = ao_two_e_integral (i, j, k, l)
!tmp2 = ao_two_e_integral_cgtos(i, j, k, l)
! acc += dif
! nrm += dabs(tmp1)
! enddo
! enddo
! enddo
! enddo
!print*, i, j, k, l
print *, ' acc (%) = ', dif * 100.d0 / nrm
!dif = abs(tmp1 - tmp2)
!!if(dif .gt. 1d-10) then
! print*, ' error on:', i, j, k, l
! print*, tmp1, tmp2, dif
! !stop
!!endif
!acc += dif
!nrm += abs(tmp1)
enddo
enddo
enddo
call wall_time(t2)
tt += t2 - t1
print*, " % done = ", 100.d0 * dble(i) / ao_num
print*, ' ellepsed time (sec) =', tt
enddo
!print *, ' acc (%) = ', dif * 100.d0 / nrm
end
@ -73,8 +85,9 @@ subroutine check_crint1()
(+1d+4, +1d+4) /)
complex*16 :: rho
complex*16 :: int_an, int_nm
double precision, external :: rint
complex*16, external :: crint_1, crint_2, crint_quad
complex*16, external :: crint_1, crint_2
n = 10
dif_thr = 1d-7
@ -92,9 +105,9 @@ subroutine check_crint1()
acc_im = 0.d0
nrm_im = 0.d0
do i = 0, n
!int_an = crint_1 (i, rho)
int_an = crint_2 (i, rho)
int_nm = crint_quad(i, rho)
!int_an = crint_1(i, rho)
int_an = crint_2(i, rho)
call crint_quad_1(i, rho, 100000000, int_nm)
dif_re = dabs(real(int_an) - real(int_nm))
dif_im = dabs(aimag(int_an) - aimag(int_nm))
@ -185,4 +198,149 @@ end
! ---
subroutine check_crint3()
implicit none
integer :: i_test, n_test
integer :: nx, ny, n, n_quad
integer :: i, seed_size, clock_time
double precision :: xr(1:4), x
double precision :: yr(1:4), y
double precision :: dif_re, dif_im, acc_re, nrm_re, acc_im, nrm_im
double precision :: delta_ref
double precision :: t1, t2, t_int1, t_int2
complex*16 :: rho
complex*16 :: int1_old, int1_ref, int2_old, int2_ref
integer, allocatable :: seed(:)
complex*16, external :: 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 = 1
acc_re = 0.d0
nrm_re = 0.d0
acc_im = 0.d0
nrm_im = 0.d0
do i_test = 1, n_test
! Re(rho)
call random_number(xr)
x = xr(1)
if(xr(2) .gt. 0.5d0) x = -x
nx = int(15.d0 * xr(3))
if(xr(4) .gt. 0.5d0) nx = -nx
x = x * 10.d0**nx
! Im(rho)
call random_number(yr)
y = yr(1)
if(yr(2) .gt. 0.5d0) y = -y
ny = int(5.d0 * yr(3))
if(yr(4) .gt. 0.5d0) ny = -ny
y = y * 10.d0**ny
rho = x + (0.d0, 1.d0) * y
call random_number(x)
x = 31.d0 * x
n = int(x)
!if(n.eq.0) cycle
n = 0
!rho = (-6.83897018210218d0, -7.24479852507338d0)
rho = (-9.83206247355480d0, 0.445269582329036d0)
print*, " n = ", n
print*, " rho = ", real(rho), aimag(rho)
call wall_time(t1)
int1_old = crint_2(n, rho)
!n_quad = 10000000
!call crint_quad_1(n, rho, n_quad, int1_old)
!!delta_ref = 1.d0
!!do while(delta_ref .gt. 1d-12)
!! n_quad = n_quad * 10
!! !print*, " delta = ", delta_ref
!! !print*, " increasing n_quad to:", n_quad
!! call crint_quad_1(n, rho, n_quad, int1_ref)
!! delta_ref = abs(int1_ref - int1_old)
!! int1_old = int1_ref
!! if(n_quad .ge. 1000000000) then
!! print*, ' convergence was not reached for crint_quad_1'
!! print*, " delta = ", delta_ref
!! exit
!! endif
!!enddo
call wall_time(t2)
t_int1 = t_int1 + t2 - t1
!print*, " n_quad for crint_quad_1:", n_quad
call wall_time(t1)
n_quad = 10000000
call crint_quad_12(n, rho, n_quad, int2_old)
!delta_ref = 1.d0
!do while(delta_ref .gt. 1d-12)
! n_quad = n_quad * 10
! !print*, " delta = ", delta_ref
! !print*, " increasing n_quad to:", n_quad
! call crint_quad_12(n, rho, n_quad, int2_ref)
! delta_ref = abs(int2_ref - int2_old)
! int2_old = int2_ref
! if(n_quad .ge. 1000000000) then
! print*, ' convergence was not reached for crint_quad_2'
! print*, " delta = ", delta_ref
! exit
! endif
!enddo
call wall_time(t2)
t_int2 = t_int2 + t2 - t1
!print*, " n_quad for crint_quad_2:", n_quad
dif_re = dabs(real(int1_old) - real(int2_old))
dif_im = dabs(aimag(int1_old) - aimag(int2_old))
if((dif_re .gt. 1d-10) .or. (dif_im .gt. 1d-10)) then
print*, ' important error found: '
print*, " n = ", n
print*, " rho = ", real(rho), aimag(rho)
print*, real(int1_old), real(int2_old), dif_re
print*, aimag(int1_old), aimag(int2_old), dif_im
!stop
endif
if((real(int1_old) /= real(int1_old)) .or. (aimag(int1_old) /= aimag(int1_old)) .or. &
(real(int2_old) /= real(int2_old)) .or. (aimag(int2_old) /= aimag(int2_old)) ) then
cycle
else
acc_re += dif_re
acc_im += dif_im
nrm_re += dabs(real(int1_old))
nrm_im += dabs(aimag(int1_old))
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*, "crint_quad_1 wall time (sec) = ", t_int1
print*, "crint_quad_2 wall time (sec) = ", t_int2
deallocate(seed)
end
! ---

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

@ -101,7 +101,7 @@ subroutine give_explicit_cpoly_and_cgaussian(P_new, P_center, p, fact_k, iorder,
! IF fact_k is too smal then: returns a "s" function centered in zero
! with an inifinite exponent and a zero polynom coef
tmp_mod = dsqrt(REAL(fact_k)*REAL(fact_k) + AIMAG(fact_k)*AIMAG(fact_k))
tmp_mod = dsqrt(real(fact_k)*real(fact_k) + aimag(fact_k)*aimag(fact_k))
if(tmp_mod < 1d-14) then
iorder = 0
p = (1.d+14, 0.d0)
@ -145,68 +145,6 @@ end
! ---
!subroutine give_explicit_poly_and_gaussian_double(P_new,P_center,p,fact_k,iorder,alpha,beta,gama,a,b,A_center,B_center,Nucl_center,dim)
! BEGIN_DOC
! ! Transforms the product of
! ! (x-x_A)^a(1) (x-x_B)^b(1) (x-x_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) exp(-(r-Nucl_center)^2 gama
! !
! ! 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 )
! ! * [ 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_z = 0,i_order(3)) P_new(l_z,3) * (z-P_center(3))^l_z ] exp (- p (z-P_center(3))^2 )
! END_DOC
! implicit none
! include 'constants.include.F'
! integer, intent(in) :: dim
! integer, intent(in) :: a(3),b(3) ! powers : (x-xa)**a_x = (x-A(1))**a(1)
! double precision, intent(in) :: alpha, beta, gama ! exponents
! double precision, intent(in) :: A_center(3) ! A center
! double precision, intent(in) :: B_center (3) ! B center
! double precision, intent(in) :: Nucl_center(3) ! B center
! double precision, intent(out) :: P_center(3) ! new center
! double precision, intent(out) :: p ! new exponent
! double precision, intent(out) :: fact_k ! constant factor
! double precision, intent(out) :: P_new(0:max_dim,3)! polynomial
! integer , intent(out) :: iorder(3) ! i_order(i) = order of the polynomials
!
! double precision :: P_center_tmp(3) ! new center
! double precision :: p_tmp ! new exponent
! double precision :: fact_k_tmp,fact_k_bis ! constant factor
! double precision :: P_new_tmp(0:max_dim,3)! polynomial
! integer :: i,j
! double precision :: binom_func
!
! ! First you transform the two primitives into a sum of primitive with the same center P_center_tmp and gaussian exponent p_tmp
! call give_explicit_cpoly_and_cgaussian(P_new_tmp,P_center_tmp,p_tmp,fact_k_tmp,iorder,alpha,beta,a,b,A_center,B_center,dim)
! ! Then you create the new gaussian from the product of the new one per the Nuclei one
! call cgaussian_product(p_tmp,P_center_tmp,gama,Nucl_center,fact_k_bis,p,P_center)
! fact_k = fact_k_bis * fact_k_tmp
!
! ! Then you build the coefficient of the new polynom
! do i = 0, iorder(1)
! P_new(i,1) = 0.d0
! do j = i,iorder(1)
! P_new(i,1) = P_new(i,1) + P_new_tmp(j,1) * binom_func(j,j-i) * (P_center(1) - P_center_tmp(1))**(j-i)
! enddo
! enddo
! do i = 0, iorder(2)
! P_new(i,2) = 0.d0
! do j = i,iorder(2)
! P_new(i,2) = P_new(i,2) + P_new_tmp(j,2) * binom_func(j,j-i) * (P_center(2) - P_center_tmp(2))**(j-i)
! enddo
! enddo
! do i = 0, iorder(3)
! P_new(i,3) = 0.d0
! do j = i,iorder(3)
! P_new(i,3) = P_new(i,3) + P_new_tmp(j,3) * binom_func(j,j-i) * (P_center(3) - P_center_tmp(3))**(j-i)
! enddo
! enddo
!
!end
! ---
subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
BEGIN_DOC
@ -237,7 +175,7 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
ab = a * b * p_inv
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))
tmp_mod = dsqrt(real(k)*real(k) + aimag(k)*aimag(k))
if(tmp_mod .gt. 40.d0) then
k = (0.d0, 0.d0)
xp(1:3) = (0.d0, 0.d0)
@ -245,9 +183,9 @@ subroutine cgaussian_product(a, xa, b, xb, k, p, xp)
endif
k = zexp(-k)
xp(1) = ( a * xa(1) + b * xb(1) ) * p_inv
xp(2) = ( a * xa(2) + b * xb(2) ) * p_inv
xp(3) = ( a * xa(3) + b * xb(3) ) * p_inv
xp(1) = (a * xa(1) + b * xb(1)) * p_inv
xp(2) = (a * xa(2) + b * xb(2)) * p_inv
xp(3) = (a * xa(3) + b * xb(3)) * p_inv
end
@ -309,8 +247,6 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd)
integer, intent(out) :: nd
integer :: ndtmp, ib, ic
double precision :: tmp_mod
complex*16 :: tmp
if(ior(nc, nb) >= 0) then ! True if nc>=0 and nb>=0
continue
@ -332,9 +268,7 @@ subroutine multiply_cpoly(b, nb, c, nc, d, nd)
enddo
do nd = ndtmp, 0, -1
tmp = d(nd)
tmp_mod = dsqrt(REAL(tmp)*REAL(tmp) + AIMAG(tmp)*AIMAG(tmp))
if(tmp_mod .lt. 1.d-15) cycle
if(abs(d(nd)) .lt. 1.d-15) cycle
exit
enddo
@ -432,47 +366,42 @@ subroutine recentered_cpoly2(P_A, x_A, x_P, a, P_B, x_B, x_Q, b)
complex*16, intent(in) :: x_A, x_P, x_B, x_Q
complex*16, intent(out) :: P_A(0:a), P_B(0:b)
integer :: i, minab, maxab
complex*16 :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
integer :: i
integer :: maxbinom
complex*16 :: pows_a(0:a+b+2), pows_b(0:a+b+2)
double precision :: binom_func
if((a<0) .or. (b<0)) return
if((a < 0) .or. (b < 0)) return
maxab = max(a, b)
minab = max(min(a, b), 0)
maxbinom = size(binom_transp, 1)
pows_a(0) = (1.d0, 0.d0)
pows_a(1) = x_P - x_A
do i = 2, a
pows_a(i) = pows_a(i-1) * pows_a(1)
enddo
pows_b(0) = (1.d0, 0.d0)
pows_b(1) = x_Q - x_B
do i = 2, maxab
pows_a(i) = pows_a(i-1) * pows_a(1)
do i = 2, b
pows_b(i) = pows_b(i-1) * pows_b(1)
enddo
P_A(0) = pows_a(a)
do i = 1, min(a, maxbinom)
P_A(i) = binom_transp(i,a) * pows_a(a-i)
enddo
do i = maxbinom+1, a
P_A(i) = binom_func(a, i) * pows_a(a-i)
enddo
P_B(0) = pows_b(b)
do i = 1, min(minab, 20)
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
do i = 1, min(b, maxbinom)
P_B(i) = binom_transp(i,b) * pows_b(b-i)
enddo
do i = minab+1, min(a, 20)
P_A(i) = binom_transp(a-i,a) * pows_a(a-i)
enddo
do i = minab+1, min(b, 20)
P_B(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = 101, a
P_A(i) = binom_func(a,a-i) * pows_a(a-i)
enddo
do i = 101, b
P_B(i) = binom_func(b,b-i) * pows_b(b-i)
do i = maxbinom+1, b
P_B(i) = binom_func(b, i) * pows_b(b-i)
enddo
end

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

@ -16,7 +16,6 @@ complex*16 function crint_1(n, rho)
complex*16 :: sq_rho, rho_inv, rho_exp
complex*16 :: crint_smallz, cpx_erf_1
complex*16 :: cpx_erf_2
rho_re = real (rho)
rho_im = aimag(rho)
@ -48,10 +47,6 @@ complex*16 function crint_1(n, rho)
rho_exp = 0.5d0 * zexp(-rho)
rho_inv = (1.d0, 0.d0) / rho
!print*, sq_rho_re, sq_rho_im
!print*, cpx_erf_1(sq_rho_re, sq_rho_im)
!print*, cpx_erf_2(sq_rho_re, sq_rho_im)
crint_1 = 0.5d0 * sqpi * cpx_erf_1(sq_rho_re, sq_rho_im) / sq_rho
mmax = n
if(mmax .gt. 0) then
@ -68,66 +63,6 @@ end
! ---
complex*16 function crint_quad(n, rho)
implicit none
integer, intent(in) :: n
complex*16, intent(in) :: rho
integer :: i_quad, n_quad
double precision :: tmp_inv, tmp
n_quad = 1000000000
tmp_inv = 1.d0 / dble(n_quad)
!crint_quad = 0.5d0 * zexp(-rho)
!do i_quad = 1, n_quad - 1
! tmp = tmp_inv * dble(i_quad)
! tmp = tmp * tmp
! crint_quad += zexp(-rho*tmp) * tmp**n
!enddo
!crint_quad = crint_quad * tmp_inv
!crint_quad = 0.5d0 * zexp(-rho)
!do i_quad = 1, n_quad - 1
! tmp = tmp_inv * dble(i_quad)
! crint_quad += zexp(-rho*tmp) * tmp**n / dsqrt(tmp)
!enddo
!crint_quad = crint_quad * 0.5d0 * tmp_inv
! Composite Boole's Rule
crint_quad = 7.d0 * zexp(-rho)
do i_quad = 1, n_quad - 1
tmp = tmp_inv * dble(i_quad)
tmp = tmp * tmp
if(modulo(i_quad, 4) .eq. 0) then
crint_quad += 14.d0 * zexp(-rho*tmp) * tmp**n
else if(modulo(i_quad, 2) .eq. 0) then
crint_quad += 12.d0 * zexp(-rho*tmp) * tmp**n
else
crint_quad += 32.d0 * zexp(-rho*tmp) * tmp**n
endif
enddo
crint_quad = crint_quad * 2.d0 * tmp_inv / 45.d0
! Composite Simpson's 3/8 rule
!crint_quad = zexp(-rho)
!do i_quad = 1, n_quad - 1
! tmp = tmp_inv * dble(i_quad)
! tmp = tmp * tmp
! if(modulo(i_quad, 3) .eq. 0) then
! crint_quad += 2.d0 * zexp(-rho*tmp) * tmp**n
! else
! crint_quad += 3.d0 * zexp(-rho*tmp) * tmp**n
! endif
!enddo
!crint_quad = crint_quad * 3.d0 * tmp_inv / 8.d0
end
! ---
complex*16 function crint_sum_1(n_pt_out, rho, d1)
implicit none
@ -245,7 +180,7 @@ complex*16 function crint_smallz(n, rho)
complex*16, intent(in) :: rho
integer, parameter :: kmax = 40
double precision, parameter :: eps = 1.d-13
double precision, parameter :: eps = 1.d-10
integer :: k
double precision :: delta_mod
@ -267,7 +202,8 @@ complex*16 function crint_smallz(n, rho)
if(delta_mod > eps) then
write(*,*) ' pb in crint_smallz !'
write(*,*) ' n, rho = ', n, rho
write(*,*) ' n, rho = ', n, rho
write(*,*) ' value = ', crint_smallz
write(*,*) ' delta_mod = ', delta_mod
!stop 1
endif
@ -283,34 +219,99 @@ complex*16 function crint_2(n, rho)
integer, intent(in) :: n
complex*16, intent(in) :: rho
double precision :: tmp
complex*16 :: rho2
complex*16 :: vals(0:n)
complex*16, external :: crint_smallz
double precision :: tmp
complex*16 :: rho2
complex*16 :: vals(0:n)
if(abs(rho) < 10.d0) then
complex*16, external :: crint_smallz
if(abs(rho) < 3.5d0) then
if(abs(rho) .lt. 0.35d0) then
select case(n)
case(0)
crint_2 = (((((((((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_2 = (((((((((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_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
case(3)
crint_2 = (((((((((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_2 = (((((((((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(abs(rho) .lt. 1d-6) then
tmp = 2.d0 * dble(n)
rho2 = rho * rho
crint_2 = 1.d0 / (tmp + 1.d0) &
- rho / (tmp + 3.d0) &
+ 0.5d0 * rho2 / (tmp + 5.d0) &
- 0.16666666666666666d0 * rho * rho2 / (tmp + 7.d0)
else
crint_2 = crint_smallz(n, rho)
endif
else
if(real(rho) .ge. 0.d0) then
call zboysfun(n, rho, vals)
crint_2 = vals(n)
else
call zboysfunnrp(n, rho, vals)
crint_2 = vals(n) * zexp(-rho)
endif
endif
endif
return
@ -328,6 +329,9 @@ subroutine zboysfun(n_max, x, vals)
! Input: x --- argument, complex*16, Re(x) >= 0
! Output: vals --- values of the Boys function, n = 0, 1, ..., n_max
!
! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
! https://doi.org/10.1063/5.0062444
!
END_DOC
implicit none
@ -363,6 +367,9 @@ subroutine zboysfunnrp(n_max, x, vals)
! Input: x --- argument, complex *16 Re(x)<=0
! Output: vals --- values of e^z F(n,z), n = 0, 1, ..., n_max
!
! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
! https://doi.org/10.1063/5.0062444
!
END_DOC
implicit none
@ -393,27 +400,21 @@ complex*16 function crint_sum_2(n_pt_out, rho, d1)
integer, intent(in) :: n_pt_out
complex*16, intent(in) :: rho, d1(0:n_pt_out)
integer :: n, i
integer :: i
integer :: n_max
complex*16, allocatable :: vals(:)
!complex*16, external :: crint_2
!crint_sum_2 = (0.d0, 0.d0)
!do i = 0, n_pt_out, 2
! n = shiftr(i, 1)
! crint_sum_2 = crint_sum_2 + d1(i) * crint_2(n, rho)
!enddo
n_max = shiftr(n_pt_out, 1)
allocate(vals(0:n_max))
call crint_2_vec(n_max, rho, vals)
!call crint_2_vec(n_max, rho, vals)
! FOR DEBUG
call crint_quad_12_vec(n_max, rho, vals)
crint_sum_2 = d1(0) * vals(0)
do i = 2, n_pt_out, 2
n = shiftr(i, 1)
crint_sum_2 += d1(i) * vals(n)
crint_sum_2 += d1(i) * vals(shiftr(i, 1))
enddo
deallocate(vals)
@ -438,23 +439,82 @@ subroutine crint_2_vec(n_max, rho, vals)
abs_rho = abs(rho)
if(abs_rho < 10.d0) then
if(abs_rho < 3.5d0) then
if(abs_rho .lt. 1d-6) then
if(abs_rho .lt. 0.35d0) then
! use finite expansion for very small rho
vals(0) = (((((((((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
! rho^2 / 2
rho2 = 0.5d0 * rho * rho
! rho^3 / 6
rho3 = 0.3333333333333333d0 * rho * rho2
if(n > 0) then
vals(0) = 1.d0 - 0.3333333333333333d0 * rho + 0.2d0 * rho2 - 0.14285714285714285d0 * rho3
do n = 1, n_max
tmp = 2.d0 * dble(n)
vals(n) = 1.d0 / (tmp + 1.d0) - rho / (tmp + 3.d0) &
+ rho2 / (tmp + 5.d0) - rho3 / (tmp + 7.d0)
enddo
vals(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
if(n > 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 > 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
endif ! n > 2
endif ! n > 1
endif ! n > 0
else
@ -497,7 +557,7 @@ subroutine crint_smallz_vec(n_max, rho, vals)
complex*16, intent(out) :: vals(0:n_max)
integer, parameter :: kmax = 40
double precision, parameter :: eps = 1.d-13
double precision, parameter :: eps = 1.d-10
integer :: k, n
complex*16 :: ct, delta_k
@ -527,11 +587,12 @@ subroutine crint_smallz_vec(n_max, rho, vals)
endif
enddo
!if(abs(delta_k) > eps) then
! write(*,*) ' pb in crint_smallz_vec !'
! write(*,*) ' n, rho = ', n, rho
! write(*,*) ' |delta_k| = ', abs(delta_k)
!endif
if(abs(delta_k) > eps) then
write(*,*) ' pb in crint_smallz_vec !'
write(*,*) ' n, rho = ', n, rho
write(*,*) ' value = ', vals(n)
write(*,*) ' |delta_k| = ', abs(delta_k)
endif
enddo
deallocate(rho_k)
@ -541,3 +602,268 @@ end
! ---
subroutine crint_quad_1(n, rho, n_quad, crint_quad)
implicit none
integer, intent(in) :: n, n_quad
complex*16, intent(in) :: rho
complex*16, intent(out) :: crint_quad
integer :: i_quad
double precision :: tmp_inv, tmp0, tmp1, tmp2
double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /)
tmp_inv = 1.d0 / dble(n_quad)
crint_quad = 7.d0 * zexp(-rho)
tmp0 = 0.d0
select case (n)
case (0)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1)
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (1)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (2)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1 * tmp1
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (3)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1 * tmp1 * tmp1
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (4)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
tmp2 = tmp1 * tmp1
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp2 * tmp2
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case default
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1) * tmp1**n
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
end select
end
! ---
subroutine crint_quad_2(n, rho, n_quad, crint_quad)
implicit none
integer, intent(in) :: n, n_quad
complex*16, intent(in) :: rho
complex*16, intent(out) :: crint_quad
integer :: i_quad
double precision :: tmp_inv, tmp0, tmp1, tmp2
double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /)
complex*16 :: rhoc, rhoe
tmp_inv = 1.d0 / dble(n_quad)
crint_quad = 7.d0 * zexp(-rho)
tmp0 = 0.d0
rhoc = zexp(-rho*tmp_inv)
rhoe = (1.d0, 0.d0)
select case (n)
case (0)
!do i_quad = 1, n_quad - 1
! tmp0 = tmp0 + tmp_inv
! rhoe = rhoe * rhoc
! tmp1 = (rhoe - 1.d0) / dsqrt(tmp0)
! crint_quad = crint_quad + coef(iand(i_quad, 3)) * tmp1
!enddo
!crint_quad = 1.d0 + crint_quad * 0.022222222222222223d0 * tmp_inv
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe / dsqrt(tmp0)
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (1)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (2)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (3)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (4)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
tmp2 = tmp1 * tmp1 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp2
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case default
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0**n / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
end select
end
! ---
subroutine crint_quad_12(n, rho, n_quad, crint_quad)
implicit none
integer, intent(in) :: n, n_quad
complex*16, intent(in) :: rho
complex*16, intent(out) :: crint_quad
integer :: i_quad
double precision :: tmp_inv, tmp0, tmp1, tmp2
double precision :: coef(0:3) = (/14.d0, 32.d0, 12.d0, 32.d0 /)
complex*16 :: rhoc, rhoe
tmp_inv = 1.d0 / dble(n_quad)
crint_quad = 7.d0 * zexp(-rho)
tmp0 = 0.d0
rhoc = zexp(-rho*tmp_inv)
rhoe = (1.d0, 0.d0)
select case (n)
case (0)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
crint_quad = crint_quad + coef(iand(i_quad, 3)) * zexp(-rho*tmp1)
enddo
crint_quad = crint_quad * 0.044444444444444446d0 * tmp_inv
case (1)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (2)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (3)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0 * tmp0 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case (4)
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0 * tmp0
tmp2 = tmp1 * tmp1 / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp2
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
case default
do i_quad = 1, n_quad - 1
tmp0 = tmp0 + tmp_inv
tmp1 = tmp0**n / dsqrt(tmp0)
rhoe = rhoe * rhoc
crint_quad = crint_quad + coef(iand(i_quad, 3)) * rhoe * tmp1
enddo
crint_quad = crint_quad * 0.022222222222222223d0 * tmp_inv
end select
end
! ---
subroutine crint_quad_12_vec(n_max, rho, vals)
implicit none
integer, intent(in) :: n_max
complex*16, intent(in) :: rho
complex*16, intent(out) :: vals(0:n_max)
integer :: n
double precision :: tmp, abs_rho
complex*16 :: rho2, rho3, erho
do n = 0, n_max
call crint_quad_12(n, rho, 10000000, vals(n))
enddo
return
end
! ---

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

@ -143,7 +143,7 @@ complex*16 function erf_G(x, yabs)
idble = dble(i)
tmp0 = 0.5d0 * idble
tmp1 = tmp0 * tmp0 + x2
tmp2 = dexp( idble * yabs - tmp1 - dlog(tmp1) - dlog_2pi) * (x - (0.d0, 1.d0)*tmp0)
tmp2 = dexp(idble * yabs - tmp1 - dlog(tmp1) - dlog_2pi) * (x - (0.d0, 1.d0)*tmp0)
erf_G = erf_G + tmp2
@ -201,49 +201,6 @@ end
! ---
complex*16 function cpx_erf_2(x, y)
BEGIN_DOC
!
! compute erf(z) for z = x + i y
!
! Beylkin & Sharma, J. Chem. Phys. 155, 174117 (2021)
! https://doi.org/10.1063/5.0062444
!
END_DOC
implicit none
double precision, intent(in) :: x, y
double precision :: yabs
complex*16 :: z
yabs = dabs(y)
if(yabs .lt. 1.d-15) then
cpx_erf_2 = (1.d0, 0.d0) * derf(x)
return
else
z = x + (0.d0, 1.d0) * y
if(x .ge. 0.d0) then
call zboysfun00(z, cpx_erf_2)
else
call zboysfun00nrp(z, cpx_erf_2)
cpx_erf_2 = cpx_erf_2 * zexp(-z)
endif
endif
return
end
! ---
subroutine zboysfun00(z, val)
BEGIN_DOC
@ -333,20 +290,18 @@ subroutine zboysfun00(z, val)
complex*16, intent(out) :: val
integer :: k
complex*16 :: z1, zz, y
zz = zexp(-z)
complex*16 :: z1, y
if(abs(z) .ge. 100.0d0) then
! large |z|
z1 = 1.0d0 / zsqrt(z)
y = 1.0d0 / z
z1 = (1.d0, 0.d0) / zsqrt(z)
y = (1.d0, 0.d0) / z
val = asymcoef(7)
do k = 6, 1, -1
val = val * y + asymcoef(k)
enddo
val = zz * val * y + z1 * sqpio2
val = zexp(-z) * val * y + z1 * sqpio2
else if(abs(z) .le. 0.35d0) then
@ -359,12 +314,7 @@ subroutine zboysfun00(z, val)
else
! intermediate |z|
val = sqpio2 / zsqrt(z) - 0.5d0 * zz * sum(ff(1:22)/(z+pp(1:22)))
!val = (0.d0, 0.d0)
!do k = 1, 22
! val += ff(k) / (z + pp(k))
!enddo
!val = sqpio2 / zsqrt(z) - 0.5d0 * zz * val
val = sqpio2 / zsqrt(z) - 0.5d0 * zexp(-z) * sum(ff(1:22)/(z+pp(1:22)))
endif
@ -389,6 +339,8 @@ subroutine zboysfun00nrp(z, val)
!
END_DOC
include 'constants.include.F'
implicit none
double precision, parameter :: asymcoef(1:7) = (/ -0.499999999999999799d0, &
@ -413,7 +365,6 @@ subroutine zboysfun00nrp(z, val)
double precision, parameter :: tol = 1.0d-03
double precision, parameter :: sqpio2 = 0.886226925452758014d0 ! sqrt(pi)/2
double precision, parameter :: pi = 3.14159265358979324d0
double precision, parameter :: etmax = 25.7903399171930621d0
double precision, parameter :: etmax1 = 26.7903399171930621d0
complex*16, parameter :: ima = (0.d0, 1.d0)
@ -452,40 +403,40 @@ subroutine zboysfun00nrp(z, val)
0.03112676196932382d0, &
0.013576229705876844d0 /)
double precision, parameter :: qq (1:16) = (/ 0.0007243228510223928d0, &
0.01980651726441906d0, &
0.11641097769229371d0, &
0.38573968881461146d0, &
0.9414671037609641d0, &
1.8939510935716377d0, &
3.3275564293459383d0, &
5.280587297262129d0, &
7.730992222360452d0, &
10.590207725831563d0, &
13.706359477128965d0, &
16.876705473663804d0, &
19.867876155236257d0, &
22.441333930203022d0, &
24.380717439613566d0, &
25.51771075067431d0 /)
double precision, parameter :: qq(1:16) = (/ 0.0007243228510223928d0, &
0.01980651726441906d0, &
0.11641097769229371d0, &
0.38573968881461146d0, &
0.9414671037609641d0, &
1.8939510935716377d0, &
3.3275564293459383d0, &
5.280587297262129d0, &
7.730992222360452d0, &
10.590207725831563d0, &
13.706359477128965d0, &
16.876705473663804d0, &
19.867876155236257d0, &
22.441333930203022d0, &
24.380717439613566d0, &
25.51771075067431d0 /)
double precision, parameter :: qq1 (1:16) = (/ 0.0007524078957852004d0,&
0.020574499281252233d0, &
0.12092472113522865d0, &
0.40069643967765295d0, &
0.9779717449089211d0, &
1.9673875468969015d0, &
3.4565797939091802d0, &
5.485337886599723d0, &
8.030755321535683d0, &
11.000834641174064d0, &
14.237812708111456d0, &
17.531086359214406d0, &
20.6382373144543d0, &
23.31147887603379d0, &
25.326060444703632d0, &
26.507139770710722d0 /)
double precision, parameter :: qq1(1:16) = (/ 0.0007524078957852004d0,&
0.020574499281252233d0, &
0.12092472113522865d0, &
0.40069643967765295d0, &
0.9779717449089211d0, &
1.9673875468969015d0, &
3.4565797939091802d0, &
5.485337886599723d0, &
8.030755321535683d0, &
11.000834641174064d0, &
14.237812708111456d0, &
17.531086359214406d0, &
20.6382373144543d0, &
23.31147887603379d0, &
25.326060444703632d0, &
26.507139770710722d0 /)
double precision, parameter :: uu(1:16) = (/ 0.9992759394074501d0, &
0.9803883431758104d0, &
@ -532,8 +483,8 @@ subroutine zboysfun00nrp(z, val)
if(abs(z) .ge. 100.0d0) then
! large |z|
z1 = 1.0d0 / zsqrt(z)
y = 1.0d0 / z
z1 = (1.d0, 0.d0) / zsqrt(z)
y = (1.d0, 0.d0) / z
val = asymcoef(7)
do k = 6, 1, -1
val = val * y + asymcoef(k)
@ -560,13 +511,13 @@ subroutine zboysfun00nrp(z, val)
zsum = zsum + ww(k) * (zz - uu(k)) / (qq(k) + z)
else
q = z + qq(k)
p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
zsum = zsum + ww(k) * p *zz
p = q * (0.041666666666666664d0*q * (q * (0.2d0*q - 1.d0) + 4.d0) - 0.5d0) + 1.d0
zsum = zsum + ww(k) * p * zz
endif
enddo
zt = ima * sqrt(z / etmax)
zt = ima * zsqrt(z / etmax)
tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt))
val = sqrt(etmax) * zsum / sqrt(pi) + zz * tmp / sqrt(pi*z)
val = dsqrt(etmax) * zsum * inv_sq_pi + zz * tmp / zsqrt(pi*z)
else
zsum = (0.d0, 0.d0)
do k = 1, 16
@ -574,13 +525,14 @@ subroutine zboysfun00nrp(z, val)
zsum = zsum + ww(k) * (zz - uu1(k)) / (qq1(k) + z)
else
q = z + qq1(k)
p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
!p = 1.0d0 - 0.5d0*q + q*q/6.0d0 - q*q*q/24.0d0 + q*q*q*q/120.0d0
p = q * (0.041666666666666664d0*q * (q * (0.2d0*q - 1.d0) + 4.d0) - 0.5d0) + 1.d0
zsum = zsum + ww(k) * p * zz
endif
enddo
zt = ima * zsqrt(z / etmax1)
tmp = 0.5d0 * ima * log((1.0d0 - zt) / (1.0d0 + zt))
val = dsqrt(etmax1) * zsum / dsqrt(pi) + zz * tmp / zsqrt(pi*z)
val = dsqrt(etmax1) * zsum * inv_sq_pi + zz * tmp / zsqrt(pi*z)
endif
return

112
src/utils/integration.irp.f
View File

@ -1008,50 +1008,78 @@ subroutine recentered_poly2_v0(P_new, lda, x_A, LD_xA, x_P, a, n_points)
deallocate(pows_a)
!deallocate(pows_b)
end subroutine recentered_poly2_v0
subroutine recentered_poly2(P_new,x_A,x_P,a,P_new2,x_B,x_Q,b)
implicit none
BEGIN_DOC
! Recenter two polynomials
END_DOC
integer, intent(in) :: a,b
double precision, intent(in) :: x_A,x_P,x_B,x_Q
double precision, intent(out) :: P_new(0:a),P_new2(0:b)
double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
double precision :: binom_func
integer :: i,j,k,l, minab, maxab
if ((a<0).or.(b<0) ) return
maxab = max(a,b)
minab = max(min(a,b),0)
pows_a(0) = 1.d0
pows_a(1) = (x_P - x_A)
pows_b(0) = 1.d0
pows_b(1) = (x_Q - x_B)
do i = 2,maxab
pows_a(i) = pows_a(i-1)*pows_a(1)
pows_b(i) = pows_b(i-1)*pows_b(1)
enddo
P_new (0) = pows_a(a)
P_new2(0) = pows_b(b)
do i = 1,min(minab,20)
P_new (i) = binom_transp(a-i,a) * pows_a(a-i)
P_new2(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = minab+1,min(a,20)
P_new (i) = binom_transp(a-i,a) * pows_a(a-i)
enddo
do i = minab+1,min(b,20)
P_new2(i) = binom_transp(b-i,b) * pows_b(b-i)
enddo
do i = 101,a
P_new(i) = binom_func(a,a-i) * pows_a(a-i)
enddo
do i = 101,b
P_new2(i) = binom_func(b,b-i) * pows_b(b-i)
enddo
end
! ---
subroutine recentered_poly2(P_new, x_A, x_P, a, Q_new, x_B, x_Q, b)
BEGIN_DOC
!
! Recenter two polynomials:
!
! (x - x_A)^a -> \sum_{i=0}^{a} \binom{a}{i} (x_A - x_P)^{a-i} (x - x_P)^i
! (x - x_B)^b -> \sum_{i=0}^{b} \binom{b}{i} (x_B - x_Q)^{b-i} (x - x_Q)^i
!
! where:
! \binom{a}{i} = \binom{a}{a-i} = a! / [i! (a-i)!]
!
END_DOC
implicit none
integer, intent(in) :: a, b
double precision, intent(in) :: x_A, x_P, x_B, x_Q
double precision, intent(out) :: P_new(0:a), Q_new(0:b)
double precision :: pows_a(-2:a+b+4), pows_b(-2:a+b+4)
integer :: i, minab, maxab
integer :: maxbinom
double precision, external :: binom_func
if((a < 0) .or. (b < 0)) return
maxbinom = size(binom_transp, 1)
maxab = max(a, b)
minab = min(min(a, b), maxbinom)
pows_a(0) = 1.d0
pows_a(1) = x_P - x_A
pows_b(0) = 1.d0
pows_b(1) = x_Q - x_B
do i = 2, maxab
pows_a(i) = pows_a(i-1) * pows_a(1)
pows_b(i) = pows_b(i-1) * pows_b(1)
enddo
P_new(0) = pows_a(a)
Q_new(0) = pows_b(b)
do i = 1, minab
P_new(i) = binom_transp(i,a) * pows_a(a-i)
Q_new(i) = binom_transp(i,b) * pows_b(b-i)
enddo
do i = minab+1, min(a, maxbinom)
P_new(i) = binom_transp(i,a) * pows_a(a-i)
enddo
do i = minab+1, min(b, maxbinom)
Q_new(i) = binom_transp(i,b) * pows_b(b-i)
enddo
do i = maxbinom+1, a
P_new(i) = binom_func(a, i) * pows_a(a-i)
enddo
do i = maxbinom+1, b
Q_new(i) = binom_func(b, i) * pows_b(b-i)
enddo
end
! ---
subroutine pol_modif_center(A_center, B_center, iorder, A_pol, B_pol)
BEGIN_DOC