LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-22 05:02:15 +02:00
Code Issues Releases Wiki Activity

Merge pull request #117 from scemama/master

Browse Source 
Accelerate selection + PT2
This commit is contained in:
Anthony Scemama 2015-11-21 00:42:34 +01:00
commit 9e24f7d6c1
9 changed files with 187 additions and 71 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
configure vendored
View File

@ -33,8 +33,8 @@ import sys
from os.path import join
if not any(i in ["--production", "--development"] for i in sys.argv):
print __doc__
sys.exit()
sys.argv += ["--development"]
if len(sys.argv) != 3:
print __doc__
sys.exit()

4
plugins/MP2/H_apply.irp.f
View File

@ -6,5 +6,9 @@ from perturbation import perturbations
s = H_apply("mp2")
s.set_perturbation("Moller_plesset")
print s
s = H_apply("mp2_selection")
s.set_selection_pt2("Moller_plesset")
print s
END_SHELL

31
plugins/MP2/mp2_wf.irp.f Normal file
View File

@ -0,0 +1,31 @@
program mp2_wf
implicit none
BEGIN_DOC
! Save the MP2 wave function
END_DOC
integer :: i,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer :: N_st, iter
N_st = N_states
allocate (pt2(N_st), norm_pert(N_st), H_pert_diag(N_st))
pt2 = 1.d0
selection_criterion = 1.e-12
selection_criterion_min = 1.e-12
TOUCH selection_criterion_min selection_criterion selection_criterion_factor
call H_apply_mp2_selection(pt2, norm_pert, H_pert_diag, N_st)
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
touch N_det psi_det psi_coef
print*,'N_det = ',N_det
print*,'-----'
print *, 'PT2 = ', pt2(1)
print *, 'E = ', HF_energy
print *, 'E_before +PT2 = ', HF_energy+pt2(1)
N_det = min(N_det,N_det_max)
call save_wavefunction
call ezfio_set_mp2_energy(HF_energy+pt2(1))
deallocate(pt2,norm_pert,H_pert_diag)
end

16
plugins/Perturbation/Moller_plesset.irp.f
View File

