added keywords for r1 grid

src/ao_many_one_e_ints/grad_lapl_jmu_modif.irp.f

@@ -195,7 +195,7 @@ END_PROVIDER
 
 ! ---
 
-BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_final_grid)]
+BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_fit, (ao_num, ao_num, n_points_final_grid)]
 
   BEGIN_DOC
   !
@@ -212,14 +212,14 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
 
   double precision, external :: overlap_gauss_r12_ao_with1s
 
-  print*, ' providing v_ij_u_cst_mu_j1b ...'
+  print*, ' providing v_ij_u_cst_mu_j1b_fit ...'
   call wall_time(wall0)
 
   provide mu_erf final_grid_points j1b_pen
   PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x
   PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
 
-  v_ij_u_cst_mu_j1b = 0.d0
+  v_ij_u_cst_mu_j1b_fit = 0.d0
 
   !$OMP PARALLEL DEFAULT (NONE)                                      &
   !$OMP PRIVATE (ipoint, i, j, i_1s, i_fit, r, coef, beta, B_center, &
@@ -228,7 +228,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
   !$OMP          final_grid_points, ng_fit_jast,                     &
   !$OMP          expo_gauss_j_mu_x, coef_gauss_j_mu_x,               &
   !$OMP          List_all_comb_b2_coef, List_all_comb_b2_expo,       & 
-  !$OMP          List_all_comb_b2_cent, v_ij_u_cst_mu_j1b)
+  !$OMP          List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_fit)
   !$OMP DO
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -278,7 +278,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
 
         enddo
 
-        v_ij_u_cst_mu_j1b(j,i,ipoint) = tmp
+        v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = tmp
       enddo
     enddo
   enddo
@@ -288,13 +288,13 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b, (ao_num, ao_num, n_points_
   do ipoint = 1, n_points_final_grid
     do i = 2, ao_num
       do j = 1, i-1
-        v_ij_u_cst_mu_j1b(j,i,ipoint) = v_ij_u_cst_mu_j1b(i,j,ipoint)
+        v_ij_u_cst_mu_j1b_fit(j,i,ipoint) = v_ij_u_cst_mu_j1b_fit(i,j,ipoint)
       enddo
     enddo
   enddo
  
   call wall_time(wall1)
-  print*, ' wall time for v_ij_u_cst_mu_j1b', wall1 - wall0
+  print*, ' wall time for v_ij_u_cst_mu_j1b_fit', wall1 - wall0
 
 END_PROVIDER 
 
@@ -327,7 +327,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin
   call wall_time(wall0)
 
   provide mu_erf final_grid_points j1b_pen
-  PROVIDE ng_fit_jast expo_gauss_j_mu_x coef_gauss_j_mu_x
   PROVIDE List_all_comb_b2_size List_all_comb_b2_coef List_all_comb_b2_expo List_all_comb_b2_cent
 
   ct = inv_sq_pi_2 / mu_erf
@@ -340,7 +339,6 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_an, (ao_num, ao_num, n_poin
   !$OMP          int_e2, int_c3, int_e3)                             &
   !$OMP SHARED  (n_points_final_grid, ao_num, List_all_comb_b2_size, & 
   !$OMP          final_grid_points, mu_erf, ct,                      &
-  !$OMP          expo_gauss_j_mu_x, coef_gauss_j_mu_x,               &
   !$OMP          List_all_comb_b2_coef, List_all_comb_b2_expo,       & 
   !$OMP          List_all_comb_b2_cent, v_ij_u_cst_mu_j1b_an)
   !$OMP DO

src/bi_ort_ints/bi_ort_ints.irp.f

@@ -9,10 +9,9 @@ program bi_ort_ints
   implicit none
 
   my_grid_becke = .True.
-  !my_n_pt_r_grid = 10
+  PROVIDE tc_grid1_a tc_grid1_r
-  !my_n_pt_a_grid = 14
+  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_r_grid = 30
+  my_n_pt_a_grid = tc_grid1_a
-  my_n_pt_a_grid = 50
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
 ! call test_3e

src/bi_ort_ints/total_twoe_pot.irp.f

@@ -261,51 +261,55 @@ END_PROVIDER
 ! ---
 
 
- BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj, (mo_num,mo_num) ]
+ BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj,          (mo_num,mo_num)]
-&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num) ]
+&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_exchange, (mo_num,mo_num)]
-&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti, (mo_num,mo_num) ]
+&BEGIN_PROVIDER [ double precision, mo_bi_ortho_tc_two_e_jj_anti,     (mo_num,mo_num)]
-  implicit none
+
   BEGIN_DOC
   ! mo_bi_ortho_tc_two_e_jj(i,j) = J_ij = <ji|W-K|ji>
   ! mo_bi_ortho_tc_two_e_jj_exchange(i,j) = K_ij = <ij|W-K|ji>
   ! mo_bi_ortho_tc_two_e_jj_anti(i,j) = J_ij - K_ij
   END_DOC
 
-  integer                        :: i,j
+  implicit none
-  double precision               :: get_two_e_integral
+  integer :: i, j
 
-  mo_bi_ortho_tc_two_e_jj = 0.d0
+  mo_bi_ortho_tc_two_e_jj          = 0.d0
   mo_bi_ortho_tc_two_e_jj_exchange = 0.d0
 
-  do i=1,mo_num
+  do i = 1, mo_num
-    do j=1,mo_num
+    do j = 1, mo_num
-      mo_bi_ortho_tc_two_e_jj(i,j) = mo_bi_ortho_tc_two_e(j,i,j,i)
+      mo_bi_ortho_tc_two_e_jj(i,j)          = mo_bi_ortho_tc_two_e(j,i,j,i)
       mo_bi_ortho_tc_two_e_jj_exchange(i,j) = mo_bi_ortho_tc_two_e(i,j,j,i)
-      mo_bi_ortho_tc_two_e_jj_anti(i,j) = mo_bi_ortho_tc_two_e_jj(i,j) - mo_bi_ortho_tc_two_e_jj_exchange(i,j)
+      mo_bi_ortho_tc_two_e_jj_anti(i,j)     = mo_bi_ortho_tc_two_e_jj(i,j) - mo_bi_ortho_tc_two_e_jj_exchange(i,j)
     enddo
   enddo
 
 END_PROVIDER
 
- BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals, (mo_num,mo_num, mo_num)]
+! ---
-&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals,(mo_num,mo_num, mo_num)]
- implicit none
- BEGIN_DOC
- ! tc_2e_3idx_coulomb_integrals(j,k,i)  = <jk|ji> 
- !
- ! tc_2e_3idx_exchange_integrals(j,k,i) = <kj|ji> 
- END_DOC
- integer :: i,j,k,l
- double precision :: get_two_e_integral
- double precision :: integral
 
