Merge pull request #13 from scemama/master

merge with scemama

configure

@@ -46,7 +46,12 @@ if len(sys.argv) != 3:
 # \_| |_ (_) |_) (_| |   | | | | (_)
 #
 
-QP_ROOT = os.getcwd()
+try:
+  QP_ROOT = os.environ["QP_ROOT"]
+except KeyError:
+  QP_ROOT = os.getcwd()
+  os.environ["QP_ROOT"] = QP_ROOT
+
 QP_ROOT_BIN = join(QP_ROOT, "bin")
 QP_ROOT_LIB = join(QP_ROOT, "lib")
 QP_ROOT_INSTALL = join(QP_ROOT, "install")

ocaml/qp_set_mo_class.ml

@@ -49,20 +49,20 @@ let t_to_string = function
   | None -> assert false
 ;;
 
-let run ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_filename =
+let set ~core ~inact ~act ~virt ~del =
 
-  Ezfio.set_file ezfio_filename ;
+  let mo_tot_num =
-  if not (Ezfio.has_mo_basis_mo_tot_num ()) then
+     Ezfio.get_mo_basis_mo_tot_num ()
-    failwith "mo_basis/mo_tot_num not found" ;
+  in
-
-  let mo_tot_num = Ezfio.get_mo_basis_mo_tot_num () in
   let n_int =
      try  N_int_number.of_int (Ezfio.get_determinants_n_int ())
      with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num 
   in
 
 
-  let mo_class = Array.init mo_tot_num ~f:(fun i -> None) in
+  let mo_class =
+     Array.init mo_tot_num ~f:(fun i -> None)
+  in
 
   (* Check input data *)
   let apply_class l = 
@@ -196,6 +196,49 @@ let run ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_fi
   |> Ezfio.set_bitmasks_cas; 
 ;;
 
+let get () =
+
+  let mo_tot_num =
+     Ezfio.get_mo_basis_mo_tot_num ()
+  in
+  let n_int =
+     try  N_int_number.of_int (Ezfio.get_determinants_n_int ())
+     with _ -> Bitlist.n_int_of_mo_tot_num mo_tot_num 
+  in
+
+  let bitmasks = 
+    match Input.Bitmasks.read () with
+    | Some x -> x
+    | None -> failwith "No data to print"
+  in
+  assert (bitmasks.Input.Bitmasks.n_mask_gen |> Bitmask_number.to_int = 1);
+  assert (bitmasks.Input.Bitmasks.n_mask_cas |> Bitmask_number.to_int = 1);
+
+  let (generators,cas) =
+    Bitlist.of_int64_array bitmasks.Input.Bitmasks.generators,
+    Bitlist.of_int64_array bitmasks.Input.Bitmasks.cas
+  in
+
+  Printf.printf "MO  : %d\n" mo_tot_num;
+  Printf.printf "n_int: %d\n" (N_int_number.to_int n_int);
+  Printf.printf "Gen : %s\nCAS : %s\n"
+    (Bitlist.to_string generators)
+    (Bitlist.to_string cas)
+  
+;;
+
+let run ~print ?(core="[]") ?(inact="[]") ?(act="[]") ?(virt="[]") ?(del="[]") ezfio_filename =
+
+  Ezfio.set_file ezfio_filename ;
+  if not (Ezfio.has_mo_basis_mo_tot_num ()) then
+    failwith "mo_basis/mo_tot_num not found" ;
+
+  if print then
+     get ()
+  else
+     set ~core ~inact ~act ~virt ~del
+;;
+
 let ezfio_file =
   let failure filename = 
         eprintf "'%s' is not an EZFIO file.\n%!" filename;
@@ -240,6 +283,7 @@ let spec =
   +> flag "act"    (optional string) ~doc:"range Range of active orbitals"
   +> flag "virt"   (optional string) ~doc:"range Range of virtual orbitals"
   +> flag "del"    (optional string) ~doc:"range Range of deleted orbitals"
+  +> flag "print"  no_arg ~doc:" Print the current masks"
   +> anon ("ezfio_filename" %: ezfio_file)
 ;;
 
@@ -251,7 +295,7 @@ let command =
       The range of MOs has the form : \"[36-53,72-107,126-131]\"
         ")
     spec
-    (fun core inact act virt del ezfio_filename () -> run ?core ?inact ?act ?virt ?del ezfio_filename )
+    (fun core inact act virt del print ezfio_filename () -> run ~print ?core ?inact ?act ?virt ?del ezfio_filename )
 ;;
 
 let () =

