Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable

external/ezfio

@@ -1 +1 @@
-Subproject commit d02132ea79217c16fd24242e8f8b8a6c3ff68091
+Subproject commit dba01c4fe0ff7b84c5ecfb1c7c77ec68781311b3

scripts/qp_exc_energy.py

@@ -157,11 +157,15 @@ A = np.array( [ [ data[-1][1], 1. ],
 B = np.array( [ [ data[-1][0] ],
                 [ data[-2][0] ] ] )
 E0 = np.linalg.solve(A,B)[1]
+E0 = E0[0]
+
 A = np.array( [ [ data[-1][4], 1. ],
                 [ data[-2][4], 1. ] ] )
 B = np.array( [ [ data[-1][3] ],
                 [ data[-2][3] ] ] )
 E1 = np.linalg.solve(A,B)[1]
+E1 = E1[0]
+
 average_2 = (E1-E0)*to_eV
 
 A = np.array( [ [ data[-1][1], 1. ],
@@ -170,14 +174,18 @@ A = np.array( [ [ data[-1][1], 1. ],
 B = np.array( [ [ data[-1][0] ],
                 [ data[-2][0] ],
                 [ data[-3][0] ] ] )
-E0 = np.linalg.lstsq(A,B,rcond=None)[0][1]
+E0 = np.linalg.lstsq(A,B,rcond=None)[0]
+E0 = E0[0][0]
+
 A = np.array( [ [ data[-1][4], 1. ],
                 [ data[-2][4], 1. ],
                 [ data[-3][4], 1. ] ] )
 B = np.array( [ [ data[-1][3] ],
                 [ data[-2][3] ],
                 [ data[-3][3] ] ] )
-E1 = np.linalg.lstsq(A,B,rcond=None)[0][1]
+E1 = np.linalg.lstsq(A,B,rcond=None)[0]
+E1 = E1[0][0]
+
 average_3 = (E1-E0)*to_eV
 
 exc = ((data[-1][3] + data[-1][4]) - (data[-1][0] + data[-1][1])) * to_eV

src/casscf_cipsi/chol_bielec.irp.f

@@ -191,10 +191,15 @@ double precision function bielec_PQxx_no(i_mo, j_mo, i_ca, j_ca)
  END_DOC
  integer, intent(in) :: i_ca, j_ca, i_mo, j_mo
  integer :: ii_ca, jj_ca
- double precision :: bielec_no_basis
  ii_ca = list_core_inact_act(i_ca)
  jj_ca = list_core_inact_act(j_ca)
- bielec_PQxx_no = bielec_no_basis(i_mo,j_mo,ii_ca,jj_ca)
+! double precision :: bielec_no_basis
+! bielec_PQxx_no = bielec_no_basis(i_mo,j_mo,ii_ca,jj_ca)
+ integer :: i
+ bielec_PQxx_no = 0.d0
+ do i = 1, cholesky_mo_num
+  bielec_PQxx_no = bielec_PQxx_no + cholesky_no_total_transp(i,i_mo, j_mo) * cholesky_no_total_transp(i,ii_ca,jj_ca)
+ enddo
 end
 
 double precision function bielec_PxxQ_no(i_mo, j_ca, i_ca, j_mo)
@@ -206,10 +211,15 @@ double precision function bielec_PxxQ_no(i_mo, j_ca, i_ca, j_mo)
   END_DOC
  integer, intent(in) :: i_ca, j_ca, i_mo, j_mo
  integer :: ii_ca, jj_ca
- double precision :: bielec_no_basis
  ii_ca = list_core_inact_act(i_ca)
  jj_ca = list_core_inact_act(j_ca)
- bielec_PxxQ_no = bielec_no_basis(i_mo, jj_ca, ii_ca, j_mo)
+ double precision :: bielec_no_basis
+! bielec_PxxQ_no = bielec_no_basis(i_mo, jj_ca, ii_ca, j_mo)
+ integer :: i
+ bielec_PxxQ_no = 0.d0
+ do i = 1, cholesky_mo_num
+  bielec_PxxQ_no = bielec_PxxQ_no + cholesky_no_total_transp(i,i_mo, jj_ca) * cholesky_no_total_transp(i,ii_ca,j_mo)
+ enddo
 
 end
 

src/ccsd/ccsd_space_orb.irp.f

@@ -1,9 +1,10 @@
-! Code
-
 program ccsd
 
   implicit none
 
+  BEGIN_DOC
+  ! Closed-shell CCSD
+  END_DOC
   read_wf = .True.
   touch read_wf
 

src/ccsd/ccsd_space_orb_sub.irp.f

@@ -26,7 +26,6 @@ subroutine run_ccsd_space_orb
 
   double precision, allocatable :: all_err(:,:), all_t(:,:)
   integer, allocatable          :: list_occ(:), list_vir(:)
-  integer(bit_kind)             :: det(N_int,2)
   integer                       :: nO, nV, nOa, nVa
 
   call set_multiple_levels_omp(.False.)
@@ -38,9 +37,8 @@ subroutine run_ccsd_space_orb
     PROVIDE all_mo_integrals
   endif
 
-  det = psi_det(:,:,cc_ref)
   print*,'Reference determinant:'
-  call print_det(det,N_int)
+  call print_det(psi_det(1,1,cc_ref),N_int)
 
   nOa = cc_nOa
   nVa = cc_nVa
@@ -57,10 +55,6 @@ subroutine run_ccsd_space_orb
   allocate(list_occ(nO),list_vir(nV))
   list_occ = cc_list_occ
   list_vir = cc_list_vir
-  ! Debug
-  !call extract_list_orb_space(det,nO,nV,list_occ,list_vir)
-  !print*,'occ',list_occ
-  !print*,'vir',list_vir
 
   ! GPU arrays
   call gpu_allocate(d_cc_space_f_oo, nO, nO)
@@ -183,10 +177,12 @@ subroutine run_ccsd_space_orb
   call guess_t2(nO,nV,cc_space_f_o,cc_space_f_v,cc_space_v_oovv,h_t2)
   call gpu_upload(h_t2, t2)
 
+  deallocate(h_t1, h_t2)
 
-  call update_tau_space(nO,nV,h_t1,t1,t2,tau)
+  call update_tau_space(nO,nV,t1%f,t1,t2,tau)
   call update_tau_x_space(nO,nV,tau,tau_x)
-  call det_energy(det,uncorr_energy)
+
+  call det_energy(psi_det(1,1,cc_ref),uncorr_energy)
   print*,'Det energy', uncorr_energy
 
   call ccsd_energy_space_x(nO,nV,d_cc_space_v_oovv,d_cc_space_f_vo,tau_x,t1,energy)
@@ -316,7 +312,6 @@ subroutine run_ccsd_space_orb
 
   call save_energy(uncorr_energy + energy, e_t)
 
-  deallocate(h_t1, h_t2)
   if (do_mo_cholesky) then
     call gpu_deallocate(d_cc_space_v_oo_chol)
     call gpu_deallocate(d_cc_space_v_ov_chol)

src/ccsd/ccsd_spin_orb.irp.f

@@ -1,5 +1,3 @@
-! Prog
-
 program ccsd
 
   implicit none

src/cipsi_utils/pt2_stoch_routines.irp.f

@@ -212,6 +212,7 @@ subroutine ZMQ_pt2(E, pt2_data, pt2_data_err, relative_error, N_in)
           ipos += 1
         endif
       enddo
+      call write_int(6,pt2_stoch_istate,'State')
       call write_int(6,sum(pt2_F),'Number of tasks')
       call write_int(6,ipos,'Number of fragmented tasks')
 
@@ -530,7 +531,7 @@ subroutine pt2_collector(zmq_socket_pull, E, relative_error, pt2_data, pt2_data_
         avg  = E0 + pt2_data_S(t) % pt2(pt2_stoch_istate) / dble(c)
         avg2 = v0 + pt2_data_S(t) % variance(pt2_stoch_istate) / dble(c)
         avg3(:) = n0(:) + pt2_data_S(t) % overlap(:,pt2_stoch_istate) / dble(c)
-        if ((avg /= 0.d0) .or. (n == N_det_generators) ) then
+        if (((c>=10).and.(avg /= 0.d0)) .or. (n == N_det_generators) ) then
           do_exit = .true.
         endif
         if (qp_stop()) then

src/davidson/u0_hs2_u0.irp.f

@@ -158,7 +158,7 @@ subroutine H_S2_u_0_nstates_openmp_work(v_t,s_t,u_t,N_st,sze,istart,iend,ishift,
   double precision, intent(out)  :: v_t(N_st,sze), s_t(N_st,sze)
 
 
-  PROVIDE ref_bitmask_energy N_int
+  PROVIDE ref_bitmask_energy N_int all_mo_integrals
 
   select case (N_int)
     case (1)
@@ -291,7 +291,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
   ASSERT (istart > 0)
   ASSERT (istep  > 0)
 
-  !$OMP DO SCHEDULE(guided,64)
+  !$OMP DO SCHEDULE(dynamic,64)
   do k_a=istart+ishift,iend,istep
 
     krow = psi_bilinear_matrix_rows(k_a)
@@ -469,7 +469,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
   enddo
   !$OMP END DO
 
-  !$OMP DO SCHEDULE(guided,64)
+  !$OMP DO SCHEDULE(dynamic,64)
   do k_a=istart+ishift,iend,istep
 
 

src/mo_basis/utils.irp.f

@@ -228,7 +228,7 @@ subroutine mo_as_svd_vectors_of_mo_matrix_eig(matrix,lda,m,n,eig,label)
   call dgemm('N','N',ao_num,m,m,1.d0,mo_coef_new,size(mo_coef_new,1),U,size(U,1),0.d0,mo_coef,size(mo_coef,1))
 
   do i=1,m
-    if (eig(i) > 1.d-20) then
+    if (D(i) > 1.d-20) then
       eig(i) = D(i)
     else
       eig(i) = 0.d0

src/mo_two_e_ints/map_integrals.irp.f

@@ -81,11 +81,15 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_s
  integer(key_kind)              :: idx
  real(integral_kind)            :: integral
  FREE ao_integrals_cache
+
  if (do_mo_cholesky) then
 
    call set_multiple_levels_omp(.False.)
-   !$OMP PARALLEL DO PRIVATE (k,l,ii)
+
+
+   !$OMP PARALLEL DO PRIVATE(k,l,ii) SCHEDULE(dynamic)
    do l=mo_integrals_cache_min,mo_integrals_cache_max
+     print *, l
      do k=mo_integrals_cache_min,mo_integrals_cache_max
          ii = int(l-mo_integrals_cache_min,8)
          ii = ior( shiftl(ii,mo_integrals_cache_shift), int(k-mo_integrals_cache_min,8))
@@ -101,7 +105,7 @@ BEGIN_PROVIDER [ double precision, mo_integrals_cache, (0_8:mo_integrals_cache_s
    !$OMP END PARALLEL DO
 
  else
-   !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral)
+   !$OMP PARALLEL DO PRIVATE (i,j,k,l,idx,ii,integral) SCHEDULE(dynamic)
    do l=mo_integrals_cache_min,mo_integrals_cache_max
      do k=mo_integrals_cache_min,mo_integrals_cache_max
        do j=mo_integrals_cache_min,mo_integrals_cache_max

src/mrci/EZFIO.cfg

@@ -0,0 +1,24 @@
+[energy]
+type: double precision
+doc: Calculated Selected CASSD energy
+interface: ezfio
+size: (determinants.n_states)
+
+[energy_pt2]
+type: double precision
+doc: Calculated CASSD energy + PT2
+interface: ezfio
+size: (determinants.n_states)
+
+
+[do_ddci]
+type: logical
+doc: If true, remove purely inactive double excitations
+interface: ezfio,provider,ocaml
+default: False
+
+[do_only_1h1p]
+type: logical
+doc: If true, do only one hole/one particle excitations
+interface: ezfio,provider,ocaml
+default: False

src/mrci/NEED

@@ -0,0 +1,4 @@
+cipsi
+generators_cas
+selectors_full
+davidson_undressed

src/mrci/README.rst

@@ -0,0 +1,17 @@
+====
+mrci
+====
+
+
+|CIPSI| algorithm in the multi-reference CI space (CAS + Singles and Doubles).
+
+
+This module is the same as the :ref:`fci` module, except for the choice of the
+generator and selector determinants.
+
+The inactive, active and virtual |MOs| will need to be set with the
+:ref:`qp_set_mo_class` program.
+
+.. seealso::
+
+    The documentation of the :ref:`fci` module.

src/mrci/class.irp.f

@@ -0,0 +1,8 @@
+BEGIN_PROVIDER [ logical, do_only_cas  ]
+ implicit none
+ BEGIN_DOC
+ ! In the CAS+SD case, always false
+ END_DOC
+ do_only_cas  = .False.
+END_PROVIDER
+

src/mrci/mrci.irp.f

@@ -0,0 +1,64 @@
+program mrci
+  implicit none
+  BEGIN_DOC
+! Selected CAS+Singles and Doubles with stochastic selection
+! and PT2.
+!
+! This program performs a |CIPSI|-like selected |CI| using a
+! stochastic scheme for both the selection of the important Slater
+! determinants and the computation of the |PT2| correction. This
+! |CIPSI|-like algorithm will be performed for the lowest states of
+! the variational space (see :option:`determinants n_states`). The
+! program will stop when reaching at least one the two following
+! conditions:
+!
+! * number of Slater determinants > :option:`determinants n_det_max`
+! * |PT2| < :option:`perturbation pt2_max`
+!
+! The following other options can be of interest:
+!
+! :option:`determinants read_wf`
+!   When set to |false|, the program starts with a ROHF-like Slater
+!   determinant as a guess wave function. When set to |true|, the
+!   program starts with the wave function(s) stored in the |EZFIO|
+!   directory as guess wave function(s).
+!
+! :option:`determinants s2_eig`
+!   When set to |true|, the selection will systematically add all the
+!   necessary Slater determinants in order to have a pure spin wave
+!   function with an |S^2| value corresponding to
+!   :option:`determinants expected_s2`.
+!
+! For excited states calculations, it is recommended to start with
+! :ref:`.cis.` or :ref:`.cisd.` guess wave functions, eventually in
+! a restricted set of |MOs|, and to set :option:`determinants s2_eig`
+! to |true|.
+!
+  END_DOC
+
+  PROVIDE all_mo_integrals
+  if (.not.is_zmq_slave) then
+    PROVIDE psi_det psi_coef
+
+    write(json_unit,json_array_open_fmt) 'fci'
+
+    double precision, allocatable :: Ev(:),PT2(:)
+    allocate(Ev(N_states), PT2(N_states))
+    if (do_pt2) then
+      call run_stochastic_cipsi(Ev,PT2)
+    else
+      call run_cipsi
+    endif
+
+    write(json_unit,json_dict_uopen_fmt)
+    write(json_unit,json_dict_close_fmtx)
+    write(json_unit,json_array_close_fmtx)
+    call json_close
+
+  else
+    PROVIDE pt2_min_parallel_tasks
+
+    call run_slave_cipsi
+
+  endif
+end

src/mrci/save_energy.irp.f

@@ -0,0 +1,10 @@
+subroutine save_energy(E,pt2)
+  implicit none
+  BEGIN_DOC
+! Saves the energy in |EZFIO|.
+  END_DOC
+  double precision, intent(in) :: E(N_states), pt2(N_states)
+  call ezfio_set_mrci_energy(E(1:N_states))
+  call ezfio_set_mrci_energy_pt2(E(1:N_states)+pt2(1:N_states))
+end
+