10
0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-03 09:55:59 +02:00

Corrected bitmask bug for Multi-reference

This commit is contained in:
Anthony Scemama 2014-06-03 19:14:12 +02:00
parent 66963eefb8
commit c4cfc0b3d6
15 changed files with 203 additions and 82 deletions

View File

@ -26,7 +26,6 @@ double precision function ao_bielec_integral(i,j,k,l)
ao_bielec_integral = 0.d0
double precision :: thresh
thresh = ao_integrals_threshold
thresh = 0.d0
if (num_i /= num_j .or. num_k /= num_l .or. num_j /= num_k)then
do p = 1, 3

View File

@ -88,7 +88,6 @@ subroutine add_integrals_to_map(mask_ijkl)
call wall_time(wall_1)
call cpu_time(cpu_1)
mo_integrals_threshold = 0.d0
!$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, &
!$OMP bielec_tmp_0_idx, bielec_tmp_0, bielec_tmp_1,bielec_tmp_2,bielec_tmp_3,&

View File

@ -57,15 +57,40 @@ Documentation
`full_ijkl_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L12>`_
Bitmask to include all possible MOs
`generators_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L91>`_
Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator).
3rd index is :
* 1 : hole for single exc
* 1 : particle for single exc
* 3 : hole for 1st exc of double
* 4 : particle for 1st exc of double
* 5 : hole for 2dn exc of double
* 6 : particle for 2dn exc of double
`hf_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L32>`_
Hartree Fock bit mask
`i_bitmask_gen <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L182>`_
Current bitmask for the generators
`i_bitmask_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L190>`_
Current bitmask for the reference
`n_generators_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L58>`_
Number of bitmasks for generators
`n_int <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L3>`_
Number of 64-bit integers needed to represent determinants as binary strings
`n_reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L126>`_
Number of bitmasks for reference
`ref_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L50>`_
Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask
`reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L159>`_
Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference)
`bitstring_to_hexa <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L95>`_
Transform a bit string to a string in hexadecimal format for printing

View File

@ -107,12 +107,17 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_
if (exists) then
call ezfio_get_bitmasks_generators(generators_bitmask)
else
generators_bitmask(:,:,s_hole ,1) = HF_bitmask
generators_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
generators_bitmask(:,:,d_hole1,1) = HF_bitmask
generators_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
generators_bitmask(:,:,d_hole2,1) = HF_bitmask
generators_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
integer :: k, ispin
do k=1,N_generators_bitmask
do ispin=1,2
generators_bitmask(:,ispin,s_hole ,k) = full_ijkl_bitmask(:,d_hole1)
generators_bitmask(:,ispin,s_part ,k) = full_ijkl_bitmask(:,d_part1)
generators_bitmask(:,ispin,d_hole1,k) = full_ijkl_bitmask(:,d_hole1)
generators_bitmask(:,ispin,d_part1,k) = full_ijkl_bitmask(:,d_part1)
generators_bitmask(:,ispin,d_hole2,k) = full_ijkl_bitmask(:,d_hole2)
generators_bitmask(:,ispin,d_part2,k) = full_ijkl_bitmask(:,d_part2)
enddo
enddo
call ezfio_set_bitmasks_generators(generators_bitmask)
endif

View File

@ -2,10 +2,10 @@ module bitmasks
integer, parameter :: bit_kind_shift = 6 ! 5: 32 bits, 6: 64 bits
integer, parameter :: bit_kind_size = 64
integer, parameter :: bit_kind = 64/8
integer, parameter :: s_hole = 1
integer, parameter :: s_part = 2
integer, parameter :: d_hole1 = 3
integer, parameter :: d_part1 = 4
integer, parameter :: d_hole2 = 5
integer, parameter :: d_part2 = 6
integer, parameter :: d_hole1 = 1
integer, parameter :: d_part1 = 2
integer, parameter :: d_hole2 = 3
integer, parameter :: d_part2 = 4
integer, parameter :: s_hole = 5
integer, parameter :: s_part = 6
end module

