Added canonical orthogonalization and accelerated Tasks

ocaml/Message.ml

@@ -280,9 +280,8 @@ type t =
 let of_string s = 
   let l =
     String.split ~on:' ' s
-    |> List.map ~f:String.strip
-    |> List.map ~f:String.lowercase
     |> List.filter ~f:(fun x -> (String.strip x) <> "")
+    |> List.map ~f:String.lowercase
   in
   match l with
   | "add_task"   :: state :: task ->

ocaml/Qpackage.ml

@@ -55,7 +55,6 @@ let executables = lazy (
     In_channel.input_lines in_channel
      |> List.map ~f:(fun x ->
          let e = String.split ~on:' ' x
-           |> List.map ~f:String.strip
            |> List.filter ~f:(fun x -> x <> "")
          in
          match e with

ocaml/Queuing_system.ml

@@ -29,12 +29,14 @@ let add_task ~task q =
     q.next_task_id 
   in
   { q with
-    queued = q.queued @ [ task_id ] ;
+    queued = task_id :: q.queued ;
     tasks  = Map.add q.tasks ~key:task_id ~data:task ;
     next_task_id = Id.Task.increment task_id ;
   }, task_id
 
 
+
+
 let add_client q =
   let client_id =
     q.next_client_id 

ocaml/TaskServer.ml

@@ -80,7 +80,7 @@ let stop ~port =
 
   Message.Terminate (Message.Terminate_msg.create ())
   |> Message.to_string
-  |> ZMQ.Socket.send req_socket ;
+  |> ZMQ.Socket.send ~block:false req_socket ;
 
   let msg = 
     ZMQ.Socket.recv req_socket
@@ -158,13 +158,13 @@ let run ~port =
     let terminate () = 
       running := false;
       Message.to_string ok
-      |> ZMQ.Socket.send rep_socket 
+      |> ZMQ.Socket.send ~block:false rep_socket 
 
     and newjob x =
       q := Queuing_system.create ();
       job := Some x;
       Message.to_string ok
-      |> ZMQ.Socket.send rep_socket 
+      |> ZMQ.Socket.send ~block:false rep_socket 
 
     and connect state msg = 
       let push_address = 
@@ -180,7 +180,7 @@ let run ~port =
       Message.ConnectReply (Message.ConnectReply_msg.create
         ~state ~client_id ~push_address)
       |> Message.to_string
-      |> ZMQ.Socket.send rep_socket
+      |> ZMQ.Socket.send ~block:false rep_socket
 
     and disconnect state msg = 
       let s, c =
@@ -199,7 +199,7 @@ let run ~port =
       Message.DisconnectReply (Message.DisconnectReply_msg.create
         ~state ~finished)
       |> Message.to_string 
-      |> ZMQ.Socket.send rep_socket
+      |> ZMQ.Socket.send ~block:false rep_socket
 
     and add_task state msg =
       let s, task =
@@ -207,12 +207,53 @@ let run ~port =
         msg.Message.AddTask_msg.task
       in
       assert (s = state);
-      let new_q, task_id = 
-        Queuing_system.add_task ~task !q
-      in
-      q := new_q;
       Message.to_string ok
-      |> ZMQ.Socket.send rep_socket
+      |> ZMQ.Socket.send ~block:false rep_socket
+      ;
+      begin
+        match 
+           String.split ~on:' ' msg.Message.AddTask_msg.task
+           |> List.filter ~f:(fun x -> x <> "")
+        with
+        | "triangle" :: str_l :: [] ->
+          begin
+            let l =
+              Int.of_string str_l
+            in
+            for j=1 to l
+            do
+              let task = 
+                Printf.sprintf "%d %s" j str_l
+              in
+              let new_q, _ = 
+                Queuing_system.add_task ~task !q
+              in
+              q := new_q
+            done
+          end
+        | "range" :: str_i :: str_j :: [] ->
+          begin
+            let i, j =
+              Int.of_string str_i,
+              Int.of_string str_j
+            in
+            for k=i to (j+1)
+            do
+              let task = 
+                Int.to_string k
+              in
+              let new_q, task_id = 
+                Queuing_system.add_task ~task !q
+              in
+              q := new_q
+            done
+          end
+        | _ ->
+            let new_q, task_id = 
+              Queuing_system.add_task ~task !q
+            in
+            q := new_q
+      end
       
     and get_task state msg =
       let s, client_id =
@@ -231,7 +272,7 @@ let run ~port =
         | _ -> Message.Terminate (Message.Terminate_msg.create ())
       in
       Message.to_string reply
-      |> ZMQ.Socket.send rep_socket
+      |> ZMQ.Socket.send ~block:false rep_socket
       
     and task_done state msg =
       let s, client_id, task_id =
@@ -245,12 +286,12 @@ let run ~port =
       in
       q := new_q;
       Message.to_string ok
-      |> ZMQ.Socket.send rep_socket
+      |> ZMQ.Socket.send ~block:false rep_socket
       
     and error msg = 
       Message.Error (Message.Error_msg.create msg)
       |> Message.to_string
-      |> ZMQ.Socket.send rep_socket
+      |> ZMQ.Socket.send ~block:false rep_socket
     in
 
     if (polling.(0) = Some ZMQ.Poll.In) then

plugins/Hartree_Fock/huckel.irp.f

@@ -8,8 +8,6 @@ subroutine huckel_guess
   double precision               :: c
   character*(64)                 :: label
 
-  mo_coef = ao_ortho_lowdin_coef
-  TOUCH mo_coef
   label = "Guess"
   call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
                                      size(mo_mono_elec_integral,1),  &

plugins/MP2/H_apply.irp.f

@@ -8,7 +8,7 @@ s.set_perturbation("Moller_plesset")
 print s
 
 s = H_apply("mp2_selection")
-s.set_selection_pt2("Moller_plesset")
+s.set_selection_pt2("Moller_Plesset")
 print s
 END_SHELL
 

plugins/Perturbation/pt2_equations.irp.f

@@ -123,8 +123,8 @@ subroutine pt2_moller_plesset ($arguments)
   call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
   if (degree == 2) then
     call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
-    delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
+    delta_e = (Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2)) - &
-            (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
+              (Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
     delta_e = 1.d0/delta_e
   else if (degree == 1) then
     call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
@@ -134,8 +134,14 @@ subroutine pt2_moller_plesset ($arguments)
     delta_e = 0.d0
   endif
 
-  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
+  if (delta_e /= 0.d0) then
-  h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  !  call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
+    call i_H_psi_minilist(det_pert,minilist,idx_minilist,N_minilist,psi_selectors_coef,Nint,N_minilist,psi_selectors_size,N_st,i_H_psi_array)
+    h = diag_H_mat_elem_fock(det_ref,det_pert,fock_diag_tmp,Nint)
+  else
+    i_H_psi_array(:) = 0.d0
+    h = 0.d0
+  endif
   do i =1,n_st
     H_pert_diag(i) = h
     c_pert(i) = i_H_psi_array(i) *delta_e

src/Determinants/H_apply.template.f

@@ -261,6 +261,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
 
       ! Build array of the non-zero integrals of second excitation
       $filter_integrals
+
       if (ispin == 1) then
         integer                        :: jjj
         
@@ -269,7 +270,7 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
           i_b = occ_hole_tmp(kk,other_spin)
           ASSERT (i_b > 0)
           ASSERT (i_b <= mo_tot_num)
-          do jjj=1,N_elec_in_key_part_2(other_spin)     ! particule
+          do jjj=1,N_elec_in_key_part_2(other_spin)     ! particle
             j_b = occ_particle_tmp(jjj,other_spin)
             ASSERT (j_b > 0)
             ASSERT (j_b <= mo_tot_num)

src/Integrals_Bielec/ao_bi_integrals.irp.f

@@ -371,19 +371,16 @@ BEGIN_PROVIDER [ logical, ao_bielec_integrals_in_map ]
   integer(ZMQ_PTR) :: zmq_to_qp_run_socket
   call new_parallel_job(zmq_to_qp_run_socket,'ao_integrals')
 
+  
   character*(32) :: task
   do l=1,ao_num
-    do j = 1, l
+    write(task,*) 'triangle', l
-      if (ao_overlap_abs(j,l) < ao_integrals_threshold) then
+    call add_task_to_taskserver(zmq_to_qp_run_socket,task)
-        cycle
-      endif
-      write(task,*) j, l
-      call add_task_to_taskserver(zmq_to_qp_run_socket,task)
-    enddo
   enddo
+
   external :: ao_bielec_integrals_in_map_slave_inproc, ao_bielec_integrals_in_map_collector
   call new_parallel_threads(ao_bielec_integrals_in_map_slave_inproc, ao_bielec_integrals_in_map_collector)
-  
+
   call end_parallel_job(zmq_to_qp_run_socket,'ao_integrals')
 
   print*, 'Sorting the map'

src/MOGuess/H_CORE_guess.irp.f

@@ -5,8 +5,6 @@ program H_CORE_guess
   END_DOC
   implicit none
   character*(64)                 :: label
-  mo_coef = ao_ortho_lowdin_coef
-  TOUCH mo_coef
   label = "Guess"
   call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
                                      size(mo_mono_elec_integral,1),  &

src/MOGuess/guess_overlap.irp.f

@@ -5,8 +5,6 @@ program guess_mimi
   implicit none
   character*(64)                 :: label
 
-  mo_coef = ao_ortho_lowdin_coef
-  TOUCH mo_coef
   label = "Guess"
   call mo_as_eigvectors_of_mo_matrix(ao_overlap,          &
                                      size(ao_overlap,1),  &

src/MOGuess/h_core_guess_routine.irp.f

@@ -4,8 +4,6 @@ subroutine hcore_guess
   END_DOC
   implicit none
   character*(64)                 :: label
-  mo_coef = ao_ortho_lowdin_coef
-  TOUCH mo_coef
   label = "Guess"
   call mo_as_eigvectors_of_mo_matrix(mo_mono_elec_integral,          &
                                      size(mo_mono_elec_integral,1),  &

src/MOGuess/mo_ortho_lowdin.irp.f

@@ -1,4 +1,3 @@
-
 BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]
   implicit none
   BEGIN_DOC
@@ -8,13 +7,9 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]
   END_DOC
   integer                        :: i,j,k,l
   double precision               :: tmp_matrix(ao_num_align,ao_num),accu
+  tmp_matrix(:,:) = 0.d0
   do j=1, ao_num
-    do i=1, ao_num
+    tmp_matrix(j,j) = 1.d0
-      tmp_matrix(i,j) = 0.d0
-    enddo
-  enddo
-  do i=1, ao_num
-    tmp_matrix(i,i) = 1.d0
   enddo
   call ortho_lowdin(ao_overlap,ao_num_align,ao_num,tmp_matrix,ao_num_align,ao_num)
   do i=1, ao_num
@@ -23,6 +18,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_coef, (ao_num_align,ao_num)]
     enddo
   enddo
 END_PROVIDER
+
 BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)]
   implicit none
   BEGIN_DOC
