LCPQ
/
qp2
mirror of https://github.com/QuantumPackage/qp2.git
9
1
Fork 0
Code Releases Activity

Merge branch 'dev-stable' into dev-stable

This commit is contained in:
AbdAmmar 2024-05-01 18:54:03 +02:00 committed by GitHub
parent 22e3e63866 2c899e6dd7
commit 88e0b5fc74
No known key found for this signature in database
GPG Key ID: B5690EEEBB952194
21 changed files with 860 additions and 228 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

68
plugins/local/bi_ort_ints/semi_num_ints_mo.irp.f
View File

@ -107,8 +107,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
integer :: i, j, ipoint
double precision :: wall0, wall1
print *, ' providing int2_grad1_u12_ao_transp ...'
call wall_time(wall0)
!print *, ' providing int2_grad1_u12_ao_transp ...'
!call wall_time(wall0)
if(test_cycle_tc) then
@ -142,15 +142,15 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
endif
call wall_time(wall1)
print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0
call print_memory_usage()
!call wall_time(wall1)
!print *, ' wall time for int2_grad1_u12_ao_transp (min) = ', (wall1 - wall0) / 60.d0
!call print_memory_usage()
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)]
BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3, n_points_final_grid)]
implicit none
integer :: ipoint
@ -159,7 +159,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num,
PROVIDE mo_l_coef mo_r_coef
PROVIDE int2_grad1_u12_ao_transp
!print *, ' providing int2_grad1_u12_bimo_transp'
!print *, ' providing int2_grad1_u12_bimo_transp ...'
!call wall_time(wall0)
!$OMP PARALLEL &
@ -167,33 +167,35 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num,
!$OMP PRIVATE (ipoint) &
!$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp)
!$OMP DO SCHEDULE (dynamic)
do ipoint = 1, n_points_final_grid
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
enddo
do ipoint = 1, n_points_final_grid
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,1,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,1,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,2,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,2,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
call ao_to_mo_bi_ortho( int2_grad1_u12_ao_transp (1,1,3,ipoint), size(int2_grad1_u12_ao_transp , 1) &
, int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
enddo
!$OMP END DO
!$OMP END PARALLEL
!FREE int2_grad1_u12_ao_transp
!call wall_time(wall1)
!print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0
!print *, ' wall time for int2_grad1_u12_bimo_transp (min) =', (wall1 - wall0) / 60.d0
!call print_memory_usage()
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)]
BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3, mo_num, mo_num)]
implicit none
integer :: i, j, ipoint
double precision :: wall0, wall1
!call wall_time(wall0)
!print *, ' Providing int2_grad1_u12_bimo_t ...'
!print *, ' providing int2_grad1_u12_bimo_t ...'
PROVIDE mo_l_coef mo_r_coef
PROVIDE int2_grad1_u12_bimo_transp
@ -211,17 +213,21 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_bimo_t, (n_points_final_grid,
FREE int2_grad1_u12_bimo_transp
!call wall_time(wall1)
!print *, ' wall time for int2_grad1_u12_bimo_t,', wall1 - wall0
!print *, ' wall time for int2_grad1_u12_bimo_t (min) =', (wall1 - wall0) / 60.d0
!call print_memory_usage()
END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_num)]
BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3, ao_num, ao_num)]
implicit none
integer :: i, j, ipoint
integer :: i, j, ipoint
double precision :: wall0, wall1
!call wall_time(wall0)
!print *, ' providing int2_grad1_u12_ao_t ...'
PROVIDE int2_grad1_u12_ao
@ -235,6 +241,10 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3,
enddo
enddo
!call wall_time(wall1)
!print *, ' wall time for int2_grad1_u12_ao_t (min) =', (wall1 - wall0) / 60.d0
!call print_memory_usage()
END_PROVIDER
! ---
@ -275,8 +285,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid,
double precision :: xyz
double precision :: wall0, wall1
print*, ' providing x_W_ki_bi_ortho_erf_rk ...'
call wall_time(wall0)
!print*, ' providing x_W_ki_bi_ortho_erf_rk ...'
!call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -300,8 +310,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid,
! FREE mo_v_ki_bi_ortho_erf_rk_cst_mu_transp
! FREE mo_x_v_ki_bi_ortho_erf_rk_cst_mu_transp
call wall_time(wall1)
print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0
!call wall_time(wall1)
!print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0
END_PROVIDER
@ -323,8 +333,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_
double precision :: xyz
double precision :: wall0, wall1
print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...'
call wall_time(wall0)
!print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...'
!call wall_time(wall0)
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
@ -343,8 +353,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_
!$OMP END DO
!$OMP END PARALLEL
call wall_time(wall1)
print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0
!call wall_time(wall1)
!print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0
END_PROVIDER

2
plugins/local/bi_ort_ints/three_body_ints_bi_ort.irp.f
View File

@ -168,7 +168,7 @@ subroutine give_integrals_3_body_bi_ort(n, l, k, m, j, i, integral)
integral = integral + tmp * final_weight_at_r_vector(ipoint)
enddo
end subroutine give_integrals_3_body_bi_ort
end
! ---

120
plugins/local/bi_ort_ints/total_twoe_pot.irp.f
View File

