mirror of
https://github.com/QuantumPackage/qp2.git
synced 2024-12-22 20:34:58 +01:00
Merge pull request #256 from AbdAmmar/dev-stable-tc-scf
Dev stable tc scf
This commit is contained in:
commit
298c91f718
@ -18,6 +18,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
|
||||
double precision :: int_gauss, dsqpi_3_2, int_j1b
|
||||
double precision :: factor_ij_1s, beta_ij, center_ij_1s(3), sq_pi_3_2
|
||||
double precision, allocatable :: int_fit_v(:)
|
||||
double precision, external :: overlap_gauss_r12_ao
|
||||
double precision, external :: overlap_gauss_r12_ao_with1s
|
||||
|
||||
print*, ' providing int2_grad1u2_grad2u2_j1b2_test ...'
|
||||
@ -39,48 +40,61 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
|
||||
!$OMP List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b, &
|
||||
!$OMP ao_overlap_abs,sq_pi_3_2)
|
||||
!$OMP DO SCHEDULE(dynamic)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
do i_1s = 1, List_comb_thr_b3_size(j,i)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(ao_overlap_abs(j,i) .lt. 1.d-12) then
|
||||
cycle
|
||||
endif
|
||||
|
||||
coef = List_comb_thr_b3_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b3_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
|
||||
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
|
||||
! --- --- ---
|
||||
! i_1s = 1
|
||||
! --- --- ---
|
||||
|
||||
int_j1b = ao_abs_comb_b3_j1b(1,j,i)
|
||||
do i_fit = 1, ng_fit_jast
|
||||
expo_fit = expo_gauss_1_erf_x_2(i_fit)
|
||||
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit)
|
||||
if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle
|
||||
int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j)
|
||||
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss
|
||||
enddo
|
||||
|
||||
! --- --- ---
|
||||
! i_1s > 1
|
||||
! --- --- ---
|
||||
|
||||
do i_fit = 1, ng_fit_jast
|
||||
do i_1s = 2, List_comb_thr_b3_size(j,i)
|
||||
|
||||
coef = List_comb_thr_b3_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b3_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
|
||||
B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
|
||||
|
||||
expo_fit = expo_gauss_1_erf_x_2(i_fit)
|
||||
!DIR$ FORCEINLINE
|
||||
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
|
||||
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef
|
||||
do i_fit = 1, ng_fit_jast
|
||||
expo_fit = expo_gauss_1_erf_x_2(i_fit)
|
||||
!DIR$ FORCEINLINE
|
||||
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
|
||||
coef_fit = -0.25d0 * coef_gauss_1_erf_x_2(i_fit) * coef
|
||||
! if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version
|
||||
if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
|
||||
|
||||
if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
|
||||
! call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, &
|
||||
! expo_fit, i, j, int_fit_v, n_points_final_grid)
|
||||
int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
|
||||
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss
|
||||
|
||||
enddo
|
||||
int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
@ -239,9 +253,27 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_comb_thr_b3_size(j,i)
|
||||
|
||||
! --- --- ---
|
||||
! i_1s = 1
|
||||
! --- --- ---
|
||||
|
||||
int_j1b = ao_abs_comb_b3_j1b(1,j,i)
|
||||
if(dabs(int_j1b).lt.1.d-10) cycle
|
||||
do i_fit = 1, ng_fit_jast
|
||||
expo_fit = expo_gauss_j_mu_x_2(i_fit)
|
||||
coef_fit = coef_gauss_j_mu_x_2(i_fit)
|
||||
if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle
|
||||
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
|
||||
tmp += coef_fit * int_fit
|
||||
enddo
|
||||
|
||||
! --- --- ---
|
||||
! i_1s > 1
|
||||
! --- --- ---
|
||||
|
||||
do i_1s = 2, List_comb_thr_b3_size(j,i)
|
||||
|
||||
coef = List_comb_thr_b3_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b3_expo (i_1s,j,i)
|
||||
@ -252,23 +284,15 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
|
||||
B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
|
||||
|
||||
do i_fit = 1, ng_fit_jast
|
||||
|
||||
expo_fit = expo_gauss_j_mu_x_2(i_fit)
|
||||
coef_fit = coef_gauss_j_mu_x_2(i_fit)
|
||||
!DIR$ FORCEINLINE
|
||||
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
|
||||
! if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version
|
||||
if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
|
||||
|
||||
! ---
|
||||
|
||||
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
|
||||
tmp += coef * coef_fit * int_fit
|
||||
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
tmp += coef * coef_fit * int_fit
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
enddo
|
||||
|
||||
int2_u2_j1b2_test(j,i,ipoint) = tmp
|
||||
@ -451,13 +475,34 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do i = 1, ao_num
|
||||
do j = i, ao_num
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
|
||||
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10) cycle
|
||||
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
r(3) = final_grid_points(3,ipoint)
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_comb_thr_b3_size(j,i)
|
||||
|
||||
! --- --- ---
|
||||
! i_1s = 1
|
||||
! --- --- ---
|
||||
|
||||
int_j1b = ao_abs_comb_b3_j1b(1,j,i)
|
||||
if(dabs(int_j1b).lt.1.d-10) cycle
|
||||
do i_fit = 1, ng_fit_jast
|
||||
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
|
||||
if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.1.d-15) cycle
|
||||
coef_fit = coef_gauss_j_mu_1_erf(i_fit)
|
||||
int_fit = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r, 1.d+9, r)
|
||||
tmp += coef_fit * int_fit
|
||||
enddo
|
||||
|
||||
! --- --- ---
|
||||
! i_1s > 1
|
||||
! --- --- ---
|
||||
|
||||
do i_1s = 2, List_comb_thr_b3_size(j,i)
|
||||
|
||||
coef = List_comb_thr_b3_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b3_expo (i_1s,j,i)
|
||||
@ -469,9 +514,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
|
||||
dist = (B_center(1) - r(1)) * (B_center(1) - r(1)) &
|
||||
+ (B_center(2) - r(2)) * (B_center(2) - r(2)) &
|
||||
+ (B_center(3) - r(3)) * (B_center(3) - r(3))
|
||||
|
||||
do i_fit = 1, ng_fit_jast
|
||||
|
||||
expo_fit = expo_gauss_j_mu_1_erf(i_fit)
|
||||
call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
|
||||
if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.1.d-15)cycle
|
||||
|
@ -192,9 +192,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
|
||||
double precision :: coef, beta, B_center(3)
|
||||
double precision :: tmp
|
||||
double precision :: wall0, wall1
|
||||
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot
|
||||
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b
|
||||
|
||||
double precision, external :: overlap_gauss_r12_ao
|
||||
double precision, external :: overlap_gauss_r12_ao_with1s
|
||||
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b
|
||||
|
||||
print*, ' providing v_ij_u_cst_mu_j1b_test ...'
|
||||
|
||||
@ -208,15 +210,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
|
||||
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, &
|
||||
!$OMP final_grid_points, ng_fit_jast, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, &
|
||||
!$OMP final_grid_points, ng_fit_jast, &
|
||||
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
|
||||
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
|
||||
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, &
|
||||
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
|
||||
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b, &
|
||||
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
@ -227,8 +228,26 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_comb_thr_b2_size(j,i)
|
||||
|
||||
! --- --- ---
|
||||
! i_1s = 1
|
||||
! --- --- ---
|
||||
|
||||
int_j1b = ao_abs_comb_b2_j1b(1,j,i)
|
||||
if(dabs(int_j1b).lt.1.d-10) cycle
|
||||
do i_fit = 1, ng_fit_jast
|
||||
expo_fit = expo_gauss_j_mu_x(i_fit)
|
||||
coef_fit = coef_gauss_j_mu_x(i_fit)
|
||||
if(ao_overlap_abs_grid(j,i).lt.1.d-15) cycle
|
||||
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
|
||||
tmp += coef_fit * int_fit
|
||||
enddo
|
||||
|
||||
! --- --- ---
|
||||
! i_1s > 1
|
||||
! --- --- ---
|
||||
|
||||
do i_1s = 2, List_comb_thr_b2_size(j,i)
|
||||
coef = List_comb_thr_b2_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b2_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
|
||||
@ -236,18 +255,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
|
||||
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
|
||||
|
||||
do i_fit = 1, ng_fit_jast
|
||||
|
||||
expo_fit = expo_gauss_j_mu_x(i_fit)
|
||||
coef_fit = coef_gauss_j_mu_x(i_fit)
|
||||
coeftot = coef * coef_fit
|
||||
if(dabs(coeftot).lt.1.d-15)cycle
|
||||
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot
|
||||
call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
|
||||
if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
|
||||
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
|
||||
tmp += coef * coef_fit * int_fit
|
||||
enddo
|
||||
enddo
|
||||
@ -288,9 +303,12 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
|
||||
double precision :: coef, beta, B_center(3)
|
||||
double precision :: tmp
|
||||
double precision :: wall0, wall1
|
||||
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3), coeftot
|
||||
double precision :: sigma_ij, dist_ij_ipoint, dsqpi_3_2, int_j1b
|
||||
|
||||
double precision, external :: overlap_gauss_r12_ao
|
||||
double precision, external :: overlap_gauss_r12_ao_with1s
|
||||
double precision :: sigma_ij,dist_ij_ipoint,dsqpi_3_2,int_j1b
|
||||
|
||||
dsqpi_3_2 = (dacos(-1.d0))**(1.5d0)
|
||||
|
||||
provide mu_erf final_grid_points j1b_pen
|
||||
@ -299,17 +317,16 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
|
||||
v_ij_u_cst_mu_j1b_ng_1_test = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
|
||||
!$OMP PRIVATE (ipoint, i, j, i_1s, r, coef, beta, B_center, &
|
||||
!$OMP beta_ij_u, factor_ij_1s_u, center_ij_1s_u, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, &
|
||||
!$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, &
|
||||
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
|
||||
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
|
||||
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, &
|
||||
!$OMP coef_fit, expo_fit, int_fit, tmp,coeftot,int_j1b) &
|
||||
!$OMP SHARED (n_points_final_grid, ao_num, &
|
||||
!$OMP final_grid_points, expo_good_j_mu_1gauss,coef_good_j_mu_1gauss, &
|
||||
!$OMP expo_gauss_j_mu_x, coef_gauss_j_mu_x, &
|
||||
!$OMP List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, &
|
||||
!$OMP List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, &
|
||||
!$OMP ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
|
||||
!$OMP DO
|
||||
!do ipoint = 1, 10
|
||||
do ipoint = 1, n_points_final_grid
|
||||
r(1) = final_grid_points(1,ipoint)
|
||||
r(2) = final_grid_points(2,ipoint)
|
||||
@ -320,8 +337,22 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
|
||||
if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
|
||||
|
||||
tmp = 0.d0
|
||||
do i_1s = 1, List_comb_thr_b2_size(j,i)
|
||||
|
||||
! --- --- ---
|
||||
! i_1s = 1
|
||||
! --- --- ---
|
||||
|
||||
int_j1b = ao_abs_comb_b2_j1b(1,j,i)
|
||||
if(dabs(int_j1b).lt.1.d-10) cycle
|
||||
expo_fit = expo_good_j_mu_1gauss
|
||||
int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
|
||||
tmp += int_fit
|
||||
|
||||
! --- --- ---
|
||||
! i_1s > 1
|
||||
! --- --- ---
|
||||
|
||||
do i_1s = 2, List_comb_thr_b2_size(j,i)
|
||||
coef = List_comb_thr_b2_coef (i_1s,j,i)
|
||||
beta = List_comb_thr_b2_expo (i_1s,j,i)
|
||||
int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
|
||||
@ -329,18 +360,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
|
||||
B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
|
||||
B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
|
||||
B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
|
||||
|
||||
! do i_fit = 1, ng_fit_jast
|
||||
|
||||
expo_fit = expo_good_j_mu_1gauss
|
||||
coef_fit = 1.d0
|
||||
coeftot = coef * coef_fit
|
||||
if(dabs(coeftot).lt.1.d-15)cycle
|
||||
double precision :: beta_ij_u, factor_ij_1s_u, center_ij_1s_u(3),coeftot
|
||||
call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
|
||||
if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
|
||||
int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
|
||||
|
||||
tmp += coef * coef_fit * int_fit
|
||||
! enddo
|
||||
enddo
|
||||
|
@ -102,11 +102,11 @@ END_PROVIDER
|
||||
List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i))
|
||||
enddo
|
||||
|
||||
print *, ' coeff, expo & cent of list b2'
|
||||
do i = 1, List_all_comb_b2_size
|
||||
print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i)
|
||||
print*, List_all_comb_b2_cent(1,i), List_all_comb_b2_cent(2,i), List_all_comb_b2_cent(3,i)
|
||||
enddo
|
||||
!print *, ' coeff, expo & cent of list b2'
|
||||
!do i = 1, List_all_comb_b2_size
|
||||
! print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i)
|
||||
! print*, List_all_comb_b2_cent(1,i), List_all_comb_b2_cent(2,i), List_all_comb_b2_cent(3,i)
|
||||
!enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -225,11 +225,11 @@ END_PROVIDER
|
||||
List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i))
|
||||
enddo
|
||||
|
||||
print *, ' coeff, expo & cent of list b3'
|
||||
do i = 1, List_all_comb_b3_size
|
||||
print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i)
|
||||
print*, List_all_comb_b3_cent(1,i), List_all_comb_b3_cent(2,i), List_all_comb_b3_cent(3,i)
|
||||
enddo
|
||||
!print *, ' coeff, expo & cent of list b3'
|
||||
!do i = 1, List_all_comb_b3_size
|
||||
! print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i)
|
||||
! print*, List_all_comb_b3_cent(1,i), List_all_comb_b3_cent(2,i), List_all_comb_b3_cent(3,i)
|
||||
!enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -1,4 +1,25 @@
|
||||
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_two_e_vartc_tot, (ao_num, ao_num, ao_num, ao_num) ]
|
||||
|
||||
integer :: i, j, k, l
|
||||
|
||||
provide j1b_type
|
||||
provide mo_r_coef mo_l_coef
|
||||
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
ao_two_e_vartc_tot(k,i,l,j) = ao_vartc_int_chemist(k,i,l,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
|
||||
END_PROVIDER
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_num) ]
|
||||
|
@ -299,7 +299,6 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ]
|
||||
@ -364,70 +363,100 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
|
||||
integer :: ipoint, i, j, k, l
|
||||
double precision :: weight1, ao_ik_r, ao_i_r
|
||||
double precision :: time0, time1
|
||||
double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:), tmp(:,:,:)
|
||||
double precision, allocatable :: b_mat(:,:,:), tmp(:,:,:)
|
||||
|
||||
print*, ' providing tc_grad_square_ao ...'
|
||||
call wall_time(time0)
|
||||
|
||||
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num), b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
|
||||
if(read_tc_integ) then
|
||||
|
||||
b_mat = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, k, ipoint) &
|
||||
!$OMP SHARED (aos_in_r_array_transp, b_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
|
||||
b_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
|
||||
|
||||
tmp = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (j, l, ipoint) &
|
||||
!$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
tmp(l,j,ipoint) = u12sq_j1bsq(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(l,j,ipoint) + 0.5d0 * grad12_j12(l,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
|
||||
ac_mat = 0.d0
|
||||
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
|
||||
, tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
|
||||
, 1.d0, ac_mat, ao_num*ao_num)
|
||||
deallocate(tmp, b_mat)
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l) &
|
||||
!$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
open(unit=11, form="unformatted", file='tc_grad_square_ao', action="read")
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
read(11) tc_grad_square_ao(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
|
||||
else
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
|
||||
|
||||
b_mat = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, k, ipoint) &
|
||||
!$OMP SHARED (aos_in_r_array_transp, b_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
|
||||
b_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
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
tmp = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (j, l, ipoint) &
|
||||
!$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
tmp(l,j,ipoint) = u12sq_j1bsq(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(l,j,ipoint) + 0.5d0 * grad12_j12(l,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
tc_grad_square_ao = 0.d0
|
||||
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
|
||||
, tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
|
||||
, 1.d0, tc_grad_square_ao, ao_num*ao_num)
|
||||
deallocate(tmp, b_mat)
|
||||
|
||||
call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)
|
||||
|
||||
!!$OMP PARALLEL &
|
||||
!!$OMP DEFAULT (NONE) &
|
||||
!!$OMP PRIVATE (i, j, k, l) &
|
||||
!!$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num)
|
||||
!!$OMP DO SCHEDULE (static)
|
||||
! do j = 1, ao_num
|
||||
! do l = 1, ao_num
|
||||
! do i = 1, ao_num
|
||||
! do k = 1, ao_num
|
||||
! tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
!!$OMP END DO
|
||||
!!$OMP END PARALLEL
|
||||
endif
|
||||
|
||||
deallocate(ac_mat)
|
||||
if(write_tc_integ) then
|
||||
open(unit=11, form="unformatted", file='tc_grad_square_ao', action="write")
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
write(11) tc_grad_square_ao(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
endif
|
||||
|
||||
call wall_time(time1)
|
||||
print*, ' Wall time for tc_grad_square_ao = ', time1 - time0
|
||||
|
@ -11,11 +11,177 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
|
||||
integer :: ipoint, i, j, k, l
|
||||
double precision :: weight1, ao_ik_r, ao_i_r,contrib,contrib2
|
||||
double precision :: time0, time1
|
||||
double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:), tmp(:,:,:)
|
||||
double precision, allocatable :: b_mat(:,:,:), tmp(:,:,:)
|
||||
|
||||
print*, ' providing tc_grad_square_ao_test ...'
|
||||
call wall_time(time0)
|
||||
|
||||
if(read_tc_integ) then
|
||||
|
||||
open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="read")
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
read(11) tc_grad_square_ao_test(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
|
||||
else
|
||||
|
||||
provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
|
||||
|
||||
b_mat = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, k, ipoint) &
|
||||
!$OMP SHARED (aos_in_r_array_transp, b_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
|
||||
b_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
|
||||
|
||||
tmp = 0.d0
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (j, l, ipoint) &
|
||||
!$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
tmp(l,j,ipoint) = u12sq_j1bsq_test(l,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b_test(l,j,ipoint) + 0.5d0 * grad12_j12_test(l,j,ipoint)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
tc_grad_square_ao_test = 0.d0
|
||||
call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
|
||||
, tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
|
||||
, 1.d0, tc_grad_square_ao_test, ao_num*ao_num)
|
||||
deallocate(tmp, b_mat)
|
||||
|
||||
call sum_A_At(tc_grad_square_ao_test(1,1,1,1), ao_num*ao_num)
|
||||
!do i = 1, ao_num
|
||||
! do j = 1, ao_num
|
||||
! do k = i, ao_num
|
||||
|
||||
! do l = max(j,k), ao_num
|
||||
! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
|
||||
! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
|
||||
! end do
|
||||
|
||||
! !if (j.eq.k) then
|
||||
! ! do l = j+1, ao_num
|
||||
! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
|
||||
! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
|
||||
! ! end do
|
||||
! !else
|
||||
! ! do l = j, ao_num
|
||||
! ! tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
|
||||
! ! tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
|
||||
! ! enddo
|
||||
! !endif
|
||||
|
||||
! enddo
|
||||
! enddo
|
||||
!enddo
|
||||
!tc_grad_square_ao_test = 2.d0 * tc_grad_square_ao_test
|
||||
! !$OMP PARALLEL &
|
||||
! !$OMP DEFAULT (NONE) &
|
||||
! !$OMP PRIVATE (i, j, k, l) &
|
||||
! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
|
||||
! !$OMP DO SCHEDULE (static)
|
||||
! integer :: ii
|
||||
! ii = 0
|
||||
! do j = 1, ao_num
|
||||
! do l = 1, ao_num
|
||||
! do i = 1, ao_num
|
||||
! do k = 1, ao_num
|
||||
! if((i.lt.j) .and. (k.lt.l)) cycle
|
||||
! ii = ii + 1
|
||||
! tc_grad_square_ao_test(k,i,l,j) = tc_grad_square_ao_test(k,i,l,j) + tc_grad_square_ao_test(l,j,k,i)
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! print *, ' ii =', ii
|
||||
! !$OMP END DO
|
||||
! !$OMP END PARALLEL
|
||||
|
||||
! !$OMP PARALLEL &
|
||||
! !$OMP DEFAULT (NONE) &
|
||||
! !$OMP PRIVATE (i, j, k, l) &
|
||||
! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
|
||||
! !$OMP DO SCHEDULE (static)
|
||||
! do j = 1, ao_num
|
||||
! do l = 1, ao_num
|
||||
! do i = 1, j-1
|
||||
! do k = 1, l-1
|
||||
! ii = ii + 1
|
||||
! tc_grad_square_ao_test(k,i,l,j) = tc_grad_square_ao_test(l,j,k,i)
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! print *, ' ii =', ii
|
||||
! print *, ao_num * ao_num * ao_num * ao_num
|
||||
! !$OMP END DO
|
||||
! !$OMP END PARALLEL
|
||||
|
||||
endif
|
||||
|
||||
if(write_tc_integ) then
|
||||
open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="write")
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
write(11) tc_grad_square_ao_test(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
endif
|
||||
|
||||
call wall_time(time1)
|
||||
print*, ' Wall time for tc_grad_square_ao_test = ', time1 - time0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, ao_num, ao_num)]
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! tc_grad_square_ao_test_ref(k,i,l,j) = -1/2 <kl | |\grad_1 u(r1,r2)|^2 + |\grad_1 u(r1,r2)|^2 | ij>
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: ipoint, i, j, k, l
|
||||
double precision :: weight1, ao_ik_r, ao_i_r,contrib,contrib2
|
||||
double precision :: time0, time1
|
||||
double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:), tmp(:,:,:)
|
||||
|
||||
print*, ' providing tc_grad_square_ao_test_ref ...'
|
||||
call wall_time(time0)
|
||||
|
||||
provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test
|
||||
|
||||
allocate(ac_mat(ao_num,ao_num,ao_num,ao_num), b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
|
||||
@ -61,13 +227,13 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l) &
|
||||
!$OMP SHARED (ac_mat, tc_grad_square_ao_test, ao_num)
|
||||
!$OMP SHARED (ac_mat, tc_grad_square_ao_test_ref, ao_num)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
tc_grad_square_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
tc_grad_square_ao_test_ref(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
@ -78,7 +244,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
|
||||
deallocate(ac_mat)
|
||||
|
||||
call wall_time(time1)
|
||||
print*, ' Wall time for tc_grad_square_ao_test = ', time1 - time0
|
||||
print*, ' Wall time for tc_grad_square_ao_test_ref = ', time1 - time0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -25,7 +25,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: ipoint, i, j
|
||||
integer :: ipoint, i, j, m
|
||||
double precision :: time0, time1
|
||||
double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
|
||||
|
||||
@ -33,54 +33,76 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_
|
||||
call wall_time(time0)
|
||||
|
||||
PROVIDE j1b_type
|
||||
|
||||
if(j1b_type .eq. 3) then
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
if(read_tc_integ) then
|
||||
|
||||
tmp0 = 0.5d0 * v_1b(ipoint)
|
||||
tmp_x = v_1b_grad(1,ipoint)
|
||||
tmp_y = v_1b_grad(2,ipoint)
|
||||
tmp_z = v_1b_grad(3,ipoint)
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
|
||||
tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
|
||||
|
||||
int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
|
||||
int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
|
||||
int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
|
||||
open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="read")
|
||||
do m = 1, 3
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
read(11) int2_grad1_u12_ao(i,j,ipoint,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
|
||||
else
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
|
||||
|
||||
int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,1)
|
||||
int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,2)
|
||||
int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,3)
|
||||
if(j1b_type .eq. 3) then
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
tmp0 = 0.5d0 * v_1b(ipoint)
|
||||
tmp_x = v_1b_grad(1,ipoint)
|
||||
tmp_y = v_1b_grad(2,ipoint)
|
||||
tmp_z = v_1b_grad(3,ipoint)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint)
|
||||
tmp2 = v_ij_u_cst_mu_j1b(i,j,ipoint)
|
||||
int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
|
||||
int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
|
||||
int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
|
||||
|
||||
int2_grad1_u12_ao *= 0.5d0
|
||||
int2_grad1_u12_ao(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,1)
|
||||
int2_grad1_u12_ao(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,2)
|
||||
int2_grad1_u12_ao(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_transp_bis(ipoint,i,j,3)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
int2_grad1_u12_ao *= 0.5d0
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
if(write_tc_integ) then
|
||||
open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="write")
|
||||
do m = 1, 3
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
write(11) int2_grad1_u12_ao(i,j,ipoint,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
endif
|
||||
|
||||
call wall_time(time1)
|
||||
print*, ' Wall time for int2_grad1_u12_ao = ', time1 - time0
|
||||
|
||||
@ -290,65 +312,95 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao, (ao_num, ao_num, ao_num,
|
||||
integer :: ipoint, i, j, k, l, m
|
||||
double precision :: weight1, ao_k_r, ao_i_r
|
||||
double precision :: time0, time1
|
||||
double precision, allocatable :: ac_mat(:,:,:,:), b_mat(:,:,:,:)
|
||||
double precision, allocatable :: b_mat(:,:,:,:)
|
||||
|
||||
print*, ' providing tc_grad_and_lapl_ao ...'
|
||||
call wall_time(time0)
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num))
|
||||
if(read_tc_integ) then
|
||||
|
||||
b_mat = 0.d0
|
||||
!$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
|
||||
|
||||
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))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
ac_mat = 0.d0
|
||||
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, ac_mat, ao_num*ao_num)
|
||||
|
||||
enddo
|
||||
deallocate(b_mat)
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l) &
|
||||
!$OMP SHARED (ac_mat, tc_grad_and_lapl_ao, ao_num)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao', action="read")
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
!tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j)
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
read(11) tc_grad_and_lapl_ao(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
|
||||
else
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3))
|
||||
|
||||
b_mat = 0.d0
|
||||
!$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
|
||||
|
||||
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))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
tc_grad_and_lapl_ao = 0.d0
|
||||
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, tc_grad_and_lapl_ao, ao_num*ao_num)
|
||||
|
||||
enddo
|
||||
deallocate(b_mat)
|
||||
|
||||
call sum_A_At(tc_grad_and_lapl_ao(1,1,1,1), ao_num*ao_num)
|
||||
! !$OMP PARALLEL &
|
||||
! !$OMP DEFAULT (NONE) &
|
||||
! !$OMP PRIVATE (i, j, k, l) &
|
||||
! !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao, ao_num)
|
||||
! !$OMP DO SCHEDULE (static)
|
||||
! do j = 1, ao_num
|
||||
! do l = 1, ao_num
|
||||
! do i = 1, ao_num
|
||||
! do k = 1, ao_num
|
||||
! tc_grad_and_lapl_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! enddo
|
||||
! !$OMP END DO
|
||||
! !$OMP END PARALLEL
|
||||
|
||||
deallocate(ac_mat)
|
||||
endif
|
||||
|
||||
if(write_tc_integ) then
|
||||
open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao', action="write")
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
write(11) tc_grad_and_lapl_ao(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
endif
|
||||
|
||||
call wall_time(time1)
|
||||
print*, ' Wall time for tc_grad_and_lapl_ao = ', time1 - time0
|
||||
|
@ -24,7 +24,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
integer :: ipoint, i, j
|
||||
integer :: ipoint, i, j, m
|
||||
double precision :: time0, time1
|
||||
double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
|
||||
|
||||
@ -32,52 +32,73 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
|
||||
call wall_time(time0)
|
||||
|
||||
PROVIDE j1b_type
|
||||
|
||||
if(j1b_type .eq. 3) then
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
if(read_tc_integ) then
|
||||
|
||||
tmp0 = 0.5d0 * v_1b(ipoint)
|
||||
tmp_x = v_1b_grad(1,ipoint)
|
||||
tmp_y = v_1b_grad(2,ipoint)
|
||||
tmp_z = v_1b_grad(3,ipoint)
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
|
||||
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint)
|
||||
tmp2 = v_ij_u_cst_mu_j1b_test(i,j,ipoint)
|
||||
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,1) - tmp2 * tmp_x
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,2) - tmp2 * tmp_y
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,3) - tmp2 * tmp_z
|
||||
open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="read")
|
||||
do m = 1, 3
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
read(11) int2_grad1_u12_ao_test(i,j,ipoint,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
|
||||
else
|
||||
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
|
||||
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,1)
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,2)
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,3)
|
||||
|
||||
if(j1b_type .eq. 3) then
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
tmp0 = 0.5d0 * v_1b(ipoint)
|
||||
tmp_x = v_1b_grad(1,ipoint)
|
||||
tmp_y = v_1b_grad(2,ipoint)
|
||||
tmp_z = v_1b_grad(3,ipoint)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
tmp1 = tmp0 * v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint)
|
||||
tmp2 = v_ij_u_cst_mu_j1b_test(i,j,ipoint)
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,1) - tmp2 * tmp_x
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,2) - tmp2 * tmp_y
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - tmp0 * x_v_ij_erf_rk_cst_mu_j1b_test(i,j,ipoint,3) - tmp2 * tmp_z
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do ipoint = 1, n_points_final_grid
|
||||
x = final_grid_points(1,ipoint)
|
||||
y = final_grid_points(2,ipoint)
|
||||
z = final_grid_points(3,ipoint)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,1) = tmp1 * x - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,1)
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,2) = tmp1 * y - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,2)
|
||||
int2_grad1_u12_ao_test(i,j,ipoint,3) = tmp1 * z - x_v_ij_erf_rk_cst_mu_tmp(i,j,ipoint,3)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
int2_grad1_u12_ao_test *= 0.5d0
|
||||
endif
|
||||
|
||||
int2_grad1_u12_ao_test *= 0.5d0
|
||||
endif
|
||||
|
||||
if(write_tc_integ) then
|
||||
open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="write")
|
||||
do m = 1, 3
|
||||
do ipoint = 1, n_points_final_grid
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
write(11) int2_grad1_u12_ao_test(i,j,ipoint,m)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
endif
|
||||
|
||||
call wall_time(time1)
|
||||
@ -109,61 +130,93 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
|
||||
print*, ' providing tc_grad_and_lapl_ao_test ...'
|
||||
call wall_time(time0)
|
||||
|
||||
provide int2_grad1_u12_ao_test
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num))
|
||||
|
||||
b_mat = 0.d0
|
||||
!$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
|
||||
|
||||
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))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
ac_mat = 0.d0
|
||||
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_test(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
|
||||
, 1.d0, ac_mat, ao_num*ao_num)
|
||||
|
||||
enddo
|
||||
deallocate(b_mat)
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l) &
|
||||
!$OMP SHARED (ac_mat, tc_grad_and_lapl_ao_test, ao_num)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
if(read_tc_integ) then
|
||||
|
||||
open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="read")
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
tc_grad_and_lapl_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
read(11) tc_grad_and_lapl_ao_test(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
|
||||
else
|
||||
|
||||
provide int2_grad1_u12_ao_test
|
||||
|
||||
allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num))
|
||||
|
||||
b_mat = 0.d0
|
||||
!$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
|
||||
|
||||
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))
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
ac_mat = 0.d0
|
||||
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_test(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
|
||||
, 1.d0, ac_mat, ao_num*ao_num)
|
||||
|
||||
enddo
|
||||
deallocate(b_mat)
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l) &
|
||||
!$OMP SHARED (ac_mat, tc_grad_and_lapl_ao_test, ao_num)
|
||||
!$OMP DO SCHEDULE (static)
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
tc_grad_and_lapl_ao_test(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
deallocate(ac_mat)
|
||||
|
||||
deallocate(ac_mat)
|
||||
endif
|
||||
|
||||
if(write_tc_integ) then
|
||||
open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="write")
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
write(11) tc_grad_and_lapl_ao_test(l,k,j,i)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
close(11)
|
||||
endif
|
||||
|
||||
call wall_time(time1)
|
||||
print*, ' Wall time for tc_grad_and_lapl_ao_test = ', time1 - time0
|
||||
|
@ -1,6 +1,44 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num, ao_num)]
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, l
|
||||
double precision :: wall1, wall0
|
||||
|
||||
print *, ' providing ao_vartc_int_chemist ...'
|
||||
call wall_time(wall0)
|
||||
|
||||
if(test_cycle_tc) then
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
ao_vartc_int_chemist(k,i,l,j) = tc_grad_square_ao_test(k,i,l,j) + ao_two_e_coul(k,i,l,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
ao_vartc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
call wall_time(wall1)
|
||||
print *, ' wall time for ao_vartc_int_chemist ', wall1 - wall0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao_num)]
|
||||
|
||||
implicit none
|
||||
@ -11,17 +49,17 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
|
||||
call wall_time(wall0)
|
||||
|
||||
if(test_cycle_tc)then
|
||||
ao_tc_int_chemist = ao_tc_int_chemist_test
|
||||
ao_tc_int_chemist = ao_tc_int_chemist_test
|
||||
else
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
ao_tc_int_chemist(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
call wall_time(wall1)
|
||||
|
@ -86,7 +86,7 @@ default: False
|
||||
type: Threshold
|
||||
doc: Threshold on the convergence of the Hartree Fock energy.
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-12
|
||||
default: 1.e-10
|
||||
|
||||
[n_it_tcscf_max]
|
||||
type: Strictly_positive_int
|
||||
@ -94,6 +94,12 @@ doc: Maximum number of SCF iterations
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 100
|
||||
|
||||
[selection_tc]
|
||||
type: integer
|
||||
doc: if +1: only positive is selected, -1: only negative is selected, :0 both positive and negative
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0
|
||||
|
||||
[j1b_pen]
|
||||
type: double precision
|
||||
doc: exponents of the 1-body Jastrow
|
||||
@ -130,12 +136,30 @@ doc: nb of Gaussians used to fit Jastrow fcts
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 20
|
||||
|
||||
[max_dim_diis_tcscf]
|
||||
type: integer
|
||||
doc: Maximum size of the DIIS extrapolation procedure
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 15
|
||||
|
||||
[level_shift_tcscf]
|
||||
type: Positive_float
|
||||
doc: Energy shift on the virtual MOs to improve TCSCF convergence
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.
|
||||
|
||||
[tcscf_algorithm]
|
||||
type: character*(32)
|
||||
doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS]
|
||||
interface: ezfio,provider,ocaml
|
||||
default: Simple
|
||||
|
||||
[im_thresh_tcscf]
|
||||
type: Threshold
|
||||
doc: Thresholds on the Imag part of energy
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-7
|
||||
|
||||
[test_cycle_tc]
|
||||
type: logical
|
||||
doc: If |true|, the integrals of the three-body jastrow are computed with cycles
|
||||
@ -154,29 +178,23 @@ doc: Threshold to determine if non-diagonal elements of L.T x R are close enouph
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-6
|
||||
|
||||
[max_dim_diis_tcscf]
|
||||
type: integer
|
||||
doc: Maximum size of the DIIS extrapolation procedure
|
||||
[var_tc]
|
||||
type: logical
|
||||
doc: If |true|, use VAR-TC
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 15
|
||||
default: False
|
||||
|
||||
[threshold_diis_tcscf]
|
||||
type: Threshold
|
||||
doc: Threshold on the convergence of the DIIS error vector during a TCSCF calculation. If 0. is chosen, the square root of thresh_tcscf will be used.
|
||||
[read_tc_integ]
|
||||
type: logical
|
||||
doc: If |true|, read integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.
|
||||
default: False
|
||||
|
||||
[level_shift_tcscf]
|
||||
type: Positive_float
|
||||
doc: Energy shift on the virtual MOs to improve TCSCF convergence
|
||||
[write_tc_integ]
|
||||
type: logical
|
||||
doc: If |true|, write integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 0.
|
||||
|
||||
[im_thresh_tcscf]
|
||||
type: Threshold
|
||||
doc: Thresholds on the Imag part of energy
|
||||
interface: ezfio,provider,ocaml
|
||||
default: 1.e-7
|
||||
default: False
|
||||
|
||||
[debug_tc_pt2]
|
||||
type: integer
|
||||
|
96
src/tc_scf/diago_vartcfock.irp.f
Normal file
96
src/tc_scf/diago_vartcfock.irp.f
Normal file
@ -0,0 +1,96 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, fock_vartc_eigvec_mo, (mo_num, mo_num)]
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: i, j
|
||||
integer :: liwork, lwork, n, info
|
||||
integer, allocatable :: iwork(:)
|
||||
double precision, allocatable :: work(:), F(:,:), F_save(:,:)
|
||||
double precision, allocatable :: diag(:)
|
||||
|
||||
PROVIDE mo_r_coef
|
||||
PROVIDE Fock_matrix_vartc_mo_tot
|
||||
|
||||
allocate( F(mo_num,mo_num), F_save(mo_num,mo_num) )
|
||||
allocate (diag(mo_num) )
|
||||
|
||||
do j = 1, mo_num
|
||||
do i = 1, mo_num
|
||||
F(i,j) = Fock_matrix_vartc_mo_tot(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Insert level shift here
|
||||
do i = elec_beta_num+1, elec_alpha_num
|
||||
F(i,i) += 0.5d0 * level_shift_tcscf
|
||||
enddo
|
||||
do i = elec_alpha_num+1, mo_num
|
||||
F(i,i) += level_shift_tcscf
|
||||
enddo
|
||||
|
||||
n = mo_num
|
||||
lwork = 1+6*n + 2*n*n
|
||||
liwork = 3 + 5*n
|
||||
|
||||
allocate(work(lwork))
|
||||
allocate(iwork(liwork) )
|
||||
|
||||
lwork = -1
|
||||
liwork = -1
|
||||
|
||||
F_save = F
|
||||
call dsyevd('V', 'U', mo_num, F, size(F, 1), diag, work, lwork, iwork, liwork, info)
|
||||
|
||||
if (info /= 0) then
|
||||
print *, irp_here//' DSYEVD failed : ', info
|
||||
stop 1
|
||||
endif
|
||||
lwork = int(work(1))
|
||||
liwork = iwork(1)
|
||||
deallocate(iwork)
|
||||
deallocate(work)
|
||||
|
||||
allocate(work(lwork))
|
||||
allocate(iwork(liwork) )
|
||||
call dsyevd('V', 'U', mo_num, F, size(F, 1), diag, work, lwork, iwork, liwork, info)
|
||||
deallocate(iwork)
|
||||
|
||||
if (info /= 0) then
|
||||
F = F_save
|
||||
call dsyev('V', 'L', mo_num, F, size(F, 1), diag, work, lwork, info)
|
||||
|
||||
if (info /= 0) then
|
||||
print *, irp_here//' DSYEV failed : ', info
|
||||
stop 1
|
||||
endif
|
||||
endif
|
||||
|
||||
do i = 1, mo_num
|
||||
do j = 1, mo_num
|
||||
fock_vartc_eigvec_mo(j,i) = F(j,i)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
deallocate(work, F, F_save, diag)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, fock_vartc_eigvec_ao, (ao_num, mo_num)]
|
||||
|
||||
implicit none
|
||||
|
||||
PROVIDE mo_r_coef
|
||||
|
||||
call dgemm( 'N', 'N', ao_num, mo_num, mo_num, 1.d0 &
|
||||
, mo_r_coef, size(mo_r_coef, 1), fock_vartc_eigvec_mo, size(fock_vartc_eigvec_mo, 1) &
|
||||
, 0.d0, fock_vartc_eigvec_ao, size(fock_vartc_eigvec_ao, 1))
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
@ -1,17 +1,3 @@
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, threshold_DIIS_nonzero_TCSCF ]
|
||||
|
||||
implicit none
|
||||
|
||||
if(threshold_DIIS_TCSCF == 0.d0) then
|
||||
threshold_DIIS_nonzero_TCSCF = dsqrt(thresh_tcscf)
|
||||
else
|
||||
threshold_DIIS_nonzero_TCSCF = threshold_DIIS_TCSCF
|
||||
endif
|
||||
ASSERT(threshold_DIIS_nonzero_TCSCF >= 0.d0)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
@ -100,13 +86,30 @@ END_PROVIDER
|
||||
BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
|
||||
|
||||
implicit none
|
||||
integer :: i, j
|
||||
double precision, allocatable :: tmp(:,:)
|
||||
double precision, allocatable :: F(:,:)
|
||||
|
||||
allocate(F(ao_num,ao_num))
|
||||
if(var_tc) then
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
F(j,i) = Fock_matrix_vartc_ao_tot(j,i)
|
||||
enddo
|
||||
enddo
|
||||
else
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
F(j,i) = Fock_matrix_tc_ao_tot(j,i)
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
allocate(tmp(ao_num,ao_num))
|
||||
|
||||
! F x Q
|
||||
call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 &
|
||||
, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), Q_matrix, size(Q_matrix, 1) &
|
||||
call dgemm( 'N', 'N', ao_num, ao_num, ao_num, 1.d0 &
|
||||
, F, size(F, 1), Q_matrix, size(Q_matrix, 1) &
|
||||
, 0.d0, tmp, size(tmp, 1) )
|
||||
|
||||
! F x Q x S
|
||||
@ -121,11 +124,12 @@ BEGIN_PROVIDER [double precision, FQS_SQF_ao, (ao_num, ao_num)]
|
||||
, 0.d0, tmp, size(tmp, 1) )
|
||||
|
||||
! F x Q x S - S x Q x F
|
||||
call dgemm( 'N', 'N', ao_num, ao_num, ao_num, -1.d0 &
|
||||
, tmp, size(tmp, 1), Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) &
|
||||
call dgemm( 'N', 'N', ao_num, ao_num, ao_num, -1.d0 &
|
||||
, tmp, size(tmp, 1), F, size(F, 1) &
|
||||
, 1.d0, FQS_SQF_ao, size(FQS_SQF_ao, 1) )
|
||||
|
||||
deallocate(tmp)
|
||||
deallocate(F)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
@ -18,6 +18,8 @@
|
||||
double precision :: density, density_a, density_b
|
||||
double precision :: t0, t1
|
||||
|
||||
!print*, ' providing two_e_tc_non_hermit_integral_seq ...'
|
||||
!call wall_time(t0)
|
||||
|
||||
two_e_tc_non_hermit_integral_seq_alpha = 0.d0
|
||||
two_e_tc_non_hermit_integral_seq_beta = 0.d0
|
||||
@ -31,6 +33,15 @@
|
||||
density_b = TCSCF_density_matrix_ao_beta (l,j)
|
||||
density = density_a + density_b
|
||||
|
||||
!! rho(l,j) * < k l| T | i j>
|
||||
!two_e_tc_non_hermit_integral_seq_alpha(k,i) += density * ao_two_e_tc_tot(l,j,k,i)
|
||||
!! rho(l,j) * < k l| T | i j>
|
||||
!two_e_tc_non_hermit_integral_seq_beta (k,i) += density * ao_two_e_tc_tot(l,j,k,i)
|
||||
!! rho_a(l,j) * < l k| T | i j>
|
||||
!two_e_tc_non_hermit_integral_seq_alpha(k,i) -= density_a * ao_two_e_tc_tot(k,j,l,i)
|
||||
!! rho_b(l,j) * < l k| T | i j>
|
||||
!two_e_tc_non_hermit_integral_seq_beta (k,i) -= density_b * ao_two_e_tc_tot(k,j,l,i)
|
||||
|
||||
! rho(l,j) * < k l| T | i j>
|
||||
two_e_tc_non_hermit_integral_seq_alpha(k,i) += density * ao_two_e_tc_tot(k,i,l,j)
|
||||
! rho(l,j) * < k l| T | i j>
|
||||
@ -45,6 +56,8 @@
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!call wall_time(t1)
|
||||
!print*, ' wall time for two_e_tc_non_hermit_integral_seq after = ', t1 - t0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -67,6 +80,8 @@ END_PROVIDER
|
||||
double precision :: t0, t1
|
||||
double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
|
||||
|
||||
!print*, ' providing two_e_tc_non_hermit_integral ...'
|
||||
!call wall_time(t0)
|
||||
|
||||
two_e_tc_non_hermit_integral_alpha = 0.d0
|
||||
two_e_tc_non_hermit_integral_beta = 0.d0
|
||||
@ -112,6 +127,8 @@ END_PROVIDER
|
||||
deallocate(tmp_a, tmp_b)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
!call wall_time(t1)
|
||||
!print*, ' wall time for two_e_tc_non_hermit_integral after = ', t1 - t0
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
@ -156,6 +173,14 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_alpha, (mo_num, mo_num) ]
|
||||
|
||||
if(bi_ortho) then
|
||||
|
||||
!allocate(tmp(ao_num,ao_num))
|
||||
!tmp = Fock_matrix_tc_ao_alpha
|
||||
!if(three_body_h_tc) then
|
||||
! tmp += fock_3e_uhf_ao_a
|
||||
!endif
|
||||
!call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1))
|
||||
!deallocate(tmp)
|
||||
|
||||
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
|
||||
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
|
||||
if(three_body_h_tc) then
|
||||
@ -184,6 +209,14 @@ BEGIN_PROVIDER [ double precision, Fock_matrix_tc_mo_beta, (mo_num,mo_num) ]
|
||||
|
||||
if(bi_ortho) then
|
||||
|
||||
!allocate(tmp(ao_num,ao_num))
|
||||
!tmp = Fock_matrix_tc_ao_beta
|
||||
!if(three_body_h_tc) then
|
||||
! tmp += fock_3e_uhf_ao_b
|
||||
!endif
|
||||
!call ao_to_mo_bi_ortho(tmp, size(tmp, 1), Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1))
|
||||
!deallocate(tmp)
|
||||
|
||||
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta, size(Fock_matrix_tc_ao_beta, 1) &
|
||||
, Fock_matrix_tc_mo_beta, size(Fock_matrix_tc_mo_beta, 1) )
|
||||
if(three_body_h_tc) then
|
||||
@ -216,10 +249,6 @@ END_PROVIDER
|
||||
do k = elec_beta_num+1, elec_alpha_num
|
||||
grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i)))
|
||||
grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k)))
|
||||
!grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i))
|
||||
!grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k))
|
||||
!grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i)
|
||||
!grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -227,10 +256,6 @@ END_PROVIDER
|
||||
do k = elec_alpha_num+1, mo_num
|
||||
grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i)))
|
||||
grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k)))
|
||||
!grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i))
|
||||
!grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k))
|
||||
grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i)
|
||||
grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
@ -238,15 +263,10 @@ END_PROVIDER
|
||||
do k = elec_alpha_num+1, mo_num
|
||||
grad_non_hermit_left = max(grad_non_hermit_left , dabs(Fock_matrix_tc_mo_tot(k,i)))
|
||||
grad_non_hermit_right = max(grad_non_hermit_right, dabs(Fock_matrix_tc_mo_tot(i,k)))
|
||||
!grad_non_hermit_left += dabs(Fock_matrix_tc_mo_tot(k,i))
|
||||
!grad_non_hermit_right += dabs(Fock_matrix_tc_mo_tot(i,k))
|
||||
grad_non_hermit_left += Fock_matrix_tc_mo_tot(k,i) * Fock_matrix_tc_mo_tot(k,i)
|
||||
grad_non_hermit_right += Fock_matrix_tc_mo_tot(i,k) * Fock_matrix_tc_mo_tot(i,k)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
!grad_non_hermit = dsqrt(grad_non_hermit_left) + dsqrt(grad_non_hermit_right)
|
||||
grad_non_hermit = grad_non_hermit_left + grad_non_hermit_right
|
||||
grad_non_hermit = max(grad_non_hermit_left, grad_non_hermit_right)
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
|
287
src/tc_scf/fock_vartc.irp.f
Normal file
287
src/tc_scf/fock_vartc.irp.f
Normal file
@ -0,0 +1,287 @@
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, two_e_vartc_integral_alpha, (ao_num, ao_num)]
|
||||
&BEGIN_PROVIDER [ double precision, two_e_vartc_integral_beta , (ao_num, ao_num)]
|
||||
|
||||
implicit none
|
||||
integer :: i, j, k, l
|
||||
double precision :: density, density_a, density_b, I_coul, I_kjli
|
||||
double precision :: t0, t1
|
||||
double precision, allocatable :: tmp_a(:,:), tmp_b(:,:)
|
||||
|
||||
two_e_vartc_integral_alpha = 0.d0
|
||||
two_e_vartc_integral_beta = 0.d0
|
||||
|
||||
!$OMP PARALLEL DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j, k, l, density_a, density_b, density, tmp_a, tmp_b, I_coul, I_kjli) &
|
||||
!$OMP SHARED (ao_num, TCSCF_density_matrix_ao_alpha, TCSCF_density_matrix_ao_beta, ao_two_e_vartc_tot, &
|
||||
!$OMP two_e_vartc_integral_alpha, two_e_vartc_integral_beta)
|
||||
|
||||
allocate(tmp_a(ao_num,ao_num), tmp_b(ao_num,ao_num))
|
||||
tmp_a = 0.d0
|
||||
tmp_b = 0.d0
|
||||
|
||||
!$OMP DO
|
||||
do j = 1, ao_num
|
||||
do l = 1, ao_num
|
||||
density_a = TCSCF_density_matrix_ao_alpha(l,j)
|
||||
density_b = TCSCF_density_matrix_ao_beta (l,j)
|
||||
density = density_a + density_b
|
||||
do i = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
|
||||
I_coul = density * ao_two_e_vartc_tot(k,i,l,j)
|
||||
I_kjli = ao_two_e_vartc_tot(k,j,l,i)
|
||||
|
||||
tmp_a(k,i) += I_coul - density_a * I_kjli
|
||||
tmp_b(k,i) += I_coul - density_b * I_kjli
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO NOWAIT
|
||||
|
||||
!$OMP CRITICAL
|
||||
do i = 1, ao_num
|
||||
do j = 1, ao_num
|
||||
two_e_vartc_integral_alpha(j,i) += tmp_a(j,i)
|
||||
two_e_vartc_integral_beta (j,i) += tmp_b(j,i)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END CRITICAL
|
||||
|
||||
deallocate(tmp_a, tmp_b)
|
||||
!$OMP END PARALLEL
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_alpha, (ao_num, ao_num)]
|
||||
|
||||
implicit none
|
||||
|
||||
Fock_matrix_vartc_ao_alpha = ao_one_e_integrals_tc_tot + two_e_vartc_integral_alpha
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_beta, (ao_num, ao_num)]
|
||||
|
||||
implicit none
|
||||
|
||||
Fock_matrix_vartc_ao_beta = ao_one_e_integrals_tc_tot + two_e_vartc_integral_beta
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_alpha, (mo_num, mo_num) ]
|
||||
|
||||
implicit none
|
||||
|
||||
call ao_to_mo_bi_ortho( Fock_matrix_vartc_ao_alpha, size(Fock_matrix_vartc_ao_alpha, 1) &
|
||||
, Fock_matrix_vartc_mo_alpha, size(Fock_matrix_vartc_mo_alpha, 1) )
|
||||
if(three_body_h_tc) then
|
||||
Fock_matrix_vartc_mo_alpha += fock_3e_uhf_mo_a
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_beta, (mo_num,mo_num) ]
|
||||
|
||||
implicit none
|
||||
|
||||
call ao_to_mo_bi_ortho( Fock_matrix_vartc_ao_beta, size(Fock_matrix_vartc_ao_beta, 1) &
|
||||
, Fock_matrix_vartc_mo_beta, size(Fock_matrix_vartc_mo_beta, 1) )
|
||||
if(three_body_h_tc) then
|
||||
Fock_matrix_vartc_mo_beta += fock_3e_uhf_mo_b
|
||||
endif
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, grad_vartc]
|
||||
|
||||
implicit none
|
||||
integer :: i, k
|
||||
double precision :: grad_left, grad_right
|
||||
|
||||
grad_left = 0.d0
|
||||
grad_right = 0.d0
|
||||
|
||||
do i = 1, elec_beta_num ! doc --> SOMO
|
||||
do k = elec_beta_num+1, elec_alpha_num
|
||||
grad_left = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i)))
|
||||
grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k)))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i = 1, elec_beta_num ! doc --> virt
|
||||
do k = elec_alpha_num+1, mo_num
|
||||
grad_left = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i)))
|
||||
grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k)))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do i = elec_beta_num+1, elec_alpha_num ! SOMO --> virt
|
||||
do k = elec_alpha_num+1, mo_num
|
||||
grad_left = max(grad_left , dabs(Fock_matrix_vartc_mo_tot(k,i)))
|
||||
grad_right = max(grad_right, dabs(Fock_matrix_vartc_mo_tot(i,k)))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
grad_vartc = grad_left + grad_right
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_ao_tot, (ao_num, ao_num) ]
|
||||
|
||||
implicit none
|
||||
|
||||
call mo_to_ao_bi_ortho( Fock_matrix_vartc_mo_tot, size(Fock_matrix_vartc_mo_tot, 1) &
|
||||
, Fock_matrix_vartc_ao_tot, size(Fock_matrix_vartc_ao_tot, 1) )
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_mo_tot, (mo_num,mo_num) ]
|
||||
&BEGIN_PROVIDER [ double precision, Fock_matrix_vartc_diag_mo_tot, (mo_num)]
|
||||
|
||||
implicit none
|
||||
integer :: i, j, n
|
||||
|
||||
if(elec_alpha_num == elec_beta_num) then
|
||||
Fock_matrix_vartc_mo_tot = Fock_matrix_vartc_mo_alpha
|
||||
else
|
||||
|
||||
do j = 1, elec_beta_num
|
||||
! F-K
|
||||
do i = 1, elec_beta_num !CC
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
|
||||
- (Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
|
||||
enddo
|
||||
! F+K/2
|
||||
do i = elec_beta_num+1, elec_alpha_num !CA
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
|
||||
+ 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
|
||||
enddo
|
||||
! F
|
||||
do i = elec_alpha_num+1, mo_num !CV
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do j = elec_beta_num+1, elec_alpha_num
|
||||
! F+K/2
|
||||
do i = 1, elec_beta_num !AC
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
|
||||
+ 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
|
||||
enddo
|
||||
! F
|
||||
do i = elec_beta_num+1, elec_alpha_num !AA
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))
|
||||
enddo
|
||||
! F-K/2
|
||||
do i = elec_alpha_num+1, mo_num !AV
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
|
||||
- 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
|
||||
enddo
|
||||
enddo
|
||||
|
||||
do j = elec_alpha_num+1, mo_num
|
||||
! F
|
||||
do i = 1, elec_beta_num !VC
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))
|
||||
enddo
|
||||
! F-K/2
|
||||
do i = elec_beta_num+1, elec_alpha_num !VA
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j))&
|
||||
- 0.5d0*(Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
|
||||
enddo
|
||||
! F+K
|
||||
do i = elec_alpha_num+1, mo_num !VV
|
||||
Fock_matrix_vartc_mo_tot(i,j) = 0.5d0*(Fock_matrix_vartc_mo_alpha(i,j)+Fock_matrix_vartc_mo_beta(i,j)) &
|
||||
+ (Fock_matrix_vartc_mo_beta(i,j) - Fock_matrix_vartc_mo_alpha(i,j))
|
||||
enddo
|
||||
enddo
|
||||
if(three_body_h_tc)then
|
||||
! C-O
|
||||
do j = 1, elec_beta_num
|
||||
do i = elec_beta_num+1, elec_alpha_num
|
||||
Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
|
||||
Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
|
||||
enddo
|
||||
enddo
|
||||
! C-V
|
||||
do j = 1, elec_beta_num
|
||||
do i = elec_alpha_num+1, mo_num
|
||||
Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
|
||||
Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
|
||||
enddo
|
||||
enddo
|
||||
! O-V
|
||||
do j = elec_beta_num+1, elec_alpha_num
|
||||
do i = elec_alpha_num+1, mo_num
|
||||
Fock_matrix_vartc_mo_tot(i,j) += 0.5d0*(fock_a_tot_3e_bi_orth(i,j) + fock_b_tot_3e_bi_orth(i,j))
|
||||
Fock_matrix_vartc_mo_tot(j,i) += 0.5d0*(fock_a_tot_3e_bi_orth(j,i) + fock_b_tot_3e_bi_orth(j,i))
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
endif
|
||||
|
||||
do i = 1, mo_num
|
||||
Fock_matrix_vartc_diag_mo_tot(i) = Fock_matrix_vartc_mo_tot(i,i)
|
||||
enddo
|
||||
|
||||
if(frozen_orb_scf)then
|
||||
integer :: iorb, jorb
|
||||
do i = 1, n_core_orb
|
||||
iorb = list_core(i)
|
||||
do j = 1, n_act_orb
|
||||
jorb = list_act(j)
|
||||
Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
|
||||
Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
if(no_oa_or_av_opt)then
|
||||
do i = 1, n_act_orb
|
||||
iorb = list_act(i)
|
||||
do j = 1, n_inact_orb
|
||||
jorb = list_inact(j)
|
||||
Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
|
||||
Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0
|
||||
enddo
|
||||
do j = 1, n_virt_orb
|
||||
jorb = list_virt(j)
|
||||
Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
|
||||
Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0
|
||||
enddo
|
||||
do j = 1, n_core_orb
|
||||
jorb = list_core(j)
|
||||
Fock_matrix_vartc_mo_tot(iorb,jorb) = 0.d0
|
||||
Fock_matrix_vartc_mo_tot(jorb,iorb) = 0.d0
|
||||
enddo
|
||||
enddo
|
||||
endif
|
||||
|
||||
!call check_sym(Fock_matrix_vartc_mo_tot, mo_num)
|
||||
!do i = 1, mo_num
|
||||
! write(*,'(100(F15.8, I4))') Fock_matrix_vartc_mo_tot(i,:)
|
||||
!enddo
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
@ -1,336 +0,0 @@
|
||||
! ---
|
||||
|
||||
subroutine rh_tcscf()
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Roothaan-Hall algorithm for TC-SCF calculation
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer :: i, j
|
||||
integer :: iteration_TCSCF, dim_DIIS, index_dim_DIIS
|
||||
double precision :: energy_TCSCF, energy_TCSCF_1e, energy_TCSCF_2e, energy_TCSCF_3e, gradie_TCSCF
|
||||
double precision :: energy_TCSCF_previous, delta_energy_TCSCF
|
||||
double precision :: gradie_TCSCF_previous, delta_gradie_TCSCF
|
||||
double precision :: max_error_DIIS_TCSCF
|
||||
double precision :: level_shift_save
|
||||
double precision :: delta_energy_tmp, delta_gradie_tmp
|
||||
double precision, allocatable :: F_DIIS(:,:,:), e_DIIS(:,:,:)
|
||||
double precision, allocatable :: mo_r_coef_save(:,:), mo_l_coef_save(:,:)
|
||||
|
||||
logical, external :: qp_stop
|
||||
|
||||
|
||||
!PROVIDE ao_md5 mo_occ
|
||||
PROVIDE level_shift_TCSCF
|
||||
|
||||
allocate( mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num) &
|
||||
, F_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF), e_DIIS(ao_num,ao_num,max_dim_DIIS_TCSCF) )
|
||||
|
||||
F_DIIS = 0.d0
|
||||
e_DIIS = 0.d0
|
||||
mo_l_coef_save = 0.d0
|
||||
mo_r_coef_save = 0.d0
|
||||
|
||||
call write_time(6)
|
||||
|
||||
! ---
|
||||
! Initialize energies and density matrices
|
||||
|
||||
energy_TCSCF_previous = TC_HF_energy
|
||||
energy_TCSCF_1e = TC_HF_one_e_energy
|
||||
energy_TCSCF_2e = TC_HF_two_e_energy
|
||||
energy_TCSCF_3e = 0.d0
|
||||
if(three_body_h_tc) then
|
||||
energy_TCSCF_3e = diag_three_elem_hf
|
||||
endif
|
||||
gradie_TCSCF_previous = grad_non_hermit
|
||||
delta_energy_TCSCF = 1.d0
|
||||
delta_gradie_TCSCF = 1.d0
|
||||
iteration_TCSCF = 0
|
||||
dim_DIIS = 0
|
||||
max_error_DIIS_TCSCF = 1.d0
|
||||
|
||||
! ---
|
||||
|
||||
! Start of main SCF loop
|
||||
|
||||
PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot
|
||||
|
||||
do while( (max_error_DIIS_TCSCF > threshold_DIIS_nonzero_TCSCF) .or. &
|
||||
!(dabs(delta_energy_TCSCF) > thresh_TCSCF) .or. &
|
||||
(dabs(gradie_TCSCF_previous) > dsqrt(thresh_TCSCF)) )
|
||||
|
||||
iteration_TCSCF += 1
|
||||
if(iteration_TCSCF > n_it_TCSCF_max) then
|
||||
print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
|
||||
stop
|
||||
endif
|
||||
|
||||
dim_DIIS = min(dim_DIIS+1, max_dim_DIIS_TCSCF)
|
||||
|
||||
! ---
|
||||
|
||||
if((tcscf_algorithm == 'DIIS') .and. (dabs(delta_energy_TCSCF) > 1.d-6)) then
|
||||
|
||||
! store Fock and error matrices at each iteration
|
||||
index_dim_DIIS = mod(dim_DIIS-1, max_dim_DIIS_TCSCF) + 1
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
F_DIIS(i,j,index_dim_DIIS) = Fock_matrix_tc_ao_tot(i,j)
|
||||
e_DIIS(i,j,index_dim_DIIS) = FQS_SQF_ao(i,j)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
call extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1), iteration_TCSCF, dim_DIIS)
|
||||
|
||||
Fock_matrix_tc_ao_alpha = 0.5d0 * Fock_matrix_tc_ao_tot
|
||||
Fock_matrix_tc_ao_beta = 0.5d0 * Fock_matrix_tc_ao_tot
|
||||
!TOUCH Fock_matrix_tc_ao_alpha Fock_matrix_tc_ao_beta
|
||||
|
||||
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_alpha, size(Fock_matrix_tc_ao_alpha, 1) &
|
||||
, Fock_matrix_tc_mo_alpha, size(Fock_matrix_tc_mo_alpha, 1) )
|
||||
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_beta , size(Fock_matrix_tc_ao_beta , 1) &
|
||||
, Fock_matrix_tc_mo_beta , size(Fock_matrix_tc_mo_beta , 1) )
|
||||
TOUCH Fock_matrix_tc_mo_alpha Fock_matrix_tc_mo_beta
|
||||
endif
|
||||
|
||||
! ---
|
||||
|
||||
mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num)
|
||||
mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num)
|
||||
TOUCH mo_l_coef mo_r_coef
|
||||
|
||||
! ---
|
||||
|
||||
! calculate error vectors
|
||||
max_error_DIIS_TCSCF = maxval(abs(FQS_SQF_mo))
|
||||
|
||||
! ---
|
||||
|
||||
delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous
|
||||
delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous
|
||||
|
||||
! ---
|
||||
|
||||
do while((delta_gradie_tmp > 1.d-7) .and. (iteration_TCSCF > 1))
|
||||
!do while((dabs(delta_energy_tmp) > 0.5d0) .and. (iteration_TCSCF > 1))
|
||||
print *, ' very big or bad step : ', delta_energy_tmp, delta_gradie_tmp
|
||||
print *, ' TC level shift = ', level_shift_TCSCF
|
||||
|
||||
mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num)
|
||||
mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num)
|
||||
|
||||
if(level_shift_TCSCF <= .1d0) then
|
||||
level_shift_TCSCF = 1.d0
|
||||
else
|
||||
level_shift_TCSCF = level_shift_TCSCF * 3.0d0
|
||||
endif
|
||||
TOUCH mo_l_coef mo_r_coef level_shift_TCSCF
|
||||
|
||||
mo_l_coef(1:ao_num,1:mo_num) = fock_tc_leigvec_ao(1:ao_num,1:mo_num)
|
||||
mo_r_coef(1:ao_num,1:mo_num) = fock_tc_reigvec_ao(1:ao_num,1:mo_num)
|
||||
TOUCH mo_l_coef mo_r_coef
|
||||
|
||||
delta_energy_tmp = TC_HF_energy - energy_TCSCF_previous
|
||||
delta_gradie_tmp = grad_non_hermit - gradie_TCSCF_previous
|
||||
|
||||
if(level_shift_TCSCF - level_shift_save > 40.d0) then
|
||||
level_shift_TCSCF = level_shift_save * 4.d0
|
||||
SOFT_TOUCH level_shift_TCSCF
|
||||
exit
|
||||
endif
|
||||
|
||||
dim_DIIS = 0
|
||||
enddo
|
||||
! print *, ' very big step : ', delta_energy_tmp
|
||||
! print *, ' TC level shift = ', level_shift_TCSCF
|
||||
|
||||
! ---
|
||||
|
||||
level_shift_TCSCF = 0.d0
|
||||
!level_shift_TCSCF = level_shift_TCSCF * 0.5d0
|
||||
SOFT_TOUCH level_shift_TCSCF
|
||||
|
||||
gradie_TCSCF = grad_non_hermit
|
||||
energy_TCSCF = TC_HF_energy
|
||||
energy_TCSCF_1e = TC_HF_one_e_energy
|
||||
energy_TCSCF_2e = TC_HF_two_e_energy
|
||||
energy_TCSCF_3e = 0.d0
|
||||
if(three_body_h_tc) then
|
||||
energy_TCSCF_3e = diag_three_elem_hf
|
||||
endif
|
||||
delta_energy_TCSCF = energy_TCSCF - energy_TCSCF_previous
|
||||
delta_gradie_TCSCF = gradie_TCSCF - gradie_TCSCF_previous
|
||||
|
||||
energy_TCSCF_previous = energy_TCSCF
|
||||
gradie_TCSCF_previous = gradie_TCSCF
|
||||
|
||||
|
||||
level_shift_save = level_shift_TCSCF
|
||||
mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num)
|
||||
mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num)
|
||||
|
||||
|
||||
print *, ' iteration = ', iteration_TCSCF
|
||||
print *, ' total TC energy = ', energy_TCSCF
|
||||
print *, ' 1-e TC energy = ', energy_TCSCF_1e
|
||||
print *, ' 2-e TC energy = ', energy_TCSCF_2e
|
||||
print *, ' 3-e TC energy = ', energy_TCSCF_3e
|
||||
print *, ' |delta TC energy| = ', dabs(delta_energy_TCSCF)
|
||||
print *, ' TC gradient = ', gradie_TCSCF
|
||||
print *, ' delta TC gradient = ', delta_gradie_TCSCF
|
||||
print *, ' max TC DIIS error = ', max_error_DIIS_TCSCF
|
||||
print *, ' TC DIIS dim = ', dim_DIIS
|
||||
print *, ' TC level shift = ', level_shift_TCSCF
|
||||
print *, ' '
|
||||
|
||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
|
||||
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
||||
|
||||
if(qp_stop()) exit
|
||||
enddo
|
||||
|
||||
! ---
|
||||
|
||||
print *, ' TCSCF DIIS converged !'
|
||||
call print_energy_and_mos()
|
||||
|
||||
call write_time(6)
|
||||
|
||||
deallocate(mo_r_coef_save, mo_l_coef_save, F_DIIS, e_DIIS)
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
||||
subroutine extrapolate_TC_Fock_matrix(e_DIIS, F_DIIS, F_ao, size_F_ao, iteration_TCSCF, dim_DIIS)
|
||||
|
||||
BEGIN_DOC
|
||||
!
|
||||
! Compute the extrapolated Fock matrix using the DIIS procedure
|
||||
!
|
||||
! e = \sum_i c_i e_i and \sum_i c_i = 1
|
||||
! ==> lagrange multiplier with L = |e|^2 - \lambda (\sum_i c_i = 1)
|
||||
!
|
||||
END_DOC
|
||||
|
||||
implicit none
|
||||
|
||||
integer, intent(in) :: iteration_TCSCF, size_F_ao
|
||||
integer, intent(inout) :: dim_DIIS
|
||||
double precision, intent(in) :: F_DIIS(ao_num,ao_num,dim_DIIS)
|
||||
double precision, intent(in) :: e_DIIS(ao_num,ao_num,dim_DIIS)
|
||||
double precision, intent(inout) :: F_ao(size_F_ao,ao_num)
|
||||
|
||||
double precision, allocatable :: B_matrix_DIIS(:,:), X_vector_DIIS(:), C_vector_DIIS(:)
|
||||
|
||||
integer :: i, j, k, l, i_DIIS, j_DIIS
|
||||
integer :: lwork
|
||||
double precision :: rcond, ferr, berr
|
||||
integer, allocatable :: iwork(:)
|
||||
double precision, allocatable :: scratch(:,:)
|
||||
|
||||
if(dim_DIIS < 1) then
|
||||
return
|
||||
endif
|
||||
|
||||
allocate( B_matrix_DIIS(dim_DIIS+1,dim_DIIS+1), X_vector_DIIS(dim_DIIS+1) &
|
||||
, C_vector_DIIS(dim_DIIS+1), scratch(ao_num,ao_num) )
|
||||
|
||||
! Compute the matrices B and X
|
||||
B_matrix_DIIS(:,:) = 0.d0
|
||||
do j = 1, dim_DIIS
|
||||
j_DIIS = min(dim_DIIS, mod(iteration_TCSCF-j, max_dim_DIIS_TCSCF)+1)
|
||||
|
||||
do i = 1, dim_DIIS
|
||||
i_DIIS = min(dim_DIIS, mod(iteration_TCSCF-i, max_dim_DIIS_TCSCF)+1)
|
||||
|
||||
! Compute product of two errors vectors
|
||||
do l = 1, ao_num
|
||||
do k = 1, ao_num
|
||||
B_matrix_DIIS(i,j) = B_matrix_DIIS(i,j) + e_DIIS(k,l,i_DIIS) * e_DIIS(k,l,j_DIIS)
|
||||
enddo
|
||||
enddo
|
||||
|
||||
enddo
|
||||
enddo
|
||||
|
||||
! Pad B matrix and build the X matrix
|
||||
|
||||
C_vector_DIIS(:) = 0.d0
|
||||
do i = 1, dim_DIIS
|
||||
B_matrix_DIIS(i,dim_DIIS+1) = -1.d0
|
||||
B_matrix_DIIS(dim_DIIS+1,i) = -1.d0
|
||||
enddo
|
||||
C_vector_DIIS(dim_DIIS+1) = -1.d0
|
||||
|
||||
deallocate(scratch)
|
||||
|
||||
! Estimate condition number of B
|
||||
integer :: info
|
||||
double precision :: anorm
|
||||
integer, allocatable :: ipiv(:)
|
||||
double precision, allocatable :: AF(:,:)
|
||||
double precision, external :: dlange
|
||||
|
||||
lwork = max((dim_DIIS+1)**2, (dim_DIIS+1)*5)
|
||||
allocate(AF(dim_DIIS+1,dim_DIIS+1))
|
||||
allocate(ipiv(2*(dim_DIIS+1)), iwork(2*(dim_DIIS+1)) )
|
||||
allocate(scratch(lwork,1))
|
||||
scratch(:,1) = 0.d0
|
||||
|
||||
anorm = dlange('1', dim_DIIS+1, dim_DIIS+1, B_matrix_DIIS, size(B_matrix_DIIS, 1), scratch(1,1))
|
||||
|
||||
AF(:,:) = B_matrix_DIIS(:,:)
|
||||
call dgetrf(dim_DIIS+1, dim_DIIS+1, AF, size(AF, 1), ipiv, info)
|
||||
if(info /= 0) then
|
||||
dim_DIIS = 0
|
||||
return
|
||||
endif
|
||||
|
||||
call dgecon('1', dim_DIIS+1, AF, size(AF, 1), anorm, rcond, scratch, iwork, info)
|
||||
if(info /= 0) then
|
||||
dim_DIIS = 0
|
||||
return
|
||||
endif
|
||||
|
||||
if(rcond < 1.d-14) then
|
||||
dim_DIIS = 0
|
||||
return
|
||||
endif
|
||||
|
||||
! solve the linear system C = B x X
|
||||
|
||||
X_vector_DIIS = C_vector_DIIS
|
||||
call dgesv(dim_DIIS+1, 1, B_matrix_DIIS, size(B_matrix_DIIS, 1), ipiv , X_vector_DIIS, size(X_vector_DIIS, 1), info)
|
||||
|
||||
deallocate(scratch, AF, iwork)
|
||||
if(info < 0) then
|
||||
stop ' bug in TC-DIIS'
|
||||
endif
|
||||
|
||||
! Compute extrapolated Fock matrix
|
||||
|
||||
!$OMP PARALLEL DO PRIVATE(i,j,k) DEFAULT(SHARED) if (ao_num > 200)
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
F_ao(i,j) = 0.d0
|
||||
enddo
|
||||
do k = 1, dim_DIIS
|
||||
if(dabs(X_vector_DIIS(k)) < 1.d-10) cycle
|
||||
do i = 1,ao_num
|
||||
! FPE here
|
||||
F_ao(i,j) = F_ao(i,j) + X_vector_DIIS(k) * F_DIIS(i,j,dim_DIIS-k+1)
|
||||
enddo
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END PARALLEL DO
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
@ -21,7 +21,7 @@ subroutine rh_tcscf_diis()
|
||||
dim_DIIS = 0
|
||||
g_delta_th = 1d0
|
||||
er_delta_th = 1d0
|
||||
rate_th = 100.d0 !0.01d0 !0.2d0
|
||||
rate_th = 0.1d0
|
||||
|
||||
allocate(mo_r_coef_save(ao_num,mo_num), mo_l_coef_save(ao_num,mo_num))
|
||||
mo_l_coef_save = 0.d0
|
||||
@ -38,17 +38,25 @@ subroutine rh_tcscf_diis()
|
||||
PROVIDE level_shift_TCSCF
|
||||
PROVIDE mo_l_coef mo_r_coef
|
||||
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
'====', '================', '================', '================', '================', '================' &
|
||||
, '================', '================', '================', '====', '========'
|
||||
!write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
! '====', '================', '================', '================', '================', '================' &
|
||||
! , '================', '================', '================', '====', '========'
|
||||
!write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
! ' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' &
|
||||
! , ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)'
|
||||
!write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
! '====', '================', '================', '================', '================', '================' &
|
||||
! , '================', '================', '================', '====', '========'
|
||||
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
'====', '================', '================', '================', '================', '================' &
|
||||
, '================', '================', '====', '========'
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' &
|
||||
, ' gradient ', ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)'
|
||||
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
, ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)'
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
'====', '================', '================', '================', '================', '================' &
|
||||
, '================', '================', '================', '====', '========'
|
||||
, '================', '================', '====', '========'
|
||||
|
||||
|
||||
! first iteration (HF orbitals)
|
||||
@ -61,23 +69,26 @@ subroutine rh_tcscf_diis()
|
||||
if(three_body_h_tc) then
|
||||
etc_3e = diag_three_elem_hf
|
||||
endif
|
||||
tc_grad = grad_non_hermit
|
||||
!tc_grad = grad_non_hermit
|
||||
er_DIIS = maxval(abs(FQS_SQF_mo))
|
||||
e_delta = dabs(etc_tot - e_save)
|
||||
|
||||
e_save = etc_tot
|
||||
g_save = tc_grad
|
||||
!g_save = tc_grad
|
||||
er_save = er_DIIS
|
||||
|
||||
call wall_time(t1)
|
||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
!write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
! it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE FQS_SQF_ao Fock_matrix_tc_ao_tot
|
||||
|
||||
do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. threshold_DIIS_nonzero_TCSCF))
|
||||
!do while((tc_grad .gt. dsqrt(thresh_tcscf)) .and. (er_DIIS .gt. dsqrt(thresh_tcscf)))
|
||||
do while(er_DIIS .gt. dsqrt(thresh_tcscf))
|
||||
|
||||
call wall_time(t0)
|
||||
|
||||
@ -118,12 +129,10 @@ subroutine rh_tcscf_diis()
|
||||
|
||||
! ---
|
||||
|
||||
g_delta = grad_non_hermit - g_save
|
||||
!g_delta = grad_non_hermit - g_save
|
||||
er_delta = maxval(abs(FQS_SQF_mo)) - er_save
|
||||
|
||||
!if((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1)) then
|
||||
if((g_delta > rate_th * g_delta_th) .and. (it > 1)) then
|
||||
!if((g_delta > 0.d0) .and. (it > 1)) then
|
||||
if((er_delta > rate_th * er_save) .and. (it > 1)) then
|
||||
|
||||
Fock_matrix_tc_ao_tot(1:ao_num,1:ao_num) = F_DIIS(1:ao_num,1:ao_num,index_dim_DIIS)
|
||||
call ao_to_mo_bi_ortho( Fock_matrix_tc_ao_tot, size(Fock_matrix_tc_ao_tot, 1) &
|
||||
@ -140,15 +149,16 @@ subroutine rh_tcscf_diis()
|
||||
|
||||
! ---
|
||||
|
||||
g_delta = grad_non_hermit - g_save
|
||||
!g_delta = grad_non_hermit - g_save
|
||||
er_delta = maxval(abs(FQS_SQF_mo)) - er_save
|
||||
|
||||
mo_l_coef_save(1:ao_num,1:mo_num) = mo_l_coef(1:ao_num,1:mo_num)
|
||||
mo_r_coef_save(1:ao_num,1:mo_num) = mo_r_coef(1:ao_num,1:mo_num)
|
||||
|
||||
!do while((g_delta > rate_th * g_delta_th) .and. (er_delta > rate_th * er_delta_th) .and. (it > 1))
|
||||
do while((g_delta > rate_th * g_delta_th) .and. (it > 1))
|
||||
print *, ' big or bad step : ', g_delta, rate_th * g_delta_th
|
||||
do while((er_delta > rate_th * er_save) .and. (it > 1))
|
||||
print *, ' big or bad step '
|
||||
!print *, g_delta , rate_th * g_save
|
||||
print *, er_delta, rate_th * er_save
|
||||
|
||||
mo_l_coef(1:ao_num,1:mo_num) = mo_l_coef_save(1:ao_num,1:mo_num)
|
||||
mo_r_coef(1:ao_num,1:mo_num) = mo_r_coef_save(1:ao_num,1:mo_num)
|
||||
@ -165,7 +175,7 @@ subroutine rh_tcscf_diis()
|
||||
!call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
||||
TOUCH mo_l_coef mo_r_coef
|
||||
|
||||
g_delta = grad_non_hermit - g_save
|
||||
!g_delta = grad_non_hermit - g_save
|
||||
er_delta = maxval(abs(FQS_SQF_mo)) - er_save
|
||||
|
||||
if(level_shift_TCSCF - level_shift_save > 40.d0) then
|
||||
@ -189,25 +199,27 @@ subroutine rh_tcscf_diis()
|
||||
if(three_body_h_tc) then
|
||||
etc_3e = diag_three_elem_hf
|
||||
endif
|
||||
tc_grad = grad_non_hermit
|
||||
!tc_grad = grad_non_hermit
|
||||
er_DIIS = maxval(abs(FQS_SQF_mo))
|
||||
e_delta = dabs(etc_tot - e_save)
|
||||
g_delta = tc_grad - g_save
|
||||
!g_delta = tc_grad - g_save
|
||||
er_delta = er_DIIS - er_save
|
||||
|
||||
e_save = etc_tot
|
||||
g_save = tc_grad
|
||||
!g_save = tc_grad
|
||||
level_shift_save = level_shift_TCSCF
|
||||
er_save = er_DIIS
|
||||
|
||||
g_delta_th = dabs(tc_grad) ! g_delta)
|
||||
!g_delta_th = dabs(tc_grad) ! g_delta)
|
||||
er_delta_th = dabs(er_DIIS) !er_delta)
|
||||
|
||||
call wall_time(t1)
|
||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
!write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
! it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, tc_grad, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
|
||||
if(g_delta .lt. 0.d0) then
|
||||
if(er_delta .lt. 0.d0) then
|
||||
call ezfio_set_tc_scf_bitc_energy(etc_tot)
|
||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
|
||||
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
||||
|
89
src/tc_scf/rh_vartcscf_simple.irp.f
Normal file
89
src/tc_scf/rh_vartcscf_simple.irp.f
Normal file
@ -0,0 +1,89 @@
|
||||
! ---
|
||||
|
||||
subroutine rh_vartcscf_simple()
|
||||
|
||||
implicit none
|
||||
integer :: i, j, it, dim_DIIS
|
||||
double precision :: t0, t1
|
||||
double precision :: e_save, e_delta, rho_delta
|
||||
double precision :: etc_tot, etc_1e, etc_2e, etc_3e
|
||||
double precision :: er_DIIS
|
||||
|
||||
|
||||
it = 0
|
||||
e_save = 0.d0
|
||||
dim_DIIS = 0
|
||||
|
||||
! ---
|
||||
|
||||
PROVIDE level_shift_tcscf
|
||||
PROVIDE mo_r_coef
|
||||
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
'====', '================', '================', '================', '================', '================' &
|
||||
, '================', '================', '====', '========'
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
' it ', ' SCF TC Energy ', ' E(1e) ', ' E(2e) ', ' E(3e) ', ' energy diff ' &
|
||||
, ' DIIS error ', ' level shift ', 'DIIS', ' WT (m)'
|
||||
write(6, '(A4,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A16,1X, A4, 1X, A8)') &
|
||||
'====', '================', '================', '================', '================', '================' &
|
||||
, '================', '================', '====', '========'
|
||||
|
||||
|
||||
! first iteration (HF orbitals)
|
||||
call wall_time(t0)
|
||||
|
||||
etc_tot = VARTC_HF_energy
|
||||
etc_1e = VARTC_HF_one_e_energy
|
||||
etc_2e = VARTC_HF_two_e_energy
|
||||
etc_3e = 0.d0
|
||||
if(three_body_h_tc) then
|
||||
etc_3e = diag_three_elem_hf
|
||||
endif
|
||||
er_DIIS = maxval(abs(FQS_SQF_mo))
|
||||
e_delta = dabs(etc_tot - e_save)
|
||||
e_save = etc_tot
|
||||
|
||||
call wall_time(t1)
|
||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
|
||||
do while(er_DIIS .gt. dsqrt(thresh_tcscf))
|
||||
call wall_time(t0)
|
||||
|
||||
it += 1
|
||||
if(it > n_it_tcscf_max) then
|
||||
print *, ' max of TCSCF iterations is reached ', n_it_TCSCF_max
|
||||
stop
|
||||
endif
|
||||
|
||||
mo_r_coef = fock_vartc_eigvec_ao
|
||||
mo_l_coef = mo_r_coef
|
||||
call ezfio_set_bi_ortho_mos_mo_l_coef(mo_l_coef)
|
||||
call ezfio_set_bi_ortho_mos_mo_r_coef(mo_r_coef)
|
||||
TOUCH mo_l_coef mo_r_coef
|
||||
|
||||
etc_tot = VARTC_HF_energy
|
||||
etc_1e = VARTC_HF_one_e_energy
|
||||
etc_2e = VARTC_HF_two_e_energy
|
||||
etc_3e = 0.d0
|
||||
if(three_body_h_tc) then
|
||||
etc_3e = diag_three_elem_hf
|
||||
endif
|
||||
er_DIIS = maxval(abs(FQS_SQF_mo))
|
||||
e_delta = dabs(etc_tot - e_save)
|
||||
e_save = etc_tot
|
||||
|
||||
call ezfio_set_tc_scf_bitc_energy(etc_tot)
|
||||
|
||||
call wall_time(t1)
|
||||
write(6, '(I4,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, F16.10,1X, I4,1X, F8.2)') &
|
||||
it, etc_tot, etc_1e, etc_2e, etc_3e, e_delta, er_DIIS, level_shift_tcscf, dim_DIIS, (t1-t0)/60.d0
|
||||
enddo
|
||||
|
||||
print *, ' VAR-TCSCF Simple converged !'
|
||||
|
||||
end
|
||||
|
||||
! ---
|
||||
|
@ -73,3 +73,4 @@ subroutine create_guess()
|
||||
end subroutine create_guess
|
||||
|
||||
! ---
|
||||
|
||||
|
@ -30,5 +30,34 @@
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
BEGIN_PROVIDER [ double precision, VARTC_HF_energy]
|
||||
&BEGIN_PROVIDER [ double precision, VARTC_HF_one_e_energy]
|
||||
&BEGIN_PROVIDER [ double precision, VARTC_HF_two_e_energy]
|
||||
|
||||
implicit none
|
||||
integer :: i, j
|
||||
|
||||
PROVIDE mo_r_coef
|
||||
|
||||
VARTC_HF_energy = nuclear_repulsion
|
||||
VARTC_HF_one_e_energy = 0.d0
|
||||
VARTC_HF_two_e_energy = 0.d0
|
||||
|
||||
do j = 1, ao_num
|
||||
do i = 1, ao_num
|
||||
VARTC_HF_two_e_energy += 0.5d0 * ( two_e_vartc_integral_alpha(i,j) * TCSCF_density_matrix_ao_alpha(i,j) &
|
||||
+ two_e_vartc_integral_beta (i,j) * TCSCF_density_matrix_ao_beta (i,j) )
|
||||
VARTC_HF_one_e_energy += ao_one_e_integrals_tc_tot(i,j) &
|
||||
* (TCSCF_density_matrix_ao_alpha(i,j) + TCSCF_density_matrix_ao_beta (i,j) )
|
||||
enddo
|
||||
enddo
|
||||
|
||||
VARTC_HF_energy += VARTC_HF_one_e_energy + VARTC_HF_two_e_energy
|
||||
VARTC_HF_energy += diag_three_elem_hf
|
||||
|
||||
END_PROVIDER
|
||||
|
||||
! ---
|
||||
|
||||
|
1059
src/tc_scf/test_int.irp.f
Normal file
1059
src/tc_scf/test_int.irp.f
Normal file
File diff suppressed because it is too large
Load Diff
@ -459,3 +459,38 @@ subroutine v2_over_x(v,x,res)
|
||||
res = 0.5d0 * (tmp - delta_E)
|
||||
|
||||
end
|
||||
|
||||
subroutine sum_A_At(A, N)
|
||||
|
||||
!BEGIN_DOC
|
||||
! useful for symmetrizing a tensor without a temporary tensor
|
||||
!END_DOC
|
||||
|
||||
implicit none
|
||||
integer, intent(in) :: N
|
||||
double precision, intent(inout) :: A(N,N)
|
||||
integer :: i, j
|
||||
|
||||
!$OMP PARALLEL &
|
||||
!$OMP DEFAULT (NONE) &
|
||||
!$OMP PRIVATE (i, j) &
|
||||
!$OMP SHARED (A, N)
|
||||
!$OMP DO
|
||||
do j = 1, N
|
||||
do i = j, N
|
||||
A(i,j) += A(j,i)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
|
||||
!$OMP DO
|
||||
do j = 2, N
|
||||
do i = 1, j-1
|
||||
A(i,j) = A(j,i)
|
||||
enddo
|
||||
enddo
|
||||
!$OMP END DO
|
||||
!$OMP END PARALLEL
|
||||
|
||||
end
|
||||
|
||||
|
Loading…
Reference in New Issue
Block a user