9
1
mirror of https://github.com/QuantumPackage/qp2.git synced 2024-10-10 09:51:31 +02:00

Merge pull request #19 from QuantumPackage/dev-stable

Dev stable
This commit is contained in:
AbdAmmar 2023-06-20 20:16:16 +02:00 committed by GitHub
commit fca309a5b2
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
36 changed files with 1150 additions and 666 deletions

View File

@ -38,7 +38,8 @@ let run slave ?prefix exe ezfio_file =
| Unix.Unix_error _ -> try_new_port (port_number+100)
in
let result =
try_new_port 41279
let port = 10*(Unix.getpid () mod 2823) + 32_769 in
try_new_port port
in
Zmq.Socket.close dummy_socket;
Zmq.Context.terminate zmq_context;

View File

@ -1,7 +1,7 @@
#!/usr/bin/env python3
"""
Save the .o from a .f90
and is the .o is asked a second time, retur it
and is the .o is asked a second time, return it
Take in argv command like:
ifort -g -openmp -I IRPF90_temp/Ezfio_files/ -c IRPF90_temp/Integrals_Monoelec/kin_ao_ints.irp.module.F90 -o IRPF90_temp/Integrals_Monoelec/kin_ao_ints.irp.module.o
"""

View File

@ -104,6 +104,9 @@ BEGIN_PROVIDER [ double precision, ao_integrals_n_e, (ao_num,ao_num)]
IF(do_pseudo) THEN
ao_integrals_n_e += ao_pseudo_integrals
ENDIF
IF(point_charges) THEN
ao_integrals_n_e += ao_integrals_pt_chrg
ENDIF
endif

View File

@ -4,7 +4,7 @@ BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
! Number of Cholesky vectors in AO basis
END_DOC
cholesky_ao_num_guess = ao_num*ao_num / 2
cholesky_ao_num_guess = ao_num*ao_num
END_PROVIDER
BEGIN_PROVIDER [ integer, cholesky_ao_num ]
@ -44,19 +44,12 @@ END_PROVIDER
do m=0,9
do l=1+m,ao_num,10
!$OMP DO SCHEDULE(dynamic)
do j=1,l
do j=1,ao_num
do k=1,ao_num
do i=1,min(k,j)
do i=1,ao_num
if (ao_two_e_integral_zero(i,j,k,l)) cycle
integral = get_ao_two_e_integral(i,j,k,l, ao_integrals_map)
ao_integrals(i,k,j,l) = integral
ao_integrals(k,i,j,l) = integral
ao_integrals(i,k,l,j) = integral
ao_integrals(k,i,l,j) = integral
ao_integrals(j,l,i,k) = integral
ao_integrals(j,l,k,i) = integral
ao_integrals(l,j,i,k) = integral
ao_integrals(l,j,k,i) = integral
enddo
enddo
enddo

View File

@ -88,7 +88,6 @@ subroutine test_5idx
! if (dabs(three_e_5_idx_direct_bi_ort(m,l,j,k,i) - three_e_5_idx_exch12_bi_ort(m,l,i,k,j)) > 1.d-10) then
! stop
! endif
new = three_e_5_idx_direct_bi_ort(m,l,j,k,i)
ref = three_e_5_idx_direct_bi_ort_old(m,l,j,k,i)
contrib = dabs(new - ref)
@ -118,6 +117,51 @@ subroutine test_5idx
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle1'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_cycle_2_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_2_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle2'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch23_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch23_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch23'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_exch13_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_exch13_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'exch13'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
! endif
!
! new = three_e_5_idx_cycle_1_bi_ort(m,l,j,k,i)
! ref = three_e_5_idx_cycle_1_bi_ort_old(m,l,j,k,i)
! contrib = dabs(new - ref)
! accu += contrib
! if(contrib .gt. 1.d-10)then
! print*,'cycle1'
! print*,i,k,j,l,m
! print*,ref,new,contrib
! stop
@ -155,7 +199,7 @@ subroutine test_5idx
! print*,ref,new,contrib
! stop
! endif
!
enddo
enddo
enddo

