Merge branch 'master' of https://github.com/garniron/quantum_package into garniron-master

Conflicts:
	ocaml/.gitignore
	plugins/CAS_SD/README.rst
	plugins/Full_CI/README.rst
	plugins/MRCC_Utils/README.rst
	plugins/Perturbation/README.rst
	src/Determinants/README.rst
	src/Ezfio_files/README.rst

config/gfortran.cfg

@@ -10,7 +10,7 @@
 #
 #
 [COMMON]
-FC           : gfortran -ffree-line-length-none -I .
+FC           : gfortran -ffree-line-length-none -I . -mavx
 LAPACK_LIB   : -llapack -lblas
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=32
@@ -35,7 +35,7 @@ OPENMP  : 1          ; Append OpenMP flags
 #                           -ffast-math and the Fortran-specific
 #                           -fno-protect-parens and -fstack-arrays.
 [OPT]
-FCFLAGS : -Ofast 
+FCFLAGS : -Ofast -g
 
 # Profiling flags
 #################

ocaml/.gitignore

@@ -3,45 +3,44 @@ ezfio.ml
 Qptypes.ml
 qptypes_generator.byte
 _build
-qp_basis_clean.native
 qp_create_ezfio_from_xyz.native
-qp_edit.native
-qp_print.native
-qp_run.native
 qp_set_ddci.native
-qp_set_mo_class.native
+qp_print.native
+qp_edit.native
+qp_set_mo_class.native
+qp_basis_clean.native
+qp_run.native
 qp_edit.native
-test_atom.byte
-test_basis.byte
-test_bitlist.byte
-test_determinants.byte
-test_elements.byte
-test_excitation.byte
-test_gto.byte
 test_mo_label.byte
-test_molecule.byte
 test_point3d.byte
-test_atom
+test_gto.byte
-test_basis
+test_excitation.byte
-test_bitlist
+test_determinants.byte
-test_determinants
+test_basis.byte
-test_elements
+test_molecule.byte
-test_excitation
+test_elements.byte
-test_gto
+test_bitlist.byte
+test_atom.byte
 test_mo_label
-test_molecule
 test_point3d
-qp_basis_clean
+test_gto
+test_excitation
+test_determinants
+test_basis
+test_molecule
+test_elements
+test_bitlist
+test_atom
 qp_create_ezfio_from_xyz
-qp_edit
-qp_print
-qp_run
 qp_set_ddci
+qp_print
+qp_edit
 qp_set_mo_class
+qp_basis_clean
+qp_run
 Input_determinants.ml
-Input_properties.ml
-Input_hartree_fock.ml
 Input_integrals_bielec.ml
+Input_pseudo.ml
 Input_perturbation.ml
 Input_pseudo.ml
 qp_edit.ml

plugins/Full_CI/full_ci.irp.f

@@ -11,6 +11,7 @@ program full_ci
   
   pt2 = 1.d0
   diag_algorithm = "Lapack"
+  
   if (N_det > N_det_max) then
     call diagonalize_CI
     call save_wavefunction

plugins/MRCC_Utils/H_apply.irp.f

@@ -9,7 +9,7 @@ s.data["declarations"] += """
     double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
     double precision, intent(in) :: delta_ii_(Ndet_ref,*)
 """
