LCPQ/qp2
9
1
Fork 0
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-07-11 05:33:34 +02:00
Code Releases Activity

Merge branch 'dev' into csf

Browse Source
This commit is contained in:
Anthony Scemama 2021-01-25 18:41:48 +01:00
commit 091bc0e13c
14 changed files with 220 additions and 126 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
RELEASE_NOTES.org
View File

@ -53,6 +53,7 @@
- Added ~print_energy~ - Added ~print_energy~
- Added ~print_hamiltonian~ - Added ~print_hamiltonian~
- Added input for two body RDM - Added input for two body RDM
- Added keyword ~save_wf_after_selection~
*** Code *** Code
@ -79,3 +80,4 @@

6
docs/source/users_guide/quickstart.rst
View File

@ -128,6 +128,12 @@ and the atomic basis set:
ao_two_e_erf_ints density_for_dft electrons mo_two_e_ints scf_utils ao_two_e_erf_ints density_for_dft electrons mo_two_e_ints scf_utils
ao_two_e_ints determinants ezfio nuclei work ao_two_e_ints determinants ezfio nuclei work
If you need to run using an already existing EZFIO database, use
.. code:: bash
qp set_file hcn
Run a Hartree-Fock calculation Run a Hartree-Fock calculation
------------------------------ ------------------------------

6
src/cipsi/EZFIO.cfg
View File

@ -4,6 +4,12 @@ doc: If true, computes the one- and two-body rdms with perturbation theory
interface: ezfio,provider,ocaml interface: ezfio,provider,ocaml
default: False default: False
[save_wf_after_selection]
type: logical
doc: If true, saves the wave function after the selection, before the diagonalization
interface: ezfio,provider,ocaml
default: False
[seniority_max] [seniority_max]
type: integer type: integer
doc: Maximum number of allowed open shells. Using -1 selects all determinants doc: Maximum number of allowed open shells. Using -1 selects all determinants

5
src/cipsi/cipsi.irp.f
View File

@ -114,7 +114,10 @@ subroutine run_cipsi
! Add selected determinants ! Add selected determinants
call copy_H_apply_buffer_to_wf() call copy_H_apply_buffer_to_wf()
! call save_wavefunction
if (save_wf_after_selection) then
call save_wavefunction
endif
PROVIDE psi_coef PROVIDE psi_coef
PROVIDE psi_det PROVIDE psi_det

72
src/cipsi/selection.irp.f
View File

