Merge branch 'master' of git://github.com/LCPQ/quantum_package

Conflicts:
	src/NEEDED_MODULES

scripts/generate_h_apply.py

@@ -220,7 +220,7 @@ class H_apply(object):
 
       self.data["printout_always"] = """
       do k=1,N_st
-        norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k)
+          norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k)
         delta_pt2(k) = pt2(k) - pt2_old(k)
       enddo
       """
@@ -275,7 +275,7 @@ class H_apply(object):
         if (select_max(i_generator) < selection_criterion_min*selection_criterion_factor) then
           !$ call omp_set_lock(lck)
           do k=1,N_st
-            norm_psi(k) = norm_psi(k) + psi_coef(i_generator,k)*psi_coef(i_generator,k)
+            norm_psi(k) = norm_psi(k) + psi_coef_generators(i_generator,k)*psi_coef_generators(i_generator,k)
 !            delta_pt2(k) = 0.d0
 !            pt2_old(k) = 0.d0
 !            pt2(k) = select_max(i_generator)

src/CAS_SD_selected/cas_sd.irp.f

@@ -54,4 +54,33 @@ program full_ci
       exit
     endif
   enddo
+
+  ! Check that it is a CAS-SD
+  logical :: in_cas
+  integer :: exc_max, degree_min
+  exc_max = 0
+  print *,  'CAS determinants : ', N_det_generators
+  do i=1,N_det_generators
+    do k=i,N_det_generators
+      call get_excitation_degree(psi_det_generators(1,1,k),psi_det_generators(1,1,i),degree,N_int)
+      exc_max = max(exc_max,degree)
+    enddo
+    call debug_det(psi_det_generators(1,1,i),N_int)
+    print *,  ''
+  enddo
+  print *,  'Max excitation degree in the CAS :', exc_max
+  do i=1,N_det
+    in_cas = .False.
+    degree_min = 1000
+    do k=1,N_det_generators
+      call get_excitation_degree(psi_det_generators(1,1,k),psi_det(1,1,i),degree,N_int)
+      degree_min = min(degree_min,degree)
+    enddo
+    if (degree_min > 2) then
+        print *,  'Error : This is not a CAS-SD : '
+        print *,  'Excited determinant:', degree_min
+        call debug_det(psi_det(1,1,k),N_int)
+        stop
+    endif
+  enddo
 end

src/Dets/H_apply_template.f

@@ -421,32 +421,32 @@ subroutine $subroutine($params_main)
       do k=1,N_int
         mask(k,ispin,s_hole) =                                      &
             iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen),  &
-            psi_generators(k,ispin,i_generator) )
+            psi_det_generators(k,ispin,i_generator) )
         mask(k,ispin,s_part) =                                      &
             iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen),  &
-            not(psi_generators(k,ispin,i_generator)) )
+            not(psi_det_generators(k,ispin,i_generator)) )
         mask(k,ispin,d_hole1) =                                      &
             iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen),  &
-            psi_generators(k,ispin,i_generator) )
+            psi_det_generators(k,ispin,i_generator) )
         mask(k,ispin,d_part1) =                                      &
             iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen),  &
-            not(psi_generators(k,ispin,i_generator)) )
+            not(psi_det_generators(k,ispin,i_generator)) )
         mask(k,ispin,d_hole2) =                                      &
             iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen),  &
-            psi_generators(k,ispin,i_generator) )
+            psi_det_generators(k,ispin,i_generator) )
         mask(k,ispin,d_part2) =                                      &
             iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen),  &
-            not(psi_generators(k,ispin,i_generator)) )
+            not(psi_det_generators(k,ispin,i_generator)) )
       enddo
     enddo
     if($do_double_excitations)then