@@ -40,7 +36,7 @@ BEGIN_PROVIDER [double precision, ao_ortho_lowdin_overlap, (ao_num_align,ao_num)
     do j=1, ao_num
       c = 0.d0
       do l=1, ao_num
-        c = ao_ortho_lowdin_coef(j,l) * ao_overlap(k,l)
+        c +=  ao_ortho_lowdin_coef(j,l) * ao_overlap(k,l)
       enddo
       do i=1, ao_num
         ao_ortho_lowdin_overlap(i,j) += ao_ortho_lowdin_coef(i,k) * c

src/MO_Basis/ao_ortho_canonical.irp.f

@@ -1,21 +1,20 @@
-
+ BEGIN_PROVIDER [ double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
-BEGIN_PROVIDER [double precision, ao_ortho_canonical_coef, (ao_num_align,ao_num)]
+&BEGIN_PROVIDER [ integer, ao_ortho_canonical_num ]
   implicit none
   BEGIN_DOC
 ! matrix of the coefficients of the mos generated by the 
 ! orthonormalization by the S^{-1/2} canonical transformation of the aos
 ! ao_ortho_canonical_coef(i,j) = coefficient of the ith ao on the jth ao_ortho_canonical orbital
   END_DOC
-  integer                        :: i,j,k,l
+  integer :: i
-  tmp_matrix(:,:) = 0.d0
+  ao_ortho_canonical_coef(:,:) = 0.d0
-  do j=1, ao_num
+  do i=1,ao_num
-    tmp_matrix(j,j) = 1.d0
+    ao_ortho_canonical_coef(i,i) = 1.d0
   enddo
-  call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,mo_tot_num)
+  call ortho_canonical(ao_overlap,ao_num_align,ao_num,ao_ortho_canonical_coef,ao_num_align,ao_ortho_canonical_num)
-  SOFT_TOUCH mo_tot_num
 END_PROVIDER
 
-BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_num_align,ao_num)]
+BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_ortho_canonical_num,ao_ortho_canonical_num)]
   implicit none
   BEGIN_DOC
 ! overlap matrix of the ao_ortho_canonical