@ -166,6 +166,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
integer :: l_a, nmax, idx integer :: l_a, nmax, idx
integer, allocatable :: indices(:), exc_degree(:), iorder(:) integer, allocatable :: indices(:), exc_degree(:), iorder(:)
double precision, parameter :: norm_thr = 1.d-16
allocate (indices(N_det), & allocate (indices(N_det), &
exc_degree(max(N_det_alpha_unique,N_det_beta_unique))) exc_degree(max(N_det_alpha_unique,N_det_beta_unique)))
@ -185,7 +186,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
i = psi_bilinear_matrix_rows(l_a) i = psi_bilinear_matrix_rows(l_a)
if (nt + exc_degree(i) <= 4) then if (nt + exc_degree(i) <= 4) then
idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a)) idx = psi_det_sorted_order(psi_bilinear_matrix_order(l_a))
if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
indices(k) = idx indices(k) = idx
k=k+1 k=k+1
endif endif
@ -212,7 +213,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
idx = psi_det_sorted_order( & idx = psi_det_sorted_order( &
psi_bilinear_matrix_order( & psi_bilinear_matrix_order( &
psi_bilinear_matrix_transp_order(l_a))) psi_bilinear_matrix_transp_order(l_a)))
if (psi_average_norm_contrib_sorted(idx) > 1.d-20) then if (psi_average_norm_contrib_sorted(idx) > norm_thr) then
indices(k) = idx indices(k) = idx
k=k+1 k=k+1
endif endif
@ -742,12 +743,9 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
alpha_h_psi = mat(istate, p1, p2) alpha_h_psi = mat(istate, p1, p2)
do jstate=1,N_states pt2_data % overlap(:,istate) = pt2_data % overlap(:,istate) + coef(:) * coef(istate)
pt2_data % overlap(jstate,istate) += coef(jstate) * coef(istate) pt2_data % variance(istate) = pt2_data % variance(istate) + alpha_h_psi * alpha_h_psi
enddo pt2_data % pt2(istate) = pt2_data % pt2(istate) + e_pert(istate)
pt2_data % variance(istate) += alpha_h_psi * alpha_h_psi
pt2_data % pt2(istate) += e_pert(istate)
!!!DEBUG !!!DEBUG
! delta_E = E0(istate) - Hii + E_shift ! delta_E = E0(istate) - Hii + E_shift
@ -1578,7 +1576,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int) phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = mo_two_e_integral(p1, p2, h1, h2) * phase hij = mo_two_e_integral(p1, p2, h1, h2) * phase
end if end if
mat(:, p1, p2) += coefs(:) * hij mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
end do end do
end do end do
else ! AA BB else ! AA BB
@ -1595,7 +1593,7 @@ subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
else else
hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int) hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj))* get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
end if end if
mat(:, puti, putj) += coefs(:) * hij mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do end do
end do end do
end if end if
@ -1654,18 +1652,18 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
do putj=1, hfix-1 do putj=1, hfix-1
if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states) tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do end do
do putj=hfix+1, mo_num do putj=hfix+1, mo_num
if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle if(lbanned(putj, ma) .or. banned(putj, puti,bant)) cycle
hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
tmp_row(1:N_states,putj) += hij * coefs(1:N_states) tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do end do
if(ma == 1) then if(ma == 1) then
mat(1:N_states,1:mo_num,puti) += tmp_row(1:N_states,1:mo_num) mat(1:N_states,1:mo_num,puti) = mat(1:N_states,1:mo_num,puti) + tmp_row(1:N_states,1:mo_num)
else else
mat(1:N_states,puti,1:mo_num) += tmp_row(1:N_states,1:mo_num) mat(1:N_states,puti,1:mo_num) = mat(1:N_states,puti,1:mo_num) + tmp_row(1:N_states,1:mo_num)
end if end if
end if end if
@ -1679,22 +1677,22 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
putj = p1 putj = p1
if(.not. banned(putj,puti,bant)) then if(.not. banned(putj,puti,bant)) then
hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int) hij = mo_two_e_integral(p2,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p2, puti, pfix, N_int)
tmp_row(:,puti) += hij * coefs(:) tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
end if end if
putj = p2 putj = p2
if(.not. banned(putj,puti,bant)) then if(.not. banned(putj,puti,bant)) then
hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int) hij = mo_two_e_integral(p1,pfix,hfix,puti) * get_phase_bi(phasemask, ma, mi, hfix, p1, puti, pfix, N_int)
tmp_row2(:,puti) += hij * coefs(:) tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
end if end if
end do end do
if(mi == 1) then if(mi == 1) then
mat(:,:,p1) += tmp_row(:,:) mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:)
mat(:,:,p2) += tmp_row2(:,:) mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:)
else else
mat(:,p1,:) += tmp_row(:,:) mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:)
mat(:,p2,:) += tmp_row2(:,:) mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:)
end if end if
else else
if(p(0,ma) == 3) then if(p(0,ma) == 3) then
@ -1707,16 +1705,16 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
do putj=1,hfix-1 do putj=1,hfix-1
if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int) hij = (mo_two_e_integral(p1, p2, putj, hfix)-mo_two_e_integral(p2,p1,putj,hfix)) * get_phase_bi(phasemask, ma, ma, putj, p1, hfix, p2, N_int)
tmp_row(:,putj) += hij * coefs(:) tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
end do end do
do putj=hfix+1,mo_num do putj=hfix+1,mo_num
if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle if(lbanned(putj,ma) .or. banned(puti,putj,1)) cycle
hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int) hij = (mo_two_e_integral(p1, p2, hfix, putj)-mo_two_e_integral(p2,p1,hfix,putj)) * get_phase_bi(phasemask, ma, ma, hfix, p1, putj, p2, N_int)
tmp_row(:,putj) += hij * coefs(:) tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
end do end do
mat(:, :puti-1, puti) += tmp_row(:,:puti-1) mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1)
mat(:, puti, puti:) += tmp_row(:,puti:) mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row(:,puti:)
end do end do
else else
hfix = h(1,mi) hfix = h(1,mi)
@ -1730,19 +1728,19 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
putj = p2 putj = p2
if(.not. banned(puti,putj,1)) then if(.not. banned(puti,putj,1)) then
hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int) hij = mo_two_e_integral(pfix, p1, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p1, N_int)
tmp_row(:,puti) += hij * coefs(:) tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
end if end if
putj = p1 putj = p1
if(.not. banned(puti,putj,1)) then if(.not. banned(puti,putj,1)) then
hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int) hij = mo_two_e_integral(pfix, p2, hfix, puti) * get_phase_bi(phasemask, mi, ma, hfix, pfix, puti, p2, N_int)
tmp_row2(:,puti) += hij * coefs(:) tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
end if end if
end do end do
mat(:,:p2-1,p2) += tmp_row(:,:p2-1) mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1)
mat(:,p2,p2:) += tmp_row(:,p2:) mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row(:,p2:)
mat(:,:p1-1,p1) += tmp_row2(:,:p1-1) mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1)
mat(:,p1,p1:) += tmp_row2(:,p1:) mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row2(:,p1:)
end if end if
end if end if
deallocate(lbanned) deallocate(lbanned)
@ -1765,7 +1763,7 @@ subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle if(bannedOrb(p1, s1) .or. bannedOrb(p2, s2) .or. banned(p1, p2, 1)) cycle
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int) call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
call i_h_j(gen, det, N_int, hij) call i_h_j(gen, det, N_int, hij)
mat(:, p1, p2) += coefs(:) * hij mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
end do end do
end do end do
end end
@ -1818,9 +1816,9 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2, N_int)
if(ma == 1) then if(ma == 1) then
mat(:, putj, puti) += coefs(:) * hij mat(:, putj, puti) = mat(:, putj, puti) + coefs(:) * hij
else else
mat(:, puti, putj) += coefs(:) * hij mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end if end if
end do end do
else else
@ -1836,7 +1834,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p1 = p(turn2(i), 1) p1 = p(turn2(i), 1)
hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2,N_int)
mat(:, puti, putj) += coefs(:) * hij mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do end do
end do end do
end if end if
@ -1856,7 +1854,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p1 = p(i1, ma) p1 = p(i1, ma)
p2 = p(i2, ma) p2 = p(i2, ma)
hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, ma, ma, h1, p1, h2, p2,N_int)
mat(:, puti, putj) += coefs(:) * hij mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do end do
end do end do
else if(tip == 3) then else if(tip == 3) then
@ -1870,7 +1868,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
p2 = p(i, ma) p2 = p(i, ma)
hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int) hij = mo_two_e_integral(p1, p2, h1, h2) * get_phase_bi(phasemask, mi, ma, h1, p1, h2, p2,N_int)
mat(:, min(puti, putj), max(puti, putj)) += coefs(:) * hij mat(:, min(puti, putj), max(puti, putj)) = mat(:, min(puti, putj), max(puti, putj)) + coefs(:) * hij
end do end do
else ! tip == 4 else ! tip == 4
puti = p(1, sp) puti = p(1, sp)
@ -1881,7 +1879,7 @@ subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p,
h1 = h(1, mi) h1 = h(1, mi)
h2 = h(2, mi) h2 = h(2, mi)
hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int) hij = (mo_two_e_integral(p1, p2, h1, h2) - mo_two_e_integral(p2,p1, h1, h2)) * get_phase_bi(phasemask, mi, mi, h1, p1, h2, p2,N_int)
mat(:, puti, putj) += coefs(:) * hij mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end if end if
end if end if
end if end if