View File

@ -65,7 +65,7 @@ end
tmp_mat = 0.d0
call print_memory_usage
!
do m = 1, mo_num
allocate(grad_mli(n_points_final_grid,mo_num))

11
src/ccsd/EZFIO.cfg Normal file
View File

@ -0,0 +1,11 @@
[energy]
type: double precision
doc: CCSD energy
interface: ezfio
[energy_t]
type: double precision
doc: CCSD(T) energy
interface: ezfio

View File

@ -135,8 +135,11 @@ subroutine run_ccsd_space_orb
write(*,'(A15,1pE10.2,A3)')' Conv = ', max_r
print*,''
call write_t1(nO,nV,t1)
call write_t2(nO,nV,t2)
if (write_amplitudes) then
call write_t1(nO,nV,t1)
call write_t2(nO,nV,t2)
call ezfio_set_utils_cc_io_amplitudes('Read')
endif
! Deallocation
if (cc_update_method == 'diis') then
@ -147,6 +150,7 @@ subroutine run_ccsd_space_orb
! CCSD(T)
double precision :: e_t
e_t = 0.d0
if (cc_par_t .and. elec_alpha_num + elec_beta_num > 2) then
@ -182,8 +186,7 @@ subroutine run_ccsd_space_orb
print*,''
endif
print*,'Reference determinant:'
call print_det(det,N_int)
call save_energy(uncorr_energy + energy, e_t)
deallocate(t1,t2)

View File

@ -269,8 +269,11 @@ subroutine run_ccsd_spin_orb
write(*,'(A15,1pE10.2,A3)')' Conv = ', max_r
print*,''
call write_t1(nO,nV,t1)
call write_t2(nO,nV,t2)
if (write_amplitudes) then
call write_t1(nO,nV,t1)
call write_t2(nO,nV,t2)
call ezfio_set_utils_cc_io_amplitudes('Read')
endif
! Deallocate
if (cc_update_method == 'diis') then
@ -284,8 +287,9 @@ subroutine run_ccsd_spin_orb
deallocate(v_ovoo,v_oovo)
deallocate(v_ovvo,v_ovov,v_oovv)
double precision :: t_corr
t_corr = 0.d0
if (cc_par_t .and. elec_alpha_num +elec_beta_num > 2) then
double precision :: t_corr
print*,'CCSD(T) calculation...'
call wall_time(ta)
!allocate(v_vvvo(nV,nV,nV,nO))
@ -307,8 +311,8 @@ subroutine run_ccsd_spin_orb
write(*,'(A15,F18.12,A3)') ' Correlation = ', energy + t_corr, ' Ha'
print*,''
endif
print*,'Reference determinant:'
call print_det(det,N_int)
call save_energy(uncorr_energy + energy, t_corr)
deallocate(f_oo,f_ov,f_vv,f_o,f_v)
deallocate(v_ooov,v_vvoo,t1,t2)

View File