@ -16,10 +16,10 @@ double precision function bi_ortho_mo_ints(l, k, j, i)
integer :: m, n, p, q
bi_ortho_mo_ints = 0.d0
do m = 1, ao_num
do p = 1, ao_num
do n = 1, ao_num
do q = 1, ao_num
do p = 1, ao_num
do m = 1, ao_num
do q = 1, ao_num
do n = 1, ao_num
! p1h1p2h2 l1 l2 r1 r2
bi_ortho_mo_ints += ao_two_e_tc_tot(n,q,m,p) * mo_l_coef(m,l) * mo_l_coef(n,k) * mo_r_coef(p,j) * mo_r_coef(q,i)
enddo
@ -27,7 +27,7 @@ double precision function bi_ortho_mo_ints(l, k, j, i)
enddo
enddo
end function bi_ortho_mo_ints
end
! ---
@ -40,38 +40,106 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num,
END_DOC
implicit none
integer :: i, j, k, l, m, n, p, q
integer :: i, j, k, l, m, n, p, q, s, r
double precision :: t1, t2, tt1, tt2
double precision, allocatable :: a1(:,:,:,:), a2(:,:,:,:)
double precision, allocatable :: a_jkp(:,:,:), a_kpq(:,:,:), ao_two_e_tc_tot_tmp(:,:,:)
print *, ' PROVIDING mo_bi_ortho_tc_two_e_chemist ...'
call wall_time(t1)
call print_memory_usage()
PROVIDE mo_r_coef mo_l_coef
allocate(a2(ao_num,ao_num,ao_num,mo_num))
if(ao_to_mo_tc_n3) then
call dgemm( 'T', 'N', ao_num*ao_num*ao_num, mo_num, ao_num, 1.d0 &
, ao_two_e_tc_tot(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num &
, 0.d0 , a2(1,1,1,1), ao_num*ao_num*ao_num)
print*, ' memory scale of TC ao -> mo: O(N3) '
allocate(a1(ao_num,ao_num,mo_num,mo_num))
if(.not.read_tc_integ) then
stop 'read_tc_integ needs to be set to true'
endif
call dgemm( 'T', 'N', ao_num*ao_num*mo_num, mo_num, ao_num, 1.d0 &
, a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num &
, 0.d0, a1(1,1,1,1), ao_num*ao_num*mo_num)
allocate(a_jkp(ao_num,ao_num,mo_num))
allocate(a_kpq(ao_num,mo_num,mo_num))
allocate(ao_two_e_tc_tot_tmp(ao_num,ao_num,ao_num))
deallocate(a2)
allocate(a2(ao_num,mo_num,mo_num,mo_num))
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="read")
call dgemm( 'T', 'N', ao_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
, a1(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num &
, 0.d0, a2(1,1,1,1), ao_num*mo_num*mo_num)
call wall_time(tt1)
deallocate(a1)
mo_bi_ortho_tc_two_e_chemist(:,:,:,:) = 0.d0
do l = 1, ao_num
read(11) ao_two_e_tc_tot_tmp(:,:,:)
call dgemm( 'T', 'N', mo_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
, a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num &
, 0.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,1), mo_num*mo_num*mo_num)
do s = 1, mo_num
deallocate(a2)
call dgemm( 'T', 'N', ao_num*ao_num, mo_num, ao_num, 1.d0 &
, ao_two_e_tc_tot_tmp(1,1,1), ao_num, mo_l_coef(1,1), ao_num &
, 0.d0, a_jkp(1,1,1), ao_num*ao_num)
call dgemm( 'T', 'N', ao_num*mo_num, mo_num, ao_num, 1.d0 &
, a_jkp(1,1,1), ao_num, mo_r_coef(1,1), ao_num &
, 0.d0, a_kpq(1,1,1), ao_num*mo_num)
call dgemm( 'T', 'N', mo_num*mo_num, mo_num, ao_num, mo_r_coef(l,s) &
, a_kpq(1,1,1), ao_num, mo_l_coef(1,1), ao_num &
, 1.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,s), mo_num*mo_num)
enddo ! s
if(l == 2) then
call wall_time(tt2)
print*, ' 1 / mo_num done in (min)', (tt2-tt1)/60.d0
print*, ' estimated time required (min)', dble(mo_num-1)*(tt2-tt1)/60.d0
elseif(l == 11) then
call wall_time(tt2)
print*, ' 10 / mo_num done in (min)', (tt2-tt1)/60.d0
print*, ' estimated time required (min)', dble(mo_num-10)*(tt2-tt1)/(60.d0*10.d0)
elseif(l == 101) then
call wall_time(tt2)
print*, ' 100 / mo_num done in (min)', (tt2-tt1)/60.d0
print*, ' estimated time required (min)', dble(mo_num-100)*(tt2-tt1)/(60.d0*100.d0)
endif
enddo ! l
close(11)
deallocate(a_jkp, a_kpq, ao_two_e_tc_tot_tmp)
else
print*, ' memory scale of TC ao -> mo: O(N4) '
allocate(a2(ao_num,ao_num,ao_num,mo_num))
call dgemm( 'T', 'N', ao_num*ao_num*ao_num, mo_num, ao_num, 1.d0 &
, ao_two_e_tc_tot(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num &
, 0.d0, a2(1,1,1,1), ao_num*ao_num*ao_num)
FREE ao_two_e_tc_tot
allocate(a1(ao_num,ao_num,mo_num,mo_num))
call dgemm( 'T', 'N', ao_num*ao_num*mo_num, mo_num, ao_num, 1.d0 &
, a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num &
, 0.d0, a1(1,1,1,1), ao_num*ao_num*mo_num)
deallocate(a2)
allocate(a2(ao_num,mo_num,mo_num,mo_num))
call dgemm( 'T', 'N', ao_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
, a1(1,1,1,1), ao_num, mo_l_coef(1,1), ao_num &
, 0.d0, a2(1,1,1,1), ao_num*mo_num*mo_num)
deallocate(a1)
call dgemm( 'T', 'N', mo_num*mo_num*mo_num, mo_num, ao_num, 1.d0 &
, a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num &
, 0.d0, mo_bi_ortho_tc_two_e_chemist(1,1,1,1), mo_num*mo_num*mo_num)
deallocate(a2)
endif
!allocate(a1(mo_num,ao_num,ao_num,ao_num))
!a1 = 0.d0
@ -135,6 +203,10 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num,
!enddo
!deallocate(a1)
call wall_time(t2)
print *, ' WALL TIME for PROVIDING mo_bi_ortho_tc_two_e_chemist (min)', (t2-t1)/60.d0
call print_memory_usage()
END_PROVIDER
! ---

8
plugins/local/bi_ortho_mos/overlap.irp.f
View File

@ -56,10 +56,10 @@
print*,'Average trace of overlap_bi_ortho is different from 1 by ', dabs(accu_d-1.d0)
print*,'And bi orthogonality is off by an average of ',accu_nd
print*,'****************'
print*,'Overlap matrix betwee mo_l_coef and mo_r_coef '
do i = 1, mo_num
write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:)
enddo
!print*,'Overlap matrix betwee mo_l_coef and mo_r_coef '
!do i = 1, mo_num
! write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:)
!enddo
endif
print*,'Average trace of overlap_bi_ortho (should be 1.)'
print*,'accu_d = ',accu_d

175
plugins/local/non_h_ints_mu/tc_integ.irp.f
View File