4
src/cipsi/stochastic_cipsi.irp.f
View File

@ -104,7 +104,9 @@ subroutine run_stochastic_cipsi
! Add selected determinants ! Add selected determinants
call copy_H_apply_buffer_to_wf() call copy_H_apply_buffer_to_wf()
! call save_wavefunction if (save_wf_after_selection) then
call save_wavefunction
endif
PROVIDE psi_coef PROVIDE psi_coef
PROVIDE psi_det PROVIDE psi_det

2
src/davidson/davidson_parallel.irp.f
View File

@ -428,7 +428,7 @@ subroutine H_S2_u_0_nstates_zmq(v_0,s_0,u_0,N_st,sze)
integer :: istep, imin, imax, ishift, ipos integer :: istep, imin, imax, ishift, ipos
integer, external :: add_task_to_taskserver integer, external :: add_task_to_taskserver
integer, parameter :: tasksize=10000 integer, parameter :: tasksize=20000
character*(100000) :: task character*(100000) :: task
istep=1 istep=1
ishift=0 ishift=0

176
src/davidson/u0_h_u0.irp.f
View File

@ -211,6 +211,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
double precision :: rss, mem, ratio double precision :: rss, mem, ratio
double precision, allocatable :: utl(:,:) double precision, allocatable :: utl(:,:)
integer, parameter :: block_size=128 integer, parameter :: block_size=128
logical :: u_is_sparse
! call resident_memory(rss) ! call resident_memory(rss)
! mem = dble(singles_beta_csc_size) / 1024.d0**3 ! mem = dble(singles_beta_csc_size) / 1024.d0**3
@ -222,6 +223,7 @@ subroutine H_S2_u_0_nstates_openmp_work_$N_int(v_t,s_t,u_t,N_st,sze,istart,iend,
! endif ! endif
compute_singles=.True. compute_singles=.True.
maxab = max(N_det_alpha_unique, N_det_beta_unique)+1 maxab = max(N_det_alpha_unique, N_det_beta_unique)+1
allocate(idx0(maxab)) allocate(idx0(maxab))
@ -249,7 +251,7 @@ compute_singles=.True.
!$OMP singles_alpha_csc,singles_alpha_csc_idx, & !$OMP singles_alpha_csc,singles_alpha_csc_idx, &
!$OMP singles_beta_csc,singles_beta_csc_idx) & !$OMP singles_beta_csc,singles_beta_csc_idx) &
!$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, & !$OMP PRIVATE(krow, kcol, tmp_det, spindet, k_a, k_b, i, &
!$OMP lcol, lrow, l_a, l_b, utl, kk, & !$OMP lcol, lrow, l_a, l_b, utl, kk, u_is_sparse, &
!$OMP buffer, doubles, n_doubles, umax, & !$OMP buffer, doubles, n_doubles, umax, &
!$OMP tmp_det2, hij, sij, idx, l, kcol_prev, & !$OMP tmp_det2, hij, sij, idx, l, kcol_prev, &
!$OMP singles_a, n_singles_a, singles_b, ratio, & !$OMP singles_a, n_singles_a, singles_b, ratio, &
@ -266,6 +268,22 @@ compute_singles=.True.
kcol_prev=-1 kcol_prev=-1
! Check if u has multiple zeros
kk=1 ! Avoid division by zero
!$OMP DO
do k=1,N_det
umax = 0.d0
do l=1,N_st
umax = max(umax, dabs(u_t(l,k)))
enddo
if (umax < 1.d-20) then
!$OMP ATOMIC
kk = kk+1
endif
enddo
!$OMP END DO
u_is_sparse = N_det / kk < 20 ! 5%
ASSERT (iend <= N_det) ASSERT (iend <= N_det)
ASSERT (istart > 0) ASSERT (istart > 0)
ASSERT (istep > 0) ASSERT (istep > 0)
@ -405,16 +423,26 @@ compute_singles=.True.
do k = 1,n_singles_a,block_size do k = 1,n_singles_a,block_size
umax = 0.d0 umax = 0.d0
! Prefetch u_t(:,l_a) ! Prefetch u_t(:,l_a)
do kk=0,block_size-1 if (u_is_sparse) then
if (k+kk > n_singles_a) exit do kk=0,block_size-1
l_a = singles_a(k+kk) if (k+kk > n_singles_a) exit
ASSERT (l_a <= N_det) l_a = singles_a(k+kk)
ASSERT (l_a <= N_det)
do l=1,N_st do l=1,N_st
utl(l,kk+1) = u_t(l,l_a) utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1))) umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo enddo
enddo else
do kk=0,block_size-1
if (k+kk > n_singles_a) exit
l_a = singles_a(k+kk)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle if (umax < 1.d-20) cycle
do kk=0,block_size-1 do kk=0,block_size-1
@ -497,16 +525,26 @@ compute_singles=.True.
do i=1,n_singles_a,block_size do i=1,n_singles_a,block_size
umax = 0.d0 umax = 0.d0
! Prefetch u_t(:,l_a) ! Prefetch u_t(:,l_a)
do kk=0,block_size-1 if (u_is_sparse) then
if (i+kk > n_singles_a) exit do kk=0,block_size-1
l_a = singles_a(i+kk) if (i+kk > n_singles_a) exit
ASSERT (l_a <= N_det) l_a = singles_a(i+kk)
ASSERT (l_a <= N_det)
do l=1,N_st do l=1,N_st
utl(l,kk+1) = u_t(l,l_a) utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1))) umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo enddo
enddo else
do kk=0,block_size-1
if (i+kk > n_singles_a) exit
l_a = singles_a(i+kk)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle if (umax < 1.d-20) cycle
do kk=0,block_size-1 do kk=0,block_size-1
@ -534,16 +572,26 @@ compute_singles=.True.
do i=1,n_doubles,block_size do i=1,n_doubles,block_size
umax = 0.d0 umax = 0.d0
! Prefetch u_t(:,l_a) ! Prefetch u_t(:,l_a)
do kk=0,block_size-1 if (u_is_sparse) then
if (i+kk > n_doubles) exit do kk=0,block_size-1
l_a = doubles(i+kk) if (i+kk > n_doubles) exit
ASSERT (l_a <= N_det) l_a = doubles(i+kk)
ASSERT (l_a <= N_det)
do l=1,N_st do l=1,N_st
utl(l,kk+1) = u_t(l,l_a) utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1))) umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo enddo
enddo else
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_a = doubles(i+kk)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle if (umax < 1.d-20) cycle
do kk=0,block_size-1 do kk=0,block_size-1
@ -611,19 +659,30 @@ compute_singles=.True.
!DIR$ LOOP COUNT avg(1000) !DIR$ LOOP COUNT avg(1000)
do i=1,n_singles_b,block_size do i=1,n_singles_b,block_size
umax = 0.d0 umax = 0.d0
do kk=0,block_size-1 if (u_is_sparse) then
if (i+kk > n_singles_b) exit do kk=0,block_size-1
l_b = singles_b(i+kk) if (i+kk > n_singles_b) exit
ASSERT (l_b <= N_det) l_b = singles_b(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
l_a = psi_bilinear_matrix_transp_order(l_b) do l=1,N_st
ASSERT (l_a <= N_det) utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
do l=1,N_st enddo
utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1)))
enddo enddo
enddo else
do kk=0,block_size-1
if (i+kk > n_singles_b) exit
l_b = singles_b(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle if (umax < 1.d-20) cycle
do kk=0,block_size-1 do kk=0,block_size-1
@ -649,18 +708,29 @@ compute_singles=.True.
!DIR$ LOOP COUNT avg(50000) !DIR$ LOOP COUNT avg(50000)
do i=1,n_doubles,block_size do i=1,n_doubles,block_size
umax = 0.d0 umax = 0.d0
do kk=0,block_size-1 if (u_is_sparse) then
if (i+kk > n_doubles) exit do kk=0,block_size-1
l_b = doubles(i+kk) if (i+kk > n_doubles) exit
ASSERT (l_b <= N_det) l_b = doubles(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b) l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_a <= N_det) ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
do l=1,N_st do l=1,N_st
utl(l,kk+1) = u_t(l,l_a) utl(l,kk+1) = u_t(l,l_a)
umax = max(umax, dabs(utl(l,kk+1))) umax = max(umax, dabs(utl(l,kk+1)))
enddo
enddo enddo
enddo else
do kk=0,block_size-1
if (i+kk > n_doubles) exit
l_b = doubles(i+kk)
l_a = psi_bilinear_matrix_transp_order(l_b)
ASSERT (l_b <= N_det)
ASSERT (l_a <= N_det)
utl(:,kk+1) = u_t(:,l_a)
enddo
umax = 1.d0
endif
if (umax < 1.d-20) cycle if (umax < 1.d-20) cycle
do kk=0,block_size-1 do kk=0,block_size-1
@ -688,10 +758,14 @@ compute_singles=.True.
! Initial determinant is at k_a in alpha-major representation ! Initial determinant is at k_a in alpha-major representation
! ----------------------------------------------------------------------- ! -----------------------------------------------------------------------
umax = 0.d0 if (u_is_sparse) then
do l=1,N_st umax = 0.d0
umax = max(umax, dabs(u_t(l,k_a))) do l=1,N_st
enddo umax = max(umax, dabs(u_t(l,k_a)))
enddo
else
umax = 1.d0
endif
if (umax < 1.d-20) cycle if (umax < 1.d-20) cycle
krow = psi_bilinear_matrix_rows(k_a) krow = psi_bilinear_matrix_rows(k_a)

