same as previous but including the files

plugins/dress_zmq/run_dress_slave.irp.f

@@ -41,6 +41,8 @@ subroutine run_dress_slave(thread,iproce,energy)
 ! double precision, external :: omp_get_wtime
   double precision :: time, time0
   integer :: ntask_tbd, task_tbd(Nproc), i_gen_tbd(Nproc), subset_tbd(Nproc)
+  logical :: interesting
+
 
   allocate(delta_det(N_states, N_det, 0:pt2_N_teeth+1, 2))
   allocate(cp(N_states, N_det, dress_N_cp, 2))
@@ -70,7 +72,7 @@ subroutine run_dress_slave(thread,iproce,energy)
   ending = dress_N_cp+1
   ntask_tbd = 0
   !$OMP PARALLEL DEFAULT(SHARED) &
-  !$OMP PRIVATE(breve_delta_m, task_id) &
+  !$OMP PRIVATE(interesting, breve_delta_m, task_id) &
   !$OMP PRIVATE(tmp,fac,m,l,t,sum_f,n_tasks) &
   !$OMP PRIVATE(i,p,will_send, i_generator, subset, iproc) &
   !$OMP PRIVATE(zmq_to_qp_run_socket, zmq_socket_push, worker_id) &
@@ -157,9 +159,12 @@ subroutine run_dress_slave(thread,iproce,energy)
       !UPDATE i_generator
 
       breve_delta_m(:,:,:) = 0d0
-      call generator_start(i_generator, iproc)
+      call generator_start(i_generator, iproc, interesting)
+
       time0 = omp_get_wtime()
+      if(interesting) then
         call alpha_callback(breve_delta_m, i_generator, subset, pt2_F(i_generator), iproc)
+      end if
       time = omp_get_wtime()
       t = dress_T(i_generator)
     

plugins/mrcc/EZFIO.cfg

@@ -0,0 +1,45 @@
+[lambda_type]
+type: Positive_int
+doc: lambda type 
+interface: ezfio,provider,ocaml
+default: 0
+
+[energy]
+type: double precision
+doc: Calculated energy
+interface: ezfio
+
+[energy_pt2]
+type: double precision
+doc: Calculated energy with PT2 contribution
+interface: ezfio
+
+[perturbative_triples]
+type: logical
+doc: Compute perturbative contribution of the Triples
+interface: ezfio,provider,ocaml
+default: true
+
+[energy]
+type: double precision
+doc: Calculated energy
+interface: ezfio
+
+[thresh_dressed_ci]
+type: Threshold
+doc: Threshold on the convergence of the dressed CI energy
+interface: ezfio,provider,ocaml
+default: 1.e-5
+
+[n_it_max_dressed_ci]
+type: Strictly_positive_int
+doc: Maximum number of dressed CI iterations
+interface: ezfio,provider,ocaml
+default: 30
+
+[dress_relative_error]
+type: Normalized_float
+doc: Stop stochastic dressing when the relative error is smaller than PT2_relative_error
+interface: ezfio,provider,ocaml
+default: 0.01
+

plugins/mrcc/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+dress_zmq DavidsonDressed Selectors_full Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils ZMQ

plugins/mrcc/mrcc.irp.f

@@ -0,0 +1,16 @@
+program shifted_bk
+  implicit none
+  BEGIN_DOC
+! TODO
+  END_DOC
+  !print *, "neu verzion"  
+  
+  PROVIDE psi_bilinear_matrix_columns_loc psi_det_alpha_unique psi_det_beta_unique
+  PROVIDE psi_bilinear_matrix_rows psi_det_sorted_gen_order psi_bilinear_matrix_order
+  PROVIDE psi_bilinear_matrix_transp_rows_loc psi_bilinear_matrix_transp_columns
+  PROVIDE psi_bilinear_matrix_transp_order
+
+
+  call dress_zmq()
+end
+

plugins/mrcc/mrcc_routines.irp.f