@ -44,14 +44,92 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
elseif(tc_integ_type .eq. "numeric") then
print *, ' Numerical integration over r1 and r2 will be performed'
! TODO combine 1shot & int2_grad1_u12_ao_num
PROVIDE int2_grad1_u12_ao_num
int2_grad1_u12_ao = int2_grad1_u12_ao_num
if(tc_save_mem) then
!PROVIDE int2_grad1_u12_ao_num_1shot
!int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot
integer :: n_blocks, n_rest, n_pass
integer :: i_blocks, i_rest, i_pass, ii
double precision :: mem, n_double
double precision, allocatable :: tmp(:,:,:), xx(:)
double precision, allocatable :: tmp_grad1_u12(:,:,:)
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num), xx(n_points_extra_final_grid))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j, i, jpoint) &
!$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
!$OMP DO COLLAPSE(2)
do j = 1, ao_num
do i = 1, ao_num
do jpoint = 1, n_points_extra_final_grid
tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call total_memory(mem)
mem = max(1.d0, qp_max_mem - mem)
n_double = mem * 1.d8
n_blocks = int(min(n_double / (n_points_extra_final_grid * 4.d0), 1.d0*n_points_final_grid))
n_rest = int(mod(n_points_final_grid, n_blocks))
n_pass = int((n_points_final_grid - n_rest) / n_blocks)
call write_int(6, n_pass, 'Number of passes')
call write_int(6, n_blocks, 'Size of the blocks')
call write_int(6, n_rest, 'Size of the last block')
allocate(tmp_grad1_u12(n_points_extra_final_grid,n_blocks,3))
do i_pass = 1, n_pass
ii = (i_pass-1)*n_blocks + 1
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_blocks, ipoint) &
!$OMP SHARED (n_blocks, n_points_extra_final_grid, ii, final_grid_points, xx, tmp_grad1_u12)
!$OMP DO
do i_blocks = 1, n_blocks
ipoint = ii - 1 + i_blocks ! r1
call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, tmp_grad1_u12(1,i_blocks,1), tmp_grad1_u12(1,i_blocks,2), tmp_grad1_u12(1,i_blocks,3), xx(1))
enddo
!$OMP END DO
!$OMP END PARALLEL
do m = 1, 3
call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, 1.d0 &
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid &
, 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num)
enddo
enddo
deallocate(tmp_grad1_u12)
if(n_rest .gt. 0) then
allocate(tmp_grad1_u12(n_points_extra_final_grid,n_rest,3))
ii = n_pass*n_blocks + 1
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_rest, ipoint) &
!$OMP SHARED (n_rest, n_points_extra_final_grid, ii, final_grid_points, xx, tmp_grad1_u12)
!$OMP DO
do i_rest = 1, n_rest
ipoint = ii - 1 + i_rest ! r1
call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, tmp_grad1_u12(1,i_rest,1), tmp_grad1_u12(1,i_rest,2), tmp_grad1_u12(1,i_rest,3), xx(1))
enddo
!$OMP END DO
!$OMP END PARALLEL
do m = 1, 3
call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, 1.d0 &
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12(1,1,m), n_points_extra_final_grid &
, 0.d0, int2_grad1_u12_ao(1,1,ii,m), ao_num*ao_num)
enddo
deallocate(tmp_grad1_u12)
endif
deallocate(tmp,xx)
else
! TODO combine 1shot & int2_grad1_u12_ao_num
PROVIDE int2_grad1_u12_ao_num
int2_grad1_u12_ao = int2_grad1_u12_ao_num
!PROVIDE int2_grad1_u12_ao_num_1shot
!int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot
endif
elseif(tc_integ_type .eq. "semi-analytic") then
@ -177,13 +255,88 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
print *, ' Numerical integration over r1 and r2 will be performed'
! TODO combine 1shot & int2_grad1_u12_square_ao_num
if(tc_save_mem) then
PROVIDE int2_grad1_u12_square_ao_num
int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num
integer :: n_blocks, n_rest, n_pass
integer :: i_blocks, i_rest, i_pass, ii
double precision :: mem, n_double
double precision, allocatable :: tmp(:,:,:), xx(:,:,:)
double precision, allocatable :: tmp_grad1_u12_squared(:,:)
!PROVIDE int2_grad1_u12_square_ao_num_1shot
!int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num_1shot
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
allocate(tmp(n_points_extra_final_grid,ao_num,ao_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (j, i, jpoint) &
!$OMP SHARED (tmp, ao_num, n_points_extra_final_grid, final_weight_at_r_vector_extra, aos_in_r_array_extra_transp)
!$OMP DO COLLAPSE(2)
do j = 1, ao_num
do i = 1, ao_num
do jpoint = 1, n_points_extra_final_grid
tmp(jpoint,i,j) = final_weight_at_r_vector_extra(jpoint) * aos_in_r_array_extra_transp(jpoint,i) * aos_in_r_array_extra_transp(jpoint,j)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call total_memory(mem)
mem = max(1.d0, qp_max_mem - mem)
n_double = mem * 1.d8
n_blocks = int(min(n_double / (n_points_extra_final_grid * 4.d0), 1.d0*n_points_final_grid))
n_rest = int(mod(n_points_final_grid, n_blocks))
n_pass = int((n_points_final_grid - n_rest) / n_blocks)
call write_int(6, n_pass, 'Number of passes')
call write_int(6, n_blocks, 'Size of the blocks')
call write_int(6, n_rest, 'Size of the last block')
allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_blocks), xx(n_points_extra_final_grid,n_blocks,3))
do i_pass = 1, n_pass
ii = (i_pass-1)*n_blocks + 1
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_blocks, ipoint) &
!$OMP SHARED (n_blocks, n_points_extra_final_grid, ii, xx, final_grid_points, tmp_grad1_u12_squared)
!$OMP DO
do i_blocks = 1, n_blocks
ipoint = ii - 1 + i_blocks ! r1
call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, xx(1,i_blocks,1), xx(1,i_blocks,2), xx(1,i_blocks,3), tmp_grad1_u12_squared(1,i_blocks))
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( "T", "N", ao_num*ao_num, n_blocks, n_points_extra_final_grid, -0.5d0 &
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid &
, 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num)
enddo
deallocate(tmp_grad1_u12_squared, xx)
if(n_rest .gt. 0) then
ii = n_pass*n_blocks + 1
allocate(tmp_grad1_u12_squared(n_points_extra_final_grid,n_rest), xx(n_points_extra_final_grid,n_rest,3))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i_rest, ipoint) &
!$OMP SHARED (n_rest, n_points_extra_final_grid, ii, xx, final_grid_points, tmp_grad1_u12_squared)
!$OMP DO
do i_rest = 1, n_rest
ipoint = ii - 1 + i_rest ! r1
call get_grad1_u12_withsq_r1_seq(ipoint, n_points_extra_final_grid, xx(1,i_rest,1), xx(1,i_rest,2), xx(1,i_rest,3), tmp_grad1_u12_squared(1,i_rest))
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( "T", "N", ao_num*ao_num, n_rest, n_points_extra_final_grid, -0.5d0 &
, tmp(1,1,1), n_points_extra_final_grid, tmp_grad1_u12_squared(1,1), n_points_extra_final_grid &
, 0.d0, int2_grad1_u12_square_ao(1,1,ii), ao_num*ao_num)
deallocate(tmp_grad1_u12_squared, xx)
endif
deallocate(tmp)
else
! TODO combine 1shot & int2_grad1_u12_square_ao_num
PROVIDE int2_grad1_u12_square_ao_num
int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num
!PROVIDE int2_grad1_u12_square_ao_num_1shot
!int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num_1shot
endif
elseif(tc_integ_type .eq. "semi-analytic") then

