diff --git a/config/gfortran.cfg b/config/gfortran.cfg index 12cd6e0c..86f625c6 100644 --- a/config/gfortran.cfg +++ b/config/gfortran.cfg @@ -10,7 +10,7 @@ # # [COMMON] -FC : gfortran -ffree-line-length-none -I . +FC : gfortran -ffree-line-length-none -I . -mavx LAPACK_LIB : -llapack -lblas IRPF90 : irpf90 IRPF90_FLAGS : --ninja --align=32 @@ -35,7 +35,7 @@ OPENMP : 1 ; Append OpenMP flags # -ffast-math and the Fortran-specific # -fno-protect-parens and -fstack-arrays. [OPT] -FCFLAGS : -Ofast +FCFLAGS : -Ofast -g # Profiling flags ################# diff --git a/ocaml/.gitignore b/ocaml/.gitignore index 235419de..0cf01652 100644 --- a/ocaml/.gitignore +++ b/ocaml/.gitignore @@ -3,45 +3,44 @@ ezfio.ml Qptypes.ml qptypes_generator.byte _build -qp_basis_clean.native qp_create_ezfio_from_xyz.native -qp_edit.native -qp_print.native -qp_run.native qp_set_ddci.native -qp_set_mo_class.native +qp_print.native +qp_edit.native +qp_set_mo_class.native +qp_basis_clean.native +qp_run.native qp_edit.native -test_atom.byte -test_basis.byte -test_bitlist.byte -test_determinants.byte -test_elements.byte -test_excitation.byte -test_gto.byte test_mo_label.byte -test_molecule.byte test_point3d.byte -test_atom -test_basis -test_bitlist -test_determinants -test_elements -test_excitation -test_gto +test_gto.byte +test_excitation.byte +test_determinants.byte +test_basis.byte +test_molecule.byte +test_elements.byte +test_bitlist.byte +test_atom.byte test_mo_label -test_molecule test_point3d -qp_basis_clean +test_gto +test_excitation +test_determinants +test_basis +test_molecule +test_elements +test_bitlist +test_atom qp_create_ezfio_from_xyz -qp_edit -qp_print -qp_run qp_set_ddci +qp_print +qp_edit qp_set_mo_class +qp_basis_clean +qp_run Input_determinants.ml -Input_properties.ml -Input_hartree_fock.ml Input_integrals_bielec.ml +Input_pseudo.ml Input_perturbation.ml Input_pseudo.ml qp_edit.ml diff --git a/plugins/Full_CI/full_ci.irp.f b/plugins/Full_CI/full_ci.irp.f index c7b8452d..074ec7e1 100644 --- a/plugins/Full_CI/full_ci.irp.f +++ b/plugins/Full_CI/full_ci.irp.f @@ -11,6 +11,7 @@ program full_ci pt2 = 1.d0 diag_algorithm = "Lapack" + if (N_det > N_det_max) then call diagonalize_CI call save_wavefunction diff --git a/plugins/MRCC_Utils/H_apply.irp.f b/plugins/MRCC_Utils/H_apply.irp.f index 455c37da..7314674c 100644 --- a/plugins/MRCC_Utils/H_apply.irp.f +++ b/plugins/MRCC_Utils/H_apply.irp.f @@ -9,7 +9,7 @@ s.data["declarations"] += """ double precision, intent(in) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) double precision, intent(in) :: delta_ii_(Ndet_ref,*) """ -s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc)" +s.data["keys_work"] = "call mrcc_dress(delta_ij_,delta_ii_,Ndet_ref,Ndet_non_ref,i_generator,key_idx,keys_out,N_int,iproc,key_mask)" s.data["params_post"] += ", delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" s.data["params_main"] += "delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref" s.data["decls_main"] += """ diff --git a/plugins/MRCC_Utils/mrcc_dress.irp.f b/plugins/MRCC_Utils/mrcc_dress.irp.f index 8c86f7fa..3a4717ed 100644 --- a/plugins/MRCC_Utils/mrcc_dress.irp.f +++ b/plugins/MRCC_Utils/mrcc_dress.irp.f @@ -12,12 +12,61 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), psi_ref_lock, (psi_det_size) ] END_PROVIDER -subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc) + +subroutine create_minilist(key_mask, fullList, miniList, idx_miniList, N_fullList, N_miniList, Nint) + use bitmasks + implicit none + + integer(bit_kind), intent(in) :: fullList(Nint, 2, N_fullList) + integer, intent(in) :: N_fullList + integer(bit_kind),intent(out) :: miniList(Nint, 2, N_fullList) + integer,intent(out) :: idx_miniList(N_fullList), N_miniList + integer, intent(in) :: Nint + integer(bit_kind) :: key_mask(Nint, 2) + integer :: ni, i, n_a, n_b, e_a, e_b + + + n_a = 0 + n_b = 0 + do ni=1,nint + n_a = n_a + popcnt(key_mask(ni,1)) + n_b = n_b + popcnt(key_mask(ni,2)) + end do + + if(n_a == 0) then + N_miniList = N_fullList + miniList(:,:,:) = fullList(:,:,:) + do i=1,N_fullList + idx_miniList(i) = i + end do + return + end if + + N_miniList = 0 + + do i=1,N_fullList + e_a = n_a + e_b = n_b + do ni=1,nint + e_a -= popcnt(iand(fullList(ni, 1, i), key_mask(ni, 1))) + e_b -= popcnt(iand(fullList(ni, 2, i), key_mask(ni, 2))) + end do + + if(e_a + e_b <= 2) then + N_miniList = N_miniList + 1 + miniList(:,:,N_miniList) = fullList(:,:,i) + idx_miniList(N_miniList) = i + end if + end do +end subroutine + + +subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n_selected,det_buffer,Nint,iproc,key_mask) use bitmasks implicit none integer, intent(in) :: i_generator,n_selected, Nint, iproc - integer, intent(in) :: Ndet_ref, Ndet_non_ref + integer, intent(in) :: Ndet_ref, Ndet_non_ref double precision, intent(inout) :: delta_ij_(Ndet_ref,Ndet_non_ref,*) double precision, intent(inout) :: delta_ii_(Ndet_ref,*) @@ -40,7 +89,12 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n integer(bit_kind) :: tmp_det(Nint,2) integer :: iint, ipos integer :: i_state, k_sd, l_sd, i_I, i_alpha - + + integer(bit_kind) :: miniList(Nint, 2, N_det_non_ref), key_mask(Nint, 2) + integer :: idx_miniList(N_det_non_ref), N_miniList + + + call find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq,N_tq) allocate (dIa_hla(N_states,Ndet_non_ref)) @@ -48,9 +102,20 @@ subroutine mrcc_dress(delta_ij_, delta_ii_, Ndet_ref, Ndet_non_ref,i_generator,n ! |I> ! |alpha> + + if(N_tq > 0) then + call create_minilist(key_mask, psi_non_ref, miniList, idx_miniList, N_det_non_ref, N_minilist, Nint) + end if + + do i_alpha=1,N_tq - call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) - +! call get_excitation_degree_vector(psi_non_ref,tq(1,1,i_alpha),degree_alpha,Nint,N_det_non_ref,idx_alpha) + call get_excitation_degree_vector(miniList,tq(1,1,i_alpha),degree_alpha,Nint,N_minilist,idx_alpha) + + do j=1,idx_alpha(0) + idx_alpha(j) = idx_miniList(idx_alpha(j)) + end do + ! |I> do i_I=1,N_det_ref ! Find triples and quadruple grand parents @@ -233,7 +298,8 @@ subroutine find_triples_and_quadruples(i_generator,n_selected,det_buffer,Nint,tq ! Select determinants that are triple or quadruple excitations ! from the ref good = .True. - call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) + call get_excitation_degree_vector(psi_ref,det_buffer(1,1,i),degree,Nint,N_det_ref,idx) + !good=(idx(0) == 0) tant que degree > 2 pas retourné par get_excitation_degree_vector do k=1,idx(0) if (degree(k) < 3) then good = .False. diff --git a/src/Determinants/H_apply.template.f b/src/Determinants/H_apply.template.f index 3a05ee0d..00d82996 100644 --- a/src/Determinants/H_apply.template.f +++ b/src/Determinants/H_apply.template.f @@ -1,4 +1,116 @@ + + subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters ) + + integer(bit_kind), intent(in) :: key_in(N_int, 2), hole_1(N_int, 2), hole_2(N_int, 2) + integer(bit_kind), intent(in) :: particl_1(N_int, 2), particl_2(N_int, 2) + integer(bit_kind) :: p1_mask(N_int, 2), p2_mask(N_int, 2) + integer,intent(in) :: i_generator,iproc_in + integer(bit_kind) :: status(N_int*bit_kind_size, 2) + integer :: highest, p1,p2,sp,ni,i,mi + $declarations + + + highest = 0 + status(:,:) = 0 + do sp=1,2 + do ni=1,N_int + do i=1,bit_kind_size + if(iand(1,ishft(key_in(ni, sp), -(i-1))) == 0) then + cycle + end if + mi = (ni-1)*bit_kind_size+i + status(mi, sp) = iand(1,ishft(hole_1(ni, sp), -(i-1))) + status(mi, sp) = status(mi, sp) + 2*iand(1,ishft(hole_2(ni, sp), -(i-1))) + if(status(mi, sp) /= 0 .and. mi > highest) then + highest = mi + end if + end do + end do + end do + + do sp=1,2 + do p1=1,highest + if(status(p1, sp) == 0) then + cycle + end if + do p2=1,highest + if(status(p2, sp) == 0) then + cycle + end if + if((status(p1, sp) == 1 .and. status(p2, sp) > 1) .or. & + (status(p1, sp) == 2 .and. status(p2, sp) == 3) .or. & + (status(p1, sp) == 3 .and. status(p2, sp) == 3 .and. p2 > p1)) then + call $subroutine_diexcP(key_in, sp, p1, particl_1, sp, p2, particl_2, i_generator, iproc_in $parameters ) + end if + end do + end do + end do + do p1=1,highest + if(status(p1, 1) == 0) then + cycle + end if + do p2=1,highest + if(status(p2, 2) == 0) then + cycle + end if + if((status(p1, 1) == 3) .or. & + (status(p1, 1) == 1 .and. status(p2, 2) >= 2) .or. & + (status(p1, 1) == 2 .and. status(p2, 2) /= 2)) then + + call $subroutine_diexcP(key_in, 1, p1, particl_1, 2, p2, particl_2, i_generator, iproc_in $parameters ) + end if + end do + end do +end subroutine + + +subroutine $subroutine_diexcP(key_in, fs1, fh1, particl_1, fs2, fh2, particl_2, i_generator, iproc_in $parameters ) + + integer(bit_kind), intent(in) :: key_in(N_int, 2), particl_1(N_int, 2), particl_2(N_int, 2) + integer(bit_kind) :: p1_mask(N_int, 2), p2_mask(N_int, 2), key_mask(N_int, 2) + integer,intent(in) :: fh1,fh2,fs1,fs2,i_generator,iproc_in + integer(bit_kind) :: miniList(N_int, 2, N_det) + integer :: n_minilist, n_alpha, n_beta, deg(2), i, ni + $declarations + + p1_mask(:,:) = 0 + p2_mask(:,:) = 0 + p1_mask(fh1/bit_kind_size + 1, fs1) = 2**(mod(fh1-1,bit_kind_size)) + p2_mask(fh2/bit_kind_size + 1, fs2) = 2**(mod(fh2-1,bit_kind_size)) + +! n_alpha = 0 +! n_beta = 0 + key_mask(:,:) = key_in(:,:) + key_mask(fh1/bit_kind_size + 1, fs1) -= 2**(mod(fh1-1,bit_kind_size)) + key_mask(fh2/bit_kind_size + 1, fs2) -= 2**(mod(fh2-1,bit_kind_size)) +! +! do i=1,N_int +! n_alpha = n_alpha + popcnt(key_mask(i, 1)) +! n_beta = n_beta + popcnt(key_mask(i, 2)) +! end do +! +! do i=1, N_det +! deg(1) = n_alpha +! deg(2) = n_beta +! +! do ni = 1, N_int +! ! deg(1) = deg(1) - popcnt(iand(key_mask(ni, 1), psi_non_ref(ni, 1, i))) +! ! deg(2) = deg(2) - popcnt(iand(key_mask(ni, 2), psi_non_ref(ni, 2, i))) +! end do +! +! +! if(deg(1) + deg(2) <= 2) then +! ! ndet_out = ndet_out + 1 +! ! idx(ndet_out) = i +! end if +! end do + + call $subroutine_diexcOrg(key_in, key_mask, p1_mask, particl_1, p2_mask, particl_2, i_generator, iproc_in $parameters ) +end subroutine + + +subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl_2, i_generator, iproc_in $parameters ) use omp_lib use bitmasks implicit none @@ -10,7 +122,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene integer,parameter :: size_max = $size_max $declarations integer ,intent(in) :: i_generator - integer(bit_kind),intent(in) :: key_in(N_int,2) + integer(bit_kind),intent(in) :: key_in(N_int,2), key_mask(N_int, 2) integer(bit_kind),allocatable :: keys_out(:,:,:) integer(bit_kind), intent(in) :: hole_1(N_int,2), particl_1(N_int,2) integer(bit_kind), intent(in) :: hole_2(N_int,2), particl_2(N_int,2) @@ -290,13 +402,18 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa double precision :: diag_H_mat_elem integer(omp_lock_kind), save :: lck, ifirst=0 integer :: iproc - + + integer(bit_kind) :: key_mask(N_int, 2) + logical :: check_double_excitation + + key_mask(:,:) = 0_8 + iproc = iproc_in check_double_excitation = .True. $check_double_excitation - + if (ifirst == 0) then ifirst=1 diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f index 49a3604a..8f594738 100644 --- a/src/Determinants/connected_to_ref.irp.f +++ b/src/Determinants/connected_to_ref.irp.f @@ -169,7 +169,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) ! output : 0 : not connected ! i : connected to determinant i of the past - ! -i : is the ith determinant of the refernce wf keys + ! -i : is the ith determinant of the reference wf keys ASSERT (Nint > 0) ASSERT (Nint == N_int) diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f index be847472..edda4362 100644 --- a/src/Determinants/slater_rules.irp.f +++ b/src/Determinants/slater_rules.irp.f @@ -787,6 +787,7 @@ subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) ASSERT (Ndet > 0) ASSERT (Ndet_max >= Ndet) i_H_psi_array = 0.d0 + call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) do ii=1,idx(0) i = idx(ii)