init mrcepa0

plugins/MRCC_Utils/mrcc_utils.irp.f

@@ -110,6 +110,7 @@ END_PROVIDER
 
  BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
 &BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
+&BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref,N_det_ref,N_states) ]
  implicit none
  BEGIN_DOC
  ! Dressing matrix in N_det basis
@@ -117,6 +118,7 @@ END_PROVIDER
  integer :: i,j,m
  delta_ij = 0.d0
  delta_ii = 0.d0
+
  call H_apply_mrcc(delta_ij,delta_ii,N_det_ref,N_det_non_ref)
  double precision :: max_delta
  double precision :: accu
@@ -157,6 +159,17 @@ BEGIN_PROVIDER [ double precision, h_matrix_dressed, (N_det,N_det,N_states) ]
        h_matrix_dressed(i,j,istate) = h_matrix_all_dets(i,j) 
      enddo
    enddo
+   
+   !!!!!!!!!!
+   do ii = 1, N_det_ref
+   do jj = 1, N_det_ref
+    i = idx_ref(ii)
+    j = idx_ref(jj)
+    h_matrix_dressed(i,j,istate) += delta_cas(ii,jj,istate)
+    h_matrix_dressed(j,i,istate) += delta_cas(ii,jj,istate)
+   end do
+   end do
+   !!!!!!!!!!!!!
    do ii = 1, N_det_ref
      i =idx_ref(ii)
      h_matrix_dressed(i,i,istate) += delta_ii(ii,istate)

plugins/mrcepa0/dressing.irp.f

