merged

.travis.yml

@@ -24,7 +24,7 @@ python:
 
 script: 
   - ./configure  --production ./config/gfortran.cfg
-  - source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles
+  - source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD mrcepa0 All_singles
   - source ./quantum_package.rc ; ninja
   - source ./quantum_package.rc ; cd ocaml ; make ; cd -
   - source ./quantum_package.rc ; cd tests ; ./run_tests.sh #-v

README.md

@@ -7,11 +7,14 @@ Set of quantum chemistry programs and libraries.
 
 For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>), or below for the installation instructions.
 
+
+
 Demo
 ====
 
 [![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_295x166.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo")
 
+[![Frozen-core Full-CI energy of Ti](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png)](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png "Convergence of Ti in cc-pv{DTQ}Z")
 
 # Installation
 
@@ -155,7 +158,7 @@ Program exited with code 139.
 
 #### Why ?
 
-It's caused when we call the DGEM routine of LAPACK. 
+It's caused when we call the DGEMM routine of LAPACK. 
 
 ##### Fix
 

plugins/FOBOCI/H_apply_dressed_autonom.irp.f

@@ -273,7 +273,7 @@ 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) :: E_ref
+    double precision, intent(inout) :: 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)
 
@@ -438,7 +438,7 @@ subroutine H_apply_dressed_pert(delta_ij_generators_,  Ndet_generators,psi_det_g
   
   
     integer, intent(in) :: Ndet_generators
-    double precision, intent(in) :: E_ref
+    double precision, intent(inout) :: 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)
 

plugins/FOBOCI/NEEDED_CHILDREN_MODULES

@@ -1 +1 @@
-Perturbation Selectors_no_sorted Hartree_Fock Davidson
+Perturbation Selectors_no_sorted Hartree_Fock Davidson CISD

plugins/FOBOCI/hcc_1h1p.irp.f

@@ -15,11 +15,10 @@ subroutine routine
   call diagonalize_CI
   call test_hcc
   call test_mulliken
