DUMP for TCHint added

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
 

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(1) = (tmp0 * grad_u12(1) + tmp2 * grad_x_j1b_nucl_num(r1)) * tmp1
-  grad(2) = tmp3 * grad_u12(2) + tmp4 * grad_y_j1b_nucl_num(r1)
+  grad(2) = (tmp0 * grad_u12(2) + tmp2 * grad_y_j1b_nucl_num(r1)) * tmp1
-  grad(3) = tmp3 * grad_u12(3) + tmp4 * grad_z_j1b_nucl_num(r1)
+  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
+
+! ---
 
 

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) &

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
 

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_r_grid = 30
-  my_n_pt_a_grid = 170
+  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

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
+
+! ---
+

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

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
+
+! ---
+

src/tc_scf/three_e_energy_bi_ortho.irp.f

@@ -1,24 +1,32 @@
 
-subroutine contrib_3e_diag_sss(i,j,k,integral)
+subroutine contrib_3e_diag_sss(i, j, k, integral)
- implicit none
+
- integer, intent(in) :: i,j,k
+  BEGIN_DOC
- BEGIN_DOC
+  ! returns the pure same spin contribution to diagonal matrix element of 3e term
- ! returns the pure same spin contribution to diagonal matrix element of 3e term
+  END_DOC
- END_DOC
+
- double precision, intent(out) :: integral
+  implicit none
- double precision :: direct_int, exch_13_int, exch_23_int, exch_12_int, c_3_int, c_minus_3_int
+  integer,          intent(in)  :: i, j, k
- call  give_integrals_3_body_bi_ort(i, k, j, i, k, j, direct_int )!!! < i k j | i k j >
+  double precision, intent(out) :: integral
- call  give_integrals_3_body_bi_ort(i, k, j, j, i, k, c_3_int)      ! < i k j | j i k >
+  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, k, j, i, c_minus_3_int)! < i k j | k j i >
+
- integral = direct_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 >
- ! negative terms :: exchange contrib
+  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, j, k, i, exch_13_int)!!! < i k j | j k i > : E_13 
+  call give_integrals_3_body_bi_ort(i, k, j, k, j, i, c_minus_3_int)! < i k j | k j i >
- call  give_integrals_3_body_bi_ort(i, k, j, i, j, k, exch_23_int)!!! < i k j | i j k > : E_23
+  integral = direct_int + c_3_int + c_minus_3_int 
- 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 
+  ! negative terms :: exchange contrib
- integral = -integral
+  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
+  implicit none
- double precision :: integral
+  integer                       :: i, j, k
- integer :: i,j,k
+  double precision              :: integral
- integral_aaa = 0.d0
+  double precision, intent(out) :: integral_aaa
- do i = 1, elec_alpha_num
+
-  do j = 1, elec_alpha_num
+  integral_aaa = 0.d0
-   do k = 1, elec_alpha_num
+  do i = 1, elec_alpha_num
-    call contrib_3e_diag_sss(i,j,k,integral)
+    do j = 1, elec_alpha_num
-    integral_aaa += integral
+      do k = 1, elec_alpha_num
-   enddo
+        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