40
src/determinants/configurations.irp.f
View File

@ -42,6 +42,7 @@ subroutine configuration_to_dets_size(o,sze,n_alpha,Nint)
amax -= popcnt( o(k,2) ) amax -= popcnt( o(k,2) )
enddo enddo
if (binom_int(bmax, amax) > huge(1)) then if (binom_int(bmax, amax) > huge(1)) then
print *, bmax, amax
print *, irp_here, ': Too many determinants to generate' print *, irp_here, ': Too many determinants to generate'
stop 1 stop 1
endif endif
@ -340,7 +341,7 @@ BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ]
! Returns the index of the configuration for each determinant ! Returns the index of the configuration for each determinant
END_DOC END_DOC
integer :: i,j,k,r,l integer :: i,j,k,r,l
integer*8 :: key integer*8 :: key, key2
integer(bit_kind) :: occ(N_int,2) integer(bit_kind) :: occ(N_int,2)
logical :: found logical :: found
integer*8, allocatable :: bit_tmp(:) integer*8, allocatable :: bit_tmp(:)
@ -361,36 +362,23 @@ BEGIN_PROVIDER [ integer, det_to_configuration, (N_det) ]
key = configuration_search_key(occ,N_int) key = configuration_search_key(occ,N_int)
! TODO: Binary search l = 0
l = 1 r = N_configuration+1
r = N_configuration j = shiftr(r-l,1)
! do while(r-l > 32) do while (j>=1)
! j = shiftr(r+l,1) j = j+l
! if (bit_tmp(j) < key) then key2 = configuration_search_key(psi_configuration(1,1,j),N_int)
! l = j if (key2 == key) then
! else
! r = j
! endif
! enddo
do j=l,r
found = .True.
do k=1,N_int
if ( (occ(k,1) /= psi_configuration(k,1,j)) &
.or. (occ(k,2) /= psi_configuration(k,2,j)) ) then
found = .False.
exit
endif
enddo
if (found) then
det_to_configuration(i) = j det_to_configuration(i) = j
exit exit
else if (key2 > key) then
r = j
else
l = j
endif endif
j = shiftr(r-l,1)
enddo enddo
if (.not.found) then
print *, '3 bug in ', irp_here
stop -1
endif
enddo enddo
!$OMP END PARALLEL DO !$OMP END PARALLEL DO
deallocate(bit_tmp) deallocate(bit_tmp)

