removed some stuffs to clean the CIPSI module

src/cipsi/EZFIO.cfg

@@ -1,9 +1,3 @@
-[pert_2rdm]
-type: logical
-doc: If true, computes the one- and two-body rdms with perturbation theory
-interface: ezfio,provider,ocaml
-default: False
-
 [save_wf_after_selection]
 type: logical
 doc: If true, saves the wave function after the selection, before the diagonalization

src/cipsi/NEED

@@ -2,5 +2,4 @@ perturbation
 zmq
 mpi
 iterations
-two_body_rdm
 csf

src/cipsi/pert_rdm_providers.irp.f

@@ -1,183 +0,0 @@
-
-use bitmasks
-use omp_lib
-
-BEGIN_PROVIDER [ integer(omp_lock_kind), pert_2rdm_lock]
-  use f77_zmq
-  implicit none
-  call omp_init_lock(pert_2rdm_lock)
-END_PROVIDER
-
-BEGIN_PROVIDER [integer, n_orb_pert_rdm]
- implicit none
- n_orb_pert_rdm = n_act_orb
-END_PROVIDER 
-
-BEGIN_PROVIDER [integer, list_orb_reverse_pert_rdm, (mo_num)]
- implicit none
- list_orb_reverse_pert_rdm = list_act_reverse
-
-END_PROVIDER 
-
-BEGIN_PROVIDER [integer, list_orb_pert_rdm, (n_orb_pert_rdm)]
- implicit none
- list_orb_pert_rdm = list_act
-
-END_PROVIDER
-
-BEGIN_PROVIDER [double precision, pert_2rdm_provider, (n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm,n_orb_pert_rdm)]
- implicit none
- pert_2rdm_provider = 0.d0
-
-END_PROVIDER
-
-subroutine fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf, psi_det_connection, psi_coef_connection_reverse, n_det_connection)
-  use bitmasks
-  use selection_types
-  implicit none
-  
-  integer, intent(in)           :: n_det_connection
-  double precision, intent(in)  :: psi_coef_connection_reverse(N_states,n_det_connection)
-  integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection)
-  integer, intent(in) :: i_generator, sp, h1, h2
-  double precision, intent(in) :: mat(N_states, mo_num, mo_num)
-  logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num)
-  double precision, intent(in)           :: fock_diag_tmp(mo_num)
-  double precision, intent(in)    :: E0(N_states)
-  type(pt2_type), intent(inout)   :: pt2_data
-  type(selection_buffer), intent(inout) :: buf
-  logical :: ok
-  integer :: s1, s2, p1, p2, ib, j, istate, jstate
-  integer(bit_kind) :: mask(N_int, 2), det(N_int, 2)
-  double precision :: e_pert, delta_E, val, Hii, sum_e_pert, tmp, alpha_h_psi, coef(N_states)
-  double precision, external :: diag_H_mat_elem_fock
-  double precision :: E_shift
-
-  logical, external :: detEq
-  double precision, allocatable :: values(:)
-  integer, allocatable          :: keys(:,:)
-  integer                       :: nkeys
-  integer :: sze_buff
-  sze_buff = 5 * mo_num ** 2 
-  allocate(keys(4,sze_buff),values(sze_buff))
-  nkeys = 0
-  if(sp == 3) then
-    s1 = 1
-    s2 = 2
-  else
-    s1 = sp
-    s2 = sp
-  end if
-  call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
-  E_shift = 0.d0
-
-  if (h0_type == 'CFG') then
-    j = det_to_configuration(i_generator)
-    E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j)
-  endif
-
-  do p1=1,mo_num
-    if(bannedOrb(p1, s1)) cycle
-    ib = 1
-    if(sp /= 3) ib = p1+1
-
-    do p2=ib,mo_num
-
-! -----
-! /!\ Generating only single excited determinants doesn't work because a
-! determinant generated by a single excitation may be doubly excited wrt
-! to a determinant of the future. In that case, the determinant will be
-! detected as already generated when generating in the future with a
-! double excitation.
-!     
-!      if (.not.do_singles) then
-!        if ((h1 == p1) .or. (h2 == p2)) then
-!          cycle
-!        endif
-!      endif
-!
-!      if (.not.do_doubles) then
-!        if ((h1 /= p1).and.(h2 /= p2)) then
-!          cycle
-!        endif
-!      endif
-! -----
-
-      if(bannedOrb(p2, s2)) cycle
-      if(banned(p1,p2)) cycle
-
-
-      if( sum(abs(mat(1:N_states, p1, p2))) == 0d0) cycle
-      call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
-
-      if (do_only_cas) then
-        integer, external :: number_of_holes, number_of_particles
-        if (number_of_particles(det)>0) then
-          cycle
-        endif
-        if (number_of_holes(det)>0) then
-          cycle
-        endif
-      endif
-
-      if (do_ddci) then
-        logical, external  :: is_a_two_holes_two_particles
-        if (is_a_two_holes_two_particles(det)) then
-          cycle
-        endif
-      endif
-
-      if (do_only_1h1p) then
-        logical, external :: is_a_1h1p
-        if (.not.is_a_1h1p(det)) cycle
-      endif
-
-
-      Hii = diag_H_mat_elem_fock(psi_det_generators(1,1,i_generator),det,fock_diag_tmp,N_int)
-
-      sum_e_pert = 0d0
-      integer :: degree
-      call get_excitation_degree(det,HF_bitmask,degree,N_int)
-      if(degree == 2)cycle
-      do istate=1,N_states
-        delta_E = E0(istate) - Hii + E_shift
-        alpha_h_psi = mat(istate, p1, p2)
-        val = alpha_h_psi + alpha_h_psi
-        tmp = dsqrt(delta_E * delta_E + val * val)
-        if (delta_E < 0.d0) then
-            tmp = -tmp
-        endif
-        e_pert = 0.5d0 * (tmp - delta_E)
-        coef(istate) = e_pert / alpha_h_psi
-        print*,e_pert,coef,alpha_h_psi
-        pt2_data % pt2(istate) += e_pert
-        pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi
-      enddo
-
-      do istate=1,N_states
-        alpha_h_psi = mat(istate, p1, p2)
-        e_pert = coef(istate) * alpha_h_psi
-        do jstate=1,N_states
-          pt2_data % overlap(jstate,jstate) = coef(istate) * coef(jstate)
-        enddo
-
-        if (weight_selection /= 5) then
-          ! Energy selection
-          sum_e_pert = sum_e_pert + e_pert * selection_weight(istate)
-
-        else
-          ! Variance selection
-          sum_e_pert = sum_e_pert - alpha_h_psi * alpha_h_psi * selection_weight(istate)
-        endif
-      end do
-      call give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff)
-
-      if(sum_e_pert <= buf%mini) then
-        call add_to_selection_buffer(buf, det, sum_e_pert)
-      end if
-    end do
-  end do
-  call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
-end
-
-

