From 1240b324ce731596f972e901e0bb617c849eac22 Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 27 Jan 2023 18:01:57 +0100 Subject: [PATCH 1/3] corrected big bug introduced before in determinants/fock_diag.irp.f --- src/determinants/fock_diag.irp.f | 34 ++++++++++++++++---------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/determinants/fock_diag.irp.f b/src/determinants/fock_diag.irp.f index c7c951b3..a8ce33b8 100644 --- a/src/determinants/fock_diag.irp.f +++ b/src/determinants/fock_diag.irp.f @@ -33,59 +33,59 @@ subroutine build_fock_tmp(fock_diag_tmp,det_ref,Nint) ! Occupied MOs do ii=1,elec_alpha_num i = occ(ii,1) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_one_e(i,i) - E0 = E0 + mo_bi_ortho_tc_one_e(i,i) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i) + E0 = E0 + mo_one_e_integrals(i,i) do jj=1,elec_alpha_num j = occ(jj,1) if (i==j) cycle - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) - E0 = E0 + 0.5d0*mo_bi_ortho_tc_two_e_jj_anti(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j) enddo do jj=1,elec_beta_num j = occ(jj,2) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj(i,j) - E0 = E0 + mo_bi_ortho_tc_two_e_jj(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j) + E0 = E0 + mo_two_e_integrals_jj(i,j) enddo enddo do ii=1,elec_beta_num i = occ(ii,2) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_one_e(i,i) - E0 = E0 + mo_bi_ortho_tc_one_e(i,i) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i) + E0 = E0 + mo_one_e_integrals(i,i) do jj=1,elec_beta_num j = occ(jj,2) if (i==j) cycle - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) - E0 = E0 + 0.5d0*mo_bi_ortho_tc_two_e_jj_anti(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j) enddo do jj=1,elec_alpha_num j = occ(jj,1) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j) enddo enddo ! Virtual MOs do i=1,mo_num if (fock_diag_tmp(1,i) /= 0.d0) cycle - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_one_e(i,i) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i) do jj=1,elec_alpha_num j = occ(jj,1) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j) enddo do jj=1,elec_beta_num j = occ(jj,2) - fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_bi_ortho_tc_two_e_jj(i,j) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j) enddo enddo do i=1,mo_num if (fock_diag_tmp(2,i) /= 0.d0) cycle - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_one_e(i,i) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i) do jj=1,elec_beta_num j = occ(jj,2) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj_anti(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j) enddo do jj=1,elec_alpha_num j = occ(jj,1) - fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_bi_ortho_tc_two_e_jj(i,j) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j) enddo enddo From 18d186228d9cc0afc1907fcbdd11774bdb41551e Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 27 Jan 2023 18:04:22 +0100 Subject: [PATCH 2/3] corrected cipsi_tc_bi_ortho/selection.irp.f with build_fock_tmp_tc --- src/cipsi_tc_bi_ortho/fock_diag.irp.f | 95 +++++++++++++++++++++++++++ src/cipsi_tc_bi_ortho/selection.irp.f | 2 +- 2 files changed, 96 insertions(+), 1 deletion(-) create mode 100644 src/cipsi_tc_bi_ortho/fock_diag.irp.f diff --git a/src/cipsi_tc_bi_ortho/fock_diag.irp.f b/src/cipsi_tc_bi_ortho/fock_diag.irp.f new file mode 100644 index 00000000..af6849ab --- /dev/null +++ b/src/cipsi_tc_bi_ortho/fock_diag.irp.f @@ -0,0 +1,95 @@ +subroutine build_fock_tmp_tc(fock_diag_tmp,det_ref,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Build the diagonal of the Fock matrix corresponding to a generator +! determinant. $F_{00}$ is $\langle i|H|i \rangle = E_0$. + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det_ref(Nint,2) + double precision, intent(out) :: fock_diag_tmp(2,mo_num+1) + + integer :: occ(Nint*bit_kind_size,2) + integer :: ne(2), i, j, ii, jj + double precision :: E0 + + ! Compute Fock matrix diagonal elements + call bitstring_to_list_ab(det_ref,occ,Ne,Nint) + + fock_diag_tmp = 0.d0 + E0 = 0.d0 + + if (Ne(1) /= elec_alpha_num) then + print *, 'Error in build_fock_tmp_tc (alpha)', Ne(1), Ne(2) + call debug_det(det_ref,N_int) + stop -1 + endif + if (Ne(2) /= elec_beta_num) then + print *, 'Error in build_fock_tmp_tc (beta)', Ne(1), Ne(2) + call debug_det(det_ref,N_int) + stop -1 + endif + + ! Occupied MOs + do ii=1,elec_alpha_num + i = occ(ii,1) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i) + E0 = E0 + mo_one_e_integrals(i,i) + do jj=1,elec_alpha_num + j = occ(jj,1) + if (i==j) cycle + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j) + enddo + do jj=1,elec_beta_num + j = occ(jj,2) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j) + E0 = E0 + mo_two_e_integrals_jj(i,j) + enddo + enddo + do ii=1,elec_beta_num + i = occ(ii,2) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i) + E0 = E0 + mo_one_e_integrals(i,i) + do jj=1,elec_beta_num + j = occ(jj,2) + if (i==j) cycle + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j) + E0 = E0 + 0.5d0*mo_two_e_integrals_jj_anti(i,j) + enddo + do jj=1,elec_alpha_num + j = occ(jj,1) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j) + enddo + enddo + + ! Virtual MOs + do i=1,mo_num + if (fock_diag_tmp(1,i) /= 0.d0) cycle + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_one_e_integrals(i,i) + do jj=1,elec_alpha_num + j = occ(jj,1) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj_anti(i,j) + enddo + do jj=1,elec_beta_num + j = occ(jj,2) + fock_diag_tmp(1,i) = fock_diag_tmp(1,i) + mo_two_e_integrals_jj(i,j) + enddo + enddo + do i=1,mo_num + if (fock_diag_tmp(2,i) /= 0.d0) cycle + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_one_e_integrals(i,i) + do jj=1,elec_beta_num + j = occ(jj,2) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj_anti(i,j) + enddo + do jj=1,elec_alpha_num + j = occ(jj,1) + fock_diag_tmp(2,i) = fock_diag_tmp(2,i) + mo_two_e_integrals_jj(i,j) + enddo + enddo + + fock_diag_tmp(1,mo_num+1) = E0 + fock_diag_tmp(2,mo_num+1) = E0 + +end diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 659e50a8..8137b922 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -19,7 +19,7 @@ subroutine select_connected(i_generator, E0, pt2_data, b, subset, csubset) allocate(fock_diag_tmp(2,mo_num+1)) - call build_fock_tmp(fock_diag_tmp, psi_det_generators(1,1,i_generator), N_int) + call build_fock_tmp_tc(fock_diag_tmp, psi_det_generators(1,1,i_generator), N_int) do k = 1, N_int hole_mask(k,1) = iand(generators_bitmask(k,1,s_hole), psi_det_generators(k,1,i_generator)) From 55750974cdec5e4cf7ec9f94de64c157255771cb Mon Sep 17 00:00:00 2001 From: eginer Date: Fri, 27 Jan 2023 18:15:32 +0100 Subject: [PATCH 3/3] rm src/ao_one_e_ints/point_charges.irp.f --- external/qp2-dependencies | 2 +- src/ao_one_e_ints/point_charges.irp.f | 272 -------------------------- 2 files changed, 1 insertion(+), 273 deletions(-) delete mode 100644 src/ao_one_e_ints/point_charges.irp.f diff --git a/external/qp2-dependencies b/external/qp2-dependencies index 242151e0..90ee61f5 160000 --- a/external/qp2-dependencies +++ b/external/qp2-dependencies @@ -1 +1 @@ -Subproject commit 242151e03d1d6bf042387226431d82d35845686a +Subproject commit 90ee61f5041c7c94a0c605625a264860292813a0 diff --git a/src/ao_one_e_ints/point_charges.irp.f b/src/ao_one_e_ints/point_charges.irp.f deleted file mode 100644 index 82388c0d..00000000 --- a/src/ao_one_e_ints/point_charges.irp.f +++ /dev/null @@ -1,272 +0,0 @@ - -! --- - - -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_ao_one_e_ints_n_pts_charge(has) - if (has) then - write(6,'(A)') '.. >>>>> [ IO READ: n_pts_charge ] <<<<< ..' - call ezfio_get_ao_one_e_ints_n_pts_charge(n_pts_charge) - else - print *, 'ao_one_e_ints/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_ao_one_e_ints_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_ao_one_e_ints_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_ao_one_e_ints_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_ao_one_e_ints_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, ao_integrals_pt_chrg, (ao_num,ao_num)] - - BEGIN_DOC - ! Point charge-electron interaction, in the |AO| basis set. - ! - ! :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle` - ! - ! These integrals also contain the pseudopotential integrals. - END_DOC - - implicit none - integer :: num_A, num_B, power_A(3), power_B(3) - integer :: i, j, k, l, n_pt_in, m - double precision :: alpha, beta - double precision :: A_center(3),B_center(3),C_center(3) - double precision :: overlap_x,overlap_y,overlap_z,overlap,dx,NAI_pol_mult - - ao_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