- do i = 1, mo_num
+ BEGIN_PROVIDER [double precision, tc_2e_3idx_coulomb_integrals , (mo_num,mo_num,mo_num)]
-  do k = 1, mo_num
+&BEGIN_PROVIDER [double precision, tc_2e_3idx_exchange_integrals, (mo_num,mo_num,mo_num)]
-   do j = 1, mo_num
+
-     tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i ) 
+  BEGIN_DOC
-     tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i ) 
+  ! tc_2e_3idx_coulomb_integrals (j,k,i) = <jk|ji> 
-   enddo
+  ! tc_2e_3idx_exchange_integrals(j,k,i) = <kj|ji> 
+  END_DOC
+
+  implicit none
+  integer :: i, j, k
+
+  do i = 1, mo_num
+    do k = 1, mo_num
+      do j = 1, mo_num
+        tc_2e_3idx_coulomb_integrals(j, k,i) = mo_bi_ortho_tc_two_e(j ,k ,j ,i ) 
+        tc_2e_3idx_exchange_integrals(j,k,i) = mo_bi_ortho_tc_two_e(k ,j ,j ,i ) 
+      enddo
+    enddo
   enddo
- enddo
 
 END_PROVIDER
+
+! ---
+

src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f

@@ -1,47 +1,54 @@
+
+! ---
+
 subroutine run_stochastic_cipsi
+
+  BEGIN_DOC
+  ! Selected Full Configuration Interaction with Stochastic selection and PT2.
+  END_DOC
+
   use selection_types
   implicit none
-  BEGIN_DOC
+  integer                       :: i, j, k, ndet
-! Selected Full Configuration Interaction with Stochastic selection and PT2.
+  integer                       :: to_select
-  END_DOC
+  logical                       :: print_pt2
-  integer                        :: i,j,k,ndet
+  logical                       :: has
-  double precision, allocatable  :: zeros(:)
+  type(pt2_type)                :: pt2_data, pt2_data_err
-  integer                        :: to_select
+  double precision              :: rss
-  type(pt2_type)                 :: pt2_data, pt2_data_err
+  double precision              :: correlation_energy_ratio, E_denom, E_tc, norm
-  logical, external              :: qp_stop
+  double precision              :: hf_energy_ref
-  logical                        :: print_pt2
+  double precision              :: relative_error
+  double precision, allocatable :: ept2(:), pt1(:), extrap_energy(:)
+  double precision, allocatable :: zeros(:)
 
-  double precision :: rss
+  logical,          external    :: qp_stop
-  double precision, external :: memory_of_double
+  double precision, external    :: memory_of_double
-  double precision :: correlation_energy_ratio,E_denom,E_tc,norm
+
-  double precision, allocatable :: ept2(:), pt1(:),extrap_energy(:)
+  PROVIDE mo_l_coef mo_r_coef
   PROVIDE H_apply_buffer_allocated distributed_davidson 
 
-  print*,'Diagonal elements of the Fock matrix '
+  print*, ' Diagonal elements of the Fock matrix '
   do i = 1, mo_num
-   write(*,*)i,Fock_matrix_tc_mo_tot(i,i)
+    write(*,*) i, Fock_matrix_tc_mo_tot(i,i)
   enddo
+
   N_iter = 1
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators
 
   rss = memory_of_double(N_states)*4.d0
-  call check_mem(rss,irp_here)
+  call check_mem(rss, irp_here)
 
-  allocate (zeros(N_states))
+  allocate(zeros(N_states))
   call pt2_alloc(pt2_data, N_states)
   call pt2_alloc(pt2_data_err, N_states)
 
-  double precision               :: hf_energy_ref
+  relative_error = PT2_relative_error
-  logical                        :: has
-  double precision               :: relative_error
 
-  relative_error=PT2_relative_error
+  zeros               = 0.d0
-
+  pt2_data % pt2      = -huge(1.e0)
-  zeros = 0.d0
+  pt2_data % rpt2     = -huge(1.e0)
-  pt2_data % pt2   = -huge(1.e0)
+  pt2_data % overlap  = 0.d0
-  pt2_data % rpt2  = -huge(1.e0)
-  pt2_data % overlap= 0.d0
   pt2_data % variance = huge(1.e0)
 
   !!!! WARNING  !!!! SEEMS TO BE PROBLEM WTH make_s2_eigenfunction !!!! THE DETERMINANTS CAN APPEAR TWICE IN THE WFT DURING SELECTION
@@ -49,7 +56,7 @@ subroutine run_stochastic_cipsi
 !    call make_s2_eigenfunction
 !  endif
   print_pt2 = .False.
-  call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
+  call diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2)
 !  call routine_save_right
 
 
@@ -74,14 +81,12 @@ subroutine run_stochastic_cipsi
 ! soft_touch thresh_it_dav
 
   print_pt2 = .True.
-  do while (                                                         &
+  do while( (N_det < N_det_max) .and. &
-        (N_det < N_det_max) .and.                                    &
+            (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max))
-        (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max)          &
-        )
-      print*,'maxval(abs(pt2_data % pt2(1:N_states)))',maxval(abs(pt2_data % pt2(1:N_states)))
-      print*,pt2_max
-      write(*,'(A)')  '--------------------------------------------------------------------------------'
 
+    print*,'maxval(abs(pt2_data % pt2(1:N_states)))',maxval(abs(pt2_data % pt2(1:N_states)))
+    print*,pt2_max
+    write(*,'(A)')  '--------------------------------------------------------------------------------'
 
     to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
     to_select = max(N_states_diag, to_select)
@@ -94,8 +99,7 @@ subroutine run_stochastic_cipsi
     call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
 !    stop
 
-    call print_summary(psi_energy_with_nucl_rep, &
+    call print_summary(psi_energy_with_nucl_rep, pt2_data, pt2_data_err, N_det, N_configuration, N_states, psi_s2)
-       pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
 
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
 
@@ -109,13 +113,13 @@ subroutine run_stochastic_cipsi
     ! Add selected determinants
     call copy_H_apply_buffer_to_wf_tc()
 
-    PROVIDE  psi_l_coef_bi_ortho psi_r_coef_bi_ortho
+    PROVIDE psi_l_coef_bi_ortho psi_r_coef_bi_ortho
-    PROVIDE  psi_det
+    PROVIDE psi_det
-    PROVIDE  psi_det_sorted_tc
+    PROVIDE psi_det_sorted_tc
 
     ept2(N_iter-1) = E_tc + nuclear_repulsion + (pt2_data % pt2(1))/norm
-    pt1(N_iter-1) = dsqrt(pt2_data % overlap(1,1))
+    pt1(N_iter-1)  = dsqrt(pt2_data % overlap(1,1))
-    call diagonalize_CI_tc_bi_ortho(ndet, E_tc,norm,pt2_data,print_pt2)
+    call diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2)
 !    stop
     if (qp_stop()) exit
   enddo

src/fci_tc_bi/diagonalize_ci.irp.f

@@ -17,6 +17,8 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2)
   integer                         :: i, j
   double precision                :: pt2_tmp, pt1_norm, rpt2_tmp, abs_pt2
 
+  PROVIDE mo_l_coef mo_r_coef
+
   pt2_tmp  = pt2_data % pt2(1)
   abs_pt2  = pt2_data % variance(1)
   pt1_norm = pt2_data % overlap(1,1)
@@ -60,7 +62,7 @@ subroutine diagonalize_CI_tc_bi_ortho(ndet, E_tc, norm, pt2_data, print_pt2)
       psi_coef(i,j)            = dabs(psi_l_coef_bi_ortho(i,j) * psi_r_coef_bi_ortho(i,j))   
     enddo
   enddo
-  SOFT_TOUCH eigval_left_tc_bi_orth  eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth 
+  SOFT_TOUCH eigval_left_tc_bi_orth eigval_right_tc_bi_orth leigvec_tc_bi_orth reigvec_tc_bi_orth norm_ground_left_right_bi_orth 
   SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho psi_coef psi_energy psi_s2
 
   call save_tc_bi_ortho_wavefunction()