src/cipsi/pt2_stoch_routines.irp.f

@@ -133,7 +133,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
   PROVIDE psi_bilinear_matrix_transp_order psi_selectors_coef_transp psi_det_sorted
   PROVIDE psi_det_hii selection_weight pseudo_sym
   PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
-  PROVIDE pert_2rdm excitation_beta_max  excitation_alpha_max excitation_max
+  PROVIDE excitation_beta_max  excitation_alpha_max excitation_max
 
   if (h0_type == 'CFG') then
     PROVIDE psi_configuration_hii det_to_configuration

src/cipsi/selection.irp.f

@@ -464,14 +464,14 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
       allocate (fullminilist (N_int, 2, fullinteresting(0)), &
                     minilist (N_int, 2,     interesting(0)) )
-      if(pert_2rdm)then
+!      if(pert_2rdm)then
-        allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
+!        allocate(coef_fullminilist_rev(N_states,fullinteresting(0)))
-        do i=1,fullinteresting(0)
+!        do i=1,fullinteresting(0)
-          do j = 1, N_states
+!          do j = 1, N_states
-            coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
+!            coef_fullminilist_rev(j,i) = psi_coef_sorted(fullinteresting(i),j)
-          enddo
+!          enddo
-        enddo
+!        enddo
-      endif
+!      endif
 
       do i=1,fullinteresting(0)
         do k=1,N_int