@@ -0,0 +1,234 @@
+use bitmasks
+
+
+ BEGIN_PROVIDER [ integer, cepa0_shortcut, (0:N_det_non_ref+1) ]
+&BEGIN_PROVIDER [ integer, det_cepa0_idx, (N_det_non_ref) ]
+&BEGIN_PROVIDER [ integer(bit_kind), det_cepa0_active, (N_det_non_ref) ]
+&BEGIN_PROVIDER [ integer(bit_kind), det_ref_active, (N_det_ref) ]
+  use bitmasks
+  implicit none
+  
+  integer(bit_kind) :: det_noactive(N_int, 2, N_det_non_ref), active_sorb(2)
+  integer i, II, j, k
+  logical, external :: detEq
+  
+  print *, "provide cepa0"
+  active_sorb = (/b'001100000000', b'001100000000'/)
+  do i=1, N_det_non_ref
+    det_noactive(1,1,i) = iand(psi_non_ref(1,1,i), not(active_sorb(1)))
+    det_noactive(1,2,i) = iand(psi_non_ref(1,2,i), not(active_sorb(2))) 	 	
+  end do
+  
+  call sort_dets_ab(det_noactive, det_cepa0_idx, cepa0_shortcut, N_det_non_ref, N_int)
+  do i=1,N_det_ref
+    det_ref_active(i) = iand(psi_ref(1,1,i), active_sorb(1))
+    det_ref_active(i) = det_ref_active(i) + iand(psi_ref(1,2,i), active_sorb(2)) * 2_8**32_8
+  end do
+  
+  cepa0_shortcut(0) = 1
+  cepa0_shortcut(1) = 1
+  det_cepa0_active(1) = iand(psi_non_ref(1,1,det_cepa0_idx(1)), active_sorb(1))
+  det_cepa0_active(1) = det_cepa0_active(1) + iand(psi_non_ref(1,2,det_cepa0_idx(1)), active_sorb(2)) * 2_8**32_8
+  
+  do i=2,N_det_non_ref
+    det_cepa0_active(i) = iand(psi_non_ref(1,1,det_cepa0_idx(i)), active_sorb(1))
+    det_cepa0_active(i) = det_cepa0_active(i) + iand(psi_non_ref(1,2,det_cepa0_idx(i)), active_sorb(2)) * 2_8**32_8
+    
+    if(.not. detEq(det_noactive(1,1,i), det_noactive(1,1,i-1), N_int)) then
+      cepa0_shortcut(0) += 1
+      cepa0_shortcut(cepa0_shortcut(0)) = i
+    end if
+  end do
+  cepa0_shortcut(0) += 1
+  cepa0_shortcut(cepa0_shortcut(0)) = N_det_non_ref+1
+  
+END_PROVIDER
+
+
+
+ BEGIN_PROVIDER [ double precision, delta_cas, (N_det_ref, N_det_ref, N_states) ]
+  use bitmasks
+  implicit none
+  integer :: i,j,k
+  double precision :: Hjk, Hki
+  integer i_state
+  
+  i_state = 1
+  
+  do i=1,N_det_ref
+    do j=1,i
+      delta_cas(i,j,i_state) = 0d0
+      do k=1,N_det_non_ref
+        call i_h_j(psi_ref(1,1,j), psi_non_ref(1,1,k),N_int,Hjk)
+        call i_h_j(psi_ref(1,1,i), psi_non_ref(1,1,k),N_int,Hki)
+        delta_cas(i,j,i_state) += Hjk * Hki * lambda_mrcc(i_state, k)
+      end do
+      delta_cas(j,i,i_state) = delta_cas(i,j,i_state)
+    end do
+  end do
+  
+  print *, "mrcepa0_cas_dressing",  delta_cas(:,:,1)
+ END_PROVIDER
+ 
+logical function detEq(a,b,Nint)
+   use bitmasks
+   implicit none
+   integer, intent(in) :: Nint
+   integer(bit_kind), intent(in) :: a(Nint,2), b(Nint,2)
+   integer :: ni, i
+ 
+   detEq = .false.
+   do i=1,2
+   do ni=1,Nint
+     if(a(ni,i) /= b(ni,i)) return
+   end do
+   end do
+   detEq = .true.
+end function
+
+
+
+ BEGIN_PROVIDER [ double precision, delta_ij, (N_det_ref,N_det_non_ref,N_states) ]
+&BEGIN_PROVIDER [ double precision, delta_ii, (N_det_ref,N_states) ]
+  use bitmasks
+  implicit none
+  
+  integer                         :: i_state, i, i_I, J, k, degree, degree2, m, l, deg, ni
+  integer                         :: p1,p2,h1,h2,s1,s2, p1_,p2_,h1_,h2_,s1_,s2_
+  logical                         :: ok
+  double precision                :: phase_iI, phase_Ik, phase_Jl, phase_IJ, phase_al, diI, hIi, hJi, delta_JI, dkI(1), HkI, ci_inv(1), dia_hla(1)
+  double precision                :: contrib
+  integer, dimension(0:2,2,2)     :: exc_iI, exc_Ik, exc_IJ
+  integer(bit_kind)               :: det_tmp(N_int, 2), made_hole, made_particle
+  integer, allocatable            :: idx_sorted_bit(:)
+  integer, external               :: get_index_in_psi_det_sorted_bit
+  logical, external               :: is_in_wavefunction
+  
+  integer :: II, blok
+  
+  provide det_cepa0_active
+  
+  if(N_int /= 1) then
+    print *, "mrcepa0 experimental N_int==1"
+    stop
+  end if
+  
+  i_state = 1
+  delta_ii(:,:) = 0
+  delta_ij(:,:,:) = 0
+  
+!   do i=1,N_det_ref
+!     delta_ii(i,i_state) = delta_cas(i,i)
+!   end do
+  
+  provide mo_bielec_integrals_in_map
+  allocate(idx_sorted_bit(N_det))
+  
+  idx_sorted_bit(:) = -1
+  do i=1,N_det_non_ref
+    idx_sorted_bit(get_index_in_psi_det_sorted_bit(psi_non_ref(1,1,i), N_int)) = i
+  enddo
+  
+  !sd $OMP PARALLEL DO schedule(dynamic,10) default(firstprivate) shared(delta_ij, delta_ii)
+  do blok=1,cepa0_shortcut(0)
+  do i=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
+    do II=1,N_det_ref
+      
+      made_hole = iand(det_cepa0_active(i), xor(det_cepa0_active(i), det_ref_active(II)))
+      made_particle = iand(det_ref_active(II), xor(det_cepa0_active(i), det_ref_active(II)))
+      
+      if(popcnt(made_hole) + popcnt(made_particle) > 2) cycle
+      
+      do k=cepa0_shortcut(blok), cepa0_shortcut(blok+1)-1
+        do J=1,N_det_ref
+          if(iand(made_hole, det_ref_active(J)) /= made_hole) cycle
+          if(iand(made_particle, det_ref_active(J)) /= 0) cycle
+          call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,det_cepa0_idx(k)),N_int,Hki)
+          contrib = Hki * lambda_mrcc(i_state, det_cepa0_idx(k)) * delta_cas(II,J,i_state)
+          delta_ij(II, det_cepa0_idx(i), i_state) += contrib
+          
+!                     if(dabs(psi_ref_coef(II,i_state)).ge.5.d-5) then
+!                       !qs$OMP CRITICAL 
+!                       delta_ii(II,i_state) -= contrib / psi_ref_coef(II, i_state) * psi_non_ref_coef(det_cepa0_idx(k),i_state)
+!                       !dsd$OMP END CRITICAL
+!                     endif
+
+        end do
+      end do
+    end do
+  end do
+  end do
+  !qsd $OMP END PARALLEL DO
+  deallocate(idx_sorted_bit)
+END_PROVIDER
+
+
+subroutine set_det_bit(det, p, s)
+  use bitmasks
+  implicit none
+  integer(bit_kind),intent(inout) :: det(N_int, 2)
+  integer, intent(in) :: p, s
+  integer :: ni, pos
+  
+  ni = (p-1)/bit_kind_size + 1
+  pos = mod(p-1, bit_kind_size)
+  det(ni,s) = ibset(det(ni,s), pos)
+end subroutine
+  
+  
+  
+subroutine apply_excitation(det, exc, res, ok, Nint)
+  use bitmasks
+  implicit none
+  
+  integer, intent(in) :: Nint
+  integer, intent(in) :: exc(0:2,2,2)
+  integer(bit_kind),intent(in) :: det(Nint, 2)
+  integer(bit_kind),intent(out) :: res(Nint, 2)
+  logical, intent(out) :: ok
+  integer :: h1,p1,h2,p2,s1,s2,degree
+  integer :: ii, pos
+  
+  
+  ok = .false.
+  degree = exc(0,1,1) + exc(0,1,2)
+  if(.not. (degree > 0 .and. degree <= 2)) then
+    print *, "apply ex"
+    STOP
+  endif
+  
+  call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
+  res = det
+  
+  ii = (h1-1)/bit_kind_size + 1
+  pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
+  if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
+  res(ii, s1) = ibclr(res(ii, s1), pos)
+  
+  ii = (p1-1)/bit_kind_size + 1
+  pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
+  if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
+  res(ii, s1) = ibset(res(ii, s1), pos)
+  
+  if(degree == 2) then
+    ii = (h2-1)/bit_kind_size + 1
+    pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
+    if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
+    res(ii, s2) = ibclr(res(ii, s2), pos)
+    
+    ii = (p2-1)/bit_kind_size + 1
+    pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
+    if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
+    res(ii, s2) = ibset(res(ii, s2), pos)
+  endif
+
+  ok = .true.
+end subroutine
+
+
+
+
+
+
+
+