-s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc)"
+s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)"
 s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
 s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref"
 s.data["decls_main"] += """

plugins/MRCC_Utils/mrcc_dress.irp.f

@@ -12,12 +12,61 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ]
 
 END_PROVIDER
 
-subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc)
+
+subroutine create_minilist(key_mask, fullList, miniList, idx_miniList, N_fullList, N_miniList, Nint)
+  use bitmasks
+  implicit none
+  
+  integer(bit_kind), intent(in)            :: fullList(Nint, 2, N_fullList)
+  integer, intent(in)                      :: N_fullList
+  integer(bit_kind),intent(out)            :: miniList(Nint, 2, N_fullList)
+  integer,intent(out)                      :: idx_miniList(N_fullList), N_miniList
+  integer, intent(in)                      :: Nint
+  integer(bit_kind)                        :: key_mask(Nint, 2)
+  integer                                  :: ni, i, n_a, n_b, e_a, e_b
+  
+  
+  n_a = 0
+  n_b = 0
+  do ni=1,nint
+    n_a = n_a + popcnt(key_mask(ni,1))
+    n_b = n_b + popcnt(key_mask(ni,2))
+  end do
+  
+  if(n_a == 0) then
+    N_miniList = N_fullList
+    miniList(:,:,:) = fullList(:,:,:)
+    do i=1,N_fullList
+      idx_miniList(i) = i
+    end do
+    return
+  end if
+  
+  N_miniList = 0
+  
+  do i=1,N_fullList
+    e_a = n_a
+    e_b = n_b
+    do ni=1,nint
+      e_a -= popcnt(iand(fullList(ni, 1, i), key_mask(ni, 1)))
+      e_b -= popcnt(iand(fullList(ni, 2, i), key_mask(ni, 2)))
+    end do
+    
+    if(e_a + e_b <= 2) then
+      N_miniList = N_miniList + 1
+      miniList(:,:,N_miniList) = fullList(:,:,i)
+      idx_miniList(N_miniList) = i
+    end if
+  end do
+end subroutine
+
+
+subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask)
  use bitmasks
  implicit none
 
   integer, intent(in)            :: i_generator,n_selected, Nint, iproc
-  integer, intent(in) :: Ndet_ref, Ndet_non_ref
+  integer, intent(in)            :: Ndet_ref, Ndet_non_ref
   double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*)
   double precision, intent(inout) :: delta_ii_(Ndet_ref,*)
 
@@ -40,7 +89,12 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
   integer(bit_kind)              :: tmp_det(Nint,2)
   integer                        :: iint, ipos
   integer                        :: i_state, k_sd, l_sd, i_I, i_alpha
-
+  
+  integer(bit_kind)              :: miniList(Nint, 2, N_det_non_ref), key_mask(Nint, 2)
+  integer                        :: idx_miniList(N_det_non_ref), N_miniList
+  
+  
+  
   call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq)
 
   allocate (dIa_hla(N_states,Ndet_non_ref))
@@ -48,9 +102,20 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n
   ! |I>
 
   ! |alpha>
+  
+   if(N_tq > 0) then
+     call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint)
+   end if
+  
+  
   do i_alpha=1,N_tq
-    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
+!    call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha)
-
+    call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha)
+    
+    do j=1,idx_alpha(0)
+      idx_alpha(j) = idx_miniList(idx_alpha(j))
+    end do
+     
     ! |I>
     do i_I=1,N_det_ref
        ! Find triples and quadruple grand parents
@@ -233,7 +298,8 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq
     ! Select determinants that are triple or quadruple excitations
     ! from the ref
     good = .True.
-    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx)
+    call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) 
+    !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector
     do k=1,idx(0)
       if (degree(k) < 3) then
         good = .False.

src/Determinants/H_apply.template.f

@@ -1,4 +1,116 @@
+
+
 subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters )
+  
+  integer(bit_kind), intent(in)         :: key_in(N_int, 2), hole_1(N_int, 2), hole_2(N_int, 2)
+  integer(bit_kind), intent(in)         :: particl_1(N_int, 2), particl_2(N_int, 2)
+  integer(bit_kind)                     :: p1_mask(N_int, 2), p2_mask(N_int, 2)
+  integer,intent(in)                    :: i_generator,iproc_in
+  integer(bit_kind)                     :: status(N_int*bit_kind_size, 2)
+  integer                               :: highest, p1,p2,sp,ni,i,mi
+  $declarations
+  
+  
+  highest = 0
+  status(:,:) = 0
+  do sp=1,2
+    do ni=1,N_int
+      do i=1,bit_kind_size
+        if(iand(1,ishft(key_in(ni, sp), -(i-1))) == 0) then
+          cycle
+        end if
+        mi = (ni-1)*bit_kind_size+i
+        status(mi, sp) = iand(1,ishft(hole_1(ni, sp), -(i-1)))
+        status(mi, sp) = status(mi, sp) + 2*iand(1,ishft(hole_2(ni, sp), -(i-1)))
+        if(status(mi, sp) /= 0 .and. mi > highest) then
+          highest = mi
+        end if
+      end do
+    end do
+  end do
+  
+  do sp=1,2
+    do p1=1,highest
+      if(status(p1, sp) == 0) then
+        cycle
+      end if
+      do p2=1,highest
+        if(status(p2, sp) == 0) then
+          cycle
+        end if
+        if((status(p1, sp) == 1 .and. status(p2, sp) > 1) .or. &
+            (status(p1, sp) == 2 .and. status(p2, sp) == 3) .or. &
+            (status(p1, sp) == 3 .and. status(p2, sp) == 3 .and. p2 > p1)) then
+          call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, i_generator, iproc_in $parameters )
+        end if
+      end do
+    end do
+  end do
+  do p1=1,highest
+    if(status(p1, 1) == 0) then
+      cycle
+    end if
+    do p2=1,highest
+      if(status(p2, 2) == 0) then
+        cycle
+      end if
+      if((status(p1, 1) == 3) .or. &
+          (status(p1, 1) == 1 .and. status(p2, 2) >= 2) .or. &
+          (status(p1, 1) == 2 .and. status(p2, 2) /= 2)) then
+          
+          call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, i_generator, iproc_in $parameters )
+      end if
+    end do
+  end do
+end subroutine
+
+
+subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, i_generator, iproc_in $parameters )
+  
+  integer(bit_kind), intent(in)         :: key_in(N_int, 2), particl_1(N_int, 2), particl_2(N_int, 2)
+  integer(bit_kind)                     :: p1_mask(N_int, 2), p2_mask(N_int, 2), key_mask(N_int, 2)
+  integer,intent(in)                    :: fh1,fh2,fs1,fs2,i_generator,iproc_in
+  integer(bit_kind)                     :: miniList(N_int, 2, N_det)
+  integer                               :: n_minilist, n_alpha, n_beta, deg(2), i, ni
+  $declarations
+  
+  p1_mask(:,:) = 0
+  p2_mask(:,:) = 0
+  p1_mask(fh1/bit_kind_size + 1, fs1) = 2**(mod(fh1-1,bit_kind_size))
+  p2_mask(fh2/bit_kind_size + 1, fs2) = 2**(mod(fh2-1,bit_kind_size))
+  
+!   n_alpha = 0
+!   n_beta = 0
+  key_mask(:,:) = key_in(:,:)
+  key_mask(fh1/bit_kind_size + 1, fs1) -= 2**(mod(fh1-1,bit_kind_size))
+  key_mask(fh2/bit_kind_size + 1, fs2) -= 2**(mod(fh2-1,bit_kind_size))
+!   
+!   do i=1,N_int
+!     n_alpha = n_alpha + popcnt(key_mask(i, 1))
+!     n_beta = n_beta + popcnt(key_mask(i, 2))
+!   end do
+!   
+!   do i=1, N_det
+!     deg(1) = n_alpha
+!     deg(2) = n_beta
+!     
+!     do ni = 1, N_int
+! !       deg(1) = deg(1) - popcnt(iand(key_mask(ni, 1), psi_non_ref(ni, 1, i)))
+! !       deg(2) = deg(2) - popcnt(iand(key_mask(ni, 2), psi_non_ref(ni, 2, i)))
+!     end do
+!     
+!     
+!     if(deg(1) + deg(2) <= 2) then
+! !       ndet_out = ndet_out + 1
+! !       idx(ndet_out) = i
+!     end if
+!   end do
+
+  call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, i_generator, iproc_in $parameters )
+end subroutine
+
+
+subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, i_generator, iproc_in $parameters )
   use omp_lib
   use bitmasks
   implicit none
@@ -10,7 +122,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
   integer,parameter              :: size_max = $size_max
   $declarations
   integer          ,intent(in)   :: i_generator
-  integer(bit_kind),intent(in)   :: key_in(N_int,2)
+  integer(bit_kind),intent(in)   :: key_in(N_int,2), key_mask(N_int, 2)
   integer(bit_kind),allocatable  :: keys_out(:,:,:)
   integer(bit_kind), intent(in)  :: hole_1(N_int,2), particl_1(N_int,2)
   integer(bit_kind), intent(in)  :: hole_2(N_int,2), particl_2(N_int,2)
@@ -290,13 +402,18 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa
   double precision               :: diag_H_mat_elem
   integer(omp_lock_kind), save   :: lck, ifirst=0
   integer                        :: iproc
-
+  
+  integer(bit_kind)              :: key_mask(N_int, 2)
+  
   logical :: check_double_excitation 
+  
+  key_mask(:,:) = 0_8
+  
   iproc = iproc_in
 
   check_double_excitation = .True.
   $check_double_excitation
-
+  
 
   if (ifirst == 0) then
     ifirst=1

src/Determinants/connected_to_ref.irp.f

@@ -169,7 +169,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
   
   ! output :   0 : not connected
   !            i : connected to determinant i of the past
-  !           -i : is the ith determinant of the refernce wf keys
+  !           -i : is the ith determinant of the reference wf keys
   
   ASSERT (Nint > 0)
   ASSERT (Nint == N_int)

src/Determinants/slater_rules.irp.f

@@ -787,6 +787,7 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array)
   ASSERT (Ndet > 0)
   ASSERT (Ndet_max >= Ndet)
   i_H_psi_array = 0.d0
+  
   call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx)
   do ii=1,idx(0)
     i = idx(ii)