mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-11-19 04:22:32 +01:00
Compare commits
No commits in common. "0465a0f4397a53daa5a3a1c8374a5e34f5b61c67" and "22e3e63866c6cac1ee85b5f5a427b0e90bebd044" have entirely different histories.
0465a0f439
...
22e3e63866
2
external/irpf90
vendored
2
external/irpf90
vendored
@ -1 +1 @@
|
|||||||
Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6
|
Subproject commit beac615343f421bd6c0571a408ba389a6d5a32ac
|
||||||
@ -107,8 +107,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
|
|||||||
integer :: i, j, ipoint
|
integer :: i, j, ipoint
|
||||||
double precision :: wall0, wall1
|
double precision :: wall0, wall1
|
||||||
|
|
||||||
!print *, ' providing int2_grad1_u12_ao_transp ...'
|
print *, ' providing int2_grad1_u12_ao_transp ...'
|
||||||
!call wall_time(wall0)
|
call wall_time(wall0)
|
||||||
|
|
||||||
if(test_cycle_tc) then
|
if(test_cycle_tc) then
|
||||||
|
|
||||||
@ -142,15 +142,15 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_transp, (ao_num, ao_num, 3,
|
|||||||
|
|
||||||
endif
|
endif
|
||||||
|
|
||||||
!call wall_time(wall1)
|
call wall_time(wall1)
|
||||||
!print *, ' wall time for int2_grad1_u12_ao_transp (min) = ', (wall1 - wall0) / 60.d0
|
print *, ' wall time for int2_grad1_u12_ao_transp ', wall1 - wall0
|
||||||
!call print_memory_usage()
|
call print_memory_usage()
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
implicit none
|
||||||
integer :: ipoint
|
integer :: ipoint
|
||||||
@ -159,7 +159,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3
|
|||||||
PROVIDE mo_l_coef mo_r_coef
|
PROVIDE mo_l_coef mo_r_coef
|
||||||
PROVIDE int2_grad1_u12_ao_transp
|
PROVIDE int2_grad1_u12_ao_transp
|
||||||
|
|
||||||
!print *, ' providing int2_grad1_u12_bimo_transp ...'
|
!print *, ' providing int2_grad1_u12_bimo_transp'
|
||||||
!call wall_time(wall0)
|
!call wall_time(wall0)
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
@ -167,35 +167,33 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_transp, (mo_num, mo_num, 3
|
|||||||
!$OMP PRIVATE (ipoint) &
|
!$OMP PRIVATE (ipoint) &
|
||||||
!$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp)
|
!$OMP SHARED (n_points_final_grid,int2_grad1_u12_ao_transp,int2_grad1_u12_bimo_transp)
|
||||||
!$OMP DO SCHEDULE (dynamic)
|
!$OMP DO SCHEDULE (dynamic)
|
||||||
do ipoint = 1, n_points_final_grid
|
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) &
|
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) )
|
, 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) &
|
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) )
|
, 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) &
|
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) )
|
, int2_grad1_u12_bimo_transp(1,1,3,ipoint), size(int2_grad1_u12_bimo_transp, 1) )
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
!FREE int2_grad1_u12_ao_transp
|
|
||||||
|
|
||||||
!call wall_time(wall1)
|
!call wall_time(wall1)
|
||||||
!print *, ' wall time for int2_grad1_u12_bimo_transp (min) =', (wall1 - wall0) / 60.d0
|
!print *, ' Wall time for providing int2_grad1_u12_bimo_transp',wall1 - wall0
|
||||||
!call print_memory_usage()
|
!call print_memory_usage()
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
implicit none
|
||||||
integer :: i, j, ipoint
|
integer :: i, j, ipoint
|
||||||
double precision :: wall0, wall1
|
double precision :: wall0, wall1
|
||||||
|
|
||||||
!call wall_time(wall0)
|
!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 mo_l_coef mo_r_coef
|
||||||
PROVIDE int2_grad1_u12_bimo_transp
|
PROVIDE int2_grad1_u12_bimo_transp
|
||||||
@ -213,21 +211,17 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_bimo_t, (n_points_final_grid, 3
|
|||||||
FREE int2_grad1_u12_bimo_transp
|
FREE int2_grad1_u12_bimo_transp
|
||||||
|
|
||||||
!call wall_time(wall1)
|
!call wall_time(wall1)
|
||||||
!print *, ' wall time for int2_grad1_u12_bimo_t (min) =', (wall1 - wall0) / 60.d0
|
!print *, ' wall time for int2_grad1_u12_bimo_t,', wall1 - wall0
|
||||||
!call print_memory_usage()
|
!call print_memory_usage()
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
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
|
PROVIDE int2_grad1_u12_ao
|
||||||
|
|
||||||
@ -241,10 +235,6 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao_t, (n_points_final_grid, 3,
|
|||||||
enddo
|
enddo
|
||||||
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
|
END_PROVIDER
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
@ -285,8 +275,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk, (n_points_final_grid,
|
|||||||
double precision :: xyz
|
double precision :: xyz
|
||||||
double precision :: wall0, wall1
|
double precision :: wall0, wall1
|
||||||
|
|
||||||
!print*, ' providing x_W_ki_bi_ortho_erf_rk ...'
|
print*, ' providing x_W_ki_bi_ortho_erf_rk ...'
|
||||||
!call wall_time(wall0)
|
call wall_time(wall0)
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
@ -310,8 +300,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_v_ki_bi_ortho_erf_rk_cst_mu_transp
|
||||||
! FREE mo_x_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)
|
call wall_time(wall1)
|
||||||
!print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0
|
print *, ' time to provide x_W_ki_bi_ortho_erf_rk = ', wall1 - wall0
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
@ -333,8 +323,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_
|
|||||||
double precision :: xyz
|
double precision :: xyz
|
||||||
double precision :: wall0, wall1
|
double precision :: wall0, wall1
|
||||||
|
|
||||||
!print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...'
|
print*,'providing x_W_ki_bi_ortho_erf_rk_diag ...'
|
||||||
!call wall_time(wall0)
|
call wall_time(wall0)
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
@ -353,8 +343,8 @@ BEGIN_PROVIDER [ double precision, x_W_ki_bi_ortho_erf_rk_diag, (n_points_final_
|
|||||||
!$OMP END DO
|
!$OMP END DO
|
||||||
!$OMP END PARALLEL
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
!call wall_time(wall1)
|
call wall_time(wall1)
|
||||||
!print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0
|
print*,'time to provide x_W_ki_bi_ortho_erf_rk_diag = ',wall1 - wall0
|
||||||
|
|
||||||
END_PROVIDER
|
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)
|
integral = integral + tmp * final_weight_at_r_vector(ipoint)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
end
|
end subroutine give_integrals_3_body_bi_ort
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
|
|||||||
@ -16,10 +16,10 @@ double precision function bi_ortho_mo_ints(l, k, j, i)
|
|||||||
integer :: m, n, p, q
|
integer :: m, n, p, q
|
||||||
|
|
||||||
bi_ortho_mo_ints = 0.d0
|
bi_ortho_mo_ints = 0.d0
|
||||||
do p = 1, ao_num
|
do m = 1, ao_num
|
||||||
do m = 1, ao_num
|
do p = 1, ao_num
|
||||||
do q = 1, ao_num
|
do n = 1, ao_num
|
||||||
do n = 1, ao_num
|
do q = 1, ao_num
|
||||||
! p1h1p2h2 l1 l2 r1 r2
|
! 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)
|
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
|
enddo
|
||||||
@ -27,7 +27,7 @@ double precision function bi_ortho_mo_ints(l, k, j, i)
|
|||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
end
|
end function bi_ortho_mo_ints
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
@ -40,106 +40,38 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num,
|
|||||||
END_DOC
|
END_DOC
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i, j, k, l, m, n, p, q, s, r
|
integer :: i, j, k, l, m, n, p, q
|
||||||
double precision :: t1, t2, tt1, tt2
|
|
||||||
double precision, allocatable :: a1(:,:,:,:), a2(:,:,:,:)
|
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
|
PROVIDE mo_r_coef mo_l_coef
|
||||||
|
|
||||||
if(ao_to_mo_tc_n3) then
|
allocate(a2(ao_num,ao_num,ao_num,mo_num))
|
||||||
|
|
||||||
print*, ' memory scale of TC ao -> mo: O(N3) '
|
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)
|
||||||
|
|
||||||
if(.not.read_tc_integ) then
|
allocate(a1(ao_num,ao_num,mo_num,mo_num))
|
||||||
stop 'read_tc_integ needs to be set to true'
|
|
||||||
endif
|
|
||||||
|
|
||||||
allocate(a_jkp(ao_num,ao_num,mo_num))
|
call dgemm( 'T', 'N', ao_num*ao_num*mo_num, mo_num, ao_num, 1.d0 &
|
||||||
allocate(a_kpq(ao_num,mo_num,mo_num))
|
, a2(1,1,1,1), ao_num, mo_r_coef(1,1), ao_num &
|
||||||
allocate(ao_two_e_tc_tot_tmp(ao_num,ao_num,ao_num))
|
, 0.d0, a1(1,1,1,1), ao_num*ao_num*mo_num)
|
||||||
|
|
||||||
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="read")
|
deallocate(a2)
|
||||||
|
allocate(a2(ao_num,mo_num,mo_num,mo_num))
|
||||||
|
|
||||||
call wall_time(tt1)
|
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)
|
||||||
|
|
||||||
mo_bi_ortho_tc_two_e_chemist(:,:,:,:) = 0.d0
|
deallocate(a1)
|
||||||
do l = 1, ao_num
|
|
||||||
read(11) ao_two_e_tc_tot_tmp(:,:,:)
|
|
||||||
|
|
||||||
do s = 1, mo_num
|
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)
|
||||||
|
|
||||||
call dgemm( 'T', 'N', ao_num*ao_num, mo_num, ao_num, 1.d0 &
|
deallocate(a2)
|
||||||
, 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))
|
!allocate(a1(mo_num,ao_num,ao_num,ao_num))
|
||||||
!a1 = 0.d0
|
!a1 = 0.d0
|
||||||
@ -203,10 +135,6 @@ BEGIN_PROVIDER [double precision, mo_bi_ortho_tc_two_e_chemist, (mo_num, mo_num,
|
|||||||
!enddo
|
!enddo
|
||||||
!deallocate(a1)
|
!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
|
END_PROVIDER
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|||||||
@ -56,10 +56,10 @@
|
|||||||
print*,'Average trace of overlap_bi_ortho is different from 1 by ', dabs(accu_d-1.d0)
|
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*,'And bi orthogonality is off by an average of ',accu_nd
|
||||||
print*,'****************'
|
print*,'****************'
|
||||||
!print*,'Overlap matrix betwee mo_l_coef and mo_r_coef '
|
print*,'Overlap matrix betwee mo_l_coef and mo_r_coef '
|
||||||
!do i = 1, mo_num
|
do i = 1, mo_num
|
||||||
! write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:)
|
write(*,'(100(F16.10,X))')overlap_bi_ortho(i,:)
|
||||||
!enddo
|
enddo
|
||||||
endif
|
endif
|
||||||
print*,'Average trace of overlap_bi_ortho (should be 1.)'
|
print*,'Average trace of overlap_bi_ortho (should be 1.)'
|
||||||
print*,'accu_d = ',accu_d
|
print*,'accu_d = ',accu_d
|
||||||
|
|||||||
@ -45,91 +45,13 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
|
|||||||
|
|
||||||
print *, ' Numerical integration over r1 and r2 will be performed'
|
print *, ' Numerical integration over r1 and r2 will be performed'
|
||||||
|
|
||||||
if(tc_save_mem) then
|
! TODO combine 1shot & int2_grad1_u12_ao_num
|
||||||
|
|
||||||
integer :: n_blocks, n_rest, n_pass
|
PROVIDE int2_grad1_u12_ao_num
|
||||||
integer :: i_blocks, i_rest, i_pass, ii
|
int2_grad1_u12_ao = int2_grad1_u12_ao_num
|
||||||
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
|
!PROVIDE int2_grad1_u12_ao_num_1shot
|
||||||
|
!int2_grad1_u12_ao = int2_grad1_u12_ao_num_1shot
|
||||||
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
|
elseif(tc_integ_type .eq. "semi-analytic") then
|
||||||
|
|
||||||
@ -255,88 +177,13 @@ 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'
|
print *, ' Numerical integration over r1 and r2 will be performed'
|
||||||
|
|
||||||
if(tc_save_mem) then
|
! TODO combine 1shot & int2_grad1_u12_square_ao_num
|
||||||
|
|
||||||
integer :: n_blocks, n_rest, n_pass
|
PROVIDE int2_grad1_u12_square_ao_num
|
||||||
integer :: i_blocks, i_rest, i_pass, ii
|
int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num
|
||||||
double precision :: mem, n_double
|
|
||||||
double precision, allocatable :: tmp(:,:,:), xx(:,:,:)
|
|
||||||
double precision, allocatable :: tmp_grad1_u12_squared(:,:)
|
|
||||||
|
|
||||||
PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
|
!PROVIDE int2_grad1_u12_square_ao_num_1shot
|
||||||
|
!int2_grad1_u12_square_ao = int2_grad1_u12_square_ao_num_1shot
|
||||||
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
|
elseif(tc_integ_type .eq. "semi-analytic") then
|
||||||
|
|
||||||
|
|||||||
@ -131,7 +131,7 @@
|
|||||||
deallocate(tmp)
|
deallocate(tmp)
|
||||||
|
|
||||||
call wall_time(time1)
|
call wall_time(time1)
|
||||||
print*, ' wall time for int2_grad1_u12_ao_num & int2_grad1_u12_square_ao_num = (min)', (time1-time0) / 60.d0
|
print*, ' wall time for int2_grad1_u12_ao_num & int2_grad1_u12_square_ao_num =', time1-time0
|
||||||
call print_memory_usage()
|
call print_memory_usage()
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|||||||
@ -33,10 +33,8 @@ 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 :: weight1, ao_k_r, ao_i_r
|
||||||
double precision :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq
|
double precision :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq
|
||||||
double precision :: time0, time1
|
double precision :: time0, time1
|
||||||
double precision, allocatable :: c_mat(:,:,:)
|
double precision, allocatable :: b_mat(:,:,:,:), c_mat(:,:,:)
|
||||||
logical, external :: ao_two_e_integral_zero
|
|
||||||
double precision, external :: get_ao_two_e_integral
|
double precision, external :: get_ao_two_e_integral
|
||||||
double precision, external :: ao_two_e_integral
|
|
||||||
|
|
||||||
PROVIDe tc_integ_type
|
PROVIDe tc_integ_type
|
||||||
PROVIDE env_type
|
PROVIDE env_type
|
||||||
@ -55,9 +53,7 @@ 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'
|
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")
|
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="read")
|
||||||
do i = 1, ao_num
|
read(11) ao_two_e_tc_tot
|
||||||
read(11) ao_two_e_tc_tot(:,:,:,i)
|
|
||||||
enddo
|
|
||||||
close(11)
|
close(11)
|
||||||
|
|
||||||
else
|
else
|
||||||
@ -69,58 +65,26 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
|||||||
|
|
||||||
PROVIDE int2_grad1_u12_square_ao
|
PROVIDE int2_grad1_u12_square_ao
|
||||||
|
|
||||||
if(tc_save_mem_loops) then
|
allocate(c_mat(n_points_final_grid,ao_num,ao_num))
|
||||||
|
|
||||||
print*, ' LOOPS are used to evaluate Hermitian part of ao_two_e_tc_tot ...'
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PARALLEL &
|
!$OMP PRIVATE (i, k, ipoint) &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP SHARED (aos_in_r_array_transp, c_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
|
||||||
!$OMP PRIVATE (i, j, k, l, ipoint, ao_i_r, ao_k_r, weight1) &
|
!$OMP DO SCHEDULE (static)
|
||||||
!$OMP SHARED (ao_num, n_points_final_grid, ao_two_e_tc_tot, &
|
do i = 1, ao_num
|
||||||
!$OMP aos_in_r_array_transp, final_weight_at_r_vector, int2_grad1_u12_square_ao)
|
do k = 1, ao_num
|
||||||
!$OMP DO COLLAPSE(4)
|
do ipoint = 1, n_points_final_grid
|
||||||
do i = 1, ao_num
|
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)
|
||||||
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
|
enddo
|
||||||
!$OMP END DO
|
enddo
|
||||||
!$OMP END PARALLEL
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
else
|
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 &
|
||||||
print*, ' DGEMM are used to evaluate Hermitian part of ao_two_e_tc_tot ...'
|
, 0.d0, ao_two_e_tc_tot, ao_num*ao_num)
|
||||||
|
|
||||||
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
|
FREE int2_grad1_u12_square_ao
|
||||||
|
|
||||||
@ -132,7 +96,6 @@ 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
|
! an additional term is added here directly instead of
|
||||||
! being added in int2_grad1_u12_square_ao for performance
|
! 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
|
PROVIDE int2_u2_env2
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
@ -164,13 +127,10 @@ 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 &
|
, 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)
|
, 1.d0, ao_two_e_tc_tot(1,1,1,1), ao_num*ao_num)
|
||||||
|
|
||||||
deallocate(c_mat)
|
|
||||||
FREE int2_u2_env2
|
FREE int2_u2_env2
|
||||||
endif ! use_ipp
|
endif ! use_ipp
|
||||||
|
|
||||||
call wall_time(time1)
|
deallocate(c_mat)
|
||||||
print*, ' done with Hermitian part after (min) ', (time1 - time0) / 60.d0
|
|
||||||
call print_memory_usage()
|
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
@ -178,71 +138,39 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
|||||||
|
|
||||||
PROVIDE int2_grad1_u12_ao
|
PROVIDE int2_grad1_u12_ao
|
||||||
|
|
||||||
if(tc_save_mem_loops) then
|
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3))
|
||||||
|
|
||||||
print*, ' LOOPS are used to evaluate non-Hermitian part of ao_two_e_tc_tot ...'
|
!$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
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
|
||||||
!$OMP DEFAULT (NONE) &
|
ao_i_r = aos_in_r_array_transp(ipoint,i)
|
||||||
!$OMP PRIVATE (i, j, k, l, ipoint, ao_i_r, ao_k_r, weight1) &
|
ao_k_r = aos_in_r_array_transp(ipoint,k)
|
||||||
!$OMP SHARED (ao_num, n_points_final_grid, ao_two_e_tc_tot, &
|
|
||||||
!$OMP aos_in_r_array_transp, final_weight_at_r_vector, &
|
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))
|
||||||
!$OMP int2_grad1_u12_ao, aos_grad_in_r_array_transp_bis)
|
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))
|
||||||
!$OMP DO COLLAPSE(4)
|
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))
|
||||||
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
|
enddo
|
||||||
!$OMP END DO
|
enddo
|
||||||
!$OMP END PARALLEL
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
else
|
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)
|
||||||
|
|
||||||
print*, ' DGEMM are used to evaluate non-Hermitian part of ao_two_e_tc_tot ...'
|
FREE int2_grad1_u12_ao
|
||||||
|
|
||||||
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
|
if(tc_integ_type .eq. "semi-analytic") then
|
||||||
FREE int2_grad1_u2e_ao
|
FREE int2_grad1_u2e_ao
|
||||||
@ -250,67 +178,30 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
|
|||||||
|
|
||||||
endif ! var_tc
|
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)
|
call sum_A_At(ao_two_e_tc_tot(1,1,1,1), ao_num*ao_num)
|
||||||
|
|
||||||
! ---
|
PROVIDE ao_integrals_map
|
||||||
|
|
||||||
logical :: integ_zero
|
!$OMP PARALLEL DEFAULT(NONE) &
|
||||||
double precision :: integ_val
|
!$OMP SHARED(ao_num, ao_two_e_tc_tot, ao_integrals_map) &
|
||||||
|
!$OMP PRIVATE(i, j, k, l)
|
||||||
print*, ' adding ERI to ao_two_e_tc_tot ...'
|
!$OMP DO
|
||||||
|
do j = 1, ao_num
|
||||||
if(tc_save_mem) then
|
do l = 1, ao_num
|
||||||
print*, ' ao_integrals_map will not be used'
|
do i = 1, ao_num
|
||||||
!$OMP PARALLEL DEFAULT(NONE) &
|
do k = 1, ao_num
|
||||||
!$OMP PRIVATE(i, j, k, l, integ_zero, integ_val) &
|
! < 1:i, 2:j | 1:k, 2:l >
|
||||||
!$OMP SHARED(ao_num, ao_two_e_tc_tot)
|
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)
|
||||||
!$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
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
enddo
|
||||||
!$OMP END PARALLEL
|
!$OMP END DO
|
||||||
else
|
!$OMP END PARALLEL
|
||||||
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") .and. (.not. tc_save_mem)) then
|
if(tc_integ_type .eq. "numeric") then
|
||||||
FREE int2_grad1_u12_ao_num int2_grad1_u12_square_ao_num
|
FREE int2_grad1_u12_ao_num int2_grad1_u12_square_ao_num
|
||||||
endif
|
endif
|
||||||
|
|
||||||
@ -320,9 +211,7 @@ 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'
|
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")
|
open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/ao_two_e_tc_tot', action="write")
|
||||||
call ezfio_set_work_empty(.False.)
|
call ezfio_set_work_empty(.False.)
|
||||||
do i = 1, ao_num
|
write(11) ao_two_e_tc_tot
|
||||||
write(11) ao_two_e_tc_tot(:,:,:,i)
|
|
||||||
enddo
|
|
||||||
close(11)
|
close(11)
|
||||||
call ezfio_set_tc_keywords_io_tc_integ('Read')
|
call ezfio_set_tc_keywords_io_tc_integ('Read')
|
||||||
endif
|
endif
|
||||||
|
|||||||
@ -2144,7 +2144,6 @@ subroutine impose_biorthog_degen_eigvec(n, deg_num, e0, L0, R0)
|
|||||||
enddo
|
enddo
|
||||||
!print*,' accu_nd after = ', accu_nd
|
!print*,' accu_nd after = ', accu_nd
|
||||||
if(accu_nd .gt. 1d-12) then
|
if(accu_nd .gt. 1d-12) then
|
||||||
print*, ' accu_nd =', accu_nd
|
|
||||||
print*, ' your strategy for degenerates orbitals failed !'
|
print*, ' your strategy for degenerates orbitals failed !'
|
||||||
print*, m, 'deg on', i
|
print*, m, 'deg on', i
|
||||||
stop
|
stop
|
||||||
|
|||||||
@ -1 +0,0 @@
|
|||||||
|
|
||||||
@ -1,4 +0,0 @@
|
|||||||
================
|
|
||||||
normal_order_old
|
|
||||||
================
|
|
||||||
|
|
||||||
@ -1,7 +0,0 @@
|
|||||||
determinants
|
|
||||||
normal_order_old
|
|
||||||
bi_ort_ints
|
|
||||||
bi_ortho_mos
|
|
||||||
tc_keywords
|
|
||||||
non_hermit_dav
|
|
||||||
dav_general_mat
|
|
||||||
@ -1,7 +0,0 @@
|
|||||||
program slater_tc
|
|
||||||
implicit none
|
|
||||||
BEGIN_DOC
|
|
||||||
! TODO : Put the documentation of the program here
|
|
||||||
END_DOC
|
|
||||||
print *, 'Hello world'
|
|
||||||
end
|
|
||||||
@ -285,22 +285,3 @@ type: logical
|
|||||||
doc: If |true|, you minimize the angle between the left and right vectors associated to degenerate orbitals
|
doc: If |true|, you minimize the angle between the left and right vectors associated to degenerate orbitals
|
||||||
interface: ezfio,provider,ocaml
|
interface: ezfio,provider,ocaml
|
||||||
default: False
|
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 :: loc_1, loc_2, loc_3
|
||||||
double precision, allocatable :: Okappa(:), Jkappa(:,:)
|
double precision, allocatable :: Okappa(:), Jkappa(:,:)
|
||||||
double precision, allocatable :: tmp_omp_d1(:), tmp_omp_d2(:,:)
|
double precision, allocatable :: tmp_omp_d1(:), tmp_omp_d2(:,:)
|
||||||
double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:), tmp_22(:,:,:)
|
double precision, allocatable :: tmp_1(:,:), tmp_2(:,:,:,:)
|
||||||
double precision, allocatable :: tmp_3(:,:,:), tmp_4(:,:,:)
|
double precision, allocatable :: tmp_3(:,:,:), tmp_4(:,:,:)
|
||||||
|
|
||||||
PROVIDE mo_l_coef mo_r_coef
|
PROVIDE mo_l_coef mo_r_coef
|
||||||
@ -63,13 +63,17 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
|
|||||||
allocate(tmp_1(n_points_final_grid,4))
|
allocate(tmp_1(n_points_final_grid,4))
|
||||||
|
|
||||||
do ipoint = 1, n_points_final_grid
|
do ipoint = 1, n_points_final_grid
|
||||||
|
|
||||||
loc_1 = 2.d0 * Okappa(ipoint)
|
loc_1 = 2.d0 * Okappa(ipoint)
|
||||||
|
|
||||||
tmp_1(ipoint,1) = loc_1 * Jkappa(ipoint,1)
|
tmp_1(ipoint,1) = loc_1 * Jkappa(ipoint,1)
|
||||||
tmp_1(ipoint,2) = loc_1 * Jkappa(ipoint,2)
|
tmp_1(ipoint,2) = loc_1 * Jkappa(ipoint,2)
|
||||||
tmp_1(ipoint,3) = loc_1 * Jkappa(ipoint,3)
|
tmp_1(ipoint,3) = loc_1 * Jkappa(ipoint,3)
|
||||||
|
|
||||||
tmp_1(ipoint,4) = Okappa(ipoint)
|
tmp_1(ipoint,4) = Okappa(ipoint)
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
|
|
||||||
!$OMP PARALLEL &
|
!$OMP PARALLEL &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PRIVATE (ipoint, i, j, loc_1, tmp_omp_d2) &
|
!$OMP PRIVATE (ipoint, i, j, loc_1, tmp_omp_d2) &
|
||||||
@ -108,81 +112,58 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
|
|||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
if(tc_save_mem) then
|
allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num))
|
||||||
|
|
||||||
allocate(tmp_22(n_points_final_grid,4,mo_num))
|
!$OMP PARALLEL &
|
||||||
do a = 1, mo_num
|
!$OMP DEFAULT (NONE) &
|
||||||
!$OMP PARALLEL &
|
!$OMP PRIVATE (ipoint, a, b) &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP SHARED (n_points_final_grid, mo_num, &
|
||||||
!$OMP PRIVATE (ipoint, b, i) &
|
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
||||||
!$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, a, &
|
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
!$OMP tmp_2)
|
||||||
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
!$OMP DO COLLAPSE(2)
|
||||||
!$OMP tmp_22)
|
do a = 1, mo_num
|
||||||
!$OMP DO
|
do b = 1, mo_num
|
||||||
do b = 1, mo_num
|
do ipoint = 1, n_points_final_grid
|
||||||
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_22(ipoint,1,b) = 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_22(ipoint,2,b) = 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)
|
||||||
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
|
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
|
enddo
|
||||||
deallocate(tmp_22)
|
enddo
|
||||||
|
!$OMP END DO
|
||||||
|
!$OMP END PARALLEL
|
||||||
|
|
||||||
else
|
!$OMP PARALLEL &
|
||||||
|
!$OMP DEFAULT (NONE) &
|
||||||
allocate(tmp_2(n_points_final_grid,4,mo_num,mo_num))
|
!$OMP PRIVATE (ipoint, a, b, i) &
|
||||||
!$OMP PARALLEL &
|
!$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, &
|
||||||
!$OMP DEFAULT (NONE) &
|
!$OMP final_weight_at_r_vector, int2_grad1_u12_bimo_t, &
|
||||||
!$OMP PRIVATE (ipoint, a, b, i) &
|
!$OMP tmp_2)
|
||||||
!$OMP SHARED (n_points_final_grid, mo_num, elec_beta_num, &
|
!$OMP DO COLLAPSE(2)
|
||||||
!$OMP mos_l_in_r_array_transp, mos_r_in_r_array_transp, &
|
do a = 1, mo_num
|
||||||
!$OMP int2_grad1_u12_bimo_t, final_weight_at_r_vector, &
|
do b = 1, mo_num
|
||||||
!$OMP tmp_2)
|
tmp_2(:,4,b,a) = 0.d0
|
||||||
!$OMP DO COLLAPSE(2)
|
do i = 1, elec_beta_num
|
||||||
do a = 1, mo_num
|
|
||||||
do b = 1, mo_num
|
|
||||||
do ipoint = 1, n_points_final_grid
|
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,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) &
|
||||||
tmp_2(ipoint,2,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,2,b,a)
|
+ int2_grad1_u12_bimo_t(ipoint,2,b,i) * int2_grad1_u12_bimo_t(ipoint,2,i,a) &
|
||||||
tmp_2(ipoint,3,b,a) = final_weight_at_r_vector(ipoint) * int2_grad1_u12_bimo_t(ipoint,3,b,a)
|
+ int2_grad1_u12_bimo_t(ipoint,3,b,i) * int2_grad1_u12_bimo_t(ipoint,3,i,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
|
||||||
enddo
|
enddo
|
||||||
!$OMP END DO
|
enddo
|
||||||
!$OMP END PARALLEL
|
!$OMP END DO
|
||||||
call dgemv( 'T', 4*n_points_final_grid, mo_num*mo_num, -2.d0 &
|
!$OMP END PARALLEL
|
||||||
, 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
|
! ---
|
||||||
|
|
||||||
deallocate(tmp_1)
|
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)
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
@ -291,7 +272,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_cs, (mo_num, mo_num)]
|
|||||||
! ---
|
! ---
|
||||||
|
|
||||||
!call wall_time(tf)
|
!call wall_time(tf)
|
||||||
!print *, ' total Wall time for fock_3e_uhf_mo_cs =', (tf - ti) / 60.d0
|
!print *, ' total Wall time for fock_3e_uhf_mo_cs =', tf - ti
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|||||||
@ -32,7 +32,7 @@ BEGIN_PROVIDER [double precision, fock_3e_uhf_mo_a, (mo_num, mo_num)]
|
|||||||
endif
|
endif
|
||||||
|
|
||||||
!call wall_time(tf)
|
!call wall_time(tf)
|
||||||
!print *, ' Wall time for fock_3e_uhf_mo_a (min) =', (tf - ti)/60.d0
|
!print *, ' Wall time for fock_3e_uhf_mo_a =', tf - ti
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
BEGIN_DOC
|
||||||
! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis
|
! Total alpha TC Fock matrix : h_c + Two-e^TC terms on the MO basis
|
||||||
|
|||||||
@ -20,7 +20,7 @@ program minimize_tc_angles
|
|||||||
! TODO
|
! TODO
|
||||||
! check if rotations of orbitals affect the TC energy
|
! check if rotations of orbitals affect the TC energy
|
||||||
! and refuse the step
|
! and refuse the step
|
||||||
call minimize_tc_orb_angles()
|
call minimize_tc_orb_angles
|
||||||
|
|
||||||
end
|
end
|
||||||
|
|
||||||
|
|||||||
@ -40,6 +40,9 @@ subroutine LTxSxR(n, m, L, S, R, C)
|
|||||||
|
|
||||||
end subroutine LTxR
|
end subroutine LTxR
|
||||||
|
|
||||||
|
! ---
|
||||||
|
|
||||||
|
|
||||||
! ---
|
! ---
|
||||||
|
|
||||||
subroutine minimize_tc_orb_angles()
|
subroutine minimize_tc_orb_angles()
|
||||||
@ -100,10 +103,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
|
|||||||
double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
|
double precision, allocatable :: stmp(:,:), T(:,:), Snew(:,:), smat2(:,:)
|
||||||
double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:)
|
double precision, allocatable :: mo_l_coef_tmp(:,:), mo_r_coef_tmp(:,:), mo_l_coef_new(:,:)
|
||||||
|
|
||||||
PROVIDE TC_HF_energy
|
E_thr = 1d-04
|
||||||
PROVIDE mo_r_coef mo_l_coef
|
|
||||||
|
|
||||||
E_thr = 1d-07
|
|
||||||
E_old = TC_HF_energy
|
E_old = TC_HF_energy
|
||||||
allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num))
|
allocate(mo_l_coef_old(ao_num,mo_num), mo_r_coef_old(ao_num,mo_num))
|
||||||
mo_r_coef_old = mo_r_coef
|
mo_r_coef_old = mo_r_coef
|
||||||
@ -111,7 +111,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
|
|||||||
|
|
||||||
good_angles = .False.
|
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 *, ' ***************************************'
|
||||||
print *, ' ***************************************'
|
print *, ' ***************************************'
|
||||||
@ -123,7 +123,7 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
|
|||||||
mo_r_coef_good = mo_r_coef
|
mo_r_coef_good = mo_r_coef
|
||||||
mo_l_coef_good = mo_l_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
|
mo_r_coef_new = mo_r_coef
|
||||||
do i = 1, mo_num
|
do i = 1, mo_num
|
||||||
norm = 1.d0/dsqrt(overlap_mo_r(i,i))
|
norm = 1.d0/dsqrt(overlap_mo_r(i,i))
|
||||||
@ -141,11 +141,10 @@ 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 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)
|
! call give_degen(fock_diag,mo_num,thr_deg,list_degen,n_degen_list)
|
||||||
if(n_core_orb.ne.0)then
|
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
|
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
|
endif
|
||||||
|
|
||||||
print *, ' fock_matrix_mo'
|
print *, ' fock_matrix_mo'
|
||||||
do i = 1, mo_num
|
do i = 1, mo_num
|
||||||
print *, i, fock_diag(i), angle_left_right(i)
|
print *, i, fock_diag(i), angle_left_right(i)
|
||||||
@ -157,52 +156,50 @@ subroutine routine_save_rotated_mos(thr_deg, good_angles)
|
|||||||
! n_degen = ilast - ifirst +1
|
! n_degen = ilast - ifirst +1
|
||||||
|
|
||||||
n_degen = list_degen(i,0)
|
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
|
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(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(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))
|
allocate(T(n_degen,n_degen), Snew(n_degen,n_degen))
|
||||||
|
|
||||||
print*,'Right orbitals before'
|
print*,'Right orbitals before'
|
||||||
do j = 1, n_degen
|
do j = 1, n_degen
|
||||||
write(*,'(1000(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j))
|
write(*,'(100(F16.10,X))') mo_r_coef_new(1:ao_num,list_degen(i,j))
|
||||||
enddo
|
enddo
|
||||||
print*,'Left orbitals before'
|
print*,'Left orbitals before'
|
||||||
do j = 1, n_degen
|
do j = 1, n_degen
|
||||||
write(*,'(1000(F16.10,X))') mo_l_coef(1:ao_num,list_degen(i,j))
|
write(*,'(100(F16.10,X))')mo_l_coef(1:ao_num,list_degen(i,j))
|
||||||
enddo
|
enddo
|
||||||
if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then
|
if(angle_left_right(list_degen(i,1)).gt.80.d0.and.n_degen==2)then
|
||||||
integer :: i_list, j_list
|
integer :: i_list, j_list
|
||||||
i_list = list_degen(i,1)
|
i_list = list_degen(i,1)
|
||||||
j_list = list_degen(i,2)
|
j_list = list_degen(i,2)
|
||||||
print*,'Huge angle !!! == ',angle_left_right(list_degen(i,1)),angle_left_right(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 = ',i_list
|
||||||
print*,'i_list = ',j_list
|
print*,'i_list = ',j_list
|
||||||
print*,'Swapping left/right orbitals'
|
print*,'Swapping left/right orbitals'
|
||||||
call print_strong_overlap(i_list, j_list)
|
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,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_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,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)
|
mo_l_coef_tmp(1:ao_num,2) = mo_r_coef_new(1:ao_num,j_list)
|
||||||
else
|
else
|
||||||
do j = 1, n_degen
|
do j = 1, n_degen
|
||||||
print*,'i_list = ',list_degen(i,j)
|
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_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))
|
mo_l_coef_tmp(1:ao_num,j) = mo_l_coef(1:ao_num,list_degen(i,j))
|
||||||
enddo
|
enddo
|
||||||
endif
|
endif
|
||||||
print*,'Right orbitals '
|
print*,'Right orbitals '
|
||||||
do j = 1, n_degen
|
do j = 1, n_degen
|
||||||
write(*,'(1000(F16.10,X))') mo_r_coef_tmp(1:ao_num,j)
|
write(*,'(100(F16.10,X))')mo_r_coef_tmp(1:ao_num,j)
|
||||||
enddo
|
enddo
|
||||||
print*,'Left orbitals '
|
print*,'Left orbitals '
|
||||||
do j = 1, n_degen
|
do j = 1, n_degen
|
||||||
write(*,'(100(F16.10,X))') mo_l_coef_tmp(1:ao_num,j)
|
write(*,'(100(F16.10,X))')mo_l_coef_tmp(1:ao_num,j)
|
||||||
enddo
|
enddo
|
||||||
! Orthogonalization of right functions
|
! Orthogonalization of right functions
|
||||||
print *, ' Orthogonalization of RIGHT functions'
|
print *, ' Orthogonalization of RIGHT functions'
|
||||||
print *, ' ------------------------------------'
|
print *, ' ------------------------------------'
|
||||||
|
|||||||
@ -1,58 +0,0 @@
|
|||||||
! ---
|
|
||||||
|
|
||||||
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,12 +47,8 @@ END_PROVIDER
|
|||||||
END_DOC
|
END_DOC
|
||||||
|
|
||||||
implicit none
|
implicit none
|
||||||
integer :: i, j, k, l, i_count
|
integer :: i,j,k,l,i_count
|
||||||
double precision :: r(3)
|
double precision :: r(3)
|
||||||
double precision :: wall0, wall1
|
|
||||||
|
|
||||||
call wall_time(wall0)
|
|
||||||
print *, ' Providing extra_final_grid_points ...'
|
|
||||||
|
|
||||||
i_count = 0
|
i_count = 0
|
||||||
do j = 1, nucl_num
|
do j = 1, nucl_num
|
||||||
@ -70,25 +66,10 @@ END_PROVIDER
|
|||||||
index_final_points_extra(2,i_count) = i
|
index_final_points_extra(2,i_count) = i
|
||||||
index_final_points_extra(3,i_count) = j
|
index_final_points_extra(3,i_count) = j
|
||||||
index_final_points_extra_reverse(k,i,j) = i_count
|
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
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|
||||||
call wall_time(wall1)
|
|
||||||
print *, ' wall time for extra_final_grid_points,', wall1 - wall0
|
|
||||||
call print_memory_usage()
|
|
||||||
|
|
||||||
END_PROVIDER
|
END_PROVIDER
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@ -72,11 +72,7 @@ END_PROVIDER
|
|||||||
print *, ' !!! WARNING !!!'
|
print *, ' !!! WARNING !!!'
|
||||||
print *, ' negative weight !!!!'
|
print *, ' negative weight !!!!'
|
||||||
print *, i_count, final_weight_at_r_vector(i_count)
|
print *, i_count, final_weight_at_r_vector(i_count)
|
||||||
if(dabs(final_weight_at_r_vector(i_count)) .lt. 1d-10) then
|
stop
|
||||||
final_weight_at_r_vector(i_count) = 0.d0
|
|
||||||
else
|
|
||||||
stop
|
|
||||||
endif
|
|
||||||
endif
|
endif
|
||||||
enddo
|
enddo
|
||||||
enddo
|
enddo
|
||||||
|
|||||||
@ -91,42 +91,3 @@ BEGIN_PROVIDER [double precision, multi_s_dipole_moment, (N_states, N_states)]
|
|||||||
enddo
|
enddo
|
||||||
|
|
||||||
END_PROVIDER
|
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
|
|
||||||
|
|
||||||
! ---
|
|
||||||
|
|
||||||
|
|
||||||
|
|||||||
@ -1,97 +0,0 @@
|
|||||||
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
|
list_degen(n_degen_list,1) = i
|
||||||
icount = 1
|
icount = 1
|
||||||
do j = i+1, n
|
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.
|
is_ok(j) = .False.
|
||||||
icount += 1
|
icount += 1
|
||||||
list_degen(n_degen_list,icount) = j
|
list_degen(n_degen_list,icount) = j
|
||||||
|
|||||||
@ -18,30 +18,3 @@ double precision, parameter :: c_4_3 = 4.d0/3.d0
|
|||||||
double precision, parameter :: c_1_3 = 1.d0/3.d0
|
double precision, parameter :: c_1_3 = 1.d0/3.d0
|
||||||
double precision, parameter :: sq_op5 = dsqrt(0.5d0)
|
double precision, parameter :: sq_op5 = dsqrt(0.5d0)
|
||||||
double precision, parameter :: dlog_2pi = dlog(2.d0*dacos(-1.d0))
|
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