diff --git a/src/ao_one_e_ints/EZFIO.cfg b/src/ao_one_e_ints/EZFIO.cfg index 262301e0..8d4fff57 100644 --- a/src/ao_one_e_ints/EZFIO.cfg +++ b/src/ao_one_e_ints/EZFIO.cfg @@ -106,26 +106,3 @@ interface: ezfio,provider,ocaml default: 1.e-15 ezfio_name: threshold_ao -[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: (ao_one_e_ints.n_pts_charge) - -[pts_charge_coord] -type: double precision -doc: Coordinate of each point charge. -interface: ezfio -size: (ao_one_e_ints.n_pts_charge,3) - -[point_charges] -type: logical -doc: If |true|, point charges (see ao_one_e_ints/write_pt_charges.py) are added to the one-electron potential -interface: ezfio,provider,ocaml -default: False 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 c038458d..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 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/cas_based_on_top/two_body_dens_rout.irp.f b/src/cas_based_on_top/two_body_dens_rout.irp.f index 4a57a868..5d066831 100644 --- a/src/cas_based_on_top/two_body_dens_rout.irp.f +++ b/src/cas_based_on_top/two_body_dens_rout.irp.f @@ -132,7 +132,7 @@ end subroutine give_n2_cas(r1,r2,istate,n2_psi) implicit none BEGIN_DOC -! returns mu(r), f_psi, n2_psi for a general cas wave function +! returns n2_psi for a general cas wave function END_DOC integer, intent(in) :: istate double precision, intent(in) :: r1(3),r2(3) diff --git a/src/cipsi_tc_bi_ortho/get_d.irp.f b/src/cipsi_tc_bi_ortho/get_d.irp.f index 58b1972a..7fdc5e12 100644 --- a/src/cipsi_tc_bi_ortho/get_d.irp.f +++ b/src/cipsi_tc_bi_ortho/get_d.irp.f @@ -523,10 +523,10 @@ subroutine get_d1(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s ! call get_mo_bi_ortho_tc_two_es(hfix,pfix,p1,mo_num,hij_cache(1,1),mo_integrals_map) ! call get_mo_bi_ortho_tc_two_es(hfix,pfix,p2,mo_num,hij_cache(1,2),mo_integrals_map) do mm = 1, mo_num - hji_cache(mm,1) = mo_bi_ortho_tc_two_e(pfix,p1,mm,hfix) - hji_cache(mm,2) = mo_bi_ortho_tc_two_e(pfix,p2,mm,hfix) hij_cache(mm,1) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p1) hij_cache(mm,2) = mo_bi_ortho_tc_two_e(mm,hfix,pfix,p2) + hji_cache(mm,1) = mo_bi_ortho_tc_two_e(pfix,p1,mm,hfix) + hji_cache(mm,2) = mo_bi_ortho_tc_two_e(pfix,p2,mm,hfix) enddo putj = p1 do puti = 1, mo_num !HOT @@ -800,7 +800,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s if(bannedOrb(p1, 1)) cycle ! call get_mo_bi_ortho_tc_two_es(p1,h2,h1,mo_num,hij_cache1,mo_integrals_map) do mm =1, mo_num - hij_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1) + hji_cache1(mm) = mo_bi_ortho_tc_two_e(mm,p1,h2,h1) hji_cache1(mm) = mo_bi_ortho_tc_two_e(h2,h1,mm,p1) enddo do p2=1, mo_num @@ -811,8 +811,8 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat_p, mat_m, mask, h, p, s call apply_particles(mask, 1,p1,2,p2, det, ok, N_int) ! call i_h_j(gen, det, N_int, hij) !!! GUESS ON THE ORDER - call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hij) - call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hji) + call htilde_mu_mat_opt_bi_ortho_no_3e(det,gen,N_int, hji) + call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,N_int, hij) else ! print*,'ELSE ' phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) diff --git a/src/cipsi_tc_bi_ortho/selection.irp.f b/src/cipsi_tc_bi_ortho/selection.irp.f index 9c695ba8..659e50a8 100644 --- a/src/cipsi_tc_bi_ortho/selection.irp.f +++ b/src/cipsi_tc_bi_ortho/selection.irp.f @@ -807,11 +807,17 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d print*,'-- bad ' print*,psi_h_alpha_tmp,alpha_h_psi_tmp print*,'-- details good' + double precision :: accu_1, accu_2 + accu_1 = 0.d0 + accu_2 = 0.d0 do iii = 1, N_det call get_excitation_degree( psi_det(1,1,iii), det, degree, N_int) call htilde_mu_mat_bi_ortho_tot(psi_det(1,1,iii), det, N_int, i_h_alpha) call htilde_mu_mat_bi_ortho_tot(det, psi_det(1,1,iii), N_int, alpha_h_i) print*,iii,degree,i_h_alpha,alpha_h_i + accu_1 += i_h_alpha + accu_2 += alpha_h_i + print*,accu_1,accu_2 enddo ! if(dabs(psi_h_alpha*alpha_h_psi).gt.1.d-10)then diff --git a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f index 01cb57dd..dc8a4c07 100644 --- a/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f +++ b/src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f @@ -44,9 +44,10 @@ subroutine run_stochastic_cipsi pt2_data % overlap= 0.d0 pt2_data % variance = huge(1.e0) - if (s2_eig) then - call make_s2_eigenfunction - endif + !!!! WARNING !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION +! if (s2_eig) then +! call make_s2_eigenfunction +! endif print_pt2 = .False. call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2) ! call routine_save_right diff --git a/src/nuclei/EZFIO.cfg b/src/nuclei/EZFIO.cfg index 34c27c46..bd25e38a 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..f765e107 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_repulsion + pt_chrg_repulsion + 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..86038742 --- /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_repulsion] + implicit none + BEGIN_DOC + ! repulsion between the point charges + END_DOC + integer :: i,j + double precision :: Z_A, z_B,A_center(3), B_center(3), dist + pt_chrg_repulsion = 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_repulsion += Z_A*Z_B/dist + enddo + enddo + print*,'Repulsion of point charges ' + print*,'pt_chrg_repulsion = ',pt_chrg_repulsion +END_PROVIDER + +BEGIN_PROVIDER [ double precision, pt_chrg_nuclei_repulsion] + 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_repulsion = 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_repulsion += Z_A*Z_B/dist + enddo + enddo + print*,'Repulsion between point charges and nuclei' + print*,'pt_chrg_nuclei_repulsion = ',pt_chrg_nuclei_repulsion +END_PROVIDER + diff --git a/src/ao_one_e_ints/write_pt_charges.py b/src/nuclei/write_pt_charges.py similarity index 94% rename from src/ao_one_e_ints/write_pt_charges.py rename to src/nuclei/write_pt_charges.py index d4b6d251..d722faa8 100755 --- a/src/ao_one_e_ints/write_pt_charges.py +++ b/src/nuclei/write_pt_charges.py @@ -13,11 +13,11 @@ def zip_in_ezfio(ezfio,tmp): cmdzip="gzip -c "+tmp+" > "+tmpzip os.system(cmdzip) os.system("rm "+tmp) - cmdmv="mv "+tmpzip+" "+EZFIO+"/ao_one_e_ints/"+tmpzip + cmdmv="mv "+tmpzip+" "+EZFIO+"/nuclei/"+tmpzip os.system(cmdmv) def mv_in_ezfio(ezfio,tmp): - cmdmv="mv "+tmp+" "+EZFIO+"/ao_one_e_ints/"+tmp + cmdmv="mv "+tmp+" "+EZFIO+"/nuclei/"+tmp os.system(cmdmv)