-! call SC2_1h1p(psi_det,psi_coef,energies, &
-!       diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
   allocate(H_matrix(N_det,N_det))
-  call SC2_1h1p_full(psi_det,psi_coef,energies, &
+  stop 'SC2_1h1p_full is not in the git!'
-        H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
+!  call SC2_1h1p_full(psi_det,psi_coef,energies, &
+!        H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
   deallocate(H_matrix)
   integer :: i,j
   double precision :: accu,coef_hf

plugins/Full_CI_ZMQ/fci_zmq.irp.f

@@ -46,14 +46,16 @@ program fci_zmq
     PROVIDE  psi_det_sorted
 
     call diagonalize_CI
+    call save_wavefunction
     
     if (N_det > N_det_max) then
       psi_det = psi_det_sorted
       psi_coef = psi_coef_sorted
       N_det = N_det_max
       soft_touch N_det psi_det psi_coef
-    endif
+      call diagonalize_CI
       call save_wavefunction
+    endif
     
     print *,  'N_det          = ', N_det
     print *,  'N_states       = ', N_states
@@ -79,17 +81,13 @@ program fci_zmq
     E_CI_before(1:N_states) = CI_energy(1:N_states)
     call ezfio_set_full_ci_energy(CI_energy)
   enddo
-  if (N_det > N_det_max) then
+
-    N_det = N_det_max
-    touch N_det psi_det psi_coef
-    call diagonalize_CI
-  endif
   if(do_pt2_end)then
     print*,'Last iteration only to compute the PT2'
     threshold_selectors = 1.d0
     threshold_generators = 0.9999d0
     E_CI_before(1:N_states) = CI_energy(1:N_states)
-    call ZMQ_selection(1, pt2)
+    call ZMQ_selection(0, pt2)
     print *,  'Final step'
     print *,  'N_det    = ', N_det
     print *,  'N_states = ', N_states
@@ -108,7 +106,7 @@ end
 
 
 
-subroutine ZMQ_selection(N, pt2)
+subroutine ZMQ_selection(N_in, pt2)
   use f77_zmq
   use selection_types
   
@@ -116,13 +114,14 @@ subroutine ZMQ_selection(N, pt2)
   
   character*(512)                :: task 
   integer(ZMQ_PTR)               :: zmq_to_qp_run_socket 
-  integer, intent(in)            :: N
+  integer, intent(in)            :: N_in
   type(selection_buffer)         :: b
-  integer                        :: i
+  integer                        :: i, N
   integer, external              :: omp_get_thread_num
   double precision, intent(out)  :: pt2(N_states)
   
   
+  N = max(N_in,1)
   provide nproc
   provide ci_electronic_energy
   call new_parallel_job(zmq_to_qp_run_socket,"selection")
@@ -147,16 +146,18 @@ subroutine ZMQ_selection(N, pt2)
       if (i==0) then
         call selection_collector(b, pt2)
       else
-        call selection_dressing_slave_inproc(i)
+        call selection_slave_inproc(i)
       endif
   !$OMP END PARALLEL
   call end_parallel_job(zmq_to_qp_run_socket, 'selection') 
+  if (N_in > 0) then
     call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN
     call copy_H_apply_buffer_to_wf()
+  endif
 end subroutine
 
 
-subroutine selection_dressing_slave_inproc(i)
+subroutine selection_slave_inproc(i)
   implicit none
   integer, intent(in)            :: i
 

plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f

@@ -0,0 +1,107 @@
+program selection_slave
+  implicit none
+  BEGIN_DOC
+! Helper program to compute the PT2 in distributed mode.
+  END_DOC
+
+  read_wf = .False.
+  SOFT_TOUCH read_wf
+  call provide_everything
+  call switch_qp_run_to_master
+  call run_wf
+end
+
+subroutine provide_everything
+  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
+!   PROVIDE ci_electronic_energy mo_tot_num N_int
+end
+
+subroutine run_wf
+  use f77_zmq
+  implicit none
+
+  integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
+  integer(ZMQ_PTR) :: zmq_to_qp_run_socket
+  double precision :: energy(N_states_diag)
+  character*(64) :: states(2)
+  integer :: rc, i
+  
+  call provide_everything
+  
+  zmq_context = f77_zmq_ctx_new ()
+  states(1) = 'selection'
+  states(2) = 'davidson'
+
+  zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
+
+  do
+
+    call wait_for_states(states,zmq_state,2)
+
+    if(trim(zmq_state) == 'Stopped') then
+
+      exit
+
+    else if (trim(zmq_state) == 'selection') then
+
+      ! Selection
+      ! ---------
+
+      print *,  'Selection'
+      call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states_diag)
+  
+      !$OMP PARALLEL PRIVATE(i)
+      i = omp_get_thread_num()
+      call selection_slave_tcp(i, energy)
+      !$OMP END PARALLEL
+      print *,  'Selection done'
+
+    else if (trim(zmq_state) == 'davidson') then
+
+      ! Davidson
+      ! --------
+
+      print *,  'Davidson'
+      call davidson_miniserver_get()
+      !$OMP PARALLEL PRIVATE(i)
+      i = omp_get_thread_num()
+      call davidson_slave_tcp(i)
+      !$OMP END PARALLEL
+      print *,  'Davidson done'
+
+    endif
+
+  end do
+end
+
+subroutine update_energy(energy)
+  implicit none
+  double precision, intent(in) :: energy(N_states_diag)
+  BEGIN_DOC
+! Update energy when it is received from ZMQ
+  END_DOC
+  integer :: j,k
+  do j=1,N_states
+    do k=1,N_det
+      CI_eigenvectors(k,j) = psi_coef(k,j)
+    enddo
+  enddo
+  call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int)
+  if (.True.) then
+    do k=1,size(ci_electronic_energy)
+      ci_electronic_energy(k) = energy(k)
+    enddo
+    TOUCH ci_electronic_energy CI_eigenvectors_s2 CI_eigenvectors
+  endif
+
+  call write_double(6,ci_energy,'Energy')
+end
+
+subroutine selection_slave_tcp(i,energy)
+  implicit none
+  double precision, intent(in) :: energy(N_states_diag)
+  integer, intent(in)            :: i
+
+  call run_selection_slave(0,i,energy)
+end
+

plugins/Full_CI_ZMQ/selection_slave.irp.f

@@ -23,20 +23,19 @@ subroutine run_wf
   integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
   integer(ZMQ_PTR) :: zmq_to_qp_run_socket
   double precision :: energy(N_states_diag)
-  character*(64) :: states(2)
+  character*(64) :: states(1)
   integer :: rc, i
   
   call provide_everything
   
   zmq_context = f77_zmq_ctx_new ()
   states(1) = 'selection'
-  states(2) = 'davidson'
 
   zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
 
   do
 
-    call wait_for_states(states,zmq_state,2)
+    call wait_for_states(states,zmq_state,1)
 
     if(trim(zmq_state) == 'Stopped') then
 
@@ -52,24 +51,10 @@ subroutine run_wf
   
       !$OMP PARALLEL PRIVATE(i)
       i = omp_get_thread_num()
-      call selection_dressing_slave_tcp(i, energy)
+      call selection_slave_tcp(i, energy)
       !$OMP END PARALLEL
       print *,  'Selection done'
 
-    else if (trim(zmq_state) == 'davidson') then
-
-      ! Davidson
-      ! --------
-
-      print *,  'Davidson'
-      call davidson_miniserver_get()
-
-      !!$OMP PARALLEL PRIVATE(i)
-      !i = omp_get_thread_num()
-      call davidson_slave_tcp(0)
-      !!$OMP END PARALLEL
-      print *,  'Davidson done'
-
     endif
 
   end do
@@ -98,7 +83,7 @@ subroutine update_energy(energy)
   call write_double(6,ci_energy,'Energy')
 end
 
-subroutine selection_dressing_slave_tcp(i,energy)
+subroutine selection_slave_tcp(i,energy)
   implicit none
   double precision, intent(in) :: energy(N_states_diag)
   integer, intent(in)            :: i

plugins/MRCC_CASSD/.gitignore

@@ -1,34 +0,0 @@
-# Automatically created by $QP_ROOT/scripts/module/module_handler.py 
-.ninja_deps
-.ninja_log
-AO_Basis
-Bitmask
-Determinants
-Electrons
-Ezfio_files
-Generators_full
-Hartree_Fock
-IRPF90_man
-IRPF90_temp
-Integrals_Bielec
-Integrals_Monoelec
-MOGuess
-MO_Basis
-MRCC_Utils
-Makefile
-Makefile.depend
-Nuclei
-Perturbation
-Properties
-Pseudo
-Psiref_CAS
-Psiref_Utils
-Selectors_full
-Utils
-ZMQ
-ezfio_interface.irp.f
-irpf90.make
-irpf90_entities
-mrcc_cassd
-mrcc_noiter
-tags

plugins/MRCC_CASSD/EZFIO.cfg

@@ -1,17 +0,0 @@
-[energy]
-type: double precision
-doc: Calculated energy
-interface: ezfio
-
-[thresh_mrcc]
-type: Threshold
-doc: Threshold on the convergence of the MRCC energy
-interface: ezfio,provider,ocaml
-default: 1.e-5
-
-[n_it_mrcc_max]
-type: Strictly_positive_int
-doc: Maximum number of MRCC iterations
-interface: ezfio,provider,ocaml
-default: 10
-

plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES

@@ -1 +0,0 @@
-Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils

plugins/MRCC_CASSD/README.rst

@@ -1,80 +0,0 @@
-===========
-MRCC Module
-===========
-
-MRCC as a coupled cluster on a CAS+SD wave function.
-
-Needed Modules
-==============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-.. image:: tree_dependency.png
-
-* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
-* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
-* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
-* `Psiref_CAS <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_CAS>`_
-* `MRCC_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils>`_
-
-Documentation
-=============
-
-.. Do not edit this section. It was auto-generated from the
-.. by the `update_README.py` script.
-
-`mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_CASSD/mrcc_cassd.irp.f#L1>`_
-  Undocumented
-
-
-`print_cas_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_CASSD/mrcc_cassd.irp.f#L11>`_
-  Undocumented
-
-Needed Modules
-==============
-.. Do not edit this section It was auto-generated
-.. by the `update_README.py` script.
-
-
-.. image:: tree_dependency.png
-
-* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation>`_
-* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
-* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Generators_full>`_
-* `Psiref_CAS <http://github.com/LCPQ/quantum_package/tree/master/plugins/Psiref_CAS>`_
-* `MRCC_Utils <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils>`_
-
-Documentation
-=============
-.. Do not edit this section It was auto-generated
-.. by the `update_README.py` script.
-
-
-`mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_cassd.irp.f#L1>`_
-  Undocumented
-
-
-`mrcc_noiter <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L1>`_
-  Undocumented
-
-
-`n_it_mrcc_max <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/ezfio_interface.irp.f#L6>`_
-  Maximum number of MRCC iterations
-
-
-`print_cas_coefs <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L78>`_
-  Undocumented
-
-
-`run <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L17>`_
-  Undocumented
-
-
-`run_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L36>`_
-  Undocumented
-
-
-`thresh_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/ezfio_interface.irp.f#L25>`_
-  Threshold on the convergence of the MRCC energy
-

plugins/MRCC_CASSD/mrcc_cassd.irp.f

@@ -1,120 +0,0 @@
-program mrcc
-  implicit none
-  double precision, allocatable  :: energy(:)
-  allocate (energy(N_states))
-  
-  read_wf = .True.
-  SOFT_TOUCH read_wf
-  call print_cas_coefs
-  call set_generators_bitmasks_as_holes_and_particles
-  call run(N_states,energy)
-  if(do_pt2_end)then
-    call run_pt2(N_states,energy)
-  endif
-  deallocate(energy)
-end
-
-subroutine run(N_st,energy)
-  implicit none
-  
-  integer, intent(in) :: N_st 
-  double precision, intent(out) :: energy(N_st) 
-
-  integer :: i
-
-  double precision :: E_new, E_old, delta_e
-  integer :: iteration
-  double precision :: E_past(4), lambda
-  E_new = 0.d0
-  delta_E = 1.d0
-  iteration = 0
-  lambda = 1.d0
-  do while (delta_E > thresh_mrcc)
-    iteration += 1
-    print *,  '===========================' 
-    print *,  'MRCC Iteration', iteration
-    print *,  '===========================' 
-    print *,  ''
-    E_old = sum(ci_energy_dressed(1:N_st))
-    call write_double(6,ci_energy_dressed(1),"MRCC energy")
-    call diagonalize_ci_dressed(lambda)
-    E_new = sum(ci_energy_dressed(1:N_st))
-    delta_E = dabs(E_new - E_old)
-    call save_wavefunction
-    call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
-    if (iteration > n_it_mrcc_max) then
-      exit
-    endif
-  enddo
-  call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
-  energy(1:N_st) = ci_energy_dressed(1:N_st)
-
-end
-
-
-subroutine run_pt2(N_st,energy) 
-  implicit none 
-  integer :: i,j,k 
-  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:) 
-  integer, intent(in)          :: N_st 
-  double precision, intent(in) :: energy(N_st) 
-  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st)) 
-  pt2 = 0.d0 
- 
-  print*,'Last iteration only to compute the PT2' 
-  threshold_selectors = 1.d0 
-  threshold_generators = 0.999d0 
- 
-  N_det_generators = lambda_mrcc_pt2(0) + N_det_cas
-  do i=1,N_det_cas
-    do k=1,N_int 
-      psi_det_generators(k,1,i) = psi_ref(k,1,i) 
-      psi_det_generators(k,2,i) = psi_ref(k,2,i) 
-    enddo 
-    do k=1,N_st 
-      psi_coef_generators(i,k) = psi_ref_coef(i,k) 
-    enddo 
-  enddo 
-  do i=N_det_cas+1,N_det_generators
-    j = lambda_mrcc_pt2(i)  
-    do k=1,N_int 
-      psi_det_generators(k,1,i) = psi_non_ref(k,1,j) 
-      psi_det_generators(k,2,i) = psi_non_ref(k,2,j) 
-    enddo 
-    do k=1,N_st 
-      psi_coef_generators(i,k) = psi_non_ref_coef(j,k) 
-    enddo 
-  enddo 
-  SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed 
- 
- 
-  call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag,  N_st) 
-  print *,  'Final step' 
-  print *,  'N_det    = ', N_det 
-  print *,  'N_states = ', N_states 
-  print *,  'PT2      = ', pt2 
-  print *,  'E        = ', energy 
-  print *,  'E+PT2    = ', energy+pt2 
-  print *,  '-----' 
- 
-
-  call ezfio_set_full_ci_energy_pt2(energy+pt2)
-  deallocate(pt2,norm_pert)
-
-end 
-
-
-subroutine print_cas_coefs
-  implicit none
-
-  integer :: i,j
-  print *,  'CAS'
-  print *,  '==='
-  do i=1,N_det_cas
-    print *,  psi_cas_coef(i,:)
-    call debug_det(psi_cas(1,1,i),N_int)
-  enddo
-  call write_double(6,ci_energy(1),"Initial CI energy")
-
-end
-

