diff --git a/scripts/generate_h_apply.py b/scripts/generate_h_apply.py index 36bd48aa..ad6a57ee 100755 --- a/scripts/generate_h_apply.py +++ b/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) diff --git a/src/CAS_SD_selected/cas_sd.irp.f b/src/CAS_SD_selected/cas_sd.irp.f index e3a8f6ec..e5e0497a 100644 --- a/src/CAS_SD_selected/cas_sd.irp.f +++ b/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 diff --git a/src/Dets/H_apply_template.f b/src/Dets/H_apply_template.f index f4486c36..e44562e7 100644 --- a/src/Dets/H_apply_template.f +++ b/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 diff --git a/src/Dets/connected_to_ref.irp.f b/src/Dets/connected_to_ref.irp.f index fd87cc3c..07331782 100644 --- a/src/Dets/connected_to_ref.irp.f +++ b/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 - endif + +! if (i > N_det) then +! 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 -! degree = 1 -! do i=1,N_det -! integer :: degree -! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) -! if (degree == 0) then -! exit -! endif -! enddo -! if (degree /=0) then -! stop 'pouet 1' -! endif -! else -! do i=1,N_det -! call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) -! if (degree == 0) then -! stop 'pouet 2' -! endif -! enddo -! endif + if (is_in_wavefunction) then + degree = 1 + do i=1,N_det + integer :: degree + call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) + if (degree == 0) then + exit + endif + enddo + if (degree /=0) then + stop 'pouet 1' + endif + else + do i=1,N_det + call get_excitation_degree(key,psi_det(1,1,i),degree,N_int) + if (degree == 0) then + stop 'pouet 2' + endif + enddo + endif ! END DEBUG is_in_wf end diff --git a/src/Generators_CAS/generators.irp.f b/src/Generators_CAS/generators.irp.f index 1ea98eee..5dcfc1b0 100644 --- a/src/Generators_CAS/generators.irp.f +++ b/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_generators(k,2,m) = psi_det(k,2,i) + psi_det_generators(k,1,m) = psi_det(k,1,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 diff --git a/src/Generators_full/generators.irp.f b/src/Generators_full/generators.irp.f index 58fefef7..4d261acd 100644 --- a/src/Generators_full/generators.irp.f +++ b/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_generators(k,2,i) = psi_det_sorted(k,2,i) + psi_det_generators(k,1,i) = psi_det_sorted(k,1,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 diff --git a/src/Generators_restart/generators.irp.f b/src/Generators_restart/generators.irp.f index 19a6a9fd..2e0bc375 100644 --- a/src/Generators_restart/generators.irp.f +++ b/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 [ double precision, psi_generators_coef, (N_det_generators,N_states) + 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 ! 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_generators(k,2,i) = psi_det(k,2,i) + psi_det_generators(k,1,i) = psi_det(k,1,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 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..4f284112 --- /dev/null +++ b/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 + 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..3711786f --- /dev/null +++ b/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 + diff --git a/src/MRCC/README.rst b/src/MRCC/README.rst new file mode 100644 index 00000000..1fbcccdb --- /dev/null +++ b/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 + diff --git a/src/MRCC/mrcc.irp.f b/src/MRCC/mrcc.irp.f new file mode 100644 index 00000000..663b9f5a --- /dev/null +++ b/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 diff --git a/src/MRCC/mrcc_dress.irp.f b/src/MRCC/mrcc_dress.irp.f new file mode 100644 index 00000000..e5c8d233 --- /dev/null +++ b/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 / (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 + diff --git a/src/MRCC/mrcc_utils.irp.f b/src/MRCC/mrcc_utils.irp.f new file mode 100644 index 00000000..2c5dc811 --- /dev/null +++ b/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/ + 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 + diff --git a/src/Makefile.common b/src/Makefile.common index 9ecc93c0..eddaa481 100644 --- a/src/Makefile.common +++ b/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 diff --git a/src/NEEDED_MODULES b/src/NEEDED_MODULES index f71f28fd..a144c42c 100644 --- a/src/NEEDED_MODULES +++ b/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 diff --git a/src/Perturbation/perturbation_template.f b/src/Perturbation/perturbation_template.f index 3627dd17..9f50b636 100644 --- a/src/Perturbation/perturbation_template.f +++ b/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 diff --git a/src/SingleRefMethod/generators.irp.f b/src/SingleRefMethod/generators.irp.f index 6ef30fa1..ce71f996 100644 --- a/src/SingleRefMethod/generators.irp.f +++ b/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_generators(i,2,1) = HF_bitmask(i,2) + psi_det_generators(i,1,1) = HF_bitmask(i,1) + 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_generators(:,:,j) = psi_det(:,:,j-1) - enddo - do j=k+1,N_det - psi_generators(:,:,j) = psi_det(:,:,j) - enddo + psi_det_generators(:,:,1) = psi_det(:,:,j) + psi_coef_generators(1,:) = psi_coef_generators(j,:) END_PROVIDER