moved the point charges definition in nuclei and added the interaction between the charges and nuclei

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

src/ao_one_e_ints/point_charges.irp.f

@@ -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

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

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)

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(det,gen,N_int, hji)
-          call htilde_mu_mat_opt_bi_ortho_no_3e(gen,det,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)

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

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
+  !!!! WARNING  !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION
-    call make_s2_eigenfunction
+!  if (s2_eig) then
-  endif
+!    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

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

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)

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 
+

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)