LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-12-22 20:35:19 +01:00
Code Issues Releases Wiki Activity

Add all the mising file...

Browse Source
This commit is contained in:
Thomas Applencourt 2015-04-20 16:45:06 +02:00
parent dcd4017cb9
commit bf997c5583
41 changed files with 7058 additions and 261 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

251
ocaml/Input_determinants.ml
View File

@ -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

1
scripts/ezfio_interface/ei_handler.py
View File

@ -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)

7
src/Determinants/ASSUMPTIONS.rst Normal file
View File

@ -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

103
src/Determinants/EZFIO.cfg Normal file
View File

@ -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

229
src/Determinants/H_apply.irp.f Normal file
View File

@ -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

542
src/Determinants/H_apply_template.f Normal file
View File

@ -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

6
src/Determinants/Makefile Normal file
View File

@ -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

1
src/Determinants/NEEDED_MODULES Normal file
View File

@ -0,0 +1 @@
AOs Bielec_integrals Bitmask Electrons Ezfio_files MonoInts MOs Nuclei Output Utils

696
src/Determinants/README.rst Normal file
View File

@ -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 <http://github.com/LCPQ/quantum_package/tree/master/src/AOs>`_
* `Bielec_integrals <http://github.com/LCPQ/quantum_package/tree/master/src/Bielec_integrals>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
* `MonoInts <http://github.com/LCPQ/quantum_package/tree/master/src/MonoInts>`_
* `MOs <http://github.com/LCPQ/quantum_package/tree/master/src/MOs>`_
* `Nuclei <http://github.com/LCPQ/quantum_package/tree/master/src/Nuclei>`_
* `Output <http://github.com/LCPQ/quantum_package/tree/master/src/Output>`_
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L100>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L187>`_
Fill the H_apply buffer with determiants for CISD
`h_apply_buffer_allocated <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L15>`_
Buffer of determinants/coefficients/perturbative energy for H_apply.
Uninitialized. Filled by H_apply subroutines.
`h_apply_buffer_lock <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L16>`_
Buffer of determinants/coefficients/perturbative energy for H_apply.
Uninitialized. Filled by H_apply subroutines.
`resize_h_apply_buffer <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/H_apply.irp.f#L48>`_
Undocumented
`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/SC2.irp.f#L1>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L155>`_
Undocumented
`connected_to_ref_by_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L253>`_
Undocumented
`det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L1>`_
Return an integer*8 corresponding to a determinant index for searching
`get_index_in_psi_det_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L48>`_
Returns the index of the determinant in the ``psi_det_sorted_bit`` array
`is_in_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L34>`_
True if the determinant ``det`` is in the wave function
`occ_pattern_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/connected_to_ref.irp.f#L17>`_
Return an integer*8 corresponding to a determinant index for searching
`do_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/create_excitations.irp.f#L1>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L382>`_
True if the Davidson algorithm is converged
`davidson_criterion <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L372>`_
Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
`davidson_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L18>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L68>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L1>`_
Max number of Davidson iterations
`davidson_sze_max <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L9>`_
Max number of Davidson sizes
`davidson_threshold <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/davidson.irp.f#L373>`_
Can be : [ energy | residual | both | wall_time | cpu_time | iterations ]
`one_body_dm_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L164>`_
One-body density matrix
`one_body_dm_mo_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L1>`_
Alpha and beta one-body density matrix for each state
`one_body_dm_mo_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L2>`_
Alpha and beta one-body density matrix for each state
`one_body_single_double_dm_mo_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L80>`_
Alpha and beta one-body density matrix for each state
`one_body_single_double_dm_mo_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L81>`_
Alpha and beta one-body density matrix for each state
`one_body_spin_density_mo <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L172>`_
rho(alpha) - rho(beta)
`save_natural_mos <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L196>`_
Save natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
`set_natural_mos <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L180>`_
Set natural orbitals, obtained by diagonalization of the one-body density matrix in the MO basis
`state_average_weight <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/density_matrix.irp.f#L207>`_
Weights in the state-average calculation of the density matrix
`det_svd <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/det_svd.irp.f#L1>`_
Computes the SVD of the Alpha x Beta determinant coefficient matrix
`filter_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L426>`_
Returns a determinant with only the 3 highest electrons
`int_of_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L391>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L32>`_
Maximum degree of excitation in the wf
`n_det <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L3>`_
Number of determinants in the wave function
`psi_average_norm_contrib <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L276>`_
Contribution of determinants to the state-averaged density
`psi_average_norm_contrib_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L306>`_
Wave function sorted by determinants contribution to the norm (state-averaged)
`psi_coef <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L230>`_
The wave function coefficients. Initialized with Hartree-Fock if the EZFIO file
is empty
`psi_coef_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L305>`_
Wave function sorted by determinants contribution to the norm (state-averaged)
`psi_coef_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L453>`_
Determinants on which we apply <i|H|j>.
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L336>`_
Determinants on which we apply <i|H|psi> 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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L65>`_
The wave function determinants. Initialized with Hartree-Fock if the EZFIO file
is empty
`psi_det_size <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L47>`_
Size of the psi_det/psi_coef arrays
`psi_det_sorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L304>`_
Wave function sorted by determinants contribution to the norm (state-averaged)
`psi_det_sorted_ab <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L452>`_
Determinants on which we apply <i|H|j>.
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L335>`_
Determinants on which we apply <i|H|psi> 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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L454>`_
Determinants on which we apply <i|H|j>.
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L583>`_
Reads the determinants from the EZFIO file
`save_wavefunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L630>`_
Save the wave function into the EZFIO file
`save_wavefunction_general <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L649>`_
Save the wave function into the EZFIO file
`save_wavefunction_unsorted <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L640>`_
Save the wave function into the EZFIO file
`sort_dets_by_3_highest_electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L474>`_
Determinants on which we apply <i|H|j>.
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants.irp.f#L349>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants_bitmasks.irp.f#L40>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants_bitmasks.irp.f#L31>`_
Number of double excitation bitmasks
`n_single_exc_bitmasks <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants_bitmasks.irp.f#L8>`_
Number of single excitation bitmasks
`single_exc_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/determinants_bitmasks.irp.f#L17>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L37>`_
Eigenvectors/values of the CI matrix
`ci_eigenvectors_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L38>`_
Eigenvectors/values of the CI matrix
`ci_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L36>`_
Eigenvectors/values of the CI matrix
`ci_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L18>`_
N_states lowest eigenvalues of the CI matrix
`diag_algorithm <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L1>`_
Diagonalization algorithm (Davidson or Lapack)
`diagonalize_ci <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI.irp.f#L96>`_
Replace the coefficients of the CI states by the coefficients of the
eigenstates of the CI matrix
`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_SC2.irp.f#L27>`_
Eigenvectors/values of the CI matrix
`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_SC2.irp.f#L26>`_
Eigenvectors/values of the CI matrix
`ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_SC2.irp.f#L1>`_
N_states_diag lowest eigenvalues of the CI matrix
`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_SC2.irp.f#L46>`_
Replace the coefficients of the CI states_diag by the coefficients of the
eigenstates of the CI matrix
`threshold_convergence_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_SC2.irp.f#L18>`_
convergence of the correlation energy of SC2 iterations
`ci_eigenvectors_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_mono.irp.f#L2>`_
Eigenvectors/values of the CI matrix
`ci_eigenvectors_s2_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_mono.irp.f#L3>`_
Eigenvectors/values of the CI matrix
`ci_electronic_energy_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_mono.irp.f#L1>`_
Eigenvectors/values of the CI matrix
`diagonalize_ci_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/diagonalize_CI_mono.irp.f#L59>`_
Replace the coefficients of the CI states by the coefficients of the
eigenstates of the CI matrix
`apply_mono <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/excitations_utils.irp.f#L1>`_
Undocumented
`filter_connected <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L2>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L163>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L293>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L392>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/filter_connected.irp.f#L101>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/guess_triplet.irp.f#L1>`_
Undocumented
`det_to_occ_pattern <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/occ_pattern.irp.f#L2>`_
Transform a determinant to an occupation pattern
`make_s2_eigenfunction <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/occ_pattern.irp.f#L251>`_
Undocumented
`n_occ_pattern <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/occ_pattern.irp.f#L143>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/occ_pattern.irp.f#L42>`_
Generate all possible determinants for a give occ_pattern
`occ_pattern_to_dets_size <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/occ_pattern.irp.f#L20>`_
Number of possible determinants for a given occ_pattern
`psi_occ_pattern <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/occ_pattern.irp.f#L142>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/occ_pattern.irp.f#L102>`_
Undocumented
`n_states_diag <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/options.irp.f#L1>`_
Number of states to consider for the diagonalization
`pouet <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/program_beginer_determinants.irp.f#L1>`_
Undocumented
`routine <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/program_beginer_determinants.irp.f#L7>`_
Undocumented
`idx_cas <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L5>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L62>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L6>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L63>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L3>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L4>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L47>`_
CAS determinants sorted to accelerate the search of a random determinant in the wave
function.
`psi_cas_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L46>`_
CAS determinants sorted to accelerate the search of a random determinant in the wave
function.
`psi_non_cas <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L60>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L61>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L100>`_
CAS determinants sorted to accelerate the search of a random determinant in the wave
function.
`psi_non_cas_sorted_bit <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/psi_cas.irp.f#L99>`_
CAS determinants sorted to accelerate the search of a random determinant in the wave
function.
`bi_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ref_bitmask.irp.f#L5>`_
Energy of the reference bitmask used in Slater rules
`kinetic_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ref_bitmask.irp.f#L3>`_
Energy of the reference bitmask used in Slater rules
`mono_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ref_bitmask.irp.f#L2>`_
Energy of the reference bitmask used in Slater rules
`nucl_elec_ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ref_bitmask.irp.f#L4>`_
Energy of the reference bitmask used in Slater rules
`ref_bitmask_energy <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/ref_bitmask.irp.f#L1>`_
Energy of the reference bitmask used in Slater rules
`expected_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L48>`_
Expected value of S2 : S*(S+1)
`get_s2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L1>`_
Returns <S^2>
`get_s2_u0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L82>`_
Undocumented
`s2_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L67>`_
array of the averaged values of the S^2 operator on the various states
`s_z <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L36>`_
z component of the Spin
`s_z2_sz <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/s2.irp.f#L37>`_
z component of the Spin
`prog_save_casino <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_casino.irp.f#L266>`_
Undocumented
`save_casino <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_casino.irp.f#L1>`_
Undocumented
`save_dets_qmcchem <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L1>`_
Undocumented
`save_for_qmc <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_for_qmcchem.irp.f#L48>`_
Undocumented
`save_natorb <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/save_natorb.irp.f#L1>`_
Undocumented
`a_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L962>`_
Needed for diag_H_mat_elem
`ac_operator <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1007>`_
Needed for diag_H_mat_elem
`decode_exc <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L76>`_
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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1139>`_
Build connection proxy between determinants
`diag_h_mat_elem <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L900>`_
Computes <i|H|i>
`get_double_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L141>`_
Returns the two excitation operators between two doubly excited determinants and the phase
`get_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L30>`_
Returns the excitation operators between two determinants and the phase
`get_excitation_degree <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1>`_
Returns the excitation degree between two determinants
`get_excitation_degree_vector <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L816>`_
Applies get_excitation_degree to an array of determinants
`get_mono_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L274>`_
Returns the excitation operator between two singly excited determinants and the phase
`get_occ_from_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1055>`_
Returns a list of occupation numbers from a bitstring
`h_u_0 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1071>`_
Computes v_0 = H|u_0>
.br
n : number of determinants
.br
H_jj : array of <j|H|j>
`i_h_j <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L355>`_
Returns <i|H|j> where i and j are determinants
`i_h_j_verbose <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L492>`_
Returns <i|H|j> where i and j are determinants
`i_h_psi <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L631>`_
<key|H|psi> for the various Nstates
`i_h_psi_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L713>`_
<key|H|psi> 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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L760>`_
<key|H|psi> 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 <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L666>`_
<key|H|psi> for the various Nstates
`n_con_int <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/slater_rules.irp.f#L1131>`_
Number of integers to represent the connections between determinants
`create_wf_of_psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L473>`_
Matrix of wf coefficients. Outer product of alpha and beta determinants
`generate_all_alpha_beta_det_products <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L528>`_
Create a wave function from all possible alpha x beta determinants
`get_index_in_psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L131>`_
Returns the index of the determinant in the ``psi_det_alpha_unique`` array
`get_index_in_psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L212>`_
Returns the index of the determinant in the ``psi_det_beta_unique`` array
`n_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L54>`_
Unique alpha determinants
`n_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L91>`_
Unique beta determinants
`psi_det_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L25>`_
List of alpha determinants of psi_det
`psi_det_alpha_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L53>`_
Unique alpha determinants
`psi_det_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L39>`_
List of beta determinants of psi_det
`psi_det_beta_unique <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L90>`_
Unique beta determinants
`psi_svd_alpha <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L568>`_
SVD wave function
`psi_svd_beta <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L569>`_
SVD wave function
`psi_svd_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L570>`_
SVD wave function
`psi_svd_matrix <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L457>`_
Matrix of wf coefficients. Outer product of alpha and beta determinants
`psi_svd_matrix_columns <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L398>`_
Matrix of wf coefficients. Outer product of alpha and beta determinants
`psi_svd_matrix_rows <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L397>`_
Matrix of wf coefficients. Outer product of alpha and beta determinants
`psi_svd_matrix_values <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L396>`_
Matrix of wf coefficients. Outer product of alpha and beta determinants
`spin_det_search_key <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L9>`_
Return an integer*8 corresponding to a determinant index for searching
`write_spindeterminants <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/spindeterminants.irp.f#L294>`_
Undocumented
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/truncate_wf.irp.f#L1>`_
Undocumented
`h_matrix_all_dets <http://github.com/LCPQ/quantum_package/tree/master/src/Determinants/utils.irp.f#L1>`_
H matrix on the basis of the slater determinants defined by psi_det