View File

@ -23,7 +23,7 @@ Needed Modules
* `BiInts <http://github.com/LCPQ/quantum_package/tree/master/src/BiInts>`_
* `Bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask>`_
* `CISD <http://github.com/LCPQ/quantum_package/tree/master/src/CISD>`_
* `CISD_SC2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2>`_
* `SC2 <http://github.com/LCPQ/quantum_package/tree/master/src/SC2>`_
* `CISD_selected <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_selected>`_
* `Dets <http://github.com/LCPQ/quantum_package/tree/master/src/Dets>`_
* `Electrons <http://github.com/LCPQ/quantum_package/tree/master/src/Electrons>`_

View File

@ -19,7 +19,7 @@ BEGIN_PROVIDER [ logical, H_apply_buffer_allocated ]
! Uninitialized. Filled by H_apply subroutines.
END_DOC
integer :: iproc, sze
sze = 100
sze = 10000
if (.not.associated(H_apply_buffer)) then
allocate(H_apply_buffer(0:nproc-1))
iproc = 0

View File

@ -47,14 +47,6 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2, i_gene
occ_hole_tmp(N_int*bit_kind_size,2))
$init_thread
!print*,'key_in !!'
!call print_key(key_in)
!print*,'hole_1, particl_1'
!call print_key(hole_1)
!call print_key(particl_1)
!print*,'hole_2, particl_2'
!call print_key(hole_2)
!call print_key(particl_2)
!!!! First couple hole particle
@ -352,9 +344,11 @@ subroutine $subroutine($params_main)
$decls_main
integer :: i_generator, k, nmax
double precision :: wall_0, wall_1, wall_2, d
integer :: i_generator, nmax
double precision :: wall_0, wall_1, wall_2
integer(omp_lock_kind) :: lck
integer(bit_kind), allocatable :: mask(:,:,:)
integer :: ispin, k
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_selectors psi_generators
PROVIDE psi_det_sorted_bit coef_hf_selector
@ -363,22 +357,45 @@ subroutine $subroutine($params_main)
call wall_time(wall_1)
!$ call omp_init_lock(lck)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(i_generator,wall_2)
!$OMP PRIVATE(i_generator,wall_2,ispin,k,mask)
allocate( mask(N_int,2,6) )
!$OMP DO SCHEDULE(guided)
do i_generator=1,nmax
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_generators(k,ispin,i_generator) )
mask(k,ispin,s_part) = &
iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), &
not(psi_generators(k,ispin,i_generator)) )
mask(k,ispin,d_hole1) = &
iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), &
psi_generators(k,ispin,i_generator) )
mask(k,ispin,d_part1) = &
iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), &
not(psi_generators(k,ispin,i_generator)) )
mask(k,ispin,d_hole2) = &
iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), &
psi_generators(k,ispin,i_generator) )
mask(k,ispin,d_part2) = &
iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), &
not (psi_generators(k,ispin,i_generator)) )
enddo
enddo
call $subroutine_diexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,d_hole1,i_bitmask_gen), &
generators_bitmask(1,1,d_part1,i_bitmask_gen), &
generators_bitmask(1,1,d_hole2,i_bitmask_gen), &
generators_bitmask(1,1,d_part2,i_bitmask_gen), &
mask(1,1,d_hole1), mask(1,1,d_part1), &
mask(1,1,d_hole2), mask(1,1,d_part2), &
i_generator $params_post)
call $subroutine_monoexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,s_hole ,i_bitmask_gen), &
generators_bitmask(1,1,s_part ,i_bitmask_gen), &
mask(1,1,s_hole ), mask(1,1,s_part ), &
i_generator $params_post)
!$ call omp_set_lock(lck)
call wall_time(wall_2)
@ -390,23 +407,46 @@ subroutine $subroutine($params_main)
!$ call omp_unset_lock(lck)
enddo
!$OMP END DO
deallocate( mask )
!$OMP END PARALLEL
!$ call omp_destroy_lock(lck)
allocate( mask(N_int,2,6) )
do i_generator=nmax+1,N_det_generators
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_generators(k,ispin,i_generator) )
mask(k,ispin,s_part) = &
iand(generators_bitmask(k,ispin,s_part,i_bitmask_gen), &
not(psi_generators(k,ispin,i_generator)) )
mask(k,ispin,d_hole1) = &
iand(generators_bitmask(k,ispin,d_hole1,i_bitmask_gen), &
psi_generators(k,ispin,i_generator) )
mask(k,ispin,d_part1) = &
iand(generators_bitmask(k,ispin,d_part1,i_bitmask_gen), &
not(psi_generators(k,ispin,i_generator)) )
mask(k,ispin,d_hole2) = &
iand(generators_bitmask(k,ispin,d_hole2,i_bitmask_gen), &
psi_generators(k,ispin,i_generator) )
mask(k,ispin,d_part2) = &
iand(generators_bitmask(k,ispin,d_part2,i_bitmask_gen), &
not(psi_generators(k,ispin,i_generator)) )
enddo
enddo
call $subroutine_diexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,d_hole1,i_bitmask_gen), &
generators_bitmask(1,1,d_part1,i_bitmask_gen), &
generators_bitmask(1,1,d_hole2,i_bitmask_gen), &
generators_bitmask(1,1,d_part2,i_bitmask_gen), &
mask(1,1,d_hole1), mask(1,1,d_part1), &
mask(1,1,d_hole2), mask(1,1,d_part2), &
i_generator $params_post)
call $subroutine_monoexc(psi_generators(1,1,i_generator), &
generators_bitmask(1,1,s_hole ,i_bitmask_gen), &
generators_bitmask(1,1,s_part ,i_bitmask_gen), &
mask(1,1,s_hole ), mask(1,1,s_part ), &
i_generator $params_post)
call wall_time(wall_2)
$printout_always
@ -419,6 +459,7 @@ subroutine $subroutine($params_main)
$copy_buffer
$generate_psi_guess
abort_here = abort_all
deallocate( mask )
end