2
plugins/local/non_h_ints_mu/tc_integ_num.irp.f
View File

@ -131,7 +131,7 @@
deallocate(tmp)
call wall_time(time1)
print*, ' wall time for int2_grad1_u12_ao_num & int2_grad1_u12_square_ao_num =', time1-time0
print*, ' wall time for int2_grad1_u12_ao_num & int2_grad1_u12_square_ao_num = (min)', (time1-time0) / 60.d0
call print_memory_usage()
END_PROVIDER

235
plugins/local/non_h_ints_mu/total_tc_int.irp.f
View File

@ -33,8 +33,10 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
double precision :: weight1, ao_k_r, ao_i_r
double precision :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq
double precision :: time0, time1
double precision, allocatable :: b_mat(:,:,:,:), c_mat(:,:,:)
double precision, allocatable :: c_mat(:,:,:)
logical, external :: ao_two_e_integral_zero
double precision, external :: get_ao_two_e_integral
double precision, external :: ao_two_e_integral
PROVIDe tc_integ_type
PROVIDE env_type
@ -53,7 +55,9 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
print*, ' Reading ao_two_e_tc_tot from ', trim(ezfio_filename) // '/work/ao_two_e_tc_tot'
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="read")
read(11) ao_two_e_tc_tot
do i = 1, ao_num
read(11) ao_two_e_tc_tot(:,:,:,i)
enddo
close(11)
else
@ -65,27 +69,59 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
PROVIDE int2_grad1_u12_square_ao
allocate(c_mat(n_points_final_grid,ao_num,ao_num))
if(tc_save_mem_loops) then
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint) &
!$OMP SHARED (aos_in_r_array_transp, c_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
!$OMP DO SCHEDULE (static)
do i = 1, ao_num
do k = 1, ao_num
do ipoint = 1, n_points_final_grid
c_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k)
print*, ' LOOPS are used to evaluate Hermitian part of ao_two_e_tc_tot ...'
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, ipoint, ao_i_r, ao_k_r, weight1) &
!$OMP SHARED (ao_num, n_points_final_grid, ao_two_e_tc_tot, &
!$OMP aos_in_r_array_transp, final_weight_at_r_vector, int2_grad1_u12_square_ao)
!$OMP DO COLLAPSE(4)
do i = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
do j = 1, ao_num
ao_two_e_tc_tot(j,l,k,i) = 0.d0
do ipoint = 1, n_points_final_grid
weight1 = final_weight_at_r_vector(ipoint)
ao_i_r = aos_in_r_array_transp(ipoint,i)
ao_k_r = aos_in_r_array_transp(ipoint,k)
ao_two_e_tc_tot(j,l,k,i) = ao_two_e_tc_tot(j,l,k,i) + int2_grad1_u12_square_ao(j,l,ipoint) * weight1 * ao_i_r * ao_k_r
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!$OMP END DO
!$OMP END PARALLEL
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, c_mat(1,1,1), n_points_final_grid &
, 0.d0, ao_two_e_tc_tot, ao_num*ao_num)
else
print*, ' DGEMM are used to evaluate Hermitian part of ao_two_e_tc_tot ...'
allocate(c_mat(n_points_final_grid,ao_num,ao_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint) &
!$OMP SHARED (aos_in_r_array_transp, c_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
!$OMP DO SCHEDULE (static)
do i = 1, ao_num
do k = 1, ao_num
do ipoint = 1, n_points_final_grid
c_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
, int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, c_mat(1,1,1), n_points_final_grid &
, 0.d0, ao_two_e_tc_tot, ao_num*ao_num)
deallocate(c_mat)
endif
FREE int2_grad1_u12_square_ao
if( (tc_integ_type .eq. "semi-analytic") .and. &
@ -96,6 +132,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
! an additional term is added here directly instead of
! being added in int2_grad1_u12_square_ao for performance
allocate(c_mat(n_points_final_grid,ao_num,ao_num))
PROVIDE int2_u2_env2
!$OMP PARALLEL &
@ -127,10 +164,13 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
, int2_u2_env2(1,1,1), ao_num*ao_num, c_mat(1,1,1), n_points_final_grid &
, 1.d0, ao_two_e_tc_tot(1,1,1,1), ao_num*ao_num)
deallocate(c_mat)
FREE int2_u2_env2
endif ! use_ipp
deallocate(c_mat)
call wall_time(time1)
print*, ' done with Hermitian part after (min) ', (time1 - time0) / 60.d0
call print_memory_usage()
! ---
@ -138,39 +178,71 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
PROVIDE int2_grad1_u12_ao
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3))
if(tc_save_mem_loops) then
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) &
!$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, &
!$OMP ao_num, n_points_final_grid, final_weight_at_r_vector)
!$OMP DO SCHEDULE (static)
do i = 1, ao_num
do k = 1, ao_num
do ipoint = 1, n_points_final_grid
print*, ' LOOPS are used to evaluate non-Hermitian part of ao_two_e_tc_tot ...'
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
ao_i_r = aos_in_r_array_transp(ipoint,i)
ao_k_r = aos_in_r_array_transp(ipoint,k)
b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1))
b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2))
b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, j, k, l, ipoint, ao_i_r, ao_k_r, weight1) &
!$OMP SHARED (ao_num, n_points_final_grid, ao_two_e_tc_tot, &
!$OMP aos_in_r_array_transp, final_weight_at_r_vector, &
!$OMP int2_grad1_u12_ao, aos_grad_in_r_array_transp_bis)
!$OMP DO COLLAPSE(4)
do i = 1, ao_num
do k = 1, ao_num
do l = 1, ao_num
do j = 1, ao_num
do ipoint = 1, n_points_final_grid
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
ao_i_r = aos_in_r_array_transp(ipoint,i)
ao_k_r = aos_in_r_array_transp(ipoint,k)
ao_two_e_tc_tot(j,l,k,i) = ao_two_e_tc_tot(j,l,k,i) &
- weight1 * int2_grad1_u12_ao(j,l,ipoint,1) * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1)) &
- weight1 * int2_grad1_u12_ao(j,l,ipoint,2) * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2)) &
- weight1 * int2_grad1_u12_ao(j,l,ipoint,3) * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3))
enddo
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!$OMP END DO
!$OMP END PARALLEL
do m = 1, 3
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, -1.d0 &
, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
, 1.d0, ao_two_e_tc_tot(1,1,1,1), ao_num*ao_num)
enddo
deallocate(b_mat)
else
FREE int2_grad1_u12_ao
print*, ' DGEMM are used to evaluate non-Hermitian part of ao_two_e_tc_tot ...'
allocate(c_mat(n_points_final_grid,ao_num,ao_num))
do m = 1, 3
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r) &
!$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, c_mat, &
!$OMP ao_num, n_points_final_grid, final_weight_at_r_vector, m)
!$OMP DO SCHEDULE (static)
do i = 1, ao_num
do k = 1, ao_num
do ipoint = 1, n_points_final_grid
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
ao_i_r = aos_in_r_array_transp(ipoint,i)
ao_k_r = aos_in_r_array_transp(ipoint,k)
c_mat(ipoint,k,i) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,m) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,m))
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, -1.d0 &
, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num, c_mat(1,1,1), n_points_final_grid &
, 1.d0, ao_two_e_tc_tot(1,1,1,1), ao_num*ao_num)
enddo
deallocate(c_mat)
end if
if(tc_integ_type .eq. "semi-analytic") then
FREE int2_grad1_u2e_ao
@ -178,30 +250,67 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
endif ! var_tc
call wall_time(time1)
print*, ' done with non-Hermitian part after (min) ', (time1 - time0) / 60.d0
call print_memory_usage()
! ---
call sum_A_At(ao_two_e_tc_tot(1,1,1,1), ao_num*ao_num)
PROVIDE ao_integrals_map
! ---
logical :: integ_zero
double precision :: integ_val
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(ao_num, ao_two_e_tc_tot, ao_integrals_map) &
!$OMP PRIVATE(i, j, k, l)
!$OMP DO
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! < 1:i, 2:j | 1:k, 2:l >
ao_two_e_tc_tot(k,i,l,j) = ao_two_e_tc_tot(k,i,l,j) + get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
print*, ' adding ERI to ao_two_e_tc_tot ...'
if(tc_save_mem) then
print*, ' ao_integrals_map will not be used'
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP PRIVATE(i, j, k, l, integ_zero, integ_val) &
!$OMP SHARED(ao_num, ao_two_e_tc_tot)
!$OMP DO COLLAPSE(4)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
integ_zero = ao_two_e_integral_zero(i,j,k,l)
if(.not. integ_zero) then
! i,k : r1 j,l : r2
integ_val = ao_two_e_integral(i,k,j,l)
ao_two_e_tc_tot(k,i,l,j) = ao_two_e_tc_tot(k,i,l,j) + integ_val
endif
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!$OMP END DO
!$OMP END PARALLEL
else
print*, ' ao_integrals_map will be used'
PROVIDE ao_integrals_map
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(ao_num, ao_two_e_tc_tot, ao_integrals_map) &
!$OMP PRIVATE(i, j, k, l)
!$OMP DO COLLAPSE(4)
do j = 1, ao_num
do l = 1, ao_num
do i = 1, ao_num
do k = 1, ao_num
! < 1:i, 2:j | 1:k, 2:l >
ao_two_e_tc_tot(k,i,l,j) = ao_two_e_tc_tot(k,i,l,j) + get_ao_two_e_integral(i, j, k, l, ao_integrals_map)
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!call clear_ao_map()
FREE ao_integrals_map
endif
if(tc_integ_type .eq. "numeric") then
if((tc_integ_type .eq. "numeric") .and. (.not. tc_save_mem)) then
FREE int2_grad1_u12_ao_num int2_grad1_u12_square_ao_num
endif
@ -211,7 +320,9 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
print*, ' Saving ao_two_e_tc_tot in ', trim(ezfio_filename) // '/work/ao_two_e_tc_tot'
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="write")
call ezfio_set_work_empty(.False.)
write(11) ao_two_e_tc_tot
do i = 1, ao_num
write(11) ao_two_e_tc_tot(:,:,:,i)
enddo
close(11)
call ezfio_set_tc_keywords_io_tc_integ('Read')
endif