@ -104,17 +104,17 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
integer*8, allocatable :: sampled(:)
! integer(omp_lock_kind), allocatable :: lock(:)
integer*2 , allocatable :: abc(:,:)
integer*8 :: Nabc, i8
integer*8 :: Nabc, i8,kiter
integer*8, allocatable :: iorder(:)
double precision :: eocc
double precision :: norm
integer :: kiter, isample
integer :: isample
! Prepare table of triplets (a,b,c)
Nabc = (int(nV,8) * int(nV+1,8) * int(nV+2,8))/6_8 - nV
allocate (memo(Nabc), sampled(Nabc), Pabc(Nabc), waccu(Nabc))
allocate (memo(Nabc), sampled(Nabc), Pabc(Nabc), waccu(0:Nabc))
allocate (abc(4,Nabc), iorder(Nabc)) !, lock(Nabc))
! eocc = 3.d0/dble(nO) * sum(f_o(1:nO))
@ -124,21 +124,21 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
do c = b+1, nV
Nabc = Nabc + 1_8
Pabc(Nabc) = -1.d0/(f_v(a) + f_v(b) + f_v(c))
abc(1,Nabc) = a
abc(2,Nabc) = b
abc(3,Nabc) = c
abc(1,Nabc) = int(a,2)
abc(2,Nabc) = int(b,2)
abc(3,Nabc) = int(c,2)
enddo
Nabc = Nabc + 1_8
abc(1,Nabc) = a
abc(2,Nabc) = b
abc(3,Nabc) = a
abc(1,Nabc) = int(a,2)
abc(2,Nabc) = int(b,2)
abc(3,Nabc) = int(a,2)
Pabc(Nabc) = -1.d0/(2.d0*f_v(a) + f_v(b))
Nabc = Nabc + 1_8
abc(1,Nabc) = b
abc(2,Nabc) = a
abc(3,Nabc) = b
abc(1,Nabc) = int(b,2)
abc(2,Nabc) = int(a,2)
abc(3,Nabc) = int(b,2)
Pabc(Nabc) = -1.d0/(f_v(a) + 2.d0*f_v(b))
enddo
enddo
@ -169,6 +169,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
waccu(i8) = waccu(i8+1) - Pabc(i8+1)
enddo
waccu(:) = waccu(:) + 1.d0
waccu(0) = 0.d0
logical :: converged, do_comp
double precision :: eta, variance, error, sample
@ -222,8 +223,12 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
do kiter=1,Nabc
!$OMP MASTER
do while ((imin <= Nabc).and.(sampled(imin)>-1_8))
imin = imin+1
do while (imin <= Nabc)
if (sampled(imin)>-1_8) then
imin = imin+1
else
exit
endif
enddo
! Deterministic part
@ -301,6 +306,7 @@ subroutine ccsd_par_t_space_stoch(nO,nV,t1,t2,f_o,f_v,v_vvvo,v_vvoo,v_vooo,energ
endif
enddo
isample = min(isample,nbuckets)
do ieta=bounds(1,isample), Nabc
w = dble(max(sampled(ieta),0_8))
tmp = w * memo(ieta) * Pabc(ieta)
@ -331,33 +337,39 @@ end
integer*8 function binary_search(arr, key, size)
integer*8 function binary_search(arr, key, sze)
implicit none
BEGIN_DOC
! Searches the key in array arr(1:size) between l_in and r_in, and returns its index
! Searches the key in array arr(1:sze) between l_in and r_in, and returns its index
END_DOC
integer*8 :: size, i, j, mid, l_in, r_in
double precision, dimension(size) :: arr(1:size)
integer*8 :: sze, i, j, mid
double precision :: arr(0:sze)
double precision :: key
i = 1_8
j = size
if ( key < arr(1) ) then
binary_search = 0_8
return
end if
do while (j >= i)
mid = i + (j - i) / 2
if (arr(mid) >= key) then
if (mid > 1 .and. arr(mid - 1) < key) then
binary_search = mid
return
end if
j = mid - 1
else if (arr(mid) < key) then
i = mid + 1
else
binary_search = mid + 1
return
end if
if ( key >= arr(sze) ) then
binary_search = sze
return
end if
i = 0_8
j = sze + 1_8
do while (.True.)
mid = (i + j) / 2_8
if ( key >= arr(mid) ) then
i = mid
else
j = mid
end if
if (j-i <= 1_8) then
binary_search = i
return
endif
end do
binary_search = i
end function binary_search

View File

@ -0,0 +1,13 @@
subroutine save_energy(E,ET)
implicit none
BEGIN_DOC
! Saves the energy in |EZFIO|.
END_DOC
double precision, intent(in) :: E, ET
call ezfio_set_ccsd_energy(E)
if (ET /= 0.d0) then
call ezfio_set_ccsd_energy_t(E+ET)
endif
end

View File

@ -88,6 +88,10 @@ subroutine select_connected(i_generator,E0,pt2_data,b,subset,csubset)
particle_mask(k,1) = iand(generators_bitmask(k,1,s_part), not(psi_det_generators(k,1,i_generator)) )
particle_mask(k,2) = iand(generators_bitmask(k,2,s_part), not(psi_det_generators(k,2,i_generator)) )
enddo
if ((subset == 1).and.(sum(hole_mask(:,2)) == 0_bit_kind)) then
! No beta electron to excite
call select_singles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,b)
endif
call select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_diag_tmp,E0,pt2_data,b,subset,csubset)
deallocate(fock_diag_tmp)
end subroutine
@ -142,7 +146,7 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
use selection_types
implicit none
BEGIN_DOC
! WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted
! WARNING /!\ : It is assumed that the generators and selectors are psi_det_sorted
END_DOC
integer, intent(in) :: i_generator, subset, csubset
@ -237,7 +241,6 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
enddo
! Iterate on 0S alpha, and find betas TQ such that exc_degree <= 4
! Remove also contributions < 1.d-20)
do j=1,N_det_alpha_unique
call get_excitation_degree_spin(psi_det_alpha_unique(1,j), &
psi_det_generators(1,1,i_generator), nt, N_int)
@ -480,7 +483,9 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
do s2=s1,2
sp = s1
if(s1 /= s2) sp = 3
if(s1 /= s2) then
sp = 3
endif
ib = 1
if(s1 == s2) ib = i1+1
@ -528,7 +533,10 @@ subroutine select_singles_and_doubles(i_generator,hole_mask,particle_mask,fock_d
deallocate(preinteresting, prefullinteresting, interesting, fullinteresting)
deallocate(banned, bannedOrb,mat)
end subroutine
subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
BEGIN_TEMPLATE
subroutine fill_buffer_$DOUBLE(i_generator, sp, h1, h2, bannedOrb, banned, fock_diag_tmp, E0, pt2_data, mat, buf)
use bitmasks
use selection_types
implicit none
@ -562,7 +570,20 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
s1 = sp
s2 = sp
end if
call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
if ($IS_DOUBLE) then
if (h2 == 0) then
print *, 'h2=0 in '//trim(irp_here)
stop
endif
call apply_holes(psi_det_generators(1,1,i_generator), s1, h1, s2, h2, mask, ok, N_int)
else
if (h2 /= 0) then
print *, 'h2 /= in '//trim(irp_here)
stop
endif
call apply_hole(psi_det_generators(1,1,i_generator), s1, h1, mask, ok, N_int)
endif
E_shift = 0.d0
if (h0_type == 'CFG') then
@ -570,12 +591,15 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
E_shift = psi_det_Hii(i_generator) - psi_configuration_Hii(j)
endif
do p1=1,mo_num
if(bannedOrb(p1, s1)) cycle
$DO_p1
! do p1=1,mo_num
if (bannedOrb(p1, s1)) cycle
ib = 1
if(sp /= 3) ib = p1+1
do p2=ib,mo_num
$DO_p2
! do p2=ib,mo_num
! -----
! /!\ Generating only single excited determinants doesn't work because a
@ -584,9 +608,10 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
! detected as already generated when generating in the future with a
! double excitation.
! -----
if(bannedOrb(p2, s2)) cycle
if(banned(p1,p2)) cycle
if ($IS_DOUBLE) then
if(bannedOrb(p2, s2)) cycle
if(banned(p1,p2)) cycle
endif
if(pseudo_sym)then
if(dabs(mat(1, p1, p2)).lt.thresh_sym)then
@ -596,7 +621,11 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
val = maxval(abs(mat(1:N_states, p1, p2)))
if( val == 0d0) cycle
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
if ($IS_DOUBLE) then
call apply_particles(mask, s1, p1, s2, p2, det, ok, N_int)
else
call apply_particle(mask, s1, p1, det, ok, N_int)
endif
if (do_only_cas) then
integer, external :: number_of_holes, number_of_particles
@ -797,7 +826,7 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
case(5)
! Variance selection
if (h0_type == 'CFG') then
w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) &
w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate)) &
/ c0_weight(istate)
else
w = min(w, - alpha_h_psi * alpha_h_psi * s_weight(istate,istate))
@ -857,10 +886,19 @@ subroutine fill_buffer_double(i_generator, sp, h1, h2, bannedOrb, banned, fock_d
if(w <= buf%mini) then
call add_to_selection_buffer(buf, det, w)
end if
end do
end do
! enddo
$ENDDO_p1
! enddo
$ENDDO_p2
end
SUBST [ DOUBLE , DO_p1 , ENDDO_p1 , DO_p2 , ENDDO_p2 , IS_DOUBLE ]
double ; do p1=1,mo_num ; enddo ; do p2=ib,mo_num ; enddo ; .True. ;;
single ; do p1=1,mo_num ; enddo ; p2=1 ; ; .False. ;;
END_TEMPLATE
subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, interesting)
use bitmasks
implicit none
@ -882,6 +920,7 @@ subroutine splash_pq(mask, sp, det, i_gen, N_sel, bannedOrb, banned, mat, intere
PROVIDE psi_selectors_coef_transp psi_det_sorted
mat = 0d0
p=0
do i=1,N_int
negMask(i,1) = not(mask(i,1))
@ -1435,7 +1474,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
p1 = p(1,sp)
p2 = p(2,sp)
do puti=1, mo_num
if(bannedOrb(puti, sp)) cycle
if (bannedOrb(puti, sp)) cycle
call get_mo_two_e_integrals(puti,p2,p1,mo_num,hij_cache1,mo_integrals_map)
call get_mo_two_e_integrals(puti,p1,p2,mo_num,hij_cache2,mo_integrals_map)
do putj=puti+1, mo_num
@ -1446,7 +1485,7 @@ subroutine get_d0(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
call i_h_j(gen, det, N_int, hij)
if (hij == 0.d0) cycle
else
hij = (mo_two_e_integral(p1, p2, puti, putj) - mo_two_e_integral(p2, p1, puti, putj))
hij = hij_cache1(putj) - hij_cache2(putj)
if (hij == 0.d0) cycle
hij = hij * get_phase_bi(phasemask, sp, sp, puti, p1 , putj, p2, N_int)
end if
@ -1506,7 +1545,7 @@ subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
use bitmasks
implicit none
BEGIN_DOC
! Identify the determinants in det which are in the internal space. These are
! Identify the determinants in det that are in the internal space. These are
! the determinants that can be produced by creating two particles on the mask.
END_DOC
@ -1534,7 +1573,7 @@ subroutine spot_isinwf(mask, det, i_gen, N, banned, fullMatch, interesting)
if(iand(det(j,2,i), mask(j,2)) /= mask(j, 2)) cycle genl
end do
! If det(i) < det(i_gen), it hs already been considered
! If det(i) < det(i_gen), it has already been considered
if(interesting(i) < i_gen) then
fullMatch = .true.
return
@ -1585,352 +1624,4 @@ end
! OLD unoptimized routines for debugging
! ======================================
subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
integer :: i, j, s, h1, h2, p1, p2, puti, putj
double precision :: hij, phase
double precision, external :: get_phase_bi, mo_two_e_integral
logical :: ok
integer :: bant
bant = 1
if(sp == 3) then ! AB
h1 = p(1,1)
h2 = p(1,2)
do p1=1, mo_num
if(bannedOrb(p1, 1)) cycle
do p2=1, mo_num
if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable?
if(p1 == h1 .or. p2 == h2) then
call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
call i_h_j(gen, det, N_int, hij)
else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = mo_two_e_integral(p1, p2, h1, h2) * phase
end if
mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
end do
end do
else ! AA BB
p1 = p(1,sp)
p2 = p(2,sp)
do puti=1, mo_num
if(bannedOrb(puti, sp)) cycle
do putj=puti+1, mo_num
if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable?
if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
call i_h_j(gen, det, N_int, hij)
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)
end if
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do
end do
end if
end
subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision :: hij, tmp_row(N_states, mo_num), tmp_row2(N_states, mo_num)
double precision, external :: get_phase_bi, mo_two_e_integral
logical :: ok
logical, allocatable :: lbanned(:,:)
integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j
integer :: hfix, pfix, h1, h2, p1, p2, ib
integer, parameter :: turn2(2) = (/2,1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
allocate (lbanned(mo_num, 2))
lbanned = bannedOrb
do i=1, p(0,1)
lbanned(p(i,1), 1) = .true.
end do
do i=1, p(0,2)
lbanned(p(i,2), 2) = .true.
end do
ma = 1
if(p(0,2) >= 2) ma = 2
mi = turn2(ma)
bant = 1
if(sp == 3) then
!move MA
if(ma == 2) bant = 2
puti = p(1,mi)
hfix = h(1,ma)
p1 = p(1,ma)
p2 = p(2,ma)
if(.not. bannedOrb(puti, mi)) then
tmp_row = 0d0
do putj=1, hfix-1
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)
tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do
do putj=hfix+1, mo_num
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)
tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do
if(ma == 1) then
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
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
!MOVE MI
pfix = p(1,mi)
tmp_row = 0d0
tmp_row2 = 0d0
do puti=1,mo_num
if(lbanned(puti,mi)) cycle
!p1 fixed
putj = p1
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)
tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
end if
putj = p2
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)
tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
end if
end do
if(mi == 1) then
mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:)
mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:)
else
mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:)
mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:)
end if
else
if(p(0,ma) == 3) then
do i=1,3
hfix = h(1,ma)
puti = p(i, ma)
p1 = p(turn3(1,i), ma)
p2 = p(turn3(2,i), ma)
tmp_row = 0d0
do putj=1,hfix-1
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)
tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
end do
do putj=hfix+1,mo_num
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)
tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
end do
mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1)
mat(:, puti, puti:) = mat(:, puti, puti:) + tmp_row(:,puti:)
end do
else
hfix = h(1,mi)
pfix = p(1,mi)
p1 = p(1,ma)
p2 = p(2,ma)
tmp_row = 0d0
tmp_row2 = 0d0
do puti=1,mo_num
if(lbanned(puti,ma)) cycle
putj = p2
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)
tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
end if
putj = p1
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)
tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
end if
end do
mat(:,:p2-1,p2) = mat(:,:p2-1,p2) + tmp_row(:,:p2-1)
mat(:,p2,p2:) = mat(:,p2,p2:) + tmp_row(:,p2:)
mat(:,:p1-1,p1) = mat(:,:p1-1,p1) + tmp_row2(:,:p1-1)
mat(:,p1,p1:) = mat(:,p1,p1:) + tmp_row2(:,p1:)
end if
end if
deallocate(lbanned)
!! MONO
if(sp == 3) then
s1 = 1
s2 = 2
else
s1 = sp
s2 = sp
end if
do i1=1,p(0,s1)
ib = 1
if(s1 == s2) ib = i1+1
do i2=ib,p(0,s2)
p1 = p(i1,s1)
p2 = p(i2,s2)
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 i_h_j(gen, det, N_int, hij)
mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
end do
end do
end
subroutine get_d2_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(2,N_int)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision, external :: get_phase_bi, mo_two_e_integral
integer :: i, j, tip, ma, mi, puti, putj
integer :: h1, h2, p1, p2, i1, i2
double precision :: hij, phase
integer, parameter:: turn2d(2,3,4) = reshape((/0,0, 0,0, 0,0, 3,4, 0,0, 0,0, 2,4, 1,4, 0,0, 2,3, 1,3, 1,2 /), (/2,3,4/))
integer, parameter :: turn2(2) = (/2, 1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
bant = 1
tip = p(0,1) * p(0,2)
ma = sp
if(p(0,1) > p(0,2)) ma = 1
if(p(0,1) < p(0,2)) ma = 2
mi = mod(ma, 2) + 1
if(sp == 3) then
if(ma == 2) bant = 2
if(tip == 3) then
puti = p(1, mi)
do i = 1, 3
putj = p(i, ma)
if(banned(putj,puti,bant)) cycle
i1 = turn3(1,i)
i2 = turn3(2,i)
p1 = p(i1, ma)
p2 = p(i2, ma)
h1 = h(1, ma)
h2 = h(2, 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)
if(ma == 1) then
mat(:, putj, puti) = mat(:, putj, puti) + coefs(:) * hij
else
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end if
end do
else
h1 = h(1,1)
h2 = h(1,2)
do j = 1,2
putj = p(j, 2)
p2 = p(turn2(j), 2)
do i = 1,2
puti = p(i, 1)
if(banned(puti,putj,bant)) cycle
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)
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do
end do
end if
else
if(tip == 0) then
h1 = h(1, ma)
h2 = h(2, ma)
do i=1,3
puti = p(i, ma)
do j=i+1,4
putj = p(j, ma)
if(banned(puti,putj,1)) cycle
i1 = turn2d(1, i, j)
i2 = turn2d(2, i, j)
p1 = p(i1, 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)
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do
end do
else if(tip == 3) then
h1 = h(1, mi)
h2 = h(1, ma)
p1 = p(1, mi)
do i=1,3
puti = p(turn3(1,i), ma)
putj = p(turn3(2,i), ma)
if(banned(puti,putj,1)) cycle
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)
mat(:, min(puti, putj), max(puti, putj)) = mat(:, min(puti, putj), max(puti, putj)) + coefs(:) * hij
end do
else ! tip == 4
puti = p(1, sp)
putj = p(2, sp)
if(.not. banned(puti,putj,1)) then
p1 = p(1, mi)
p2 = p(2, mi)
h1 = h(1, 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)
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end if
end if
end if
end