@@ -531,19 +531,19 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
 
             call splash_pq(mask, sp, minilist, i_generator, interesting(0), bannedOrb, banned, mat, interesting)
 
-            if(.not.pert_2rdm)then
+!            if(.not.pert_2rdm)then
              call fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
-            else
+!            else
-             call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
+!             call fill_buffer_double_rdm(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf,fullminilist, coef_fullminilist_rev, fullinteresting(0))
-            endif
+!            endif
           end if
         enddo
         if(s1 /= s2) monoBdo = .false.
       enddo
       deallocate(fullminilist,minilist)
-      if(pert_2rdm)then
+!      if(pert_2rdm)then
-        deallocate(coef_fullminilist_rev)
+!        deallocate(coef_fullminilist_rev)
-      endif
+!      endif
     enddo
   enddo
   deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)

src/cipsi/update_2rdm.irp.f

@@ -1,223 +0,0 @@
-use bitmasks
-
-subroutine give_2rdm_pert_contrib(det,coef,psi_det_connection,psi_coef_connection_reverse,n_det_connection,nkeys,keys,values,sze_buff)
- implicit none
- integer, intent(in) :: n_det_connection,sze_buff
- double precision, intent(in) :: coef(N_states)
- integer(bit_kind), intent(in) :: det(N_int,2)
- integer(bit_kind), intent(in) :: psi_det_connection(N_int,2,n_det_connection) 
- double precision, intent(in)  :: psi_coef_connection_reverse(N_states,n_det_connection)
- integer,           intent(inout) :: keys(4,sze_buff),nkeys
- double precision,  intent(inout) :: values(sze_buff)
- integer :: i,j
- integer                        :: exc(0:2,2,2)
- integer                        :: degree
- double precision               :: phase, contrib
- do i = 1, n_det_connection
-  call get_excitation(det,psi_det_connection(1,1,i),exc,degree,phase,N_int) 
-  if(degree.gt.2)cycle
-  contrib = 0.d0
-  do j = 1, N_states
-   contrib += state_average_weight(j) * psi_coef_connection_reverse(j,i) * phase * coef(j)
-  enddo
-  ! case of single excitations 
-  if(degree == 1)then
-   if (nkeys + 6 * elec_alpha_num .ge. sze_buff)then
-    call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
-    nkeys = 0
-   endif
-   call update_buffer_single_exc_rdm(det,psi_det_connection(1,1,i),exc,phase,contrib,nkeys,keys,values,sze_buff)
-  else 
- !! case of double excitations 
- ! if (nkeys + 4 .ge. sze_buff)then
- !  call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
- !  nkeys = 0
- ! endif
- ! call update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
-  endif
- enddo
-!call update_keys_values(keys,values,nkeys,n_orb_pert_rdm,pert_2rdm_provider,pert_2rdm_lock)
-!nkeys = 0
-
-end
-
-subroutine update_buffer_single_exc_rdm(det1,det2,exc,phase,contrib,nkeys,keys,values,sze_buff)
- implicit none
- integer, intent(in) :: sze_buff
- integer(bit_kind), intent(in) :: det1(N_int,2)
- integer(bit_kind), intent(in) :: det2(N_int,2)
- integer,intent(in)             :: exc(0:2,2,2)
- double precision,intent(in)    :: phase, contrib
- integer, intent(inout)         :: nkeys, keys(4,sze_buff)
- double precision, intent(inout):: values(sze_buff)
- 
- integer                        :: occ(N_int*bit_kind_size,2)
- integer                        :: n_occ_ab(2),ispin,other_spin
- integer :: h1,h2,p1,p2,i
- call bitstring_to_list_ab(det1, occ, n_occ_ab, N_int)
-  
- if (exc(0,1,1) == 1) then
-  ! Mono alpha
-  h1 = exc(1,1,1)
-  p1 = exc(1,2,1)
-  ispin = 1
-  other_spin = 2
- else
-  ! Mono beta
-  h1 = exc(1,1,2)
-  p1 = exc(1,2,2)
-  ispin = 2
-  other_spin = 1
- endif
- if(list_orb_reverse_pert_rdm(h1).lt.0)return
- h1 = list_orb_reverse_pert_rdm(h1)
- if(list_orb_reverse_pert_rdm(p1).lt.0)return
- p1 = list_orb_reverse_pert_rdm(p1)
- !update the alpha/beta part
- do i = 1, n_occ_ab(other_spin)
-  h2 = occ(i,other_spin)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
-  h2 = list_orb_reverse_pert_rdm(h2)
-
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = h2
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h2
-  keys(2,nkeys) = h1
-  keys(3,nkeys) = h2
-  keys(4,nkeys) = p1
- enddo 
- !update the same spin part 
-!do i = 1, n_occ_ab(ispin)
-! h2 = occ(i,ispin)
-! if(list_orb_reverse_pert_rdm(h2).lt.0)return
-! h2 = list_orb_reverse_pert_rdm(h2)
-
-! nkeys += 1
-! values(nkeys) = 0.5d0 * contrib * phase
-! keys(1,nkeys) = h1
-! keys(2,nkeys) = h2
-! keys(3,nkeys) = p1
-! keys(4,nkeys) = h2
-
-! nkeys += 1
-! values(nkeys) = - 0.5d0 * contrib * phase
-! keys(1,nkeys) = h1
-! keys(2,nkeys) = h2
-! keys(3,nkeys) = h2
-! keys(4,nkeys) = p1
-!
-! nkeys += 1
-! values(nkeys) = 0.5d0 * contrib * phase
-! keys(1,nkeys) = h2
-! keys(2,nkeys) = h1
-! keys(3,nkeys) = h2
-! keys(4,nkeys) = p1
-
-! nkeys += 1
-! values(nkeys) = - 0.5d0 * contrib * phase
-! keys(1,nkeys) = h2
-! keys(2,nkeys) = h1
-! keys(3,nkeys) = p1
-! keys(4,nkeys) = h2
-!enddo 
-
-end
-
-subroutine update_buffer_double_exc_rdm(exc,phase,contrib,nkeys,keys,values,sze_buff)
- implicit none
- integer, intent(in) :: sze_buff
- integer,intent(in)             :: exc(0:2,2,2)
- double precision,intent(in)    :: phase, contrib
- integer, intent(inout)         :: nkeys, keys(4,sze_buff)
- double precision, intent(inout):: values(sze_buff)
- integer :: h1,h2,p1,p2
-
- if (exc(0,1,1) == 1) then
-  ! Double alpha/beta
-  h1 = exc(1,1,1) 
-  h2 = exc(1,1,2) 
-  p1 = exc(1,2,1)
-  p2 = exc(1,2,2)
-  ! check if the orbitals involved are within the orbital range
-  if(list_orb_reverse_pert_rdm(h1).lt.0)return
-  h1 = list_orb_reverse_pert_rdm(h1)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
-  h2 = list_orb_reverse_pert_rdm(h2)
-  if(list_orb_reverse_pert_rdm(p1).lt.0)return
-  p1 = list_orb_reverse_pert_rdm(p1)
-  if(list_orb_reverse_pert_rdm(p2).lt.0)return
-  p2 = list_orb_reverse_pert_rdm(p2)
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = p2
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = p1
-  keys(2,nkeys) = p2
-  keys(3,nkeys) = h1
-  keys(4,nkeys) = h2 
-
- else 
-  if (exc(0,1,1) == 2) then
-    ! Double alpha/alpha
-   h1 =  exc(1,1,1)
-   h2 =  exc(2,1,1)
-   p1 =  exc(1,2,1)
-   p2 =  exc(2,2,1)
-  else if (exc(0,1,2) == 2) then
-   ! Double beta
-   h1 = exc(1,1,2)
-   h2 = exc(2,1,2)
-   p1 = exc(1,2,2)
-   p2 = exc(2,2,2)
-  endif
-  ! check if the orbitals involved are within the orbital range
-  if(list_orb_reverse_pert_rdm(h1).lt.0)return
-  h1 = list_orb_reverse_pert_rdm(h1)
-  if(list_orb_reverse_pert_rdm(h2).lt.0)return
-  h2 = list_orb_reverse_pert_rdm(h2)
-  if(list_orb_reverse_pert_rdm(p1).lt.0)return
-  p1 = list_orb_reverse_pert_rdm(p1)
-  if(list_orb_reverse_pert_rdm(p2).lt.0)return
-  p2 = list_orb_reverse_pert_rdm(p2)
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = p2
- 
-  nkeys += 1
-  values(nkeys) = - 0.5d0 * contrib * phase
-  keys(1,nkeys) = h1
-  keys(2,nkeys) = h2
-  keys(3,nkeys) = p2
-  keys(4,nkeys) = p1
-                                        
-  nkeys += 1
-  values(nkeys) = 0.5d0 * contrib * phase
-  keys(1,nkeys) = h2
-  keys(2,nkeys) = h1
-  keys(3,nkeys) = p2
-  keys(4,nkeys) = p1
- 
-  nkeys += 1
-  values(nkeys) = - 0.5d0 * contrib * phase
-  keys(1,nkeys) = h2
-  keys(2,nkeys) = h1
-  keys(3,nkeys) = p1
-  keys(4,nkeys) = p2
- endif
- 
-end
-
-