@@ -0,0 +1,315 @@
+subroutine generator_start(i_gen, iproc, interesting)
+  implicit none
+  integer, intent(in) :: i_gen, iproc
+  logical, intent(inout) :: interesting
+  integer :: i
+  logical, external :: deteq
+  PROVIDE dij
+  interesting = .true.
+  do i=1,N_det_ref
+    if(deteq(psi_det_generators(1,1,i_gen), psi_ref(1,1,i), N_int)) then
+      interesting = .false.
+      exit
+    end if
+  end do
+end subroutine
+
+ BEGIN_PROVIDER [ double precision, hij_cache_, (N_det,Nproc) ]
+&BEGIN_PROVIDER [ double precision, sij_cache_, (N_det,Nproc) ]
+&BEGIN_PROVIDER [ double precision, dIa_hla_, (N_states,N_det,Nproc) ]
+&BEGIN_PROVIDER [ double precision, dIa_sla_, (N_states,N_det,Nproc) ]
+&BEGIN_PROVIDER [ integer, excs_ , (0:2,2,2,N_det,Nproc) ]
+&BEGIN_PROVIDER [ double precision, phases_, (N_det, Nproc) ]
+BEGIN_DOC
+ ! temporay arrays for dress_with_alpha_buffer. Avoids reallocation.
+END_DOC
+END_PROVIDER
+
+
+subroutine dress_with_alpha_buffer(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minilist, det_minilist, n_minilist, alpha, iproc)
+  use bitmasks
+  implicit none
+  BEGIN_DOC
+  !delta_ij_loc(:,:,1) : dressing column for H
+  !delta_ij_loc(:,:,2) : dressing column for S2
+  !i_gen : generator index in psi_det_generators
+  !minilist : indices of determinants connected to alpha ( in psi_det )
+  !n_minilist : size of minilist
+  !alpha : alpha determinant
+  END_DOC
+  integer, intent(in)             :: Nint, Ndet, Nstates, n_minilist, iproc, i_gen
+  integer(bit_kind), intent(in)   :: alpha(Nint,2), det_minilist(Nint, 2, n_minilist)
+  integer,intent(in)              :: minilist(n_minilist)
+  integer(bit_kind)               :: dettmp(Nint,2)
+  double precision, intent(inout) :: delta_ij_loc(Nstates,N_det,2)
+  double precision :: hij, sij
+  double precision, external :: diag_H_mat_elem_fock
+  double precision :: c_alpha(N_states)
+  double precision :: hdress, sdress
+  integer :: i, l_sd, j, k, i_I, s, ni
+  logical :: ok
+  double precision               :: phase, phase2
+  integer                        :: degree, exc(0:2,2,2)
+  if(n_minilist == 1) return
+  !chekc actutally not linked to reference
+  do i=1,n_minilist
+    if(idx_non_ref_rev(minilist(i)) == 0) then
+      return
+    end if
+  end do
+  
+  c_alpha = 0d0
+  
+  do i_I=1,N_det_ref
+    do i=1,n_minilist
+    diamondloop : do j=i+1,n_minilist
+      do s=1,2
+      do ni=1,nint
+        dettmp(ni,s) = ieor(det_minilist(ni,s,i), det_minilist(ni,s,j))
+        dettmp(ni,s) = ieor(dettmp(ni,s), psi_ref(ni,s,i_I))
+        dettmp(ni,s) = ieor(dettmp(ni,s), alpha(ni,s))
+        if(dettmp(ni,s) /= 0_8) cycle diamondloop
+      end do 
+      end do
+      !diamond found
+
+      call get_excitation(psi_ref(1,1,i_I),det_minilist(1,1,j),exc,degree,phase,Nint)
+      call get_excitation(alpha,det_minilist(1,1,i),exc,degree,phase2,Nint)
+
+      do s=1,Nstates
+        c_alpha(s) += psi_ref_coef(i_I, s) * dij(i_I, idx_non_ref_rev(minilist(i)), s) &
+           * dij(i_I, idx_non_ref_rev(minilist(j)), s) * phase * phase2
+      end do
+    end do diamondloop
+    end do
+  end do
+  
+  if(maxval(c_alpha) == 0d0 .and. minval(c_alpha) == 0d0) return
+  
+  do i=1,n_minilist
+    call i_h_j_s2(alpha,det_minilist(1,1,i),N_int,hij, sij)
+    do s=1,Nstates
+      hdress = c_alpha(s) * hij
+      sdress = c_alpha(s) * sij
+      delta_ij_loc(s, minilist(i), 1) += hdress
+      delta_ij_loc(s, minilist(i), 2) += sdress
+    end do
+  end do
+end subroutine
+
+
+
+subroutine dress_with_alpha_buffer_old(Nstates,Ndet,Nint,delta_ij_loc, i_gen, minilist, det_minilist, n_minilist, alpha, iproc)
+ use bitmasks
+ implicit none
+  BEGIN_DOC
+  !delta_ij_loc(:,:,1) : dressing column for H
+  !delta_ij_loc(:,:,2) : dressing column for S2
+  !minilist : indices of determinants connected to alpha ( in psi_det_sorted )
+  !n_minilist : size of minilist
+  !alpha : alpha determinant
+  END_DOC
+  integer(bit_kind), intent(in)   :: alpha(Nint,2), det_minilist(Nint, 2, n_minilist)
+  integer,intent(in)              :: minilist(n_minilist), n_minilist, iproc, i_gen, Nstates, Ndet, Nint
+  double precision, intent(inout) :: delta_ij_loc(Nstates,Ndet,2)
+
+
+  integer                        :: i,j,k,l,m
+  integer                        :: degree1, degree2, degree
+
+  double precision               :: hIk, hla, hIl, sla, dIk(Nstates), dka(Nstates), dIa(Nstates), hka
+  double precision               :: phase, phase2
+  integer                        :: exc(0:2,2,2)
+  integer                        :: h1,h2,p1,p2,s1,s2
+  integer(bit_kind)              :: tmp_det(Nint,2), ctrl
+  integer                        :: i_state, k_sd, l_sd, m_sd, ll_sd, i_I
+  double precision :: Delta_E_inv(Nstates)
+  double precision :: sdress, hdress
+  logical :: ok, ok2
+  integer :: canbediamond
+
+  PROVIDE mo_class dij N_int N_states elec_num n_act_orb 
+  
+  if(n_minilist == 1) return
+  
+  do i=1,n_minilist
+    if(idx_non_ref_rev(minilist(i)) == 0) return
+  end do
+
+  if (perturbative_triples) then
+    PROVIDE one_anhil fock_virt_total fock_core_inactive_total one_creat
+  endif
+
+  canbediamond = 0
+  do l_sd=1,n_minilist
+    call get_excitation(det_minilist(1,1,l_sd),alpha,exc,degree1,phase,Nint)
+    call decode_exc(exc,degree1,h1,p1,h2,p2,s1,s2)
+   
+    ok = (mo_class(h1)(1:1) == 'A' .or. mo_class(h1)(1:1) == 'I') .and. &
+         (mo_class(p1)(1:1) == 'A' .or. mo_class(p1)(1:1) == 'V') 
+    if(ok .and. degree1 == 2) then
+            ok = (mo_class(h2)(1:1) == 'A' .or. mo_class(h2)(1:1) == 'I') .and. &
+                 (mo_class(p2)(1:1) == 'A' .or. mo_class(p2)(1:1) == 'V') 
+    end if
+    
+    if(ok) then
+      canbediamond += 1
+      excs_(:,:,:,l_sd,iproc) = exc(:,:,:)
+      phases_(l_sd, iproc) = phase
+    else
+      phases_(l_sd, iproc) = 0d0
+    end if
+    call i_h_j_s2(alpha,det_minilist(1,1,l_sd),Nint,hij_cache_(l_sd,iproc), sij_cache_(l_sd,iproc))
+  enddo
+  if(canbediamond <= 1) return
+
+  do i_I=1,N_det_ref
+    call get_excitation_degree(alpha,psi_ref(1,1,i_I),degree1,Nint)
+    if (degree1 > 4) then
+      cycle
+    endif
+    
+    do i_state=1,Nstates
+      dIa(i_state) = 0.d0
+    enddo
+
+    do k_sd=1,n_minilist
+      if(phases_(k_sd,iproc) == 0d0) cycle
+      call get_excitation_degree(psi_ref(1,1,i_I),det_minilist(1,1,k_sd),degree,Nint)
+      if (degree > 2) then
+        cycle
+      endif
+      
+      phase = phases_(k_sd, iproc)
+      exc(:,:,:) = excs_(:,:,:,k_sd,iproc)
+      degree2 = exc(0,1,1) + exc(0,1,2)
+      call apply_excitation(psi_ref(1,1,i_I), exc, tmp_det, ok, Nint)
+      if((.not. ok) .and. (.not. perturbative_triples)) cycle
+
+      do i_state=1,Nstates
+        dka(i_state) = 0.d0
+      enddo
+      
+      ok2 = .false.
+      !do i_state=1,Nstates
+      !  !if(dka(i_state) == 0) cycle
+      !  dIk(i_state) = dij(i_I, idx_non_ref_rev(minilist(k_sd)), i_state)
+      !  if(dIk(i_state) /= 0d0) then
+      !    ok2 = .true.
+      !  endif
+      !enddo
+      !if(.not. ok2) cycle
+
+      if (ok) then
+        phase2 = 0d0
+        do l_sd=k_sd+1,n_minilist
+          if(phases_(l_sd, iproc) == 0d0) cycle
+          call get_excitation_degree(tmp_det,det_minilist(1,1,l_sd),degree,Nint)
+          if (degree == 0) then
+            do i_state=1,Nstates
+              dIk(i_state) = dij(i_I, idx_non_ref_rev(minilist(k_sd)), i_state)
+              if(dIk(i_state) /= 0d0) then
+                if(phase2 == 0d0) call get_excitation(psi_ref(1,1,i_I),det_minilist(1,1,l_sd),exc,degree,phase2,Nint)
+                dka(i_state) = dij(i_I, idx_non_ref_rev(minilist(l_sd)), i_state) * phase * phase2
+              end if
+            end do
+
+            exit
+
+          endif
+        enddo
+      else if (perturbative_triples) then
+        hka = hij_cache_(k_sd,iproc)
+        if (dabs(hka) > 1.d-12) then
+          call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv)
+
+          do i_state=1,Nstates
+            ASSERT (Delta_E_inv(i_state) < 0.d0)
+            dka(i_state) = hka / Delta_E_inv(i_state)
+          enddo
+        endif
+      endif
+      
+           
+      if (perturbative_triples.and. (degree2 == 1) ) then
+          if(sum(popcnt(tmp_det(:,1))) /= elec_alpha_num) stop "STOP 1"
+  if(sum(popcnt(tmp_det(:,2))) /= elec_beta_num) stop "STOP 2"
+  if(sum(popcnt(tmp_det(:,1))) /= elec_alpha_num) stop "STOP 3"
+  if(sum(popcnt(tmp_det(:,2))) /= elec_beta_num) stop "STOP 4"
+  
+
+          call i_h_j(psi_ref(1,1,i_I),tmp_det,Nint,hka)
+          hka = hij_cache_(k_sd,iproc) - hka
+          if (dabs(hka) > 1.d-12) then
+            call get_delta_e_dyall_general_mp(psi_ref(1,1,i_I),alpha,Delta_E_inv)
+            do i_state=1,Nstates
+              ASSERT (Delta_E_inv(i_state) < 0.d0)
+              dka(i_state) = hka / Delta_E_inv(i_state)
+            enddo
+          endif
+      endif
+      do i_state=1,Nstates
+        dIa(i_state) = dIa(i_state) + dIk(i_state) * dka(i_state)
+      enddo
+    enddo
+    
+    ok2 = .false.
+    do i_state=1,Nstates
+      if(dIa(i_state) /= 0d0) ok2 = .true.
+    enddo
+    if(.not. ok2) cycle
+
+    do l_sd=1,n_minilist
+      k_sd = minilist(l_sd)
+      hla = hij_cache_(l_sd,iproc)
+      sla = sij_cache_(l_sd,iproc)
+      do i_state=1,Nstates
+        hdress =  dIa(i_state) * hla * psi_ref_coef(i_I,i_state)
+        sdress =  dIa(i_state) * sla * psi_ref_coef(i_I,i_state)
+        !!!$OMP ATOMIC
+        delta_ij_loc(i_state,k_sd,1) += hdress
+        !!!$OMP ATOMIC
+        delta_ij_loc(i_state,k_sd,2) += sdress
+      enddo
+    enddo
+  enddo
+end subroutine
+
+
+
+
+
+!! TESTS MINILIST
+subroutine test_minilist(minilist, n_minilist, alpha)
+  use bitmasks
+  implicit none
+  integer, intent(in)            :: n_minilist
+  integer(bit_kind),intent(in)  :: alpha(N_int, 2)
+  integer, intent(in) :: minilist(n_minilist)
+  integer :: a, i, deg
+  integer :: refc(N_det), testc(N_det)
+  
+  refc = 0
+  testc = 0
+  do i=1,N_det
+    call get_excitation_degree(psi_det(1,1,i), alpha, deg, N_int) 
+    if(deg <= 2) refc(i) = refc(i) + 1
+  end do
+  do i=1,n_minilist
+    call get_excitation_degree(psi_det(1,1,minilist(i)), alpha, deg, N_int) 
+    if(deg <= 2) then
+      testc(minilist(i)) += 1
+    else
+      stop "NON LINKED IN MINILIST"
+    end if
+  end do
+  
+  do i=1,N_det
+    if(refc(i) /= testc(i)) then
+      print *, "MINILIST FAIL ", sum(refc), sum(testc), n_minilist
+      exit
+    end if
+  end do
+end subroutine
+
+