plugins/mrcepa0/mrcepa0.irp.f

@@ -0,0 +1,24 @@
+program mrcepa0
+  implicit none
+  if (.not.read_wf) then
+    print *,  'read_wf has to be true.'
+    stop 1
+  endif
+  call print_cas_coefs
+  call run_mrcepa0
+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/mrcepa0/mrcepa0_general.irp.f

@@ -0,0 +1,96 @@
+subroutine run_mrcepa0
+  implicit none
+  call set_generators_bitmasks_as_holes_and_particles
+  call mrcepa0_iterations
+end
+
+subroutine mrcepa0_iterations
+  implicit none
+  
+  integer :: i,j
+
+  double precision :: E_new, E_old, delta_e
+  integer :: iteration,i_oscillations
+  double precision :: E_past(4)
+  E_new = 0.d0
+  delta_E = 1.d0
+  iteration = 0
+  j = 1
+  i_oscillations = 0
+  do while (delta_E > 1.d-7)
+    iteration += 1
+    print *,  '===========================' 
+    print *,  'MRCEPA0 Iteration', iteration
+    print *,  '===========================' 
+    print *,  ''
+    E_old = sum(ci_energy_dressed)
+    call write_double(6,ci_energy_dressed(1),"MRCEPA0 energy")
+    call diagonalize_ci_dressed
+    E_new = sum(ci_energy_dressed)
+    delta_E = dabs(E_new - E_old)
+    E_past(j) = E_new
+    j +=1
+    if(j>4)then
+     j=1
+    endif
+    if(iteration > 4) then 
+     if(delta_E > 1.d-10)then
+      if(dabs(E_past(1) - E_past(3)) .le. delta_E .and. dabs(E_past(2) - E_past(4)).le. delta_E)then 
+        print*,'OSCILLATIONS !!!'
+!        oscillations = .True.
+       i_oscillations +=1
+       lambda_mrcc_tmp = lambda_mrcc
+      endif
+     endif
+    endif
+    call save_wavefunction
+!   if (i_oscillations > 5) then
+!    exit
+!   endif
+    if (iteration > 100) then
+      exit
+    endif
+    print*,'------------'
+    print*,'VECTOR'
+    do i = 1, N_det_ref
+     print*,''
+     print*,'psi_ref_coef(i,1) = ',psi_ref_coef(i,1)
+     print*,'delta_ii(i,1)     = ',delta_ii(i,1)
+    enddo
+    print*,'------------'
+  enddo
+  call write_double(6,ci_energy_dressed(1),"Final MRCEPA0 energy")
+  call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
+  call save_wavefunction
+
+end
+
+subroutine set_generators_bitmasks_as_holes_and_particles
+ implicit none
+ integer :: i,k
+ do k = 1, N_generators_bitmask
+  do i = 1, N_int
+   ! Pure single part 
+   generators_bitmask(i,1,1,k) = holes_operators(i,1)   ! holes for pure single exc alpha 
+   generators_bitmask(i,1,2,k) = particles_operators(i,1) ! particles for pure single exc alpha 
+   generators_bitmask(i,2,1,k) = holes_operators(i,2)   ! holes for pure single exc beta 
+   generators_bitmask(i,2,2,k) = particles_operators(i,2) ! particles for pure single exc beta 
+
+   ! Double excitation 
+   generators_bitmask(i,1,3,k) = holes_operators(i,1)   ! holes for first single exc alpha 
+   generators_bitmask(i,1,4,k) = particles_operators(i,1) ! particles for first single exc alpha 
+   generators_bitmask(i,2,3,k) = holes_operators(i,2)   ! holes for first single exc beta 
+   generators_bitmask(i,2,4,k) = particles_operators(i,2) ! particles for first single exc beta 
+
+   generators_bitmask(i,1,5,k) = holes_operators(i,1)   ! holes for second single exc alpha 
+   generators_bitmask(i,1,6,k) = particles_operators(i,1) ! particles for second single exc alpha 
+   generators_bitmask(i,2,5,k) = holes_operators(i,2)   ! holes for second single exc beta 
+   generators_bitmask(i,2,6,k) = particles_operators(i,2) ! particles for second single exc beta 
+
+  enddo
+ enddo
+ touch generators_bitmask
+
+
+
+end