plugins/FOBOCI/H_apply_dressed_autonom.irp.f

@@ -273,7 +273,8 @@ subroutine H_apply_dressed_pert_monoexc(key_in, hole_1,particl_1,i_generator,ipr
   integer,parameter              :: size_max = 3072
   
     integer, intent(in) :: Ndet_generators
-    double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators),E_ref
+    double precision, intent(in) :: E_ref
+    double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
     integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
 
   integer          ,intent(in)   :: i_generator
@@ -437,8 +438,9 @@ subroutine H_apply_dressed_pert(delta_ij_generators_,  Ndet_generators,psi_det_g
   
   
     integer, intent(in) :: Ndet_generators
-    integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators),E_ref
+    double precision, intent(in) :: E_ref
-    double precision, intent(in) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
+    double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
+    integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
 
   
   integer                        :: i_generator, nmax

plugins/Full_CI/H_apply.irp.f

@@ -9,7 +9,7 @@ print s
 
 s = H_apply_zmq("FCI_PT2")
 s.set_perturbation("epstein_nesbet_2x2")
-#s.unset_openmp()
+s.unset_openmp()
 print s
 
 s = H_apply("FCI_no_skip")

plugins/MRCC_Utils/mrcc_dress.irp.f

@@ -274,7 +274,7 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
       enddo
       
       do i_state=1,Nstates
-        ci_inv(i_state) = 1.d0/psi_ref_coef(i_I,i_state)
+        ci_inv(i_state) = psi_ref_coef_inv(i_I,i_state)
       enddo
       do l_sd=1,idx_alpha(0)
         k_sd = idx_alpha(l_sd)
@@ -285,16 +285,20 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Nstates, Ndet_non_ref, Ndet_ref,i_ge
         enddo
       enddo
       call omp_set_lock( psi_ref_lock(i_I) )
-      do l_sd=1,idx_alpha(0)
+      do i_state=1,Nstates
-        k_sd = idx_alpha(l_sd)
+        if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
-        do i_state=1,Nstates
+          do l_sd=1,idx_alpha(0)
-          delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
+            k_sd = idx_alpha(l_sd)
-          if(dabs(psi_ref_coef(i_I,i_state)).ge.5.d-5)then
+              delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
-            delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
+              delta_ii_(i_state,i_I) = delta_ii_(i_state,i_I) - dIa_hla(i_state,k_sd) * ci_inv(i_state) * psi_non_ref_coef_transp(i_state,k_sd)
-          else
+          enddo
-            delta_ii_(i_state,i_I)  = 0.d0
+        else
-          endif
+          delta_ii_(i_state,i_I)  = 0.d0
-        enddo
+          do l_sd=1,idx_alpha(0)
+            k_sd = idx_alpha(l_sd)
+              delta_ij_(i_state,k_sd,i_I) = delta_ij_(i_state,k_sd,i_I) + dIa_hla(i_state,k_sd)
+          enddo
+        endif
       enddo
       call omp_unset_lock( psi_ref_lock(i_I) )
     enddo

plugins/Psiref_CAS/overwrite_with_cas.irp.f

@@ -0,0 +1,3 @@
+program overwrite_w_cas
+  call extract_ref
+end

plugins/Psiref_CAS/psi_ref.irp.f

@@ -26,21 +26,21 @@ use bitmasks
 
 END_PROVIDER
 
-BEGIN_PROVIDER [ double precision, psi_ref_coef_normalized,  (psi_det_size,n_states) ]
+BEGIN_PROVIDER [ double precision, psi_ref_coef_inv, (psi_det_size,n_states) ]
  implicit none
  BEGIN_DOC
-! Normalized coefficients of the reference
+ ! 1/psi_ref_coef
  END_DOC
- integer :: i,j,k
+ integer :: i, i_state
- do k=1,N_states
+ do i_state=1,N_states
-   do j=1,N_det_ref
+  do i=1,N_det_ref
-     psi_ref_coef_normalized(j,k) = psi_ref_coef(j,k)
+    psi_ref_coef_inv(i,i_state) = 1.d0/psi_ref_coef(i,i_state)
-   enddo
+  enddo
-   call normalize(psi_ref_coef_normalized(1,k), N_det_ref)
  enddo
 
 END_PROVIDER
 
+
  BEGIN_PROVIDER [ integer(bit_kind), psi_ref_restart, (N_int,2,psi_det_size) ]
 &BEGIN_PROVIDER [ double precision, psi_ref_coef_restart,  (psi_det_size,n_states) ]
   implicit none

plugins/Psiref_Utils/NEEDED_CHILDREN_MODULES

