From e6eb789ab30653c398746bb28649896d773177d7 Mon Sep 17 00:00:00 2001
From: Anthony Scemama <scemama@irsamc.ups-tlse.fr>
Date: Tue, 18 Jun 2019 00:12:17 +0200
Subject: [PATCH] Dev lcpq (#49)

* Add energy components

* Fixed beta_rs

* Update do_single_excitation
---
 ocaml/qp_tunnel.ml                        |  6 +++
 src/determinants/create_excitations.irp.f | 60 ++++++++++++++++++-----
 src/dft_utils_one_e/ec_scan.irp.f         |  2 +-
 src/tools/print_ci_vectors.irp.f          |  1 +
 src/zmq/utils.irp.f                       |  3 +-
 5 files changed, 58 insertions(+), 14 deletions(-)

diff --git a/ocaml/qp_tunnel.ml b/ocaml/qp_tunnel.ml
index c35a2bac..dee01980 100644
--- a/ocaml/qp_tunnel.ml
+++ b/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
 
diff --git a/src/determinants/create_excitations.irp.f b/src/determinants/create_excitations.irp.f
index ddb9ae0f..cec87901 100644
--- a/src/determinants/create_excitations.irp.f
+++ b/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
diff --git a/src/dft_utils_one_e/ec_scan.irp.f b/src/dft_utils_one_e/ec_scan.irp.f
index 7a4b587b..4807b89f 100644
--- a/src/dft_utils_one_e/ec_scan.irp.f
+++ b/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
diff --git a/src/tools/print_ci_vectors.irp.f b/src/tools/print_ci_vectors.irp.f
index 9ba06d9a..97dfdc0b 100644
--- a/src/tools/print_ci_vectors.irp.f
+++ b/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)
diff --git a/src/zmq/utils.irp.f b/src/zmq/utils.irp.f
index 2a0c1d2e..70f0830b 100644
--- a/src/zmq/utils.irp.f
+++ b/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)