-     call $subroutine_diexc(psi_generators(1,1,i_generator),          &
+     call $subroutine_diexc(psi_det_generators(1,1,i_generator),      &
          mask(1,1,d_hole1), mask(1,1,d_part1),                        &
          mask(1,1,d_hole2), mask(1,1,d_part2),                        &
          i_generator, 0 $params_post)
     endif
     if($do_mono_excitations)then
-     call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
+     call $subroutine_monoexc(psi_det_generators(1,1,i_generator),    &
          mask(1,1,s_hole ), mask(1,1,s_part ),                        &
          i_generator, 0 $params_post)
     endif
@@ -481,33 +481,33 @@ subroutine $subroutine($params_main)
       do k=1,N_int
         mask(k,ispin,s_hole) =                                      &
             iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen),  &
-            psi_generators(k,ispin,i_generator) )
+            psi_det_generators(k,ispin,i_generator) )
         mask(k,ispin,s_part) =                                      &
             iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen),  &
-            not(psi_generators(k,ispin,i_generator)) )
+            not(psi_det_generators(k,ispin,i_generator)) )
         mask(k,ispin,d_hole1) =                                      &
             iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen),  &
-            psi_generators(k,ispin,i_generator) )
+            psi_det_generators(k,ispin,i_generator) )
         mask(k,ispin,d_part1) =                                      &
             iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen),  &
-            not(psi_generators(k,ispin,i_generator)) )
+            not(psi_det_generators(k,ispin,i_generator)) )
         mask(k,ispin,d_hole2) =                                      &
             iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen),  &
-            psi_generators(k,ispin,i_generator) )
+            psi_det_generators(k,ispin,i_generator) )
         mask(k,ispin,d_part2) =                                      &
             iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen),  &
-            not (psi_generators(k,ispin,i_generator)) )
+            not (psi_det_generators(k,ispin,i_generator)) )
       enddo
     enddo
 
     if($do_double_excitations)then
-    call $subroutine_diexc(psi_generators(1,1,i_generator),          &
+      call $subroutine_diexc(psi_det_generators(1,1,i_generator),    &
         mask(1,1,d_hole1), mask(1,1,d_part1),                        &
         mask(1,1,d_hole2), mask(1,1,d_part2),                        &
         i_generator, iproc $params_post)
     endif
     if($do_mono_excitations)then
-    call $subroutine_monoexc(psi_generators(1,1,i_generator),        &
+      call $subroutine_monoexc(psi_det_generators(1,1,i_generator),  &
         mask(1,1,s_hole ), mask(1,1,s_part ),                        &
         i_generator, iproc $params_post)
     endif

src/Dets/connected_to_ref.irp.f

@@ -60,6 +60,7 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint)
   integer*8, external            :: det_search_key
   logical                        :: is_in_wavefunction
   
+  is_in_wavefunction = .False.
   get_index_in_psi_det_sorted_bit = 0
   ibegin = 1
   iend   = N_det+1
@@ -93,9 +94,10 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint)
     endif
   enddo
   i += 1
-  if (i > N_det) then
+
-    return
+!  if (i > N_det) then
-  endif
+!    return
+!  endif
 
   !DIR$ FORCEINLINE
   do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref)
@@ -104,50 +106,49 @@ integer function get_index_in_psi_det_sorted_bit(key,Nint)
       continue
     else
       is_in_wavefunction = .True.
-      !DEC$ LOOP COUNT MIN(3)
+      !DIR$ IVDEP
+      !DIR$ LOOP COUNT MIN(3)
       do l=2,Nint
         if ( (key(l,1) /= psi_det_sorted_bit(l,1,i)).or.                           &
               (key(l,2) /= psi_det_sorted_bit(l,2,i)) ) then
           is_in_wavefunction = .False.
-          exit
         endif
       enddo
       if (is_in_wavefunction) then
+        get_index_in_psi_det_sorted_bit = i
         exit
+!        return
       endif
     endif
     i += 1
     if (i > N_det) then
-      return
+      exit
- !    exit
+!      return
     endif
     
   enddo
-  if (is_in_wavefunction) then
-    get_index_in_psi_det_sorted_bit = i
-  endif
 
 ! DEBUG is_in_wf