@@ -0,0 +1 @@
+Bitmask Determinants

plugins/Psiref_Utils/extract_ref.irp.f

@@ -0,0 +1,24 @@
+subroutine extract_ref
+  implicit none
+  BEGIN_DOC
+  ! Replaces the total wave function by the normalized projection on the reference
+  END_DOC
+  
+  integer                        :: i,j,k
+  do k=1,N_states
+    do j=1,N_det_ref
+      psi_coef(j,k) = psi_ref_coef_normalized(j,k)
+    enddo
+  enddo
+
+  do j=1,N_det_ref
+    do k=1,N_int
+      psi_det(k,1,j) = psi_ref(k,1,j) 
+      psi_det(k,2,j) = psi_ref(k,2,j) 
+    enddo
+  enddo
+  N_det = N_det_ref
+
+  call save_wavefunction
+
+end

plugins/Psiref_Utils/psi_ref_utils.irp.f

@@ -27,6 +27,22 @@ BEGIN_PROVIDER [ double precision, psi_ref_coef_transp, (n_states,psi_det_size)
  enddo
 END_PROVIDER
 
+BEGIN_PROVIDER [ double precision, psi_ref_coef_normalized,  (psi_det_size,n_states) ]
+ implicit none
+ BEGIN_DOC
+! Normalized coefficients of the reference
+ END_DOC
+ integer :: i,j,k
+ do k=1,N_states
+   do j=1,N_det_ref
+     psi_ref_coef_normalized(j,k) = psi_ref_coef(j,k)
+   enddo
+   call normalize(psi_ref_coef_normalized(1,k), N_det_ref)
+ enddo
+
+END_PROVIDER
+
+
 BEGIN_PROVIDER [ double precision, psi_non_ref_coef_transp, (n_states,psi_det_size) ]
  implicit none
  BEGIN_DOC

scripts/compilation/qp_create_ninja.py

@@ -184,7 +184,7 @@ def ninja_ezfio_config_rule():
 
 def get_children_of_ezfio_cfg(l_module_with_ezfio_cfg):
     """
-    From a module list of ezfio_cfg return all the stuff create by him
+    From a module list of ezfio_cfg return all the stuff created by it
     """
     config_folder = join(QP_EZFIO, "config")
 

src/Determinants/davidson.irp.f

@@ -315,6 +315,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
   double precision, intent(inout) :: u_in(dim_in,N_st)
   double precision, intent(out)  :: energies(N_st)
   
+  integer                        :: sze_8
   integer                        :: iter
   integer                        :: i,j,k,l,m
   logical                        :: converged
@@ -334,6 +335,7 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
   double precision               :: to_print(2,N_st)
   double precision               :: cpu, wall
   
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: U, W, R, Wt, y, h, lambda
 
 
   call write_time(iunit)
@@ -362,12 +364,15 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
   enddo
   write(iunit,'(A)') trim(write_buffer)
 
+  integer, external :: align_double
+  sze_8 = align_double(sze)
+
   allocate(                                                          &
       kl_pairs(2,N_st*(N_st+1)/2),                                   &
-      W(sze,N_st,davidson_sze_max),                                                   &
+      W(sze_8,N_st,davidson_sze_max),                                                   &
       Wt(sze),                                                        &
-      U(sze,N_st,davidson_sze_max),                                  &
+      U(sze_8,N_st,davidson_sze_max),                                  &
-      R(sze,N_st),                                                   &
+      R(sze_8,N_st),                                                   &
       h(N_st,davidson_sze_max,N_st,davidson_sze_max),                &
       y(N_st,davidson_sze_max,N_st,davidson_sze_max),                &
       lambda(N_st*davidson_sze_max))
@@ -473,7 +478,10 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
       ! Express eigenvectors of h in the determinant basis
       ! --------------------------------------------------
       
+     !$OMP PARALLEL DEFAULT(NONE)                                     &
+        !$OMP PRIVATE(k,i,l,iter2) SHARED(U,W,R,y,iter,lambda,N_st,sze)
       do k=1,N_st
+        !$OMP DO
         do i=1,sze
           U(i,k,iter+1) = 0.d0
           W(i,k,iter+1) = 0.d0
@@ -484,7 +492,9 @@ subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iun
             enddo
           enddo
         enddo
+        !$OMP END DO
       enddo
+      !$OMP END PARALLEL
       
       ! Compute residual vector
       ! -----------------------

src/Determinants/slater_rules.irp.f

@@ -1608,12 +1608,11 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
   integer                        :: i,j,k,l, jj,ii
   integer                        :: i0, j0
   
-  integer, allocatable           :: shortcut(:), sort_idx(:)
+  integer, allocatable           :: shortcut(:,:), sort_idx(:,:)
-  integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
+  integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
   integer(bit_kind)              :: sorted_i(Nint)
   
   integer                        :: sh, sh2, ni, exa, ext, org_i, org_j, endi
-  double precision               :: local_threshold
   
 
   ASSERT (Nint > 0)
@@ -1621,104 +1620,83 @@ subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint)
   ASSERT (n>0)
   PROVIDE ref_bitmask_energy davidson_criterion
 
