From 210179e8a04454c6a458e2b3bd114c82ce7215ce Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 6 Dec 2024 14:55:44 +0100 Subject: [PATCH] pot_ao_extra ints work --- plugins/local/ao_extra_basis/EZFIO.cfg | 18 ++-- .../local/ao_extra_basis/extra_basis.irp.f | 1 + plugins/local/extra_basis_int/NEED | 2 +- .../local/extra_basis_int/ao_overlap.irp.f | 53 +++++++++++ .../extra_basis_int/extra_basis_int.irp.f | 45 ++++++++- .../local/extra_basis_int/pot_ao_ints.irp.f | 93 +++++++++++++++++++ .../extra_basis_int/ref_extra_basis.irp.f | 35 +++++++ 7 files changed, 236 insertions(+), 11 deletions(-) create mode 100644 plugins/local/extra_basis_int/pot_ao_ints.irp.f create mode 100644 plugins/local/extra_basis_int/ref_extra_basis.irp.f diff --git a/plugins/local/ao_extra_basis/EZFIO.cfg b/plugins/local/ao_extra_basis/EZFIO.cfg index bddf4d37..8b3a8667 100644 --- a/plugins/local/ao_extra_basis/EZFIO.cfg +++ b/plugins/local/ao_extra_basis/EZFIO.cfg @@ -11,37 +11,37 @@ interface: ezfio, provider [ao_extra_prim_num] type: integer doc: Number of primitives per |ao_extra| -size: (extra_basis.ao_extra_num) +size: (ao_extra_basis.ao_extra_num) interface: ezfio, provider [ao_extra_prim_num_max] type: integer doc: Maximum number of primitives -default: =maxval(extra_basis.ao_extra_prim_num) +default: =maxval(ao_extra_basis.ao_extra_prim_num) interface: ezfio [ao_extra_nucl] type: integer doc: Index of the nucleus on which the |ao_extra| is centered -size: (extra_basis.ao_extra_num) +size: (ao_extra_basis.ao_extra_num) interface: ezfio, provider [ao_extra_power] type: integer doc: Powers of x, y and z for each |ao_extra| -size: (extra_basis.ao_extra_num,3) +size: (ao_extra_basis.ao_extra_num,3) interface: ezfio, provider [ao_extra_coef] type: double precision doc: Primitive coefficients, read from input. Those should not be used directly, as the MOs are expressed on the basis of **normalized** ao_extras. -size: (extra_basis.ao_extra_num,extra_basis.ao_extra_prim_num_max) +size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max) interface: ezfio, provider [ao_extra_expo] type: double precision doc: Exponents for each primitive of each |ao_extra| -size: (extra_basis.ao_extra_num,extra_basis.ao_extra_prim_num_max) +size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max) interface: ezfio, provider [ao_extra_md5] @@ -70,18 +70,18 @@ default: true [ao_extra_expo_im] type: double precision doc: imag part for Exponents for each primitive of each cGTOs |ao_extra| -size: (extra_basis.ao_extra_num,extra_basis.ao_extra_prim_num_max) +size: (ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max) interface: ezfio, provider [ao_extra_expo_pw] type: double precision doc: plane wave part for each primitive GTOs |ao_extra| -size: (3,extra_basis.ao_extra_num,extra_basis.ao_extra_prim_num_max) +size: (3,ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max) interface: ezfio, provider [ao_extra_expo_phase] type: double precision doc: phase shift for each primitive GTOs |ao_extra| -size: (3,extra_basis.ao_extra_num,extra_basis.ao_extra_prim_num_max) +size: (3,ao_extra_basis.ao_extra_num,ao_extra_basis.ao_extra_prim_num_max) interface: ezfio, provider diff --git a/plugins/local/ao_extra_basis/extra_basis.irp.f b/plugins/local/ao_extra_basis/extra_basis.irp.f index 19b36349..9773af2f 100644 --- a/plugins/local/ao_extra_basis/extra_basis.irp.f +++ b/plugins/local/ao_extra_basis/extra_basis.irp.f @@ -12,4 +12,5 @@ program extra_basis print*,'extra_nucl_coord = ' print*,extra_nucl_coord(i,1:3) enddo + print*,ao_extra_num end diff --git a/plugins/local/extra_basis_int/NEED b/plugins/local/extra_basis_int/NEED index 1dc12e7b..7f39af7c 100644 --- a/plugins/local/extra_basis_int/NEED +++ b/plugins/local/extra_basis_int/NEED @@ -1,2 +1,2 @@ -extra_basis +ao_extra_basis ao_one_e_ints diff --git a/plugins/local/extra_basis_int/ao_overlap.irp.f b/plugins/local/extra_basis_int/ao_overlap.irp.f index a821e161..4f8debb6 100644 --- a/plugins/local/extra_basis_int/ao_overlap.irp.f +++ b/plugins/local/extra_basis_int/ao_overlap.irp.f @@ -134,3 +134,56 @@ END_PROVIDER ! --- +subroutine get_ao_mixed_overlap(r_nucl,ao_mixed_overlap) + implicit none + BEGIN_DOC +! returns the overlap integrals between the AOs and the extra_AOs located at r_nucl + END_DOC + double precision, intent(in) :: r_nucl(extra_nucl_num,3) + double precision, intent(out):: ao_mixed_overlap(ao_extra_num,ao_num) + integer :: j,i,l,n, power_A(3), power_B(3), dim1 + double precision :: A_center(3), B_center(3), alpha, beta + double precision :: overlap_x,overlap_y,overlap_z,overlap,c + dim1=100 + ao_mixed_overlap = 0.d0 + +!$OMP PARALLEL DO SCHEDULE(GUIDED) & +!$OMP DEFAULT(NONE) & +!$OMP PRIVATE(A_center,B_center,power_A,power_B,& +!$OMP overlap_x,overlap_y, overlap_z, overlap, & +!$OMP alpha, beta,i,j,n,l,c) & +!$OMP SHARED(r_nucl,ao_extra_power,ao_extra_prim_num, & +!$OMP ao_mixed_overlap,ao_extra_num,ao_extra_coef_normalized_ordered_transp,ao_extra_nucl, & +!$OMP ao_extra_expo_ordered_transp,dim1, & +!$OMP nucl_coord,ao_power,ao_prim_num, & +!$OMP ao_num,ao_coef_normalized_ordered_transp,ao_nucl, & +!$OMP ao_expo_ordered_transp) + do i = 1, ao_num + B_center(1) = nucl_coord( ao_nucl(i), 1 ) + B_center(2) = nucl_coord( ao_nucl(i), 2 ) + B_center(3) = nucl_coord( ao_nucl(i), 3 ) + power_B(1) = ao_power( i, 1 ) + power_B(2) = ao_power( i, 2 ) + power_B(3) = ao_power( i, 3 ) + do l = 1, ao_prim_num(i) + beta = ao_expo_ordered_transp(l,i) + do j=1,ao_extra_num + A_center(1) = r_nucl( ao_extra_nucl(j), 1 ) + A_center(2) = r_nucl( ao_extra_nucl(j), 2 ) + A_center(3) = r_nucl( ao_extra_nucl(j), 3 ) + power_A(1) = ao_extra_power( j, 1 ) + power_A(2) = ao_extra_power( j, 2 ) + power_A(3) = ao_extra_power( j, 3 ) + do n = 1,ao_extra_prim_num(j) + alpha = ao_extra_expo_ordered_transp(n,j) + call overlap_gaussian_xyz(A_center,B_center,alpha,beta,power_A,power_B,overlap_x,overlap_y,overlap_z,overlap,dim1) + c = ao_extra_coef_normalized_ordered_transp(n,j) * ao_coef_normalized_ordered_transp(l,i) + ao_mixed_overlap(j,i) += c * overlap + enddo + enddo + enddo + enddo +!$OMP END PARALLEL DO + + +end diff --git a/plugins/local/extra_basis_int/extra_basis_int.irp.f b/plugins/local/extra_basis_int/extra_basis_int.irp.f index 9c1cf351..09838c4e 100644 --- a/plugins/local/extra_basis_int/extra_basis_int.irp.f +++ b/plugins/local/extra_basis_int/extra_basis_int.irp.f @@ -3,5 +3,48 @@ program extra_basis_int BEGIN_DOC ! TODO : Put the documentation of the program here END_DOC - print *, 'Hello world' +! call test_overlap + call routine_test_pot_ne +end + +subroutine test_overlap + implicit none + integer :: i,j + do i = 1, ao_extra_num + do j = 1, ao_extra_num + write(33,*)ao_extra_overlap(j,i) + enddo + enddo +end + +subroutine test_overlap_mixed + implicit none + integer :: i,j + double precision, allocatable :: ao_mixed_overlap(:,:) + allocate(ao_mixed_overlap(ao_extra_num,ao_num)) + call get_ao_mixed_overlap(extra_nucl_coord,ao_mixed_overlap) + do i = 1, ao_extra_num + do j = 1, ao_num + write(33,*)dabs(ao_extra_overlap_mixed(j,i)-ao_mixed_overlap(i,j)) + write(*,*)ao_extra_overlap_mixed(j,i),ao_mixed_overlap(i,j),dabs(ao_extra_overlap_mixed(j,i)-ao_mixed_overlap(i,j)) + enddo + enddo +end + +subroutine routine_test_pot_ne + implicit none + integer :: i,j + double precision :: integral, C_center(3), mu_in + double precision :: NAI_pol_mult_erf_ao_extra + C_center(1) = 0.1d0 + C_center(2) = -0.3d0 + C_center(3) = 0.8d0 + mu_in = 1.d10 + do i = 1, ao_extra_num + do j = 1, ao_extra_num + integral = NAI_pol_mult_erf_ao_extra(i, j, mu_in, C_center) + write(33,*)j,i,integral + enddo + enddo + end diff --git a/plugins/local/extra_basis_int/pot_ao_ints.irp.f b/plugins/local/extra_basis_int/pot_ao_ints.irp.f new file mode 100644 index 00000000..5f3af244 --- /dev/null +++ b/plugins/local/extra_basis_int/pot_ao_ints.irp.f @@ -0,0 +1,93 @@ + +double precision function NAI_pol_mult_erf_ao_extra(i_ao, j_ao, mu_in, C_center) + + BEGIN_DOC + ! + ! Computes the following integral : + ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + ! + ! where $\chi_i(r)$ AND $\chi_j(r)$ belongs to the extra basis + END_DOC + + implicit none + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: mu_in, C_center(3) + + integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in + double precision :: A_center(3), B_center(3), integral, alpha, beta + + double precision :: NAI_pol_mult_erf + + num_A = ao_extra_nucl(i_ao) + power_A(1:3) = ao_extra_power(i_ao,1:3) + A_center(1:3) = extra_nucl_coord(num_A,1:3) + num_B = ao_extra_nucl(j_ao) + power_B(1:3) = ao_extra_power(j_ao,1:3) + B_center(1:3) = extra_nucl_coord(num_B,1:3) + + n_pt_in = n_pt_max_extra_basis_integrals + + NAI_pol_mult_erf_ao_extra = 0.d0 + do i = 1, ao_extra_prim_num(i_ao) + alpha = ao_extra_expo_ordered_transp(i,i_ao) + do j = 1, ao_extra_prim_num(j_ao) + beta = ao_extra_expo_ordered_transp(j,j_ao) + + integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in,mu_in) + + NAI_pol_mult_erf_ao_extra += integral * ao_extra_coef_normalized_ordered_transp(j,j_ao) * ao_extra_coef_normalized_ordered_transp(i,i_ao) + enddo + enddo + +end function NAI_pol_mult_erf_ao_extra + +! --- + +double precision function NAI_pol_mult_erf_ao_extra_mixed(i_ao, j_ao, mu_in, C_center) + + BEGIN_DOC + ! + ! Computes the following integral : + ! $\int_{-\infty}^{infty} dr \chi_i(r) \chi_j(r) \frac{\erf(\mu |r - R_C|)}{|r - R_C|}$. + ! + ! + ! where $\chi_i(r)$ belongs to the extra basis and $\chi_j(r)$ to the regular basis + END_DOC + + implicit none + integer, intent(in) :: i_ao, j_ao + double precision, intent(in) :: mu_in, C_center(3) + + integer :: i, j, num_A, num_B, power_A(3), power_B(3), n_pt_in + double precision :: A_center(3), B_center(3), integral, alpha, beta + + double precision :: NAI_pol_mult_erf + + ! A = chi_i == extra basis + num_A = ao_extra_nucl(i_ao) + power_A(1:3) = ao_extra_power(i_ao,1:3) + A_center(1:3) = extra_nucl_coord(num_A,1:3) + ! B = chi_j == regular basis + num_B = ao_nucl(j_ao) + power_B(1:3) = ao_power(j_ao,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + n_pt_in = max(n_pt_max_integrals,n_pt_max_extra_basis_integrals) + + NAI_pol_mult_erf_ao_extra_mixed = 0.d0 + do i = 1, ao_extra_prim_num(i_ao) + alpha = ao_extra_expo_ordered_transp(i,i_ao) + do j = 1, ao_prim_num(j_ao) + beta = ao_expo_ordered_transp(j,j_ao) + + integral = NAI_pol_mult_erf(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in,mu_in) + + NAI_pol_mult_erf_ao_extra_mixed += integral * ao_coef_normalized_ordered_transp(j,j_ao) * ao_extra_coef_normalized_ordered_transp(i,i_ao) + enddo + enddo + +end + +! --- + diff --git a/plugins/local/extra_basis_int/ref_extra_basis.irp.f b/plugins/local/extra_basis_int/ref_extra_basis.irp.f new file mode 100644 index 00000000..8a50be2c --- /dev/null +++ b/plugins/local/extra_basis_int/ref_extra_basis.irp.f @@ -0,0 +1,35 @@ +program pouet + implicit none +! call ref_overlap + call ref_pot + +end + +subroutine ref_overlap + implicit none + integer :: i,j + do i = 1, ao_num + do j = 1, ao_num + write(34,*)ao_overlap(j,i) + enddo + enddo + +end + +subroutine ref_pot + implicit none + integer :: i,j + double precision :: integral, C_center(3), mu_in + double precision :: NAI_pol_mult_erf_ao + C_center(1) = 0.1d0 + C_center(2) = -0.3d0 + C_center(3) = 0.8d0 + mu_in = 1.d10 + do i = 1, ao_num + do j = 1, ao_num + integral = NAI_pol_mult_erf_ao(i, j, mu_in, C_center) + write(34,*)j,i,integral + enddo + enddo + +end