View File

@ -0,0 +1,350 @@
! OLD unoptimized routines for debugging
! ======================================
subroutine get_d0_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: gen(N_int, 2), mask(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
integer :: i, j, s, h1, h2, p1, p2, puti, putj
double precision :: hij, phase
double precision, external :: get_phase_bi, mo_two_e_integral
logical :: ok
integer :: bant
bant = 1
if(sp == 3) then ! AB
h1 = p(1,1)
h2 = p(1,2)
do p1=1, mo_num
if(bannedOrb(p1, 1)) cycle
do p2=1, mo_num
if(bannedOrb(p2,2)) cycle
if(banned(p1, p2, bant)) cycle ! rentable?
if(p1 == h1 .or. p2 == h2) then
call apply_particles(mask, 1,p1,2,p2, det, ok, N_int)
call i_h_j(gen, det, N_int, hij)
else
phase = get_phase_bi(phasemask, 1, 2, h1, p1, h2, p2, N_int)
hij = mo_two_e_integral(p1, p2, h1, h2) * phase
end if
mat(:, p1, p2) = mat(:, p1, p2) + coefs(:) * hij
end do
end do
else ! AA BB
p1 = p(1,sp)
p2 = p(2,sp)
do puti=1, mo_num
! do not cycle here? otherwise singles will be missed??
if(bannedOrb(puti, sp)) cycle
do putj=puti+1, mo_num
if(bannedOrb(putj, sp)) cycle
if(banned(puti, putj, bant)) cycle ! rentable?
if(puti == p1 .or. putj == p2 .or. puti == p2 .or. putj == p1) then
call apply_particles(mask, sp,puti,sp,putj, det, ok, N_int)
call i_h_j(gen, det, N_int, hij)
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)
end if
mat(:, puti, putj) = mat(:, puti, putj) + coefs(:) * hij
end do
end do
end if
end
subroutine get_d1_reference(gen, phasemask, bannedOrb, banned, mat, mask, h, p, sp, coefs)
use bitmasks
implicit none
integer(bit_kind), intent(in) :: mask(N_int, 2), gen(N_int, 2)
integer(bit_kind), intent(in) :: phasemask(N_int,2)
logical, intent(in) :: bannedOrb(mo_num, 2), banned(mo_num, mo_num,2)
integer(bit_kind) :: det(N_int, 2)
double precision, intent(in) :: coefs(N_states)
double precision, intent(inout) :: mat(N_states, mo_num, mo_num)
integer, intent(in) :: h(0:2,2), p(0:4,2), sp
double precision :: hij, tmp_row(N_states, mo_num), tmp_row2(N_states, mo_num)
double precision, external :: get_phase_bi, mo_two_e_integral
logical :: ok
logical, allocatable :: lbanned(:,:)
integer :: puti, putj, ma, mi, s1, s2, i, i1, i2, j
integer :: hfix, pfix, h1, h2, p1, p2, ib
integer, parameter :: turn2(2) = (/2,1/)
integer, parameter :: turn3(2,3) = reshape((/2,3, 1,3, 1,2/), (/2,3/))
integer :: bant
allocate (lbanned(mo_num, 2))
lbanned = bannedOrb
do i=1, p(0,1)
lbanned(p(i,1), 1) = .true.
end do
do i=1, p(0,2)
lbanned(p(i,2), 2) = .true.
end do
ma = 1
if(p(0,2) >= 2) ma = 2
mi = turn2(ma)
bant = 1
if(sp == 3) then
!move MA
if(ma == 2) bant = 2
puti = p(1,mi)
hfix = h(1,ma)
p1 = p(1,ma)
p2 = p(2,ma)
if(.not. bannedOrb(puti, mi)) then
tmp_row = 0d0
do putj=1, hfix-1
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)
tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do
do putj=hfix+1, mo_num
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)
tmp_row(1:N_states,putj) = tmp_row(1:N_states,putj) + hij * coefs(1:N_states)
end do
if(ma == 1) then
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
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
!MOVE MI
pfix = p(1,mi)
tmp_row = 0d0
tmp_row2 = 0d0
do puti=1,mo_num
if(lbanned(puti,mi)) cycle
!p1 fixed
putj = p1
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)
tmp_row(:,puti) = tmp_row(:,puti) + hij * coefs(:)
end if
putj = p2
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)
tmp_row2(:,puti) = tmp_row2(:,puti) + hij * coefs(:)
end if
end do
if(mi == 1) then
mat(:,:,p1) = mat(:,:,p1) + tmp_row(:,:)
mat(:,:,p2) = mat(:,:,p2) + tmp_row2(:,:)
else
mat(:,p1,:) = mat(:,p1,:) + tmp_row(:,:)
mat(:,p2,:) = mat(:,p2,:) + tmp_row2(:,:)
end if
else
if(p(0,ma) == 3) then
do i=1,3
hfix = h(1,ma)
puti = p(i, ma)
p1 = p(turn3(1,i), ma)
p2 = p(turn3(2,i), ma)
tmp_row = 0d0
do putj=1,hfix-1
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)
tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
end do
do putj=hfix+1,mo_num
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)
tmp_row(:,putj) = tmp_row(:,putj) + hij * coefs(:)
end do
mat(:, :puti-1, puti) = mat(:, :puti-1, puti) + tmp_row(:,:puti-1)