src/fci_tc_bi/fci_tc_bi_ortho.irp.f

@@ -1,5 +1,8 @@
-program fci
+
-  implicit none
+! ---
+
+program fci_tc_bi
+
   BEGIN_DOC
   ! Selected Full Configuration Interaction with stochastic selection
   ! and PT2.
@@ -36,10 +39,12 @@ program fci
   !
   END_DOC
 
+  implicit none
 
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid 
 
   pruning = -1.d0
@@ -62,18 +67,18 @@ subroutine run_cipsi_tc()
 
     PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e
 
-    if(elec_alpha_num+elec_beta_num .ge. 3) then
+    if((elec_alpha_num+elec_beta_num) .ge. 3) then
       if(three_body_h_tc) then
         call provide_all_three_ints_bi_ortho()
       endif
     endif
 
-    FREE int2_grad1_u12_ao
+    FREE int2_grad1_u12_ao int2_grad1_u12_ao_t int2_grad1_u12_ao_transp
-    FREE int2_grad1_u12_bimo_transp int2_grad1_u12_ao_transp
+    FREE int2_grad1_u12_bimo_transp 
 
     write(json_unit,json_array_open_fmt) 'fci_tc'
 
-    if (do_pt2) then
+    if(do_pt2) then
       call run_stochastic_cipsi()
     else
       call run_cipsi()
@@ -88,14 +93,14 @@ subroutine run_cipsi_tc()
 
     PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks
 
-    if(elec_alpha_num+elec_beta_num .ge. 3) then
+    if((elec_alpha_num+elec_beta_num) .ge. 3) then
       if(three_body_h_tc) then
         call provide_all_three_ints_bi_ortho()
       endif
     endif
 
-    FREE int2_grad1_u12_ao
+    FREE int2_grad1_u12_ao int2_grad1_u12_ao_t int2_grad1_u12_ao_transp
-    FREE int2_grad1_u12_bimo_transp int2_grad1_u12_ao_transp
+    FREE int2_grad1_u12_bimo_transp 
 
     call run_slave_cipsi
 

src/fci_tc_bi/pt2_tc.irp.f

@@ -1,31 +1,42 @@
+
+! ---
+
 program tc_pt2_prog
+
   implicit none
+
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid 
+
   pruning = -1.d0
   touch pruning
+
 !  pt2_relative_error = 0.01d0
 !  touch pt2_relative_error
-  call run_pt2_tc
+  call run_pt2_tc()
 
 end
 
+! ---
 
-subroutine run_pt2_tc
+subroutine run_pt2_tc()
 
- implicit none
+  implicit none
 
   PROVIDE psi_det psi_coef mo_bi_ortho_tc_two_e mo_bi_ortho_tc_one_e
-  if(elec_alpha_num+elec_beta_num.ge.3)then
+
+  if(elec_alpha_num+elec_beta_num.ge.3) then
     if(three_body_h_tc)then
-      call provide_all_three_ints_bi_ortho
+      call provide_all_three_ints_bi_ortho()
     endif
   endif
-  ! ---
-
- call tc_pt2
 
+  call tc_pt2()
 
 end
+
+! ---
+

src/non_h_ints_mu/debug_fit.irp.f

@@ -6,13 +6,9 @@ program debug_fit
   implicit none
 
   my_grid_becke  = .True.
-
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_r_grid = 30
+  my_n_pt_r_grid = tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_a_grid = tc_grid1_a
-  !my_n_pt_r_grid = 100
-  !my_n_pt_a_grid = 170
-  !my_n_pt_r_grid = 150
-  !my_n_pt_a_grid = 194
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   PROVIDE mu_erf j1b_pen

src/non_h_ints_mu/debug_integ_jmu_modif.irp.f

@@ -6,13 +6,9 @@ program debug_integ_jmu_modif
   implicit none
 
   my_grid_becke  = .True.
-
+  PROVIDE tc_grid1_a tc_grid1_r
-  !my_n_pt_r_grid = 30
+  my_n_pt_r_grid = tc_grid1_r
-  !my_n_pt_a_grid = 50
+  my_n_pt_a_grid = tc_grid1_a
-  !my_n_pt_r_grid = 100
-  !my_n_pt_a_grid = 170
-  my_n_pt_r_grid = 150
-  my_n_pt_a_grid = 194
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   PROVIDE mu_erf j1b_pen
@@ -48,22 +44,21 @@ subroutine test_v_ij_u_cst_mu_j1b()
 
   print*, ' test_v_ij_u_cst_mu_j1b ...'
 
-  PROVIDE v_ij_u_cst_mu_j1b
+  PROVIDE v_ij_u_cst_mu_j1b_fit
 
   eps_ij  = 1d-3
   acc_tot = 0.d0
   normalz = 0.d0
 
-  !do ipoint = 1, 10
   do ipoint = 1, n_points_final_grid
     do j = 1, ao_num
       do i = 1, ao_num
 
-        i_exc  =     v_ij_u_cst_mu_j1b(i,j,ipoint) 
+        i_exc  =     v_ij_u_cst_mu_j1b_fit(i,j,ipoint) 
-        i_num  = num_v_ij_u_cst_mu_j1b(i,j,ipoint)
+        i_num  = num_v_ij_u_cst_mu_j1b    (i,j,ipoint)
         acc_ij = dabs(i_exc - i_num)
         if(acc_ij .gt. eps_ij) then
-          print *, ' problem in v_ij_u_cst_mu_j1b on', i, j, ipoint
+          print *, ' problem in v_ij_u_cst_mu_j1b_fit on', i, j, ipoint
           print *, ' analyt integ = ', i_exc
           print *, ' numeri integ = ', i_num
           print *, ' diff         = ', acc_ij

src/non_h_ints_mu/new_grad_tc.irp.f

@@ -1,68 +1,3 @@
-! ---
-
-!BEGIN_PROVIDER [ double precision, int1_grad2_u12_ao, (3, ao_num, ao_num, n_points_final_grid)]
-!
-!  BEGIN_DOC
-!  !
-!  ! int1_grad2_u12_ao(:,i,j,ipoint) = \int dr1 [-1 * \grad_r2 J(r1,r2)] \phi_i(r1) \phi_j(r1) 
-!  !
-!  ! where r1 = r(ipoint)
-!  !
-!  ! if J(r1,r2) = u12:
-!  !
-!  ! int1_grad2_u12_ao(:,i,j,ipoint) = +0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] \phi_i(r1) \phi_j(r1)
-!  !                                 = -0.5 * [ v_ij_erf_rk_cst_mu(i,j,ipoint) * r(:) - x_v_ij_erf_rk_cst_mu(i,j,ipoint,:) ]
-!  !                                 = -int2_grad1_u12_ao(i,j,ipoint,:)
-!  !
-!  ! if J(r1,r2) = u12 x v1 x v2
-!  !
-!  ! int1_grad2_u12_ao(:,i,j,ipoint) =      v2    x [ 0.5 x \int dr1 [-(r1 - r2) (erf(mu * r12)-1)r_12] v1 \phi_i(r1) \phi_j(r1) ]
-!  !                                 - \grad_2 v2 x [       \int dr1                   u12              v1 \phi_i(r1) \phi_j(r1) ] 
-!  !                                 =   -0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) 
-!  !                                 +    0.5 v_1b(ipoint) * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,:) 
-!  !                                 - v_1b_grad[:,ipoint] * v_ij_u_cst_mu_j1b(i,j,ipoint)
-!  !
-!  !
-!  END_DOC
-!
-!  implicit none
-!  integer          :: ipoint, i, j
-!  double precision :: x, y, z, tmp_x, tmp_y, tmp_z, tmp0, tmp1, tmp2
-!
-!  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)
-!
-!      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)
-!
-!          int1_grad2_u12_ao(1,i,j,ipoint) = -tmp1 * x + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,1) - tmp2 * tmp_x
-!          int1_grad2_u12_ao(2,i,j,ipoint) = -tmp1 * y + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,2) - tmp2 * tmp_y
-!          int1_grad2_u12_ao(3,i,j,ipoint) = -tmp1 * z + tmp0 * x_v_ij_erf_rk_cst_mu_j1b(i,j,ipoint,3) - tmp2 * tmp_z
-!        enddo
-!      enddo
-!    enddo
-!
-!  else
-!
-!    int1_grad2_u12_ao = -1.d0 * int2_grad1_u12_ao
-!
-!  endif
-!
-!END_PROVIDER 
 
 ! ---
 