@ -26,10 +26,18 @@ subroutine pt2_moller_plesset(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
ASSERT (Nint == N_int)
ASSERT (Nint > 0)
call get_excitation(ref_bitmask,det_pert,exc,degree,phase,Nint)
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
(Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
delta_e = 1.d0/delta_e
if (degree == 2) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) + Fock_matrix_diag_mo(h2) - &
(Fock_matrix_diag_mo(p1) + Fock_matrix_diag_mo(p2))
delta_e = 1.d0/delta_e
else if (degree == 1) then
call decode_exc(exc,degree,h1,p1,h2,p2,s1,s2)
delta_e = Fock_matrix_diag_mo(h1) - Fock_matrix_diag_mo(p1)
delta_e = 1.d0/delta_e
else
delta_e = 0.d0
endif
call i_H_psi(det_pert,psi_selectors,psi_selectors_coef,Nint,N_det,psi_selectors_size,n_st,i_H_psi_array)
h = diag_H_mat_elem(det_pert,Nint)

6
scripts/generate_h_apply.py
View File

@ -60,7 +60,7 @@ class H_apply(object):
s["omp_master"] = "!$OMP MASTER"
s["omp_end_master"] = "!$OMP END MASTER"
s["omp_barrier"] = "!$OMP BARRIER"
s["omp_do"] = "!$OMP DO SCHEDULE (static)"
s["omp_do"] = "!$OMP DO SCHEDULE (static,1)"
s["omp_enddo"] = "!$OMP ENDDO NOWAIT"
d = { True : '.True.', False : '.False.'}
@ -201,7 +201,7 @@ class H_apply(object):
"""
self.data["size_max"] = "256"
self.data["initialization"] = """
PROVIDE CI_electronic_energy psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
PROVIDE psi_selectors_coef psi_selectors E_corr_per_selectors psi_det_sorted_bit
"""
self.data["keys_work"] = """
call perturb_buffer_%s(i_generator,keys_out,key_idx,e_2_pert_buffer,coef_pert_buffer,sum_e_2_pert, &
@ -219,7 +219,7 @@ class H_apply(object):
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag(N_st)
double precision :: delta_pt2(N_st), norm_psi(N_st), pt2_old(N_st)
PROVIDE CI_electronic_energy N_det_generators
PROVIDE N_det_generators
do k=1,N_st
pt2(k) = 0.d0
norm_pert(k) = 0.d0

21
src/Determinants/H_apply.template.f
View File

@ -10,6 +10,9 @@ subroutine $subroutine_diexc(key_in, key_prev, hole_1,particl_1, hole_2, particl
integer :: highest, p1,p2,sp,ni,i,mi,nt,ns
integer(bit_kind), intent(in) :: key_prev(N_int, 2, *)
PROVIDE N_int
PROVIDE N_det
$declarations
@ -183,10 +186,8 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
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)
call bitstring_to_list_ab(particle,occ_particle,N_elec_in_key_part_1,N_int)
call bitstring_to_list_ab(hole,occ_hole,N_elec_in_key_hole_1,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))
@ -249,10 +250,8 @@ subroutine $subroutine_diexcOrg(key_in,key_mask,hole_1,particl_1,hole_2, particl
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)
call bitstring_to_list_ab(particle_tmp,occ_particle_tmp,N_elec_in_key_part_2,N_int)
call bitstring_to_list_ab(hole_tmp,occ_hole_tmp,N_elec_in_key_hole_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
@ -444,10 +443,8 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1,i_generator,iproc_in $pa
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)
call bitstring_to_list_ab(particle,occ_particle,N_elec_in_key_part_1,N_int)
call bitstring_to_list_ab(hole,occ_hole,N_elec_in_key_hole_1,N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2))
do ispin=1,2

2
src/Determinants/density_matrix.irp.f
View File

@ -71,7 +71,6 @@
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
@ -157,7 +156,6 @@ END_PROVIDER
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

29
src/Determinants/filter_connected.irp.f
View File

@ -132,7 +132,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
popcnt(xor( key1(1,2,i), key2(1,2)))
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
else if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
@ -148,7 +148,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
popcnt(xor( key1(2,2,i), key2(2,2)))
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
else if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
@ -166,7 +166,7 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
popcnt(xor( key1(3,2,i), key2(3,2)))
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
else if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
@ -175,23 +175,28 @@ subroutine filter_connected_i_H_psi0(key1,key2,Nint,sze,idx)
else
!DIR$ LOOP COUNT (1000)
do i=1,sze
outer: 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
if ( key1(m,1,i) /= key2(m,1)) then
degree_x2 = degree_x2+ popcnt(xor( key1(m,1,i), key2(m,1)))
if (degree_x2 > 4) then
cycle outer
endif
endif
if ( key1(m,2,i) /= key2(m,2)) then
degree_x2 = degree_x2+ popcnt(xor( key1(m,2,i), key2(m,2)))
if (degree_x2 > 4) then
cycle outer
endif
endif
enddo
if (degree_x2 > 4) then
cycle
else if(degree_x2 .ne. 0)then
if(degree_x2 /= 0)then
idx(l) = i
l = l+1
endif
enddo
enddo outer
endif
idx(0) = l-1

145
src/Determinants/slater_rules.irp.f
View File

