From 766fabf1d2f8b6c6bf65537fc61df88a82da22b5 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 7 Feb 2023 12:50:33 +0100 Subject: [PATCH 1/3] point charges work --- external/qp2-dependencies | 2 +- src/ao_one_e_ints/pot_ao_ints.irp.f | 4 + src/ao_one_e_ints/pot_pt_charges.irp.f | 108 +++++++++++++ src/nuclei/EZFIO.cfg | 24 +++ src/nuclei/nuclei.irp.f | 3 + src/nuclei/point_charges.irp.f | 209 +++++++++++++++++++++++++ 6 files changed, 349 insertions(+), 1 deletion(-) create mode 100644 src/ao_one_e_ints/pot_pt_charges.irp.f create mode 100644 src/nuclei/point_charges.irp.f diff --git a/external/qp2-dependencies b/external/qp2-dependencies index 242151e0..f40bde09 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 242151e03d1d6bf042387226431d82d35845686a +Subproject commit f40bde0925808bbec0424b57bfcef1b26473a1c8 diff --git a/src/ao_one_e_ints/pot_ao_ints.irp.f b/src/ao_one_e_ints/pot_ao_ints.irp.f index 4108ce71..1d92dc7d 100644 --- a/src/ao_one_e_ints/pot_ao_ints.irp.f +++ b/src/ao_one_e_ints/pot_ao_ints.irp.f @@ -80,6 +80,10 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)] IF (DO_PSEUDO) THEN ao_integrals_n_e += ao_pseudo_integrals ENDIF + IF(point_charges) THEN + ao_integrals_n_e += ao_integrals_pt_chrg + ENDIF + endif diff --git a/src/ao_one_e_ints/pot_pt_charges.irp.f b/src/ao_one_e_ints/pot_pt_charges.irp.f new file mode 100644 index 00000000..93f1acff --- /dev/null +++ b/src/ao_one_e_ints/pot_pt_charges.irp.f @@ -0,0 +1,108 @@ + +BEGIN_PROVIDER [ double precision, ao_integrals_pt_chrg, (ao_num,ao_num)] + + BEGIN_DOC + ! Point charge-electron interaction, in the |AO| basis set. + ! + ! :math:`\langle \chi_i | -\sum_charge charge * \frac{1}{|r-R_charge|} | \chi_j \rangle` + ! + ! Notice the minus sign convention as it is supposed to be for electrons. + END_DOC + + implicit none + integer :: num_A, num_B, power_A(3), power_B(3) + integer :: i, j, k, l, n_pt_in, m + double precision :: alpha, beta + double precision :: A_center(3),B_center(3),C_center(3) + double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult + + ao_integrals_pt_chrg = 0.d0 + +! if (read_ao_integrals_pt_chrg) then +! +! call ezfio_get_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) +! print *, 'AO N-e integrals read from disk' +! +! else + +! if(use_cosgtos) then +! !print *, " use_cosgtos for ao_integrals_pt_chrg ?", use_cosgtos +! +! do j = 1, ao_num +! do i = 1, ao_num +! ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg_cosgtos(i,j) +! enddo +! enddo +! +! else + + !$OMP PARALLEL & + !$OMP DEFAULT (NONE) & + !$OMP PRIVATE (i,j,k,l,m,alpha,beta,A_center,B_center,C_center,power_A,power_B,& + !$OMP num_A,num_B,Z,c,c1,n_pt_in) & + !$OMP SHARED (ao_num,ao_prim_num,ao_expo_ordered_transp,ao_power,ao_nucl,pts_charge_coord,ao_coef_normalized_ordered_transp,nucl_coord,& + !$OMP n_pt_max_integrals,ao_integrals_pt_chrg,n_pts_charge,pts_charge_z) + + n_pt_in = n_pt_max_integrals + + !$OMP DO SCHEDULE (dynamic) + + do j = 1, ao_num + num_A = ao_nucl(j) + power_A(1:3)= ao_power(j,1:3) + A_center(1:3) = nucl_coord(num_A,1:3) + + do i = 1, ao_num + + num_B = ao_nucl(i) + power_B(1:3)= ao_power(i,1:3) + B_center(1:3) = nucl_coord(num_B,1:3) + + do l=1,ao_prim_num(j) + alpha = ao_expo_ordered_transp(l,j) + + do m=1,ao_prim_num(i) + beta = ao_expo_ordered_transp(m,i) + + double precision :: c, c1 + c = 0.d0 + + do k = 1, n_pts_charge + double precision :: Z + Z = pts_charge_z(k) + + C_center(1:3) = pts_charge_coord(k,1:3) + + c1 = NAI_pol_mult( A_center, B_center, power_A, power_B & + , alpha, beta, C_center, n_pt_in ) + + c = c - Z * c1 + + enddo + ao_integrals_pt_chrg(i,j) = ao_integrals_pt_chrg(i,j) & + + ao_coef_normalized_ordered_transp(l,j) & + * ao_coef_normalized_ordered_transp(m,i) * c + enddo + enddo + enddo + enddo + + !$OMP END DO + !$OMP END PARALLEL + +! endif + + +! IF(do_pseudo) THEN +! ao_integrals_pt_chrg += ao_pseudo_integrals +! ENDIF + +! endif + + +! if (write_ao_integrals_pt_chrg) then +! call ezfio_set_ao_one_e_ints_ao_integrals_pt_chrg(ao_integrals_pt_chrg) +! print *, 'AO N-e integrals written to disk' +! endif + +END_PROVIDER diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index 34c27c46..060eede6 100644 --- a/src/nuclei/EZFIO.cfg +++ b/src/nuclei/EZFIO.cfg @@ -37,3 +37,27 @@ type: logical doc: If true, the calculation uses periodic boundary conditions interface: ezfio, provider, ocaml default: false +[n_pts_charge] +type: integer +doc: Number of point charges to be added to the potential +interface: ezfio +default: 0 + +[pts_charge_z] +type: double precision +doc: Charge associated to each point charge +interface: ezfio +size: (nuclei.n_pts_charge) + +[pts_charge_coord] +type: double precision +doc: Coordinate of each point charge. +interface: ezfio +size: (nuclei.n_pts_charge,3) + +[point_charges] +type: logical +doc: If |true|, point charges (see nuclei/write_pt_charges.py) are added to the one-electron potential +interface: ezfio,provider,ocaml +default: False + diff --git a/src/nuclei/nuclei.irp.f b/src/nuclei/nuclei.irp.f index c1b5f52f..3c04316f 100644 --- a/src/nuclei/nuclei.irp.f +++ b/src/nuclei/nuclei.irp.f @@ -205,6 +205,9 @@ BEGIN_PROVIDER [ double precision, nuclear_repulsion ] enddo enddo nuclear_repulsion *= 0.5d0 + if(point_charges)then + nuclear_repulsion += pt_chrg_nuclei_interaction + pt_chrg_interaction + endif end if call write_time(6) diff --git a/src/nuclei/point_charges.irp.f b/src/nuclei/point_charges.irp.f new file mode 100644 index 00000000..b955537f --- /dev/null +++ b/src/nuclei/point_charges.irp.f @@ -0,0 +1,209 @@ +! --- + + +BEGIN_PROVIDER [ integer, n_pts_charge ] + implicit none + BEGIN_DOC +! Number of point charges to be added to the potential + END_DOC + + logical :: has + PROVIDE ezfio_filename + if (mpi_master) then + + call ezfio_has_nuclei_n_pts_charge(has) + if (has) then + write(6,'(A)') '.. >>>>> [ IO READ: n_pts_charge ] <<<<< ..' + call ezfio_get_nuclei_n_pts_charge(n_pts_charge) + else + print *, 'nuclei/n_pts_charge not found in EZFIO file' + stop 1 + 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( n_pts_charge, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read n_pts_charge with MPI' + endif + IRP_ENDIF + + call write_time(6) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pts_charge_z, (n_pts_charge) ] + + BEGIN_DOC + ! Charge associated to each point charge. + END_DOC + + implicit none + logical :: exists + + PROVIDE ezfio_filename + + if (mpi_master) then + call ezfio_has_nuclei_pts_charge_z(exists) + 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(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_z with MPI' + endif + IRP_ENDIF + + if (exists) then + + if (mpi_master) then + write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_z ] <<<<< ..' + call ezfio_get_nuclei_pts_charge_z(pts_charge_z) + IRP_IF MPI + call MPI_BCAST(pts_charge_z, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_z with MPI' + endif + IRP_ENDIF + endif + + else + + integer :: i + do i = 1, n_pts_charge + pts_charge_z(i) = 0.d0 + enddo + + endif + print*,'Point charges ' + do i = 1, n_pts_charge + print*,'i,pts_charge_z(i)',i,pts_charge_z(i) + enddo + +END_PROVIDER + + +BEGIN_PROVIDER [ double precision, pts_charge_coord, (n_pts_charge,3) ] + + BEGIN_DOC + ! Coordinates of each point charge. + END_DOC + + implicit none + logical :: exists + + PROVIDE ezfio_filename + + if (mpi_master) then + call ezfio_has_nuclei_pts_charge_coord(exists) + 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(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_coord with MPI' + endif + IRP_ENDIF + + if (exists) then + + if (mpi_master) then + double precision, allocatable :: buffer(:,:) + allocate (buffer(n_pts_charge,3)) + write(6,'(A)') '.. >>>>> [ IO READ: pts_charge_coord ] <<<<< ..' + call ezfio_get_nuclei_pts_charge_coord(buffer) + integer :: i,j + do i=1,3 + do j=1,n_pts_charge + pts_charge_coord(j,i) = buffer(j,i) + enddo + enddo + deallocate(buffer) + IRP_IF MPI + call MPI_BCAST(pts_charge_coord, (n_pts_charge), MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr) + if (ierr /= MPI_SUCCESS) then + stop 'Unable to read pts_charge_coord with MPI' + endif + IRP_ENDIF + endif + + else + + do i = 1, n_pts_charge + pts_charge_coord(i,:) = 0.d0 + enddo + + endif + print*,'Coordinates for the point charges ' + do i = 1, n_pts_charge + write(*,'(I3,X,3(F16.8,X))') i,pts_charge_coord(i,1:3) + enddo + +END_PROVIDER + +! --- +BEGIN_PROVIDER [ double precision, pt_chrg_interaction] + implicit none + BEGIN_DOC + ! Interaction between the point charges + END_DOC + integer :: i,j + double precision :: Z_A, z_B,A_center(3), B_center(3), dist + pt_chrg_interaction = 0.d0 + do i = 1, n_pts_charge + Z_A = pts_charge_z(i) + A_center(1:3) = pts_charge_coord(i,1:3) + do j = i+1, n_pts_charge + Z_B = pts_charge_z(j) + B_center(1:3) = pts_charge_coord(j,1:3) + dist = (A_center(1)-B_center(1))**2 + (A_center(2)-B_center(2))**2 + (A_center(3)-B_center(3))**2 + dist = dsqrt(dist) + pt_chrg_interaction += Z_A*Z_B/dist + enddo + enddo + print*,'Interaction between the point charges ' + print*,'pt_chrg_interaction = ',pt_chrg_interaction +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pt_chrg_nuclei_interaction] + implicit none + BEGIN_DOC + ! repulsion between the point charges and the nuclei + END_DOC + integer :: i,j + double precision :: Z_A, z_B,A_center(3), B_center(3), dist + pt_chrg_nuclei_interaction = 0.d0 + do i = 1, n_pts_charge + Z_A = pts_charge_z(i) + A_center(1:3) = pts_charge_coord(i,1:3) + do j = 1, nucl_num + Z_B = nucl_charge(j) + B_center(1:3) = nucl_coord(j,1:3) + dist = (A_center(1)-B_center(1))**2 + (A_center(2)-B_center(2))**2 + (A_center(3)-B_center(3))**2 + dist = dsqrt(dist) + pt_chrg_nuclei_interaction += Z_A*Z_B/dist + enddo + enddo + print*,'Interaction between point charges and nuclei' + print*,'pt_chrg_nuclei_interaction = ',pt_chrg_nuclei_interaction +END_PROVIDER + From cf0c8e75ae2bdcc1a94b42a1c4a0171677dd2ce3 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 7 Feb 2023 13:11:42 +0100 Subject: [PATCH 2/3] added point charges test --- src/hartree_fock/10.hf.bats | 32 ++++++++++++++++++++++++++++++++ 1 file changed, 32 insertions(+) diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index 65117b76..a52ce075 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -18,6 +18,38 @@ function run() { } +function run_pt_charges() { + thresh=1.e-5 + cp ${QP_ROOT}/src/nuclei/write_pt_charges.py . + cat > hcn.xyz << EOF +3 +HCN molecule +C 0.0 0.0 0.0 +H 0.0 0.0 1.064 +N 0.0 0.0 -1.156 +EOF + +cat > hcn_charges.xyz << EOF +0.5 2.0 0.0 0.0 +0.5 -2.0 0.0 0.0 +EOF + +rm -rf hcn.ezfio +qp create_ezfio -b def2-svp hcn.xyz +qp run scf +mv hcn_charges.xyz hcn.ezfio_point_charges.xyz +python write_pt_charges.py hcn.ezfio +qp set nuclei point_charges True +qp run scf | tee hcn.ezfio.pt_charges.out + energy="$(ezfio get hartree_fock energy)" +good=-92.76613324421798 + eq $energy $good $thresh +} + +@test "point charges" { + run_pt_charges +} + @test "B-B" { # 3s run b2_stretched.ezfio -48.9950585434279 } From c34171898240886d820f80d5d7d2ba7aaacf2f53 Mon Sep 17 00:00:00 2001 From: eginer Date: Tue, 7 Feb 2023 13:23:00 +0100 Subject: [PATCH 3/3] fixed the 10.hf.bats --- src/hartree_fock/10.hf.bats | 1 + 1 file changed, 1 insertion(+) diff --git a/src/hartree_fock/10.hf.bats b/src/hartree_fock/10.hf.bats index a52ce075..20b59603 100644 --- a/src/hartree_fock/10.hf.bats +++ b/src/hartree_fock/10.hf.bats @@ -42,6 +42,7 @@ python write_pt_charges.py hcn.ezfio qp set nuclei point_charges True qp run scf | tee hcn.ezfio.pt_charges.out energy="$(ezfio get hartree_fock energy)" +rm -rf hcn.ezfio good=-92.76613324421798 eq $energy $good $thresh }