-! if (is_in_wavefunction) then
+ if (is_in_wavefunction) then
-!   degree = 1
+   degree = 1
-!   do i=1,N_det
+   do i=1,N_det
-!     integer                        :: degree
+     integer                        :: degree
-!     call get_excitation_degree(key,psi_det(1,1,i),degree,N_int)
+     call get_excitation_degree(key,psi_det(1,1,i),degree,N_int)
-!     if (degree == 0) then
+     if (degree == 0) then
-!       exit
+       exit
-!     endif
+     endif
-!   enddo
+   enddo
-!   if (degree /=0) then
+   if (degree /=0) then
-!     stop 'pouet 1'
+     stop 'pouet 1'
-!   endif
+   endif
-! else
+ else
-!   do i=1,N_det
+   do i=1,N_det
-!     call get_excitation_degree(key,psi_det(1,1,i),degree,N_int)
+     call get_excitation_degree(key,psi_det(1,1,i),degree,N_int)
-!     if (degree == 0) then
+     if (degree == 0) then
-!       stop 'pouet 2'
+       stop 'pouet 2'
-!     endif
+     endif
-!   enddo
+   enddo
-! endif
+ endif
 ! END DEBUG is_in_wf
 end
 

src/Generators_CAS/generators.irp.f

@@ -15,9 +15,9 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
       do k=1,N_int
         good = good .and. (                                          &
             iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) ==         &
-            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. (  &
+            iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) ) .and. ( &
             iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) ==         &
-            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
+            iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2)) )
       enddo
       if (good) then
         exit
@@ -31,7 +31,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
   call write_int(output_dets,N_det_generators,'Number of generators')
 END_PROVIDER
 
-BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
   implicit none
   BEGIN_DOC
   ! For Single reference wave functions, the generator is the
@@ -44,11 +45,11 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
     do l=1,n_cas_bitmask
       good = .True.
       do k=1,N_int
-        good = good .and. (                                          &
+        good = good .and. (                                         &
             iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) ==         &
-            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. (  &
+            iand(not(cas_bitmask(k,1,l)), HF_bitmask(k,1)) .and. (   &
             iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) ==         &
-            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
+            iand(not(cas_bitmask(k,2,l)), HF_bitmask(k,2) )) )
       enddo
       if (good) then
         exit
@@ -57,10 +58,11 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
     if (good) then
       m = m+1
       do k=1,N_int
-        psi_generators(k,1,m) = psi_det(k,1,i)
+        psi_det_generators(k,1,m) = psi_det(k,1,i)
-        psi_generators(k,2,m) = psi_det(k,2,i)
+        psi_det_generators(k,2,m) = psi_det(k,2,i)
       enddo
-!     call debug_det(psi_generators(1,1,m),N_int)
+      psi_coef_generators(m,:) = psi_coef(m,:)
+!     call debug_det(psi_det_generators(1,1,m),N_int)
     endif
   enddo
   

src/Generators_full/generators.irp.f

@@ -34,7 +34,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  call write_int(output_dets,N_det_generators,'Number of generators')
 END_PROVIDER
 
-BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
  implicit none
  BEGIN_DOC
  ! For Single reference wave functions, the generator is the
@@ -43,9 +44,10 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
  integer                        :: i, k
  do i=1,N_det
    do k=1,N_int
-     psi_generators(k,1,i) = psi_det_sorted(k,1,i)
+     psi_det_generators(k,1,i) = psi_det_sorted(k,1,i)
-     psi_generators(k,2,i) = psi_det_sorted(k,2,i)
+     psi_det_generators(k,2,i) = psi_det_sorted(k,2,i)
    enddo
+   psi_coef_generators(i,:) = psi_coef_sorted(i,:)
  enddo
 
 END_PROVIDER
@@ -58,7 +60,7 @@ BEGIN_PROVIDER [integer, degree_max_generators]
  integer :: i,degree
   degree_max_generators = 0
   do i = 1, N_det_generators