1
plugins/local/non_hermit_dav/lapack_diag_non_hermit.irp.f
View File

@ -2144,6 +2144,7 @@ subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0)
enddo
!print*,' accu_nd after = ', accu_nd
if(accu_nd .gt. 1d-12) then
print*, ' accu_nd =', accu_nd
print*, ' your strategy for degenerates orbitals failed !'
print*, m, 'deg on', i
stop

19
plugins/local/tc_keywords/EZFIO.cfg
View File

@ -285,3 +285,22 @@ type: logical
doc: If |true|, you minimize the angle between the left and right vectors associated to degenerate orbitals
interface: ezfio,provider,ocaml
default: False
[ao_to_mo_tc_n3]
type: logical
doc: If |true|, memory scale of TC ao -> mo: O(N3)
interface: ezfio,provider,ocaml
default: False
[tc_save_mem_loops]
type: logical
doc: If |true|, use loops to save memory TC
interface: ezfio,provider,ocaml
default: False
[tc_save_mem]
type: logical
doc: If |true|, more calc but less mem
interface: ezfio,provider,ocaml
default: False

121
plugins/local/tc_scf/fock_3e_bi_ortho_cs.irp.f
View File

@ -9,7 +9,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
double precision :: loc_1, loc_2, loc_3
double precision, allocatable :: Okappa(:), Jkappa(:,:)
double precision, allocatable :: tmp_omp_d1(:), tmp_omp_d2(:,:)
double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:)
double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:), tmp_22(:,:,:)
double precision, allocatable :: tmp_3(:,:,:), tmp_4(:,:,:)
PROVIDE mo_l_coef mo_r_coef
@ -63,17 +63,13 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
allocate(tmp_1(n_points_final_grid,4))
do ipoint = 1, n_points_final_grid
loc_1 = 2.d0 * Okappa(ipoint)
tmp_1(ipoint,1) = loc_1 * Jkappa(ipoint,1)
tmp_1(ipoint,2) = loc_1 * Jkappa(ipoint,2)
tmp_1(ipoint,3) = loc_1 * Jkappa(ipoint,3)
tmp_1(ipoint,4) = Okappa(ipoint)
enddo
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, i, j, loc_1, tmp_omp_d2) &
@ -112,58 +108,81 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
! ---
allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num))
if(tc_save_mem) then
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, a, b) &
!$OMP SHARED (n_points_final_grid, mo_num, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP tmp_2)
!$OMP DO COLLAPSE(2)
do a = 1, mo_num
do b = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp_2(ipoint,1,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a)
tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a)
tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a)
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, a, b, i) &
!$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, &
!$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t, &
!$OMP tmp_2)
!$OMP DO COLLAPSE(2)
do a = 1, mo_num
do b = 1, mo_num
tmp_2(:,4,b,a) = 0.d0
do i = 1, elec_beta_num
allocate(tmp_22(n_points_final_grid,4,mo_num))
do a = 1, mo_num
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, b, i) &
!$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, a, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP tmp_22)
!$OMP DO
do b = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp_2(ipoint,4,b,a) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) &
+ int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
+ int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) )
tmp_22(ipoint,1,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a)
tmp_22(ipoint,2,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a)
tmp_22(ipoint,3,b) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a)
enddo
tmp_22(:,4,b) = 0.d0
do i = 1, elec_beta_num
do ipoint = 1, n_points_final_grid
tmp_22(ipoint,4,b) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) &
+ int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
+ int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) )
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
call dgemv( 'T', 4*n_points_final_grid, mo_num, -2.d0 &
, tmp_22(1,1,1), size(tmp_22, 1) * size(tmp_22, 2) &
, tmp_1(1,1), 1 &
, 0.d0, fock_3e_uhf_mo_cs(1,a), 1)
enddo
deallocate(tmp_22)
else
allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num))
!$OMP PARALLEL &
!$OMP DEFAULT (NONE) &
!$OMP PRIVATE (ipoint, a, b, i) &
!$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, &
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
!$OMP tmp_2)
!$OMP DO COLLAPSE(2)
do a = 1, mo_num
do b = 1, mo_num
do ipoint = 1, n_points_final_grid
tmp_2(ipoint,1,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,1,b,a)
tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a)
tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a)
enddo
tmp_2(:,4,b,a) = 0.d0
do i = 1, elec_beta_num
do ipoint = 1, n_points_final_grid
tmp_2(ipoint,4,b,a) -= final_weight_at_r_vector(ipoint) * ( int2_grad1_u12_bimo_t(ipoint,1,b,i) * int2_grad1_u12_bimo_t(ipoint,1,i,a) &
+ int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
+ int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,a) )
enddo
enddo
enddo
enddo
enddo
!$OMP END DO
!$OMP END PARALLEL
!$OMP END DO
!$OMP END PARALLEL
call dgemv( 'T', 4*n_points_final_grid, mo_num*mo_num, -2.d0 &
, tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2) &
, tmp_1(1,1), 1 &
, 0.d0, fock_3e_uhf_mo_cs(1,1), 1)
deallocate(tmp_2)
! ---
endif
call dgemv( 'T', 4*n_points_final_grid, mo_num*mo_num, -2.d0 &
, tmp_2(1,1,1,1), size(tmp_2, 1) * size(tmp_2, 2) &
, tmp_1(1,1), 1 &
, 0.d0, fock_3e_uhf_mo_cs(1,1), 1)
deallocate(tmp_1, tmp_2)
deallocate(tmp_1)
! ---
@ -272,7 +291,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
! ---
!call wall_time(tf)
!print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti
!print *, ' total Wall time for fock_3e_uhf_mo_cs =', (tf - ti) / 60.d0
END_PROVIDER

