Dev lcpq (#49)

* Add energy components

* Fixed beta_rs

* Update do_single_excitation

ocaml/qp_tunnel.ml

@@ -363,6 +363,12 @@ let () =
                |> Zmq.Socket.send socket_in
       in
 
+      Printf.printf "On remote hosts, create ssh tunnel using:
+ssh -L %d:%s:%d -L %d:%s:%d -L %d:%s:%d %s\n%!" 
+                (port  ) localhost (localport  )
+                (port+1) localhost (localport+1)
+                (port+9) localhost (localport+9)
+                (Unix.gethostname ());
       Printf.printf "Ready\n%!";
       while !run_status do
 

src/determinants/create_excitations.irp.f

@@ -12,6 +12,7 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok)
   integer(bit_kind), intent(inout) :: key_in(N_int,2)
   integer, intent(out)           :: i_ok
   integer                        :: k,j,i
+  integer(bit_kind)              :: mask
   use bitmasks
   ASSERT (i_hole > 0 )
   ASSERT (i_particle <= mo_num)
@@ -19,31 +20,66 @@ subroutine do_single_excitation(key_in,i_hole,i_particle,ispin,i_ok)
   ! hole
   k = shiftr(i_hole-1,bit_kind_shift)+1
   j = i_hole-shiftl(k-1,bit_kind_shift)-1
+  mask = ibset(0_bit_kind,j)
 ! check whether position j is occupied
-  if (ibits(key_in(k,ispin),j,1).eq.1) then
+  if (iand(key_in(k,ispin),mask) /= 0_bit_kind) then
    key_in(k,ispin) = ibclr(key_in(k,ispin),j)
   else 
    i_ok= -1
+   return
   end if
 
   ! particle
   k = shiftr(i_particle-1,bit_kind_shift)+1
   j = i_particle-shiftl(k-1,bit_kind_shift)-1
-  key_in(k,ispin) = ibset(key_in(k,ispin),j)
+  mask = ibset(0_bit_kind,j)
+  if (iand(key_in(k,ispin),mask) == 0_bit_kind) then
+   key_in(k,ispin) = ibset(key_in(k,ispin),j)
+  else 
+   i_ok= -1
+   return
+  end if
 
-  integer                        :: n_elec_tmp
-  n_elec_tmp = 0
-  do i = 1, N_int
-    n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
-  enddo
-  if(n_elec_tmp .ne. elec_num)then
-    !print*, n_elec_tmp,elec_num
-    !call debug_det(key_in,N_int)
-    i_ok = -1
-  endif
+!  integer                        :: n_elec_tmp
+!  n_elec_tmp = 0
+!  do i = 1, N_int
+!    n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2))
+!  enddo
+!  if(n_elec_tmp .ne. elec_num)then
+!    print*, n_elec_tmp,elec_num
+!    call debug_det(key_in,N_int)
+!    stop -1
+!  endif
 end
 
 
+subroutine build_singly_excited_wavefunction(i_hole,i_particle,ispin,det_out,coef_out)
+  implicit none
+  BEGIN_DOC
+  ! Applies the single excitation operator : a^{dager}_(i_particle) a_(i_hole) of
+  ! spin = ispin to the current wave function (psi_det, psi_coef)
+  END_DOC
+  integer, intent(in)            :: i_hole,i_particle,ispin
+  integer(bit_kind), intent(out) :: det_out(N_int,2,N_det)
+  double precision, intent(out)  :: coef_out(N_det,N_states)
+
+  integer :: k
+  integer :: i_ok
+  double precision :: phase
+  do k=1,N_det
+    coef_out(k,:) = psi_coef(k,:)
+    det_out(:,:,k) = psi_det(:,:,k)
+    call do_single_excitation(det_out(1,1,k),i_hole,i_particle,ispin,i_ok)
+    if (i_ok == 1) then
+      call get_phase(psi_det(1,1,k), det_out(1,1,k),phase,N_int)
+      coef_out(k,:) = phase * coef_out(k,:)
+    else
+      coef_out(k,:) = 0.d0
+      det_out(:,:,k) = psi_det(:,:,k)
+    endif
+  enddo
+end
+
 logical function is_spin_flip_possible(key_in,i_flip,ispin)
   implicit none
   BEGIN_DOC

src/dft_utils_one_e/ec_scan.irp.f

@@ -95,6 +95,6 @@ end
 double precision function beta_rs(rs)
  implicit none
  double precision, intent(in) ::rs
- beta_rs(rs) = 0.066725d0 * (1.d0 + 0.1d0 * rs)/(1.d0 + 0.1778d0 * rs)
+ beta_rs = 0.066725d0 * (1.d0 + 0.1d0 * rs)/(1.d0 + 0.1778d0 * rs)
 
 end

src/tools/print_ci_vectors.irp.f

@@ -24,6 +24,7 @@ subroutine routine
  implicit none
  integer :: i,k
  integer :: degree
+ call print_energy_components
  do i = 1, N_det
   print *, 'Determinant ', i
   call debug_det(psi_det(1,1,i),N_int)

src/zmq/utils.irp.f

@@ -748,10 +748,11 @@ integer function add_task_to_taskserver(zmq_to_qp_run_socket,task)
   character*(*), intent(in)      :: task
 
   integer                        :: rc, sze
-  character(len=:), allocatable :: message
+  character(len=:), allocatable  :: message
 
   add_task_to_taskserver = 0
 
+  allocate(character(len=len(task)+10+len(zmq_state)) :: message)
   message='add_task '//trim(zmq_state)//' '//trim(task)
   sze = len(message)
   rc = f77_zmq_send(zmq_to_qp_run_socket, message, sze, 0)