View File

@ -56,16 +56,40 @@ logical function is_in_wavefunction(key,Nint,Ndet)
endif
enddo
if (is_in_wavefunction) then
return
exit
endif
endif
i += 1
if (i > N_det) then
exit
return
! exit
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)
@ -92,14 +116,10 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
N_past = max(1,N_past_in)
if (Nint == 1) then
do i=N_past-1,1,-1
do i=1,N_past-1
degree_x2 = popcnt(xor( key(1,1), keys(1,1,i))) + &
popcnt(xor( key(1,2), keys(1,2,i)))
if(degree_x2 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
if (degree_x2 > 4) then
cycle
else
connected_to_ref = i
@ -112,16 +132,12 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
else if (Nint==2) then
do i=N_past-1,1,-1
do i=1,N_past-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 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
if (degree_x2 > 4) then
cycle
else
connected_to_ref = i
@ -133,18 +149,14 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
else if (Nint==3) then
do i=N_past-1,1,-1
do i=1,N_past-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 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
if (degree_x2 > 4) then
cycle
else
connected_to_ref = i
@ -156,7 +168,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
else
do i=N_past-1,1,-1
do i=1,N_past-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)
@ -164,11 +176,7 @@ integer function connected_to_ref(key,keys,Nint,N_past_in,Ndet)
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 == 0)then
connected_to_ref = -i
return
endif
if (degree_x2 > 5) then
if (degree_x2 > 4) then
cycle
else
connected_to_ref = i

View File

@ -6,5 +6,9 @@ s = H_apply("FCI",openmp=True)
s.set_selection_pt2("epstein_nesbet_2x2")
print s
s = H_apply("FCI_PT2",openmp=True)
s.set_perturbation("epstein_nesbet_2x2")
print s
END_SHELL

View File