2
plugins/local/tc_scf/fock_3e_bi_ortho_uhf.irp.f
View File

@ -32,7 +32,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
endif
!call wall_time(tf)
!print *, ' Wall time for fock_3e_uhf_mo_a =', tf - ti
!print *, ' Wall time for fock_3e_uhf_mo_a (min) =', (tf - ti)/60.d0
END_PROVIDER

2
plugins/local/tc_scf/fock_tc.irp.f
View File

@ -175,7 +175,7 @@ END_PROVIDER
! ---
BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
BEGIN_PROVIDER [double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num)]
BEGIN_DOC
! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis

2
plugins/local/tc_scf/minimize_tc_angles.irp.f
View File

@ -20,7 +20,7 @@ program minimize_tc_angles
! TODO
! check if rotations of orbitals affect the TC energy
! and refuse the step
call minimize_tc_orb_angles
call minimize_tc_orb_angles()
end

79
plugins/local/tc_scf/routines_rotates.irp.f
View File

@ -40,9 +40,6 @@ subroutine LTxSxR(n, m, L, S, R, C)
end subroutine LTxR
! ---
! ---
subroutine minimize_tc_orb_angles()
@ -103,7 +100,10 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:)
E_thr = 1d-04
PROVIDE TC_HF_energy
PROVIDE mo_r_coef mo_l_coef
E_thr = 1d-07
E_old = TC_HF_energy
allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num))
mo_r_coef_old = mo_r_coef
@ -111,7 +111,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
good_angles = .False.
allocate(mo_l_coef_good(ao_num, mo_num), mo_r_coef_good(ao_num,mo_num))
allocate(mo_l_coef_good(ao_num,mo_num), mo_r_coef_good(ao_num,mo_num))
print *, ' ***************************************'
print *, ' ***************************************'
@ -123,7 +123,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
mo_r_coef_good = mo_r_coef
mo_l_coef_good = mo_l_coef
allocate(mo_r_coef_new(ao_num, mo_num))
allocate(mo_r_coef_new(ao_num,mo_num))
mo_r_coef_new = mo_r_coef
do i = 1, mo_num
norm = 1.d0/dsqrt(overlap_mo_r(i,i))
@ -141,10 +141,11 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
call build_s_matrix(ao_num, mo_num, mo_r_coef_new, mo_r_coef_new, ao_overlap, s_mat)
! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
if(n_core_orb.ne.0)then
call give_degen_full_listcore(fock_diag, mo_num, list_core, n_core_orb, thr_deg, list_degen, n_degen_list)
call give_degen_full_listcore(fock_diag, mo_num, list_core, n_core_orb, thr_deg, list_degen, n_degen_list)
else
call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
call give_degen_full_list(fock_diag, mo_num, thr_deg, list_degen, n_degen_list)
endif
print *, ' fock_matrix_mo'
do i = 1, mo_num
print *, i, fock_diag(i), angle_left_right(i)
@ -156,50 +157,52 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
! n_degen = ilast - ifirst +1
n_degen = list_degen(i,0)
if(n_degen .ge. 1000)n_degen = 1 ! convention for core orbitals
if(n_degen .ge. 1000) n_degen = 1 ! convention for core orbitals
if(n_degen .eq. 1) cycle
print*, ' working on orbital', i
print*, ' multiplicity =', n_degen
allocate(stmp(n_degen,n_degen), smat2(n_degen,n_degen))
allocate(mo_r_coef_tmp(ao_num,n_degen), mo_l_coef_tmp(ao_num,n_degen), mo_l_coef_new(ao_num,n_degen))
allocate(T(n_degen,n_degen), Snew(n_degen,n_degen))
print*,'Right orbitals before'
do j = 1, n_degen
write(*,'(100(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j))
enddo
do j = 1, n_degen
write(*,'(1000(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j))
enddo
print*,'Left orbitals before'
do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_l_coef(1:ao_num,list_degen(i,j))
enddo
do j = 1, n_degen
write(*,'(1000(F16.10,X))') mo_l_coef(1:ao_num,list_degen(i,j))
enddo
if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then
integer :: i_list, j_list
i_list = list_degen(i,1)
j_list = list_degen(i,2)
print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(list_degen(i,2))
print*,'i_list = ',i_list
print*,'i_list = ',j_list
print*,'Swapping left/right orbitals'
call print_strong_overlap(i_list, j_list)
mo_r_coef_tmp(1:ao_num,1) = mo_r_coef_new(1:ao_num,i_list)
mo_r_coef_tmp(1:ao_num,2) = mo_l_coef(1:ao_num,i_list)
mo_l_coef_tmp(1:ao_num,1) = mo_l_coef(1:ao_num,j_list)
mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list)
integer :: i_list, j_list
i_list = list_degen(i,1)
j_list = list_degen(i,2)
print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(list_degen(i,2))
print*,'i_list = ',i_list
print*,'i_list = ',j_list
print*,'Swapping left/right orbitals'
call print_strong_overlap(i_list, j_list)
mo_r_coef_tmp(1:ao_num,1) = mo_r_coef_new(1:ao_num,i_list)
mo_r_coef_tmp(1:ao_num,2) = mo_l_coef(1:ao_num,i_list)
mo_l_coef_tmp(1:ao_num,1) = mo_l_coef(1:ao_num,j_list)
mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list)
else
do j = 1, n_degen
print*,'i_list = ',list_degen(i,j)
mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
enddo
do j = 1, n_degen
print*,'i_list = ',list_degen(i,j)
mo_r_coef_tmp(1:ao_num,j) = mo_r_coef_new(1:ao_num,list_degen(i,j))
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
enddo
endif
print*,'Right orbitals '
do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_r_coef_tmp(1:ao_num,j)
enddo
do j = 1, n_degen
write(*,'(1000(F16.10,X))') mo_r_coef_tmp(1:ao_num,j)
enddo
print*,'Left orbitals '
do j = 1, n_degen
write(*,'(100(F16.10,X))')mo_l_coef_tmp(1:ao_num,j)
enddo
do j = 1, n_degen
write(*,'(100(F16.10,X))') mo_l_coef_tmp(1:ao_num,j)
enddo
! Orthogonalization of right functions
print *, ' Orthogonalization of RIGHT functions'
print *, ' ------------------------------------'