plugins/mrcc/mrcc_slave.irp.f

@@ -0,0 +1,170 @@
+program shifted_bk_slave
+  implicit none
+  BEGIN_DOC
+! Helper program to compute the dress in distributed mode.
+  END_DOC
+
+  read_wf = .False.
+  distributed_davidson = .False.
+  SOFT_TOUCH read_wf distributed_davidson
+  call provide_all
+  call switch_qp_run_to_master
+  call run_w
+end
+
+subroutine provide_all
+  PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context N_states_diag
+  PROVIDE dress_e0_denominator mo_tot_num N_int ci_energy mpi_master zmq_state zmq_context
+  PROVIDE psi_det psi_coef threshold_generators threshold_selectors state_average_weight 
+  PROVIDE N_det_selectors dress_stoch_istate N_det 
+end
+
+subroutine run_w
+  use f77_zmq
+  
+  implicit none
+  IRP_IF MPI
+    include 'mpif.h'
+  IRP_ENDIF
+
+  integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
+  integer(ZMQ_PTR) :: zmq_to_qp_run_socket
+  double precision :: energy(N_states)
+  character*(64) :: states(3)
+  character*(64) :: old_state
+  integer :: rc, i, ierr
+  double precision :: t0, t1
+  
+  integer, external              :: zmq_get_dvector, zmq_get_N_det_generators 
+  integer, external              :: zmq_get_ivector
+  integer, external              :: zmq_get_psi, zmq_get_N_det_selectors, zmq_get_int
+  integer, external              :: zmq_get_N_states_diag
+
+  zmq_context = f77_zmq_ctx_new ()
+  states(1) = 'selection'
+  states(2) = 'davidson'
+  states(3) = 'dress'
+  old_state = 'Waiting'
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+
+  PROVIDE psi_det psi_coef threshold_generators threshold_selectors state_average_weight mpi_master
+  PROVIDE zmq_state N_det_selectors dress_stoch_istate N_det dress_e0_denominator
+  PROVIDE N_det_generators N_states N_states_diag
+  IRP_IF MPI
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+  do
+
+    if (mpi_master) then
+      call wait_for_states(states,zmq_state,size(states))
+      if (zmq_state(1:64) == old_state(1:64)) then
+        call sleep(1)
+        cycle
+      else
+        old_state(1:64) = zmq_state(1:64)
+      endif
+      print *,  trim(zmq_state)
+    endif
+
+    IRP_IF MPI_DEBUG
+      print *,  irp_here, mpi_rank
+      call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+    IRP_ENDIF
+    IRP_IF MPI
+      call MPI_BCAST (zmq_state, 128, MPI_CHARACTER, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        print *,  irp_here, 'error in broadcast of zmq_state'
+      endif
+    IRP_ENDIF
+
+    if(zmq_state(1:7) == 'Stopped') then
+      exit
+    endif
+
+
+    if (zmq_state(1:8) == 'davidson') then
+
+      ! Davidson
+      ! --------
+
+      call wall_time(t0)
+      if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
+      if (zmq_get_N_states_diag(zmq_to_qp_run_socket,1) == -1) cycle
+      if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states_diag) == -1) cycle
+
+      call wall_time(t1)
+      if (mpi_master) then
+        call write_double(6,(t1-t0),'Broadcast time')
+      endif
+
+      call omp_set_nested(.True.)
+      call davidson_slave_tcp(0)
+      call omp_set_nested(.False.)
+      print *,  'Davidson done'
+      IRP_IF MPI
+        call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+        if (ierr /= MPI_SUCCESS) then
+          print *,  irp_here, 'error in barrier'
+        endif
+      IRP_ENDIF
+      print *,  'All Davidson done'
+
+    else if (zmq_state(1:5) == 'dress') then
+
+      ! Dress
+      ! ---
+
+      call wall_time(t0)
+      if (zmq_get_psi(zmq_to_qp_run_socket,1) == -1) cycle
+
+      if (zmq_get_N_det_generators (zmq_to_qp_run_socket, 1) == -1) cycle
+
+      if (zmq_get_N_det_selectors(zmq_to_qp_run_socket, 1) == -1) cycle
+
+      if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_generators',threshold_generators,1) == -1) cycle
+
+      if (zmq_get_dvector(zmq_to_qp_run_socket,1,'threshold_selectors',threshold_selectors,1) == -1) cycle
+
+      if (zmq_get_dvector(zmq_to_qp_run_socket,1,'energy',energy,N_states) == -1) cycle
+
+      if (zmq_get_int(zmq_to_qp_run_socket,1,'dress_stoch_istate',dress_stoch_istate) == -1) cycle
+
+      if (zmq_get_dvector(zmq_to_qp_run_socket,1,'state_average_weight',state_average_weight,N_states) == -1) cycle
+
+      psi_energy(1:N_states) = energy(1:N_states)
+      TOUCH psi_energy state_average_weight dress_stoch_istate threshold_selectors threshold_generators
+      if (mpi_master) then
+        print *,  'N_det', N_det
+        print *,  'N_det_generators', N_det_generators
+        print *,  'N_det_selectors', N_det_selectors
+        print *,  'psi_energy', psi_energy
+        print *,  'dress_stoch_istate', dress_stoch_istate
+        print *,  'state_average_weight', state_average_weight
+      endif
+
+      call wall_time(t1)
+      call write_double(6,(t1-t0),'Broadcast time')
+
+      call dress_slave_tcp(0, energy)
+
+
+      IRP_IF MPI
+        call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+        if (ierr /= MPI_SUCCESS) then
+          print *,  irp_here, 'error in barrier'
+        endif
+      IRP_ENDIF
+      print *,  'All dress done'
+
+    endif
+
+  end do
+  IRP_IF MPI
+    call MPI_finalize(ierr)
+  IRP_ENDIF
+end
+
+
+
+

plugins/shiftedbk/shifted_bk_routines.irp.f

@@ -10,11 +10,12 @@
 END_PROVIDER
 
 
-subroutine generator_start(i_gen, iproc)
+subroutine generator_start(i_gen, iproc, interesting)
   implicit none
   integer, intent(in) :: i_gen, iproc
+  logical, intent(inout) :: interesting
   integer :: i
-  
+  interesting = .true.
   call build_fock_tmp(fock_diag_tmp_(1,1,iproc),psi_det_generators(1,1,i_gen),N_int)
 end subroutine