diff --git a/src/MRCC/ASSUMPTIONS.rst b/src/MRCC/ASSUMPTIONS.rst new file mode 100644 index 00000000..e69de29b diff --git a/src/MRCC/H_apply.irp.f b/src/MRCC/H_apply.irp.f new file mode 100644 index 00000000..757627e1 --- /dev/null +++ b/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 + diff --git a/src/MRCC/Makefile b/src/MRCC/Makefile new file mode 100644 index 00000000..06dc50ff --- /dev/null +++ b/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 diff --git a/src/MRCC/NEEDED_MODULES b/src/MRCC/NEEDED_MODULES new file mode 100644 index 00000000..31a27888 --- /dev/null +++ b/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 diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst new file mode 100644 index 00000000..ebb762c1 --- /dev/null +++ b/src/MRCC/README.rst @@ -0,0 +1,4 @@ +======= + Module +======= + diff --git a/src/MRCC/cas_sd.irp.f b/src/MRCC/cas_sd.irp.f new file mode 100644 index 00000000..e3a8f6ec --- /dev/null +++ b/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 diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f new file mode 100644 index 00000000..4f2a6858 --- /dev/null +++ b/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 diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f new file mode 100644 index 00000000..86ffc9f5 --- /dev/null +++ b/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 + diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f new file mode 100644 index 00000000..327c220e --- /dev/null +++ b/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/ + 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