-   call get_excitation_degree(HF_bitmask,psi_generators(1,1,i),degree,N_int)
+   call get_excitation_degree(HF_bitmask,psi_det_generators(1,1,i),degree,N_int)
    if(degree .gt. degree_max_generators)then
     degree_max_generators = degree
    endif

src/Generators_restart/generators.irp.f

@@ -17,8 +17,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
 END_PROVIDER
 
 
- BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,N_det_generators) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
-&BEGIN_PROVIDER [ double precision, psi_generators_coef, (N_det_generators,N_states)
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
  implicit none
  BEGIN_DOC
  ! read wf
@@ -29,11 +29,11 @@ END_PROVIDER
  if(ifirst == 0)then
   do i=1,N_det_generators
     do k=1,N_int
-      psi_generators(k,1,i) = psi_det(k,1,i)
+      psi_det_generators(k,1,i) = psi_det(k,1,i)
-      psi_generators(k,2,i) = psi_det(k,2,i)
+      psi_det_generators(k,2,i) = psi_det(k,2,i)
     enddo
    do k = 1, N_states
-    psi_generators_coef(i,k) = psi_coef(i,k)
+    psi_coef_generators(i,k) = psi_coef(i,k)
    enddo
   enddo
   ifirst = 1

src/MRCC/H_apply.irp.f

@@ -0,0 +1,26 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import *
+
+s = H_apply("mrcc")
+s.data["parameters"] = ", delta_ij_sd_, Ndet_sd"
+s.data["declarations"] += """
+    integer, intent(in) :: Ndet_sd
+    double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*)
+"""
+s.data["keys_work"] = "call mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,key_idx,keys_out,N_int,iproc)"
+s.data["params_post"] += ", delta_ij_sd_, Ndet_sd"
+s.data["params_main"] += "delta_ij_sd_, Ndet_sd"
+s.data["decls_main"] += """
+    integer, intent(in) :: Ndet_sd
+    double precision, intent(in) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*)
+"""
+s.data["finalization"] = ""
+s.data["copy_buffer"] = ""
+s.data["generate_psi_guess"] = ""
+s.data["size_max"] = "256"
+
+print s
+
+END_SHELL
+

src/MRCC/Makefile

@@ -0,0 +1,6 @@
+# Define here all new external source files and objects.Don't forget to prefix the
+# object files with IRPF90_temp/
+SRC=
+OBJ=
+
+include $(QPACKAGE_ROOT)/src/Makefile.common

src/MRCC/NEEDED_MODULES

@@ -0,0 +1,2 @@
+AOs BiInts Bitmask CAS_SD_selected Dets Electrons Ezfio_files Generators_CAS Generators_full Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils
+

src/MRCC/README.rst

@@ -0,0 +1,14 @@
+-4.142795384334731
+
+4.695183071437694E-002
+ Determinant           64
+ ---------------------------------------------
+ 000000000000002E|000000000000002E
+ |-+++-+----------------------------------------------------------|
+ |-+++-+----------------------------------------------------------|
+
+CAS-SD: -4.14214374069306
+      : -4.14230904320551
+
+E0 = -11.5634986758976
+

src/MRCC/mrcc.irp.f

@@ -0,0 +1,43 @@
+program mrcc
+  implicit none
+  read_wf = .True.
+  TOUCH read_wf
+  call run
+end
+
+subroutine run
+  implicit none
+
+  print *, N_det
+  print *, N_det_cas
+  print *, N_det_sd
+
+!  call update_generators
+  integer :: i
+  print *,  'CAS'
+  print *,  '==='
+  do i=1,N_det_cas
+    print *,  psi_cas_coefs(i,:)
+    call debug_det(psi_cas(1,1,i),N_int)
+  enddo
+
+!  print *,  'SD'
+!  print *,  '=='
+!  do i=1,N_det_sd
+!    print *,  psi_sd_coefs(i,:)
+!    call debug_det(psi_sd(1,1,i),N_int)
+!  enddo
+!
+  print *,  'MRCC'
+  print *,  '===='
+  print *,  ci_energy(:)
+  print *,  h_matrix_all_dets(3,3), delta_ij(3,3,1)
+  print *,  h_matrix_all_dets(3,3), delta_ij(3,3,1)
+  print *,  ci_energy_dressed(:)
+!  print *,  'max', maxval(delta_ij_sd)
+!  print *,  'min', minval(delta_ij_sd)
+!
+!  do i=1,N_det
+!     print '(10(F10.6,X))',  delta_ij(i,1:N_det,1)
+!  enddo
+end

src/MRCC/mrcc_dress.irp.f

@@ -0,0 +1,190 @@
+subroutine mrcc_dress(delta_ij_sd_,Ndet_sd,i_generator,n_selected,det_buffer,Nint,iproc)
+ use bitmasks
+ implicit none
+
+  integer, intent(in)            :: i_generator,n_selected, Nint, iproc
+  integer, intent(in) :: Ndet_sd
+  double precision, intent(inout) :: delta_ij_sd_(Ndet_sd,Ndet_sd,*)
+
+  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
+  integer                        :: i,j,k,m
+  integer                        :: new_size
+  logical                        :: is_in_wavefunction
+  integer                        :: degree(psi_det_size)
+  integer                        :: idx(0:psi_det_size)
+  logical                        :: good
+
+  integer(bit_kind)              :: tq(Nint,2,n_selected)
+  integer                        :: N_tq, c_ref
+  integer                        :: connected_to_ref
+
+  N_tq = 0
+  do i=1,N_selected
+    c_ref = connected_to_ref(det_buffer(1,1,i),psi_det_generators,Nint, &
+       i_generator,N_det_generators)
+
+    if (c_ref /= 0) then
+      cycle
+    endif
+
+    ! Select determinants that are triple or quadruple excitations
+    ! from the CAS
+    good = .True.
+    call get_excitation_degree_vector(psi_cas,det_buffer(1,1,i),degree,Nint,N_det_cas,idx)
+    do k=1,idx(0)
+      if (degree(k) < 3) then
+        good = .False.
+        exit
+      endif
+    enddo
+    if (good) then
+      if (.not. is_in_wavefunction(det_buffer(1,1,i),Nint,N_det)) then
+        N_tq += 1
+        do k=1,N_int
+          tq(k,1,N_tq) = det_buffer(k,1,i)
+          tq(k,2,N_tq) = det_buffer(k,2,i)
+        enddo
+      endif
+    endif
+  enddo
+
+  ! Compute <k|H|a><a|H|j> / (E0 - Haa)
+  double precision :: hka, haa
+  double precision :: haj
+  double precision :: f(N_states)
+
+  do i=1,N_tq
+    call get_excitation_degree_vector(psi_sd,tq(1,1,i),degree,Nint,Ndet_sd,idx)
+    call i_h_j(tq(1,1,i),tq(1,1,i),Nint,haa)
+    do m=1,N_states
+      f(m) = 1.d0/(ci_electronic_energy(m)-haa)
+    enddo
+    do k=1,idx(0)
+      call i_h_j(tq(1,1,i),psi_sd(1,1,idx(k)),Nint,hka)
+      do j=k,idx(0)
+        call i_h_j(tq(1,1,i),psi_sd(1,1,idx(j)),Nint,haj)
+        do m=1,N_states
+          delta_ij_sd_(idx(k), idx(j),m) += haj*hka* f(m)
+          delta_ij_sd_(idx(j), idx(k),m) += haj*hka* f(m)
+        enddo
+      enddo 
+    enddo
+  enddo
+end
+
+BEGIN_PROVIDER [ double precision, delta_ij_sd, (N_det_sd, N_det_sd,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Dressing matrix in SD basis
+ END_DOC
+ delta_ij_sd = 0.d0
+ call H_apply_mrcc(delta_ij_sd,N_det_sd)
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, delta_ij, (N_det,N_det,N_states) ]
+ implicit none
+ BEGIN_DOC
+ ! Dressing matrix in N_det basis
+ END_DOC
+ integer :: i,j,m
+ delta_ij = 0.d0
+ do m=1,N_states
+  do j=1,N_det_sd
+   do i=1,N_det_sd
+     delta_ij(idx_sd(i),idx_sd(j),m) = delta_ij_sd(i,j,m) 
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det) ]
+ implicit none
+ BEGIN_DOC
+ ! Dressed H with Delta_ij
+ END_DOC
+ integer                        :: i, j
+ do j=1,N_det
+   do i=1,N_det
+     h_matrix_dressed(i,j) = h_matrix_all_dets(i,j) + delta_ij(i,j,1)
+     if (i==j) then
+       print *,  i, delta_ij(i,j,1), h_matrix_all_dets(i,j) 
+     endif
+   enddo
+ enddo
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, CI_electronic_energy_dressed, (N_states_diag) ]
+&BEGIN_PROVIDER [ double precision, CI_eigenvectors_dressed, (N_det,N_states_diag) ]
+&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_dressed, (N_states_diag) ]
+  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_eigenvectors_dressed(i,j) = psi_coef(i,j)
+    enddo
+  enddo
+  
+  if (diag_algorithm == "Davidson") then
+    
+    stop 'use Lapack'
+!    call davidson_diag(psi_det,CI_eigenvectors_dressed,CI_electronic_energy_dressed, &
+!        size(CI_eigenvectors_dressed,1),N_det,N_states_diag,N_int,output_Dets)
+    
+  else if (diag_algorithm == "Lapack") then
+    
+    double precision, allocatable  :: eigenvectors(:,:), eigenvalues(:)
+    allocate (eigenvectors(size(H_matrix_dressed,1),N_det))
+    allocate (eigenvalues(N_det))
+    call lapack_diag(eigenvalues,eigenvectors,                       &
+        H_matrix_dressed,size(H_matrix_dressed,1),N_det)
+    CI_electronic_energy_dressed(:) = 0.d0
+    do i=1,N_det
+       CI_eigenvectors_dressed(i,1) = eigenvectors(i,1)
+    enddo
+    integer :: i_state
+    double precision :: s2
+    i_state = 0
+    do j=1,N_det
+      call get_s2_u0(psi_det,eigenvectors(1,j),N_det,N_det,s2)
+      if(dabs(s2-expected_s2).le.0.3d0)then
+       i_state += 1
+       do i=1,N_det
+         CI_eigenvectors_dressed(i,i_state) = eigenvectors(i,j)
+       enddo
+       CI_electronic_energy_dressed(i_state) = eigenvalues(j)
+       CI_eigenvectors_s2_dressed(i_state) = s2
+      endif
+      if (i_state.ge.N_states_diag) then
+        exit
+      endif
+    enddo
+    deallocate(eigenvectors,eigenvalues)
+  endif
+  
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
+  implicit none
+  BEGIN_DOC
+  ! N_states lowest eigenvalues of the dressed CI matrix
+  END_DOC
+  
+  integer                        :: j
+  character*(8)                  :: st
+  call write_time(output_Dets)
+  do j=1,N_states_diag
+    CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion
+    write(st,'(I4)') j
+    call write_double(output_Dets,CI_energy(j),'Energy of state '//trim(st))
+    call write_double(output_Dets,CI_eigenvectors_s2(j),'S^2 of state '//trim(st))
+  enddo
+
+END_PROVIDER
+

src/MRCC/mrcc_utils.irp.f

@@ -0,0 +1,150 @@
+ use bitmasks
+
+BEGIN_PROVIDER [ integer(bit_kind), cas_bitmask, (N_int,2,N_cas_bitmask) ]
+ implicit none
+ BEGIN_DOC
+ ! Bitmasks for CAS reference determinants. (N_int, alpha/beta, CAS reference)
+ END_DOC
+ logical                        :: exists
+ integer                        :: i
+ PROVIDE ezfio_filename
+  
+ call ezfio_has_bitmasks_cas(exists)
+ if (exists) then                
+   call ezfio_get_bitmasks_cas(cas_bitmask)
+ else
+   do i=1,N_cas_bitmask
+     cas_bitmask(:,:,i) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
+   enddo 
+ endif  
+            
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, N_det_cas ]
+  implicit none
+  BEGIN_DOC
+  ! Number of generator detetrminants
+  END_DOC
+  integer                        :: i,k,l
+  logical                        :: good
+  call write_time(output_dets)
+  N_det_cas = 0
+  do i=1,N_det
+    do l=1,n_cas_bitmask
+      good = .True.
+      do k=1,N_int
+        good = good .and. (                                          &
+            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) ==         &
+            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. (  &
+            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) ==         &
+            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
+      enddo
+      if (good) then
+        exit
+      endif
+    enddo
+    if (good) then
+      N_det_cas += 1
+    endif
+  enddo
+  N_det_cas = max(N_det_cas, 1)
+  call write_int(output_dets,N_det_cas, 'Number of determinants in the CAS')
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_cas) ]
+&BEGIN_PROVIDER [ double precision, psi_cas_coefs,  (N_det_cas,n_states) ]
+&BEGIN_PROVIDER [ integer, idx_cas, (N_det_cas) ]
+  implicit none
+  BEGIN_DOC
+  ! For Single reference wave functions, the generator is the
+  ! Hartree-Fock determinant
+  END_DOC
+  integer                        :: i, k, l, m
+  logical                        :: good
+  m=0
+  do i=1,N_det
+    do l=1,n_cas_bitmask
+      good = .True.
+      do k=1,N_int
+        good = good .and. (                                          &
+            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,i)) ==         &
+            iand(not(cas_bitmask(k,1,l)), psi_det(k,1,1)) ) .and. (  &
+            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,i)) ==         &
+            iand(not(cas_bitmask(k,2,l)), psi_det(k,2,1)) )
+      enddo
+      if (good) then
+        exit
+      endif
+    enddo
+    if (good) then
+      m = m+1
+      do k=1,N_int
+        psi_cas(k,1,m) = psi_det(k,1,i)
+        psi_cas(k,2,m) = psi_det(k,2,i)
+      enddo
+      idx_cas(m) = i
+      do k=1,N_states
+        psi_cas_coefs(m,k) = psi_coef(i,k)
+      enddo
+    endif
+  enddo
+
+END_PROVIDER
+
+
+
+
+ BEGIN_PROVIDER [ integer(bit_kind), psi_sd,  (N_int,2,N_det) ]
+&BEGIN_PROVIDER [ double precision, psi_sd_coefs, (N_det,n_states) ]
+&BEGIN_PROVIDER [ integer, idx_sd,  (N_det) ]
+&BEGIN_PROVIDER [ integer, N_det_sd]
+ implicit none
+ BEGIN_DOC
+ ! SD
+ END_DOC
+ integer                        :: i_sd,j,k
+ integer                        :: degree
+ logical                        :: in_cas
+ i_sd =0
+ do k=1,N_det
+   in_cas = .False.
+   do j=1,N_det_cas
+     call get_excitation_degree(psi_cas(1,1,j), psi_det(1,1,k), degree, N_int)
+     if (degree == 0) then
+       in_cas = .True.
+       exit
+     endif
+   enddo
+   if (.not.in_cas) then
+     double precision :: hij
+     i_sd += 1
+     psi_sd(1:N_int,1:2,i_sd) = psi_det(1:N_int,1:2,k)
+     psi_sd_coefs(i_sd,1:N_states) = psi_coef(k,1:N_states)
+     idx_sd(i_sd) = k
+   endif
+ enddo
+ N_det_sd = i_sd
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, lambda_mrcc, (psi_det_size,n_states) ]
+ implicit none
+ BEGIN_DOC
+ ! cm/<Psi_0|H|D_m>
+ END_DOC
+ integer :: i,k
+ double precision :: ihpsi(N_states)
+ do i=1,N_det_sd
+   call i_h_psi(psi_sd(1,1,i), psi_cas, psi_cas_coefs, N_int, N_det_cas, &
+     size(psi_cas_coefs,1), n_states, ihpsi)
+   double precision :: hij
+   do k=1,N_states
+     if (dabs(ihpsi(k)) < 1.d-6) then
+       lambda_mrcc(i,k) = 0.d0
+     else
+       lambda_mrcc(i,k) = psi_sd_coefs(i,k)/ihpsi(k)
+       lambda_mrcc(i,k) = min( lambda_mrcc (i,k),0.d0 )
+     endif
+   enddo
+ enddo
+END_PROVIDER
+