-  allocate (shortcut(0:n+1), sort_idx(n), sorted(Nint,n), version(Nint,n))
+  allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
   v_0 = 0.d0
 
-  call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
+  call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
+  call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
   
   !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold,sorted_i)&
+      !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,sorted_i)&
-      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version)
+      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,sorted,shortcut,sort_idx,version)
   allocate(vt(n))
   Vt = 0.d0
   
   !$OMP DO SCHEDULE(dynamic)
-  do sh=1,shortcut(0)
+  do sh=1,shortcut(0,1)
-    do sh2=1,sh
+    do sh2=sh,shortcut(0,1)
       exa = 0
       do ni=1,Nint
-        exa = exa + popcnt(xor(version(ni,sh), version(ni,sh2)))
+        exa = exa + popcnt(xor(version(ni,sh,1), version(ni,sh2,1)))
       end do
       if(exa > 2) then
         cycle
       end if
       
-      do i=shortcut(sh),shortcut(sh+1)-1
+      do i=shortcut(sh,1),shortcut(sh+1,1)-1
-        org_i = sort_idx(i)
+        org_i = sort_idx(i,1)
-        local_threshold = threshold_davidson - dabs(u_0(org_i))
         if(sh==sh2) then
           endi = i-1
         else
-          endi = shortcut(sh2+1)-1
+          endi = shortcut(sh2+1,1)-1
         end if
         do ni=1,Nint
-          sorted_i(ni) = sorted(ni,i)
+          sorted_i(ni) = sorted(ni,i,1)
         enddo
         
-        do j=shortcut(sh2),endi
+        do j=shortcut(sh2,1),endi
-          org_j = sort_idx(j)
+          org_j = sort_idx(j,1)
-          if ( dabs(u_0(org_j)) > local_threshold ) then
+          ext = exa
-            ext = exa
+          do ni=1,Nint
-            do ni=1,Nint
+            ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j,1)))
-              ext = ext + popcnt(xor(sorted_i(ni), sorted(ni,j)))
+          end do
-            end do
+          if(ext <= 4) then
-            if(ext <= 4) then
+            call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
-              call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
+            vt (org_i) = vt (org_i) + hij*u_0(org_j)
-              vt (org_i) = vt (org_i) + hij*u_0(org_j)
+            vt (org_j) = vt (org_j) + hij*u_0(org_i)
-              vt (org_j) = vt (org_j) + hij*u_0(org_i)
-            endif
           endif
         enddo
       enddo
     enddo
   enddo
-  !$OMP END DO
+  !$OMP END DO NOWAIT
-  
-  !$OMP CRITICAL
-  do i=1,n
-    v_0(i) = v_0(i) + vt(i)
-  enddo
-  !$OMP END CRITICAL
-  
-  deallocate(vt)
-  !$OMP END PARALLEL
-  
-  call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut, version, n, Nint)
-  
-  !$OMP PARALLEL DEFAULT(NONE)                                       &
-      !$OMP PRIVATE(i,hij,j,k,jj,vt,ii,sh,sh2,ni,exa,ext,org_i,org_j,endi,local_threshold)&
-      !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0,threshold_davidson,sorted,shortcut,sort_idx,version)
-  allocate(vt(n))
-  Vt = 0.d0
   
   !$OMP DO SCHEDULE(dynamic)
