mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-08-29 20:23:40 +02:00
Merge branch 'dev-stable' of https://github.com/QuantumPackage/qp2 into dev-stable
This commit is contained in:
commit
c150f2a2a2
@ -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
|
||||
|
||||
|
@ -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
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -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
|
||||
|
||||
! ---
|
||||
|
@ -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
|
||||
|
@ -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
|
||||
|
||||
|
@ -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
|
||||
|
@ -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
|
||||
|
@ -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
|
||||
|
@ -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
|
||||
|
||||
|
@ -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
|
||||
|
||||
|
@ -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
|
||||
|
||||
|
@ -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
|
||||
|
@ -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
|
||||
|
||||
|
@ -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
58
plugins/local/tc_scf/write_ao_2e_tc_integ.irp.f
Normal 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
|
||||
|
||||
! ---
|
||||
|
@ -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
|
||||
|
||||
|
||||
|
@ -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
|
||||
|
@ -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
97
src/tools/print_detweights.irp.f
Normal 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
|
||||
|
||||
|
||||
|
@ -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
|
||||
|
@ -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
|
||||
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user