58
plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f Normal file
View File

@ -0,0 +1,58 @@
! ---
program write_ao_2e_tc_integ
implicit none
PROVIDE j1e_type
PROVIDE j2e_type
print *, ' j1e_type = ', j1e_type
print *, ' j2e_type = ', j2e_type
my_grid_becke = .True.
PROVIDE tc_grid1_a tc_grid1_r
my_n_pt_r_grid = tc_grid1_r
my_n_pt_a_grid = tc_grid1_a
touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
call write_int(6, my_n_pt_r_grid, 'radial external grid over')
call write_int(6, my_n_pt_a_grid, 'angular external grid over')
if(tc_integ_type .eq. "numeric") then
my_extra_grid_becke = .True.
PROVIDE tc_grid2_a tc_grid2_r
my_n_pt_r_extra_grid = tc_grid2_r
my_n_pt_a_extra_grid = tc_grid2_a
touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
call write_int(6, my_n_pt_r_extra_grid, 'radial internal grid over')
call write_int(6, my_n_pt_a_extra_grid, 'angular internal grid over')
endif
call main()
end
! ---
subroutine main()
implicit none
PROVIDE io_tc_integ
print*, 'io_tc_integ = ', io_tc_integ
if(io_tc_integ .ne. "Write") then
print*, 'io_tc_integ != Write'
print*, io_tc_integ
stop
endif
PROVIDE ao_two_e_tc_tot
end
! ---

23
src/becke_numerical_grid/extra_grid_vector.irp.f
View File

@ -47,8 +47,12 @@ END_PROVIDER
END_DOC
implicit none
integer :: i,j,k,l,i_count
double precision :: r(3)
integer :: i, j, k, l, i_count
double precision :: r(3)
double precision :: wall0, wall1
call wall_time(wall0)
print *, ' Providing extra_final_grid_points ...'
i_count = 0
do j = 1, nucl_num
@ -66,10 +70,25 @@ END_PROVIDER
index_final_points_extra(2,i_count) = i
index_final_points_extra(3,i_count) = j
index_final_points_extra_reverse(k,i,j) = i_count
if(final_weight_at_r_vector_extra(i_count) .lt. 0.d0) then
print *, ' !!! WARNING !!!'
print *, ' negative weight !!!!'
print *, i_count, final_weight_at_r_vector_extra(i_count)
if(dabs(final_weight_at_r_vector_extra(i_count)) .lt. 1d-10) then
final_weight_at_r_vector_extra(i_count) = 0.d0
else
stop
endif
endif
enddo
enddo
enddo
call wall_time(wall1)
print *, ' wall time for extra_final_grid_points,', wall1 - wall0
call print_memory_usage()
END_PROVIDER

6
src/becke_numerical_grid/grid_becke_vector.irp.f
View File

@ -72,7 +72,11 @@ END_PROVIDER
print *, ' !!! WARNING !!!'
print *, ' negative weight !!!!'
print *, i_count, final_weight_at_r_vector(i_count)
stop
if(dabs(final_weight_at_r_vector(i_count)) .lt. 1d-10) then
final_weight_at_r_vector(i_count) = 0.d0
else
stop
endif
endif
enddo
enddo