-  do sh=1,shortcut(0)
+  do sh=1,shortcut(0,2)
-    do i=shortcut(sh),shortcut(sh+1)-1
+    do i=shortcut(sh,2),shortcut(sh+1,2)-1
-      org_i = sort_idx(i)
+      org_i = sort_idx(i,2)
-      local_threshold = threshold_davidson - dabs(u_0(org_i))
+      do j=shortcut(sh,2),i-1
-      do j=shortcut(sh),i-1
+        org_j = sort_idx(j,2)
-        org_j = sort_idx(j)
+        ext = 0
-        if ( dabs(u_0(org_j)) > local_threshold ) then
+        do ni=1,Nint
-          ext = 0
+          ext = ext + popcnt(xor(sorted(ni,i,2), sorted(ni,j,2)))
-          do ni=1,Nint
+        end do
-            ext = ext + popcnt(xor(sorted(ni,i), sorted(ni,j)))
+        if(ext == 4) then
-          end do
+          call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
-          if(ext == 4) then
+          vt (org_i) = vt (org_i) + hij*u_0(org_j)
-            call i_H_j(keys_tmp(1,1,org_j),keys_tmp(1,1,org_i),Nint,hij)
+          vt (org_j) = vt (org_j) + hij*u_0(org_i)
-            vt (org_i) = vt (org_i) + hij*u_0(org_j)
-            vt (org_j) = vt (org_j) + hij*u_0(org_i)
-          end if
         end if
       end do
     end do
   enddo
-  !$OMP END DO
+  !$OMP END DO NOWAIT
   
   !$OMP CRITICAL
-  do i=1,n
+  do i=n,1,-1
     v_0(i) = v_0(i) + vt(i)
   enddo
   !$OMP END CRITICAL
+
   deallocate(vt)
   !$OMP END PARALLEL
   

src/Integrals_Bielec/ao_bielec_integrals_in_map_slave.irp.f

@@ -67,6 +67,8 @@ end
 
 
 
+
+
 subroutine ao_bielec_integrals_in_map_slave(thread,iproc)
   use map_module
   use f77_zmq
@@ -107,7 +109,7 @@ subroutine ao_bielec_integrals_in_map_slave(thread,iproc)
       call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, 0)
     enddo
     call compute_ao_integrals_jl(l,l,n_integrals,buffer_i,buffer_value) 
-    call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id,n_integrals)
+    call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
     call push_integrals(zmq_socket_push, n_integrals, buffer_i, buffer_value, task_id)
   enddo
 
@@ -127,7 +129,7 @@ subroutine pull_integrals(zmq_socket_pull, n_integrals, buffer_i, buffer_value,
   BEGIN_DOC
   ! How the collector pulls the computed integrals
   END_DOC
-  integer(ZMQ_PTR), intent(out)  :: zmq_socket_pull
+  integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
   integer, intent(out)           :: n_integrals
   integer(key_kind), intent(out) :: buffer_i(*)
   real(integral_kind), intent(out) :: buffer_value(*)

src/Integrals_Bielec/qp_ao_ints.irp.f

@@ -8,7 +8,8 @@ program qp_ao_ints
 
   call switch_qp_run_to_master
 
-  PROVIDE zmq_context
+  zmq_context = f77_zmq_ctx_new ()
+
   ! Set the state of the ZMQ
   zmq_state = 'ao_integrals'
 

src/MO_Basis/swap_mos.irp.f

@@ -0,0 +1,14 @@
+program swap_mos
+  implicit none
+  integer :: i,j, i1, i2
+  double precision :: x
+  print *,  'MOs to swap?'
+  read(*,*) i1, i2
+  do i=1,ao_num_align
+    x = mo_coef(i,i1)
+    mo_coef(i,i1) = mo_coef(i,i2)
+    mo_coef(i,i2) = x
+  enddo
+  call save_mos
+
+end

src/Utils/util.irp.f

@@ -324,6 +324,7 @@ double precision function u_dot_v(u,v,sze)
   t3 = t2+t2
   t4 = t3+t2
   u_dot_v = 0.d0
+  !DIR$ VECTOR ALWAYS
   do i=1,t2
     u_dot_v = u_dot_v + u(t1+i)*v(t1+i) + u(t2+i)*v(t2+i) +          &
         u(t3+i)*v(t3+i) + u(t4+i)*v(t4+i)
@@ -359,6 +360,7 @@ double precision function u_dot_u(u,sze)
 !   u_dot_u = u_dot_u+u(i)*u(i)
 ! enddo
   
+  !DIR$ VECTOR ALWAYS
   do i=1,sze
     u_dot_u = u_dot_u + u(i)*u(i)
   enddo

src/ZMQ/utils.irp.f

@@ -361,6 +361,8 @@ subroutine end_zmq_pull_socket(zmq_socket_pull)
     stop 'error'
   endif
   
+  call sleep(1) ! see https://github.com/zeromq/libzmq/issues/1922
+
   rc = f77_zmq_setsockopt(zmq_socket_pull,ZMQ_LINGER,0,4)
   if (rc /= 0) then
     stop 'Unable to set ZMQ_LINGER on zmq_socket_pull'