215
src/Determinants/SC2.irp.f Normal file
View File

@ -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

357
src/Determinants/connected_to_ref.irp.f Normal file
View File

@ -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

36
src/Determinants/create_excitations.irp.f Normal file
View File

@ -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

418
src/Determinants/davidson.irp.f Normal file
View File

@ -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 = <u_k | W_l> = <u_k| H |u_l>
! -------------------------------------------
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

214
src/Determinants/density_matrix.irp.f Normal file
View File

@ -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

61
src/Determinants/det_svd.irp.f Normal file
View File

@ -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

9
src/Determinants/determinants.irp.f
View File

@ -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

57
src/Determinants/determinants_bitmasks.irp.f Normal file
View File

@ -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

109
src/Determinants/diagonalize_CI.irp.f Normal file
View File

@ -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

59
src/Determinants/diagonalize_CI_SC2.irp.f Normal file
View File

@ -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

72
src/Determinants/diagonalize_CI_mono.irp.f Normal file
View File

@ -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

16
src/Determinants/excitations_utils.irp.f Normal file
View File

@ -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

611
src/Determinants/filter_connected.irp.f Normal file
View File

@ -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

79
src/Determinants/guess_doublet.irp.f Normal file
View File

@ -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

44
src/Determinants/guess_singlet.irp.f Normal file
View File