39
src/mol_properties/multi_s_dipole_moment.irp.f
View File

@ -91,3 +91,42 @@ BEGIN_PROVIDER [double precision, multi_s_dipole_moment, (N_states, N_states)]
enddo
END_PROVIDER
! ---
BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment_eigenvec, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment_eigenvec, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment_eigenvec, (N_states, N_states)]
&BEGIN_PROVIDER [double precision, multi_s_x_dipole_moment_eigenval, (N_states)]
&BEGIN_PROVIDER [double precision, multi_s_y_dipole_moment_eigenval, (N_states)]
&BEGIN_PROVIDER [double precision, multi_s_z_dipole_moment_eigenval, (N_states)]
implicit none
double precision, allocatable :: eigval(:), eigvec(:,:), A(:,:)
PROVIDE multi_s_x_dipole_moment multi_s_y_dipole_moment multi_s_z_dipole_moment
allocate(A(N_states,N_states), eigvec(N_states,N_states), eigval(N_states))
A = multi_s_x_dipole_moment
call lapack_diag(eigval(1), eigvec(1,1), A(1,1), N_states, N_states)
multi_s_x_dipole_moment_eigenval = eigval
multi_s_x_dipole_moment_eigenvec = eigvec
A = multi_s_y_dipole_moment
call lapack_diag(eigval(1), eigvec(1,1), A(1,1), N_states, N_states)
multi_s_y_dipole_moment_eigenval = eigval
multi_s_y_dipole_moment_eigenvec = eigvec
A = multi_s_z_dipole_moment
call lapack_diag(eigval(1), eigvec(1,1), A(1,1), N_states, N_states)
multi_s_z_dipole_moment_eigenval = eigval
multi_s_z_dipole_moment_eigenvec = eigvec
deallocate(A, eigvec, eigval)
END_PROVIDER
! ---

97
src/tools/print_detweights.irp.f Normal file
View File

@ -0,0 +1,97 @@
program print_detweights
implicit none
read_wf = .True.
touch read_wf
call print_exc()
!call main()
end
! ---
subroutine main()
implicit none
integer :: i
integer :: degree
integer :: ios
integer, allocatable :: deg(:), ii(:), deg_sorted(:)
double precision, allocatable :: c(:)
PROVIDE N_int
PROVIDE N_det
PROVIDE psi_det
PROVIDe psi_coef
allocate(deg(N_det), ii(N_det), deg_sorted(N_det), c(N_det))
do i = 1, N_det
call debug_det(psi_det(1,1,i), N_int)
call get_excitation_degree(psi_det(1,1,i), psi_det(1,1,1), degree, N_int)
ii (i) = i
deg(i) = degree
c (i) = dabs(psi_coef(i,1))
enddo
call dsort(c, ii, N_det)
do i = 1, N_det
deg_sorted(i) = deg(ii(i))
print *, deg_sorted(i), c(i)
enddo
print *, ' saving psi'
! Writing output in binary format
open(unit=10, file="coef.bin", status="replace", action="write", iostat=ios, form="unformatted")
if(ios /= 0) then
print *, ' Error opening file!'
stop
endif
write(10) N_det
write(10) deg_sorted
write(10) c
close(10)
deallocate(deg, ii, deg_sorted, c)
end
! ---
subroutine print_exc()
implicit none
integer :: i
integer, allocatable :: deg(:)
PROVIDE N_int
PROVIDE N_det
PROVIDE psi_det
allocate(deg(N_det))
do i = 1, N_det
call get_excitation_degree(psi_det(1,1,1), psi_det(1,1,i), deg(i), N_int)
enddo
open(unit=10, file="exc.dat", action="write")
write(10,*) N_det
write(10,*) deg
close(10)
deallocate(deg)
end

2
src/utils/block_diag_degen.irp.f
View File

@ -191,7 +191,7 @@ subroutine give_degen_full_list(A, n, thr, list_degen, n_degen_list)
list_degen(n_degen_list,1) = i
icount = 1
do j = i+1, n
if(dabs(A(i)-A(j)).lt.thr.and.is_ok(j)) then
if(dabs(A(i)-A(j)).lt.thr .and. is_ok(j)) then
is_ok(j) = .False.
icount += 1
list_degen(n_degen_list,icount) = j

27
src/utils/constants.include.F
View File

@ -18,3 +18,30 @@ double precision, parameter :: c_4_3 = 4.d0/3.d0
double precision, parameter :: c_1_3 = 1.d0/3.d0
double precision, parameter :: sq_op5 = dsqrt(0.5d0)
double precision, parameter :: dlog_2pi = dlog(2.d0*dacos(-1.d0))
! physical constants and units conversion factors
double precision, parameter :: k_boltzman_si = 1.38066d-23 ! K k^-1
double precision, parameter :: k_boltzman_au = 3.1667d-6 ! Hartree k^-1
double precision, parameter :: k_boltzman_m1_au = 315795.26d0 ! Hartree^-1 k
double precision, parameter :: bohr_radius_si = 0.529177d-10 ! m
double precision, parameter :: bohr_radius_cm = 0.529177d-8 ! cm
double precision, parameter :: bohr_radius_angs = 0.529177d0 ! Angstrom
double precision, parameter :: electronmass_si = 9.10953d-31 ! Kg
double precision, parameter :: electronmass_uma = 5.4858d-4 ! uma
double precision, parameter :: electronvolt_si = 1.6021892d-19 ! J
double precision, parameter :: uma_si = 1.66057d-27 ! Kg
double precision, parameter :: debye_si = 3.33564d-30 ! coulomb meter
double precision, parameter :: debye_au = 0.393427228d0 ! e * Bohr
double precision, parameter :: angstrom_to_au = 1.889727d0 ! au
double precision, parameter :: au_to_ohmcmm1 = 46000.0d0 ! (ohm cm)^-1
double precision, parameter :: au_to_kb = 294210.0d0 ! kbar
double precision, parameter :: au_to_eV = 27.211652d0
double precision, parameter :: uma_to_au = 1822.89d0
double precision, parameter :: au_to_terahertz = 2.4189d-5
double precision, parameter :: au_to_sec = 2.4189d-17
double precision, parameter :: au_to_fsec = 2.4189d-2
double precision, parameter :: Wcm2 = 3.5d16
double precision, parameter :: amconv = 1.66042d-24/9.1095d-28*0.5d0 ! mass conversion: a.m.u to a.u. (ry)
double precision, parameter :: uakbar = 147105.d0 ! pressure conversion from ry/(a.u)^3 to k