LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2025-01-08 20:33:26 +01:00
Code Issues Releases Wiki Activity

Added guess_lowest_state

Browse Source
This commit is contained in:
Anthony Scemama 2015-11-12 01:08:04 +01:00
parent 2e27f8bd37
commit 1e15823494
6 changed files with 190 additions and 147 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/Determinants/H_apply.irp.f
View File

@ -41,6 +41,15 @@ type(H_apply_buffer_type), pointer :: H_apply_buffer(:)
call omp_init_lock(H_apply_buffer_lock(1,iproc))
!$OMP END PARALLEL
endif
do iproc=2,nproc-1
if (.not.allocated(H_apply_buffer(iproc)%det)) then
print *, ' ===================== Error =================== '
print *, 'H_apply_buffer_allocated should be provided outside'
print *, 'of an OpenMP section'
print *, ' =============================================== '
stop
endif
enddo
END_PROVIDER
@ -111,7 +120,6 @@ subroutine copy_H_apply_buffer_to_wf
double precision, allocatable :: buffer_coef(:,:)
integer :: i,j,k
integer :: N_det_old
integer :: iproc
PROVIDE H_apply_buffer_allocated
@ -158,7 +166,7 @@ subroutine copy_H_apply_buffer_to_wf
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)
!$OMP SHARED(N_int,H_apply_buffer,psi_det,psi_coef,N_states,psi_det_size)
j=0
!$ j=omp_get_thread_num()
do k=0,j-1

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

@ -8,6 +8,7 @@ BEGIN_PROVIDER [ integer, N_det ]
logical :: exists
character*64 :: label
PROVIDE ezfio_filename
PROVIDE nproc
if (read_wf) then
call ezfio_has_determinants_n_det(exists)
if (exists) then

170
src/Determinants/guess_lowest_state.irp.f Normal file
View File

