From bf997c558372b02ae1b87fd78c3f628675060697 Mon Sep 17 00:00:00 2001 From: Thomas Applencourt Date: Mon, 20 Apr 2015 16:45:06 +0200 Subject: [PATCH] Add all the mising file... --- ocaml/Input_determinants.ml | 251 ---- scripts/ezfio_interface/ei_handler.py | 1 - src/Determinants/ASSUMPTIONS.rst | 7 + src/Determinants/EZFIO.cfg | 103 ++ src/Determinants/H_apply.irp.f | 229 +++ src/Determinants/H_apply_template.f | 542 +++++++ src/Determinants/Makefile | 6 + src/Determinants/NEEDED_MODULES | 1 + src/Determinants/README.rst | 696 +++++++++ src/Determinants/SC2.irp.f | 215 +++ src/Determinants/connected_to_ref.irp.f | 357 +++++ src/Determinants/create_excitations.irp.f | 36 + src/Determinants/davidson.irp.f | 418 ++++++ src/Determinants/density_matrix.irp.f | 214 +++ src/Determinants/det_svd.irp.f | 61 + src/Determinants/determinants.irp.f | 9 - src/Determinants/determinants_bitmasks.irp.f | 57 + src/Determinants/diagonalize_CI.irp.f | 109 ++ src/Determinants/diagonalize_CI_SC2.irp.f | 59 + src/Determinants/diagonalize_CI_mono.irp.f | 72 + src/Determinants/excitations_utils.irp.f | 16 + src/Determinants/filter_connected.irp.f | 611 ++++++++ src/Determinants/guess_doublet.irp.f | 79 + src/Determinants/guess_singlet.irp.f | 44 + src/Determinants/guess_triplet.irp.f | 48 + src/Determinants/occ_pattern.irp.f | 339 +++++ src/Determinants/options.irp.f | 22 + .../program_beginer_determinants.irp.f | 138 ++ src/Determinants/psi_cas.irp.f | 114 ++ src/Determinants/ref_bitmask.irp.f | 57 + src/Determinants/s2.irp.f | 106 ++ src/Determinants/save_for_casino.irp.f | 268 ++++ src/Determinants/save_for_qmcchem.irp.f | 51 + src/Determinants/save_natorb.irp.f | 6 + src/Determinants/slater_rules.irp.f | 1301 +++++++++++++++++ .../spindeterminants.ezfio_config | 17 + src/Determinants/spindeterminants.irp.f | 615 ++++++++ src/Determinants/truncate_wf.irp.f | 18 + src/Determinants/utils.irp.f | 20 + src/Output/README.rst | 1 + src/Properties/EZFIO.cfg | 5 + 41 files changed, 7058 insertions(+), 261 deletions(-) delete mode 100644 ocaml/Input_determinants.ml create mode 100644 src/Determinants/ASSUMPTIONS.rst create mode 100644 src/Determinants/EZFIO.cfg create mode 100644 src/Determinants/H_apply.irp.f create mode 100644 src/Determinants/H_apply_template.f create mode 100644 src/Determinants/Makefile create mode 100644 src/Determinants/NEEDED_MODULES create mode 100644 src/Determinants/README.rst create mode 100644 src/Determinants/SC2.irp.f create mode 100644 src/Determinants/connected_to_ref.irp.f create mode 100644 src/Determinants/create_excitations.irp.f create mode 100644 src/Determinants/davidson.irp.f create mode 100644 src/Determinants/density_matrix.irp.f create mode 100644 src/Determinants/det_svd.irp.f create mode 100644 src/Determinants/determinants_bitmasks.irp.f create mode 100644 src/Determinants/diagonalize_CI.irp.f create mode 100644 src/Determinants/diagonalize_CI_SC2.irp.f create mode 100644 src/Determinants/diagonalize_CI_mono.irp.f create mode 100644 src/Determinants/excitations_utils.irp.f create mode 100644 src/Determinants/filter_connected.irp.f create mode 100644 src/Determinants/guess_doublet.irp.f create mode 100644 src/Determinants/guess_singlet.irp.f create mode 100644 src/Determinants/guess_triplet.irp.f create mode 100644 src/Determinants/occ_pattern.irp.f create mode 100644 src/Determinants/options.irp.f create mode 100644 src/Determinants/program_beginer_determinants.irp.f create mode 100644 src/Determinants/psi_cas.irp.f create mode 100644 src/Determinants/ref_bitmask.irp.f create mode 100644 src/Determinants/s2.irp.f create mode 100644 src/Determinants/save_for_casino.irp.f create mode 100644 src/Determinants/save_for_qmcchem.irp.f create mode 100644 src/Determinants/save_natorb.irp.f create mode 100644 src/Determinants/slater_rules.irp.f create mode 100644 src/Determinants/spindeterminants.ezfio_config create mode 100644 src/Determinants/spindeterminants.irp.f create mode 100644 src/Determinants/truncate_wf.irp.f create mode 100644 src/Determinants/utils.irp.f create mode 100644 src/Properties/EZFIO.cfg diff --git a/ocaml/Input_determinants.ml b/ocaml/Input_determinants.ml deleted file mode 100644 index df046231..00000000 --- a/ocaml/Input_determinants.ml +++ /dev/null @@ -1,251 +0,0 @@ -(* =~=~ *) -(* Init *) -(* =~=~ *) - -open Qptypes;; -open Qputils;; -open Core.Std;; - -module Determinants : sig -(* Generate type *) - type t = - { - n_det_max_jacobi : Strictly_positive_int.t; - threshold_generators : Threshold.t; - threshold_selectors : Threshold.t; - n_states : Strictly_positive_int.t; - s2_eig : bool; - read_wf : bool; - only_single_double_dm : bool; - } with sexp - ;; - val read : unit -> t option - val write : t-> unit - val to_string : t -> string - val to_rst : t -> Rst_string.t - val of_rst : Rst_string.t -> t option -end = struct -(* Generate type *) - type t = - { - n_det_max_jacobi : Strictly_positive_int.t; - threshold_generators : Threshold.t; - threshold_selectors : Threshold.t; - n_states : Strictly_positive_int.t; - s2_eig : bool; - read_wf : bool; - only_single_double_dm : bool; - } with sexp - ;; - - let get_default = Qpackage.get_ezfio_default "determinants";; - -(* =~=~=~=~=~=~==~=~=~=~=~=~ *) -(* Generate Special Function *) -(* =~=~=~==~=~~=~=~=~=~=~=~=~ *) - -(* Read snippet for n_det_max_jacobi *) - let read_n_det_max_jacobi () = - if not (Ezfio.has_determinants_n_det_max_jacobi ()) then - get_default "n_det_max_jacobi" - |> Int.of_string - |> Ezfio.set_determinants_n_det_max_jacobi - ; - Ezfio.get_determinants_n_det_max_jacobi () - |> Strictly_positive_int.of_int - ;; -(* Write snippet for n_det_max_jacobi *) - let write_n_det_max_jacobi var = - Strictly_positive_int.to_int var - |> Ezfio.set_determinants_n_det_max_jacobi - ;; - -(* Read snippet for n_states *) - let read_n_states () = - if not (Ezfio.has_determinants_n_states ()) then - get_default "n_states" - |> Int.of_string - |> Ezfio.set_determinants_n_states - ; - Ezfio.get_determinants_n_states () - |> Strictly_positive_int.of_int - ;; -(* Write snippet for n_states *) - let write_n_states var = - Strictly_positive_int.to_int var - |> Ezfio.set_determinants_n_states - ;; - -(* Read snippet for only_single_double_dm *) - let read_only_single_double_dm () = - if not (Ezfio.has_determinants_only_single_double_dm ()) then - get_default "only_single_double_dm" - |> Bool.of_string - |> Ezfio.set_determinants_only_single_double_dm - ; - Ezfio.get_determinants_only_single_double_dm () - ;; -(* Write snippet for only_single_double_dm *) - let write_only_single_double_dm = - Ezfio.set_determinants_only_single_double_dm - ;; - -(* Read snippet for read_wf *) - let read_read_wf () = - if not (Ezfio.has_determinants_read_wf ()) then - get_default "read_wf" - |> Bool.of_string - |> Ezfio.set_determinants_read_wf - ; - Ezfio.get_determinants_read_wf () - ;; -(* Write snippet for read_wf *) - let write_read_wf = - Ezfio.set_determinants_read_wf - ;; - -(* Read snippet for s2_eig *) - let read_s2_eig () = - if not (Ezfio.has_determinants_s2_eig ()) then - get_default "s2_eig" - |> Bool.of_string - |> Ezfio.set_determinants_s2_eig - ; - Ezfio.get_determinants_s2_eig () - ;; -(* Write snippet for s2_eig *) - let write_s2_eig = - Ezfio.set_determinants_s2_eig - ;; - -(* Read snippet for threshold_generators *) - let read_threshold_generators () = - if not (Ezfio.has_determinants_threshold_generators ()) then - get_default "threshold_generators" - |> Float.of_string - |> Ezfio.set_determinants_threshold_generators - ; - Ezfio.get_determinants_threshold_generators () - |> Threshold.of_float - ;; -(* Write snippet for threshold_generators *) - let write_threshold_generators var = - Threshold.to_float var - |> Ezfio.set_determinants_threshold_generators - ;; - -(* Read snippet for threshold_selectors *) - let read_threshold_selectors () = - if not (Ezfio.has_determinants_threshold_selectors ()) then - get_default "threshold_selectors" - |> Float.of_string - |> Ezfio.set_determinants_threshold_selectors - ; - Ezfio.get_determinants_threshold_selectors () - |> Threshold.of_float - ;; -(* Write snippet for threshold_selectors *) - let write_threshold_selectors var = - Threshold.to_float var - |> Ezfio.set_determinants_threshold_selectors - ;; - -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) -(* Generate Global Function *) -(* =~=~=~=~=~=~=~=~=~=~=~=~ *) - -(* Read all *) - let read() = - Some - { - n_det_max_jacobi = read_n_det_max_jacobi (); - threshold_generators = read_threshold_generators (); - threshold_selectors = read_threshold_selectors (); - n_states = read_n_states (); - s2_eig = read_s2_eig (); - read_wf = read_read_wf (); - only_single_double_dm = read_only_single_double_dm (); - } - ;; -(* Write all *) - let write{ - n_det_max_jacobi; - threshold_generators; - threshold_selectors; - n_states; - s2_eig; - read_wf; - only_single_double_dm; - } = - write_n_det_max_jacobi n_det_max_jacobi; - write_threshold_generators threshold_generators; - write_threshold_selectors threshold_selectors; - write_n_states n_states; - write_s2_eig s2_eig; - write_read_wf read_wf; - write_only_single_double_dm only_single_double_dm; - ;; -(* to_string*) - let to_string b = - Printf.sprintf " - n_det_max_jacobi = %s - threshold_generators = %s - threshold_selectors = %s - n_states = %s - s2_eig = %s - read_wf = %s - only_single_double_dm = %s - " - (Strictly_positive_int.to_string b.n_det_max_jacobi) - (Threshold.to_string b.threshold_generators) - (Threshold.to_string b.threshold_selectors) - (Strictly_positive_int.to_string b.n_states) - (Bool.to_string b.s2_eig) - (Bool.to_string b.read_wf) - (Bool.to_string b.only_single_double_dm) - ;; -(* to_rst*) - let to_rst b = - Printf.sprintf " - Maximum number of determinants diagonalized by Jacobi :: - - n_det_max_jacobi = %s - - Percentage of the norm of the state-averaged wave function to consider for the generators :: - - threshold_generators = %s - - Percentage of the norm of the state-averaged wave function to consider for the selectors :: - - threshold_selectors = %s - - Number of states to consider :: - - n_states = %s - - Force the wave function to be an eigenfunction of S^2 :: - - s2_eig = %s - - If true, read the wave function from the EZFIO file :: - - read_wf = %s - - If true, The One body DM is calculated with ignoing the Double <-> Doubles extra diag elements :: - - only_single_double_dm = %s - - " - (Strictly_positive_int.to_string b.n_det_max_jacobi) - (Threshold.to_string b.threshold_generators) - (Threshold.to_string b.threshold_selectors) - (Strictly_positive_int.to_string b.n_states) - (Bool.to_string b.s2_eig) - (Bool.to_string b.read_wf) - (Bool.to_string b.only_single_double_dm) - |> Rst_string.of_string - ;; - include Generic_input_of_rst;; - let of_rst = of_rst t_of_sexp;; - -end \ No newline at end of file diff --git a/scripts/ezfio_interface/ei_handler.py b/scripts/ezfio_interface/ei_handler.py index e5c08895..6d18d071 100755 --- a/scripts/ezfio_interface/ei_handler.py +++ b/scripts/ezfio_interface/ei_handler.py @@ -278,7 +278,6 @@ def get_dict_config_file(config_file_path, module_lower): try: d[pvd]["default"] = is_bool(default_raw) - print is_bool(default_raw) except TypeError: d[pvd]["default"] = Type(None, default_raw, default_raw) diff --git a/src/Determinants/ASSUMPTIONS.rst b/src/Determinants/ASSUMPTIONS.rst new file mode 100644 index 00000000..e9e24d09 --- /dev/null +++ b/src/Determinants/ASSUMPTIONS.rst @@ -0,0 +1,7 @@ +* The MOs are orthonormal +* All the determinants have the same number of electrons +* The determinants are orthonormal +* The number of generator determinants <= the number of determinants +* All the determinants in the H_apply buffer are supposed to be different from the + wave function determinants +* All the determinants in the H_apply buffer are supposed to be unique diff --git a/src/Determinants/EZFIO.cfg b/src/Determinants/EZFIO.cfg new file mode 100644 index 00000000..5f63404b --- /dev/null +++ b/src/Determinants/EZFIO.cfg @@ -0,0 +1,103 @@ +[N_states] +type: States_number +doc: Number of states to consider +interface: input +default: 1 + +[N_det_max_jacobi] +type: Strictly_positive_int +doc: Maximum number of determinants diagonalized by Jacobi +interface: input +default: 1000 + +[read_wf] +type: logical +doc: If true, read the wave function from the EZFIO file +interface: input +default: False + +[only_single_double_dm] +type: logical +doc: If true, The One body DM is calculated with ignoring the Double<->Doubles extra diag elements +interface: input +default: False + +[s2_eig] +type: logical +doc: Force the wave function to be an eigenfunction of S^2 +interface: input +default: False + +[threshold_generators] +type: Threshold +doc: Thresholds on generators (fraction of the norm) +interface: input +default: 0.99 + +[threshold_selectors] +type: Threshold +doc: Thresholds on selectors (fraction of the norm) +interface: input +default: 0.999 + + +# Only create the ezfio_config, (no Input_* and no PROVIDER) + +[n_states_diag] +type: integer +doc: n_states_diag +interface: Ocaml + +[n_int] +interface: OCaml +doc: n_int +type: N_int_number + +[bit_kind] +interface: OCaml +doc: bit_kind +type: Bit_kind + +[mo_label] +interface: OCaml +doc: o_label +type: character*(64) + +[n_det] +interface: OCaml +doc: n_det +type: integer + +[psi_coef] +interface: OCaml +doc: psi_coef +type: double precision +size: (determinants_n_det,determinants_n_states) + +[psi_det] +interface: OCaml +doc: psi_det +type: integer*8 +size: (determinants_n_int*determinants_bit_kind/8,2,determinants_n_det) + +[det_num] +interface: OCaml +doc: det_num +type: integer + +[det_occ] +interface: OCaml +doc: det_occ +type: integer +size: (electrons_elec_alpha_num,determinants_det_num,2) + +[det_coef] +interface: OCaml +doc: det_coef +type: double precision +size: (determinants_det_num) + +[expected_s2] +interface: OCaml +doc: expcted_s2 +type: double precision \ No newline at end of file diff --git a/src/Determinants/H_apply.irp.f b/src/Determinants/H_apply.irp.f new file mode 100644 index 00000000..801d00a5 --- /dev/null +++ b/src/Determinants/H_apply.irp.f @@ -0,0 +1,229 @@ +use bitmasks +use omp_lib + +type H_apply_buffer_type +integer :: N_det +integer :: sze +integer(bit_kind), pointer :: det(:,:,:) +double precision , pointer :: coef(:,:) +double precision , pointer :: e2(:,:) +end type H_apply_buffer_type + +type(H_apply_buffer_type), pointer :: H_apply_buffer(:) + + + BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ] +&BEGIN_PROVIDER [ integer(omp_lock_kind), H_apply_buffer_lock, (64,0:nproc-1) ] + use omp_lib + implicit none + BEGIN_DOC + ! Buffer of determinants/coefficients/perturbative energy for H_apply. + ! Uninitialized. Filled by H_apply subroutines. + END_DOC + integer :: iproc, sze + sze = 10000 + if (.not.associated(H_apply_buffer)) then + allocate(H_apply_buffer(0:nproc-1)) + iproc = 0 + !$OMP PARALLEL PRIVATE(iproc) DEFAULT(NONE) & + !$OMP SHARED(H_apply_buffer,N_int,sze,N_states,H_apply_buffer_lock) + !$ iproc = omp_get_thread_num() + H_apply_buffer(iproc)%N_det = 0 + H_apply_buffer(iproc)%sze = sze + allocate ( & + H_apply_buffer(iproc)%det(N_int,2,sze), & + H_apply_buffer(iproc)%coef(sze,N_states), & + H_apply_buffer(iproc)%e2(sze,N_states) & + ) + H_apply_buffer(iproc)%det = 0_bit_kind + H_apply_buffer(iproc)%coef = 0.d0 + H_apply_buffer(iproc)%e2 = 0.d0 + call omp_init_lock(H_apply_buffer_lock(1,iproc)) + !$OMP END PARALLEL + endif + +END_PROVIDER + + +subroutine resize_H_apply_buffer(new_size,iproc) + implicit none + integer, intent(in) :: new_size, iproc + integer(bit_kind), pointer :: buffer_det(:,:,:) + double precision, pointer :: buffer_coef(:,:) + double precision, pointer :: buffer_e2(:,:) + integer :: i,j,k + integer :: Ndet + PROVIDE H_apply_buffer_allocated + + ASSERT (new_size > 0) + ASSERT (iproc >= 0) + ASSERT (iproc < nproc) + + call omp_set_lock(H_apply_buffer_lock(1,iproc)) + allocate ( buffer_det(N_int,2,new_size), & + buffer_coef(new_size,N_states), & + buffer_e2(new_size,N_states) ) + + do i=1,min(new_size,H_apply_buffer(iproc)%N_det) + do k=1,N_int + buffer_det(k,1,i) = H_apply_buffer(iproc)%det(k,1,i) + buffer_det(k,2,i) = H_apply_buffer(iproc)%det(k,2,i) + enddo + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i))) == elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num ) + enddo + deallocate(H_apply_buffer(iproc)%det) + H_apply_buffer(iproc)%det => buffer_det + + do k=1,N_states + do i=1,min(new_size,H_apply_buffer(iproc)%N_det) + buffer_coef(i,k) = H_apply_buffer(iproc)%coef(i,k) + enddo + enddo + deallocate(H_apply_buffer(iproc)%coef) + H_apply_buffer(iproc)%coef => buffer_coef + + do k=1,N_states + do i=1,min(new_size,H_apply_buffer(iproc)%N_det) + buffer_e2(i,k) = H_apply_buffer(iproc)%e2(i,k) + enddo + enddo + deallocate(H_apply_buffer(iproc)%e2) + H_apply_buffer(iproc)%e2 => buffer_e2 + + H_apply_buffer(iproc)%sze = new_size + H_apply_buffer(iproc)%N_det = min(new_size,H_apply_buffer(iproc)%N_det) + call omp_unset_lock(H_apply_buffer_lock(1,iproc)) + +end + +subroutine copy_H_apply_buffer_to_wf + use omp_lib + implicit none + BEGIN_DOC +! Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det +! after calling this function. +! After calling this subroutine, N_det, psi_det and psi_coef need to be touched + END_DOC + integer(bit_kind), allocatable :: buffer_det(:,:,:) + double precision, allocatable :: buffer_coef(:,:) + integer :: i,j,k + integer :: N_det_old + integer :: iproc + + PROVIDE H_apply_buffer_allocated + + ASSERT (N_int > 0) + ASSERT (N_det > 0) + + allocate ( buffer_det(N_int,2,N_det), buffer_coef(N_det,N_states) ) + + do i=1,N_det + do k=1,N_int + ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num) + ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num) + buffer_det(k,1,i) = psi_det(k,1,i) + buffer_det(k,2,i) = psi_det(k,2,i) + enddo + enddo + do k=1,N_states + do i=1,N_det + buffer_coef(i,k) = psi_coef(i,k) + enddo + enddo + + N_det_old = N_det + do j=0,nproc-1 + N_det = N_det + H_apply_buffer(j)%N_det + enddo + + if (psi_det_size < N_det) then + psi_det_size = N_det + TOUCH psi_det_size + endif + do i=1,N_det_old + do k=1,N_int + psi_det(k,1,i) = buffer_det(k,1,i) + psi_det(k,2,i) = buffer_det(k,2,i) + enddo + ASSERT (sum(popcnt(psi_det(:,1,i))) == elec_alpha_num) + ASSERT (sum(popcnt(psi_det(:,2,i))) == elec_beta_num ) + enddo + do k=1,N_states + do i=1,N_det_old + psi_coef(i,k) = buffer_coef(i,k) + enddo + enddo + !$OMP PARALLEL DEFAULT(SHARED) & + !$OMP PRIVATE(j,k,i) FIRSTPRIVATE(N_det_old) & + !$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_coef,N_states) + j=0 + !$ j=omp_get_thread_num() + do k=0,j-1 + N_det_old += H_apply_buffer(k)%N_det + enddo + do i=1,H_apply_buffer(j)%N_det + do k=1,N_int + psi_det(k,1,i+N_det_old) = H_apply_buffer(j)%det(k,1,i) + psi_det(k,2,i+N_det_old) = H_apply_buffer(j)%det(k,2,i) + enddo + ASSERT (sum(popcnt(psi_det(:,1,i+N_det_old))) == elec_alpha_num) + ASSERT (sum(popcnt(psi_det(:,2,i+N_det_old))) == elec_beta_num ) + enddo + do k=1,N_states + do i=1,H_apply_buffer(j)%N_det + psi_coef(i+N_det_old,k) = H_apply_buffer(j)%coef(i,k) + enddo + enddo + !$OMP BARRIER + H_apply_buffer(j)%N_det = 0 + !$OMP END PARALLEL + call normalize(psi_coef,N_det) + SOFT_TOUCH N_det psi_det psi_coef + +end + + +subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc) + use bitmasks + implicit none + BEGIN_DOC + ! Fill the H_apply buffer with determiants for CISD + END_DOC + + integer, intent(in) :: n_selected, Nint, iproc + integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected) + integer :: i,j,k + integer :: new_size + PROVIDE H_apply_buffer_allocated + new_size = H_apply_buffer(iproc)%N_det + n_selected + if (new_size > H_apply_buffer(iproc)%sze) then + call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc) + endif + call omp_set_lock(H_apply_buffer_lock(1,iproc)) + do i=1,H_apply_buffer(iproc)%N_det + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) + enddo + do i=1,n_selected + do j=1,N_int + H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i) + H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i) + enddo + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num) + enddo + do j=1,N_states + do i=1,N_selected + H_apply_buffer(iproc)%coef(i,j) = 0.d0 + enddo + enddo + H_apply_buffer(iproc)%N_det = new_size + do i=1,H_apply_buffer(iproc)%N_det + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num) + ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num) + enddo + call omp_unset_lock(H_apply_buffer_lock(1,iproc)) +end + + diff --git a/src/Determinants/H_apply_template.f b/src/Determinants/H_apply_template.f new file mode 100644 index 00000000..a9a282ae --- /dev/null +++ b/src/Determinants/H_apply_template.f @@ -0,0 +1,542 @@ +subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_generator, iproc_in $parameters ) + use omp_lib + use bitmasks + implicit none + BEGIN_DOC + ! Generate all double excitations of key_in using the bit masks of holes and + ! particles. + ! Assume N_int is already provided. + END_DOC + 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),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) + integer, intent(in) :: iproc_in + integer(bit_kind), allocatable :: hole_save(:,:) + integer(bit_kind), allocatable :: key(:,:),hole(:,:), particle(:,:) + integer(bit_kind), allocatable :: hole_tmp(:,:), particle_tmp(:,:) + integer :: ii,i,jj,j,k,ispin,l + integer, allocatable :: occ_particle(:,:), occ_hole(:,:) + integer, allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:) + integer :: kk,pp,other_spin,key_idx + integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) + integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) + + double precision :: mo_bielec_integral + logical :: is_a_two_holes_two_particles + integer, allocatable :: ia_ja_pairs(:,:,:) + integer, allocatable :: ib_jb_pairs(:,:) + double precision :: diag_H_mat_elem + integer :: iproc + integer(omp_lock_kind), save :: lck, ifirst=0 + if (ifirst == 0) then +!$ call omp_init_lock(lck) + ifirst=1 + endif + + logical :: check_double_excitation + check_double_excitation = .True. + iproc = iproc_in + + + $initialization + + $omp_parallel +!$ iproc = omp_get_thread_num() + allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & + key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& + particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & + occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),& + occ_hole_tmp(N_int*bit_kind_size,2)) + $init_thread + + + + !!!! First couple hole particle + do j = 1, N_int + hole(j,1) = iand(hole_1(j,1),key_in(j,1)) + hole(j,2) = iand(hole_1(j,2),key_in(j,2)) + particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1)) + particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2)) + enddo + call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int) + call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int) + call bitstring_to_list(hole(1,1),occ_hole(1,1),N_elec_in_key_hole_1(1),N_int) + call bitstring_to_list(hole(1,2),occ_hole(1,2),N_elec_in_key_hole_1(2),N_int) + allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2), & + ib_jb_pairs(2,0:(elec_alpha_num)*mo_tot_num)) + + do ispin=1,2 + i=0 + do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole + i_a = occ_hole(ii,ispin) + ASSERT (i_a > 0) + ASSERT (i_a <= mo_tot_num) + + do jj=1,N_elec_in_key_part_1(ispin) !particle + j_a = occ_particle(jj,ispin) + ASSERT (j_a > 0) + ASSERT (j_a <= mo_tot_num) + i += 1 + ia_ja_pairs(1,i,ispin) = i_a + ia_ja_pairs(2,i,ispin) = j_a + enddo + enddo + ia_ja_pairs(1,0,ispin) = i + enddo + + key_idx = 0 + + integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b + integer(bit_kind) :: test(N_int,2) + double precision :: accu + logical, allocatable :: array_pairs(:,:) + allocate(array_pairs(mo_tot_num,mo_tot_num)) + accu = 0.d0 + do ispin=1,2 + other_spin = iand(ispin,1)+1 + if (abort_here) then + exit + endif + $omp_do + do ii=1,ia_ja_pairs(1,0,ispin) + if (abort_here) then + cycle + endif + i_a = ia_ja_pairs(1,ii,ispin) + ASSERT (i_a > 0) + ASSERT (i_a <= mo_tot_num) + j_a = ia_ja_pairs(2,ii,ispin) + ASSERT (j_a > 0) + ASSERT (j_a <= mo_tot_num) + hole = key_in + k = ishft(i_a-1,-bit_kind_shift)+1 + j = i_a-ishft(k-1,bit_kind_shift)-1 + hole(k,ispin) = ibclr(hole(k,ispin),j) + k_a = ishft(j_a-1,-bit_kind_shift)+1 + l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 + hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) + + !!!! Second couple hole particle + do j = 1, N_int + hole_tmp(j,1) = iand(hole_2(j,1),hole(j,1)) + hole_tmp(j,2) = iand(hole_2(j,2),hole(j,2)) + particle_tmp(j,1) = iand(xor(particl_2(j,1),hole(j,1)),particl_2(j,1)) + particle_tmp(j,2) = iand(xor(particl_2(j,2),hole(j,2)),particl_2(j,2)) + enddo + + call bitstring_to_list(particle_tmp(1,1),occ_particle_tmp(1,1),N_elec_in_key_part_2(1),N_int) + call bitstring_to_list(particle_tmp(1,2),occ_particle_tmp(1,2),N_elec_in_key_part_2(2),N_int) + call bitstring_to_list(hole_tmp (1,1),occ_hole_tmp (1,1),N_elec_in_key_hole_2(1),N_int) + call bitstring_to_list(hole_tmp (1,2),occ_hole_tmp (1,2),N_elec_in_key_hole_2(2),N_int) + + ! hole = a^(+)_j_a(ispin) a_i_a(ispin)|key_in> : mono exc :: orb(i_a,ispin) --> orb(j_a,ispin) + hole_save = hole + + ! Build array of the non-zero integrals of second excitation + $filter_integrals + if (ispin == 1) then + integer :: jjj + + i=0 + do kk = 1,N_elec_in_key_hole_2(other_spin) + i_b = occ_hole_tmp(kk,other_spin) + ASSERT (i_b > 0) + ASSERT (i_b <= mo_tot_num) + do jjj=1,N_elec_in_key_part_2(other_spin) ! particule + j_b = occ_particle_tmp(jjj,other_spin) + ASSERT (j_b > 0) + ASSERT (j_b <= mo_tot_num) + if (array_pairs(i_b,j_b)) then + i+= 1 + ib_jb_pairs(1,i) = i_b + ib_jb_pairs(2,i) = j_b + endif + enddo + enddo + ib_jb_pairs(1,0) = i + + do kk = 1,ib_jb_pairs(1,0) + hole = hole_save + i_b = ib_jb_pairs(1,kk) + j_b = ib_jb_pairs(2,kk) + k = ishft(i_b-1,-bit_kind_shift)+1 + j = i_b-ishft(k-1,bit_kind_shift)-1 + hole(k,other_spin) = ibclr(hole(k,other_spin),j) + key = hole + k = ishft(j_b-1,-bit_kind_shift)+1 + l = j_b-ishft(k-1,bit_kind_shift)-1 + key(k,other_spin) = ibset(key(k,other_spin),l) + $filter2h2p + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = key(k,1) + keys_out(k,2,key_idx) = key(k,2) + enddo + ASSERT (key_idx <= size_max) + if (key_idx == size_max) then + $keys_work + key_idx = 0 + endif + if (abort_here) then + exit + endif + enddo + endif + + ! does all the mono excitations of the same spin + i=0 + do kk = 1,N_elec_in_key_hole_2(ispin) + i_b = occ_hole_tmp(kk,ispin) + if (i_b <= i_a.or.i_b == j_a) cycle + ASSERT (i_b > 0) + ASSERT (i_b <= mo_tot_num) + do jjj=1,N_elec_in_key_part_2(ispin) ! particule + j_b = occ_particle_tmp(jjj,ispin) + ASSERT (j_b > 0) + ASSERT (j_b <= mo_tot_num) + if (j_b <= j_a) cycle + if (array_pairs(i_b,j_b)) then + i+= 1 + ib_jb_pairs(1,i) = i_b + ib_jb_pairs(2,i) = j_b + endif + enddo + enddo + ib_jb_pairs(1,0) = i + + do kk = 1,ib_jb_pairs(1,0) + hole = hole_save + i_b = ib_jb_pairs(1,kk) + j_b = ib_jb_pairs(2,kk) + k = ishft(i_b-1,-bit_kind_shift)+1 + j = i_b-ishft(k-1,bit_kind_shift)-1 + hole(k,ispin) = ibclr(hole(k,ispin),j) + key = hole + k = ishft(j_b-1,-bit_kind_shift)+1 + l = j_b-ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibset(key(k,ispin),l) + $filter2h2p + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = key(k,1) + keys_out(k,2,key_idx) = key(k,2) + enddo + ASSERT (key_idx <= size_max) + if (key_idx == size_max) then + $keys_work + key_idx = 0 + endif + if (abort_here) then + exit + endif + enddo ! kk + + enddo ! ii + $omp_enddo + enddo ! ispin + $keys_work + $deinit_thread + deallocate (ia_ja_pairs, ib_jb_pairs, & + keys_out, hole_save, & + key,hole, particle, hole_tmp,& + particle_tmp, occ_particle, & + occ_hole, occ_particle_tmp,& + occ_hole_tmp,array_pairs) + $omp_end_parallel + $finalization +end + +subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $parameters ) + use omp_lib + use bitmasks + implicit none + BEGIN_DOC + ! Generate all single excitations of key_in using the bit masks of holes and + ! particles. + ! Assume N_int is already provided. + END_DOC + 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) :: hole_1(N_int,2), particl_1(N_int,2) + integer, intent(in) :: iproc_in + integer(bit_kind),allocatable :: keys_out(:,:,:) + integer(bit_kind),allocatable :: hole_save(:,:) + integer(bit_kind),allocatable :: key(:,:),hole(:,:), particle(:,:) + integer(bit_kind),allocatable :: hole_tmp(:,:), particle_tmp(:,:) + integer(bit_kind),allocatable :: hole_2(:,:), particl_2(:,:) + integer :: ii,i,jj,j,k,ispin,l + integer,allocatable :: occ_particle(:,:), occ_hole(:,:) + integer,allocatable :: occ_particle_tmp(:,:), occ_hole_tmp(:,:) + integer,allocatable :: ib_jb_pairs(:,:) + integer :: kk,pp,other_spin,key_idx + integer :: N_elec_in_key_hole_1(2),N_elec_in_key_part_1(2) + integer :: N_elec_in_key_hole_2(2),N_elec_in_key_part_2(2) + logical :: is_a_two_holes_two_particles + + integer, allocatable :: ia_ja_pairs(:,:,:) + logical, allocatable :: array_pairs(:,:) + double precision :: diag_H_mat_elem + integer(omp_lock_kind), save :: lck, ifirst=0 + integer :: iproc + + logical :: check_double_excitation + iproc = iproc_in + + check_double_excitation = .True. + $check_double_excitation + + + if (ifirst == 0) then + ifirst=1 +!$ call omp_init_lock(lck) + endif + + $initialization + + $omp_parallel +!$ iproc = omp_get_thread_num() + allocate (keys_out(N_int,2,size_max), hole_save(N_int,2), & + key(N_int,2),hole(N_int,2), particle(N_int,2), hole_tmp(N_int,2),& + particle_tmp(N_int,2), occ_particle(N_int*bit_kind_size,2), & + occ_hole(N_int*bit_kind_size,2), occ_particle_tmp(N_int*bit_kind_size,2),& + occ_hole_tmp(N_int*bit_kind_size,2)) + $init_thread + !!!! First couple hole particle + do j = 1, N_int + hole(j,1) = iand(hole_1(j,1),key_in(j,1)) + hole(j,2) = iand(hole_1(j,2),key_in(j,2)) + particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1)) + particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2)) + enddo + + call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int) + call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int) + call bitstring_to_list(hole (1,1),occ_hole (1,1),N_elec_in_key_hole_1(1),N_int) + call bitstring_to_list(hole (1,2),occ_hole (1,2),N_elec_in_key_hole_1(2),N_int) + allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2)) + + do ispin=1,2 + i=0 + do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole + i_a = occ_hole(ii,ispin) + do jj=1,N_elec_in_key_part_1(ispin) !particule + j_a = occ_particle(jj,ispin) + i += 1 + ia_ja_pairs(1,i,ispin) = i_a + ia_ja_pairs(2,i,ispin) = j_a + enddo + enddo + ia_ja_pairs(1,0,ispin) = i + enddo + + key_idx = 0 + + integer :: i_a,j_a,i_b,j_b,k_a,l_a,k_b,l_b + integer(bit_kind) :: test(N_int,2) + double precision :: accu + accu = 0.d0 + do ispin=1,2 + other_spin = iand(ispin,1)+1 + $omp_do + do ii=1,ia_ja_pairs(1,0,ispin) + i_a = ia_ja_pairs(1,ii,ispin) + j_a = ia_ja_pairs(2,ii,ispin) + hole = key_in + k = ishft(i_a-1,-bit_kind_shift)+1 + j = i_a-ishft(k-1,bit_kind_shift)-1 + $filterhole + hole(k,ispin) = ibclr(hole(k,ispin),j) + k_a = ishft(j_a-1,-bit_kind_shift)+1 + l_a = j_a-ishft(k_a-1,bit_kind_shift)-1 + $filterparticle + hole(k_a,ispin) = ibset(hole(k_a,ispin),l_a) + $filter2h2p + key_idx += 1 + do k=1,N_int + keys_out(k,1,key_idx) = hole(k,1) + keys_out(k,2,key_idx) = hole(k,2) + enddo + if (key_idx == size_max) then + $keys_work + key_idx = 0 + endif + enddo ! ii + $omp_enddo + enddo ! ispin + $keys_work + $deinit_thread + deallocate (ia_ja_pairs, & + keys_out, hole_save, & + key,hole, particle, hole_tmp,& + particle_tmp, occ_particle, & + occ_hole, occ_particle_tmp,& + occ_hole_tmp) + $omp_end_parallel + $finalization + +end + + +subroutine $subroutine($params_main) + implicit none + use omp_lib + use bitmasks + BEGIN_DOC + ! Calls H_apply on the HF determinant and selects all connected single and double + ! excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script. + END_DOC + + $decls_main + + integer :: i_generator, nmax + double precision :: wall_0, wall_1 + integer(omp_lock_kind) :: lck + integer(bit_kind), allocatable :: mask(:,:,:) + integer :: ispin, k + integer :: iproc + + $initialization + PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators + + + nmax = mod( N_det_generators,nproc ) + + + !$ call omp_init_lock(lck) + call start_progress(N_det_generators,'Selection (norm)',0.d0) + + call wall_time(wall_0) + + iproc = 0 + allocate( mask(N_int,2,6) ) + do i_generator=1,nmax + + progress_bar(1) = i_generator + + if (abort_here) then + exit + endif + $skip + + ! Create bit masks for holes and particles + do ispin=1,2 + do k=1,N_int + mask(k,ispin,s_hole) = & + iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & + psi_det_generators(k,ispin,i_generator) ) + mask(k,ispin,s_part) = & + iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & + not(psi_det_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole1) = & + iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & + psi_det_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part1) = & + iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & + not(psi_det_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole2) = & + iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & + psi_det_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part2) = & + iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & + not(psi_det_generators(k,ispin,i_generator)) ) + enddo + enddo + if($do_double_excitations)then + 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_det_generators(1,1,i_generator), & + mask(1,1,s_hole ), mask(1,1,s_part ), & + i_generator, iproc $params_post) + endif + call wall_time(wall_1) + $printout_always + if (wall_1 - wall_0 > 2.d0) then + $printout_now + wall_0 = wall_1 + endif + enddo + + deallocate( mask ) + + !$OMP PARALLEL DEFAULT(SHARED) & + !$OMP PRIVATE(i_generator,wall_1,wall_0,ispin,k,mask,iproc) + call wall_time(wall_0) + !$ iproc = omp_get_thread_num() + allocate( mask(N_int,2,6) ) + !$OMP DO SCHEDULE(dynamic,1) + do i_generator=nmax+1,N_det_generators + if (iproc == 0) then + progress_bar(1) = i_generator + endif + if (abort_here) then + cycle + endif + $skip + + ! Create bit masks for holes and particles + do ispin=1,2 + do k=1,N_int + mask(k,ispin,s_hole) = & + iand(generators_bitmask(k,ispin,s_hole,i_bitmask_gen), & + psi_det_generators(k,ispin,i_generator) ) + mask(k,ispin,s_part) = & + iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), & + not(psi_det_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole1) = & + iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), & + psi_det_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part1) = & + iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), & + not(psi_det_generators(k,ispin,i_generator)) ) + mask(k,ispin,d_hole2) = & + iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), & + psi_det_generators(k,ispin,i_generator) ) + mask(k,ispin,d_part2) = & + iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), & + not (psi_det_generators(k,ispin,i_generator)) ) + enddo + enddo + + if($do_double_excitations)then + 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_det_generators(1,1,i_generator), & + mask(1,1,s_hole ), mask(1,1,s_part ), & + i_generator, iproc $params_post) + endif + !$ call omp_set_lock(lck) + call wall_time(wall_1) + $printout_always + if (wall_1 - wall_0 > 2.d0) then + $printout_now + wall_0 = wall_1 + endif + !$ call omp_unset_lock(lck) + enddo + !$OMP END DO + deallocate( mask ) + !$OMP END PARALLEL + !$ call omp_destroy_lock(lck) + + abort_here = abort_all + call stop_progress + + $copy_buffer + $generate_psi_guess + +end + diff --git a/src/Determinants/Makefile b/src/Determinants/Makefile new file mode 100644 index 00000000..092d879d --- /dev/null +++ b/src/Determinants/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=H_apply_template.f +OBJ= + +include $(QPACKAGE_ROOT)/src/Makefile.common diff --git a/src/Determinants/NEEDED_MODULES b/src/Determinants/NEEDED_MODULES new file mode 100644 index 00000000..824c75ed --- /dev/null +++ b/src/Determinants/NEEDED_MODULES @@ -0,0 +1 @@ +AOs Bielec_integrals Bitmask Electrons Ezfio_files MonoInts MOs Nuclei Output Utils diff --git a/src/Determinants/README.rst b/src/Determinants/README.rst new file mode 100644 index 00000000..445c8b5e --- /dev/null +++ b/src/Determinants/README.rst @@ -0,0 +1,696 @@ +=========== +Dets Module +=========== + +This module contains the determinants of the CI wave function. + +H is applied on the list of generator determinants. Selected determinants +are added into the *H_apply buffer*. Then the new wave function is +constructred as the concatenation of the odl wave function and +some determinants of the H_apply buffer. Generator determinants are built +as a subset of the determinants of the wave function. + + +Assumptions +=========== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* The MOs are orthonormal +* All the determinants have the same number of electrons +* The determinants are orthonormal +* The number of generator determinants <= the number of determinants +* All the determinants in the H_apply buffer are supposed to be different from the + wave function determinants +* All the determinants in the H_apply buffer are supposed to be unique + + +Needed Modules +============== + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +* `AOs `_ +* `Bielec_integrals `_ +* `Bitmask `_ +* `Electrons `_ +* `Ezfio_files `_ +* `MonoInts `_ +* `MOs `_ +* `Nuclei `_ +* `Output `_ +* `Utils `_ + +Documentation +============= + +.. Do not edit this section. It was auto-generated from the +.. NEEDED_MODULES file. + +`copy_h_apply_buffer_to_wf `_ + Copies the H_apply buffer to psi_coef. You need to touch psi_det, psi_coef and N_det + after calling this function. + After calling this subroutine, N_det, psi_det and psi_coef need to be touched + +`fill_h_apply_buffer_no_selection `_ + Fill the H_apply buffer with determiants for CISD + +`h_apply_buffer_allocated `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. + +`h_apply_buffer_lock `_ + Buffer of determinants/coefficients/perturbative energy for H_apply. + Uninitialized. Filled by H_apply subroutines. + +`resize_h_apply_buffer `_ + Undocumented + +`cisd_sc2 `_ + CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + Initial guess vectors are not necessarily orthonormal + +`connected_to_ref `_ + Undocumented + +`connected_to_ref_by_mono `_ + Undocumented + +`det_search_key `_ + Return an integer*8 corresponding to a determinant index for searching + +`get_index_in_psi_det_sorted_bit `_ + Returns the index of the determinant in the ``psi_det_sorted_bit`` array + +`is_in_wavefunction `_ + True if the determinant ``det`` is in the wave function + +`occ_pattern_search_key `_ + Return an integer*8 corresponding to a determinant index for searching + +`do_mono_excitation `_ + Apply the mono excitation operator : a^{dager}_(i_particle) a_(i_hole) of spin = ispin + on key_in + ispin = 1 == alpha + ispin = 2 == beta + i_ok = 1 == the excitation is possible + i_ok = -1 == the excitation is not possible + +`davidson_converged `_ + True if the Davidson algorithm is converged + +`davidson_criterion `_ + Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] + +`davidson_diag `_ + Davidson diagonalization. + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + iunit : Unit number for the I/O + .br + Initial guess vectors are not necessarily orthonormal + +`davidson_diag_hjj `_ + Davidson diagonalization with specific diagonal elements of the H matrix + .br + H_jj : specific diagonal H matrix elements to diagonalize de Davidson + .br + dets_in : bitmasks corresponding to determinants + .br + u_in : guess coefficients on the various states. Overwritten + on exit + .br + dim_in : leftmost dimension of u_in + .br + sze : Number of determinants + .br + N_st : Number of eigenstates + .br + iunit : Unit for the I/O + .br + Initial guess vectors are not necessarily orthonormal + +`davidson_iter_max `_ + Max number of Davidson iterations + +`davidson_sze_max `_ + Max number of Davidson sizes + +`davidson_threshold `_ + Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] + +`one_body_dm_mo `_ + One-body density matrix + +`one_body_dm_mo_alpha `_ + Alpha and beta one-body density matrix for each state + +`one_body_dm_mo_beta `_ + Alpha and beta one-body density matrix for each state + +`one_body_single_double_dm_mo_alpha `_ + Alpha and beta one-body density matrix for each state + +`one_body_single_double_dm_mo_beta `_ + Alpha and beta one-body density matrix for each state + +`one_body_spin_density_mo `_ + rho(alpha) - rho(beta) + +`save_natural_mos `_ + Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + +`set_natural_mos `_ + Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + +`state_average_weight `_ + Weights in the state-average calculation of the density matrix + +`det_svd `_ + Computes the SVD of the Alpha x Beta determinant coefficient matrix + +`filter_3_highest_electrons `_ + Returns a determinant with only the 3 highest electrons + +`int_of_3_highest_electrons `_ + Returns an integer*8 as : + .br + |_<--- 21 bits ---><--- 21 bits ---><--- 21 bits --->| + .br + |0<--- i1 ---><--- i2 ---><--- i3 --->| + .br + It encodes the value of the indices of the 3 highest MOs + in descending order + .br + +`max_degree_exc `_ + Maximum degree of excitation in the wf + +`n_det `_ + Number of determinants in the wave function + +`psi_average_norm_contrib `_ + Contribution of determinants to the state-averaged density + +`psi_average_norm_contrib_sorted `_ + Wave function sorted by determinants contribution to the norm (state-averaged) + +`psi_coef `_ + The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file + is empty + +`psi_coef_sorted `_ + Wave function sorted by determinants contribution to the norm (state-averaged) + +`psi_coef_sorted_ab `_ + Determinants on which we apply . + They are sorted by the 3 highest electrons in the alpha part, + then by the 3 highest electrons in the beta part to accelerate + the research of connected determinants. + +`psi_coef_sorted_bit `_ + Determinants on which we apply for perturbation. + They are sorted by determinants interpreted as integers. Useful + to accelerate the search of a random determinant in the wave + function. + +`psi_det `_ + The wave function determinants. Initialized with Hartree-Fock if the EZFIO file + is empty + +`psi_det_size `_ + Size of the psi_det/psi_coef arrays + +`psi_det_sorted `_ + Wave function sorted by determinants contribution to the norm (state-averaged) + +`psi_det_sorted_ab `_ + Determinants on which we apply . + They are sorted by the 3 highest electrons in the alpha part, + then by the 3 highest electrons in the beta part to accelerate + the research of connected determinants. + +`psi_det_sorted_bit `_ + Determinants on which we apply for perturbation. + They are sorted by determinants interpreted as integers. Useful + to accelerate the search of a random determinant in the wave + function. + +`psi_det_sorted_next_ab `_ + Determinants on which we apply . + They are sorted by the 3 highest electrons in the alpha part, + then by the 3 highest electrons in the beta part to accelerate + the research of connected determinants. + +`read_dets `_ + Reads the determinants from the EZFIO file + +`save_wavefunction `_ + Save the wave function into the EZFIO file + +`save_wavefunction_general `_ + Save the wave function into the EZFIO file + +`save_wavefunction_unsorted `_ + Save the wave function into the EZFIO file + +`sort_dets_by_3_highest_electrons `_ + Determinants on which we apply . + They are sorted by the 3 highest electrons in the alpha part, + then by the 3 highest electrons in the beta part to accelerate + the research of connected determinants. + +`sort_dets_by_det_search_key `_ + Determinants are sorted are sorted according to their det_search_key. + Useful to accelerate the search of a random determinant in the wave + function. + +`double_exc_bitmask `_ + double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1 + double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1 + double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2 + double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2 + for a given couple of hole/particle excitations i. + +`n_double_exc_bitmasks `_ + Number of double excitation bitmasks + +`n_single_exc_bitmasks `_ + Number of single excitation bitmasks + +`single_exc_bitmask `_ + single_exc_bitmask(:,1,i) is the bitmask for holes + single_exc_bitmask(:,2,i) is the bitmask for particles + for a given couple of hole/particle excitations i. + +`ci_eigenvectors `_ + Eigenvectors/values of the CI matrix + +`ci_eigenvectors_s2 `_ + Eigenvectors/values of the CI matrix + +`ci_electronic_energy `_ + Eigenvectors/values of the CI matrix + +`ci_energy `_ + N_states lowest eigenvalues of the CI matrix + +`diag_algorithm `_ + Diagonalization algorithm (Davidson or Lapack) + +`diagonalize_ci `_ + Replace the coefficients of the CI states by the coefficients of the + eigenstates of the CI matrix + +`ci_sc2_eigenvectors `_ + Eigenvectors/values of the CI matrix + +`ci_sc2_electronic_energy `_ + Eigenvectors/values of the CI matrix + +`ci_sc2_energy `_ + N_states_diag lowest eigenvalues of the CI matrix + +`diagonalize_ci_sc2 `_ + Replace the coefficients of the CI states_diag by the coefficients of the + eigenstates of the CI matrix + +`threshold_convergence_sc2 `_ + convergence of the correlation energy of SC2 iterations + +`ci_eigenvectors_mono `_ + Eigenvectors/values of the CI matrix + +`ci_eigenvectors_s2_mono `_ + Eigenvectors/values of the CI matrix + +`ci_electronic_energy_mono `_ + Eigenvectors/values of the CI matrix + +`diagonalize_ci_mono `_ + Replace the coefficients of the CI states by the coefficients of the + eigenstates of the CI matrix + +`apply_mono `_ + Undocumented + +`filter_connected `_ + Filters out the determinants that are not connected by H + .br + returns the array idx which contains the index of the + .br + determinants in the array key1 that interact + .br + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + +`filter_connected_davidson `_ + Filters out the determinants that are not connected by H + returns the array idx which contains the index of the + determinants in the array key1 that interact + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + key1 should come from psi_det_sorted_ab. + +`filter_connected_i_h_psi0 `_ + returns the array idx which contains the index of the + .br + determinants in the array key1 that interact + .br + via the H operator with key2. + .br + idx(0) is the number of determinants that interact with key1 + +`filter_connected_i_h_psi0_sc2 `_ + standard filter_connected_i_H_psi but returns in addition + .br + the array of the index of the non connected determinants to key1 + .br + in order to know what double excitation can be repeated on key1 + .br + idx_repeat(0) is the number of determinants that can be used + .br + to repeat the excitations + +`filter_connected_sorted_ab `_ + Filters out the determinants that are not connected by H + returns the array idx which contains the index of the + determinants in the array key1 that interact + via the H operator with key2. + idx(0) is the number of determinants that interact with key1 + .br + Determinants are taken from the psi_det_sorted_ab array + +`put_gess `_ + Undocumented + +`det_to_occ_pattern `_ + Transform a determinant to an occupation pattern + +`make_s2_eigenfunction `_ + Undocumented + +`n_occ_pattern `_ + array of the occ_pattern present in the wf + psi_occ_pattern(:,1,j) = jth occ_pattern of the wave function : represent all the single occupation + psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation + +`occ_pattern_to_dets `_ + Generate all possible determinants for a give occ_pattern + +`occ_pattern_to_dets_size `_ + Number of possible determinants for a given occ_pattern + +`psi_occ_pattern `_ + array of the occ_pattern present in the wf + psi_occ_pattern(:,1,j) = jth occ_pattern of the wave function : represent all the single occupation + psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation + +`rec_occ_pattern_to_dets `_ + Undocumented + +`n_states_diag `_ + Number of states to consider for the diagonalization + +`pouet `_ + Undocumented + +`routine `_ + Undocumented + +`idx_cas `_ + CAS wave function, defined from the application of the CAS bitmask on the + determinants. idx_cas gives the indice of the CAS determinant in psi_det. + +`idx_non_cas `_ + Set of determinants which are not part of the CAS, defined from the application + of the CAS bitmask on the determinants. + idx_non_cas gives the indice of the determinant in psi_det. + +`n_det_cas `_ + CAS wave function, defined from the application of the CAS bitmask on the + determinants. idx_cas gives the indice of the CAS determinant in psi_det. + +`n_det_non_cas `_ + Set of determinants which are not part of the CAS, defined from the application + of the CAS bitmask on the determinants. + idx_non_cas gives the indice of the determinant in psi_det. + +`psi_cas `_ + CAS wave function, defined from the application of the CAS bitmask on the + determinants. idx_cas gives the indice of the CAS determinant in psi_det. + +`psi_cas_coef `_ + CAS wave function, defined from the application of the CAS bitmask on the + determinants. idx_cas gives the indice of the CAS determinant in psi_det. + +`psi_cas_coef_sorted_bit `_ + CAS determinants sorted to accelerate the search of a random determinant in the wave + function. + +`psi_cas_sorted_bit `_ + CAS determinants sorted to accelerate the search of a random determinant in the wave + function. + +`psi_non_cas `_ + Set of determinants which are not part of the CAS, defined from the application + of the CAS bitmask on the determinants. + idx_non_cas gives the indice of the determinant in psi_det. + +`psi_non_cas_coef `_ + Set of determinants which are not part of the CAS, defined from the application + of the CAS bitmask on the determinants. + idx_non_cas gives the indice of the determinant in psi_det. + +`psi_non_cas_coef_sorted_bit `_ + CAS determinants sorted to accelerate the search of a random determinant in the wave + function. + +`psi_non_cas_sorted_bit `_ + CAS determinants sorted to accelerate the search of a random determinant in the wave + function. + +`bi_elec_ref_bitmask_energy `_ + Energy of the reference bitmask used in Slater rules + +`kinetic_ref_bitmask_energy `_ + Energy of the reference bitmask used in Slater rules + +`mono_elec_ref_bitmask_energy `_ + Energy of the reference bitmask used in Slater rules + +`nucl_elec_ref_bitmask_energy `_ + Energy of the reference bitmask used in Slater rules + +`ref_bitmask_energy `_ + Energy of the reference bitmask used in Slater rules + +`expected_s2 `_ + Expected value of S2 : S*(S+1) + +`get_s2 `_ + Returns + +`get_s2_u0 `_ + Undocumented + +`s2_values `_ + array of the averaged values of the S^2 operator on the various states + +`s_z `_ + z component of the Spin + +`s_z2_sz `_ + z component of the Spin + +`prog_save_casino `_ + Undocumented + +`save_casino `_ + Undocumented + +`save_dets_qmcchem `_ + Undocumented + +`save_for_qmc `_ + Undocumented + +`save_natorb `_ + Undocumented + +`a_operator `_ + Needed for diag_H_mat_elem + +`ac_operator `_ + Needed for diag_H_mat_elem + +`decode_exc `_ + Decodes the exc arrays returned by get_excitation. + h1,h2 : Holes + p1,p2 : Particles + s1,s2 : Spins (1:alpha, 2:beta) + degree : Degree of excitation + +`det_connections `_ + Build connection proxy between determinants + +`diag_h_mat_elem `_ + Computes + +`get_double_excitation `_ + Returns the two excitation operators between two doubly excited determinants and the phase + +`get_excitation `_ + Returns the excitation operators between two determinants and the phase + +`get_excitation_degree `_ + Returns the excitation degree between two determinants + +`get_excitation_degree_vector `_ + Applies get_excitation_degree to an array of determinants + +`get_mono_excitation `_ + Returns the excitation operator between two singly excited determinants and the phase + +`get_occ_from_key `_ + Returns a list of occupation numbers from a bitstring + +`h_u_0 `_ + Computes v_0 = H|u_0> + .br + n : number of determinants + .br + H_jj : array of + +`i_h_j `_ + Returns where i and j are determinants + +`i_h_j_verbose `_ + Returns where i and j are determinants + +`i_h_psi `_ + for the various Nstates + +`i_h_psi_sc2 `_ + for the various Nstate + .br + returns in addition + .br + the array of the index of the non connected determinants to key1 + .br + in order to know what double excitation can be repeated on key1 + .br + idx_repeat(0) is the number of determinants that can be used + .br + to repeat the excitations + +`i_h_psi_sc2_verbose `_ + for the various Nstate + .br + returns in addition + .br + the array of the index of the non connected determinants to key1 + .br + in order to know what double excitation can be repeated on key1 + .br + idx_repeat(0) is the number of determinants that can be used + .br + to repeat the excitations + +`i_h_psi_sec_ord `_ + for the various Nstates + +`n_con_int `_ + Number of integers to represent the connections between determinants + +`create_wf_of_psi_svd_matrix `_ + Matrix of wf coefficients. Outer product of alpha and beta determinants + +`generate_all_alpha_beta_det_products `_ + Create a wave function from all possible alpha x beta determinants + +`get_index_in_psi_det_alpha_unique `_ + Returns the index of the determinant in the ``psi_det_alpha_unique`` array + +`get_index_in_psi_det_beta_unique `_ + Returns the index of the determinant in the ``psi_det_beta_unique`` array + +`n_det_alpha_unique `_ + Unique alpha determinants + +`n_det_beta_unique `_ + Unique beta determinants + +`psi_det_alpha `_ + List of alpha determinants of psi_det + +`psi_det_alpha_unique `_ + Unique alpha determinants + +`psi_det_beta `_ + List of beta determinants of psi_det + +`psi_det_beta_unique `_ + Unique beta determinants + +`psi_svd_alpha `_ + SVD wave function + +`psi_svd_beta `_ + SVD wave function + +`psi_svd_coefs `_ + SVD wave function + +`psi_svd_matrix `_ + Matrix of wf coefficients. Outer product of alpha and beta determinants + +`psi_svd_matrix_columns `_ + Matrix of wf coefficients. Outer product of alpha and beta determinants + +`psi_svd_matrix_rows `_ + Matrix of wf coefficients. Outer product of alpha and beta determinants + +`psi_svd_matrix_values `_ + Matrix of wf coefficients. Outer product of alpha and beta determinants + +`spin_det_search_key `_ + Return an integer*8 corresponding to a determinant index for searching + +`write_spindeterminants `_ + Undocumented + +`cisd `_ + Undocumented + +`h_matrix_all_dets `_ + H matrix on the basis of the slater determinants defined by psi_det + + + diff --git a/src/Determinants/SC2.irp.f b/src/Determinants/SC2.irp.f new file mode 100644 index 00000000..440b2870 --- /dev/null +++ b/src/Determinants/SC2.irp.f @@ -0,0 +1,215 @@ +subroutine CISD_SC2(dets_in,u_in,energies,dim_in,sze,N_st,Nint,convergence) + use bitmasks + implicit none + BEGIN_DOC + ! CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not) + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: energies(N_st) + double precision, intent(in) :: convergence + ASSERT (N_st > 0) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + integer :: iter + integer :: i,j,k,l,m + logical :: converged + double precision :: overlap(N_st,N_st) + double precision :: u_dot_v, u_dot_u + + integer :: degree,N_double,index_hf + double precision :: hij_elec, e_corr_double,e_corr,diag_h_mat_elem,inv_c0 + double precision :: e_corr_double_before,accu,cpu_2,cpu_1 + integer,allocatable :: degree_exc(:), index_double(:) + integer :: i_ok + double precision,allocatable :: e_corr_array(:),H_jj_ref(:),H_jj_dressed(:),hij_double(:) + integer(bit_kind), allocatable :: doubles(:,:,:) + + + allocate (doubles(Nint,2,sze),e_corr_array(sze),H_jj_ref(sze),H_jj_dressed(sze),& + index_double(sze), degree_exc(sze), hij_double(sze)) + call write_time(output_determinants) + write(output_determinants,'(A)') '' + write(output_determinants,'(A)') 'CISD SC2' + write(output_determinants,'(A)') '========' + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(sze,N_st, & + !$OMP H_jj_ref,Nint,dets_in,u_in) & + !$OMP PRIVATE(i) + + !$OMP DO SCHEDULE(guided) + do i=1,sze + H_jj_ref(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) + enddo + !$OMP END DO NOWAIT + !$OMP END PARALLEL + + N_double = 0 + e_corr = 0.d0 + e_corr_double = 0.d0 + do i = 1, sze + call get_excitation_degree(ref_bitmask,dets_in(1,1,i),degree,Nint) + degree_exc(i) = degree+1 + if(degree==0)then + index_hf=i + else if (degree == 2)then + N_double += 1 + index_double(N_double) = i + doubles(:,:,N_double) = dets_in(:,:,i) + call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec) + hij_double(N_double) = hij_elec + e_corr_array(N_double) = u_in(i,1)* hij_elec + e_corr_double += e_corr_array(N_double) + e_corr += e_corr_array(N_double) + else if (degree == 1)then + call i_H_j(ref_bitmask,dets_in(1,1,i),Nint,hij_elec) + e_corr += u_in(i,1)* hij_elec + endif + enddo + inv_c0 = 1.d0/u_in(index_hf,1) + do i = 1, N_double + e_corr_array(i) = e_corr_array(i) * inv_c0 + enddo + e_corr = e_corr * inv_c0 + e_corr_double = e_corr_double * inv_c0 + converged = .False. + e_corr_double_before = e_corr_double + iter = 0 + do while (.not.converged) + if (abort_here) then + exit + endif + iter +=1 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,j,degree,accu) & + !$OMP SHARED(H_jj_dressed,sze,H_jj_ref,index_hf,N_int,N_double,& + !$OMP dets_in,doubles,degree_exc,e_corr_array,e_corr_double) + !$OMP DO SCHEDULE(STATIC) + do i=1,sze + H_jj_dressed(i) = H_jj_ref(i) + if (i==index_hf)cycle + accu = -e_corr_double + select case (N_int) + case (1) + do j=1,N_double + degree = & + popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + & + popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) + + if (degree<=ishft(degree_exc(i),1)) then + accu += e_corr_array(j) + endif + enddo + case (2) + do j=1,N_double + degree = & + popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + & + popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) + & + popcnt(xor( dets_in(2,1,i),doubles(2,1,j))) + & + popcnt(xor( dets_in(2,2,i),doubles(2,2,j))) + + if (degree<=ishft(degree_exc(i),1)) then + accu += e_corr_array(j) + endif + enddo + case (3) + do j=1,N_double + degree = & + popcnt(xor( dets_in(1,1,i),doubles(1,1,j))) + & + popcnt(xor( dets_in(1,2,i),doubles(1,2,j))) + & + popcnt(xor( dets_in(2,1,i),doubles(2,1,j))) + & + popcnt(xor( dets_in(2,2,i),doubles(2,2,j))) + & + popcnt(xor( dets_in(3,1,i),doubles(3,1,j))) + & + popcnt(xor( dets_in(3,2,i),doubles(3,2,j))) + + if (degree<=ishft(degree_exc(i),1)) then + accu += e_corr_array(j) + endif + enddo + case default + do j=1,N_double + call get_excitation_degree(dets_in(1,1,i),doubles(1,1,j),degree,N_int) + if (degree<=degree_exc(i)) then + accu += e_corr_array(j) + endif + enddo + end select + H_jj_dressed(i) -= accu + enddo + !$OMP END DO + !$OMP END PARALLEL + + if(sze<=N_det_max_jacobi)then + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:),H_matrix_tmp(:,:) + allocate (H_matrix_tmp(size(H_matrix_all_dets,1),sze),eigenvalues(sze),eigenvectors(size(H_matrix_all_dets,1),sze)) + do j=1,sze + do i=1,sze + H_matrix_tmp(i,j) = H_matrix_all_dets(i,j) + enddo + enddo + do i = 1,sze + H_matrix_tmp(i,i) = H_jj_dressed(i) + enddo + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_tmp,size(H_matrix_all_dets,1),sze) + do j=1,min(N_states_diag,sze) + do i=1,sze + u_in(i,j) = eigenvectors(i,j) + enddo + energies(j) = eigenvalues(j) + enddo + deallocate (H_matrix_tmp, eigenvalues, eigenvectors) + else + call davidson_diag_hjj(dets_in,u_in,H_jj_dressed,energies,dim_in,sze,N_st,Nint,output_determinants) + endif + + e_corr_double = 0.d0 + inv_c0 = 1.d0/u_in(index_hf,1) + do i = 1, N_double + e_corr_array(i) = u_in(index_double(i),1)*inv_c0 * hij_double(i) + e_corr_double += e_corr_array(i) + enddo + write(output_determinants,'(A,I3)') 'SC2 Iteration ', iter + write(output_determinants,'(A)') '------------------' + write(output_determinants,'(A)') '' + write(output_determinants,'(A)') '===== ================' + write(output_determinants,'(A)') 'State Energy ' + write(output_determinants,'(A)') '===== ================' + do i=1,N_st + write(output_determinants,'(I5,X,F16.10)') i, energies(i)+nuclear_repulsion + enddo + write(output_determinants,'(A)') '===== ================' + write(output_determinants,'(A)') '' + call write_double(output_determinants,(e_corr_double - e_corr_double_before),& + 'Delta(E_corr)') + converged = dabs(e_corr_double - e_corr_double_before) < convergence + converged = converged .or. abort_here + if (converged) then + exit + endif + e_corr_double_before = e_corr_double + + enddo + + call write_time(output_determinants) + deallocate (doubles,e_corr_array,H_jj_ref,H_jj_dressed, & + index_double, degree_exc, hij_double) + +end + + diff --git a/src/Determinants/connected_to_ref.irp.f b/src/Determinants/connected_to_ref.irp.f new file mode 100644 index 00000000..2d40b621 --- /dev/null +++ b/src/Determinants/connected_to_ref.irp.f @@ -0,0 +1,357 @@ +integer*8 function det_search_key(det,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Return an integer*8 corresponding to a determinant index for searching + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det(Nint,2) + integer :: i + det_search_key = iand(det(1,1),det(1,2)) + do i=2,Nint + det_search_key = ieor(det_search_key,iand(det(i,1),det(i,2))) + enddo +end + + +integer*8 function occ_pattern_search_key(det,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Return an integer*8 corresponding to a determinant index for searching + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det(Nint,2) + integer :: i + occ_pattern_search_key = ieor(det(1,1),det(1,2)) + do i=2,Nint + occ_pattern_search_key = ieor(occ_pattern_search_key,iand(det(i,1),det(i,2))) + enddo +end + + + +logical function is_in_wavefunction(key,Nint,Ndet) + use bitmasks + implicit none + BEGIN_DOC +! True if the determinant ``det`` is in the wave function + END_DOC + integer, intent(in) :: Nint, Ndet + integer(bit_kind), intent(in) :: key(Nint,2) + integer, external :: get_index_in_psi_det_sorted_bit + + !DIR$ FORCEINLINE + is_in_wavefunction = get_index_in_psi_det_sorted_bit(key,Nint) > 0 +end + +integer function get_index_in_psi_det_sorted_bit(key,Nint) + use bitmasks + BEGIN_DOC +! Returns the index of the determinant in the ``psi_det_sorted_bit`` array + END_DOC + implicit none + + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key(Nint,2) + + integer :: i, ibegin, iend, istep, l + integer*8 :: det_ref, det_search + 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 + + !DIR$ FORCEINLINE + det_ref = det_search_key(key,Nint) + !DIR$ FORCEINLINE + det_search = det_search_key(psi_det_sorted_bit(1,1,1),Nint) + + istep = ishft(iend-ibegin,-1) + i=ibegin+istep + do while (istep > 0) + !DIR$ FORCEINLINE + det_search = det_search_key(psi_det_sorted_bit(1,1,i),Nint) + if ( det_search > det_ref ) then + iend = i + else if ( det_search == det_ref ) then + exit + else + ibegin = i + endif + istep = ishft(iend-ibegin,-1) + i = ibegin + istep + end do + + !DIR$ FORCEINLINE + do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref) + i = i-1 + if (i == 0) then + exit + endif + enddo + i += 1 + + if (i > N_det) then + return + endif + + !DIR$ FORCEINLINE + do while (det_search_key(psi_det_sorted_bit(1,1,i),Nint) == det_ref) + if ( (key(1,1) /= psi_det_sorted_bit(1,1,i)).or. & + (key(1,2) /= psi_det_sorted_bit(1,2,i)) ) then + continue + else + is_in_wavefunction = .True. + !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. + 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 +! exit + return + endif + + enddo + +! 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 +! END DEBUG is_in_wf +end + +integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet) + use bitmasks + implicit none + integer, intent(in) :: Nint, N_past_in, Ndet + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + + integer :: N_past + integer :: i, l + integer :: degree_x2 + logical :: t + double precision :: hij_elec + + ! output : 0 : not connected + ! i : connected to determinant i of the past + ! -i : is the ith determinant of the refernce wf keys + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + connected_to_ref = 0 + N_past = max(1,N_past_in) + if (Nint == 1) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + if (degree_x2 > 4) then + cycle + else + connected_to_ref = i + return + endif + enddo + + return + + + else if (Nint==2) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + if (degree_x2 > 4) then + cycle + else + connected_to_ref = i + return + endif + enddo + + return + + else if (Nint==3) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + & + popcnt(xor( key(3,1), keys(3,1,i))) + & + popcnt(xor( key(3,2), keys(3,2,i))) + if (degree_x2 > 4) then + cycle + else + connected_to_ref = i + return + endif + enddo + + return + + else + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +& + popcnt(xor( key(l,2), keys(l,2,i))) + enddo + if (degree_x2 > 4) then + cycle + else + connected_to_ref = i + return + endif + enddo + + endif + +end + + + +integer function connected_to_ref_by_mono(key,keys,Nint,N_past_in,Ndet) + use bitmasks + implicit none + integer, intent(in) :: Nint, N_past_in, Ndet + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + + integer :: N_past + integer :: i, l + integer :: degree_x2 + logical :: t + double precision :: hij_elec + + ! output : 0 : not connected + ! i : connected to determinant i of the past + ! -i : is the ith determinant of the refernce wf keys + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + connected_to_ref_by_mono = 0 + N_past = max(1,N_past_in) + if (Nint == 1) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + return + + + else if (Nint==2) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + return + + else if (Nint==3) then + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + & + popcnt(xor( key(2,1), keys(2,1,i))) + & + popcnt(xor( key(2,2), keys(2,2,i))) + & + popcnt(xor( key(3,1), keys(3,1,i))) + & + popcnt(xor( key(3,2), keys(3,2,i))) + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + return + + else + + do i=N_past-1,1,-1 + degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + & + popcnt(xor( key(1,2), keys(1,2,i))) + !DEC$ LOOP COUNT MIN(3) + do l=2,Nint + degree_x2 = degree_x2 + popcnt(xor( key(l,1), keys(l,1,i))) +& + popcnt(xor( key(l,2), keys(l,2,i))) + enddo + if (degree_x2 > 3.and. degree_x2 <5) then + cycle + else if (degree_x2 == 4)then + cycle + else if(degree_x2 == 2)then + connected_to_ref_by_mono = i + return + endif + enddo + + endif + +end + + diff --git a/src/Determinants/create_excitations.irp.f b/src/Determinants/create_excitations.irp.f new file mode 100644 index 00000000..a33525c7 --- /dev/null +++ b/src/Determinants/create_excitations.irp.f @@ -0,0 +1,36 @@ +subroutine do_mono_excitation(key_in,i_hole,i_particle,ispin,i_ok) + implicit none + BEGIN_DOC + ! Apply the mono excitation operator : a^{dager}_(i_particle) a_(i_hole) of spin = ispin + ! on key_in + ! ispin = 1 == alpha + ! ispin = 2 == beta + ! i_ok = 1 == the excitation is possible + ! i_ok = -1 == the excitation is not possible + END_DOC + integer, intent(in) :: i_hole,i_particle,ispin + integer(bit_kind), intent(inout) :: key_in(N_int,2) + integer, intent(out) :: i_ok + integer :: k,j,i + use bitmasks + ASSERT (i_hole > 0 ) + ASSERT (i_particle <= mo_tot_num) + i_ok = 1 + ! hole + k = ishft(i_hole-1,-bit_kind_shift)+1 + j = i_hole-ishft(k-1,bit_kind_shift)-1 + key_in(k,ispin) = ibclr(key_in(k,ispin),j) + + ! particle + k = ishft(i_particle-1,-bit_kind_shift)+1 + j = i_particle-ishft(k-1,bit_kind_shift)-1 + key_in(k,ispin) = ibset(key_in(k,ispin),j) + integer :: n_elec_tmp + n_elec_tmp = 0 + do i = 1, N_int + n_elec_tmp += popcnt(key_in(i,1)) + popcnt(key_in(i,2)) + enddo + if(n_elec_tmp .ne. elec_num)then + i_ok = -1 + endif +end diff --git a/src/Determinants/davidson.irp.f b/src/Determinants/davidson.irp.f new file mode 100644 index 00000000..bdc979c4 --- /dev/null +++ b/src/Determinants/davidson.irp.f @@ -0,0 +1,418 @@ +BEGIN_PROVIDER [ integer, davidson_iter_max ] + implicit none + BEGIN_DOC + ! Max number of Davidson iterations + END_DOC + davidson_iter_max = 100 +END_PROVIDER + +BEGIN_PROVIDER [ integer, davidson_sze_max ] + implicit none + BEGIN_DOC + ! Max number of Davidson sizes + END_DOC + ASSERT (davidson_sze_max <= davidson_iter_max) + davidson_sze_max = 8*N_states_diag +END_PROVIDER + +subroutine davidson_diag(dets_in,u_in,energies,dim_in,sze,N_st,Nint,iunit) + use bitmasks + implicit none + BEGIN_DOC + ! Davidson diagonalization. + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! iunit : Unit number for the I/O + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint, iunit + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: energies(N_st) + double precision, allocatable :: H_jj(:) + + double precision :: diag_h_mat_elem + integer :: i + ASSERT (N_st > 0) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + PROVIDE mo_bielec_integrals_in_map + allocate(H_jj(sze)) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(sze,H_jj,dets_in,Nint) & + !$OMP PRIVATE(i) + !$OMP DO SCHEDULE(guided) + do i=1,sze + H_jj(i) = diag_h_mat_elem(dets_in(1,1,i),Nint) + enddo + !$OMP END DO + !$OMP END PARALLEL + + call davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit) + deallocate (H_jj) +end + +subroutine davidson_diag_hjj(dets_in,u_in,H_jj,energies,dim_in,sze,N_st,Nint,iunit) + use bitmasks + implicit none + BEGIN_DOC + ! Davidson diagonalization with specific diagonal elements of the H matrix + ! + ! H_jj : specific diagonal H matrix elements to diagonalize de Davidson + ! + ! dets_in : bitmasks corresponding to determinants + ! + ! u_in : guess coefficients on the various states. Overwritten + ! on exit + ! + ! dim_in : leftmost dimension of u_in + ! + ! sze : Number of determinants + ! + ! N_st : Number of eigenstates + ! + ! iunit : Unit for the I/O + ! + ! Initial guess vectors are not necessarily orthonormal + END_DOC + integer, intent(in) :: dim_in, sze, N_st, Nint + integer(bit_kind), intent(in) :: dets_in(Nint,2,sze) + double precision, intent(in) :: H_jj(sze) + integer, intent(in) :: iunit + double precision, intent(inout) :: u_in(dim_in,N_st) + double precision, intent(out) :: energies(N_st) + + integer :: iter + integer :: i,j,k,l,m + logical :: converged + + double precision :: overlap(N_st,N_st) + double precision :: u_dot_v, u_dot_u + + integer, allocatable :: kl_pairs(:,:) + integer :: k_pairs, kl + + integer :: iter2 + double precision, allocatable :: W(:,:,:), U(:,:,:), R(:,:) + double precision, allocatable :: y(:,:,:,:), h(:,:,:,:), lambda(:) + double precision :: diag_h_mat_elem + double precision :: residual_norm(N_st) + character*(16384) :: write_buffer + double precision :: to_print(2,N_st) + double precision :: cpu, wall + + PROVIDE det_connections + + call write_time(iunit) + call wall_time(wall) + call cpu_time(cpu) + write(iunit,'(A)') '' + write(iunit,'(A)') 'Davidson Diagonalization' + write(iunit,'(A)') '------------------------' + write(iunit,'(A)') '' + call write_int(iunit,N_st,'Number of states') + call write_int(iunit,sze,'Number of determinants') + write(iunit,'(A)') '' + write_buffer = '===== ' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ================' + enddo + write(iunit,'(A)') trim(write_buffer) + write_buffer = ' Iter' + do i=1,N_st + write_buffer = trim(write_buffer)//' Energy Residual' + enddo + write(iunit,'(A)') trim(write_buffer) + write_buffer = '===== ' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ================' + enddo + write(iunit,'(A)') trim(write_buffer) + + allocate( & + kl_pairs(2,N_st*(N_st+1)/2), & + W(sze,N_st,davidson_sze_max), & + U(sze,N_st,davidson_sze_max), & + R(sze,N_st), & + h(N_st,davidson_sze_max,N_st,davidson_sze_max), & + y(N_st,davidson_sze_max,N_st,davidson_sze_max), & + lambda(N_st*davidson_sze_max)) + + ASSERT (N_st > 0) + ASSERT (sze > 0) + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + + ! Initialization + ! ============== + + k_pairs=0 + do l=1,N_st + do k=1,l + k_pairs+=1 + kl_pairs(1,k_pairs) = k + kl_pairs(2,k_pairs) = l + enddo + enddo + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP SHARED(U,sze,N_st,overlap,kl_pairs,k_pairs, & + !$OMP Nint,dets_in,u_in) & + !$OMP PRIVATE(k,l,kl,i) + + + ! Orthonormalize initial guess + ! ============================ + + !$OMP DO + do kl=1,k_pairs + k = kl_pairs(1,kl) + l = kl_pairs(2,kl) + if (k/=l) then + overlap(k,l) = u_dot_v(U_in(1,k),U_in(1,l),sze) + overlap(l,k) = overlap(k,l) + else + overlap(k,k) = u_dot_u(U_in(1,k),sze) + endif + enddo + !$OMP END DO + !$OMP END PARALLEL + + call ortho_lowdin(overlap,size(overlap,1),N_st,U_in,size(U_in,1),sze) + + ! Davidson iterations + ! =================== + + converged = .False. + + do while (.not.converged) + + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(k,i) SHARED(U,u_in,sze,N_st) + do k=1,N_st + !$OMP DO + do i=1,sze + U(i,k,1) = u_in(i,k) + enddo + !$OMP END DO + enddo + !$OMP END PARALLEL + + do iter=1,davidson_sze_max-1 + + ! Compute W_k = H |u_k> + ! ---------------------- + + do k=1,N_st + call H_u_0(W(1,k,iter),U(1,k,iter),H_jj,sze,dets_in,Nint) + enddo + + ! Compute h_kl = = + ! ------------------------------------------- + + do l=1,N_st + do k=1,N_st + do iter2=1,iter-1 + h(k,iter2,l,iter) = u_dot_v(U(1,k,iter2),W(1,l,iter),sze) + h(k,iter,l,iter2) = h(k,iter2,l,iter) + enddo + enddo + do k=1,l + h(k,iter,l,iter) = u_dot_v(U(1,k,iter),W(1,l,iter),sze) + h(l,iter,k,iter) = h(k,iter,l,iter) + enddo + enddo + + !DEBUG H MATRIX + !do i=1,iter + ! print '(10(x,F16.10))', h(1,i,1,1:i) + !enddo + !print *, '' + !END + + ! Diagonalize h + ! ------------- + call lapack_diag(lambda,y,h,N_st*davidson_sze_max,N_st*iter) + + ! Express eigenvectors of h in the determinant basis + ! -------------------------------------------------- + + do k=1,N_st + do i=1,sze + U(i,k,iter+1) = 0.d0 + W(i,k,iter+1) = 0.d0 + do l=1,N_st + do iter2=1,iter + U(i,k,iter+1) = U(i,k,iter+1) + U(i,l,iter2)*y(l,iter2,k,1) + W(i,k,iter+1) = W(i,k,iter+1) + W(i,l,iter2)*y(l,iter2,k,1) + enddo + enddo + enddo + enddo + + ! Compute residual vector + ! ----------------------- + + do k=1,N_st + do i=1,sze + R(i,k) = lambda(k) * U(i,k,iter+1) - W(i,k,iter+1) + enddo + residual_norm(k) = u_dot_u(R(1,k),sze) + to_print(1,k) = lambda(k) + nuclear_repulsion + to_print(2,k) = residual_norm(k) + enddo + + write(iunit,'(X,I3,X,100(X,F16.10,X,E16.6))'), iter, to_print(:,1:N_st) + call davidson_converged(lambda,residual_norm,wall,iter,cpu,N_st,converged) + if (converged) then + exit + endif + + + ! Davidson step + ! ------------- + + do k=1,N_st + do i=1,sze + U(i,k,iter+1) = -1.d0/max(H_jj(i) - lambda(k),1.d-2) * R(i,k) + enddo + enddo + + ! Gram-Schmidt + ! ------------ + + double precision :: c + do k=1,N_st + do iter2=1,iter + do l=1,N_st + c = u_dot_v(U(1,k,iter+1),U(1,l,iter2),sze) + do i=1,sze + U(i,k,iter+1) -= c * U(i,l,iter2) + enddo + enddo + enddo + do l=1,k-1 + c = u_dot_v(U(1,k,iter+1),U(1,l,iter+1),sze) + do i=1,sze + U(i,k,iter+1) -= c * U(i,l,iter+1) + enddo + enddo + call normalize( U(1,k,iter+1), sze ) + enddo + + !DEBUG : CHECK OVERLAP + !print *, '===' + !do k=1,iter+1 + ! do l=1,k + ! c = u_dot_v(U(1,1,k),U(1,1,l),sze) + ! print *, k,l, c + ! enddo + !enddo + !print *, '===' + !pause + !END DEBUG + + + enddo + + if (.not.converged) then + iter = davidson_sze_max-1 + endif + + ! Re-contract to u_in + ! ----------- + + do k=1,N_st + energies(k) = lambda(k) + do i=1,sze + u_in(i,k) = 0.d0 + do iter2=1,iter + do l=1,N_st + u_in(i,k) += U(i,l,iter2)*y(l,iter2,k,1) + enddo + enddo + enddo + enddo + + enddo + + write_buffer = '===== ' + do i=1,N_st + write_buffer = trim(write_buffer)//' ================ ================' + enddo + write(iunit,'(A)') trim(write_buffer) + write(iunit,'(A)') '' + call write_time(iunit) + + deallocate ( & + kl_pairs, & + W, & + U, & + R, & + h, & + y, & + lambda & + ) + abort_here = abort_all +end + + BEGIN_PROVIDER [ character(64), davidson_criterion ] +&BEGIN_PROVIDER [ double precision, davidson_threshold ] + implicit none + BEGIN_DOC + ! Can be : [ energy | residual | both | wall_time | cpu_time | iterations ] + END_DOC + davidson_criterion = 'residual' + davidson_threshold = 1.d-6 +END_PROVIDER + +subroutine davidson_converged(energy,residual,wall,iterations,cpu,N_st,converged) + implicit none + BEGIN_DOC +! True if the Davidson algorithm is converged + END_DOC + integer, intent(in) :: N_st, iterations + logical, intent(out) :: converged + double precision, intent(in) :: energy(N_st), residual(N_st) + double precision, intent(in) :: wall, cpu + double precision :: E(N_st), time + double precision, allocatable, save :: energy_old(:) + + if (.not.allocated(energy_old)) then + allocate(energy_old(N_st)) + energy_old = 0.d0 + endif + + E = energy - energy_old + energy_old = energy + if (davidson_criterion == 'energy') then + converged = dabs(maxval(E(1:N_st))) < davidson_threshold + else if (davidson_criterion == 'residual') then + converged = dabs(maxval(residual(1:N_st))) < davidson_threshold + else if (davidson_criterion == 'both') then + converged = dabs(maxval(residual(1:N_st))) + dabs(maxval(E(1:N_st)) ) & + < davidson_threshold + else if (davidson_criterion == 'wall_time') then + call wall_time(time) + converged = time - wall > davidson_threshold + else if (davidson_criterion == 'cpu_time') then + call cpu_time(time) + converged = time - cpu > davidson_threshold + else if (davidson_criterion == 'iterations') then + converged = iterations >= int(davidson_threshold) + endif + converged = converged.or.abort_here +end diff --git a/src/Determinants/density_matrix.irp.f b/src/Determinants/density_matrix.irp.f new file mode 100644 index 00000000..f72b337c --- /dev/null +++ b/src/Determinants/density_matrix.irp.f @@ -0,0 +1,214 @@ + BEGIN_PROVIDER [ double precision, one_body_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix for each state + END_DOC + + integer :: j,k,l,m + integer :: occ(N_int*bit_kind_size,2) + double precision :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: exc(0:2,2,2),n_occ_alpha + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + + if(only_single_double_dm)then + print*,'ONLY DOUBLE DM' + one_body_dm_mo_alpha = one_body_single_double_dm_mo_alpha + one_body_dm_mo_beta = one_body_single_double_dm_mo_beta + else + one_body_dm_mo_alpha = 0.d0 + one_body_dm_mo_beta = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & + !$OMP tmp_a, tmp_b, n_occ_alpha)& + !$OMP SHARED(psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& + !$OMP elec_beta_num,one_body_dm_mo_alpha,one_body_dm_mo_beta,N_det,mo_tot_num_align,& + !$OMP mo_tot_num) + allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) + tmp_a = 0.d0 + tmp_b = 0.d0 + !$OMP DO SCHEDULE(dynamic) + do k=1,N_det + call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int) + call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int) + do m=1,N_states + ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j) += ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j) += ck + enddo + enddo + do l=1,k-1 + call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m) + if (s1==1) then + tmp_a(h1,p1) += ckl + tmp_a(p1,h1) += ckl + else + tmp_b(h1,p1) += ckl + tmp_b(p1,h1) += ckl + endif + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_body_dm_mo_alpha = one_body_dm_mo_alpha + tmp_a + !$OMP END CRITICAL + !$OMP CRITICAL + one_body_dm_mo_beta = one_body_dm_mo_beta + tmp_b + !$OMP END CRITICAL + deallocate(tmp_a,tmp_b) + !$OMP BARRIER + !$OMP END PARALLEL + + endif +END_PROVIDER + + BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_alpha, (mo_tot_num_align,mo_tot_num) ] +&BEGIN_PROVIDER [ double precision, one_body_single_double_dm_mo_beta, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! Alpha and beta one-body density matrix for each state + END_DOC + + integer :: j,k,l,m + integer :: occ(N_int*bit_kind_size,2) + double precision :: ck, cl, ckl + double precision :: phase + integer :: h1,h2,p1,p2,s1,s2, degree + integer :: exc(0:2,2,2),n_occ_alpha + double precision, allocatable :: tmp_a(:,:), tmp_b(:,:) + integer :: degree_respect_to_HF_k + integer :: degree_respect_to_HF_l + + PROVIDE elec_alpha_num elec_beta_num + + one_body_single_double_dm_mo_alpha = 0.d0 + one_body_single_double_dm_mo_beta = 0.d0 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(j,k,l,m,occ,ck, cl, ckl,phase,h1,h2,p1,p2,s1,s2, degree,exc, & + !$OMP tmp_a, tmp_b, n_occ_alpha,degree_respect_to_HF_k,degree_respect_to_HF_l)& + !$OMP SHARED(ref_bitmask,psi_det,psi_coef,N_int,N_states,state_average_weight,elec_alpha_num,& + !$OMP elec_beta_num,one_body_single_double_dm_mo_alpha,one_body_single_double_dm_mo_beta,N_det,mo_tot_num_align,& + !$OMP mo_tot_num) + allocate(tmp_a(mo_tot_num_align,mo_tot_num), tmp_b(mo_tot_num_align,mo_tot_num) ) + tmp_a = 0.d0 + tmp_b = 0.d0 + !$OMP DO SCHEDULE(dynamic) + do k=1,N_det + call bitstring_to_list(psi_det(1,1,k), occ(1,1), n_occ_alpha, N_int) + call bitstring_to_list(psi_det(1,2,k), occ(1,2), n_occ_alpha, N_int) + call get_excitation_degree(ref_bitmask,psi_det(1,1,k),degree_respect_to_HF_k,N_int) + + do m=1,N_states + ck = psi_coef(k,m)*psi_coef(k,m) * state_average_weight(m) + call get_excitation_degree(ref_bitmask,psi_det(1,1,k),degree_respect_to_HF_l,N_int) + if(degree_respect_to_HF_l.le.0)then + do l=1,elec_alpha_num + j = occ(l,1) + tmp_a(j,j) += ck + enddo + do l=1,elec_beta_num + j = occ(l,2) + tmp_b(j,j) += ck + enddo + endif + enddo + do l=1,k-1 + call get_excitation_degree(ref_bitmask,psi_det(1,1,l),degree_respect_to_HF_l,N_int) + if(degree_respect_to_HF_k.ne.0)cycle + if(degree_respect_to_HF_l.eq.2.and.degree_respect_to_HF_k.ne.2)cycle + call get_excitation_degree(psi_det(1,1,k),psi_det(1,1,l),degree,N_int) + if (degree /= 1) then + cycle + endif + call get_mono_excitation(psi_det(1,1,k),psi_det(1,1,l),exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + do m=1,N_states + ckl = psi_coef(k,m) * psi_coef(l,m) * phase * state_average_weight(m) + if (s1==1) then + tmp_a(h1,p1) += ckl + tmp_a(p1,h1) += ckl + else + tmp_b(h1,p1) += ckl + tmp_b(p1,h1) += ckl + endif + enddo + enddo + enddo + !$OMP END DO NOWAIT + !$OMP CRITICAL + one_body_single_double_dm_mo_alpha = one_body_single_double_dm_mo_alpha + tmp_a + !$OMP END CRITICAL + !$OMP CRITICAL + one_body_single_double_dm_mo_beta = one_body_single_double_dm_mo_beta + tmp_b + !$OMP END CRITICAL + deallocate(tmp_a,tmp_b) + !$OMP BARRIER + !$OMP END PARALLEL +END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_body_dm_mo, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! One-body density matrix + END_DOC + one_body_dm_mo = one_body_dm_mo_alpha + one_body_dm_mo_beta +END_PROVIDER + +BEGIN_PROVIDER [ double precision, one_body_spin_density_mo, (mo_tot_num_align,mo_tot_num) ] + implicit none + BEGIN_DOC + ! rho(alpha) - rho(beta) + END_DOC + one_body_spin_density_mo = one_body_dm_mo_alpha - one_body_dm_mo_beta +END_PROVIDER + +subroutine set_natural_mos + implicit none + BEGIN_DOC + ! Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + END_DOC + character*(64) :: label + double precision, allocatable :: tmp(:,:) + allocate(tmp(size(one_body_dm_mo,1),size(one_body_dm_mo,2))) + + ! Negation to have the occupied MOs first after the diagonalization + tmp = -one_body_dm_mo + label = "Natural" + call mo_as_eigvectors_of_mo_matrix(tmp,size(tmp,1),size(tmp,2),label) + deallocate(tmp) + +end +subroutine save_natural_mos + implicit none + BEGIN_DOC + ! Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis + END_DOC + call set_natural_mos + call save_mos + +end + + +BEGIN_PROVIDER [ double precision, state_average_weight, (N_states) ] + implicit none + BEGIN_DOC + ! Weights in the state-average calculation of the density matrix + END_DOC + state_average_weight = 1.d0/dble(N_states) +END_PROVIDER + diff --git a/src/Determinants/det_svd.irp.f b/src/Determinants/det_svd.irp.f new file mode 100644 index 00000000..0a57acf3 --- /dev/null +++ b/src/Determinants/det_svd.irp.f @@ -0,0 +1,61 @@ +program det_svd + implicit none + BEGIN_DOC +! Computes the SVD of the Alpha x Beta determinant coefficient matrix + END_DOC + integer :: i,j,k + + read_wf = .True. + TOUCH read_wf + + print *, 'SVD matrix before filling' + print *, '=========================' + print *, '' + print *, 'N_det = ', N_det + print *, 'N_det_alpha = ', N_det_alpha_unique + print *, 'N_det_beta = ', N_det_beta_unique + print *, '' + +! do i=1,N_det_alpha_unique +! do j=1,N_det_beta_unique +! print *, i,j,psi_svd_matrix(i,j,:) +! enddo +! enddo + + print *, '' + print *, 'Energy = ', ci_energy + print *, '' + + print *, psi_svd_coefs(1:20,1) + + call generate_all_alpha_beta_det_products + print *, '' + print *, 'Energy = ', ci_energy + print *, '' + + print *, 'SVD matrix after filling' + print *, '========================' + print *, '' + print *, 'N_det = ', N_det + print *, 'N_det_alpha = ', N_det_alpha_unique + print *, 'N_det_beta = ', N_det_beta_unique + print *, '' + print *, '' + call diagonalize_ci + print *, 'Energy = ', ci_energy + + do i=1,N_det_alpha_unique + do j=1,N_det_beta_unique + do k=1,N_states + if (dabs(psi_svd_matrix(i,j,k)) < 1.d-15) then + psi_svd_matrix(i,j,k) = 0.d0 + endif + enddo + enddo + enddo + + print *, '' + print *, psi_svd_coefs(1:20,1) + call save_wavefunction + +end diff --git a/src/Determinants/determinants.irp.f b/src/Determinants/determinants.irp.f index 03315836..a70d0fe8 100644 --- a/src/Determinants/determinants.irp.f +++ b/src/Determinants/determinants.irp.f @@ -68,9 +68,6 @@ BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ] ! The wave function determinants. Initialized with Hartree-Fock if the EZFIO file ! is empty END_DOC - - PROVIDE ezfio_filename - integer :: i logical :: exists character*64 :: label @@ -237,8 +234,6 @@ BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states_diag) ] ! is empty END_DOC - PROVIDE ezfio_filename - integer :: i,k, N_int2 logical :: exists double precision, allocatable :: psi_coef_read(:,:) @@ -602,8 +597,6 @@ subroutine read_dets(det,Nint,Ndet) integer :: i,k equivalence (det_8, det_bk) - PROVIDE ezfio_filename - call ezfio_get_determinants_N_int(N_int2) ASSERT (N_int2 == Nint) call ezfio_get_determinants_bit_kind(k) @@ -672,8 +665,6 @@ subroutine save_wavefunction_general(ndet,nstates,psidet,dim_psicoef,psicoef) integer :: i,k PROVIDE progress_bar - PROVIDE ezfio_filename - call start_progress(7,'Saving wfunction',0.d0) progress_bar(1) = 1 diff --git a/src/Determinants/determinants_bitmasks.irp.f b/src/Determinants/determinants_bitmasks.irp.f new file mode 100644 index 00000000..8343fa84 --- /dev/null +++ b/src/Determinants/determinants_bitmasks.irp.f @@ -0,0 +1,57 @@ +use bitmasks + +integer, parameter :: hole_ = 1 +integer, parameter :: particle_ = 2 +integer, parameter :: hole2_ = 3 +integer, parameter :: particle2_= 4 + +BEGIN_PROVIDER [ integer, N_single_exc_bitmasks ] + implicit none + BEGIN_DOC + ! Number of single excitation bitmasks + END_DOC + N_single_exc_bitmasks = 1 + !TODO : Read from input! +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), single_exc_bitmask, (N_int, 2, N_single_exc_bitmasks) ] + implicit none + BEGIN_DOC + ! single_exc_bitmask(:,1,i) is the bitmask for holes + ! single_exc_bitmask(:,2,i) is the bitmask for particles + ! for a given couple of hole/particle excitations i. + END_DOC + + single_exc_bitmask(:,hole_,1) = HF_bitmask(:,1) + single_exc_bitmask(:,particle_,1) = not(HF_bitmask(:,2)) + !TODO : Read from input! +END_PROVIDER + + +BEGIN_PROVIDER [ integer, N_double_exc_bitmasks ] + implicit none + BEGIN_DOC + ! Number of double excitation bitmasks + END_DOC + N_double_exc_bitmasks = 1 + !TODO : Read from input! +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), double_exc_bitmask, (N_int, 4, N_double_exc_bitmasks) ] + implicit none + BEGIN_DOC + ! double_exc_bitmask(:,1,i) is the bitmask for holes of excitation 1 + ! double_exc_bitmask(:,2,i) is the bitmask for particles of excitation 1 + ! double_exc_bitmask(:,3,i) is the bitmask for holes of excitation 2 + ! double_exc_bitmask(:,4,i) is the bitmask for particles of excitation 2 + ! for a given couple of hole/particle excitations i. + END_DOC + + double_exc_bitmask(:,hole_,1) = HF_bitmask(:,1) + double_exc_bitmask(:,particle_,1) = not(HF_bitmask(:,2)) + double_exc_bitmask(:,hole2_,1) = HF_bitmask(:,1) + double_exc_bitmask(:,particle2_,1) = not(HF_bitmask(:,2)) + + !TODO : Read from input! +END_PROVIDER + diff --git a/src/Determinants/diagonalize_CI.irp.f b/src/Determinants/diagonalize_CI.irp.f new file mode 100644 index 00000000..0e697ab3 --- /dev/null +++ b/src/Determinants/diagonalize_CI.irp.f @@ -0,0 +1,109 @@ +BEGIN_PROVIDER [ character*(64), diag_algorithm ] + implicit none + BEGIN_DOC + ! Diagonalization algorithm (Davidson or Lapack) + END_DOC + if (N_det > N_det_max_jacobi) then + diag_algorithm = "Davidson" + else + diag_algorithm = "Lapack" + endif + + if (N_det < N_states_diag) then + diag_algorithm = "Lapack" + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, CI_energy, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states lowest eigenvalues of the CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(output_determinants) + do j=1,N_states_diag + CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion + write(st,'(I4)') j + call write_double(output_determinants,CI_energy(j),'Energy of state '//trim(st)) + call write_double(output_determinants,CI_eigenvectors_s2(j),'S^2 of state '//trim(st)) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2, (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(i,j) = psi_coef(i,j) + enddo + enddo + + if (diag_algorithm == "Davidson") then + + call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, & + size(CI_eigenvectors,1),N_det,N_states_diag,N_int,output_determinants) + + else if (diag_algorithm == "Lapack") then + + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) + CI_electronic_energy(:) = 0.d0 + do i=1,N_det + CI_eigenvectors(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(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy(i_state) = eigenvalues(j) + CI_eigenvectors_s2(i_state) = s2 + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo +! if(i_state < min(N_states_diag,N_det))then +! print *, 'pb with the number of states' +! print *, 'i_state = ',i_state +! print *, 'N_states_diag ',N_states_diag +! print *,'stopping ...' +! stop +! endif + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + +subroutine diagonalize_CI + implicit none + BEGIN_DOC +! Replace the coefficients of the CI states by the coefficients of the +! eigenstates of the CI matrix + END_DOC + integer :: i,j + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors(i,j) + enddo + enddo + SOFT_TOUCH psi_coef CI_electronic_energy CI_energy CI_eigenvectors CI_eigenvectors_s2 +end diff --git a/src/Determinants/diagonalize_CI_SC2.irp.f b/src/Determinants/diagonalize_CI_SC2.irp.f new file mode 100644 index 00000000..3b0d7904 --- /dev/null +++ b/src/Determinants/diagonalize_CI_SC2.irp.f @@ -0,0 +1,59 @@ +BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states_diag) ] + implicit none + BEGIN_DOC + ! N_states_diag lowest eigenvalues of the CI matrix + END_DOC + + integer :: j + character*(8) :: st + call write_time(output_determinants) + do j=1,N_states_diag + CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion + write(st,'(I4)') j + call write_double(output_determinants,CI_SC2_energy(j),'Energy of state '//trim(st)) + enddo + +END_PROVIDER + + BEGIN_PROVIDER [ double precision, threshold_convergence_SC2] + implicit none + BEGIN_DOC + ! convergence of the correlation energy of SC2 iterations + END_DOC + threshold_convergence_SC2 = 1.d-10 + + END_PROVIDER + BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,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_SC2_eigenvectors(i,j) = psi_coef(i,j) + enddo +! TODO : check comment +! CI_SC2_electronic_energy(j) = CI_electronic_energy(j) + enddo + + call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, & + size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2) +END_PROVIDER + +subroutine diagonalize_CI_SC2 + implicit none + BEGIN_DOC +! Replace the coefficients of the CI states_diag by the coefficients of the +! eigenstates of the CI matrix + END_DOC + integer :: i,j + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_SC2_eigenvectors(i,j) + enddo + enddo + SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors +end diff --git a/src/Determinants/diagonalize_CI_mono.irp.f b/src/Determinants/diagonalize_CI_mono.irp.f new file mode 100644 index 00000000..1c9a4de3 --- /dev/null +++ b/src/Determinants/diagonalize_CI_mono.irp.f @@ -0,0 +1,72 @@ + BEGIN_PROVIDER [ double precision, CI_electronic_energy_mono, (N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_mono, (N_det,N_states_diag) ] +&BEGIN_PROVIDER [ double precision, CI_eigenvectors_s2_mono, (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_mono(i,j) = psi_coef(i,j) + enddo + enddo + + if (diag_algorithm == "Davidson") then + + call davidson_diag(psi_det,CI_eigenvectors_mono,CI_electronic_energy, & + size(CI_eigenvectors_mono,1),N_det,N_states_diag,N_int,output_determinants) + + else if (diag_algorithm == "Lapack") then + + double precision, allocatable :: eigenvectors(:,:), eigenvalues(:) + allocate (eigenvectors(size(H_matrix_all_dets,1),N_det)) + allocate (eigenvalues(N_det)) + call lapack_diag(eigenvalues,eigenvectors, & + H_matrix_all_dets,size(H_matrix_all_dets,1),N_det) + CI_electronic_energy_mono(:) = 0.d0 + do i=1,N_det + CI_eigenvectors_mono(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 + print*,'j = ',j + print*,'e = ',eigenvalues(j) + print*,'c = ',dabs(eigenvectors(1,j)) + if(dabs(eigenvectors(1,j)).gt.0.9d0)then + i_state += 1 + do i=1,N_det + CI_eigenvectors_mono(i,i_state) = eigenvectors(i,j) + enddo + CI_electronic_energy_mono(i_state) = eigenvalues(j) + CI_eigenvectors_s2_mono(i_state) = s2 + endif + endif + if (i_state.ge.N_states_diag) then + exit + endif + enddo + deallocate(eigenvectors,eigenvalues) + endif + +END_PROVIDER + +subroutine diagonalize_CI_mono + implicit none + BEGIN_DOC +! Replace the coefficients of the CI states by the coefficients of the +! eigenstates of the CI matrix + END_DOC + integer :: i,j + do j=1,N_states_diag + do i=1,N_det + psi_coef(i,j) = CI_eigenvectors_mono(i,j) + enddo + enddo + SOFT_TOUCH psi_coef CI_electronic_energy_mono CI_eigenvectors_mono CI_eigenvectors_s2_mono +end diff --git a/src/Determinants/excitations_utils.irp.f b/src/Determinants/excitations_utils.irp.f new file mode 100644 index 00000000..46e38b08 --- /dev/null +++ b/src/Determinants/excitations_utils.irp.f @@ -0,0 +1,16 @@ +subroutine apply_mono(i_hole,i_particle,ispin_excit,key_in,Nint) + implicit none + integer, intent(in) :: i_hole,i_particle,ispin_excit,Nint + integer(bit_kind), intent(inout) :: key_in(Nint,2) + integer :: k,j + use bitmasks + ! hole + k = ishft(i_hole-1,-bit_kind_shift)+1 + j = i_hole-ishft(k-1,bit_kind_shift)-1 + key_in(k,ispin_excit) = ibclr(key_in(k,ispin_excit),j) + + k = ishft(i_particle-1,-bit_kind_shift)+1 + j = i_particle-ishft(k-1,bit_kind_shift)-1 + key_in(k,ispin_excit) = ibset(key_in(k,ispin_excit),j) + +end diff --git a/src/Determinants/filter_connected.irp.f b/src/Determinants/filter_connected.irp.f new file mode 100644 index 00000000..93a6ee7b --- /dev/null +++ b/src/Determinants/filter_connected.irp.f @@ -0,0 +1,611 @@ + +subroutine filter_connected(key1,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Filters out the determinants that are not connected by H + ! + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,l + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1))) & + + popcnt( xor( key1(1,2,i), key2(1,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do j=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(j,1,i), key2(j,1))) +& + popcnt(xor( key1(j,2,i), key2(j,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +end + + +subroutine filter_connected_sorted_ab(key1,key2,next,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Filters out the determinants that are not connected by H + ! returns the array idx which contains the index of the + ! determinants in the array key1 that interact + ! via the H operator with key2. + ! idx(0) is the number of determinants that interact with key1 + ! + ! Determinants are taken from the psi_det_sorted_ab array + END_DOC + integer, intent(in) :: Nint, sze + integer, intent(in) :: next(2,N_det) + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,l + integer :: degree_x2 + integer(bit_kind) :: det3_1(Nint,2), det3_2(Nint,2) + + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + call filter_3_highest_electrons( key2(1,1), det3_2(1,1), Nint) + if (Nint==1) then + + i = 1 + do while ( i<= sze ) + call filter_3_highest_electrons( key1(1,1,i), det3_1(1,1), Nint) + degree_x2 = popcnt( xor( det3_1(1,1), det3_2(1,1))) + if (degree_x2 > 4) then + i = next(1,i) + cycle + else + degree_x2 = popcnt( xor( key1(1,1,i), key2(1,1)) ) + if (degree_x2 <= 4) then + degree_x2 += popcnt( xor( key1(1,2,i), key2(1,2)) ) + if (degree_x2 <= 4) then + idx(l) = i + l += 1 + endif + endif + i += 1 + endif + enddo + + else + + print *, 'Not implemented', irp_here + stop 1 + + endif + idx(0) = l-1 +end + + + + +subroutine filter_connected_davidson(key1,key2,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Filters out the determinants that are not connected by H + ! returns the array idx which contains the index of the + ! determinants in the array key1 that interact + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + ! key1 should come from psi_det_sorted_ab. + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,j,k,l + integer :: degree_x2 + integer :: j_int, j_start + integer*8 :: itmp + + PROVIDE N_con_int det_connections + ASSERT (Nint > 0) + ASSERT (sze >= 0) + + l=1 + + if (Nint==1) then + + i = idx(0) + do j_int=1,N_con_int + itmp = det_connections(j_int,i) + do while (itmp /= 0_8) + j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5) + do j = j_start+1, min(j_start+32,i-1) + degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) + & + popcnt(xor( key1(1,2,j), key2(1,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = j + l = l+1 + endif + enddo + itmp = iand(itmp-1_8,itmp) + enddo + enddo + + else if (Nint==2) then + + + i = idx(0) + do j_int=1,N_con_int + itmp = det_connections(j_int,i) + do while (itmp /= 0_8) + j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5) + do j = j_start+1, min(j_start+32,i-1) + degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) + & + popcnt(xor( key1(2,1,j), key2(2,1))) + & + popcnt(xor( key1(1,2,j), key2(1,2))) + & + popcnt(xor( key1(2,2,j), key2(2,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = j + l = l+1 + endif + enddo + itmp = iand(itmp-1_8,itmp) + enddo + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + i = idx(0) + do j_int=1,N_con_int + itmp = det_connections(j_int,i) + do while (itmp /= 0_8) + j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5) + do j = j_start+1, min(j_start+32,i-1) + degree_x2 = popcnt(xor( key1(1,1,j), key2(1,1))) + & + popcnt(xor( key1(1,2,j), key2(1,2))) + & + popcnt(xor( key1(2,1,j), key2(2,1))) + & + popcnt(xor( key1(2,2,j), key2(2,2))) + & + popcnt(xor( key1(3,1,j), key2(3,1))) + & + popcnt(xor( key1(3,2,j), key2(3,2))) + if (degree_x2 > 4) then + cycle + else + idx(l) = j + l = l+1 + endif + enddo + itmp = iand(itmp-1_8,itmp) + enddo + enddo + + else + + !DIR$ LOOP COUNT (1000) + i = idx(0) + do j_int=1,N_con_int + itmp = det_connections(j_int,i) + do while (itmp /= 0_8) + j_start = ishft(j_int-1,11) + ishft(trailz(itmp),5) + do j = j_start+1, min(j_start+32,i-1) + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do k=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(k,1,j), key2(k,1))) +& + popcnt(xor( key1(k,2,j), key2(k,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + idx(l) = j + l = l+1 + endif + enddo + itmp = iand(itmp-1_8,itmp) + enddo + enddo + + endif + idx(0) = l-1 +end + +subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx) + use bitmasks + BEGIN_DOC + ! returns the array idx which contains the index of the + ! + ! determinants in the array key1 that interact + ! + ! via the H operator with key2. + ! + ! idx(0) is the number of determinants that interact with key1 + END_DOC + implicit none + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + + integer :: i,l,m + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sze > 0) + + l=1 + + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 > 4) then + cycle + else if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 > 4) then + cycle + else if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 > 4) then + cycle + else if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do m=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) +& + popcnt(xor( key1(m,2,i), key2(m,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 > 4) then + cycle + else if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +end + +subroutine filter_connected_i_H_psi0_SC2(key1,key2,Nint,sze,idx,idx_repeat) + use bitmasks + BEGIN_DOC + ! standard filter_connected_i_H_psi but returns in addition + ! + ! the array of the index of the non connected determinants to key1 + ! + ! in order to know what double excitation can be repeated on key1 + ! + ! idx_repeat(0) is the number of determinants that can be used + ! + ! to repeat the excitations + END_DOC + implicit none + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: idx(0:sze) + integer, intent(out) :: idx_repeat(0:sze) + + integer :: i,l,l_repeat,m + integer :: degree_x2 + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sze > 0) + + integer :: degree + degree = popcnt(xor( ref_bitmask(1,1), key2(1,1))) + & + popcnt(xor( ref_bitmask(1,2), key2(1,2))) + !DEC$ NOUNROLL + do m=2,Nint + degree = degree+ popcnt(xor( ref_bitmask(m,1), key2(m,1))) + & + popcnt(xor( ref_bitmask(m,2), key2(m,2))) + enddo + degree = ishft(degree,-1) + + l_repeat=1 + l=1 + if(degree == 2)then + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree_x2>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree_x2>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if(degree_x2>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + else if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do m=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) +& + popcnt(xor( key1(m,2,i), key2(m,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + elseif(degree_x2>6)then + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + endif + elseif(degree==1)then + if (Nint==1) then + + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))) + if (degree_x2 < 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree_x2 = 0 + !DEC$ LOOP COUNT MIN(4) + do m=1,Nint + degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1))) +& + popcnt(xor( key1(m,2,i), key2(m,2))) + if (degree_x2 > 4) then + exit + endif + enddo + if (degree_x2 <= 5) then + if(degree_x2 .ne. 0)then + idx(l) = i + l = l+1 + endif + else + idx_repeat(l_repeat) = i + l_repeat = l_repeat + 1 + endif + enddo + + endif + + else +! print*,'more than a double excitation, can not apply the ' +! print*,'SC2 dressing of the diagonal element .....' +! print*,'stop !!' +! print*,'degree = ',degree +! stop + idx(0) = 0 + idx_repeat(0) = 0 + endif + idx(0) = l-1 + idx_repeat(0) = l_repeat-1 +end + diff --git a/src/Determinants/guess_doublet.irp.f b/src/Determinants/guess_doublet.irp.f new file mode 100644 index 00000000..a44697c1 --- /dev/null +++ b/src/Determinants/guess_doublet.irp.f @@ -0,0 +1,79 @@ +program put_gess + use bitmasks + implicit none + integer :: i,j,N_det_tmp,N_states_tmp + integer :: list(N_int*bit_kind_size,2) + integer(bit_kind) :: string(N_int,2) + integer(bit_kind) :: psi_det_tmp(N_int,2,3) + double precision :: psi_coef_tmp(3,1) + + integer :: iorb,jorb,korb + print*,'which open shells ?' + read(5,*)iorb,jorb,korb + print*,iorb,jorb,korb + N_states= 1 + N_det= 3 + + + list = 0 + list(1,1) = 1 + list(1,2) = 1 + list(2,1) = 2 + list(2,2) = 2 + list(3,1) = iorb + list(4,1) = jorb + list(3,2) = korb + print*,'passed' + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + print*,'passed' + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + print*,'passed' + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,1) = string(i,j) + enddo + enddo + psi_coef(1,1) = 1.d0/dsqrt(3.d0) + + print*,'passed 1' + list = 0 + list(1,1) = 1 + list(1,2) = 1 + list(2,1) = 2 + list(2,2) = 2 + list(3,1) = iorb + list(4,1) = korb + list(3,2) = jorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,2) = string(i,j) + enddo + enddo + psi_coef(2,1) = 1.d0/dsqrt(3.d0) + + print*,'passed 2' + list = 0 + list(1,1) = 1 + list(1,2) = 1 + list(2,1) = 2 + list(2,2) = 2 + list(3,1) = korb + list(4,1) = jorb + list(3,2) = iorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,3) = string(i,j) + enddo + enddo + psi_coef(3,1) = 1.d0/dsqrt(3.d0) + print*,'passed 3' + + call save_wavefunction +end diff --git a/src/Determinants/guess_singlet.irp.f b/src/Determinants/guess_singlet.irp.f new file mode 100644 index 00000000..50f8dc4e --- /dev/null +++ b/src/Determinants/guess_singlet.irp.f @@ -0,0 +1,44 @@ +program put_gess + use bitmasks + implicit none + integer :: i,j,N_det_tmp,N_states_tmp + integer :: list(N_int*bit_kind_size,2) + integer(bit_kind) :: string(N_int,2) + integer(bit_kind) :: psi_det_tmp(N_int,2,2) + double precision :: psi_coef_tmp(2,1) + + integer :: iorb,jorb + print*,'which open shells ?' + read(5,*)iorb,jorb + N_states= 1 + N_det= 2 + + + list = 0 + list(1,1) = iorb + list(1,2) = jorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,1) = string(i,j) + enddo + enddo + psi_coef(1,1) = 1.d0/dsqrt(2.d0) + + list = 0 + list(1,1) = jorb + list(1,2) = iorb + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,2) = string(i,j) + enddo + enddo + psi_coef(2,1) = 1.d0/dsqrt(2.d0) + + call save_wavefunction +end diff --git a/src/Determinants/guess_triplet.irp.f b/src/Determinants/guess_triplet.irp.f new file mode 100644 index 00000000..77f88c3e --- /dev/null +++ b/src/Determinants/guess_triplet.irp.f @@ -0,0 +1,48 @@ +program put_gess + use bitmasks + implicit none + integer :: i,j,N_det_tmp,N_states_tmp + integer :: list(N_int*bit_kind_size,2) + integer(bit_kind) :: string(N_int,2) + integer(bit_kind) :: psi_det_tmp(N_int,2,2) + double precision :: psi_coef_tmp(2,1) + + integer :: iorb,jorb + print*,'which open shells ?' + read(5,*)iorb,jorb + N_states= 1 + N_det= 2 + print*,'iorb = ',iorb + print*,'jorb = ',jorb + + + list = 0 + list(1,1) = iorb + list(1,2) = jorb + string = 0 + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,1) = string(i,j) + enddo + enddo + psi_coef(1,1) = 1.d0/dsqrt(2.d0) + + list = 0 + list(1,1) = jorb + list(1,2) = iorb + string = 0 + call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int) + call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int) + call print_det(string,N_int) + do j = 1,2 + do i = 1, N_int + psi_det(i,j,2) = string(i,j) + enddo + enddo + psi_coef(2,1) = -1.d0/dsqrt(2.d0) + + call save_wavefunction +end diff --git a/src/Determinants/occ_pattern.irp.f b/src/Determinants/occ_pattern.irp.f new file mode 100644 index 00000000..a0fd4a3c --- /dev/null +++ b/src/Determinants/occ_pattern.irp.f @@ -0,0 +1,339 @@ +use bitmasks +subroutine det_to_occ_pattern(d,o,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Transform a determinant to an occupation pattern + END_DOC + integer ,intent(in) :: Nint + integer(bit_kind),intent(in) :: d(Nint,2) + integer(bit_kind),intent(out) :: o(Nint,2) + + integer :: k + + do k=1,Nint + o(k,1) = ieor(d(k,1),d(k,2)) + o(k,2) = iand(d(k,1),d(k,2)) + enddo +end + +subroutine occ_pattern_to_dets_size(o,sze,n_alpha,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Number of possible determinants for a given occ_pattern + END_DOC + integer ,intent(in) :: Nint, n_alpha + integer(bit_kind),intent(in) :: o(Nint,2) + integer, intent(out) :: sze + integer :: amax,bmax,k + double precision, external :: binom_func + + amax = n_alpha + bmax = 0 + do k=1,Nint + bmax += popcnt( o(k,1) ) + amax -= popcnt( o(k,2) ) + enddo + sze = int( min(binom_func(bmax, amax), 1.d8) ) + +end + +subroutine occ_pattern_to_dets(o,d,sze,n_alpha,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Generate all possible determinants for a give occ_pattern + END_DOC + integer ,intent(in) :: Nint, n_alpha + integer ,intent(inout) :: sze + integer(bit_kind),intent(in) :: o(Nint,2) + integer(bit_kind),intent(out) :: d(Nint,2,sze) + + integer :: i, k, nt, na, nd, amax + integer :: list_todo(n_alpha) + integer :: list_a(n_alpha) + + amax = n_alpha + do k=1,Nint + amax -= popcnt( o(k,2) ) + enddo + + call bitstring_to_list(o(1,1), list_todo, nt, Nint) + + na = 0 + nd = 0 + d = 0 + call rec_occ_pattern_to_dets(list_todo,nt,list_a,na,d,nd,sze,amax,Nint) + + sze = nd + + do i=1,nd + ! Doubly occupied orbitals + do k=1,Nint + d(k,1,i) = ior(d(k,1,i),o(k,2)) + d(k,2,i) = ior(d(k,2,i),o(k,2)) + enddo + enddo + +! !TODO DEBUG +! integer :: j,s +! do i=1,nd +! do j=1,i-1 +! na=0 +! do k=1,Nint +! if((d(k,1,j) /= d(k,1,i)).or. & +! (d(k,2,j) /= d(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( j== 0 ) then +! print *, 'det ',i,' and ',j,' equal:' +! call debug_det(d(1,1,j),Nint) +! call debug_det(d(1,1,i),Nint) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG +end + +recursive subroutine rec_occ_pattern_to_dets(list_todo,nt,list_a,na,d,nd,sze,amax,Nint) + use bitmasks + implicit none + + integer, intent(in) :: nt, sze, amax, Nint,na + integer,intent(inout) :: list_todo(nt) + integer, intent(inout) :: list_a(na+1),nd + integer(bit_kind),intent(inout) :: d(Nint,2,sze) + + if (na == amax) then + nd += 1 + if (na > 0) then + call list_to_bitstring( d(1,1,nd), list_a, na, Nint) + endif + if (nt > 0) then + call list_to_bitstring( d(1,2,nd), list_todo, nt, Nint) + endif + else + integer :: i, j, k + integer :: list_todo_tmp(nt) + do i=1,nt + if (na > 0) then + if (list_todo(i) < list_a(na)) then + cycle + endif + endif + list_a(na+1) = list_todo(i) + k=1 + do j=1,nt + if (i/=j) then + list_todo_tmp(k) = list_todo(j) + k += 1 + endif + enddo + call rec_occ_pattern_to_dets(list_todo_tmp,nt-1,list_a,na+1,d,nd,sze,amax,Nint) + enddo + endif + +end + + BEGIN_PROVIDER [ integer(bit_kind), psi_occ_pattern, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_occ_pattern ] + implicit none + BEGIN_DOC + ! array of the occ_pattern present in the wf + ! psi_occ_pattern(:,1,j) = jth occ_pattern of the wave function : represent all the single occupation + ! psi_occ_pattern(:,2,j) = jth occ_pattern of the wave function : represent all the double occupation + END_DOC + integer :: i,j,k + + ! create + do i = 1, N_det + do k = 1, N_int + psi_occ_pattern(k,1,i) = ieor(psi_det(k,1,i),psi_det(k,2,i)) + psi_occ_pattern(k,2,i) = iand(psi_det(k,1,i),psi_det(k,2,i)) + enddo + enddo + + ! Sort + integer, allocatable :: iorder(:) + integer*8, allocatable :: bit_tmp(:) + integer*8, external :: occ_pattern_search_key + integer(bit_kind), allocatable :: tmp_array(:,:,:) + logical,allocatable :: duplicate(:) + + + allocate ( iorder(N_det), duplicate(N_det), bit_tmp(N_det), tmp_array(N_int,2,psi_det_size) ) + + do i=1,N_det + iorder(i) = i + !$DIR FORCEINLINE + bit_tmp(i) = occ_pattern_search_key(psi_occ_pattern(1,1,i),N_int) + enddo + call i8sort(bit_tmp,iorder,N_det) + !DIR$ IVDEP + do i=1,N_det + do k=1,N_int + tmp_array(k,1,i) = psi_occ_pattern(k,1,iorder(i)) + tmp_array(k,2,i) = psi_occ_pattern(k,2,iorder(i)) + enddo + duplicate(i) = .False. + enddo + + i=1 + integer (bit_kind) :: occ_pattern_tmp + do i=1,N_det + duplicate(i) = .False. + enddo + + do i=1,N_det-1 + if (duplicate(i)) then + cycle + endif + j = i+1 + do while (bit_tmp(j)==bit_tmp(i)) + if (duplicate(j)) then + j+=1 + cycle + endif + duplicate(j) = .True. + do k=1,N_int + if ( (tmp_array(k,1,i) /= tmp_array(k,1,j)) & + .or. (tmp_array(k,2,i) /= tmp_array(k,2,j)) ) then + duplicate(j) = .False. + exit + endif + enddo + j+=1 + if (j>N_det) then + exit + endif + enddo + enddo + + N_occ_pattern=0 + do i=1,N_det + if (duplicate(i)) then + cycle + endif + N_occ_pattern += 1 + do k=1,N_int + psi_occ_pattern(k,1,N_occ_pattern) = tmp_array(k,1,i) + psi_occ_pattern(k,2,N_occ_pattern) = tmp_array(k,2,i) + enddo + enddo + + deallocate(iorder,duplicate,bit_tmp,tmp_array) +! !TODO DEBUG +! integer :: s +! do i=1,N_occ_pattern +! do j=i+1,N_occ_pattern +! s = 0 +! do k=1,N_int +! if((psi_occ_pattern(k,1,j) /= psi_occ_pattern(k,1,i)).or. & +! (psi_occ_pattern(k,2,j) /= psi_occ_pattern(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( s == 0 ) then +! print *, 'Error : occ ', j, 'already in wf' +! call debug_det(psi_occ_pattern(1,1,j),N_int) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG +END_PROVIDER + +subroutine make_s2_eigenfunction + implicit none + integer :: i,j,k + integer :: smax, s + integer(bit_kind), allocatable :: d(:,:,:), det_buffer(:,:,:) + integer :: N_det_new + integer, parameter :: bufsze = 1000 + logical, external :: is_in_wavefunction + +! !TODO DEBUG +! do i=1,N_det +! do j=i+1,N_det +! s = 0 +! do k=1,N_int +! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. & +! (psi_det(k,2,j) /= psi_det(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( s == 0 ) then +! print *, 'Error0: det ', j, 'already in wf' +! call debug_det(psi_det(1,1,j),N_int) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG + + allocate (d(N_int,2,1), det_buffer(N_int,2,bufsze) ) + smax = 1 + N_det_new = 0 + + do i=1,N_occ_pattern + call occ_pattern_to_dets_size(psi_occ_pattern(1,1,i),s,elec_alpha_num,N_int) + s += 1 + if (s > smax) then + deallocate(d) + allocate ( d(N_int,2,s) ) + smax = s + endif + call occ_pattern_to_dets(psi_occ_pattern(1,1,i),d,s,elec_alpha_num,N_int) + do j=1,s + if (.not. is_in_wavefunction( d(1,1,j), N_int, N_det)) then + N_det_new += 1 + do k=1,N_int + det_buffer(k,1,N_det_new) = d(k,1,j) + det_buffer(k,2,N_det_new) = d(k,2,j) + enddo + if (N_det_new == bufsze) then + call fill_H_apply_buffer_no_selection(bufsze,det_buffer,N_int,0) + N_det_new = 0 + endif + endif + enddo + enddo + + if (N_det_new > 0) then + call fill_H_apply_buffer_no_selection(N_det_new,det_buffer,N_int,0) + call copy_H_apply_buffer_to_wf + SOFT_TOUCH N_det psi_coef psi_det + endif + + deallocate(d,det_buffer) + + +! !TODO DEBUG +! do i=1,N_det +! do j=i+1,N_det +! s = 0 +! do k=1,N_int +! if((psi_det(k,1,j) /= psi_det(k,1,i)).or. & +! (psi_det(k,2,j) /= psi_det(k,2,i))) then +! s=1 +! exit +! endif +! enddo +! if ( s == 0 ) then +! print *, 'Error : det ', j, 'already in wf at ', i +! call debug_det(psi_det(1,1,j),N_int) +! stop +! endif +! enddo +! enddo +! !TODO DEBUG + call write_int(output_determinants,N_det_new, 'Added deteminants for S^2') + +end + diff --git a/src/Determinants/options.irp.f b/src/Determinants/options.irp.f new file mode 100644 index 00000000..d4283128 --- /dev/null +++ b/src/Determinants/options.irp.f @@ -0,0 +1,22 @@ +BEGIN_PROVIDER [ integer, N_states_diag ] + implicit none + BEGIN_DOC +! Number of states to consider for the diagonalization + END_DOC + + logical :: has + PROVIDE ezfio_filename + call ezfio_has_determinants_n_states_diag(has) + if (has) then + call ezfio_get_determinants_n_states_diag(N_states_diag) + else + N_states_diag = N_states + endif + + call write_time(output_determinants) + call write_int(output_determinants, N_states_diag, & + 'N_states_diag') + + +END_PROVIDER + diff --git a/src/Determinants/program_beginer_determinants.irp.f b/src/Determinants/program_beginer_determinants.irp.f new file mode 100644 index 00000000..6375af22 --- /dev/null +++ b/src/Determinants/program_beginer_determinants.irp.f @@ -0,0 +1,138 @@ +program pouet + implicit none + print*,'HF energy = ',ref_bitmask_energy + nuclear_repulsion + call routine + +end +subroutine routine + use bitmasks + implicit none + integer :: i,j,k,l + double precision :: hij,get_mo_bielec_integral + double precision :: hmono,h_bi_ispin,h_bi_other_spin + integer(bit_kind),allocatable :: key_tmp(:,:) + integer, allocatable :: occ(:,:) + integer :: n_occ_alpha, n_occ_beta + ! First checks + print*,'N_int = ',N_int + print*,'mo_tot_num = ',mo_tot_num + print*,'mo_tot_num / 64+1= ',mo_tot_num/64+1 + ! We print the HF determinant + do i = 1, N_int + print*,'ref_bitmask(i,1) = ',ref_bitmask(i,1) + print*,'ref_bitmask(i,2) = ',ref_bitmask(i,2) + enddo + print*,'' + print*,'Hartree Fock determinant ...' + call debug_det(ref_bitmask,N_int) + allocate(key_tmp(N_int,2)) + ! We initialize key_tmp to the Hartree Fock one + key_tmp = ref_bitmask + integer :: i_hole,i_particle,ispin,i_ok,other_spin + ! We do a mono excitation on the top of the HF determinant + write(*,*)'Enter the (hole, particle) couple for the mono excitation ...' + read(5,*)i_hole,i_particle +!!i_hole = 4 +!!i_particle = 20 + write(*,*)'Enter the ispin variable ...' + write(*,*)'ispin = 1 ==> alpha ' + write(*,*)'ispin = 2 ==> beta ' + read(5,*)ispin + if(ispin == 1)then + other_spin = 2 + else if(ispin == 2)then + other_spin = 1 + else + print*,'PB !! ' + print*,'ispin must be 1 or 2 !' + stop + endif +!!ispin = 1 + call do_mono_excitation(key_tmp,i_hole,i_particle,ispin,i_ok) + ! We check if it the excitation was possible with "i_ok" + if(i_ok == -1)then + print*,'i_ok = ',i_ok + print*,'You can not do this excitation because of Pauli principle ...' + print*,'check your hole particle couple, there must be something wrong ...' + stop + + endif + print*,'New det = ' + call debug_det(key_tmp,N_int) + call i_H_j(key_tmp,ref_bitmask,N_int,hij) + ! We calculate the H matrix element between the new determinant and HF + print*,' = ',hij + print*,'' + print*,'' + print*,'Recalculating it old school style ....' + print*,'' + print*,'' + ! We recalculate this old school style !!! + ! Mono electronic part + hmono = mo_mono_elec_integral(i_hole,i_particle) + print*,'' + print*,'Mono electronic part ' + print*,'' + print*,' = ',hmono + h_bi_ispin = 0.d0 + h_bi_other_spin = 0.d0 + print*,'' + print*,'Getting all the info for the calculation of the bi electronic part ...' + print*,'' + allocate (occ(N_int*bit_kind_size,2)) + ! We get the occupation of the alpha electrons in occ(:,1) + call bitstring_to_list(key_tmp(1,1), occ(1,1), n_occ_alpha, N_int) + print*,'n_occ_alpha = ',n_occ_alpha + print*,'elec_alpha_num = ',elec_alpha_num + ! We get the occupation of the beta electrons in occ(:,2) + call bitstring_to_list(key_tmp(1,2), occ(1,2), n_occ_beta, N_int) + print*,'n_occ_beta = ',n_occ_beta + print*,'elec_beta_num = ',elec_beta_num + ! We print the occupation of the alpha electrons + print*,'Alpha electrons !' + do i = 1, n_occ_alpha + print*,'i = ',i + print*,'occ(i,1) = ',occ(i,1) + enddo + ! We print the occupation of the beta electrons + print*,'Alpha electrons !' + do i = 1, n_occ_beta + print*,'i = ',i + print*,'occ(i,2) = ',occ(i,2) + enddo + integer :: exc(0:2,2,2),degree,h1,p1,h2,p2,s1,s2 + double precision :: phase + + call get_excitation_degree(key_tmp,ref_bitmask,degree,N_int) + print*,'degree = ',degree + call get_mono_excitation(ref_bitmask,key_tmp,exc,phase,N_int) + call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + print*,'h1 = ',h1 + print*,'p1 = ',p1 + print*,'s1 = ',s1 + print*,'phase = ',phase + do i = 1, elec_num_tab(ispin) + integer :: orb_occupied + orb_occupied = occ(i,ispin) + h_bi_ispin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) & + -get_mo_bielec_integral(i_hole,i_particle,orb_occupied,orb_occupied,mo_integrals_map) + enddo + print*,'h_bi_ispin = ',h_bi_ispin + + do i = 1, elec_num_tab(other_spin) + orb_occupied = occ(i,other_spin) + h_bi_other_spin += get_mo_bielec_integral(i_hole,orb_occupied,i_particle,orb_occupied,mo_integrals_map) + enddo + print*,'h_bi_other_spin = ',h_bi_other_spin + print*,'h_bi_ispin + h_bi_other_spin = ',h_bi_ispin + h_bi_other_spin + + print*,'Total matrix element = ',phase*(h_bi_ispin + h_bi_other_spin + hmono) +!i = 1 +!j = 1 +!k = 1 +!l = 1 +!hij = get_mo_bielec_integral(i,j,k,l,mo_integrals_map) +!print*,' = ',hij + + +end diff --git a/src/Determinants/psi_cas.irp.f b/src/Determinants/psi_cas.irp.f new file mode 100644 index 00000000..8ca081d6 --- /dev/null +++ b/src/Determinants/psi_cas.irp.f @@ -0,0 +1,114 @@ +use bitmasks + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_cas, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_cas ] + implicit none + BEGIN_DOC + ! CAS wave function, defined from the application of the CAS bitmask on the + ! determinants. idx_cas gives the indice of the CAS determinant in psi_det. + END_DOC + integer :: i, k, l + logical :: good + 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 = N_det_cas+1 + do k=1,N_int + psi_cas(k,1,N_det_cas) = psi_det(k,1,i) + psi_cas(k,2,N_det_cas) = psi_det(k,2,i) + enddo + idx_cas(N_det_cas) = i + do k=1,N_states + psi_cas_coef(N_det_cas,k) = psi_coef(i,k) + enddo + endif + enddo + call write_int(output_determinants,N_det_cas, 'Number of determinants in the CAS') + +END_PROVIDER + + + BEGIN_PROVIDER [ integer(bit_kind), psi_cas_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_cas_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! CAS determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_cas, psi_cas, psi_cas_coef, & + psi_cas_sorted_bit, psi_cas_coef_sorted_bit) + +END_PROVIDER + + + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_cas_coef, (psi_det_size,n_states) ] +&BEGIN_PROVIDER [ integer, idx_non_cas, (psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_non_cas ] + implicit none + BEGIN_DOC + ! Set of determinants which are not part of the CAS, defined from the application + ! of the CAS bitmask on the determinants. + ! idx_non_cas gives the indice of the determinant in psi_det. + END_DOC + integer :: i_non_cas,j,k + integer :: degree + logical :: in_cas + i_non_cas =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_non_cas += 1 + do j=1,N_int + psi_non_cas(j,1,i_non_cas) = psi_det(j,1,k) + psi_non_cas(j,2,i_non_cas) = psi_det(j,2,k) + enddo + do j=1,N_states + psi_non_cas_coef(i_non_cas,j) = psi_coef(k,j) + enddo + idx_non_cas(i_non_cas) = k + endif + enddo + N_det_non_cas = i_non_cas +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_non_cas_sorted_bit, (N_int,2,psi_det_size) ] +&BEGIN_PROVIDER [ double precision, psi_non_cas_coef_sorted_bit, (psi_det_size,N_states) ] + implicit none + BEGIN_DOC + ! CAS determinants sorted to accelerate the search of a random determinant in the wave + ! function. + END_DOC + call sort_dets_by_det_search_key(N_det_cas, psi_non_cas, psi_non_cas_coef, & + psi_non_cas_sorted_bit, psi_non_cas_coef_sorted_bit) + +END_PROVIDER + + + + + diff --git a/src/Determinants/ref_bitmask.irp.f b/src/Determinants/ref_bitmask.irp.f new file mode 100644 index 00000000..7f760562 --- /dev/null +++ b/src/Determinants/ref_bitmask.irp.f @@ -0,0 +1,57 @@ + BEGIN_PROVIDER [ double precision, ref_bitmask_energy ] +&BEGIN_PROVIDER [ double precision, mono_elec_ref_bitmask_energy ] +&BEGIN_PROVIDER [ double precision, kinetic_ref_bitmask_energy ] +&BEGIN_PROVIDER [ double precision, nucl_elec_ref_bitmask_energy ] +&BEGIN_PROVIDER [ double precision, bi_elec_ref_bitmask_energy ] + use bitmasks + implicit none + BEGIN_DOC + ! Energy of the reference bitmask used in Slater rules + END_DOC + + integer :: occ(N_int*bit_kind_size,2) + integer :: i,j + + call bitstring_to_list(ref_bitmask(1,1), occ(1,1), i, N_int) + call bitstring_to_list(ref_bitmask(1,2), occ(1,2), i, N_int) + + + ref_bitmask_energy = 0.d0 + mono_elec_ref_bitmask_energy = 0.d0 + kinetic_ref_bitmask_energy = 0.d0 + nucl_elec_ref_bitmask_energy = 0.d0 + bi_elec_ref_bitmask_energy = 0.d0 + + do i = 1, elec_beta_num + ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1)) + mo_mono_elec_integral(occ(i,2),occ(i,2)) + kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1)) + mo_kinetic_integral(occ(i,2),occ(i,2)) + nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1)) + mo_nucl_elec_integral(occ(i,2),occ(i,2)) + enddo + + do i = elec_beta_num+1,elec_alpha_num + ref_bitmask_energy += mo_mono_elec_integral(occ(i,1),occ(i,1)) + kinetic_ref_bitmask_energy += mo_kinetic_integral(occ(i,1),occ(i,1)) + nucl_elec_ref_bitmask_energy += mo_nucl_elec_integral(occ(i,1),occ(i,1)) + enddo + + do j= 1, elec_alpha_num + do i = j+1, elec_alpha_num + bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1)) + ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,1),occ(j,1)) + enddo + enddo + + do j= 1, elec_beta_num + do i = j+1, elec_beta_num + bi_elec_ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2)) + ref_bitmask_energy += mo_bielec_integral_jj_anti(occ(i,2),occ(j,2)) + enddo + do i= 1, elec_alpha_num + bi_elec_ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2)) + ref_bitmask_energy += mo_bielec_integral_jj(occ(i,1),occ(j,2)) + enddo + enddo + mono_elec_ref_bitmask_energy = kinetic_ref_bitmask_energy + nucl_elec_ref_bitmask_energy + +END_PROVIDER + diff --git a/src/Determinants/s2.irp.f b/src/Determinants/s2.irp.f new file mode 100644 index 00000000..cd1d9fda --- /dev/null +++ b/src/Determinants/s2.irp.f @@ -0,0 +1,106 @@ +subroutine get_s2(key_i,key_j,phase,Nint) + implicit none + use bitmasks + BEGIN_DOC +! Returns + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2) + integer(bit_kind), intent(in) :: key_j(Nint,2) + double precision, intent(out) :: phase + integer :: exc(0:2,2,2) + integer :: degree + double precision :: phase_spsm + integer :: nup, i + + phase = 0.d0 + !$FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case(2) + call get_double_excitation(key_i,key_j,exc,phase_spsm,Nint) + if (exc(0,1,1) == 1) then ! Mono alpha + mono-beta + if ( (exc(1,1,1) == exc(1,2,2)).and.(exc(1,1,2) == exc(1,2,1)) ) then + phase = -phase_spsm + endif + endif + case(0) + nup = 0 + do i=1,Nint + nup += popcnt(iand(xor(key_i(i,1),key_i(i,2)),key_i(i,1))) + enddo + phase = dble(nup) + end select +end + +BEGIN_PROVIDER [ double precision, S_z ] +&BEGIN_PROVIDER [ double precision, S_z2_Sz ] + implicit none + BEGIN_DOC +! z component of the Spin + END_DOC + + S_z = 0.5d0*dble(elec_alpha_num-elec_beta_num) + S_z2_Sz = S_z*(S_z-1.d0) + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, expected_s2] + implicit none + BEGIN_DOC +! Expected value of S2 : S*(S+1) + END_DOC + logical :: has_expected_s2 + + call ezfio_has_determinants_expected_s2(has_expected_s2) + if (has_expected_s2) then + call ezfio_get_determinants_expected_s2(expected_s2) + else + double precision :: S + S = (elec_alpha_num-elec_beta_num)*0.5d0 + expected_s2 = S * (S+1.d0) +! expected_s2 = elec_alpha_num - elec_beta_num + 0.5d0 * ((elec_alpha_num - elec_beta_num)**2*0.5d0 - (elec_alpha_num-elec_beta_num)) + endif + +END_PROVIDER + +BEGIN_PROVIDER [ double precision, s2_values, (N_states) ] + implicit none + BEGIN_DOC +! array of the averaged values of the S^2 operator on the various states + END_DOC + integer :: i + double precision :: s2 + do i = 1, N_states + call get_s2_u0(psi_det,psi_coef(1,i),n_det,psi_det_size,s2) + s2_values(i) = s2 + enddo + +END_PROVIDER + + +subroutine get_s2_u0(psi_keys_tmp,psi_coefs_tmp,n,nmax,s2) + implicit none + use bitmasks + integer(bit_kind), intent(in) :: psi_keys_tmp(N_int,2,nmax) + integer, intent(in) :: n,nmax + double precision, intent(in) :: psi_coefs_tmp(nmax) + double precision, intent(out) :: s2 + integer :: i,j,l + double precision :: s2_tmp + s2 = S_z2_Sz + !$OMP PARALLEL DO DEFAULT(NONE) & + !$OMP PRIVATE(i,j,s2_tmp) SHARED(n,psi_coefs_tmp,psi_keys_tmp,N_int) & + !$OMP REDUCTION(+:s2) SCHEDULE(dynamic) + do i = 1, n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,i),s2_tmp,N_int) +! print*,'s2_tmp = ',s2_tmp + do j = 1, n + call get_s2(psi_keys_tmp(1,1,i),psi_keys_tmp(1,1,j),s2_tmp,N_int) + if (s2_tmp == 0.d0) cycle + s2 += psi_coefs_tmp(i)*psi_coefs_tmp(j)*s2_tmp + enddo + enddo + !$OMP END PARALLEL DO +end + diff --git a/src/Determinants/save_for_casino.irp.f b/src/Determinants/save_for_casino.irp.f new file mode 100644 index 00000000..631f79bd --- /dev/null +++ b/src/Determinants/save_for_casino.irp.f @@ -0,0 +1,268 @@ +subroutine save_casino + use bitmasks + implicit none + character*(128) :: message + integer :: getUnitAndOpen, iunit + integer, allocatable :: itmp(:) + integer :: n_ao_new + real, allocatable :: rtmp(:) + PROVIDE ezfio_filename + + iunit = getUnitAndOpen('gwfn.data','w') + print *, 'Title?' + read(*,*) message + write(iunit,'(A)') trim(message) + write(iunit,'(A)') '' + write(iunit,'(A)') 'BASIC_INFO' + write(iunit,'(A)') '----------' + write(iunit,'(A)') 'Generated by:' + write(iunit,'(A)') 'Quantum package' + write(iunit,'(A)') 'Method:' + print *, 'Method?' + read(*,*) message + write(iunit,'(A)') trim(message) + write(iunit,'(A)') 'DFT Functional:' + write(iunit,'(A)') 'none' + write(iunit,'(A)') 'Periodicity:' + write(iunit,'(A)') '0' + write(iunit,'(A)') 'Spin unrestricted:' + write(iunit,'(A)') '.false.' + write(iunit,'(A)') 'nuclear-nuclear repulsion energy (au/atom):' + write(iunit,*) nuclear_repulsion + write(iunit,'(A)') 'Number of electrons per primitive cell:' + write(iunit,*) elec_num + write(iunit,*) '' + + + write(iunit,*) 'GEOMETRY' + write(iunit,'(A)') '--------' + write(iunit,'(A)') 'Number of atoms:' + write(iunit,*) nucl_num + write(iunit,'(A)') 'Atomic positions (au):' + integer :: i + do i=1,nucl_num + write(iunit,'(3(1PE20.13))') nucl_coord(i,1:3) + enddo + write(iunit,'(A)') 'Atomic numbers for each atom:' + ! Add 200 if pseudopotential + allocate(itmp(nucl_num)) + do i=1,nucl_num + itmp(i) = int(nucl_charge(i)) + enddo + write(iunit,'(8(I10))') itmp(1:nucl_num) + deallocate(itmp) + write(iunit,'(A)') 'Valence charges for each atom:' + write(iunit,'(4(1PE20.13))') nucl_charge(1:nucl_num) + write(iunit,'(A)') '' + + + write(iunit,'(A)') 'BASIS SET' + write(iunit,'(A)') '---------' + write(iunit,'(A)') 'Number of Gaussian centres' + write(iunit,*) nucl_num + write(iunit,'(A)') 'Number of shells per primitive cell' + integer :: icount + icount = 0 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + icount += 1 + endif + enddo + write(iunit,*) icount + write(iunit,'(A)') 'Number of basis functions (''AO'') per primitive cell' + icount = 0 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + icount += 2*ao_l(i)+1 + endif + enddo + n_ao_new = icount + write(iunit,*) n_ao_new + write(iunit,'(A)') 'Number of Gaussian primitives per primitive cell' + allocate(itmp(ao_num)) + integer :: l + l=0 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + l += 1 + itmp(l) = ao_prim_num(i) + endif + enddo + write(iunit,'(8(I10))') sum(itmp(1:l)) + write(iunit,'(A)') 'Highest shell angular momentum (s/p/d/f... 1/2/3/4...)' + write(iunit,*) maxval(ao_l(1:ao_num))+1 + write(iunit,'(A)') 'Code for shell types (s/sp/p/d/f... 1/2/3/4/5...)' + l=0 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + l += 1 + if (ao_l(i) > 0) then + itmp(l) = ao_l(i)+2 + else + itmp(l) = ao_l(i)+1 + endif + endif + enddo + write(iunit,'(8(I10))') itmp(1:l) + write(iunit,'(A)') 'Number of primitive Gaussians in each shell' + l=0 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + l += 1 + itmp(l) = ao_prim_num(i) + endif + enddo + write(iunit,'(8(I10))') itmp(1:l) + deallocate(itmp) + write(iunit,'(A)') 'Sequence number of first shell on each centre' + allocate(itmp(nucl_num)) + l=0 + icount = 1 + itmp(icount) = 1 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + l = l+1 + if (ao_nucl(i) == icount) then + continue + else if (ao_nucl(i) == icount+1) then + icount += 1 + itmp(icount) = l + else + print *, 'Problem in order of centers of basis functions' + stop 1 + endif + endif + enddo + ! Check + if (icount /= nucl_num) then + print *, 'Error :' + print *, ' icount :', icount + print *, ' nucl_num:', nucl_num + stop 2 + endif + write(iunit,'(8(I10))') itmp(1:nucl_num) + deallocate(itmp) + write(iunit,'(A)') 'Exponents of Gaussian primitives' + allocate(rtmp(ao_num)) + l=0 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + do j=1,ao_prim_num(i) + l+=1 + rtmp(l) = ao_expo(i,ao_prim_num(i)-j+1) + enddo + endif + enddo + write(iunit,'(4(1PE20.13))') rtmp(1:l) + write(iunit,'(A)') 'Normalized contraction coefficients' + l=0 + integer :: j + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + do j=1,ao_prim_num(i) + l+=1 + rtmp(l) = ao_coef(i,ao_prim_num(i)-j+1) + enddo + endif + enddo + write(iunit,'(4(1PE20.13))') rtmp(1:l) + deallocate(rtmp) + write(iunit,'(A)') 'Position of each shell (au)' + l=0 + do i=1,ao_num + if (ao_l(i) == ao_power(i,1)) then + write(iunit,'(3(1PE20.13))') nucl_coord( ao_nucl(i), 1:3 ) + endif + enddo + write(iunit,'(A)') + + + write(iunit,'(A)') 'MULTIDETERMINANT INFORMATION' + write(iunit,'(A)') '----------------------------' + write(iunit,'(A)') 'GS' + write(iunit,'(A)') 'ORBITAL COEFFICIENTS' + write(iunit,'(A)') '------------------------' + + ! Transformation cartesian -> spherical + double precision :: tf2(6,5), tf3(10,7), tf4(15,9) + integer :: check2(3,6), check3(3,10), check4(3,15) + check2(:,1) = (/ 2, 0, 0 /) + check2(:,2) = (/ 1, 1, 0 /) + check2(:,3) = (/ 1, 0, 1 /) + check2(:,4) = (/ 0, 2, 0 /) + check2(:,5) = (/ 0, 1, 1 /) + check2(:,6) = (/ 0, 0, 2 /) + + check3(:,1) = (/ 3, 0, 0 /) + check3(:,2) = (/ 2, 1, 0 /) + check3(:,3) = (/ 2, 0, 1 /) + check3(:,4) = (/ 1, 2, 0 /) + check3(:,5) = (/ 1, 1, 1 /) + check3(:,6) = (/ 1, 0, 2 /) + check3(:,7) = (/ 0, 3, 0 /) + check3(:,8) = (/ 0, 2, 1 /) + check3(:,9) = (/ 0, 1, 2 /) + check3(:,10) = (/ 0, 0, 3 /) + + check4(:,1) = (/ 4, 0, 0 /) + check4(:,2) = (/ 3, 1, 0 /) + check4(:,3) = (/ 3, 0, 1 /) + check4(:,4) = (/ 2, 2, 0 /) + check4(:,5) = (/ 2, 1, 1 /) + check4(:,6) = (/ 2, 0, 2 /) + check4(:,7) = (/ 1, 3, 0 /) + check4(:,8) = (/ 1, 2, 1 /) + check4(:,9) = (/ 1, 1, 2 /) + check4(:,10) = (/ 1, 0, 3 /) + check4(:,11) = (/ 0, 4, 0 /) + check4(:,12) = (/ 0, 3, 1 /) + check4(:,13) = (/ 0, 2, 2 /) + check4(:,14) = (/ 0, 1, 3 /) + check4(:,15) = (/ 0, 0, 4 /) + +! tf2 = (/ +! -0.5, 0, 0, -0.5, 0, 1.0, & +! 0, 0, 1.0, 0, 0, 0, & +! 0, 0, 0, 0, 1.0, 0, & +! 0.86602540378443864676, 0, 0, -0.86602540378443864676, 0, 0, & +! 0, 1.0, 0, 0, 0, 0, & +! /) +! tf3 = (/ +! 0, 0, -0.67082039324993690892, 0, 0, 0, 0, -0.67082039324993690892, 0, 1.0, & +! -0.61237243569579452455, 0, 0, -0.27386127875258305673, 0, 1.0954451150103322269, 0, 0, 0, 0, & +! 0, -0.27386127875258305673, 0, 0, 0, 0, -0.61237243569579452455, 0, 1.0954451150103322269, 0, & +! 0, 0, 0.86602540378443864676, 0, 0, 0, 0, -0.86602540378443864676, 0, 0, & +! 0, 0, 0, 0, 1.0, 0, 0, 0, 0, 0, & +! 0.790569415042094833, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, & +! 0, 1.0606601717798212866, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0, & +! /) +! tf4 = (/ +! 0.375, 0, 0, 0.21957751641341996535, 0, -0.87831006565367986142, 0, 0, 0, 0, 0.375, 0, -0.87831006565367986142, 0, 1.0, & +! 0, 0, -0.89642145700079522998, 0, 0, 0, 0, -0.40089186286863657703, 0, 1.19522860933439364, 0, 0, 0, 0, 0, & +! 0, 0, 0, 0, -0.40089186286863657703, 0, 0, 0, 0, 0, 0, -0.89642145700079522998, 0, 1.19522860933439364, 0, & +! -0.5590169943749474241, 0, 0, 0, 0, 0.9819805060619657157, 0, 0, 0, 0, 0.5590169943749474241, 0, -0.9819805060619657157, 0, 0, & +! 0, -0.42257712736425828875, 0, 0, 0, 0, -0.42257712736425828875, 0, 1.1338934190276816816, 0, 0, 0, 0, 0, 0, & +! 0, 0, 0.790569415042094833, 0, 0, 0, 0, -1.0606601717798212866, 0, 0, 0, 0, 0, 0, 0, & +! 0, 0, 0, 0, 1.0606601717798212866, 0, 0, 0, 0, 0, 0, -0.790569415042094833, 0, 0, 0, & +! 0.73950997288745200532, 0, 0, -1.2990381056766579701, 0, 0, 0, 0, 0, 0, 0.73950997288745200532, 0, 0, 0, 0, & +! 0, 1.1180339887498948482, 0, 0, 0, 0, -1.1180339887498948482, 0, 0, 0, 0, 0, 0, 0, 0, & +! /) +! + + + allocate(rtmp(ao_num*mo_tot_num)) + l=0 + do i=1,mo_tot_num + do j=1,ao_num + l += 1 + rtmp(l) = mo_coef(j,i) + enddo + enddo + write(iunit,'(4(1PE20.13))') rtmp(1:l) + deallocate(rtmp) + close(iunit) +end + +program prog_save_casino + call save_casino +end diff --git a/src/Determinants/save_for_qmcchem.irp.f b/src/Determinants/save_for_qmcchem.irp.f new file mode 100644 index 00000000..b707ff7c --- /dev/null +++ b/src/Determinants/save_for_qmcchem.irp.f @@ -0,0 +1,51 @@ +subroutine save_dets_qmcchem + use bitmasks + implicit none + character :: c(mo_tot_num) + integer :: i,k + + integer, allocatable :: occ(:,:,:), occ_tmp(:,:) + !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: occ, occ_tmp + + read_wf = .True. + TOUCH read_wf + call ezfio_set_determinants_det_num(N_det) + call ezfio_set_determinants_det_coef(psi_coef_sorted(1,1)) + + allocate (occ(elec_alpha_num,N_det,2)) + ! OMP PARALLEL DEFAULT(NONE) & + ! OMP PRIVATE(occ_tmp,i,k)& + ! OMP SHARED(N_det,psi_det_sorted,elec_alpha_num, & + ! OMP occ,elec_beta_num,N_int) + allocate (occ_tmp(N_int*bit_kind_size,2)) + occ_tmp = 0 + ! OMP DO + do i=1,N_det + call bitstring_to_list(psi_det_sorted(1,1,i), occ_tmp(1,1), elec_alpha_num, N_int ) + call bitstring_to_list(psi_det_sorted(1,2,i), occ_tmp(1,2), elec_beta_num, N_int ) + do k=1,elec_alpha_num + occ(k,i,1) = occ_tmp(k,1) + occ(k,i,2) = occ_tmp(k,2) + enddo + enddo + ! OMP END DO + deallocate(occ_tmp) + ! OMP END PARALLEL + call ezfio_set_determinants_det_occ(occ) + call write_int(output_determinants,N_det,'Determinants saved for QMC') + deallocate(occ) + open(unit=31,file=trim(ezfio_filename)//'/mo_basis/mo_classif') + write(31,'(I1)') 1 + write(31,*) mo_tot_num + do i=1,mo_tot_num + write(31,'(A)') 'a' + enddo + close(31) + call system('gzip -f '//trim(ezfio_filename)//'/mo_basis/mo_classif') + +end + +program save_for_qmc + call save_dets_qmcchem + call write_spindeterminants +end diff --git a/src/Determinants/save_natorb.irp.f b/src/Determinants/save_natorb.irp.f new file mode 100644 index 00000000..e56f9821 --- /dev/null +++ b/src/Determinants/save_natorb.irp.f @@ -0,0 +1,6 @@ +program save_natorb + read_wf = .True. + touch read_wf + call save_natural_mos +end + diff --git a/src/Determinants/slater_rules.irp.f b/src/Determinants/slater_rules.irp.f new file mode 100644 index 00000000..7d431879 --- /dev/null +++ b/src/Determinants/slater_rules.irp.f @@ -0,0 +1,1301 @@ +subroutine get_excitation_degree(key1,key2,degree,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the excitation degree between two determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key1(Nint,2) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: degree + + integer :: l + + ASSERT (Nint > 0) + + degree = popcnt(xor( key1(1,1), key2(1,1))) + & + popcnt(xor( key1(1,2), key2(1,2))) + !DEC$ NOUNROLL + do l=2,Nint + degree = degree+ popcnt(xor( key1(l,1), key2(l,1))) + & + popcnt(xor( key1(l,2), key2(l,2))) + enddo + ASSERT (degree >= 0) + degree = ishft(degree,-1) + +end + + + +subroutine get_excitation(det1,det2,exc,degree,phase,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the excitation operators between two determinants and the phase + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det1(Nint,2) + integer(bit_kind), intent(in) :: det2(Nint,2) + integer, intent(out) :: exc(0:2,2,2) + integer, intent(out) :: degree + double precision, intent(out) :: phase + ! exc(number,hole/particle,spin) + ! ex : + ! exc(0,1,1) = number of holes alpha + ! exc(0,2,1) = number of particle alpha + ! exc(0,2,2) = number of particle beta + ! exc(1,2,1) = first particle alpha + ! exc(1,1,1) = first hole alpha + ! exc(1,2,2) = first particle beta + ! exc(1,1,2) = first hole beta + + ASSERT (Nint > 0) + + !DIR$ FORCEINLINE + call get_excitation_degree(det1,det2,degree,Nint) + select case (degree) + + case (3:) + degree = -1 + return + + case (2) + call get_double_excitation(det1,det2,exc,phase,Nint) + return + + case (1) + call get_mono_excitation(det1,det2,exc,phase,Nint) + return + + case(0) + return + + end select +end + +subroutine decode_exc(exc,degree,h1,p1,h2,p2,s1,s2) + use bitmasks + implicit none + BEGIN_DOC + ! Decodes the exc arrays returned by get_excitation. + ! h1,h2 : Holes + ! p1,p2 : Particles + ! s1,s2 : Spins (1:alpha, 2:beta) + ! degree : Degree of excitation + END_DOC + integer, intent(in) :: exc(0:2,2,2),degree + integer, intent(out) :: h1,h2,p1,p2,s1,s2 + ASSERT (degree > 0) + ASSERT (degree < 3) + + select case(degree) + case(2) + if (exc(0,1,1) == 2) then + h1 = exc(1,1,1) + h2 = exc(2,1,1) + p1 = exc(1,2,1) + p2 = exc(2,2,1) + s1 = 1 + s2 = 1 + else if (exc(0,1,2) == 2) then + h1 = exc(1,1,2) + h2 = exc(2,1,2) + p1 = exc(1,2,2) + p2 = exc(2,2,2) + s1 = 2 + s2 = 2 + else + h1 = exc(1,1,1) + h2 = exc(1,1,2) + p1 = exc(1,2,1) + p2 = exc(1,2,2) + s1 = 1 + s2 = 2 + endif + case(1) + if (exc(0,1,1) == 1) then + h1 = exc(1,1,1) + h2 = 0 + p1 = exc(1,2,1) + p2 = 0 + s1 = 1 + s2 = 0 + else + h1 = exc(1,1,2) + h2 = 0 + p1 = exc(1,2,2) + p2 = 0 + s1 = 2 + s2 = 0 + endif + case(0) + h1 = 0 + p1 = 0 + h2 = 0 + p2 = 0 + s1 = 0 + s2 = 0 + end select +end + +subroutine get_double_excitation(det1,det2,exc,phase,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the two excitation operators between two doubly excited determinants and the phase + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det1(Nint,2) + integer(bit_kind), intent(in) :: det2(Nint,2) + integer, intent(out) :: exc(0:2,2,2) + double precision, intent(out) :: phase + integer :: tz + integer :: l, ispin, idx_hole, idx_particle, ishift + integer :: nperm + integer :: i,j,k,m,n + integer :: high, low + integer :: a,b,c,d + integer(bit_kind) :: hole, particle, tmp + double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /) + + ASSERT (Nint > 0) + nperm = 0 + exc(0,1,1) = 0 + exc(0,2,1) = 0 + exc(0,1,2) = 0 + exc(0,2,2) = 0 + do ispin = 1,2 + idx_particle = 0 + idx_hole = 0 + ishift = 1-bit_kind_size + do l=1,Nint + ishift = ishift + bit_kind_size + if (det1(l,ispin) == det2(l,ispin)) then + cycle + endif + tmp = xor( det1(l,ispin), det2(l,ispin) ) + particle = iand(tmp, det2(l,ispin)) + hole = iand(tmp, det1(l,ispin)) + do while (particle /= 0_bit_kind) + tz = trailz(particle) + idx_particle = idx_particle + 1 + exc(0,2,ispin) = exc(0,2,ispin) + 1 + exc(idx_particle,2,ispin) = tz+ishift + particle = iand(particle,particle-1_bit_kind) + enddo + if (iand(exc(0,1,ispin),exc(0,2,ispin))==2) then ! exc(0,1,ispin)==2 or exc(0,2,ispin)==2 + exit + endif + do while (hole /= 0_bit_kind) + tz = trailz(hole) + idx_hole = idx_hole + 1 + exc(0,1,ispin) = exc(0,1,ispin) + 1 + exc(idx_hole,1,ispin) = tz+ishift + hole = iand(hole,hole-1_bit_kind) + enddo + if (iand(exc(0,1,ispin),exc(0,2,ispin))==2) then ! exc(0,1,ispin)==2 or exc(0,2,ispin) + exit + endif + enddo + + ! TODO : Voir si il faut sortir i,n,k,m du case. + + select case (exc(0,1,ispin)) + case(0) + cycle + + case(1) + low = min(exc(1,1,ispin), exc(1,2,ispin)) + high = max(exc(1,1,ispin), exc(1,2,ispin)) + + ASSERT (low > 0) + j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint) + n = iand(low,bit_kind_size-1) ! mod(low,bit_kind_size) + ASSERT (high > 0) + k = ishft(high-1,-bit_kind_shift)+1 + m = iand(high,bit_kind_size-1) + + if (j==k) then + nperm = nperm + popcnt(iand(det1(j,ispin), & + iand( ibset(0_bit_kind,m-1)-1_bit_kind, & + ibclr(-1_bit_kind,n)+1_bit_kind ) )) + else + nperm = nperm + popcnt(iand(det1(k,ispin), & + ibset(0_bit_kind,m-1)-1_bit_kind)) + & + popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind)) + do i=j+1,k-1 + nperm = nperm + popcnt(det1(i,ispin)) + end do + endif + + case (2) + + do i=1,2 + low = min(exc(i,1,ispin), exc(i,2,ispin)) + high = max(exc(i,1,ispin), exc(i,2,ispin)) + + ASSERT (low > 0) + j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint) + n = iand(low,bit_kind_size-1) ! mod(low,bit_kind_size) + ASSERT (high > 0) + k = ishft(high-1,-bit_kind_shift)+1 + m = iand(high,bit_kind_size-1) + + if (j==k) then + nperm = nperm + popcnt(iand(det1(j,ispin), & + iand( ibset(0_bit_kind,m-1)-1_bit_kind, & + ibclr(-1_bit_kind,n)+1_bit_kind ) )) + else + nperm = nperm + popcnt(iand(det1(k,ispin), & + ibset(0_bit_kind,m-1)-1_bit_kind)) + & + popcnt(iand(det1(j,ispin), ibclr(-1_bit_kind,n) +1_bit_kind)) + do l=j+1,k-1 + nperm = nperm + popcnt(det1(l,ispin)) + end do + endif + + enddo + + a = min(exc(1,1,ispin), exc(1,2,ispin)) + b = max(exc(1,1,ispin), exc(1,2,ispin)) + c = min(exc(2,1,ispin), exc(2,2,ispin)) + d = max(exc(2,1,ispin), exc(2,2,ispin)) + if (c>a .and. cb) then + nperm = nperm + 1 + endif + exit + end select + + enddo + phase = phase_dble(iand(nperm,1)) + +end + +subroutine get_mono_excitation(det1,det2,exc,phase,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Returns the excitation operator between two singly excited determinants and the phase + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det1(Nint,2) + integer(bit_kind), intent(in) :: det2(Nint,2) + integer, intent(out) :: exc(0:2,2,2) + double precision, intent(out) :: phase + integer :: tz + integer :: l, ispin, idx_hole, idx_particle, ishift + integer :: nperm + integer :: i,j,k,m,n + integer :: high, low + integer :: a,b,c,d + integer(bit_kind) :: hole, particle, tmp + double precision, parameter :: phase_dble(0:1) = (/ 1.d0, -1.d0 /) + + ASSERT (Nint > 0) + nperm = 0 + exc(0,1,1) = 0 + exc(0,2,1) = 0 + exc(0,1,2) = 0 + exc(0,2,2) = 0 + do ispin = 1,2 + ishift = 1-bit_kind_size + do l=1,Nint + ishift = ishift + bit_kind_size + if (det1(l,ispin) == det2(l,ispin)) then + cycle + endif + tmp = xor( det1(l,ispin), det2(l,ispin) ) + particle = iand(tmp, det2(l,ispin)) + hole = iand(tmp, det1(l,ispin)) + if (particle /= 0_bit_kind) then + tz = trailz(particle) + exc(0,2,ispin) = 1 + exc(1,2,ispin) = tz+ishift + endif + if (hole /= 0_bit_kind) then + tz = trailz(hole) + exc(0,1,ispin) = 1 + exc(1,1,ispin) = tz+ishift + endif + + if ( iand(exc(0,1,ispin),exc(0,2,ispin)) /= 1) then ! exc(0,1,ispin)/=1 and exc(0,2,ispin) /= 1 + cycle + endif + + low = min(exc(1,1,ispin),exc(1,2,ispin)) + high = max(exc(1,1,ispin),exc(1,2,ispin)) + + ASSERT (low > 0) + j = ishft(low-1,-bit_kind_shift)+1 ! Find integer in array(Nint) + n = iand(low,bit_kind_size-1) ! mod(low,bit_kind_size) + ASSERT (high > 0) + k = ishft(high-1,-bit_kind_shift)+1 + m = iand(high,bit_kind_size-1) + if (j==k) then + nperm = popcnt(iand(det1(j,ispin), & + iand(ibset(0_bit_kind,m-1)-1_bit_kind,ibclr(-1_bit_kind,n)+1_bit_kind))) + else + nperm = nperm + popcnt(iand(det1(k,ispin),ibset(0_bit_kind,m-1)-1_bit_kind)) +& + popcnt(iand(det1(j,ispin),ibclr(-1_bit_kind,n)+1_bit_kind)) + do i=j+1,k-1 + nperm = nperm + popcnt(det1(i,ispin)) + end do + endif + phase = phase_dble(iand(nperm,1)) + return + + enddo + enddo +end + + + + + +subroutine i_H_j(key_i,key_j,Nint,hij) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase,phase_2 + integer :: n_occ_alpha, n_occ_beta + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + !DEC$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint) + call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint) + has_mipi = .False. + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hij = hij + mipi(occ(k,1)) - miip(occ(k,1)) + enddo + do k = 1, elec_beta_num + hij = hij + mipi(occ(k,2)) + enddo + + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hij = hij + mipi(occ(k,1)) + enddo + do k = 1, elec_beta_num + hij = hij + mipi(occ(k,2)) - miip(occ(k,2)) + enddo + + endif + hij = phase*(hij + mo_mono_elec_integral(m,p)) + + case (0) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + + + +subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble) + use bitmasks + implicit none + BEGIN_DOC + ! Returns where i and j are determinants + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2) + double precision, intent(out) :: hij,hmono,hdouble + + integer :: exc(0:2,2,2) + integer :: degree + double precision :: get_mo_bielec_integral + integer :: m,n,p,q + integer :: i,j,k + integer :: occ(Nint*bit_kind_size,2) + double precision :: diag_H_mat_elem, phase,phase_2 + integer :: n_occ_alpha, n_occ_beta + logical :: has_mipi(Nint*bit_kind_size) + double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size) + PROVIDE mo_bielec_integrals_in_map mo_integrals_map + + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (sum(popcnt(key_i(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_i(:,2))) == elec_beta_num) + ASSERT (sum(popcnt(key_j(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(key_j(:,2))) == elec_beta_num) + + hij = 0.d0 + hmono = 0.d0 + hdouble = 0.d0 + !DEC$ FORCEINLINE + call get_excitation_degree(key_i,key_j,degree,Nint) + select case (degree) + case (2) + call get_double_excitation(key_i,key_j,exc,phase,Nint) + if (exc(0,1,1) == 1) then + ! Mono alpha, mono beta + hij = phase*get_mo_bielec_integral( & + exc(1,1,1), & + exc(1,1,2), & + exc(1,2,1), & + exc(1,2,2) ,mo_integrals_map) + else if (exc(0,1,1) == 2) then + ! Double alpha + hij = phase*(get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(1,2,1), & + exc(2,2,1) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,1), & + exc(2,1,1), & + exc(2,2,1), & + exc(1,2,1) ,mo_integrals_map) ) + else if (exc(0,1,2) == 2) then + ! Double beta + hij = phase*(get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(1,2,2), & + exc(2,2,2) ,mo_integrals_map) - & + get_mo_bielec_integral( & + exc(1,1,2), & + exc(2,1,2), & + exc(2,2,2), & + exc(1,2,2) ,mo_integrals_map) ) + endif + case (1) + call get_mono_excitation(key_i,key_j,exc,phase,Nint) + call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint) + call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint) + has_mipi = .False. + if (exc(0,1,1) == 1) then + ! Mono alpha + m = exc(1,1,1) + p = exc(1,2,1) + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hdouble = hdouble + mipi(occ(k,1)) - miip(occ(k,1)) + enddo + do k = 1, elec_beta_num + hdouble = hdouble + mipi(occ(k,2)) + enddo + + else + ! Mono beta + m = exc(1,1,2) + p = exc(1,2,2) + do k = 1, elec_beta_num + i = occ(k,2) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + miip(i) = get_mo_bielec_integral(m,i,i,p,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + do k = 1, elec_alpha_num + i = occ(k,1) + if (.not.has_mipi(i)) then + mipi(i) = get_mo_bielec_integral(m,i,p,i,mo_integrals_map) + has_mipi(i) = .True. + endif + enddo + + do k = 1, elec_alpha_num + hdouble = hdouble + mipi(occ(k,1)) + enddo + do k = 1, elec_beta_num + hdouble = hdouble + mipi(occ(k,2)) - miip(occ(k,2)) + enddo + + endif + hmono = mo_mono_elec_integral(m,p) + hij = phase*(hdouble + hmono) + + case (0) + hij = diag_H_mat_elem(key_i,Nint) + end select +end + + + +subroutine i_H_psi(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array) + use bitmasks + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet) + BEGIN_DOC + ! for the various Nstates + END_DOC + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + 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) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + enddo +end + +subroutine i_H_psi_sec_ord(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_interaction,interactions) + use bitmasks + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + double precision, intent(out) :: interactions(Ndet) + integer,intent(out) :: idx_interaction(0:Ndet) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet),n_interact + BEGIN_DOC + ! for the various Nstates + END_DOC + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_i_H_psi0(keys,key,Nint,Ndet,idx) + n_interact = 0 + do ii=1,idx(0) + i = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + if(dabs(hij).ge.1.d-8)then + if(i.ne.1)then + n_interact += 1 + interactions(n_interact) = hij + idx_interaction(n_interact) = i + endif + endif + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + enddo + idx_interaction(0) = n_interact +end + + +subroutine i_H_psi_SC2(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat) + use bitmasks + BEGIN_DOC + ! for the various Nstate + ! + ! returns in addition + ! + ! the array of the index of the non connected determinants to key1 + ! + ! in order to know what double excitation can be repeated on key1 + ! + ! idx_repeat(0) is the number of determinants that can be used + ! + ! to repeat the excitations + END_DOC + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + integer , intent(out) :: idx_repeat(0:Ndet) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet) + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat) + do ii=1,idx(0) + i = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + enddo +end + + +subroutine i_H_psi_SC2_verbose(key,keys,coef,Nint,Ndet,Ndet_max,Nstate,i_H_psi_array,idx_repeat) + use bitmasks + BEGIN_DOC + ! for the various Nstate + ! + ! returns in addition + ! + ! the array of the index of the non connected determinants to key1 + ! + ! in order to know what double excitation can be repeated on key1 + ! + ! idx_repeat(0) is the number of determinants that can be used + ! + ! to repeat the excitations + END_DOC + implicit none + integer, intent(in) :: Nint, Ndet,Ndet_max,Nstate + integer(bit_kind), intent(in) :: keys(Nint,2,Ndet) + integer(bit_kind), intent(in) :: key(Nint,2) + double precision, intent(in) :: coef(Ndet_max,Nstate) + double precision, intent(out) :: i_H_psi_array(Nstate) + integer , intent(out) :: idx_repeat(0:Ndet) + + integer :: i, ii,j + double precision :: phase + integer :: exc(0:2,2,2) + double precision :: hij + integer :: idx(0:Ndet) + + ASSERT (Nint > 0) + ASSERT (N_int == Nint) + ASSERT (Nstate > 0) + ASSERT (Ndet > 0) + ASSERT (Ndet_max >= Ndet) + i_H_psi_array = 0.d0 + call filter_connected_i_H_psi0_SC2(keys,key,Nint,Ndet,idx,idx_repeat) + print*,'--------' + do ii=1,idx(0) + print*,'--' + i = idx(ii) + !DEC$ FORCEINLINE + call i_H_j(keys(1,1,i),key,Nint,hij) + if (i==1)then + print*,'i==1 !!' + endif + print*,coef(i,1) * hij,coef(i,1),hij + do j = 1, Nstate + i_H_psi_array(j) = i_H_psi_array(j) + coef(i,j)*hij + enddo + print*,i_H_psi_array(1) + enddo + print*,'------' +end + + + +subroutine get_excitation_degree_vector(key1,key2,degree,Nint,sze,idx) + use bitmasks + implicit none + BEGIN_DOC + ! Applies get_excitation_degree to an array of determinants + END_DOC + integer, intent(in) :: Nint, sze + integer(bit_kind), intent(in) :: key1(Nint,2,sze) + integer(bit_kind), intent(in) :: key2(Nint,2) + integer, intent(out) :: degree(sze) + integer, intent(out) :: idx(0:sze) + + integer :: i,l + + ASSERT (Nint > 0) + ASSERT (sze > 0) + + l=1 + if (Nint==1) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree(l) = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))),-1) + if (degree(l) < 3) then + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==2) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree(l) = ishft(popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))),-1) + if (degree(l) < 3) then + idx(l) = i + l = l+1 + endif + enddo + + else if (Nint==3) then + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree(l) = ishft( popcnt(xor( key1(1,1,i), key2(1,1))) + & + popcnt(xor( key1(1,2,i), key2(1,2))) + & + popcnt(xor( key1(2,1,i), key2(2,1))) + & + popcnt(xor( key1(2,2,i), key2(2,2))) + & + popcnt(xor( key1(3,1,i), key2(3,1))) + & + popcnt(xor( key1(3,2,i), key2(3,2))),-1) + if (degree(l) < 3) then + idx(l) = i + l = l+1 + endif + enddo + + else + + !DIR$ LOOP COUNT (1000) + do i=1,sze + degree(l) = 0 + !DEC$ LOOP COUNT MIN(4) + do l=1,Nint + degree(l) = degree(l)+ popcnt(xor( key1(l,1,i), key2(l,1))) +& + popcnt(xor( key1(l,2,i), key2(l,2))) + enddo + degree(l) = ishft(degree(l),-1) + if (degree(l) < 3) then + idx(l) = i + l = l+1 + endif + enddo + + endif + idx(0) = l-1 +end + + + + +double precision function diag_H_mat_elem(det_in,Nint) + implicit none + BEGIN_DOC + ! Computes + END_DOC + integer,intent(in) :: Nint + integer(bit_kind),intent(in) :: det_in(Nint,2) + + integer(bit_kind) :: hole(Nint,2) + integer(bit_kind) :: particle(Nint,2) + integer :: i, nexc(2), ispin + integer :: occ_particle(Nint*bit_kind_size,2) + integer :: occ_hole(Nint*bit_kind_size,2) + integer(bit_kind) :: det_tmp(Nint,2) + integer :: na, nb + + ASSERT (Nint > 0) + ASSERT (sum(popcnt(det_in(:,1))) == elec_alpha_num) + ASSERT (sum(popcnt(det_in(:,2))) == elec_beta_num) + + nexc(1) = 0 + nexc(2) = 0 + do i=1,Nint + hole(i,1) = xor(det_in(i,1),ref_bitmask(i,1)) + hole(i,2) = xor(det_in(i,2),ref_bitmask(i,2)) + particle(i,1) = iand(hole(i,1),det_in(i,1)) + particle(i,2) = iand(hole(i,2),det_in(i,2)) + hole(i,1) = iand(hole(i,1),ref_bitmask(i,1)) + hole(i,2) = iand(hole(i,2),ref_bitmask(i,2)) + nexc(1) += popcnt(hole(i,1)) + nexc(2) += popcnt(hole(i,2)) + enddo + + diag_H_mat_elem = ref_bitmask_energy + if (nexc(1)+nexc(2) == 0) then + return + endif + + !call debug_det(det_in,Nint) + integer :: tmp + call bitstring_to_list(particle(1,1), occ_particle(1,1), tmp, Nint) + ASSERT (tmp == nexc(1)) + call bitstring_to_list(particle(1,2), occ_particle(1,2), tmp, Nint) + ASSERT (tmp == nexc(2)) + call bitstring_to_list(hole(1,1), occ_hole(1,1), tmp, Nint) + ASSERT (tmp == nexc(1)) + call bitstring_to_list(hole(1,2), occ_hole(1,2), tmp, Nint) + ASSERT (tmp == nexc(2)) + + det_tmp = ref_bitmask + do ispin=1,2 + na = elec_num_tab(ispin) + nb = elec_num_tab(iand(ispin,1)+1) + do i=1,nexc(ispin) + !DIR$ FORCEINLINE + call ac_operator( occ_particle(i,ispin), ispin, det_tmp, diag_H_mat_elem, Nint,na,nb) + !DIR$ FORCEINLINE + call a_operator ( occ_hole (i,ispin), ispin, det_tmp, diag_H_mat_elem, Nint,na,nb) + enddo + enddo +end + +subroutine a_operator(iorb,ispin,key,hjj,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Needed for diag_H_mat_elem + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hjj + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k > 0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibclr(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + !DIR$ FORCEINLINE + call get_occ_from_key(key,occ,Nint) + na -= 1 + + hjj -= mo_mono_elec_integral(iorb,iorb) + + ! Same spin + do i=1,na + hjj -= mo_bielec_integral_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + hjj -= mo_bielec_integral_jj(occ(i,other_spin),iorb) + enddo + +end + + +subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb) + use bitmasks + implicit none + BEGIN_DOC + ! Needed for diag_H_mat_elem + END_DOC + integer, intent(in) :: iorb, ispin, Nint + integer, intent(inout) :: na, nb + integer(bit_kind), intent(inout) :: key(Nint,2) + double precision, intent(inout) :: hjj + + integer :: occ(Nint*bit_kind_size,2) + integer :: other_spin + integer :: k,l,i + + ASSERT (iorb > 0) + ASSERT (ispin > 0) + ASSERT (ispin < 3) + ASSERT (Nint > 0) + + integer :: tmp + !DIR$ FORCEINLINE + call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint) + ASSERT (tmp == elec_alpha_num) + !DIR$ FORCEINLINE + call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint) + ASSERT (tmp == elec_beta_num) + + k = ishft(iorb-1,-bit_kind_shift)+1 + ASSERT (k > 0) + l = iorb - ishft(k-1,bit_kind_shift)-1 + key(k,ispin) = ibset(key(k,ispin),l) + other_spin = iand(ispin,1)+1 + + hjj += mo_mono_elec_integral(iorb,iorb) + + ! Same spin + do i=1,na + hjj += mo_bielec_integral_jj_anti(occ(i,ispin),iorb) + enddo + + ! Opposite spin + do i=1,nb + hjj += mo_bielec_integral_jj(occ(i,other_spin),iorb) + enddo + na += 1 +end + +subroutine get_occ_from_key(key,occ,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Returns a list of occupation numbers from a bitstring + END_DOC + integer(bit_kind), intent(in) :: key(Nint,2) + integer , intent(in) :: Nint + integer , intent(out) :: occ(Nint*bit_kind_size,2) + integer :: tmp + + call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint) + call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint) + +end + +subroutine H_u_0(v_0,u_0,H_jj,n,keys_tmp,Nint) + use bitmasks + implicit none + BEGIN_DOC + ! Computes v_0 = H|u_0> + ! + ! n : number of determinants + ! + ! H_jj : array of + END_DOC + integer, intent(in) :: n,Nint + double precision, intent(out) :: v_0(n) + double precision, intent(in) :: u_0(n) + double precision, intent(in) :: H_jj(n) + integer(bit_kind),intent(in) :: keys_tmp(Nint,2,n) + integer, allocatable :: idx(:) + double precision :: hij + double precision, allocatable :: vt(:) + integer :: i,j,k,l, jj + integer :: i0, j0 + ASSERT (Nint > 0) + ASSERT (Nint == N_int) + ASSERT (n>0) + PROVIDE ref_bitmask_energy + integer, parameter :: block_size = 157 + !$OMP PARALLEL DEFAULT(NONE) & + !$OMP PRIVATE(i,hij,j,k,idx,jj,vt) & + !$OMP SHARED(n,H_jj,u_0,keys_tmp,Nint,v_0) + !$OMP DO SCHEDULE(static) + do i=1,n + v_0(i) = H_jj(i) * u_0(i) + enddo + !$OMP END DO + allocate(idx(0:n), vt(n)) + Vt = 0.d0 + !$OMP DO SCHEDULE(guided) + do i=1,n + idx(0) = i + call filter_connected_davidson(keys_tmp,keys_tmp(1,1,i),Nint,i-1,idx) + do jj=1,idx(0) + j = idx(jj) + if ( (dabs(u_0(j)) > 1.d-7).or.((dabs(u_0(i)) > 1.d-7)) ) then + call i_H_j(keys_tmp(1,1,j),keys_tmp(1,1,i),Nint,hij) + vt (i) = vt (i) + hij*u_0(j) + vt (j) = vt (j) + hij*u_0(i) + endif + enddo + enddo + !$OMP END DO + !$OMP CRITICAL + do i=1,n + v_0(i) = v_0(i) + vt(i) + enddo + !$OMP END CRITICAL + deallocate(idx,vt) + !$OMP END PARALLEL +end + + + +BEGIN_PROVIDER [ integer, N_con_int ] + implicit none + BEGIN_DOC + ! Number of integers to represent the connections between determinants + END_DOC + N_con_int = 1 + ishft(N_det-1,-11) +END_PROVIDER + +BEGIN_PROVIDER [ integer*8, det_connections, (N_con_int,N_det) ] + implicit none + BEGIN_DOC + ! Build connection proxy between determinants + END_DOC + integer :: i,j + integer :: degree + integer :: j_int, j_k, j_l + integer, allocatable :: idx(:) + integer :: thread_num + integer :: omp_get_thread_num + + PROVIDE progress_bar + call start_progress(N_det,'Det connections',0.d0) + + select case(N_int) + + case(1) + + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections, & + !$OMP progress_bar,progress_value)& + !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num) + + !$ thread_num = omp_get_thread_num() + allocate (idx(0:N_det)) + !$OMP DO SCHEDULE(guided) + do i=1,N_det + if (thread_num == 0) then + progress_bar(1) = i + progress_value = dble(i) + endif + do j_int=1,N_con_int + det_connections(j_int,i) = 0_8 + j_k = ishft(j_int-1,11) + do j_l = j_k,min(j_k+2047,N_det), 32 + do j = j_l+1,min(j_l+32,i) + degree = popcnt(xor( psi_det(1,1,i),psi_det(1,1,j))) + & + popcnt(xor( psi_det(1,2,i),psi_det(1,2,j))) + if (degree < 5) then + det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) ) + exit + endif + enddo + enddo + enddo + enddo + !$OMP ENDDO + deallocate(idx) + !$OMP END PARALLEL + + case(2) + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections,& + !$OMP progress_bar,progress_value)& + !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num) + !$ thread_num = omp_get_thread_num() + allocate (idx(0:N_det)) + !$OMP DO SCHEDULE(guided) + do i=1,N_det + if (thread_num == 0) then + progress_bar(1) = i + progress_value = dble(i) + endif + do j_int=1,N_con_int + det_connections(j_int,i) = 0_8 + j_k = ishft(j_int-1,11) + do j_l = j_k,min(j_k+2047,N_det), 32 + do j = j_l+1,min(j_l+32,i) + degree = popcnt(xor( psi_det(1,1,i),psi_det(1,1,j))) + & + popcnt(xor( psi_det(1,2,i),psi_det(1,2,j))) + & + popcnt(xor( psi_det(2,1,i),psi_det(2,1,j))) + & + popcnt(xor( psi_det(2,2,i),psi_det(2,2,j))) + if (degree < 5) then + det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) ) + exit + endif + enddo + enddo + enddo + enddo + !$OMP ENDDO + deallocate(idx) + !$OMP END PARALLEL + + case(3) + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections,& + !$OMP progress_bar,progress_value)& + !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num) + !$ thread_num = omp_get_thread_num() + allocate (idx(0:N_det)) + !$OMP DO SCHEDULE(guided) + do i=1,N_det + if (thread_num == 0) then + progress_bar(1) = i + progress_value = dble(i) + endif + do j_int=1,N_con_int + det_connections(j_int,i) = 0_8 + j_k = ishft(j_int-1,11) + do j_l = j_k,min(j_k+2047,N_det), 32 + do j = j_l+1,min(j_l+32,i) + degree = popcnt(xor( psi_det(1,1,i),psi_det(1,1,j))) + & + popcnt(xor( psi_det(1,2,i),psi_det(1,2,j))) + & + popcnt(xor( psi_det(2,1,i),psi_det(2,1,j))) + & + popcnt(xor( psi_det(2,2,i),psi_det(2,2,j))) + & + popcnt(xor( psi_det(3,1,i),psi_det(3,1,j))) + & + popcnt(xor( psi_det(3,2,i),psi_det(3,2,j))) + if (degree < 5) then + det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) ) + exit + endif + enddo + enddo + enddo + enddo + !$OMP ENDDO + deallocate(idx) + !$OMP END PARALLEL + + case default + + + !$OMP PARALLEL DEFAULT (NONE) & + !$OMP SHARED(N_det, N_con_int, psi_det,N_int, det_connections,& + !$OMP progress_bar,progress_value)& + !$OMP PRIVATE(i,j_int,j_k,j_l,j,degree,idx,thread_num) + !$ thread_num = omp_get_thread_num() + allocate (idx(0:N_det)) + !$OMP DO SCHEDULE(guided) + do i=1,N_det + if (thread_num == 0) then + progress_bar(1) = i + progress_value = dble(i) + endif + do j_int=1,N_con_int + det_connections(j_int,i) = 0_8 + j_k = ishft(j_int-1,11) + do j_l = j_k,min(j_k+2047,N_det), 32 + do j = j_l+1,min(j_l+32,i) + !DIR$ FORCEINLINE + call get_excitation_degree(psi_det(1,1,i),psi_det(1,1,j),degree,N_int) + if (degree < 3) then + det_connections(j_int,i) = ibset( det_connections(j_int,i), iand(63,ishft(j_l,-5)) ) + exit + endif + enddo + enddo + enddo + enddo + !$OMP ENDDO + deallocate(idx) + !$OMP END PARALLEL + + end select + call stop_progress + +END_PROVIDER + diff --git a/src/Determinants/spindeterminants.ezfio_config b/src/Determinants/spindeterminants.ezfio_config new file mode 100644 index 00000000..39ccb82b --- /dev/null +++ b/src/Determinants/spindeterminants.ezfio_config @@ -0,0 +1,17 @@ +spindeterminants + n_det_alpha integer + n_det_beta integer + n_det integer + n_int integer + bit_kind integer + n_states integer + psi_det_alpha integer*8 (spindeterminants_n_int*spindeterminants_bit_kind/8,spindeterminants_n_det_alpha) + psi_det_beta integer*8 (spindeterminants_n_int*spindeterminants_bit_kind/8,spindeterminants_n_det_beta) + psi_coef_matrix_rows integer (spindeterminants_n_det) + psi_coef_matrix_columns integer (spindeterminants_n_det) + psi_coef_matrix_values double precision (spindeterminants_n_det,spindeterminants_n_states) + n_svd_coefs integer + psi_svd_alpha double precision (spindeterminants_n_det_alpha,spindeterminants_n_svd_coefs,spindeterminants_n_states) + psi_svd_beta double precision (spindeterminants_n_det_beta,spindeterminants_n_svd_coefs,spindeterminants_n_states) + psi_svd_coefs double precision (spindeterminants_n_svd_coefs,spindeterminants_n_states) + diff --git a/src/Determinants/spindeterminants.irp.f b/src/Determinants/spindeterminants.irp.f new file mode 100644 index 00000000..ffd28f85 --- /dev/null +++ b/src/Determinants/spindeterminants.irp.f @@ -0,0 +1,615 @@ +!==============================================================================! +! ! +! Independent alpha/beta parts ! +! ! +!==============================================================================! + +use bitmasks + +integer*8 function spin_det_search_key(det,Nint) + use bitmasks + implicit none + BEGIN_DOC +! Return an integer*8 corresponding to a determinant index for searching + END_DOC + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: det(Nint) + integer :: i + spin_det_search_key = det(1) + do i=2,Nint + spin_det_search_key = ieor(spin_det_search_key,det(i)) + enddo +end + + +BEGIN_PROVIDER [ integer(bit_kind), psi_det_alpha, (N_int,psi_det_size) ] + implicit none + BEGIN_DOC +! List of alpha determinants of psi_det + END_DOC + integer :: i,k + + do i=1,N_det + do k=1,N_int + psi_det_alpha(k,i) = psi_det(k,1,i) + enddo + enddo +END_PROVIDER + +BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta, (N_int,psi_det_size) ] + implicit none + BEGIN_DOC +! List of beta determinants of psi_det + END_DOC + integer :: i,k + + do i=1,N_det + do k=1,N_int + psi_det_beta(k,i) = psi_det(k,2,i) + enddo + enddo +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_alpha_unique, (N_int,psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_alpha_unique ] + implicit none + BEGIN_DOC + ! Unique alpha determinants + END_DOC + + integer :: i,k + integer, allocatable :: iorder(:) + integer*8, allocatable :: bit_tmp(:) + integer*8 :: last_key + integer*8, external :: spin_det_search_key + + allocate ( iorder(N_det), bit_tmp(N_det)) + + do i=1,N_det + iorder(i) = i + bit_tmp(i) = spin_det_search_key(psi_det_alpha(1,i),N_int) + enddo + + call i8sort(bit_tmp,iorder,N_det) + + N_det_alpha_unique = 0 + last_key = 0_8 + do i=1,N_det + if (bit_tmp(i) /= last_key) then + last_key = bit_tmp(i) + N_det_alpha_unique += 1 + do k=1,N_int + psi_det_alpha_unique(k,N_det_alpha_unique) = psi_det_alpha(k,iorder(i)) + enddo + endif + enddo + + deallocate (iorder, bit_tmp) +END_PROVIDER + + BEGIN_PROVIDER [ integer(bit_kind), psi_det_beta_unique, (N_int,psi_det_size) ] +&BEGIN_PROVIDER [ integer, N_det_beta_unique ] + implicit none + BEGIN_DOC + ! Unique beta determinants + END_DOC + + integer :: i,k + integer, allocatable :: iorder(:) + integer*8, allocatable :: bit_tmp(:) + integer*8 :: last_key + integer*8, external :: spin_det_search_key + + allocate ( iorder(N_det), bit_tmp(N_det)) + + do i=1,N_det + iorder(i) = i + bit_tmp(i) = spin_det_search_key(psi_det_beta(1,i),N_int) + enddo + + call i8sort(bit_tmp,iorder,N_det) + + N_det_beta_unique = 0 + last_key = 0_8 + do i=1,N_det + if (bit_tmp(i) /= last_key) then + last_key = bit_tmp(i) + N_det_beta_unique += 1 + do k=1,N_int + psi_det_beta_unique(k,N_det_beta_unique) = psi_det_beta(k,iorder(i)) + enddo + endif + enddo + + deallocate (iorder, bit_tmp) +END_PROVIDER + + + + + +integer function get_index_in_psi_det_alpha_unique(key,Nint) + use bitmasks + BEGIN_DOC +! Returns the index of the determinant in the ``psi_det_alpha_unique`` array + END_DOC + implicit none + + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key(Nint) + + integer :: i, ibegin, iend, istep, l + integer*8 :: det_ref, det_search + integer*8, external :: spin_det_search_key + logical :: is_in_wavefunction + + is_in_wavefunction = .False. + get_index_in_psi_det_alpha_unique = 0 + ibegin = 1 + iend = N_det_alpha_unique + 1 + + !DIR$ FORCEINLINE + det_ref = spin_det_search_key(key,Nint) + !DIR$ FORCEINLINE + det_search = spin_det_search_key(psi_det_alpha_unique(1,1),Nint) + + istep = ishft(iend-ibegin,-1) + i=ibegin+istep + do while (istep > 0) + !DIR$ FORCEINLINE + det_search = spin_det_search_key(psi_det_alpha_unique(1,i),Nint) + if ( det_search > det_ref ) then + iend = i + else if ( det_search == det_ref ) then + exit + else + ibegin = i + endif + istep = ishft(iend-ibegin,-1) + i = ibegin + istep + end do + + !DIR$ FORCEINLINE + do while (spin_det_search_key(psi_det_alpha_unique(1,i),Nint) == det_ref) + i = i-1 + if (i == 0) then + exit + endif + enddo + i += 1 + + if (i > N_det_alpha_unique) then + return + endif + + !DIR$ FORCEINLINE + do while (spin_det_search_key(psi_det_alpha_unique(1,i),Nint) == det_ref) + if (key(1) /= psi_det_alpha_unique(1,i)) then + continue + else + is_in_wavefunction = .True. + !DIR$ IVDEP + !DIR$ LOOP COUNT MIN(3) + do l=2,Nint + if (key(l) /= psi_det_alpha_unique(l,i)) then + is_in_wavefunction = .False. + endif + enddo + if (is_in_wavefunction) then + get_index_in_psi_det_alpha_unique = i + return + endif + endif + i += 1 + if (i > N_det_alpha_unique) then + return + endif + + enddo + +end + +integer function get_index_in_psi_det_beta_unique(key,Nint) + use bitmasks + BEGIN_DOC +! Returns the index of the determinant in the ``psi_det_beta_unique`` array + END_DOC + implicit none + + integer, intent(in) :: Nint + integer(bit_kind), intent(in) :: key(Nint) + + integer :: i, ibegin, iend, istep, l + integer*8 :: det_ref, det_search + integer*8, external :: spin_det_search_key + logical :: is_in_wavefunction + + is_in_wavefunction = .False. + get_index_in_psi_det_beta_unique = 0 + ibegin = 1 + iend = N_det_beta_unique + 1 + + !DIR$ FORCEINLINE + det_ref = spin_det_search_key(key,Nint) + !DIR$ FORCEINLINE + det_search = spin_det_search_key(psi_det_beta_unique(1,1),Nint) + + istep = ishft(iend-ibegin,-1) + i=ibegin+istep + do while (istep > 0) + !DIR$ FORCEINLINE + det_search = spin_det_search_key(psi_det_beta_unique(1,i),Nint) + if ( det_search > det_ref ) then + iend = i + else if ( det_search == det_ref ) then + exit + else + ibegin = i + endif + istep = ishft(iend-ibegin,-1) + i = ibegin + istep + end do + + !DIR$ FORCEINLINE + do while (spin_det_search_key(psi_det_beta_unique(1,i),Nint) == det_ref) + i = i-1 + if (i == 0) then + exit + endif + enddo + i += 1 + + if (i > N_det_beta_unique) then + return + endif + + !DIR$ FORCEINLINE + do while (spin_det_search_key(psi_det_beta_unique(1,i),Nint) == det_ref) + if (key(1) /= psi_det_beta_unique(1,i)) then + continue + else + is_in_wavefunction = .True. + !DIR$ IVDEP + !DIR$ LOOP COUNT MIN(3) + do l=2,Nint + if (key(l) /= psi_det_beta_unique(l,i)) then + is_in_wavefunction = .False. + endif + enddo + if (is_in_wavefunction) then + get_index_in_psi_det_beta_unique = i + return + endif + endif + i += 1 + if (i > N_det_beta_unique) then + return + endif + + enddo + +end + + +subroutine write_spindeterminants + use bitmasks + implicit none + integer*8, allocatable :: tmpdet(:,:) + integer :: N_int2 + integer :: i,j,k + integer*8 :: det_8(100) + integer(bit_kind) :: det_bk((100*8)/bit_kind) + equivalence (det_8, det_bk) + + N_int2 = (N_int*bit_kind)/8 + call ezfio_set_spindeterminants_n_det_alpha(N_det_alpha_unique) + call ezfio_set_spindeterminants_n_det_beta(N_det_beta_unique) + call ezfio_set_spindeterminants_n_det(N_det) + call ezfio_set_spindeterminants_n_int(N_int) + call ezfio_set_spindeterminants_bit_kind(bit_kind) + call ezfio_set_spindeterminants_n_states(N_states) + + allocate(tmpdet(N_int2,N_det_alpha_unique)) + do i=1,N_det_alpha_unique + do k=1,N_int + det_bk(k) = psi_det_alpha_unique(k,i) + enddo + do k=1,N_int2 + tmpdet(k,i) = det_8(k) + enddo + enddo + call ezfio_set_spindeterminants_psi_det_alpha(psi_det_alpha_unique) + deallocate(tmpdet) + + allocate(tmpdet(N_int2,N_det_beta_unique)) + do i=1,N_det_beta_unique + do k=1,N_int + det_bk(k) = psi_det_beta_unique(k,i) + enddo + do k=1,N_int2 + tmpdet(k,i) = det_8(k) + enddo + enddo + call ezfio_set_spindeterminants_psi_det_beta(psi_det_beta_unique) + deallocate(tmpdet) + + call ezfio_set_spindeterminants_psi_coef_matrix_values(psi_svd_matrix_values) + call ezfio_set_spindeterminants_psi_coef_matrix_rows(psi_svd_matrix_rows) + call ezfio_set_spindeterminants_psi_coef_matrix_columns(psi_svd_matrix_columns) + + integer :: n_svd_coefs + double precision :: norm, f + f = 1.d0/dble(N_states) + norm = 1.d0 + do n_svd_coefs=1,N_det_alpha_unique + do k=1,N_states + norm -= psi_svd_coefs(n_svd_coefs,k)*psi_svd_coefs(n_svd_coefs,k) + enddo + if (norm < 1.d-4) then + exit + endif + enddo + n_svd_coefs -= 1 + call ezfio_set_spindeterminants_n_svd_coefs(n_svd_coefs) + + double precision, allocatable :: dtmp(:,:,:) + allocate(dtmp(N_det_alpha_unique,n_svd_coefs,N_states)) + do k=1,N_states + do j=1,n_svd_coefs + do i=1,N_det_alpha_unique + dtmp(i,j,k) = psi_svd_alpha(i,j,k) + enddo + enddo + enddo + call ezfio_set_spindeterminants_psi_svd_alpha(dtmp) + deallocate(dtmp) + + allocate(dtmp(N_det_beta_unique,n_svd_coefs,N_states)) + do k=1,N_states + do j=1,n_svd_coefs + do i=1,N_det_beta_unique + dtmp(i,j,k) = psi_svd_beta(i,j,k) + enddo + enddo + enddo + call ezfio_set_spindeterminants_psi_svd_beta(dtmp) + deallocate(dtmp) + + allocate(dtmp(n_svd_coefs,N_states,1)) + do k=1,N_states + do j=1,n_svd_coefs + dtmp(j,k,1) = psi_svd_coefs(j,k) + enddo + enddo + call ezfio_set_spindeterminants_psi_svd_coefs(dtmp) + deallocate(dtmp) + +end + + +!==============================================================================! +! ! +! Alpha x Beta Matrix ! +! ! +!==============================================================================! + +BEGIN_PROVIDER [ double precision, psi_svd_matrix_values, (N_det,N_states) ] +&BEGIN_PROVIDER [ integer, psi_svd_matrix_rows, (N_det) ] +&BEGIN_PROVIDER [ integer, psi_svd_matrix_columns, (N_det) ] + use bitmasks + implicit none + BEGIN_DOC +! Matrix of wf coefficients. Outer product of alpha and beta determinants + END_DOC + integer :: i,j,k, l + integer(bit_kind) :: tmp_det(N_int,2) + integer :: idx + integer, external :: get_index_in_psi_det_sorted_bit + logical, external :: is_in_wavefunction + + + PROVIDE psi_coef_sorted_bit + +! l=0 +! do j=1,N_det_beta_unique +! do k=1,N_int +! tmp_det(k,2) = psi_det_beta_unique(k,j) +! enddo +! do i=1,N_det_alpha_unique +! do k=1,N_int +! tmp_det(k,1) = psi_det_alpha_unique(k,i) +! enddo +! idx = get_index_in_psi_det_sorted_bit(tmp_det,N_int) +! if (idx > 0) then +! l += 1 +! psi_svd_matrix_rows(l) = i +! psi_svd_matrix_columns(l) = j +! do k=1,N_states +! psi_svd_matrix_values(l,k) = psi_coef_sorted_bit(idx,k) +! enddo +! endif +! enddo +! enddo +! ASSERT (l == N_det) + + integer, allocatable :: iorder(:), to_sort(:) + integer, external :: get_index_in_psi_det_alpha_unique + integer, external :: get_index_in_psi_det_beta_unique + allocate(iorder(N_det), to_sort(N_det)) + do k=1,N_det + i = get_index_in_psi_det_alpha_unique(psi_det(1,1,k),N_int) + j = get_index_in_psi_det_beta_unique (psi_det(1,2,k),N_int) + do l=1,N_states + psi_svd_matrix_values(k,l) = psi_coef(k,l) + enddo + psi_svd_matrix_rows(k) = i + psi_svd_matrix_columns(k) = j + to_sort(k) = N_det_alpha_unique * (j-1) + i + iorder(k) = k + enddo + call isort(to_sort, iorder, N_det) + call iset_order(psi_svd_matrix_rows,iorder,N_det) + call iset_order(psi_svd_matrix_columns,iorder,N_det) + call dset_order(psi_svd_matrix_values,iorder,N_det) + deallocate(iorder,to_sort) +END_PROVIDER + +BEGIN_PROVIDER [ double precision, psi_svd_matrix, (N_det_alpha_unique,N_det_beta_unique,N_states) ] + implicit none + BEGIN_DOC +! Matrix of wf coefficients. Outer product of alpha and beta determinants + END_DOC + integer :: i,j,k,istate + psi_svd_matrix = 0.d0 + do k=1,N_det + i = psi_svd_matrix_rows(k) + j = psi_svd_matrix_columns(k) + do istate=1,N_states + psi_svd_matrix(i,j,istate) = psi_svd_matrix_values(k,istate) + enddo + enddo +END_PROVIDER + +subroutine create_wf_of_psi_svd_matrix + use bitmasks + implicit none + BEGIN_DOC +! Matrix of wf coefficients. Outer product of alpha and beta determinants + END_DOC + integer :: i,j,k + integer(bit_kind) :: tmp_det(N_int,2) + integer :: idx + integer, external :: get_index_in_psi_det_sorted_bit + logical, external :: is_in_wavefunction + double precision :: norm(N_states) + + call generate_all_alpha_beta_det_products + norm = 0.d0 + do j=1,N_det_beta_unique + do k=1,N_int + tmp_det(k,2) = psi_det_beta_unique(k,j) + enddo + do i=1,N_det_alpha_unique + do k=1,N_int + tmp_det(k,1) = psi_det_alpha_unique(k,i) + enddo + idx = get_index_in_psi_det_sorted_bit(tmp_det,N_int) + if (idx > 0) then + do k=1,N_states + psi_coef_sorted_bit(idx,k) = psi_svd_matrix(i,j,k) + norm(k) += psi_svd_matrix(i,j,k) + enddo + endif + enddo + enddo + do k=1,N_states + norm(k) = 1.d0/dsqrt(norm(k)) + do i=1,N_det + psi_coef_sorted_bit(i,k) = psi_coef_sorted_bit(i,k)*norm(k) + enddo + enddo + psi_det = psi_det_sorted_bit + psi_coef = psi_coef_sorted_bit + TOUCH psi_det psi_coef + psi_det = psi_det_sorted + psi_coef = psi_coef_sorted + norm(1) = 0.d0 + do i=1,N_det + norm(1) += psi_average_norm_contrib_sorted(i) + if (norm(1) >= 0.999999d0) then + exit + endif + enddo + N_det = min(i,N_det) + SOFT_TOUCH psi_det psi_coef N_det + +end + +subroutine generate_all_alpha_beta_det_products + implicit none + BEGIN_DOC +! Create a wave function from all possible alpha x beta determinants + END_DOC + integer :: i,j,k,l + integer :: idx, iproc + integer, external :: get_index_in_psi_det_sorted_bit + integer(bit_kind), allocatable :: tmp_det(:,:,:) + logical, external :: is_in_wavefunction + integer, external :: omp_get_thread_num + + !$OMP PARALLEL DEFAULT(NONE) SHARED(psi_coef_sorted_bit,N_det_beta_unique,& + !$OMP N_det_alpha_unique, N_int, psi_det_alpha_unique, psi_det_beta_unique,& + !$OMP N_det) & + !$OMP PRIVATE(i,j,k,l,tmp_det,idx,iproc) + !$ iproc = omp_get_thread_num() + allocate (tmp_det(N_int,2,N_det_alpha_unique)) + !$OMP DO + do j=1,N_det_beta_unique + l = 1 + do i=1,N_det_alpha_unique + do k=1,N_int + tmp_det(k,1,l) = psi_det_alpha_unique(k,i) + tmp_det(k,2,l) = psi_det_beta_unique (k,j) + enddo + if (.not.is_in_wavefunction(tmp_det(1,1,l),N_int,N_det)) then + l = l+1 + endif + enddo + call fill_H_apply_buffer_no_selection(l-1, tmp_det, N_int, iproc) + enddo + !$OMP END DO NOWAIT + deallocate(tmp_det) + !$OMP END PARALLEL + deallocate (tmp_det) + call copy_H_apply_buffer_to_wf + SOFT_TOUCH psi_det psi_coef N_det +end + + BEGIN_PROVIDER [ double precision, psi_svd_alpha, (N_det_alpha_unique,N_det_alpha_unique,N_states) ] +&BEGIN_PROVIDER [ double precision, psi_svd_beta , (N_det_beta_unique,N_det_beta_unique,N_states) ] +&BEGIN_PROVIDER [ double precision, psi_svd_coefs, (N_det_beta_unique,N_states) ] + implicit none + BEGIN_DOC + ! SVD wave function + END_DOC + + integer :: lwork, info, istate + double precision, allocatable :: work(:), tmp(:,:), copy(:,:) + allocate (work(1),tmp(N_det_beta_unique,N_det_beta_unique), & + copy(size(psi_svd_matrix,1),size(psi_svd_matrix,2))) + + do istate = 1,N_states + copy(:,:) = psi_svd_matrix(:,:,istate) + lwork=-1 + call dgesvd('A','A', N_det_alpha_unique, N_det_beta_unique, & + copy, size(copy,1), & + psi_svd_coefs(1,istate), psi_svd_alpha(1,1,istate), & + size(psi_svd_alpha,1), & + tmp, size(psi_svd_beta,2), & + work, lwork, info) + lwork = work(1) + deallocate(work) + allocate(work(lwork)) + call dgesvd('A','A', N_det_alpha_unique, N_det_beta_unique, & + copy, size(copy,1), & + psi_svd_coefs(1,istate), psi_svd_alpha(1,1,istate), & + size(psi_svd_alpha,1), & + tmp, size(psi_svd_beta,2), & + work, lwork, info) + deallocate(work) + if (info /= 0) then + print *, irp_here//': error in det SVD' + stop 1 + endif + integer :: i,j + do j=1,N_det_beta_unique + do i=1,N_det_beta_unique + psi_svd_beta(i,j,istate) = tmp(j,i) + enddo + enddo + deallocate(tmp,copy) + enddo + +END_PROVIDER + + diff --git a/src/Determinants/truncate_wf.irp.f b/src/Determinants/truncate_wf.irp.f new file mode 100644 index 00000000..f867ad7e --- /dev/null +++ b/src/Determinants/truncate_wf.irp.f @@ -0,0 +1,18 @@ +program cisd + implicit none + integer :: i,k + + + double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:) + integer :: N_st, degree + N_det=10000 + do i=1,N_det + do k=1,N_int + psi_det(k,1,i) = psi_det_sorted(k,1,i) + psi_det(k,2,i) = psi_det_sorted(k,2,i) + enddo + psi_coef(k,:) = psi_coef_sorted(k,:) + enddo + TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det + call save_wavefunction +end diff --git a/src/Determinants/utils.irp.f b/src/Determinants/utils.irp.f new file mode 100644 index 00000000..22faee83 --- /dev/null +++ b/src/Determinants/utils.irp.f @@ -0,0 +1,20 @@ +BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ] + implicit none + BEGIN_DOC + ! H matrix on the basis of the slater determinants defined by psi_det + END_DOC + integer :: i,j + double precision :: hij + call i_H_j(psi_det(1,1,1),psi_det(1,1,1),N_int,hij) + !$OMP PARALLEL DO SCHEDULE(GUIDED) PRIVATE(i,j,hij) & + !$OMP SHARED (N_det, psi_det, N_int,H_matrix_all_dets) + do i =1,N_det + do j =i,N_det + call i_H_j(psi_det(1,1,i),psi_det(1,1,j),N_int,hij) + H_matrix_all_dets(i,j) = hij + H_matrix_all_dets(j,i) = hij + enddo + enddo + !$OMP END PARALLEL DO +END_PROVIDER + diff --git a/src/Output/README.rst b/src/Output/README.rst index adcae302..7b510fc1 100644 --- a/src/Output/README.rst +++ b/src/Output/README.rst @@ -32,6 +32,7 @@ Needed Modules .. NEEDED_MODULES file. * `Utils `_ +* `Ezfio_files `_ Documentation ============= diff --git a/src/Properties/EZFIO.cfg b/src/Properties/EZFIO.cfg new file mode 100644 index 00000000..d230011d --- /dev/null +++ b/src/Properties/EZFIO.cfg @@ -0,0 +1,5 @@ +[z_one_point] +type: double precision +doc: z point on which the integrated delta rho is calculated +interface: input +default: 3.9 \ No newline at end of file