plugins/MRCC_CASSD/mrcc_noiter.irp.f

@@ -1,91 +0,0 @@
-program mrcc_noiter
-  implicit none
-  double precision, allocatable :: energy(:)
-  allocate (energy(N_states))
-  read_wf = .True.
-  threshold_generators = .9999d0
-  SOFT_TOUCH read_wf threshold_generators
-  call print_cas_coefs
-  call set_generators_bitmasks_as_holes_and_particles
-  call run(N_states,energy)
-  if(do_pt2_end)then
-    call run_pt2(N_states,energy)
-  endif
-  deallocate(energy)
-end
-
-subroutine run(N_st,energy)
-  implicit none
-  
-  integer, intent(in) :: N_st
-  double precision, intent(out) :: energy(N_st)
-  integer :: i,j
-  do j=1,N_states_diag
-    do i=1,N_det
-      psi_coef(i,j) = CI_eigenvectors_dressed(i,j) 
-    enddo
-  enddo
-  SOFT_TOUCH psi_coef ci_energy_dressed
-  call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
-  call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
-  call save_wavefunction
-  energy(:) = ci_energy_dressed(:)
-end
-
-
-subroutine run_pt2(N_st,energy)
-  implicit none
-  integer :: i,j,k
-  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
-  integer, intent(in)          :: N_st
-  double precision, intent(in) :: energy(N_st)
-  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
-  pt2 = 0.d0
-
-  print*,'Last iteration only to compute the PT2'
-  threshold_selectors = 1.d0
-  threshold_generators = 0.999d0
-
-  N_det_generators = lambda_mrcc_pt2(0)
-  do i=1,N_det_generators
-    j = lambda_mrcc_pt2(i) 
-    do k=1,N_int
-      psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
-      psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
-    enddo
-    do k=1,N_st
-      psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
-    enddo
-  enddo
-  SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
-
-
-  call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag,  N_st)
-  print *,  'Final step'
-  print *,  'N_det    = ', N_det
-  print *,  'N_states = ', N_states
-  print *,  'PT2      = ', pt2
-  print *,  'E        = ', energy
-  print *,  'E+PT2    = ', energy+pt2
-  print *,  '-----'
-
-  call ezfio_set_full_ci_energy_pt2(energy+pt2)
-  deallocate(pt2,norm_pert)
-
-end
-
-
-subroutine print_cas_coefs
-  implicit none
-
-  integer :: i,j
-  print *,  'CAS'
-  print *,  '==='
-  do i=1,N_det_cas
-    print *,  psi_cas_coef(i,:)
-    call debug_det(psi_cas(1,1,i),N_int)
-  enddo
-  call write_double(6,ci_energy(1),"Initial CI energy")
-
-end
-

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
   integer                        :: j
   character*(8)                  :: st
   call write_time(output_determinants)
-  do j=1,min(N_det,N_states_diag)
+  do j=1,min(N_det,N_states)
     write(st,'(I4)') j
     CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion
     call write_double(output_determinants,CI_energy_dressed(j),'Energy of state '//trim(st))

plugins/QmcChem/dressed_dmc.irp.f

@@ -57,7 +57,7 @@ program dressed_dmc
   enddo
 
 
-  call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,,N_int,6)
+  call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,N_int,6)
 
   call save_wavefunction
   call write_spindeterminants

plugins/Selectors_full/zmq.irp.f

@@ -114,7 +114,7 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id, energy, size_energy)
   if (N_det_selectors_read > 0) then
     N_det_selectors  = N_det_selectors_read
   endif
