diff --git a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f index 5e7ef7e9..f3e93360 100644 --- a/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f +++ b/src/non_h_ints_mu/debug_integ_jmu_modif.irp.f @@ -27,12 +27,13 @@ program debug_integ_jmu_modif ! call test_int2_grad1_u12_ao() ! ! call test_grad12_j12() + call test_tchint_rsdft() ! call test_u12sq_j1bsq() ! call test_u12_grad1_u12_j1b_grad1_j1b() ! !call test_gradu_squared_u_ij_mu() !call test_vect_overlap_gauss_r12_ao() - call test_vect_overlap_gauss_r12_ao_with1s() + !call test_vect_overlap_gauss_r12_ao_with1s() end @@ -473,6 +474,65 @@ end subroutine test_gradu_squared_u_ij_mu ! --- +subroutine test_tchint_rsdft() + + implicit none + integer :: i, j, m, ipoint, jpoint + double precision :: acc_ij, acc_tot, eps_ij, i_exc, i_num, normalz + double precision :: x(3), y(3), dj_1(3), dj_2(3), dj_3(3) + + print*, ' test rsdft_jastrow ...' + + PROVIDE grad1_u12_num + + eps_ij = 1d-4 + acc_tot = 0.d0 + normalz = 0.d0 + + do ipoint = 1, n_points_final_grid + x(1) = final_grid_points(1,ipoint) + x(2) = final_grid_points(2,ipoint) + x(3) = final_grid_points(3,ipoint) + + do jpoint = 1, n_points_extra_final_grid + y(1) = final_grid_points_extra(1,jpoint) + y(2) = final_grid_points_extra(2,jpoint) + y(3) = final_grid_points_extra(3,jpoint) + + dj_1(1) = grad1_u12_num(jpoint,ipoint,1) + dj_1(2) = grad1_u12_num(jpoint,ipoint,2) + dj_1(3) = grad1_u12_num(jpoint,ipoint,3) + + call get_tchint_rsdft_jastrow(x, y, dj_2) + + do m = 1, 3 + i_exc = dj_1(m) + i_num = dj_2(m) + acc_ij = dabs(i_exc - i_num) + if(acc_ij .gt. eps_ij) then + print *, ' problem on', ipoint, jpoint, m + print *, ' x = ', x + print *, ' y = ', y + print *, ' exc, num, diff = ', i_exc, i_num, acc_ij + call grad1_jmu_modif_num(x, y, dj_3) + print *, ' check = ', dj_3(m) + stop + endif + + acc_tot += acc_ij + normalz += dabs(i_exc) + enddo + enddo + enddo + + print*, ' acc_tot = ', acc_tot + print*, ' normalz = ', normalz + + return +end subroutine test_tchint_rsdft + +! --- + subroutine test_grad12_j12() implicit none @@ -484,7 +544,7 @@ subroutine test_grad12_j12() PROVIDE grad12_j12 - eps_ij = 1d-3 + eps_ij = 1d-6 acc_tot = 0.d0 normalz = 0.d0 diff --git a/src/non_h_ints_mu/j12_nucl_utils.irp.f b/src/non_h_ints_mu/j12_nucl_utils.irp.f index ac077fe0..7dd13f14 100644 --- a/src/non_h_ints_mu/j12_nucl_utils.irp.f +++ b/src/non_h_ints_mu/j12_nucl_utils.irp.f @@ -668,7 +668,7 @@ subroutine grad1_jmu_modif_num(r1, r2, grad) double precision, intent(in) :: r1(3), r2(3) double precision, intent(out) :: grad(3) - double precision :: tmp0, tmp1, tmp2, tmp3, tmp4, grad_u12(3) + double precision :: tmp0, tmp1, tmp2, grad_u12(3) double precision, external :: j12_mu double precision, external :: j1b_nucl @@ -681,18 +681,93 @@ subroutine grad1_jmu_modif_num(r1, r2, grad) tmp0 = j1b_nucl(r1) tmp1 = j1b_nucl(r2) tmp2 = j12_mu(r1, r2) - tmp3 = tmp0 * tmp1 - tmp4 = tmp2 * tmp1 - grad(1) = tmp3 * grad_u12(1) + tmp4 * grad_x_j1b_nucl_num(r1) - grad(2) = tmp3 * grad_u12(2) + tmp4 * grad_y_j1b_nucl_num(r1) - grad(3) = tmp3 * grad_u12(3) + tmp4 * grad_z_j1b_nucl_num(r1) + grad(1) = (tmp0 * grad_u12(1) + tmp2 * grad_x_j1b_nucl_num(r1)) * tmp1 + grad(2) = (tmp0 * grad_u12(2) + tmp2 * grad_y_j1b_nucl_num(r1)) * tmp1 + grad(3) = (tmp0 * grad_u12(3) + tmp2 * grad_z_j1b_nucl_num(r1)) * tmp1 return end subroutine grad1_jmu_modif_num ! --- - +subroutine get_tchint_rsdft_jastrow(x, y, dj) + + implicit none + double precision, intent(in) :: x(3), y(3) + double precision, intent(out) :: dj(3) + integer :: at + double precision :: a, mu_tmp, inv_sq_pi_2 + double precision :: tmp_x, tmp_y, tmp_z, tmp + double precision :: dx2, dy2, pos(3), dxy, dxy2 + double precision :: v1b_x, v1b_y + double precision :: u2b, grad1_u2b(3), grad1_v1b(3) + + PROVIDE mu_erf + + inv_sq_pi_2 = 0.5d0 / dsqrt(dacos(-1.d0)) + + dj = 0.d0 + +! double precision, external :: j12_mu, j1b_nucl +! v1b_x = j1b_nucl(x) +! v1b_y = j1b_nucl(y) +! call grad1_j1b_nucl(x, grad1_v1b) +! u2b = j12_mu(x, y) +! call grad1_j12_mu(x, y, grad1_u2b) + + ! 1b terms + v1b_x = 1.d0 + v1b_y = 1.d0 + tmp_x = 0.d0 + tmp_y = 0.d0 + tmp_z = 0.d0 + do at = 1, nucl_num + + a = j1b_pen(at) + pos(1) = nucl_coord(at,1) + pos(2) = nucl_coord(at,2) + pos(3) = nucl_coord(at,3) + + dx2 = sum((x-pos)**2) + dy2 = sum((y-pos)**2) + tmp = dexp(-a*dx2) * a + + v1b_x = v1b_x - dexp(-a*dx2) + v1b_y = v1b_y - dexp(-a*dy2) + + tmp_x = tmp_x + tmp * (x(1) - pos(1)) + tmp_y = tmp_y + tmp * (x(2) - pos(2)) + tmp_z = tmp_z + tmp * (x(3) - pos(3)) + end do + grad1_v1b(1) = 2.d0 * tmp_x + grad1_v1b(2) = 2.d0 * tmp_y + grad1_v1b(3) = 2.d0 * tmp_z + + ! 2b terms + dxy2 = sum((x-y)**2) + dxy = dsqrt(dxy2) + mu_tmp = mu_erf * dxy + u2b = 0.5d0 * dxy * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf + + if(dxy .lt. 1d-8) then + grad1_u2b(1) = 0.d0 + grad1_u2b(2) = 0.d0 + grad1_u2b(3) = 0.d0 + else + tmp = 0.5d0 * (1.d0 - derf(mu_tmp)) / dxy + grad1_u2b(1) = tmp * (x(1) - y(1)) + grad1_u2b(2) = tmp * (x(2) - y(2)) + grad1_u2b(3) = tmp * (x(3) - y(3)) + endif + + dj(1) = (grad1_u2b(1) * v1b_x + u2b * grad1_v1b(1)) * v1b_y + dj(2) = (grad1_u2b(2) * v1b_x + u2b * grad1_v1b(2)) * v1b_y + dj(3) = (grad1_u2b(3) * v1b_x + u2b * grad1_v1b(3)) * v1b_y + + return +end subroutine get_tchint_rsdft_jastrow + +! --- diff --git a/src/non_h_ints_mu/jast_deriv.irp.f b/src/non_h_ints_mu/jast_deriv.irp.f index bd7ff6b7..859f2aa5 100644 --- a/src/non_h_ints_mu/jast_deriv.irp.f +++ b/src/non_h_ints_mu/jast_deriv.irp.f @@ -70,6 +70,8 @@ elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then + PROVIDE final_grid_points + !$OMP PARALLEL & !$OMP DEFAULT (NONE) & !$OMP PRIVATE (ipoint, jpoint, r1, r2, v1b_r1, v1b_r2, u2b_r12, grad1_v1b, grad1_u2b, dx, dy, dz) & diff --git a/src/tc_bi_ortho/print_tc_dump.irp.f b/src/tc_bi_ortho/print_tc_dump.irp.f index 55df20a2..0a7e08d2 100644 --- a/src/tc_bi_ortho/print_tc_dump.irp.f +++ b/src/tc_bi_ortho/print_tc_dump.irp.f @@ -8,6 +8,8 @@ program tc_bi_ortho my_grid_becke = .True. my_n_pt_r_grid = 30 my_n_pt_a_grid = 50 + !my_n_pt_r_grid = 100 + !my_n_pt_a_grid = 170 touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid call ERI_dump() @@ -21,26 +23,66 @@ end subroutine KMat_tilde_dump() implicit none - integer :: i, j, k, l + integer :: i, j, k, l + integer :: isym, ms2, st, iii + character(16) :: corb + double precision :: t1, t2 + integer, allocatable :: orbsym(:) + + print *, ' generating FCIDUMP' + call wall_time(t1) PROVIDE mo_bi_ortho_tc_two_e_chemist + PROVIDE mo_bi_ortho_tc_one_e - print *, ' Kmat_tilde in chem notation' + isym = 1 + ms2 = elec_alpha_num - elec_beta_num + st = 0 + iii = 0 + + allocate(orbsym(mo_num)) + orbsym(1:mo_num) = 1 + + open(33, file='FCIDUMP', action='write') + + write(33,'("&",a)') 'FCI' + write(33,'(1x,a,"=",i0,",")') 'NORB', mo_num + write(33,'(1x,a,"=",i0,",")') 'NELEC', elec_num + write(33,'(1x,a,"=",i0,",")') 'MS2', ms2 + write(33,'(1x,a,"=",i0,",")') 'ISYM', isym + write(corb,'(i0)') mo_num + write(33,'(1x,a,"=",'//corb//'(i0,","))') 'ORBSYM', orbsym + write(33,'(1x,a,"=",i0,",")') 'ST', st + write(33,'(1x,a,"=",i0,",")') 'III', iii + write(33,'(1x,a,"=",i0,",")') 'OCC', (elec_num-ms2)/2+ms2 + write(33,'(1x,a,"=",i0,",")') 'CLOSED', 2*elec_alpha_num + write(33,'(1x,"/")') - open(33, file='Kmat_tilde.dat', action='write') do l = 1, mo_num do k = 1, mo_num do j = 1, mo_num do i = 1, mo_num - write(33, '(4(I4, 2X), 4X, E15.7)') i, j, k, l, mo_bi_ortho_tc_two_e_chemist(i,j,k,l) ! TCHint convention - !write(33, '(4(I4, 2X), 4X, E15.7)') i, j, k, l, mo_bi_ortho_tc_two_e_chemist(j,i,l,k) + write(33, '(E15.7, 4X, 4(I4, 2X))') mo_bi_ortho_tc_two_e_chemist(j,i,l,k), i, j, k, l enddo enddo enddo enddo + + do j = 1, mo_num + do i = 1, mo_num + ! TCHint convention + write(33, '(E15.7, 4X, 4(I4, 2X))') mo_bi_ortho_tc_one_e(i,j), i, j, 0, 0 + enddo + enddo + close(33) + deallocate(orbsym) + + call wall_time(t2) + print *, ' end after (min)', (t2-t1)/60.d0 + return end subroutine KMat_tilde_dump @@ -106,12 +148,15 @@ subroutine LMat_tilde_dump() implicit none integer :: i, j, k, l, m, n double precision :: integral + double precision :: t1, t2 - PROVIDE mo_bi_ortho_tc_two_e_chemist + print *, ' generating TCDUMP' + call wall_time(t1) - print *, ' Lmat_tilde in phys notation' + PROVIDE mo_l_coef mo_r_coef - open(33, file='Lmat_tilde.dat', action='write') + open(33, file='TCDUMP', action='write') + write(33, '(4X, I4)') mo_num do n = 1, mo_num do m = 1, mo_num do l = 1, mo_num @@ -120,7 +165,12 @@ subroutine LMat_tilde_dump() do i = 1, mo_num ! < i j k | -L | l m n > with a BI-ORTHONORMAL MOLECULAR ORBITALS call give_integrals_3_body_bi_ort(i, j, k, l, m, n, integral) - write(33, '(6(I4, 2X), 4X, E15.7)') i, j, k, l, m, n, integral + !write(33, '(6(I4, 2X), 4X, E15.7)') i, j, k, l, m, n, integral + ! TCHint convention + if(dabs(integral).gt.1d-10) then + write(33, '(E15.7, 4X, 6(I4, 2X))') -integral/3.d0, i, j, k, l, m, n + !write(33, '(E15.7, 4X, 6(I4, 2X))') -integral/3.d0, l, m, n, i, j, k + endif enddo enddo enddo @@ -129,6 +179,9 @@ subroutine LMat_tilde_dump() enddo close(33) + call wall_time(t2) + print *, ' end after (min)', (t2-t1)/60.d0 + return end subroutine LMat_tilde_dump diff --git a/src/tc_bi_ortho/print_tc_energy.irp.f b/src/tc_bi_ortho/print_tc_energy.irp.f index b9f23a8a..2f667a48 100644 --- a/src/tc_bi_ortho/print_tc_energy.irp.f +++ b/src/tc_bi_ortho/print_tc_energy.irp.f @@ -8,11 +8,12 @@ program print_tc_energy print *, 'Hello world' my_grid_becke = .True. - !my_n_pt_r_grid = 30 - !my_n_pt_a_grid = 50 - my_n_pt_r_grid = 100 - my_n_pt_a_grid = 170 + my_n_pt_r_grid = 30 + my_n_pt_a_grid = 50 + + !my_n_pt_r_grid = 100 + !my_n_pt_a_grid = 170 !my_n_pt_r_grid = 100 !my_n_pt_a_grid = 266 diff --git a/src/tc_scf/print_tcscf_energy.irp.f b/src/tc_scf/print_tcscf_energy.irp.f new file mode 100644 index 00000000..96512cb7 --- /dev/null +++ b/src/tc_scf/print_tcscf_energy.irp.f @@ -0,0 +1,58 @@ +program print_tcscf_energy + + BEGIN_DOC + ! TODO : Put the documentation of the program here + END_DOC + + implicit none + + print *, 'Hello world' + my_grid_becke = .True. + + !my_n_pt_r_grid = 30 + !my_n_pt_a_grid = 50 + + my_n_pt_r_grid = 100 + my_n_pt_a_grid = 170 + + !my_n_pt_r_grid = 100 + !my_n_pt_a_grid = 266 + + touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid + + call main() + +end + +! --- + +subroutine main() + + implicit none + double precision :: etc_tot, etc_1e, etc_2e, etc_3e + + PROVIDE mu_erf + PROVIDE j1b_type + + print*, ' mu_erf = ', mu_erf + print*, ' j1b_type = ', j1b_type + + etc_tot = TC_HF_energy + etc_1e = TC_HF_one_e_energy + etc_2e = TC_HF_two_e_energy + etc_3e = 0.d0 + if(three_body_h_tc) then + !etc_3e = diag_three_elem_hf + etc_3e = tcscf_energy_3e_naive + endif + + print *, " E_TC = ", etc_tot + print *, " E_1e = ", etc_1e + print *, " E_2e = ", etc_2e + print *, " E_3e = ", etc_3e + + return +end subroutine main + +! --- + diff --git a/src/tc_scf/tc_scf.irp.f b/src/tc_scf/tc_scf.irp.f index 04c4f92d..2f2d803f 100644 --- a/src/tc_scf/tc_scf.irp.f +++ b/src/tc_scf/tc_scf.irp.f @@ -13,8 +13,13 @@ program tc_scf print *, ' starting ...' my_grid_becke = .True. - my_n_pt_r_grid = 30 - my_n_pt_a_grid = 50 + + !my_n_pt_r_grid = 30 + !my_n_pt_a_grid = 50 + + my_n_pt_r_grid = 100 + my_n_pt_a_grid = 170 + ! my_n_pt_r_grid = 10 ! small grid for quick debug ! my_n_pt_a_grid = 26 ! small grid for quick debug touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid diff --git a/src/tc_scf/tcscf_energy_naive.irp.f b/src/tc_scf/tcscf_energy_naive.irp.f new file mode 100644 index 00000000..82bb8799 --- /dev/null +++ b/src/tc_scf/tcscf_energy_naive.irp.f @@ -0,0 +1,80 @@ + +! --- + +BEGIN_PROVIDER [double precision, tcscf_energy_3e_naive] + + implicit none + integer :: i, j, k + integer :: neu, ned, D(elec_num) + integer :: ii, jj, kk + integer :: si, sj, sk + double precision :: I_ijk, I_jki, I_kij, I_jik, I_ikj, I_kji + double precision :: I_tot + + PROVIDE mo_l_coef mo_r_coef + + neu = elec_alpha_num + ned = elec_beta_num + if (neu > 0) D(1:neu) = [(2*i-1, i = 1, neu)] + if (ned > 0) D(neu+1:neu+ned) = [(2*i, i = 1, ned)] + + !print*, "D = " + !do i = 1, elec_num + ! ii = (D(i) - 1) / 2 + 1 + ! si = mod(D(i), 2) + ! print*, i, D(i), ii, si + !enddo + + tcscf_energy_3e_naive = 0.d0 + + do i = 1, elec_num - 2 + ii = (D(i) - 1) / 2 + 1 + si = mod(D(i), 2) + + do j = i + 1, elec_num - 1 + jj = (D(j) - 1) / 2 + 1 + sj = mod(D(j), 2) + + do k = j + 1, elec_num + kk = (D(k) - 1) / 2 + 1 + sk = mod(D(k), 2) + + call give_integrals_3_body_bi_ort(ii, jj, kk, ii, jj, kk, I_ijk) + I_tot = I_ijk + + if(sj==si .and. sk==sj) then + call give_integrals_3_body_bi_ort(ii, jj, kk, jj, kk, ii, I_jki) + I_tot += I_jki + endif + + if(sk==si .and. si==sj) then + call give_integrals_3_body_bi_ort(ii, jj, kk, kk, ii, jj, I_kij) + I_tot += I_kij + endif + + if(sj==si) then + call give_integrals_3_body_bi_ort(ii, jj, kk, jj, ii, kk, I_jik) + I_tot -= I_jik + endif + + if(sk==sj) then + call give_integrals_3_body_bi_ort(ii, jj, kk, ii, kk, jj, I_ikj) + I_tot -= I_ikj + endif + + if(sk==si) then + call give_integrals_3_body_bi_ort(ii, jj, kk, kk, jj, ii, I_kji) + I_tot -= I_kji + endif + + tcscf_energy_3e_naive += I_tot + enddo + enddo + enddo + + tcscf_energy_3e_naive = -tcscf_energy_3e_naive + +END_PROVIDER + +! --- + diff --git a/src/tc_scf/three_e_energy_bi_ortho.irp.f b/src/tc_scf/three_e_energy_bi_ortho.irp.f index 64212da8..0c9ebbd7 100644 --- a/src/tc_scf/three_e_energy_bi_ortho.irp.f +++ b/src/tc_scf/three_e_energy_bi_ortho.irp.f @@ -1,24 +1,32 @@ -subroutine contrib_3e_diag_sss(i,j,k,integral) - implicit none - integer, intent(in) :: i,j,k - BEGIN_DOC - ! returns the pure same spin contribution to diagonal matrix element of 3e term - END_DOC - double precision, intent(out) :: integral - double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int - call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j > - call give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int) ! < i k j | j i k > - call give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i > - integral = direct_int + c_3_int + c_minus_3_int - ! negative terms :: exchange contrib - call give_integrals_3_body_bi_ort(i, k, j, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13 - call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23 - call give_integrals_3_body_bi_ort(i, k, j, k, i, j, exch_12_int)!!! < i k j | k i j > : E_12 - integral += - exch_13_int - exch_23_int - exch_12_int - integral = -integral +subroutine contrib_3e_diag_sss(i, j, k, integral) + + BEGIN_DOC + ! returns the pure same spin contribution to diagonal matrix element of 3e term + END_DOC + + implicit none + integer, intent(in) :: i, j, k + double precision, intent(out) :: integral + double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int + + call give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j > + call give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int) ! < i k j | j i k > + call give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i > + integral = direct_int + c_3_int + c_minus_3_int + + ! negative terms :: exchange contrib + call give_integrals_3_body_bi_ort(i, k, j, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13 + call give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23 + call give_integrals_3_body_bi_ort(i, k, j, k, i, j, exch_12_int)!!! < i k j | k i j > : E_12 + + integral += - exch_13_int - exch_23_int - exch_12_int + integral = -integral + end +! --- + subroutine contrib_3e_diag_soo(i,j,k,integral) implicit none integer, intent(in) :: i,j,k @@ -51,23 +59,30 @@ subroutine give_aaa_contrib_bis(integral_aaa) end +! --- + subroutine give_aaa_contrib(integral_aaa) - implicit none - double precision, intent(out) :: integral_aaa - double precision :: integral - integer :: i,j,k - integral_aaa = 0.d0 - do i = 1, elec_alpha_num - do j = 1, elec_alpha_num - do k = 1, elec_alpha_num - call contrib_3e_diag_sss(i,j,k,integral) - integral_aaa += integral - enddo + + implicit none + integer :: i, j, k + double precision :: integral + double precision, intent(out) :: integral_aaa + + integral_aaa = 0.d0 + do i = 1, elec_alpha_num + do j = 1, elec_alpha_num + do k = 1, elec_alpha_num + call contrib_3e_diag_sss(i, j, k, integral) + integral_aaa += integral + enddo + enddo enddo - enddo - integral_aaa *= 1.d0/6.d0 + integral_aaa *= 1.d0/6.d0 + + return end +! --- subroutine give_aab_contrib(integral_aab) implicit none