@ -349,6 +349,80 @@ subroutine get_mono_excitation(det1,det2,exc,phase,Nint)
enddo
end
subroutine bitstring_to_list_ab( string, list, n_elements, Nint)
use bitmasks
implicit none
BEGIN_DOC
! Gives the inidices(+1) of the bits set to 1 in the bit string
! For alpha/beta determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: string(Nint,2)
integer, intent(out) :: list(Nint*bit_kind_size,2)
integer, intent(out) :: n_elements(2)
integer :: i, j, ishift
integer(bit_kind) :: l
n_elements(1) = 0
n_elements(2) = 0
ishift = 1
do i=1,Nint
l = string(i,1)
do while (l /= 0_bit_kind)
j = trailz(l)
n_elements(1) = n_elements(1)+1
l = ibclr(l,j)
list(n_elements(1),1) = ishift+j
enddo
l = string(i,2)
do while (l /= 0_bit_kind)
j = trailz(l)
n_elements(2) = n_elements(2)+1
l = ibclr(l,j)
list(n_elements(2),2) = ishift+j
enddo
ishift = ishift + bit_kind_size
enddo
end
subroutine bitstring_to_list_ab_old( string, list, n_elements, Nint)
use bitmasks
implicit none
BEGIN_DOC
! Gives the inidices(+1) of the bits set to 1 in the bit string
! For alpha/beta determinants
END_DOC
integer, intent(in) :: Nint
integer(bit_kind), intent(in) :: string(Nint,2)
integer, intent(out) :: list(Nint*bit_kind_size,2)
integer, intent(out) :: n_elements(2)
integer :: i, ishift
integer(bit_kind) :: l
n_elements(1) = 0
n_elements(2) = 0
ishift = 2
do i=1,Nint
l = string(i,1)
do while (l /= 0_bit_kind)
n_elements(1) = n_elements(1)+1
list(n_elements(1),1) = ishift+popcnt(l-1_bit_kind) - popcnt(l)
l = iand(l,l-1_bit_kind)
enddo
l = string(i,2)
do while (l /= 0_bit_kind)
n_elements(2) = n_elements(2)+1
list(n_elements(2),2) = ishift+popcnt(l-1_bit_kind) - popcnt(l)
l = iand(l,l-1_bit_kind)
enddo
ishift = ishift + bit_kind_size
enddo
end
@ -370,7 +444,7 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem, phase,phase_2
integer :: n_occ_alpha, n_occ_beta
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
@ -422,8 +496,8 @@ subroutine i_H_j(key_i,key_j,Nint,hij)
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
@ -506,7 +580,7 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem
integer :: n_occ_alpha, n_occ_beta
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
@ -558,8 +632,8 @@ subroutine i_H_j_phase_out(key_i,key_j,Nint,hij,phase,exc,degree)
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
@ -642,7 +716,7 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
integer :: i,j,k
integer :: occ(Nint*bit_kind_size,2)
double precision :: diag_H_mat_elem, phase,phase_2
integer :: n_occ_alpha, n_occ_beta
integer :: n_occ_ab(2)
logical :: has_mipi(Nint*bit_kind_size)
double precision :: mipi(Nint*bit_kind_size), miip(Nint*bit_kind_size)
PROVIDE mo_bielec_integrals_in_map mo_integrals_map
@ -696,8 +770,8 @@ subroutine i_H_j_verbose(key_i,key_j,Nint,hij,hmono,hdouble)
endif
case (1)
call get_mono_excitation(key_i,key_j,exc,phase,Nint)
call bitstring_to_list(key_i(1,1), occ(1,1), n_occ_alpha, Nint)
call bitstring_to_list(key_i(1,2), occ(1,2), n_occ_beta, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key_i, occ, n_occ_ab, Nint)
has_mipi = .False.
if (exc(0,1,1) == 1) then
! Mono alpha
@ -1229,15 +1303,15 @@ double precision function diag_H_mat_elem(det_in,Nint)
endif
!call debug_det(det_in,Nint)
integer :: tmp
call bitstring_to_list(particle(1,1), occ_particle(1,1), tmp, Nint)
ASSERT (tmp == nexc(1))
call bitstring_to_list(particle(1,2), occ_particle(1,2), tmp, Nint)
ASSERT (tmp == nexc(2))
call bitstring_to_list(hole(1,1), occ_hole(1,1), tmp, Nint)
ASSERT (tmp == nexc(1))
call bitstring_to_list(hole(1,2), occ_hole(1,2), tmp, Nint)
ASSERT (tmp == nexc(2))
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(particle, occ_particle, tmp, Nint)
ASSERT (tmp(1) == nexc(1))
ASSERT (tmp(2) == nexc(2))
!DIR$ FORCEINLINE
call bitstring_to_list_ab(hole, occ_hole, tmp, Nint)
ASSERT (tmp(1) == nexc(1))
ASSERT (tmp(2) == nexc(2))
det_tmp = ref_bitmask
do ispin=1,2
@ -1266,6 +1340,7 @@ subroutine a_operator(iorb,ispin,key,hjj,Nint,na,nb)
integer :: occ(Nint*bit_kind_size,2)
integer :: other_spin
integer :: k,l,i
integer :: tmp(2)
ASSERT (iorb > 0)
ASSERT (ispin > 0)
@ -1279,19 +1354,19 @@ subroutine a_operator(iorb,ispin,key,hjj,Nint,na,nb)
other_spin = iand(ispin,1)+1
!DIR$ FORCEINLINE
call get_occ_from_key(key,occ,Nint)
na -= 1
call bitstring_to_list_ab(key, occ, tmp, Nint)
na = na-1
hjj -= mo_mono_elec_integral(iorb,iorb)
hjj = hjj - mo_mono_elec_integral(iorb,iorb)
! Same spin
do i=1,na
hjj -= mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
hjj = hjj - mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
hjj -= mo_bielec_integral_jj(occ(i,other_spin),iorb)
hjj = hjj - mo_bielec_integral_jj(occ(i,other_spin),iorb)
enddo
end
@ -1317,13 +1392,11 @@ subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb)
ASSERT (ispin < 3)
ASSERT (Nint > 0)
integer :: tmp
integer :: tmp(2)
!DIR$ FORCEINLINE
call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint)
ASSERT (tmp == elec_alpha_num)
!DIR$ FORCEINLINE
call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint)
ASSERT (tmp == elec_beta_num)
call bitstring_to_list_ab(key, occ, tmp, Nint)
ASSERT (tmp(1) == elec_alpha_num)
ASSERT (tmp(2) == elec_beta_num)
k = ishft(iorb-1,-bit_kind_shift)+1
ASSERT (k > 0)
@ -1331,18 +1404,18 @@ subroutine ac_operator(iorb,ispin,key,hjj,Nint,na,nb)
key(k,ispin) = ibset(key(k,ispin),l)
other_spin = iand(ispin,1)+1
hjj += mo_mono_elec_integral(iorb,iorb)
hjj = hjj + mo_mono_elec_integral(iorb,iorb)
! Same spin
do i=1,na
hjj += mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
hjj = hjj + mo_bielec_integral_jj_anti(occ(i,ispin),iorb)
enddo
! Opposite spin
do i=1,nb
hjj += mo_bielec_integral_jj(occ(i,other_spin),iorb)
hjj = hjj + mo_bielec_integral_jj(occ(i,other_spin),iorb)
enddo
na += 1
na = na+1
end
subroutine get_occ_from_key(key,occ,Nint)
@ -1354,10 +1427,10 @@ subroutine get_occ_from_key(key,occ,Nint)
integer(bit_kind), intent(in) :: key(Nint,2)
integer , intent(in) :: Nint
integer , intent(out) :: occ(Nint*bit_kind_size,2)
integer :: tmp
integer :: tmp(2)
call bitstring_to_list(key(1,1), occ(1,1), tmp, Nint)
call bitstring_to_list(key(1,2), occ(1,2), tmp, Nint)
!DIR$ FORCEINLINE
call bitstring_to_list_ab(key, occ, tmp, Nint)
end