-  SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators
+  SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators psi_coef_generators psi_det_generators
 
 end
 

plugins/mrcepa0/EZFIO.cfg

@@ -3,3 +3,31 @@ 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
+
+[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-4
+
+[n_it_max_dressed_ci]
+type: Strictly_positive_int
+doc: Maximum number of dressed CI iterations
+interface: ezfio,provider,ocaml
+default: 10
+

plugins/mrcepa0/mrcepa0_general.irp.f

@@ -17,8 +17,8 @@ subroutine run(N_st,energy)
   
 
   
-  thresh_mrcc = 1d-7
+  thresh_mrcc = thresh_dressed_ci
-  n_it_mrcc_max = 10
+  n_it_mrcc_max = n_it_max_dressed_ci
 
   if(n_it_mrcc_max == 1) then
     do j=1,N_states_diag
@@ -48,8 +48,8 @@ subroutine run(N_st,energy)
       E_new = sum(ci_energy_dressed)
       delta_E = dabs(E_new - E_old)
       call save_wavefunction
-      call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
+      call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
-      if (iteration > n_it_mrcc_max) then
+      if (iteration >= n_it_mrcc_max) then
         exit
       endif
     enddo
@@ -184,7 +184,7 @@ subroutine run_pt2_old(N_st,energy)
   print *,  '-----' 
  
 
-!  call ezfio_set_full_ci_energy_pt2(energy+pt2)
+  call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1))
 
 end 
 
@@ -238,5 +238,7 @@ subroutine run_pt2(N_st,energy)
   print *,  'E+PT2    = ', energy+pt2 
   print *,  '-----' 
 
+  call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1))
+
 end 
 

src/Davidson/davidson_parallel.irp.f

@@ -48,8 +48,8 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
           ext = ext + popcnt(xor(sorted_i(ni), sorted_(ni,j,1)))
         end do
         if(ext <= 4) then
-          call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
           call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2) 
+          call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
           if(.not. wrotten(ii)) then
             wrotten(ii) = .true.
             idx(ii) = org_i
@@ -60,6 +60,8 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
             vt (istate,ii) += hij*dav_ut(istate,org_j)
             st (istate,ii) += s2*dav_ut(istate,org_j)
           enddo
+!          call daxpy(N_states_diag,hij,dav_ut(1,org_j),1,vt(1,org_i),1)
+!          call daxpy(N_states_diag,s2, dav_ut(1,org_j),1,st(1,org_i),1)
         endif
       enddo
     enddo
@@ -68,7 +70,7 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
   
 
   do sh=blockb,min(blocke, shortcut_(0,2))
-  do sh2=sh, shortcut_(0,2), shortcut_(0,1)*51
+  do sh2=sh, shortcut_(0,2), shortcut_(0,1)
     do i=shortcut_(sh2,2),shortcut_(sh2+1,2)-1
       ii += 1
       org_i = sort_idx_(i,2)
@@ -88,6 +90,8 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
            vt (:,ii) = 0d0
            st (:,ii) = 0d0
          end if
+!         call daxpy(N_states_diag,hij,dav_ut(1,org_j),1,vt(1,org_i),1)
+!         call daxpy(N_states_diag,s2, dav_ut(1,org_j),1,st(1,org_i),1)
          do istate=1,N_states_diag
            vt (istate,ii) += hij*dav_ut(istate,org_j)
            st (istate,ii) += s2*dav_ut(istate,org_j)
@@ -107,6 +111,7 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
       st(:,N) = st(:,i)
     end if
   end do
+
 end subroutine
 
 
@@ -138,13 +143,33 @@ subroutine davidson_collect(blockb, blocke, N, idx, vt, st , v0t, s0t)
 end subroutine
 
 
-subroutine davidson_init(zmq_to_qp_run_socket)
+subroutine davidson_init(zmq_to_qp_run_socket,n,n_st_8,ut)
   use f77_zmq
   implicit none
   
   integer(ZMQ_PTR), intent(out)  :: zmq_to_qp_run_socket
+  integer, intent(in) :: n, n_st_8
+  double precision, intent(in) :: ut(n_st_8,n)
+  integer :: i,k
+
+  
+  dav_size = n
+  touch dav_size
+
+  do i=1,n
+    do k=1,N_int
+      dav_det(k,1,i) = psi_det(k,1,i)
+      dav_det(k,2,i) = psi_det(k,2,i)
+    enddo
+  enddo
+  do i=1,n
+    do k=1,N_states_diag
+      dav_ut(k,i) = ut(k,i)
+    enddo
+  enddo
+
+  touch dav_det dav_ut
 
-  touch dav_size dav_det dav_ut
   call new_parallel_job(zmq_to_qp_run_socket,"davidson")
 end subroutine
 
@@ -260,6 +285,7 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id)
     call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
     call davidson_push_results(zmq_socket_push, blockb, blocke, N, idx, vt, st, task_id)
   end do
+
 end subroutine
 
 
@@ -340,18 +366,18 @@ end subroutine
 
 
 
-subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0)
+subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0, LDA)
   use f77_zmq
   implicit none
 
+  integer                        :: LDA
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer(ZMQ_PTR), intent(in)   :: zmq_socket_pull
   
-  double precision    ,intent(inout)    :: v0(dav_size, N_states_diag)
+  double precision    ,intent(inout) :: v0(LDA, N_states_diag)
-  double precision    ,intent(inout)    :: s0(dav_size, N_states_diag)
+  double precision    ,intent(inout) :: s0(LDA, N_states_diag)
-
-  integer                          :: more, task_id
   
+  integer                          :: more, task_id, taskn
   
   integer                        :: blockb, blocke
   integer                        :: N
@@ -381,16 +407,25 @@ subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0)
   end do
   deallocate(idx,vt,st)
 
-  call dtranspose(v0t,size(v0t,1), v0, size(v0,1), N_states_diag, dav_size)
+  integer :: i,j
-  call dtranspose(s0t,size(s0t,1), s0, size(s0,1), N_states_diag, dav_size)
+  !DIR$ IVDEP
+  do j=1,N_states_diag
+   !DIR$ IVDEP
+   do i=1,dav_size
+     v0(i,j) = v0t(j,i)
+     s0(i,j) = s0t(j,i)
+   enddo
+  enddo
+
   deallocate(v0t,s0t)
 end subroutine
 
 
-subroutine davidson_run(zmq_to_qp_run_socket , v0, s0)
+subroutine davidson_run(zmq_to_qp_run_socket , v0, s0, LDA)
   use f77_zmq
   implicit none
   
+  integer                        :: LDA
   integer(ZMQ_PTR), intent(in)   :: zmq_to_qp_run_socket
   integer(ZMQ_PTR),external      :: new_zmq_to_qp_run_socket
   integer(ZMQ_PTR)               :: zmq_collector