@ -5,9 +5,9 @@ program cisd
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, degree
character*(64) :: perturbation
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
character*(64) :: perturbation
pt2 = 1.d0
diag_algorithm = "Lapack"

View File

@ -33,8 +33,8 @@ BEGIN_PROVIDER [ integer, N_det_generators ]
N_det_generators = N_det
do i=1,N_det
norm = norm + psi_average_norm_contrib_sorted(i)
if (norm > threshold_generators) then
N_det_generators = i-1
if (norm >= threshold_generators) then
N_det_generators = i
exit
endif
enddo

View File

@ -66,14 +66,24 @@ subroutine pt2_epstein_nesbet_2x2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det_selectors,psi_selectors_size,N_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)
do i =1,N_st
delta_e = h - CI_electronic_energy(i)
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
if (i_H_psi_array(i) /= 0.d0) then
delta_e = h - CI_electronic_energy(i)
if (delta_e > 0.d0) then
e_2_pert(i) = 0.5d0 * (delta_e - dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
else
e_2_pert(i) = 0.5d0 * (delta_e + dsqrt(delta_e * delta_e + 4.d0 * i_H_psi_array(i) * i_H_psi_array(i)))
endif
if (dabs(i_H_psi_array(i)) > 1.d-6) then
c_pert(i) = e_2_pert(i)/i_H_psi_array(i)
else
c_pert(i) = 0.d0
endif
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
else
c_pert(i) = -1.d0
e_2_pert(i) = 0.d0
c_pert(i) = 0.d0
H_pert_diag(i) = 0.d0
endif
H_pert_diag(i) = h*c_pert(i)*c_pert(i)
enddo
end

View File

@ -8,7 +8,7 @@ Documentation
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/SC2.irp.f#L1>`_
`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/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
@ -24,25 +24,51 @@ Documentation
.br
Initial guess vectors are not necessarily orthonormal
`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/SC2.irp.f#L169>`_
`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/SC2.irp.f#L169>`_
Undocumented
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/cisd_SC2.irp.f#L1>`_
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/cisd_SC2.irp.f#L1>`_
Undocumented
`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L19>`_
`ci_sc2_eigenvectors <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L19>`_
Eigenvectors/values of the CI matrix
`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L18>`_
`ci_sc2_electronic_energy <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L18>`_
Eigenvectors/values of the CI matrix
`ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L1>`_
`ci_sc2_energy <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L1>`_
N_states lowest eigenvalues of the CI matrix
`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD_SC2/diagonalize_CI_SC2.irp.f#L38>`_
`diagonalize_ci_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/diagonalize_CI_SC2.irp.f#L38>`_
Replace the coefficients of the CI states by the coefficients of the
eigenstates of the CI matrix
`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/pert_sc2.irp.f#L2>`_
compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
.br
for the various N_st states,
.br
but with the correction in the denominator
.br
comming from the interaction of that determinant with all the others determinants
.br
that can be repeated by repeating all the double excitations
.br
: you repeat all the correlation energy already taken into account in CI_electronic_energy(1)
.br
that could be repeated to this determinant.
.br
In addition, for the perturbative energetic contribution you have the standard second order
.br
e_2_pert = <psi_i|H|det_pert>^2/(Delta_E)
.br
and also the purely projected contribution
.br
H_pert_diag = <HF|H|det_pert> c_pert
`repeat_all_e_corr <http://github.com/LCPQ/quantum_package/tree/master/src/SC2/pert_sc2.irp.f#L82>`_
Undocumented
Needed Modules

View File

@ -51,7 +51,11 @@ END_PROVIDER
E_corr_per_selectors(i) = -1000.d0
endif
enddo
inv_selectors_coef_hf = 1.d0/coef_hf_selector
if (dabs(coef_hf_selector) > 1.d-8) then
inv_selectors_coef_hf = 1.d0/coef_hf_selector
else
inv_selectors_coef_hf = 0.d0
endif
do i = 1,n_double_selectors
E_corr_per_selectors(double_index_selectors(i)) *=inv_selectors_coef_hf
enddo