@@ -98,22 +33,14 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_loop, (ao_num, ao_num, ao_
     weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
 
     do i = 1, ao_num
-      !ao_i_r  = weight1 * aos_in_r_array_transp         (ipoint,i)
-      !ao_i_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,1)
-      !ao_i_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,2)
-      !ao_i_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,i,3)
       ao_i_r  = weight1 * aos_in_r_array     (i,ipoint)
       ao_i_dx = weight1 * aos_grad_in_r_array(i,ipoint,1)
       ao_i_dy = weight1 * aos_grad_in_r_array(i,ipoint,2)
       ao_i_dz = weight1 * aos_grad_in_r_array(i,ipoint,3)
 
       do k = 1, ao_num
-        !ao_k_r = aos_in_r_array_transp(ipoint,k)
         ao_k_r = aos_in_r_array(k,ipoint)
 
-        !tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1) 
-        !tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2) 
-        !tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3) 
         tmp_x = ao_k_r * ao_i_dx - ao_i_r * aos_grad_in_r_array(k,ipoint,1) 
         tmp_y = ao_k_r * ao_i_dy - ao_i_r * aos_grad_in_r_array(k,ipoint,2) 
         tmp_z = ao_k_r * ao_i_dz - ao_i_r * aos_grad_in_r_array(k,ipoint,3) 