@ -0,0 +1,170 @@
program first_guess
use bitmasks
implicit none
BEGIN_DOC
! Select all the determinants with the lowest energy as a starting point.
END_DOC
integer :: i,j
double precision, allocatable :: orb_energy(:)
double precision :: E
integer, allocatable :: kept(:)
integer :: nelec_kept(2)
PROVIDE H_apply_buffer_allocated
allocate (orb_energy(mo_tot_num), kept(0:mo_tot_num))
nelec_kept(1:2) = 0
kept(0) = 0
do i=1,mo_tot_num
orb_energy(i) = mo_mono_elec_integral(i,i)
do j=1,elec_beta_num
if (i==j) cycle
orb_energy(i) += mo_bielec_integral_jj_anti(i,j)
enddo
do j=1,elec_alpha_num
orb_energy(i) += mo_bielec_integral_jj(i,j)
enddo
if ( (orb_energy(i) > -1.d0).and.(orb_energy(i) < .5d0) ) then
kept(0) += 1
kept( kept(0) ) = i
if (i <= elec_beta_num) then
nelec_kept(2) += 1
endif
if (i <= elec_alpha_num) then
nelec_kept(1) += 1
endif
endif
enddo
integer, allocatable :: list (:,:)
integer(bit_kind), allocatable :: string(:,:)
allocate ( list(N_int*bit_kind_size,2), string(N_int,2) )
string = ref_bitmask
call bitstring_to_list( string(1,1), list(1,1), elec_alpha_num, N_int)
call bitstring_to_list( string(1,2), list(1,2), elec_beta_num , N_int)
psi_det_alpha_unique(:,1) = string(:,1)
psi_det_beta_unique (:,1) = string(:,2)
N_det_alpha_unique = 1
N_det_beta_unique = 1
integer :: i1,i2,i3,i4,i5,i6,i7,i8,i9
! select case (nelec_kept(1))
!
! case(0)
! continue
!
! case(1)
! do i1=kept(0),1,-1
! list(elec_alpha_num-0,1) = kept(i1)
! call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
! N_det_alpha_unique += 1
! psi_det_alpha_unique(:,N_det_alpha_unique) = string(:,1)
! enddo
!
! case(2)
! do i1=kept(0),1,-1
! list(elec_alpha_num-0,1) = kept(i1)
! call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
! N_det_alpha_unique += 1
! psi_det_alpha_unique(:,N_det_alpha_unique) = string(:,1)
! do i2=i1-1,1,-1
! list(elec_alpha_num-1,1) = kept(i2)
! call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
! N_det_alpha_unique += 1
! psi_det_alpha_unique(:,N_det_alpha_unique) = string(:,1)
! enddo
! enddo
psi_det_size = kept(0)**(nelec_kept(1)+nelec_kept(2))
TOUCH psi_det_size
BEGIN_SHELL [ /usr/bin/python ]
template_alpha_ext = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_alpha_num-%(i)d,1) = kept(%(i2)s)
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
"""
template_alpha = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_alpha_num-%(i)d,1) = kept(%(i2)s)
call list_to_bitstring( string(1,1), list(1,1), elec_alpha_num, N_int)
N_det_alpha_unique += 1
psi_det_alpha_unique(:,N_det_alpha_unique) = string(:,1)
"""
template_beta_ext = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_beta_num-%(i)d,2) = kept(%(i2)s)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
"""
template_beta = """
do %(i2)s = %(i1)s-1,1,-1
list(elec_beta_num-%(i)d,2) = kept(%(i2)s)
call list_to_bitstring( string(1,2), list(1,2), elec_beta_num, N_int)
N_det_beta_unique += 1
psi_det_beta_unique(:,N_det_beta_unique) = string(:,2)
"""
def write(template_ext,template,imax):
print "case(%d)"%(imax)
def aux(i2,i1,i,j):
if (i==imax-1):
print template%locals()
else:
print template_ext%locals()
i += 1
j -= 1
if (i != imax):
i1 = "i%d"%(i)
i2 = "i%d"%(i+1)
aux(i2,i1,i,j)
print "enddo"
i2 = "i1"
i1 = "kept(0)+1"
i = 0
aux(i2,i1,i,imax)
def main():
print """
select case (nelec_kept(1))
case(0)
continue
"""
for imax in range(1,10):
write(template_alpha_ext,template_alpha,imax)
print """
end select
select case (nelec_kept(2))
case(0)
continue
"""
for imax in range(1,10):
write(template_beta_ext,template_beta,imax)
print "end select"
main()
END_SHELL
TOUCH N_det_alpha_unique N_det_beta_unique psi_det_alpha_unique psi_det_beta_unique
call create_wf_of_psi_bilinear_matrix(.False.)
call diagonalize_ci
j= N_det
do i=1,N_det
if (psi_average_norm_contrib_sorted(i) < 1.d-6) then
j = i-1
exit
endif
! call debug_det(psi_det_sorted(1,1,i),N_int)
enddo
call save_wavefunction_general(j,N_states,psi_det_sorted,size(psi_coef_sorted,1),psi_coef_sorted)
deallocate(orb_energy, kept, list, string)
end

138
src/Determinants/program_initial_determinants.irp.f
View File

@ -1,138 +0,0 @@
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

10
src/Determinants/spindeterminants.irp.f
View File

@ -442,13 +442,14 @@ BEGIN_PROVIDER [ double precision, psi_bilinear_matrix, (N_det_alpha_unique,N_de
enddo
END_PROVIDER
subroutine create_wf_of_psi_bilinear_matrix
subroutine create_wf_of_psi_bilinear_matrix(truncate)
use bitmasks
implicit none
BEGIN_DOC
! Generate a wave function containing all possible products
! of alpha and beta determinants
END_DOC
logical, intent(in) :: truncate
integer :: i,j,k
integer(bit_kind) :: tmp_det(N_int,2)
integer :: idx
@ -488,8 +489,10 @@ subroutine create_wf_of_psi_bilinear_matrix
norm(1) = 0.d0
do i=1,N_det
norm(1) += psi_average_norm_contrib_sorted(i)
if (norm(1) >= 0.999999d0) then
exit
if (truncate) then
if (norm(1) >= 0.999999d0) then
exit
endif
endif
enddo
N_det = min(i,N_det)
@ -532,7 +535,6 @@ subroutine generate_all_alpha_beta_det_products
!$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

6
src/Determinants/truncate_wf.irp.f
View File

@ -8,10 +8,10 @@ program cisd
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)
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,:)
psi_coef(i,:) = psi_coef_sorted(i,:)
enddo
TOUCH psi_det psi_coef psi_det_sorted psi_coef_sorted psi_average_norm_contrib_sorted N_det
call save_wavefunction