@@ -23,19 +22,19 @@ BEGIN_PROVIDER [double precision, ao_ortho_canonical_overlap, (ao_num_align,ao_n
   END_DOC
   integer                        :: i,j,k,l
   double precision               :: c
-  do j=1, ao_num
+  do j=1, ao_ortho_canonical_num
-    do i=1, ao_num
+    do i=1, ao_ortho_canonical_num
       ao_ortho_canonical_overlap(i,j) = 0.d0
     enddo
   enddo
-  do k=1, ao_num
+  do j=1, ao_ortho_canonical_num
-    do j=1, ao_num
+    do k=1, ao_num
       c = 0.d0
       do l=1, ao_num
-        c +=  ao_ortho_canonical_coef(j,l) * ao_overlap(k,l)
+        c +=  ao_ortho_canonical_coef(l,j) * ao_overlap(l,k)
       enddo
-      do i=1, ao_num
+      do i=1, ao_ortho_canonical_num
-        ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(i,k) * c
+        ao_ortho_canonical_overlap(i,j) += ao_ortho_canonical_coef(k,i) * c
       enddo
     enddo
   enddo

src/MO_Basis/mos.irp.f

@@ -9,8 +9,9 @@ BEGIN_PROVIDER [ integer, mo_tot_num ]
   if (exists) then
     call ezfio_get_mo_basis_mo_tot_num(mo_tot_num)
   else
-    mo_tot_num = ao_num
+    mo_tot_num = ao_ortho_canonical_num
   endif
+  call write_int(6,mo_tot_num,'mo_tot_num')
   ASSERT (mo_tot_num > 0)
 END_PROVIDER
 
@@ -56,7 +57,14 @@ END_PROVIDER
     deallocate(buffer)
   else
     ! Orthonormalized AO basis
-    mo_coef = 0.
+    do i=1,mo_tot_num
+      do j=1,ao_num
+        mo_coef(j,i) = ao_ortho_canonical_coef(j,i)
+      enddo
+      do j=ao_num+1,ao_num_align
+        mo_coef(j,i) = 0.d0
+      enddo
+    enddo
   endif
 END_PROVIDER
 

src/Utils/LinearAlgebra.irp.f

@@ -42,11 +42,81 @@ subroutine svd(A,LDA,U,LDU,D,Vt,LDVt,m,n)
 end
 
 
+subroutine ortho_canonical(overlap,LDA,N,C,LDC,m)
+  implicit none
+  BEGIN_DOC
+  ! Compute C_new=C_old.U.s^-1/2 canonical orthogonalization.
+  !
+  ! overlap : overlap matrix 
+  !
+  ! LDA : leftmost dimension of overlap array
+  !
+  ! N : Overlap matrix is NxN (array is (LDA,N) )
+  !
+  ! C : Coefficients of the vectors to orthogonalize. On exit,
+  !     orthogonal vectors
+  !
+  ! LDC : leftmost dimension of C
+  !
+  ! m : Coefficients matrix is MxN, ( array is (LDC,N) )
+  !
+  END_DOC
+  
+  integer, intent(in)            :: lda, ldc, n
+  integer, intent(out)           :: m
+  double precision, intent(in)   :: overlap(lda,n)
+  double precision, intent(inout) :: C(ldc,n)
+  double precision, allocatable  :: U(:,:)
+  double precision, allocatable  :: Vt(:,:)
+  double precision, allocatable  :: D(:)
+  double precision, allocatable  :: S_half(:,:)
+  !DEC$ ATTRIBUTES ALIGN : 64    :: U, Vt, D
+  integer                        :: info, i, j
+  
+  allocate (U(ldc,n), Vt(lda,n), D(n), S_half(lda,n))
+
+  call svd(overlap,lda,U,ldc,D,Vt,lda,n,n)
+
+  m=n
+  do i=1,n
+    if ( D(i) >= 1.d-4 ) then
+      D(i) = 1.d0/dsqrt(D(i))
+    else
+      m = i-1
+      exit
+    endif
+  enddo
+  do i=m+1,n
+    D(i) = 0.d0
+  enddo
+
+  !$OMP PARALLEL DEFAULT(NONE) &
+  !$OMP SHARED(S_half,U,D,Vt,n,C,m) &
+  !$OMP PRIVATE(i,j)
+
+  !$OMP DO
+  do j=1,n
+    do i=1,n
+      S_half(i,j) = U(i,j)*D(j)
+    enddo
+    do i=1,n
+      U(i,j) = C(i,j)
+    enddo
+  enddo
+  !$OMP END DO
+  
+  !$OMP END PARALLEL
+
+  call dgemm('N','N',n,m,n,1.d0,U,size(U,1),S_half,size(S_half,1),0.d0,C,size(C,1))
+  deallocate (U, Vt, D, S_half)
+  
+end
+
 
 subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
   implicit none
   BEGIN_DOC
-  ! Compute C_new=C_old.S^-1/2 canonical orthogonalization.
+  ! Compute C_new=C_old.S^-1/2 orthogonalization.
   !
   ! overlap : overlap matrix 
   !
@@ -81,7 +151,7 @@ subroutine ortho_lowdin(overlap,LDA,N,C,LDC,m)
 
   !$OMP DO
   do i=1,n
-    if ( D(i) < 1.d-4 ) then
+    if ( D(i) < 1.d-5 ) then
       D(i) = 0.d0
     else
       D(i) = 1.d0/dsqrt(D(i))