@@ -134,44 +61,11 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_loop, (ao_num, ao_num, ao_
 
   ! ---
 
-  !do ipoint = 1, n_points_final_grid
-  !  weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
-
-  !  do l = 1, ao_num
-  !    ao_l_r  = weight1 * aos_in_r_array_transp         (ipoint,l)
-  !    ao_l_dx = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,1)
-  !    ao_l_dy = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,2)
-  !    ao_l_dz = weight1 * aos_grad_in_r_array_transp_bis(ipoint,l,3)
-
-  !    do j = 1, ao_num
-  !      ao_j_r = aos_in_r_array_transp(ipoint,j)
-
-  !      tmp_x = ao_j_r * ao_l_dx - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,1) 
-  !      tmp_y = ao_j_r * ao_l_dy - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,2) 
-  !      tmp_z = ao_j_r * ao_l_dz - ao_l_r * aos_grad_in_r_array_transp_bis(ipoint,j,3) 
-
-  !      do i = 1, ao_num
-  !        do k = 1, ao_num
-
-  !          contrib_x = int2_grad1_u12_ao(k,i,ipoint,1) * tmp_x 
-  !          contrib_y = int2_grad1_u12_ao(k,i,ipoint,2) * tmp_y 
-  !          contrib_z = int2_grad1_u12_ao(k,i,ipoint,3) * tmp_z 
-
-  !          ac_mat(k,i,l,j) += contrib_x + contrib_y + contrib_z
-  !        enddo
-  !      enddo
-  !    enddo
-  !  enddo
-  !enddo
-
-  ! ---
- 
   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_loop(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
-          !tc_grad_and_lapl_ao_loop(k,i,l,j) = ac_mat(k,i,l,j)
         enddo
       enddo
     enddo

src/non_h_ints_mu/test_non_h_ints.irp.f

@@ -1,19 +1,18 @@
+
+! ---
+
 program test_non_h
 
   implicit none
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 50
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 74
+  my_n_pt_r_grid = tc_grid1_r
-  !my_n_pt_r_grid = 400
+  my_n_pt_a_grid = tc_grid1_a
-  !my_n_pt_a_grid = 974
-
-!  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
 
-  !call routine_grad_squared
+  !call routine_grad_squared()
-  !call routine_fit
+  !call routine_fit()
   
   call test_ipp()
 end

src/tc_bi_ortho/compute_deltamu_right.irp.f

@@ -3,8 +3,9 @@ program compute_deltamu_right
   implicit none
 
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   read_wf = .True.

src/tc_bi_ortho/print_tc_dump.irp.f

@@ -6,10 +6,9 @@ program tc_bi_ortho
   implicit none
 
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
-  !my_n_pt_r_grid = 100
+  my_n_pt_a_grid = tc_grid1_a
-  !my_n_pt_a_grid = 170
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   call ERI_dump()

src/tc_bi_ortho/print_tc_energy.irp.f

@@ -7,16 +7,12 @@ program print_tc_energy
   implicit none
 
   print *, 'Hello world'
+
   my_grid_becke = .True.
-
+  PROVIDE tc_grid1_a tc_grid1_r
-  !my_n_pt_r_grid = 30
+  my_n_pt_r_grid = tc_grid1_r
-  !my_n_pt_a_grid = 50
+  my_n_pt_a_grid = tc_grid1_a
-
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-  my_n_pt_r_grid = 100
-  my_n_pt_a_grid = 170
-
-  !my_n_pt_r_grid = 100
-  !my_n_pt_a_grid = 266
 
   read_wf = .True.
   touch read_wf
@@ -24,8 +20,6 @@ program print_tc_energy
   PROVIDE j1b_type
   print*, 'j1b_type = ', j1b_type
 
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-
   call write_tc_energy()
 
 end

src/tc_bi_ortho/print_tc_spin_dens.irp.f

@@ -1,16 +1,26 @@
+
+! ---
+
 program test_spin_dens
-  implicit none
+
   BEGIN_DOC
-! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
+  ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
   END_DOC
+
+  implicit none
+
   print *, 'Hello world'
+
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
- call tc_print_mulliken_sd
+  call tc_print_mulliken_sd()
-! call test
+  !call test
 
 end

src/tc_bi_ortho/print_tc_var.irp.f

@@ -7,12 +7,15 @@ program print_tc_var
   implicit none
 
   print *, 'Hello world'
+
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   call write_tc_var()
 

src/tc_bi_ortho/print_tc_wf.irp.f

@@ -1,20 +1,31 @@
+
+! ---
+
 program print_tc_bi_ortho
-  implicit none
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! 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
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
 !  if(three_body_h_tc)then
 !   call provide_all_three_ints_bi_ortho
 !  endif
 !  call routine
- call write_l_r_wf
+  call write_l_r_wf
+
 end
 
 subroutine write_l_r_wf

src/tc_bi_ortho/pt2_tc_cisd.irp.f

@@ -7,12 +7,16 @@ program pt2_tc_cisd
   !
   END_DOC
 
+  implicit none
+
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   print*, ' nb of states = ', N_states
   print*, ' nb of det    = ', N_det

src/tc_bi_ortho/save_tc_bi_ortho_nat.irp.f

@@ -1,35 +1,59 @@
- program tc_natorb_bi_ortho
-   implicit none
-   BEGIN_DOC
- ! TODO : Put the documentation of the program here
-   END_DOC
-   print *, 'Hello world'
-   my_grid_becke = .True.
-   my_n_pt_r_grid = 30
-   my_n_pt_a_grid = 50
-   read_wf = .True.
-   touch read_wf
-   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-   call print_energy_and_mos
-   call save_tc_natorb
-!   call minimize_tc_orb_angles
- end
 
- subroutine save_tc_natorb 
+! ---
+
+program tc_natorb_bi_ortho
+
+  BEGIN_DOC
+  ! TODO : Put the documentation of the program here
+  END_DOC
+
   implicit none
+
+  print *, 'Hello world'
+
+  my_grid_becke = .True.
+  PROVIDE tc_grid1_a tc_grid1_r
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
+  read_wf = .True.
+  touch read_wf
+
+  call print_energy_and_mos()
+  call save_tc_natorb()
+  !call minimize_tc_orb_angles()
+
+end
+
+! ---
+ 
+subroutine save_tc_natorb()
+
+  implicit none
+
   print*,'Saving the natorbs '
+
   provide natorb_tc_leigvec_ao natorb_tc_reigvec_ao
+
   call ezfio_set_bi_ortho_mos_mo_l_coef(natorb_tc_leigvec_ao)
   call ezfio_set_bi_ortho_mos_mo_r_coef(natorb_tc_reigvec_ao)
-  call save_ref_determinant_nstates_1
+  call save_ref_determinant_nstates_1()
   call ezfio_set_determinants_read_wf(.False.)
- end
 
- subroutine save_ref_determinant_nstates_1
+end
-   implicit none
+
+! ---
+ 
+subroutine save_ref_determinant_nstates_1()
+
    use bitmasks
-   double precision               :: buffer(1,N_states)
+   implicit none
+   double precision :: buffer(1,N_states)
+
    buffer = 0.d0
    buffer(1,1) = 1.d0
-   call save_wavefunction_general(1,1,ref_bitmask,1,buffer)
+   call save_wavefunction_general(1, 1, ref_bitmask, 1, buffer)
- end
+
+end
+

src/tc_bi_ortho/select_dets_bi_ortho.irp.f

@@ -1,15 +1,24 @@
-program tc_bi_ortho
+
-  implicit none
+! ---
+
+program select_dets_bi_ortho()
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! 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
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   !!!!!!!!!!!!!!! WARNING NO 3-BODY
   !!!!!!!!!!!!!!! WARNING NO 3-BODY
@@ -22,6 +31,8 @@ program tc_bi_ortho
 ! call test
 end
 
+! ---
+
 subroutine routine_test
  implicit none
  use bitmasks ! you need to include the bitmasks_module.f90 features
@@ -57,5 +68,7 @@ subroutine routine_test
  enddo
  call save_wavefunction_general(n_good,n_states,dets,n_good,coef_new)
 
-
 end
+
+! ---
+

src/tc_bi_ortho/slater_tc_3e_slow.irp.f

@@ -1,4 +1,6 @@
 
+! ---
+
 subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
 
   BEGIN_DOC
@@ -22,12 +24,12 @@ subroutine diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
 
   if(core_tc_op) then
     do i = 1, Nint
-      key_i_core(i,1) = xor(key_i(i,1),core_bitmask(i,1))
+      key_i_core(i,1) = xor(key_i(i,1), core_bitmask(i,1))
-      key_i_core(i,2) = xor(key_i(i,2),core_bitmask(i,2))
+      key_i_core(i,2) = xor(key_i(i,2), core_bitmask(i,2))
     enddo
-    call bitstring_to_list_ab(key_i_core,occ,Ne,Nint)
+    call bitstring_to_list_ab(key_i_core, occ, Ne, Nint)
   else
-    call bitstring_to_list_ab(key_i,occ,Ne,Nint)
+    call bitstring_to_list_ab(key_i, occ, Ne, Nint)
   endif
 
   hthree = 0.d0

src/tc_bi_ortho/slater_tc_opt.irp.f

@@ -47,13 +47,19 @@ subroutine htilde_mu_mat_opt_bi_ortho_tot(key_j, key_i, Nint, htot)
   END_DOC
 
   use bitmasks
- integer,           intent(in) :: Nint
+  integer,           intent(in) :: Nint
- integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
+  integer(bit_kind), intent(in) :: key_i(Nint,2), key_j(Nint,2)
- double precision, intent(out) :: htot
+  double precision, intent(out) :: htot
- double precision :: hmono, htwoe, hthree
+  double precision :: hmono, htwoe, hthree
- call htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
+
+  call htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
+
 end
+
+! ---
+
 subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
+
   BEGIN_DOC
   !
   ! <key_j | H_tilde | key_i> where |key_j> is developed on the LEFT basis and |key_i> is developed on the RIGHT basis
@@ -81,11 +87,11 @@ subroutine htilde_mu_mat_opt_bi_ortho(key_j, key_i, Nint, hmono, htwoe, hthree,
   call get_excitation_degree(key_i, key_j, degree, Nint)
   if(degree.gt.2) return
 
-  if(degree == 0)then
+  if(degree == 0) then
     call diag_htilde_mu_mat_fock_bi_ortho  (Nint, key_i, hmono, htwoe, hthree, htot)
-  else if (degree == 1)then
+  else if (degree == 1) then
-    call single_htilde_mu_mat_fock_bi_ortho(Nint,key_j, key_i , hmono, htwoe, hthree, htot)
+    call single_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i , hmono, htwoe, hthree, htot)
-  else if(degree == 2)then
+  else if(degree == 2) then
     call double_htilde_mu_mat_fock_bi_ortho(Nint, key_j, key_i, hmono, htwoe, hthree, htot)
   endif
 

src/tc_bi_ortho/slater_tc_opt_diag.irp.f

@@ -111,7 +111,7 @@ end
 
 ! ---
 
-subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
+subroutine ac_tc_operator(iorb, ispin, key, hmono, htwoe, hthree, Nint, na, nb)
 
   BEGIN_DOC
   ! Routine that computes one- and two-body energy corresponding 
@@ -127,16 +127,16 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
 
   use bitmasks
   implicit none
-  integer, intent(in)            :: iorb, ispin, Nint
+  integer,           intent(in)    :: iorb, ispin, Nint
-  integer, intent(inout)         :: na, nb
+  integer,           intent(inout) :: na, nb
   integer(bit_kind), intent(inout) :: key(Nint,2)
-  double precision, intent(inout) :: hmono,htwoe,hthree
+  double precision,  intent(inout) :: hmono, htwoe, hthree
-
-  integer                        :: occ(Nint*bit_kind_size,2)
-  integer                        :: other_spin
-  integer                        :: k,l,i,jj,mm,j,m
-  double precision ::  direct_int, exchange_int
 
+  integer                          :: occ(Nint*bit_kind_size,2)
+  integer                          :: other_spin
+  integer                          :: k, l, i, jj, mm, j, m
+  integer                          :: tmp(2)
+  double precision                 :: direct_int, exchange_int
   
   if (iorb < 1) then
     print *,  irp_here, ': iorb < 1'
@@ -153,7 +153,6 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
   ASSERT (ispin < 3)
   ASSERT (Nint > 0)
 
-  integer                        :: tmp(2)
   !DIR$ FORCEINLINE
   call bitstring_to_list_ab(key, occ, tmp, Nint)
   ASSERT (tmp(1) == elec_alpha_num)
@@ -169,50 +168,54 @@ subroutine ac_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
   hmono = hmono + mo_bi_ortho_tc_one_e(iorb,iorb)
 
   ! Same spin
-  do i=1,na
+  do i = 1, na
     htwoe = htwoe + mo_bi_ortho_tc_two_e_jj_anti(occ(i,ispin),iorb)
   enddo
 
   ! Opposite spin
-  do i=1,nb
+  do i = 1, nb
     htwoe = htwoe + mo_bi_ortho_tc_two_e_jj(occ(i,other_spin),iorb)
   enddo
 
-  if(three_body_h_tc.and.elec_num.gt.2.and.three_e_3_idx_term)then
+  if(three_body_h_tc .and. (elec_num.gt.2) .and. three_e_3_idx_term) then
-   !!!!! 3-e part 
+
-   !! same-spin/same-spin
+    !!!!! 3-e part 
-   do j = 1, na
+    !! same-spin/same-spin
-    jj = occ(j,ispin)
+    do j = 1, na
-    do m = j+1, na
+      jj = occ(j,ispin)
-     mm = occ(m,ispin)
+      do m = j+1, na
-     hthree += three_e_diag_parrallel_spin_prov(mm,jj,iorb)
+        mm = occ(m,ispin)
+        hthree += three_e_diag_parrallel_spin_prov(mm,jj,iorb)
+      enddo
     enddo
-   enddo
+    !! same-spin/oposite-spin
-   !! same-spin/oposite-spin
+    do j = 1, na
-   do j = 1, na
+      jj = occ(j,ispin)
-    jj = occ(j,ispin)
+      do m = 1, nb
-    do m = 1, nb
+        mm = occ(m,other_spin)
-     mm = occ(m,other_spin)
+        direct_int   = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
-     direct_int   = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+        exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
-     exchange_int = three_e_3_idx_exch12_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+        hthree += direct_int - exchange_int
-     hthree += direct_int - exchange_int
+      enddo
     enddo
-   enddo
+    !! oposite-spin/opposite-spin
-   !! oposite-spin/opposite-spin
     do j = 1, nb
-     jj = occ(j,other_spin) 
+      jj = occ(j,other_spin) 
-     do m = j+1, nb 
+      do m = j+1, nb 
-      mm = occ(m,other_spin) 
+        mm = occ(m,other_spin) 
-      direct_int = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+        direct_int   = three_e_3_idx_direct_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
-      exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
+        exchange_int = three_e_3_idx_exch23_bi_ort(mm,jj,iorb) ! USES 3-IDX TENSOR 
-      hthree += direct_int - exchange_int
+        hthree += direct_int - exchange_int
-     enddo
+      enddo
     enddo
   endif
 
   na = na+1
+
 end
 
+! ---
+
 subroutine a_tc_operator(iorb,ispin,key,hmono,htwoe,hthree,Nint,na,nb)
   use bitmasks
   implicit none

src/tc_bi_ortho/slater_tc_slow.irp.f

@@ -21,7 +21,7 @@ subroutine htilde_mu_mat_bi_ortho_tot_slow(key_j, key_i, Nint, htot)
   integer :: degree
 
   call get_excitation_degree(key_j, key_i, degree, Nint)
-  if(degree.gt.2)then
+  if(degree.gt.2) then
     htot = 0.d0
   else
     call htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree, htot)
