mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-25 05:31:48 +02:00
Fast 4-idx transformation
This commit is contained in:
parent
0336738109
commit
a744bc30d4
@ -153,4 +153,3 @@ BEGIN_PROVIDER [real*8, bielecCI, (n_act_orb,n_act_orb,n_act_orb, mo_num)]
|
|||||||
!$OMP END PARALLEL DO
|
!$OMP END PARALLEL DO
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
@ -4,10 +4,8 @@ program casscf
|
|||||||
! TODO : Put the documentation of the program here
|
! TODO : Put the documentation of the program here
|
||||||
END_DOC
|
END_DOC
|
||||||
no_vvvv_integrals = .True.
|
no_vvvv_integrals = .True.
|
||||||
no_ivvv_integrals = .True.
|
|
||||||
no_vvv_integrals = .True.
|
|
||||||
pt2_max = 0.02
|
pt2_max = 0.02
|
||||||
SOFT_TOUCH no_vvvv_integrals no_vvv_integrals pt2_max
|
SOFT_TOUCH no_vvvv_integrals pt2_max
|
||||||
call run
|
call run
|
||||||
end
|
end
|
||||||
|
|
||||||
@ -41,8 +39,7 @@ subroutine run
|
|||||||
psi_det = psi_det_sorted
|
psi_det = psi_det_sorted
|
||||||
psi_coef = psi_coef_sorted
|
psi_coef = psi_coef_sorted
|
||||||
read_wf = .True.
|
read_wf = .True.
|
||||||
call map_deinit(mo_integrals_map)
|
call clear_mo_map
|
||||||
FREE mo_two_e_integrals_in_map mo_integrals_map
|
|
||||||
SOFT_TOUCH mo_coef N_det pt2_max psi_det psi_coef
|
SOFT_TOUCH mo_coef N_det pt2_max psi_det psi_coef
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
|
@ -11,24 +11,3 @@ interface: ezfio,provider,ocaml
|
|||||||
default: 1.e-15
|
default: 1.e-15
|
||||||
ezfio_name: threshold_mo
|
ezfio_name: threshold_mo
|
||||||
|
|
||||||
[no_vvvv_integrals]
|
|
||||||
type: logical
|
|
||||||
doc: If `True`, computes all integrals except for the integrals having 4 virtual indices
|
|
||||||
interface: ezfio,provider,ocaml
|
|
||||||
default: False
|
|
||||||
ezfio_name: no_vvvv_integrals
|
|
||||||
|
|
||||||
[no_ivvv_integrals]
|
|
||||||
type: logical
|
|
||||||
doc: Can be switched on only if `no_vvvv_integrals` is `True`, then does not compute the integrals with 3 virtual indices and 1 belonging to the core inactive active orbitals
|
|
||||||
interface: ezfio,provider,ocaml
|
|
||||||
default: False
|
|
||||||
ezfio_name: no_ivvv_integrals
|
|
||||||
|
|
||||||
[no_vvv_integrals]
|
|
||||||
type: logical
|
|
||||||
doc: Can be switched on only if `no_vvvv_integrals` is `True`, then does not compute the integrals with 3 virtual orbitals
|
|
||||||
interface: ezfio,provider,ocaml
|
|
||||||
default: False
|
|
||||||
ezfio_name: no_vvv_integrals
|
|
||||||
|
|
||||||
|
@ -1,3 +1,12 @@
|
|||||||
|
BEGIN_PROVIDER [ logical, no_vvvv_integrals ]
|
||||||
|
implicit none
|
||||||
|
BEGIN_DOC
|
||||||
|
! If `True`, computes all integrals except for the integrals having 3 or 4 virtual indices
|
||||||
|
END_DOC
|
||||||
|
|
||||||
|
no_vvvv_integrals = .False.
|
||||||
|
END_PROVIDER
|
||||||
|
|
||||||
BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ]
|
BEGIN_PROVIDER [ double precision, mo_coef_novirt, (ao_num,n_core_inact_act_orb) ]
|
||||||
implicit none
|
implicit none
|
||||||
BEGIN_DOC
|
BEGIN_DOC
|
||||||
@ -48,134 +57,66 @@ subroutine four_idx_novvvv
|
|||||||
! Retransform MO integrals for next CAS-SCF step
|
! Retransform MO integrals for next CAS-SCF step
|
||||||
END_DOC
|
END_DOC
|
||||||
integer :: i,j,k,l,n_integrals
|
integer :: i,j,k,l,n_integrals
|
||||||
double precision, allocatable :: f(:,:,:), d(:,:), T(:,:,:,:)
|
double precision, allocatable :: f(:,:,:), f2(:,:,:), d(:,:), T(:,:,:,:), T2(:,:,:,:)
|
||||||
double precision, external :: get_ao_two_e_integral
|
double precision, external :: get_ao_two_e_integral
|
||||||
integer(key_kind), allocatable :: idx(:)
|
integer(key_kind), allocatable :: idx(:)
|
||||||
real(integral_kind), allocatable :: values(:)
|
real(integral_kind), allocatable :: values(:)
|
||||||
|
|
||||||
integer :: p,q,r,s
|
integer :: p,q,r,s
|
||||||
double precision, allocatable :: ijij(:), ijji(:), jqjs(:,:,:), jqrj(:,:,:)
|
|
||||||
double precision :: c
|
double precision :: c
|
||||||
allocate (jqjs(mo_num,ao_num,ao_num), jqrj(mo_num,ao_num,ao_num))
|
allocate( T(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) , &
|
||||||
|
T2(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) )
|
||||||
allocate( T(n_core_inact_act_orb,n_core_inact_act_orb,ao_num,ao_num) )
|
|
||||||
|
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
!$OMP SHARED(mo_num,ao_num,T,n_core_inact_act_orb, mo_coef_transp, &
|
!$OMP SHARED(mo_num,ao_num,T,n_core_inact_act_orb, mo_coef_transp, &
|
||||||
!$OMP mo_integrals_threshold,mo_coef,mo_integrals_map, &
|
!$OMP mo_integrals_threshold,mo_coef,mo_integrals_map, &
|
||||||
!$OMP list_core_inact_act,jqjs,jqrj,ao_integrals_map) &
|
!$OMP list_core_inact_act,T2,ao_integrals_map) &
|
||||||
!$OMP PRIVATE(i,j,k,l,p,q,r,s,idx,values,n_integrals, &
|
!$OMP PRIVATE(i,j,k,l,p,q,r,s,idx,values,n_integrals, &
|
||||||
!$OMP f,d,c,ijij,ijji)
|
!$OMP f,f2,d,c)
|
||||||
allocate(f(ao_num,ao_num,ao_num), d(mo_num,mo_num), &
|
allocate(f(ao_num,ao_num,ao_num), f2(ao_num,ao_num,ao_num), d(mo_num,mo_num), &
|
||||||
idx(mo_num*mo_num), values(mo_num*mo_num) )
|
idx(mo_num*mo_num), values(mo_num*mo_num) )
|
||||||
allocate(ijij(mo_num), ijji(mo_num))
|
|
||||||
|
|
||||||
! <aa|vv>
|
! <aa|vv>
|
||||||
!$OMP DO
|
|
||||||
do s=1,ao_num
|
|
||||||
T(:,:,:,s) = 0.d0
|
|
||||||
jqjs(:,:,s) = 0.d0
|
|
||||||
jqrj(:,:,s) = 0.d0
|
|
||||||
enddo
|
|
||||||
!$OMP END DO
|
|
||||||
|
|
||||||
!$OMP DO
|
!$OMP DO
|
||||||
do s=1,ao_num
|
do s=1,ao_num
|
||||||
do r=1,ao_num
|
do r=1,ao_num
|
||||||
do q=1,ao_num
|
do q=1,ao_num
|
||||||
do p=1,r
|
do p=1,r
|
||||||
f(p,q,r) = get_ao_two_e_integral(p,q,r,s,ao_integrals_map)
|
f (p,q,r) = get_ao_two_e_integral(p,q,r,s,ao_integrals_map)
|
||||||
f(r,q,p) = f(p,q,r)
|
f (r,q,p) = f(p,q,r)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
! f(p,q,r) = <pq|rs>
|
|
||||||
|
|
||||||
do r=1,ao_num
|
|
||||||
call ao_to_mo_novirt(f(1,1,r),size(f,1),T(1,1,r,s),size(T,1))
|
|
||||||
enddo
|
|
||||||
! T(i,j,p,q) = <ij|rs>
|
|
||||||
|
|
||||||
! Diagonal
|
|
||||||
do r=1,ao_num
|
do r=1,ao_num
|
||||||
do q=1,ao_num
|
do q=1,ao_num
|
||||||
do p=1,ao_num
|
do p=1,ao_num
|
||||||
if (dabs(f(p,q,r)) >= mo_integrals_threshold) then
|
f2(p,q,r) = f(p,r,q)
|
||||||
do i=1,mo_num
|
|
||||||
jqjs(i,q,s) = jqjs(i,q,s) + mo_coef_transp(i,p) * f(p,q,r) * mo_coef_transp(i,r)
|
|
||||||
enddo
|
|
||||||
endif
|
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
! f (p,q,r) = <pq|rs>
|
||||||
enddo
|
! f2(p,q,r) = <pr|qs>
|
||||||
!$OMP END DO NOWAIT
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do s=1,ao_num
|
|
||||||
do r=1,ao_num
|
do r=1,ao_num
|
||||||
do q=1,r
|
call ao_to_mo_novirt(f (1,1,r),size(f ,1),T (1,1,r,s),size(T,1))
|
||||||
do p=1,ao_num
|
call ao_to_mo_novirt(f2(1,1,r),size(f2,1),T2(1,1,r,s),size(T,1))
|
||||||
f(p,q,r) = get_ao_two_e_integral(p,q,s,r,ao_integrals_map)
|
|
||||||
f(p,r,q) = f(p,q,r)
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
enddo
|
||||||
! f(p,q,r) = <pq|sr>
|
! T (i,j,p,q) = <ij|rs>
|
||||||
|
! T2(i,j,p,q) = <ir|js>
|
||||||
|
|
||||||
! Diagonal
|
|
||||||
do r=1,ao_num
|
|
||||||
do q=1,ao_num
|
|
||||||
do p=1,ao_num
|
|
||||||
if (dabs(f(p,q,r)) >= mo_integrals_threshold) then
|
|
||||||
do i=1,mo_num
|
|
||||||
jqrj(i,q,s) = jqrj(i,q,s) + mo_coef_transp(i,p) * f(p,q,r) * mo_coef_transp(i,r)
|
|
||||||
enddo
|
|
||||||
endif
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
|
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO NOWAIT
|
!$OMP END DO
|
||||||
|
|
||||||
!$OMP BARRIER
|
|
||||||
|
|
||||||
!$OMP DO
|
|
||||||
do i=1,mo_num
|
|
||||||
ijij(:) = 0.d0
|
|
||||||
ijji(:) = 0.d0
|
|
||||||
do s=1,ao_num
|
|
||||||
do q=1,ao_num
|
|
||||||
do j=1,mo_num
|
|
||||||
c = mo_coef_transp(j,q) * mo_coef_transp(j,s)
|
|
||||||
ijij(j) = ijij(j) + jqjs(i,q,s) * c
|
|
||||||
ijji(j) = ijji(j) + jqrj(i,q,s) * c
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
enddo
|
|
||||||
do j=1,mo_num
|
|
||||||
call two_e_integrals_index(i,j,i,j,idx(j))
|
|
||||||
values(j) = ijij(j)
|
|
||||||
enddo
|
|
||||||
do j=1,mo_num
|
|
||||||
call two_e_integrals_index(i,j,j,i,idx(mo_num+j))
|
|
||||||
values(mo_num+j) = ijji(j)
|
|
||||||
enddo
|
|
||||||
call map_append(mo_integrals_map, idx, values, 2*mo_num)
|
|
||||||
enddo
|
|
||||||
!$OMP END DO
|
|
||||||
|
|
||||||
!$OMP DO
|
!$OMP DO
|
||||||
do j=1,n_core_inact_act_orb
|
do j=1,n_core_inact_act_orb
|
||||||
do i=1,n_core_inact_act_orb
|
do i=1,n_core_inact_act_orb
|
||||||
do s=1,ao_num
|
do s=1,ao_num
|
||||||
do r=1,ao_num
|
do r=1,ao_num
|
||||||
f(r,s,1) = T(i,j,r,s)
|
f (r,s,1) = T (i,j,r,s)
|
||||||
|
f2(r,s,1) = T2(i,j,r,s)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
call ao_to_mo(f,size(f,1),d,size(d,1))
|
call ao_to_mo(f ,size(f ,1),d,size(d,1))
|
||||||
n_integrals = 0
|
n_integrals = 0
|
||||||
do l=1,mo_num
|
do l=1,mo_num
|
||||||
do k=1,mo_num
|
do k=1,mo_num
|
||||||
@ -185,14 +126,25 @@ subroutine four_idx_novvvv
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
call map_append(mo_integrals_map, idx, values, n_integrals)
|
call map_append(mo_integrals_map, idx, values, n_integrals)
|
||||||
|
|
||||||
|
call ao_to_mo(f2,size(f2,1),d,size(d,1))
|
||||||
|
n_integrals = 0
|
||||||
|
do l=1,mo_num
|
||||||
|
do k=1,mo_num
|
||||||
|
n_integrals+=1
|
||||||
|
call two_e_integrals_index(list_core_inact_act(i),k,list_core_inact_act(j),l,idx(n_integrals))
|
||||||
|
values(n_integrals) = d(k,l)
|
||||||
|
enddo
|
||||||
|
enddo
|
||||||
|
call map_append(mo_integrals_map, idx, values, n_integrals)
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
deallocate(f,d,ijij,ijji,idx,values)
|
deallocate(f,f2,d,idx,values)
|
||||||
|
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
deallocate(T,jqrj,jqjs)
|
deallocate(T,T2)
|
||||||
|
|
||||||
|
|
||||||
call map_sort(mo_integrals_map)
|
call map_sort(mo_integrals_map)
|
||||||
@ -200,3 +152,29 @@ subroutine four_idx_novvvv
|
|||||||
call map_shrink(mo_integrals_map,real(mo_integrals_threshold,integral_kind))
|
call map_shrink(mo_integrals_map,real(mo_integrals_threshold,integral_kind))
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
subroutine four_idx_novvvv2
|
||||||
|
use bitmasks
|
||||||
|
implicit none
|
||||||
|
integer :: i
|
||||||
|
integer(bit_kind) :: mask_ijkl(N_int,4)
|
||||||
|
|
||||||
|
print*, '<ix|ix>'
|
||||||
|
do i = 1,N_int
|
||||||
|
mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
|
||||||
|
mask_ijkl(i,2) = full_ijkl_bitmask_4(i,1)
|
||||||
|
mask_ijkl(i,3) = core_inact_act_bitmask_4(i,1)
|
||||||
|
mask_ijkl(i,4) = full_ijkl_bitmask_4(i,1)
|
||||||
|
enddo
|
||||||
|
call add_integrals_to_map(mask_ijkl)
|
||||||
|
|
||||||
|
print*, '<ii|vv>'
|
||||||
|
do i = 1,N_int
|
||||||
|
mask_ijkl(i,1) = core_inact_act_bitmask_4(i,1)
|
||||||
|
mask_ijkl(i,2) = core_inact_act_bitmask_4(i,1)
|
||||||
|
mask_ijkl(i,3) = virt_bitmask(i,1)
|
||||||
|
mask_ijkl(i,4) = virt_bitmask(i,1)
|
||||||
|
enddo
|
||||||
|
call add_integrals_to_map(mask_ijkl)
|
||||||
|
|
||||||
|
end
|
||||||
|
@ -120,10 +120,6 @@ subroutine add_integrals_to_map(mask_ijkl)
|
|||||||
call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int )
|
call bitstring_to_list( mask_ijkl(1,2), list_ijkl(1,2), n_j, N_int )
|
||||||
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
|
call bitstring_to_list( mask_ijkl(1,3), list_ijkl(1,3), n_k, N_int )
|
||||||
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
|
call bitstring_to_list( mask_ijkl(1,4), list_ijkl(1,4), n_l, N_int )
|
||||||
character*(2048) :: output(1)
|
|
||||||
print *, 'i'
|
|
||||||
call bitstring_to_str( output(1), mask_ijkl(1,1), N_int )
|
|
||||||
print *, trim(output(1))
|
|
||||||
j = 0
|
j = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
j += popcnt(mask_ijkl(i,1))
|
j += popcnt(mask_ijkl(i,1))
|
||||||
@ -132,9 +128,6 @@ subroutine add_integrals_to_map(mask_ijkl)
|
|||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
|
||||||
print*, 'j'
|
|
||||||
call bitstring_to_str( output(1), mask_ijkl(1,2), N_int )
|
|
||||||
print *, trim(output(1))
|
|
||||||
j = 0
|
j = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
j += popcnt(mask_ijkl(i,2))
|
j += popcnt(mask_ijkl(i,2))
|
||||||
@ -143,9 +136,6 @@ subroutine add_integrals_to_map(mask_ijkl)
|
|||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
|
||||||
print*, 'k'
|
|
||||||
call bitstring_to_str( output(1), mask_ijkl(1,3), N_int )
|
|
||||||
print *, trim(output(1))
|
|
||||||
j = 0
|
j = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
j += popcnt(mask_ijkl(i,3))
|
j += popcnt(mask_ijkl(i,3))
|
||||||
@ -154,9 +144,6 @@ subroutine add_integrals_to_map(mask_ijkl)
|
|||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
|
||||||
print*, 'l'
|
|
||||||
call bitstring_to_str( output(1), mask_ijkl(1,4), N_int )
|
|
||||||
print *, trim(output(1))
|
|
||||||
j = 0
|
j = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
j += popcnt(mask_ijkl(i,4))
|
j += popcnt(mask_ijkl(i,4))
|
||||||
@ -171,6 +158,7 @@ subroutine add_integrals_to_map(mask_ijkl)
|
|||||||
|
|
||||||
double precision :: accu_bis
|
double precision :: accu_bis
|
||||||
accu_bis = 0.d0
|
accu_bis = 0.d0
|
||||||
|
call wall_time(wall_1)
|
||||||
|
|
||||||
!$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, &
|
!$OMP PARALLEL PRIVATE(l1,k1,j1,i1,i2,i3,i4,i,j,k,l,c, ii1,kmax, &
|
||||||
!$OMP two_e_tmp_0_idx, two_e_tmp_0, two_e_tmp_1,two_e_tmp_2,two_e_tmp_3,&
|
!$OMP two_e_tmp_0_idx, two_e_tmp_0, two_e_tmp_1,two_e_tmp_2,two_e_tmp_3,&
|
||||||
@ -414,10 +402,6 @@ subroutine add_integrals_to_map_three_indices(mask_ijk)
|
|||||||
call bitstring_to_list( mask_ijk(1,1), list_ijkl(1,1), n_i, N_int )
|
call bitstring_to_list( mask_ijk(1,1), list_ijkl(1,1), n_i, N_int )
|
||||||
call bitstring_to_list( mask_ijk(1,2), list_ijkl(1,2), n_j, N_int )
|
call bitstring_to_list( mask_ijk(1,2), list_ijkl(1,2), n_j, N_int )
|
||||||
call bitstring_to_list( mask_ijk(1,3), list_ijkl(1,3), n_k, N_int )
|
call bitstring_to_list( mask_ijk(1,3), list_ijkl(1,3), n_k, N_int )
|
||||||
character*(2048) :: output(1)
|
|
||||||
print*, 'i'
|
|
||||||
call bitstring_to_str( output(1), mask_ijk(1,1), N_int )
|
|
||||||
print *, trim(output(1))
|
|
||||||
j = 0
|
j = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
j += popcnt(mask_ijk(i,1))
|
j += popcnt(mask_ijk(i,1))
|
||||||
@ -426,9 +410,6 @@ subroutine add_integrals_to_map_three_indices(mask_ijk)
|
|||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
|
||||||
print*, 'j'
|
|
||||||
call bitstring_to_str( output(1), mask_ijk(1,2), N_int )
|
|
||||||
print *, trim(output(1))
|
|
||||||
j = 0
|
j = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
j += popcnt(mask_ijk(i,2))
|
j += popcnt(mask_ijk(i,2))
|
||||||
@ -437,9 +418,6 @@ subroutine add_integrals_to_map_three_indices(mask_ijk)
|
|||||||
return
|
return
|
||||||
endif
|
endif
|
||||||
|
|
||||||
print*, 'k'
|
|
||||||
call bitstring_to_str( output(1), mask_ijk(1,3), N_int )
|
|
||||||
print *, trim(output(1))
|
|
||||||
j = 0
|
j = 0
|
||||||
do i = 1, N_int
|
do i = 1, N_int
|
||||||
j += popcnt(mask_ijk(i,3))
|
j += popcnt(mask_ijk(i,3))
|
||||||
|
Loading…
Reference in New Issue
Block a user