@ -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

48
src/Determinants/guess_triplet.irp.f Normal file
View File

@ -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

339
src/Determinants/occ_pattern.irp.f Normal file
View File

@ -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

22
src/Determinants/options.irp.f Normal file
View File

@ -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

138
src/Determinants/program_beginer_determinants.irp.f Normal file
View File

@ -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*,'<D_i|H|D_j> = ',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*,'<D_i|h(1)|D_j> = ',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*,'<ij|kl> = ',hij
end

114
src/Determinants/psi_cas.irp.f Normal file
View File

@ -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

57
src/Determinants/ref_bitmask.irp.f Normal file
View File

@ -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

106
src/Determinants/s2.irp.f Normal file
View File

@ -0,0 +1,106 @@
subroutine get_s2(key_i,key_j,phase,Nint)
implicit none
use bitmasks
BEGIN_DOC
! Returns <S^2>
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

268
src/Determinants/save_for_casino.irp.f Normal file
View File

@ -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

51
src/Determinants/save_for_qmcchem.irp.f Normal file
View File

@ -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

6
src/Determinants/save_natorb.irp.f Normal file
View File

@ -0,0 +1,6 @@
program save_natorb
read_wf = .True.
touch read_wf
call save_natural_mos
end

1301
src/Determinants/slater_rules.irp.f Normal file
View File

File diff suppressed because it is too large Load Diff

17
src/Determinants/spindeterminants.ezfio_config Normal file
View File

@ -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)

615
src/Determinants/spindeterminants.irp.f Normal file
View File

@ -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

18
src/Determinants/truncate_wf.irp.f Normal file
View File

@ -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

20
src/Determinants/utils.irp.f Normal file
View File

@ -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

1
src/Output/README.rst
View File

@ -32,6 +32,7 @@ Needed Modules
.. NEEDED_MODULES file.
* `Utils <http://github.com/LCPQ/quantum_package/tree/master/src/Utils>`_
* `Ezfio_files <http://github.com/LCPQ/quantum_package/tree/master/src/Ezfio_files>`_
Documentation
=============

5
src/Properties/EZFIO.cfg Normal file
View File

@ -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