@@ -60,22 +60,22 @@ subroutine htilde_mu_mat_bi_ortho_slow(key_j, key_i, Nint, hmono, htwoe, hthree,
   call get_excitation_degree(key_i, key_j, degree, Nint)
   if(degree.gt.2) return
 
-  if(degree == 0)then
+  if(degree == 0) then
     call diag_htilde_mu_mat_bi_ortho_slow(Nint, key_i, hmono, htwoe, htot)
-  else if (degree == 1)then
+  else if (degree == 1) then
     call single_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
-  else if(degree == 2)then
+  else if(degree == 2) then
     call double_htilde_mu_mat_bi_ortho_slow(Nint, key_j, key_i, hmono, htwoe, htot)
   endif
 
   if(three_body_h_tc) then
     if(degree == 2) then
-      if(.not.double_normal_ord.and.elec_num.gt.2.and.three_e_5_idx_term) then
+      if((.not.double_normal_ord) .and. (elec_num .gt. 2) .and. three_e_5_idx_term) then
         call double_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
       endif
-    else if(degree == 1.and.elec_num.gt.2.and.three_e_4_idx_term) then
+    else if((degree == 1) .and. (elec_num .gt. 2) .and. three_e_4_idx_term) then
       call single_htilde_three_body_ints_bi_ort_slow(Nint, key_j, key_i, hthree)
-    else if(degree == 0.and.elec_num.gt.2.and.three_e_3_idx_term) then
+    else if((degree == 0) .and. (elec_num .gt. 2) .and. three_e_3_idx_term) then
       call diag_htilde_three_body_ints_bi_ort_slow(Nint, key_i, hthree)
     endif
   endif

src/tc_bi_ortho/tc_bi_ortho.irp.f

@@ -8,11 +8,13 @@ program tc_bi_ortho
   END_DOC
 
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   print*, ' nb of states = ', N_states
   print*, ' nb of det    = ', N_det
@@ -20,22 +22,29 @@ program tc_bi_ortho
   call routine_diag()
   call write_tc_energy()
   call save_tc_bi_ortho_wavefunction()
-end
-
-subroutine test
- implicit none
- integer :: i,j 
- double precision :: hmono,htwoe,hthree,htot
- use bitmasks
- print*,'reading the wave function '
- do i = 1, N_det
-  call debug_det(psi_det(1,1,i),N_int)
-  print*,i,psi_l_coef_bi_ortho(i,1)*psi_r_coef_bi_ortho(i,1)
-  print*,i,psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1)
- enddo
 
 end
 
+! ---
+
+subroutine test()
+
+  use bitmasks
+  implicit none
+  integer          :: i, j 
+  double precision :: hmono, htwoe, hthree, htot
+
+  print*, 'reading the wave function '
+  do i = 1, N_det
+    call debug_det(psi_det(1,1,i), N_int)
+    print*, i, psi_l_coef_bi_ortho(i,1)*psi_r_coef_bi_ortho(i,1)
+    print*, i, psi_l_coef_bi_ortho(i,1),psi_r_coef_bi_ortho(i,1)
+  enddo
+
+end
+
+! ---
+
 subroutine routine_diag()
 
   implicit none

src/tc_bi_ortho/tc_bi_ortho_prop.irp.f