@@ -400,8 +435,8 @@ subroutine davidson_run(zmq_to_qp_run_socket , v0, s0)
   integer :: i
   integer, external              :: omp_get_thread_num
 
-  double precision    , intent(inout)     :: v0(dav_size, N_states_diag)
+  double precision    , intent(inout)     :: v0(LDA, N_states_diag)
-  double precision    , intent(inout)     :: s0(dav_size, N_states_diag)
+  double precision    , intent(inout)     :: s0(LDA, N_states_diag)
   
   call zmq_set_running(zmq_to_qp_run_socket)
   
@@ -412,19 +447,20 @@ subroutine davidson_run(zmq_to_qp_run_socket , v0, s0)
   
   PROVIDE nproc
   
-  !$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+2)
+  !$OMP PARALLEL NUM_THREADS(nproc+2) PRIVATE(i)
   i = omp_get_thread_num()
-    if (i==0) then
+  if (i == 0 ) then
-      call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0)
+    call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0, LDA)
     call end_zmq_to_qp_run_socket(zmq_collector)
     call end_zmq_pull_socket(zmq_socket_pull)
     call davidson_miniserver_end()
-    else if(i==1) then
+  else if (i == 1 ) then
-      call davidson_miniserver_run()
+    call davidson_miniserver_run ()
   else
     call davidson_slave_inproc(i)
   endif
   !$OMP END PARALLEL
+
   call end_parallel_job(zmq_to_qp_run_socket, 'davidson')
 end subroutine
 
@@ -505,15 +541,28 @@ end subroutine
 
 
 
-BEGIN_PROVIDER [ integer(bit_kind), dav_det, (N_int, 2, dav_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), dav_det, (N_int, 2, dav_size) ]
-END_PROVIDER
+&BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ]
-
+ use bitmasks
-
+ implicit none
-BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ]
+ BEGIN_DOC
+! Temporary arrays for parallel davidson
+!
+! Touched in davidson_miniserver_get
+ END_DOC
+ dav_det = 0_bit_kind
+ dav_ut = -huge(1.d0)
 END_PROVIDER
 
 
 BEGIN_PROVIDER [ integer, dav_size ]
+ implicit none
+ BEGIN_DOC
+! Size of the arrays for Davidson
+!
+! Touched in davidson_miniserver_get
+ END_DOC
+ dav_size = 1
 END_PROVIDER
 
 

src/Davidson/diagonalize_CI_SC2.irp.f

@@ -1,62 +0,0 @@
-BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states_diag) ]
-  implicit none
-  BEGIN_DOC
-  ! N_states_diag lowest eigenvalues of the CI matrix
-  END_DOC
-  
-  integer                        :: j
-  character*(8)                  :: st
-  call write_time(output_determinants)
-  do j=1,N_states_diag
-    CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion
-    write(st,'(I4)') j
-    call write_double(output_determinants,CI_SC2_energy(j),'Energy of state '//trim(st))
-  enddo
-
-END_PROVIDER
-
- BEGIN_PROVIDER [ double precision, threshold_convergence_SC2]
- implicit none
- BEGIN_DOC
- ! convergence of the correlation energy of SC2 iterations
- END_DOC
- threshold_convergence_SC2 = 1.d-10
-
- END_PROVIDER
-
- BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ]
-&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ]
-&BEGIN_PROVIDER [ double precision, Diag_H_elements_SC2, (N_det) ]
-  implicit none
-  BEGIN_DOC
-  ! Eigenvectors/values of the CI matrix
-  END_DOC
-  integer                        :: i,j
-  
-  do j=1,N_states_diag
-    do i=1,N_det
-      CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
-    enddo
-    CI_SC2_electronic_energy(j) = CI_electronic_energy(j) 
-  enddo
-  
-  call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
-!       size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
-        diag_H_elements_SC2,size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
-END_PROVIDER
-
-subroutine diagonalize_CI_SC2
-  implicit none
-  BEGIN_DOC
-!  Replace the coefficients of the CI states_diag by the coefficients of the 
-!  eigenstates of the CI matrix
-  END_DOC
-  integer :: i,j
-  do j=1,N_states_diag
-    do i=1,N_det
-      psi_coef(i,j) = CI_SC2_eigenvectors(i,j)
-    enddo
-  enddo
-  SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors diag_h_elements_sc2
-! SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors 
-end

src/Davidson/u0Hu0.irp.f

@@ -24,6 +24,7 @@ subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
   do i=1,N_st
     e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
   enddo
+  deallocate (H_jj, v_0)
 end
 
 
@@ -56,10 +57,9 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
   integer                        :: N_st_8
   
   integer, external              :: align_double
-  !!!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut
 
-  if(N_st /= N_states_diag) stop "H_u_0_nstates N_st /= N_states_diag"
+  N_st_8 = align_double(N_st)
-  N_st_8 = N_states_diag ! align_double(N_st)
 
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
@@ -165,7 +165,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
       v_0(i,istate) += H_jj(i) * u_0(i,istate)
     enddo
   enddo
-  !deallocate (shortcut, sort_idx, sorted, version, ut)
+  deallocate (shortcut, sort_idx, sorted, version, ut)
 end
 
 BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]
@@ -200,8 +200,8 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
   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, istate
@@ -213,38 +213,31 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
   integer(ZMQ_PTR) :: handler
   
   if(N_st /= N_states_diag .or. sze_8 < N_det) stop "assert fail in H_S2_u_0_nstates"
-  N_st_8 = N_st !! align_double(N_st)
+  N_st_8 = N_st ! align_double(N_st)
 
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)
   ASSERT (n>0)
   PROVIDE ref_bitmask_energy 
 
-  allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
+  allocate (shortcut(0:n+1,2), sort_idx(n), sorted(Nint,n), version(Nint,n))
   allocate(ut(N_st_8,n))
 
   v_0 = 0.d0
   s_0 = 0.d0
   
-  if(n /= N_det) stop "n /= N_det"
-  
   do i=1,n
     do istate=1,N_st
       ut(istate,i) =  u_0(i,istate)
     enddo
   enddo
-   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_ab_v(keys_tmp, sorted, sort_idx, shortcut(0,1), version, 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)
+   call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut(0,2), version, n, Nint)
-    
-  dav_size = n
-  touch dav_size
-  dav_det = psi_det
-  dav_ut = ut
     
   workload = 0
   blockb = shortcut(0,1)
   blocke = blockb
-  call davidson_init(handler)
+  call davidson_init(handler,n,N_st_8,ut)
   do sh=shortcut(0,1),1,-1
     workload += (shortcut(sh+1,1) - shortcut(sh,1))**2
     if(workload > max_workload) then
@@ -256,8 +249,7 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
   enddo
   
   if(blockb > 0) call davidson_add_task(handler, 1, blockb)
-  call davidson_run(handler, v_0, s_0)
+  call davidson_run(handler, v_0, s_0, size(v_0,1))
-
 
   do istate=1,N_st
     do i=1,n