src/Makefile.common

@@ -8,6 +8,9 @@ default:
 veryclean:
 	 clean_modules.sh
 
+clean:
+	IN_MAKE=1 make clean
+
 else  # Called by scripts/build_module.sh
 
 default: all .gitignore

src/NEEDED_MODULES

@@ -1 +1 @@
-AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected
+AOs Bielec_integrals Bitmask CID CID_SC2_selected CID_selected CIS CISD CISD_selected CISD_SC2_selected Dets Electrons Ezfio_files Full_CI Generators_full Hartree_Fock MOGuess MonoInts MOs MP2 Nuclei Output Selectors_full Utils Molden FCIdump Generators_CAS CAS_SD_selected DDCI_selected MRCC

src/Perturbation/perturbation_template.f

@@ -25,7 +25,7 @@ subroutine perturb_buffer_$PERT(i_generator,buffer,buffer_size,e_2_pert_buffer,c
   ASSERT (N_st > 0)
   do i = 1,buffer_size
 
-    c_ref = connected_to_ref(buffer(1,1,i),psi_generators,Nint,i_generator,N_det)
+    c_ref = connected_to_ref(buffer(1,1,i),psi_det_generators,Nint,i_generator,N_det_generators)
 
     if (c_ref /= 0) then
       cycle
@@ -76,7 +76,7 @@ subroutine perturb_buffer_by_mono_$PERT(i_generator,buffer,buffer_size,e_2_pert_
   ASSERT (N_st > 0)
   do i = 1,buffer_size
 
-    c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_generators,Nint,i_generator,N_det)
+    c_ref = connected_to_ref_by_mono(buffer(1,1,i),psi_det_generators,Nint,i_generator,N_det)
 
     if (c_ref /= 0) then
       cycle

src/SingleRefMethod/generators.irp.f

@@ -9,19 +9,20 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
  N_det_generators = 1
 END_PROVIDER
 
-BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
+ BEGIN_PROVIDER [ integer(bit_kind), psi_det_generators, (N_int,2,psi_det_size) ]
+&BEGIN_PROVIDER [ double precision, psi_coef_generators, (psi_det_size,N_states) ]
  implicit none
  BEGIN_DOC
  ! For Single reference wave functions, the generator is the
  ! Hartree-Fock determinant
  END_DOC
- psi_generators = 0_bit_kind
+ psi_det_generators = 0_bit_kind
  integer :: i,j,k
  integer :: degree
 
  do i=1,N_int
-   psi_generators(i,1,1) = HF_bitmask(i,1)
+   psi_det_generators(i,1,1) = HF_bitmask(i,1)
-   psi_generators(i,2,1) = HF_bitmask(i,2)
+   psi_det_generators(i,2,1) = HF_bitmask(i,2)
  enddo
 
  do j=1,N_det
@@ -32,12 +33,8 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_generators, (N_int,2,psi_det_size) ]
    endif
  end do
 
- do j=2,k
+ psi_det_generators(:,:,1) = psi_det(:,:,j)
-   psi_generators(:,:,j) = psi_det(:,:,j-1)
+ psi_coef_generators(1,:) = psi_coef_generators(j,:)
- enddo
- do j=k+1,N_det
-   psi_generators(:,:,j) = psi_det(:,:,j)
- enddo
 
 END_PROVIDER