@@ -1,19 +1,32 @@
+
+! ---
+
 program tc_bi_ortho_prop
-  implicit none
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! 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
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
-!  call routine_diag
+  !call routine_diag
- call test
+  call test
+
 end
 
+! ---
+
 subroutine test
  implicit none
  integer :: i

src/tc_bi_ortho/tc_cisd_sc2.irp.f

@@ -1,20 +1,32 @@
-program tc_bi_ortho
+
-  implicit none
+! ---
+
+program tc_cisd_sc2
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! 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
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   call test
+
 end
 
-subroutine test
+! ---
+
+subroutine test()
  implicit none
 ! double precision, allocatable :: dressing_dets(:),e_corr_dets(:)
 ! allocate(dressing_dets(N_det),e_corr_dets(N_det))

src/tc_bi_ortho/tc_h_eigvectors.irp.f

@@ -1,42 +1,56 @@
+
+! ---
+
   use bitmasks
 
- BEGIN_PROVIDER [ integer, index_HF_psi_det]                                                                                                            
+! ---
- implicit none
+
- integer :: i,degree
+BEGIN_PROVIDER [integer, index_HF_psi_det]                                                                                                            
- do i = 1, N_det
-   call get_excitation_degree(HF_bitmask,psi_det(1,1,i),degree,N_int)
-   if(degree == 0)then
-    index_HF_psi_det = i
-    exit
-   endif
- enddo
- END_PROVIDER
 
-subroutine diagonalize_CI_tc
   implicit none
+  integer :: i, degree
+
+  do i = 1, N_det
+    call get_excitation_degree(HF_bitmask, psi_det(1,1,i), degree, N_int)
+    if(degree == 0) then
+      index_HF_psi_det = i
+      exit
+    endif
+  enddo
+
+END_PROVIDER
+
+! ---
+
+subroutine diagonalize_CI_tc()
+
   BEGIN_DOC
-!  Replace the coefficients of the |CI| states by the coefficients of the
+  !  Replace the coefficients of the |CI| states by the coefficients of the
-!  eigenstates of the |CI| matrix.
+  !  eigenstates of the |CI| matrix.
   END_DOC
-  integer                        :: i,j
+
-  do j=1,N_states
+  implicit none
-    do i=1,N_det
+  integer :: i, j
+
+  do j = 1, N_states
+    do i = 1, N_det
       psi_l_coef_bi_ortho(i,j) = leigvec_tc_bi_orth(i,j)
       psi_r_coef_bi_ortho(i,j) = reigvec_tc_bi_orth(i,j)
     enddo
   enddo
 
   SOFT_TOUCH psi_l_coef_bi_ortho psi_r_coef_bi_ortho
+
 end
 
+! ---
 
-
+ BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states)      ]
- BEGIN_PROVIDER [double precision, eigval_right_tc_bi_orth, (N_states)]
+&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth , (N_states)      ]
-&BEGIN_PROVIDER [double precision, eigval_left_tc_bi_orth, (N_states)]
+&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth     , (N_det,N_states)]
-&BEGIN_PROVIDER [double precision, reigvec_tc_bi_orth, (N_det,N_states)]
+&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth     , (N_det,N_states)]
-&BEGIN_PROVIDER [double precision, leigvec_tc_bi_orth, (N_det,N_states)]
+&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth   , (N_states)      ]
-&BEGIN_PROVIDER [double precision, s2_eigvec_tc_bi_orth, (N_states)]
+&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth           ]
-&BEGIN_PROVIDER [double precision, norm_ground_left_right_bi_orth ]
 
   BEGIN_DOC
   ! eigenvalues, right and left eigenvectors of the transcorrelated Hamiltonian on the BI-ORTHO basis 
@@ -44,29 +58,29 @@ end
 
   implicit none
   integer                       :: i, idx_dress, j, istate, k
+  integer                       :: i_good_state, i_other_state, i_state
+  integer                       :: n_real_tc_bi_orth_eigval_right, igood_r, igood_l
   logical                       :: converged, dagger
-  integer                       :: n_real_tc_bi_orth_eigval_right,igood_r,igood_l
+  double precision, parameter   :: alpha = 0.1d0
-  double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:),leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
+  integer,          allocatable :: index_good_state_array(:)
+  integer,          allocatable :: iorder(:)
+  logical,          allocatable :: good_state_array(:)
+  double precision, allocatable :: reigvec_tc_bi_orth_tmp(:,:), leigvec_tc_bi_orth_tmp(:,:),eigval_right_tmp(:)
   double precision, allocatable :: s2_values_tmp(:), H_prime(:,:), expect_e(:)
-  double precision, parameter :: alpha = 0.1d0
+  double precision, allocatable :: coef_hf_r(:),coef_hf_l(:)
-  integer                        :: i_good_state,i_other_state, i_state
+  double precision, allocatable :: Stmp(:,:)
-  integer, allocatable           :: index_good_state_array(:)
-  logical, allocatable           :: good_state_array(:)
-  double precision, allocatable  :: coef_hf_r(:),coef_hf_l(:)
-  double precision, allocatable  :: Stmp(:,:)
-  integer, allocatable :: iorder(:)
 
   PROVIDE N_det N_int
 
-  if(n_det .le. N_det_max_full) then
+  if(N_det .le. N_det_max_full) then
 
-    allocate(reigvec_tc_bi_orth_tmp(N_det,N_det),leigvec_tc_bi_orth_tmp(N_det,N_det),eigval_right_tmp(N_det),expect_e(N_det))
+    allocate(reigvec_tc_bi_orth_tmp(N_det,N_det), leigvec_tc_bi_orth_tmp(N_det,N_det), eigval_right_tmp(N_det), expect_e(N_det))
-    allocate (H_prime(N_det,N_det),s2_values_tmp(N_det))
+    allocate(H_prime(N_det,N_det), s2_values_tmp(N_det))
 
     H_prime(1:N_det,1:N_det) = htilde_matrix_elmt_bi_ortho(1:N_det,1:N_det)
     if(s2_eig) then
       H_prime(1:N_det,1:N_det) += alpha * S2_matrix_all_dets(1:N_det,1:N_det)
-      do j=1,N_det
+      do j = 1, N_det
         H_prime(j,j) = H_prime(j,j) - alpha*expected_s2
       enddo
     endif

src/tc_bi_ortho/tc_hmat.irp.f

@@ -31,7 +31,9 @@
 
 END_PROVIDER 
 
- BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho_tranp, (N_det,N_det)]
+! ---
+
+BEGIN_PROVIDER [double precision, htilde_matrix_elmt_bi_ortho_tranp, (N_det,N_det)]
  implicit none
  integer ::i,j
   do i = 1, N_det

src/tc_bi_ortho/tc_som.irp.f

@@ -12,10 +12,9 @@ program tc_som
   print *, ' do not forget to do tc-scf first'
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
-!  my_n_pt_r_grid = 10 ! small grid for quick debug
+  my_n_pt_a_grid = tc_grid1_a
-!  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
 
   PROVIDE mu_erf 

src/tc_bi_ortho/test_natorb.irp.f

@@ -1,21 +1,34 @@
+
+! ---
+
 program test_natorb
-  implicit none
+
   BEGIN_DOC
-! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
+  ! TODO : Reads psi_det in the EZFIO folder and prints out the left- and right-eigenvectors together with the energy. Saves the left-right wave functions at the end. 
   END_DOC