src/cipsi/zmq_selection.irp.f

@@ -22,7 +22,7 @@ subroutine ZMQ_selection(N_in, pt2_data)
     PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
     PROVIDE psi_bilinear_matrix_transp_order selection_weight pseudo_sym
     PROVIDE n_act_orb n_inact_orb n_core_orb n_virt_orb n_del_orb seniority_max
-    PROVIDE pert_2rdm excitation_beta_max  excitation_alpha_max excitation_max
+    PROVIDE excitation_beta_max  excitation_alpha_max excitation_max
 
     call new_parallel_job(zmq_to_qp_run_socket,zmq_socket_pull,'selection')
 

src/davidson_dressed/diagonalize_ci.irp.f

@@ -12,7 +12,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
   enddo
   do j=1,min(N_det,N_states)
     write(st,'(I4)') j
-    call write_double(6,CI_energy_dressed(j),'Energy of state '//trim(st))
+    call write_double(6,CI_energy_dressed(j),'Energy dressed of state '//trim(st))
     call write_double(6,CI_eigenvectors_s2_dressed(j),'S^2 of state '//trim(st))
   enddo
 

src/dft_one_e/mu_erf_dft.irp.f

@@ -58,6 +58,7 @@ BEGIN_PROVIDER [double precision, mu_of_r_dft_average]
   r(3) = final_grid_points(3,i)
   call dm_dft_alpha_beta_at_r(r,dm_a,dm_b)
   rho = dm_a + dm_b
+  if(mu_of_r_dft(i).gt.1.d+3)cycle
   mu_of_r_dft_average += rho * mu_of_r_dft(i) * final_weight_at_r_vector(i)
  enddo
  mu_of_r_dft_average = mu_of_r_dft_average / dble(elec_alpha_num + elec_beta_num)

src/iterations/print_summary.irp.f

@@ -98,7 +98,7 @@ subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s
     enddo
   endif
 
-  call print_energy_components()
+!  call print_energy_components()
 
 end subroutine
 

src/tools/NEED

@@ -2,3 +2,4 @@ fci
 mo_two_e_erf_ints
 aux_quantities
 hartree_fock
+two_body_rdm