Working on MRCC

src/MRCC/H_apply.irp.f

@@ -0,0 +1,10 @@
+use bitmasks
+BEGIN_SHELL [ /usr/bin/env python ]
+from generate_h_apply import *
+
+s = H_apply("mrcc")
+s.data["keys_work"] = "call mrcc_dress(i_generator,key_idx,keys_out,N_int,iproc)"
+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 @@
+AOs BiInts Bitmask CAS_SD_selected Dets Electrons Ezfio_files Generators_CAS Hartree_Fock MOGuess MonoInts MOs Nuclei Output Perturbation Properties Selectors_full Utils

src/MRCC/README.rst

@@ -0,0 +1,4 @@
+=======
+ Module
+=======
+

src/MRCC/cas_sd.irp.f

@@ -0,0 +1,57 @@
+program full_ci
+  implicit none
+  integer                        :: i,k
+
+  
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
+  character*(64)                 :: perturbation
+  
+  pt2 = 1.d0
+  diag_algorithm = "Lapack"
+  if (N_det > n_det_max_fci) then
+    call diagonalize_CI
+    call save_wavefunction
+    psi_det = psi_det_sorted
+    psi_coef = psi_coef_sorted
+    N_det = n_det_max_fci
+    soft_touch N_det psi_det psi_coef
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+  endif
+
+  do while (N_det < n_det_max_fci.and.maxval(abs(pt2(1:N_st))) > pt2_max)
+    call H_apply_FCI(pt2, norm_pert, H_pert_diag,  N_st)
+
+    PROVIDE  psi_coef
+    PROVIDE  psi_det
+    PROVIDE  psi_det_sorted
+
+    if (N_det > n_det_max_fci) then
+       psi_det = psi_det_sorted
+       psi_coef = psi_coef_sorted
+       N_det = n_det_max_fci
+       soft_touch N_det psi_det psi_coef
+    endif
+    call diagonalize_CI
+    call save_wavefunction
+    print *,  'N_det    = ', N_det
+    print *,  'N_states = ', N_states
+    print *,  'PT2      = ', pt2
+    print *,  'E        = ', CI_energy
+    print *,  'E+PT2    = ', CI_energy+pt2
+    print *,  '-----'
+    call ezfio_set_full_ci_energy(CI_energy)
+    if (abort_all) then
+      exit
+    endif
+  enddo
+end

src/MRCC/mrcc.irp.f

@@ -0,0 +1,39 @@
+program mrcc
+  implicit none
+  read_wf = .True.
+  TOUCH read_wf
+
+
+  print *, N_det
+  print *, N_det_cas
+  print *, N_det_sd
+!  psi_cas, (N_int,2,N_det_generators) ]
+!psi_cas_coefs,  (N_det_generators,n_states) ]
+!psi_sd,  (N_int,2,psi_det_size) ]
+!psi_sd_coefs, (psi_det_size,n_states) ]
+
+  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
+
+  double precision, allocatable  :: pt2(:), norm_pert(:), H_pert_diag(:)
+  integer                        :: N_st, degree
+  N_st = N_states
+  allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
+
+  print *,  'MRCC'
+  print *,  '===='
+  call H_apply_mrcc(pt2, norm_pert, H_pert_diag,  N_st)
+end

src/MRCC/mrcc_dress.irp.f

@@ -0,0 +1,54 @@
+subroutine mrcc_dress(i_generator,n_selected,det_buffer,Nint,iproc)
+ use bitmasks
+ implicit none
+
+  integer, intent(in)            :: i_generator,n_selected, Nint, iproc
+  integer(bit_kind), intent(in)  :: det_buffer(Nint,2,n_selected)
+  integer                        :: i,j,k
+  integer                        :: new_size
+  logical                        :: is_in_wavefunction
+  double precision               :: degree(N_det_cas)
+  integer                        :: idx(0:N_det_cas)
+  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_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
+
+  print *, N_tq
+  do i=1,N_tq
+    call debug_det(det_buffer(1,1,i),Nint)
+  enddo
+end
+

src/MRCC/mrcc_utils.irp.f

@@ -0,0 +1,77 @@
+ use bitmasks
+ BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,N_det_generators) ]
+&BEGIN_PROVIDER [ double precision, psi_cas_coefs,  (N_det_generators,n_states) ]
+&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_cas, (N_det_generators) ]
+&BEGIN_PROVIDER [ integer, idx_sd,  (N_det) ]
+&BEGIN_PROVIDER [ integer, N_det_sd]
+&BEGIN_PROVIDER [ integer, N_det_cas]
+ implicit none
+ BEGIN_DOC
+ ! SD
+ END_DOC
+ integer                        :: i_cas,i_sd,j,k
+ integer                        :: degree
+ logical                        :: in_cas
+ i_cas=0
+ i_sd =0
+ do k=1,N_det
+   in_cas = .False.
+   do j=1,n_det_generators
+     call get_excitation_degree(psi_generators(1,1,j), psi_det(1,1,k), degree, N_int)
+     if (degree == 0) then
+       i_cas  += 1
+       psi_cas(1:N_int,1:2,i_cas) = psi_det(1:N_int,1:2,k)
+       psi_cas_coefs(i_cas,1:N_states) = psi_coef(k,1:N_states)
+       in_cas = .True.
+       idx_cas(i_cas) = k
+       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
+ N_det_cas = i_cas
+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
+
+subroutine update_generators
+  implicit none
+  integer :: i,j,k
+  n_det_generators = N_det_sd
+  do k=1,N_det_sd
+    do j=1,2
+      do i=1,N_int
+        psi_generators(i,j,k) = psi_sd(i,j,k)
+      enddo
+    enddo
+  enddo
+end