26
src/mo_two_e_ints/map_integrals.irp.f
View File

@ -99,9 +99,15 @@ double precision function get_two_e_integral(i,j,k,l,map)
type(map_type), intent(inout) :: map type(map_type), intent(inout) :: map
real(integral_kind) :: tmp real(integral_kind) :: tmp
PROVIDE mo_two_e_integrals_in_map mo_integrals_cache PROVIDE mo_two_e_integrals_in_map mo_integrals_cache
if (banned_excitation(i,k) .or. banned_excitation(j,l)) then if (use_banned_excitation) then
get_two_e_integral = 0.d0 if (banned_excitation(i,k)) then
return get_two_e_integral = 0.d0
return
endif
if (banned_excitation(j,l)) then
get_two_e_integral = 0.d0
return
endif
endif endif
ii = l-mo_integrals_cache_min ii = l-mo_integrals_cache_min
ii = ior(ii, k-mo_integrals_cache_min) ii = ior(ii, k-mo_integrals_cache_min)
@ -282,17 +288,19 @@ subroutine get_mo_two_e_integrals_exch_ii(k,l,sze,out_val,map)
end end
BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ] BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ]
&BEGIN_PROVIDER [ logical, use_banned_excitation ]
implicit none implicit none
use map_module use map_module
BEGIN_DOC BEGIN_DOC
! If true, the excitation is banned in the selection. Useful with local MOs. ! If true, the excitation is banned in the selection. Useful with local MOs.
END_DOC END_DOC
banned_excitation = .False. banned_excitation = .False.
integer :: i,j integer :: i,j, icount
integer(key_kind) :: idx integer(key_kind) :: idx
double precision :: tmp double precision :: tmp
! double precision :: buffer(mo_num)
icount = 1 ! Avoid division by zero
do j=1,mo_num do j=1,mo_num
do i=1,j-1 do i=1,j-1
call two_e_integrals_index(i,j,j,i,idx) call two_e_integrals_index(i,j,j,i,idx)
@ -300,8 +308,14 @@ BEGIN_PROVIDER [ logical, banned_excitation, (mo_num,mo_num) ]
call map_get(mo_integrals_map,idx,tmp) call map_get(mo_integrals_map,idx,tmp)
banned_excitation(i,j) = dabs(tmp) < 1.d-14 banned_excitation(i,j) = dabs(tmp) < 1.d-14
banned_excitation(j,i) = banned_excitation(i,j) banned_excitation(j,i) = banned_excitation(i,j)
if (banned_excitation(i,j)) icount = icount+2
enddo enddo
enddo enddo
use_banned_excitation = (mo_num*mo_num) / icount <= 100 !1%
if (use_banned_excitation) then
print *, 'Using sparsity of exchange integrals'
endif
END_PROVIDER END_PROVIDER

1
src/two_body_rdm/example.irp.f
View File

@ -284,3 +284,4 @@ subroutine routine_full_mos
print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3 print*,'wee_tot_st_av_3 = ',wee_tot_st_av_3
end end

2
travis/compilation.sh
View File

@ -8,7 +8,7 @@ tar -zxf $HOME/cache/config.tgz
# Configure QP2 # Configure QP2
cd qp2 cd qp2
source ./quantum_package.rc source ./quantum_package.rc
ninja -j 1 -v ninja -j 1 -v || exit -1
# Create cache # Create cache
cd .. cd ..

2
travis/configuration.sh
View File

@ -2,7 +2,7 @@
# Stage 1 # Stage 1
# Configure QP2 # Configure QP2
./configure --install all --config ./config/travis.cfg ./configure --install all --config ./config/travis.cfg || exit -1
# Create cache # Create cache
cd ../ cd ../

2
travis/testing.sh
View File

@ -8,7 +8,7 @@ tar -zxf $HOME/cache/compil.tgz
# Configure QP2 # Configure QP2
cd qp2 cd qp2
source ./quantum_package.rc source ./quantum_package.rc
qp_test -a && rm $HOME/cache/compil.tgz exec qp_test -a && rm $HOME/cache/compil.tgz