+
+  implicit none
+
   print *, 'Hello world'
+
   my_grid_becke = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
- call routine
+  call routine()
-! call test
+  ! call test()
 
 end
 
-subroutine routine
+! ---
+
+subroutine routine()
+
  implicit none
  double precision, allocatable :: fock_diag(:),eigval(:),leigvec(:,:),reigvec(:,:),mat_ref(:,:)
  allocate(eigval(mo_num),leigvec(mo_num,mo_num),reigvec(mo_num,mo_num),fock_diag(mo_num),mat_ref(mo_num, mo_num))

src/tc_bi_ortho/test_normal_order.irp.f

@@ -1,19 +1,32 @@
+
+! ---
+
 program test_normal_order
-  implicit none
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! 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
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
-  call provide_all_three_ints_bi_ortho
+  call provide_all_three_ints_bi_ortho()
-  call test
+  call test()
+
 end
 
+! ---
+
 subroutine test
  implicit none
   use bitmasks ! you need to include the bitmasks_module.f90 features

src/tc_bi_ortho/test_s2_tc.irp.f

@@ -1,14 +1,22 @@
+
+! ---
+
 program test_tc 
- implicit none
+
- read_wf = .True.
+  implicit none
- my_grid_becke = .True.
+ 
- my_n_pt_r_grid = 30
+  my_grid_becke = .True.
- my_n_pt_a_grid = 50
+  PROVIDE tc_grid1_a tc_grid1_r
- read_wf = .True.
+  my_n_pt_r_grid = tc_grid1_r
- touch read_wf
+  my_n_pt_a_grid = tc_grid1_a
- touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
- call routine_test_s2
+
- call routine_test_s2_davidson
+  read_wf = .True.
+  touch read_wf
+
+  call routine_test_s2
+  call routine_test_s2_davidson
+
 end
 
 subroutine routine_test_s2

src/tc_bi_ortho/test_tc_bi_ortho.irp.f

@@ -1,15 +1,24 @@
+
+! ---
+
 program tc_bi_ortho
-  implicit none
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! 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
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
 ! call test_h_u0
 ! call test_slater_tc_opt

src/tc_bi_ortho/test_tc_fock.irp.f

@@ -1,22 +1,32 @@
+
+! ---
+
 program test_tc_fock
- implicit none
+
   BEGIN_DOC
-! TODO : Put the documentation of the program here
+  ! 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
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
+  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
+
   read_wf = .True.
   touch read_wf
-  touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   !call routine_1
   !call routine_2
 !  call routine_3()
 
 ! call test_3e
- call routine_tot
+  call routine_tot
+
 end
 
 ! ---

src/tc_keywords/EZFIO.cfg

@@ -262,3 +262,16 @@ doc: If |true|, use Manu IPP
 interface: ezfio,provider,ocaml
 default: True
 
+[tc_grid1_a]
+type: integer
+doc: size of angular grid over r1
+interface: ezfio,provider,ocaml
+default: 50
+
+[tc_grid1_r]
+type: integer
+doc: size of radial grid over r1
+interface: ezfio,provider,ocaml
+default: 30
+
+

src/tc_scf/combine_lr_tcscf.irp.f

@@ -10,8 +10,9 @@ program combine_lr_tcscf
   implicit none
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   bi_ortho = .True.

src/tc_scf/minimize_tc_angles.irp.f

@@ -1,17 +1,26 @@
-program print_angles
+
- implicit none
+! ---
- BEGIN_DOC
+
-! program that minimizes the angle between left- and right-orbitals when degeneracies are found in the TC-Fock matrix
+program minimize_tc_angles
- END_DOC
+
+  BEGIN_DOC
+  ! program that minimizes the angle between left- and right-orbitals when degeneracies are found in the TC-Fock matrix
+  END_DOC
+
+  implicit none
+
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
   touch my_n_pt_r_grid my_n_pt_a_grid
-!  call sort_by_tc_fock
+
+  !  call sort_by_tc_fock
 
   ! TODO
   ! check if rotations of orbitals affect the TC energy
   ! and refuse the step
   call minimize_tc_orb_angles
+
 end
 

src/tc_scf/molden_lr_mos.irp.f

@@ -11,10 +11,9 @@ program molden_lr_mos
   print *, 'starting ...'
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
-!  my_n_pt_r_grid = 10 ! small grid for quick debug
+  my_n_pt_a_grid = tc_grid1_a
-!  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
 
   !call molden_lr

src/tc_scf/print_fit_param.irp.f

@@ -7,10 +7,9 @@ program print_fit_param
   implicit none
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
-!  my_n_pt_r_grid = 10 ! small grid for quick debug
+  my_n_pt_a_grid = tc_grid1_a
-!  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
 
   !call create_guess

src/tc_scf/print_tcscf_energy.irp.f

@@ -8,16 +8,9 @@ program print_tcscf_energy
 
   print *, 'Hello world'
   my_grid_becke = .True.
-
+  PROVIDE tc_grid1_a tc_grid1_r
-  !my_n_pt_r_grid = 30
+  my_n_pt_r_grid = tc_grid1_r
-  !my_n_pt_a_grid = 50
+  my_n_pt_a_grid = tc_grid1_a
-
-  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()

src/tc_scf/rotate_tcscf_orbitals.irp.f

@@ -10,8 +10,9 @@ program rotate_tcscf_orbitals
   implicit none
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
 
   bi_ortho = .True.

src/tc_scf/tc_petermann_factor.irp.f

@@ -10,10 +10,9 @@ program tc_petermann_factor
   implicit none
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
-!  my_n_pt_r_grid = 10 ! small grid for quick debug
+  my_n_pt_a_grid = tc_grid1_a
-!  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
 
   call main()

src/tc_scf/tc_scf.irp.f

@@ -14,14 +14,10 @@ program tc_scf
 
   my_grid_becke  = .True.
 
-  !my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  !my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
+  my_n_pt_a_grid = tc_grid1_a
 
-  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
 
   PROVIDE mu_erf 

src/tc_scf/test_int.irp.f

@@ -9,10 +9,9 @@ program test_ints
   print *, ' starting test_ints ...'
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  PROVIDE tc_grid1_a tc_grid1_r
-  my_n_pt_a_grid = 50
+  my_n_pt_r_grid = tc_grid1_r
-!  my_n_pt_r_grid = 15 ! small grid for quick debug
+  my_n_pt_a_grid = tc_grid1_a
-!  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
 
   my_extra_grid_becke = .True.
@@ -280,7 +279,7 @@ subroutine routine_v_ij_u_cst_mu_j1b_test
     do i = 1, ao_num
      do j = 1, ao_num
       array(j,i,l,k)     += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
+      array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b_fit (j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
      enddo
     enddo
    enddo
@@ -506,7 +505,7 @@ subroutine routine_v_ij_u_cst_mu_j1b
     do i = 1, ao_num
      do j = 1, ao_num
       array(j,i,l,k)     += v_ij_u_cst_mu_j1b_test(j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
-      array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b(j,i,ipoint)      * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
+      array_ref(j,i,l,k) += v_ij_u_cst_mu_j1b_fit (j,i,ipoint) * aos_in_r_array(k,ipoint) * aos_in_r_array(l,ipoint) * weight
      enddo
     enddo
    enddo