@@ -265,6 +257,8 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
       s_0(i,istate) = s_0(i,istate) + s2_jj(i)* u_0(i,istate)
     enddo
   enddo
+  deallocate(shortcut, sort_idx, sorted, version)
+  deallocate(ut)
 end
 
 

tests/bats/cassd.bats

@@ -0,0 +1,17 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+
+@test "CAS_SD H2O cc-pVDZ" {
+  test_exe cas_sd_selected || skip
+  INPUT=h2o.ezfio
+  qp_edit -c $INPUT
+  ezfio set_file $INPUT
+  ezfio set perturbation do_pt2_end False
+  ezfio set determinants n_det_max 1000
+  qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" 
+  qp_run cas_sd_selected $INPUT  
+  energy="$(ezfio get cas_sd energy)"
+  eq $energy -76.22213389282479 1.E-5
+}
+

tests/bats/common.bats.sh

@@ -0,0 +1,44 @@
+#!/usr/bin/env bats
+
+# floating point number comparison
+# Compare two numbers ($1, $2) with a given precision ($3)
+# If the numbers are not equal, the exit code is 1 else it is 0
+# So we strip the "-", is the abs value of the poor
+function eq() {
+    declare -a diff
+    diff=($(awk -v d1=$1 -v d2=$2 -v n1=${1#-} -v n2=${2#-} -v p=$3 'BEGIN{ if ((n1-n2)^2 < p^2) print 0; print 1 " " (d1-d2) " " d1 " " d2 }'))
+    if [[ "${diff[0]}" == "0" ]]
+    then
+       return 0
+    else
+       echo "Test      : " ${BATS_TEST_DESCRIPTION}
+       echo "Error     : " ${diff[1]}
+       echo "Reference : " ${diff[3]}
+       echo "Computed  : " ${diff[2]}
+       exit 1
+    fi
+}
+
+
+#   ___           
+#    |  ._  o _|_ 
+#   _|_ | | |  |_ 
+#                 
+source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh
+TEST_DIR=${QP_ROOT}/tests/work/
+
+mkdir -p "${TEST_DIR}"
+
+cd "${TEST_DIR}" || exit 1
+
+function test_exe() {
+  l_EXE=$(awk "/^$1 / { print \$2 }" < "${QP_ROOT}"/data/executables)
+  l_EXE=$(echo $l_EXE | sed "s|\$QP_ROOT|$QP_ROOT|")
+  if [[ -x "$l_EXE" ]]
+  then
+    return 0
+  else
+    return 127
+  fi
+}
+

tests/bats/convert.bats

@@ -0,0 +1,27 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+
+#=== Convert
+@test "gamess convert HBO.out" {
+  cp ${QP_ROOT}/tests/input/HBO.out .
+  qp_convert_output_to_ezfio.py HBO.out
+  qp_edit -c HBO.out.ezfio 
+  ezfio set_file HBO.out.ezfio
+  qp_run SCF HBO.out.ezfio 
+  # Check energy
+  energy="$(ezfio get hartree_fock energy)"
+  eq $energy -100.0185822590964 1.e-10
+}
+
+@test "g09 convert H2O.log" {
+  cp ${QP_ROOT}/tests/input/h2o.log .
+  qp_convert_output_to_ezfio.py h2o.log
+  qp_edit -c h2o.log.ezfio 
+  ezfio set_file h2o.log.ezfio
+  qp_run SCF h2o.log.ezfio 
+  # Check energy
+  energy="$(ezfio get hartree_fock energy)"
+  eq $energy -76.0270218704265 1E-10
+}
+

tests/bats/fci.bats

@@ -0,0 +1,52 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+
+function run_FCI() {
+  thresh=5.e-5
+  test_exe full_ci || skip
+  qp_edit -c $1
+  ezfio set_file $1
+  ezfio set perturbation do_pt2_end True
+  ezfio set determinants n_det_max $2
+  ezfio set davidson threshold_davidson 1.e-10
+
+  qp_run full_ci $1
+  energy="$(ezfio get full_ci energy)"
+  eq $energy $3 $thresh
+  energy_pt2="$(ezfio get full_ci energy_pt2)"
+  eq $energy_pt2 $4 $thresh
+}
+
+function run_FCI_ZMQ() {
+  thresh=5.e-5
+  test_exe full_ci || skip
+  qp_edit -c $1
+  ezfio set_file $1
+  ezfio set perturbation do_pt2_end True
+  ezfio set determinants n_det_max $2
+  ezfio set davidson threshold_davidson 1.e-10
+
+  qp_run fci_zmq $1
+  energy="$(ezfio get full_ci energy)"
+  eq $energy $3 $thresh
+  energy_pt2="$(ezfio get full_ci energy_pt2)"
+  eq $energy_pt2 $4 $thresh
+}
+
+
+
+#=== H2O
+
+@test "qp_set_mo_class H2O cc-pVDZ" {
+  qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]"
+}
+@test "FCI H2O cc-pVDZ" {
+  run_FCI h2o.ezfio 2000  -0.761255633582109E+02 -0.761258377850042E+02
+}
+
+@test "FCI-ZMQ H2O cc-pVDZ" {
+  run_FCI_ZMQ h2o.ezfio 2000  -0.761255633582109E+02 -0.761258377850042E+02
+}
+
+

tests/bats/foboci.bats

@@ -0,0 +1,27 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+
+function run_all_1h_1p() {
+  thresh=1.e-6
+  test_exe all_1h_1p || skip
+  qp_edit -c $1
+  ezfio set_file $1
+  ezfio set determinants n_det_max $2
+  ezfio set perturbation pt2_max $3
+  ezfio set davidson threshold_davidson 1.e-10
+
+  qp_run all_1h_1p $1 | tee $1.F1h1p.out
+  energy="$(ezfio get all_singles energy)"
+  eq $energy $4 $thresh
+}
+
+
+#=== DHNO
+
+@test "all_1h_1p DHNO chipman-dzp" {
+  qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio 
+  run_all_1h_1p dhno.ezfio 10000 0.0000000001  -130.4466283766202
+}
+
+

tests/bats/hf.bats

@@ -0,0 +1,52 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+
+function run_init() {
+  cp "${QP_ROOT}/tests/input/$1" .
+  qp_create_ezfio_from_xyz $1 -o $3 $2
+  qp_edit -c $3
+}
+
+
+function run_HF() {
+  thresh=1.e-7
+  test_exe SCF || skip
+  qp_edit -c $1
+  ezfio set_file $1
+  ezfio set hartree_fock thresh_scf 1.e-11
+  qp_run SCF $1
+  energy="$(ezfio get hartree_fock energy)"
+  eq $energy $2 $thresh
+}
+
+
+
+#=== DHNO
+@test "init DHNO chipman-dzp" {
+  run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio
+}
+
+@test "SCF DHNO chipman-dzp" {
+  run_HF  dhno.ezfio  -130.4278777822   
+}
+
+#=== HBO
+@test "init HBO STO-3G" {
+  run_init HBO.xyz "-b STO-3G" hbo.ezfio
+}
+
+@test "SCF HBO STO-3G" {
+  run_HF  hbo.ezfio  -98.8251985678084 
+}
+
+
+#=== H2O
+@test "init H2O cc-pVDZ" {
+  run_init h2o.xyz "-b cc-pvdz" h2o.ezfio
+}
+
+@test "SCF H2O cc-pVDZ" {
+  run_HF  h2o.ezfio  -0.760270218692179E+02
+}
+

tests/bats/mrcepa0.bats

@@ -0,0 +1,70 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+
+#=== H2O
+@test "MRCC-lambda H2O cc-pVDZ" {
+  INPUT=h2o.ezfio
+  EXE=mrcc
+  test_exe $EXE || skip
+  qp_edit -c $INPUT  
+  ezfio set_file $INPUT
+  ezfio set determinants threshold_generators 1.
+  ezfio set determinants threshold_selectors  1.
+  ezfio set determinants read_wf True
+  ezfio set mrcepa0 lambda_type 1
+  ezfio set mrcepa0 n_it_max_dressed_ci 3
+  qp_run $EXE $INPUT  
+  energy="$(ezfio get mrcepa0 energy)"
+  eq $energy -76.22903276183061 1.e-4
+}
+
+@test "MRCC H2O cc-pVDZ" {
+  INPUT=h2o.ezfio
+  EXE=mrcc
+  test_exe $EXE || skip
+  qp_edit -c $INPUT  
+  ezfio set_file $INPUT
+  ezfio set determinants threshold_generators 1.
+  ezfio set determinants threshold_selectors  1.
+  ezfio set determinants read_wf True
+  ezfio set determinants read_wf True
+  ezfio set mrcepa0 lambda_type 0
+  ezfio set mrcepa0 n_it_max_dressed_ci 3
+  qp_run $EXE $INPUT  
+  energy="$(ezfio get mrcepa0 energy)"
+  eq $energy -76.22899302846875 1.e-4
+}
+
+@test "MRSC2 H2O cc-pVDZ" {
+  INPUT=h2o.ezfio
+  EXE=mrsc2
+  test_exe $EXE || skip
+  qp_edit -c $INPUT  
+  ezfio set_file $INPUT
+  ezfio set determinants threshold_generators 1.
+  ezfio set determinants threshold_selectors  1.
+  ezfio set determinants read_wf True
+  ezfio set mrcepa0 lambda_type 0
+  ezfio set mrcepa0 n_it_max_dressed_ci 3
+  qp_run $EXE $INPUT  
+  energy="$(ezfio get mrcepa0 energy)"
+  eq $energy -76.22647345292708 1.e-4
+}
+
+@test "MRCEPA0 H2O cc-pVDZ" {
+  INPUT=h2o.ezfio
+  EXE=mrcepa0
+  test_exe $EXE || skip
+  qp_edit -c $INPUT  
+  ezfio set_file $INPUT
+  ezfio set determinants threshold_generators 1.
+  ezfio set determinants threshold_selectors  1.
+  ezfio set determinants read_wf True
+  ezfio set mrcepa0 lambda_type 0
+  ezfio set mrcepa0 n_it_max_dressed_ci 3
+  qp_run $EXE $INPUT  
+  energy="$(ezfio get mrcepa0 energy)"
+  eq $energy -76.23199784430074 1.e-4
+}
+

tests/bats/pseudo.bats

@@ -0,0 +1,53 @@
+#!/usr/bin/env bats
+
+source $QP_ROOT/tests/bats/common.bats.sh
+
+function run_init() {
+  cp "${QP_ROOT}/tests/input/$1" .
+  qp_create_ezfio_from_xyz $1 -o $3 $2
+  qp_edit -c $3
+}
+
+
+function run_HF() {
+  thresh=1.e-7
+  test_exe SCF || skip
+  qp_edit -c $1
+  ezfio set_file $1
+  ezfio set hartree_fock thresh_scf 1.e-11
+  qp_run SCF $1
+  energy="$(ezfio get hartree_fock energy)"
+  eq $energy $2 $thresh
+}
+
+
+function run_FCI_ZMQ() {
+  thresh=5.e-5
+  test_exe full_ci || skip
+  qp_edit -c $1
+  ezfio set_file $1
+  ezfio set perturbation do_pt2_end True
+  ezfio set determinants n_det_max $2
+  ezfio set davidson threshold_davidson 1.e-10
+
+  qp_run fci_zmq $1
+  energy="$(ezfio get full_ci energy)"
+  eq $energy $3 $thresh
+  energy_pt2="$(ezfio get full_ci energy_pt2)"
+  eq $energy_pt2 $4 $thresh
+}
+
+#=== H2O Pseudo
+@test "init H2O VDZ pseudo" {
+  run_init h2o.xyz "-p bfd -b vdz-bfd" h2o_pseudo.ezfio
+}
+
+@test "SCF H2O VDZ pseudo" {
+  run_HF  h2o_pseudo.ezfio  -16.9483703905461
+}
+
+@test "FCI H2O VDZ pseudo" {
+  qp_set_mo_class h2o_pseudo.ezfio -core "[1]" -act "[2-12]" -del "[13-23]"
+  run_FCI_ZMQ h2o_pseudo.ezfio 2000    -0.170399597228904E+02 -0.170400168816800E+02
+}
+

tests/bats/qp.bats

@@ -1,205 +0,0 @@
-#!/usr/bin/env bats
-
-#                
-#   |\/| o  _  _ 
-#   |  | | _> (_ 
-#                
-# floating point number comparison
-# Compare two numbers ($1, $2) with a given precision ($3)
-# If the numbers are not equal, the exit code is 1 else it is 0
-# So we strip the "-", is the abs value of the poor
-function eq() {
-    declare -a diff
-    diff=($(awk -v d1=$1 -v d2=$2 -v n1=${1#-} -v n2=${2#-} -v p=$3 'BEGIN{ if ((n1-n2)^2 < p^2) print 0; print 1 " " (d1-d2) " " d1 " " d2 }'))
-    if [[ "${diff[0]}" == "0" ]]
-    then
-       return 0
-    else
-       echo "Test      : " ${BATS_TEST_DESCRIPTION}
-       echo "Error     : " ${diff[1]}
-       echo "Reference : " ${diff[3]}
-       echo "Computed  : " ${diff[2]}
-       exit 1
-    fi
-}
-
-
-#   ___           
-#    |  ._  o _|_ 
-#   _|_ | | |  |_ 
-#                 
-source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh
-TEST_DIR=${QP_ROOT}/tests/work/
-
-mkdir -p "${TEST_DIR}"
-
-cd "${TEST_DIR}" || exit 1
-
-function run_init() {
-  cp "${QP_ROOT}/tests/input/$1" .
-  qp_create_ezfio_from_xyz $1 -o $3 $2 
-  qp_edit -c $3
-}
-
-function test_exe() {
-  EXE=$(awk "/^$1 / { print \$2 }" < "${QP_ROOT}"/data/executables)
-  EXE=$(echo $EXE | sed "s|\$QP_ROOT|$QP_ROOT|")
-  if [[ -x "$EXE" ]]
-  then
-    return 0
-  else
-    return 127
-  fi
-}
-
-function run_HF() {
-  thresh=1.e-7
-  test_exe SCF || skip
-  ezfio set_file $1
-  ezfio set hartree_fock thresh_scf 1.e-11
-  qp_run SCF $1 
-  energy="$(ezfio get hartree_fock energy)"
-  eq $energy $2 $thresh
-}
-
-function run_FCI() {
-  thresh=5.e-5
-  test_exe full_ci || skip
-  ezfio set_file $1
-  ezfio set perturbation do_pt2_end True
-  ezfio set determinants n_det_max $2
-  ezfio set davidson threshold_davidson 1.e-10
-
-  qp_run full_ci $1 
-  energy="$(ezfio get full_ci energy)"
-  eq $energy $3 $thresh
-  energy_pt2="$(ezfio get full_ci energy_pt2)"
-  eq $energy_pt2 $4 $thresh
-}
-
-function run_all_1h_1p() {
-  thresh=1.e-6
-  test_exe all_1h_1p || skip
-  ezfio set_file $1
-  ezfio set determinants n_det_max $2
-  ezfio set perturbation pt2_max $3
-  ezfio set davidson threshold_davidson 1.e-10
-
-  qp_run all_1h_1p $1 | tee $1.F1h1p.out
-  energy="$(ezfio get all_singles energy)"
-  eq $energy $4 $thresh
-}
-
-#   ___             
-#    |  _   _ _|_   
-#    | (/_ _>  |_   
-#                   
-
-
-#=== DHNO
-@test "init DHNO chipman-dzp" {
-  run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio
-}
-
-@test "SCF DHNO chipman-dzp" {
-  run_HF  dhno.ezfio  -130.4278777822   
-}
-
-@test "all_1h_1p DHNO chipman-dzp" {
-  qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio 
-  run_all_1h_1p dhno.ezfio 10000 0.0000000001  -130.4466283766202
-}
-
-#=== HBO
-@test "init HBO STO-3G" {
-  run_init HBO.xyz "-b STO-3G" hbo.ezfio
-}
-
-@test "SCF HBO STO-3G" {
-  run_HF  hbo.ezfio  -98.8251985678084 
-}
-
-
-#=== H2O
-@test "init H2O cc-pVDZ" {
-  run_init h2o.xyz "-b cc-pvdz" h2o.ezfio
-}
-
-@test "SCF H2O cc-pVDZ" {
-  run_HF  h2o.ezfio  -0.760270218692179E+02
-}
-
-@test "FCI H2O cc-pVDZ" {
-  qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]"
-  run_FCI h2o.ezfio 2000  -0.761255633582109E+02 -0.761258377850042E+02
-}
-
-@test "CAS_SD H2O cc-pVDZ" {
-  test_exe cas_sd_selected || skip
-  INPUT=h2o.ezfio
-  ezfio set_file $INPUT
-  ezfio set perturbation do_pt2_end False
-  ezfio set determinants n_det_max 1000
-  qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]" 
-  qp_run cas_sd_selected $INPUT  
-  energy="$(ezfio get cas_sd energy)"
-  eq $energy -0.762219854008117E+02 1.E-5
-}
-
-@test "MRCC H2O cc-pVDZ" {
-  test_exe mrcc_cassd || skip
-  INPUT=h2o.ezfio
-  ezfio set_file $INPUT
-  ezfio set determinants threshold_generators 1.
-  ezfio set determinants threshold_selectors  1.
-  ezfio set determinants read_wf True
-  qp_run mrcc_cassd $INPUT  
-  energy="$(ezfio get mrcc_cassd energy)"
-  eq $energy -76.2288648023833  1.e-4
-  
-}
-
-
-#=== H2O Pseudo
-@test "init H2O VDZ pseudo" {
-  run_init h2o.xyz "-p bfd -b vdz-bfd" h2o_pseudo.ezfio
-}
-
-@test "SCF H2O VDZ pseudo" {
-  run_HF  h2o_pseudo.ezfio  -16.9483703905461
-}
-
-@test "FCI H2O VDZ pseudo" {
-  qp_set_mo_class h2o_pseudo.ezfio -core "[1]" -act "[2-12]" -del "[13-23]"
-  run_FCI h2o_pseudo.ezfio 2000    -0.170399597228904E+02 -0.170400168816800E+02
-}
-
-#=== Convert
-@test "gamess convert HBO.out" {
-  cp ${QP_ROOT}/tests/input/HBO.out .
-  qp_convert_output_to_ezfio.py HBO.out
-  ezfio set_file HBO.out.ezfio
-  qp_run SCF HBO.out.ezfio 
-  # Check energy
-  energy="$(ezfio get hartree_fock energy)"
-  eq $energy -100.0185822590964 1.e-10
-}
-
-@test "g09 convert H2O.log" {
-  cp ${QP_ROOT}/tests/input/h2o.log .
-  qp_convert_output_to_ezfio.py h2o.log
-  ezfio set_file h2o.log.ezfio
-  qp_run SCF h2o.log.ezfio 
-  # Check energy
-  energy="$(ezfio get hartree_fock energy)"
-  eq $energy -76.0270218704265 1E-10
-}
-
-
-# TODO N_int = 1,2,3,4,5
-# TODO mod(64) MOs
-# TODO All G2 SCF energies
-# TODO Long and short tests
-# TODO MP2
-# TODO CISD_selected
-

tests/bats_to_sh.py

@@ -1,6 +1,8 @@
 #!/usr/bin/env python
 
-with open('bats/qp.bats','r') as f:
+import sys
+
+with open(sys.argv[1],'r') as f:
   raw_data = f.read()
 
 output = []

tests/run_tests.sh

@@ -1,18 +1,39 @@
 #!/bin/bash
 
+LIST="
+
+convert.bats
+hf.bats
+foboci.bats
+pseudo.bats
+fci.bats
+cassd.bats
+mrcepa0.bats
+
+"
+
+
 export QP_PREFIX="timeout -s 9 300"
 export QP_TASK_DEBUG=1
 
-BATS_FILE=bats/qp.bats 
-
 rm -rf work output
 
-if [[ "$1" == "-v" ]]
+
-then
+for BATS_FILE in $LIST
+do
+  echo
+  echo "-~-~-~-~-~-~"
+  echo
+  echo "Running tests for ${BATS_FILE%.bats}"
+  echo
+  BATS_FILE=bats/$BATS_FILE
+  if [[ "$1" == "-v" ]]
+  then
     echo "Verbose mode"
     ./bats_to_sh.py $BATS_FILE | bash
-else
+  else
     bats $BATS_FILE
-fi
+  fi
+done