diff --git a/external/irpf90 b/external/irpf90
index 0007f72f..4ab1b175 160000
--- a/external/irpf90
+++ b/external/irpf90
@@ -1 +1 @@
-Subproject commit 0007f72f677fe7d61c5e1ed461882cb239517102
+Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6
diff --git a/src/ao_two_e_ints/providers_ao_erf.irp.f b/src/ao_two_e_ints/providers_ao_erf.irp.f
index 293df29f..ff8c31a2 100644
--- a/src/ao_two_e_ints/providers_ao_erf.irp.f
+++ b/src/ao_two_e_ints/providers_ao_erf.irp.f
@@ -90,7 +90,7 @@ BEGIN_PROVIDER [ logical, ao_two_e_integrals_erf_in_map ]
   if (write_ao_two_e_integrals_erf) then
     call ezfio_set_work_empty(.False.)
     call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
-    call ezfio_set_ao_two_e_erf_ints_io_ao_two_e_integrals_erf("Read")
+    call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals_erf('Read')
   endif
 
 END_PROVIDER
diff --git a/src/ao_two_e_ints/routines_save_integrals_erf.irp.f b/src/ao_two_e_ints/routines_save_integrals_erf.irp.f
index 4b0cfad0..d980bc05 100644
--- a/src/ao_two_e_ints/routines_save_integrals_erf.irp.f
+++ b/src/ao_two_e_ints/routines_save_integrals_erf.irp.f
@@ -4,7 +4,7 @@ subroutine save_erf_two_e_integrals_ao
  PROVIDE ao_two_e_integrals_erf_in_map
  call ezfio_set_work_empty(.False.)
  call map_save_to_disk(trim(ezfio_filename)//'/work/ao_ints_erf',ao_integrals_erf_map)
- call ezfio_set_ao_two_e_erf_ints_io_ao_two_e_integrals_erf('Read')
+ call ezfio_set_ao_two_e_ints_io_ao_two_e_integrals_erf('Read')
 end
 
 subroutine save_erf_two_e_ints_ao_into_ints_ao
diff --git a/src/basis_correction/51.basis_c.bats b/src/basis_correction/51.basis_c.bats
index 2682361b..914b482b 100644
--- a/src/basis_correction/51.basis_c.bats
+++ b/src/basis_correction/51.basis_c.bats
@@ -10,8 +10,8 @@ function run() {
   qp set perturbation do_pt2 False
   qp set determinants n_det_max 8000
   qp set determinants n_states  1
-  qp set davidson threshold_davidson 1.e-10
-  qp set davidson n_states_diag 8
+  qp set davidson_keywords threshold_davidson 1.e-10
+  qp set davidson_keywords n_states_diag 8
   qp run fci 
   energy1="$(ezfio get fci energy | tr '[]' ' ' | cut -d ',' -f 1)"
   eq $energy1 $1 $thresh
diff --git a/src/casscf_cipsi/50.casscf.bats b/src/casscf_cipsi/50.casscf.bats
index a0db725d..9f63dfe2 100644
--- a/src/casscf_cipsi/50.casscf.bats
+++ b/src/casscf_cipsi/50.casscf.bats
@@ -9,8 +9,8 @@ function run_stoch() {
   test_exe casscf || skip
   qp set perturbation do_pt2 True
   qp set determinants n_det_max $3
-  qp set davidson threshold_davidson 1.e-10
-  qp set davidson n_states_diag 4
+  qp set davidson_keywords threshold_davidson 1.e-10
+  qp set davidson_keywords n_states_diag 4
   qp run casscf | tee casscf.out 
   energy1="$(ezfio get casscf energy_pt2 | tr '[]' ' ' | cut -d ',' -f 1)"
   eq $energy1 $1 $thresh
diff --git a/src/cis/20.cis.bats b/src/cis/20.cis.bats
index 4f255c7b..4a5c6e45 100644
--- a/src/cis/20.cis.bats
+++ b/src/cis/20.cis.bats
@@ -9,7 +9,7 @@ function run() {
   qp set_file $1
   qp edit --check
   qp set determinants n_states  3
-  qp set davidson threshold_davidson 1.e-12
+  qp set davidson_keywords threshold_davidson 1.e-12
   qp set mo_two_e_ints io_mo_two_e_integrals Write
   qp set_frozen_core 
   qp run cis 
@@ -59,7 +59,7 @@ function run() {
 
 @test "ClO" { # 1.65582s 2.06465s
   [[ -n $TRAVIS ]] && skip
-  run clo.ezfio -534.263560525680 -534.256601571199 -534.062020844428
+  run clo.ezfio -534.2635737789097 -534.2566081298855 -534.0620070783308
 }
 
 @test "SO" { # 1.9667s 2.91234s
@@ -69,7 +69,7 @@ function run() {
 
 @test "OH" { # 2.201s 2.65573s
   [[ -n $TRAVIS ]] && skip
-  run oh.ezfio -75.4314648243896 -75.4254639668256 -75.2707675632313
+  run oh.ezfio -75.4314822573358 -75.4254733392003 -75.2707586997333
 }
 
 @test "H2O2" { # 2.27079s 3.07875s
@@ -109,7 +109,7 @@ function run() {
 
 @test "DHNO" { # 6.42976s 12.9899s
   [[ -n $TRAVIS ]] && skip
-  run dhno.ezfio -130.4472288472718 -130.3571808164850 -130.2196257046987
+  run dhno.ezfio -130.447238897118 -130.357186843611 -130.219626716369
 }
 
 @test "CH4" { # 6.4969s 10.9157s
@@ -129,7 +129,7 @@ function run() {
 
 @test "[Cu(NH3)4]2+" { # 29.7711s  3.45478m
   [[ -n ${TRAVIS} ]] && skip
-  run  cu_nh3_4_2plus.ezfio -1862.97958885180 -1862.92457657404 -1862.91134959451
+  run  cu_nh3_4_2plus.ezfio -1862.97958844302 -1862.92454785007 -1862.91130869967
 
 }
 
diff --git a/src/cisd/30.cisd.bats b/src/cisd/30.cisd.bats
index 6b8fddb6..5ec11e4b 100644
--- a/src/cisd/30.cisd.bats
+++ b/src/cisd/30.cisd.bats
@@ -8,10 +8,9 @@ function run() {
   test_exe cisd || skip
   qp edit --check
   qp set determinants n_states  2
-  qp set davidson threshold_davidson 1.e-12
-  qp set davidson n_states_diag 24
-  qp run cis
-  qp run cisd 
+  qp set davidson_keywords threshold_davidson 1.e-12
+  qp set davidson_keywords n_states_diag 24
+  qp run cisd
   energy1="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 1)"
   energy2="$(qp get cisd energy | tr '[]' ' ' | cut -d ',' -f 2)"
   eq $energy1 $1 $thresh
@@ -19,7 +18,7 @@ function run() {
 }
 
 
-@test "B-B" { # 
+@test "B-B" { #
   qp set_file b2_stretched.ezfio
   qp set_frozen_core
   run  -49.120607088648597 -49.055152453388231
@@ -34,7 +33,7 @@ function run() {
 @test "HBO" { # 4.42968s  19.6099s
   qp set_file hbo.ezfio
   qp set_frozen_core
-  run -100.2019254455993 -99.79484127741013 
+  run -100.2019254455993 -99.79484127741013
 }
 
 @test "HCO" { # 6.6077s 28.6801s
@@ -46,7 +45,7 @@ function run() {
 @test "H2O" { # 7.0651s 30.6642s
   qp set_file h2o.ezfio
   qp set_frozen_core
-  run -76.22975602077072 -75.80609108747208 
+  run -76.22975602077072 -75.80609108747208
 }
 
 
@@ -78,7 +77,7 @@ function run() {
   [[ -n $TRAVIS ]] && skip
   qp set_file oh.ezfio
   qp set_frozen_core
-  run -75.6087472926588 -75.5370393736601
+  run -75.6088105201621  -75.5370802925698
 }
 
 @test "CH4" { # 19.821s 1.38648m
@@ -105,8 +104,9 @@ function run() {
 @test "DHNO" { # 24.7077s 1.46487m
   [[ -n $TRAVIS ]] && skip
   qp set_file dhno.ezfio
-  qp set_mo_class --core="[1-7]" --act="[8-64]" 
-  run -130.458814562403 -130.356308303681
+  qp set_mo_class --core="[1-7]" --act="[8-64]"
+  run -130.4659881027444 -130.2692384198501
+#  run -130.458814562403 -130.356308303681
 }
 
 @test "H3COH" { # 24.7248s 1.85043m
@@ -120,7 +120,7 @@ function run() {
   [[ -n $TRAVIS ]] && skip
   qp set_file cu_nh3_4_2plus.ezfio
   qp set_mo_class --core="[1-24]" --act="[25-45]" --del="[46-87]"
-  run -1862.98689579931  -1862.6883044626563
+  run -1862.98310702274  -1862.88506319755
 
 }
 
@@ -135,14 +135,14 @@ function run() {
   [[ -n $TRAVIS ]] && skip
   qp set_file c2h2.ezfio
   qp set_mo_class --act="[1-30]" --del="[31-36]"
-  run -12.3566731164213 -11.9495394759914 
+  run -12.3566731164213 -11.9495394759914
 }
 
 @test "ClO" { # 37.6949s
   [[ -n $TRAVIS ]] && skip
   qp set_file clo.ezfio
   qp set_frozen_core
-  run -534.5404021326773 -534.3818725793897 
+  run -534.540464615019  -534.381904487587
 }
 
 @test "F2" { # 45.2078s
@@ -155,7 +155,7 @@ function run() {
 @test "SO2" { # 47.6922s
   [[ -n $TRAVIS ]] && skip
   qp set_file so2.ezfio
-  qp set_mo_class --core="[1-8]" --act="[9-87]" 
+  qp set_mo_class --core="[1-8]" --act="[9-87]"
   run -41.5746738710350 -41.3800467740750
 }
 
@@ -177,7 +177,7 @@ function run() {
   [[ -n $TRAVIS ]] && skip
   qp set_file n2.ezfio
   qp set_mo_class --core="[1,2]" --act="[3-40]" --del="[41-60]"
-  run -109.275693633982 -108.757794570948 
+  run -109.275693633982 -108.757794570948
 }
 
 @test "HCN" { # 133.8696s
diff --git a/src/dft_one_e/mu_erf_dft.irp.f b/src/dft_one_e/mu_erf_dft.irp.f
index 0b870564..08779f0e 100644
--- a/src/dft_one_e/mu_erf_dft.irp.f
+++ b/src/dft_one_e/mu_erf_dft.irp.f
@@ -3,7 +3,7 @@ BEGIN_PROVIDER [double precision, mu_erf_dft]
  BEGIN_DOC
 ! range separation parameter used in RS-DFT. 
 ! 
-! It is set to mu_erf in order to be consistent with the module "ao_two_e_erf_ints"
+! It is set to mu_erf in order to be consistent with the module "hamiltonian"
  END_DOC
  mu_erf_dft = mu_erf
 
diff --git a/src/kohn_sham_rs/61.rsks.bats b/src/kohn_sham_rs/61.rsks.bats
index 90b82142..29d43c3b 100644
--- a/src/kohn_sham_rs/61.rsks.bats
+++ b/src/kohn_sham_rs/61.rsks.bats
@@ -13,7 +13,7 @@ function run() {
   qp set scf_utils thresh_scf 1.e-10
   qp set dft_keywords exchange_functional $functional
   qp set dft_keywords correlation_functional $functional
-  qp set ao_two_e_erf_ints mu_erf 0.5 
+  qp set hamiltonian mu_erf 0.5 
   qp set becke_numerical_grid grid_type_sgn 1 
   qp_reset --mos $1 
   qp run rs_ks_scf 
diff --git a/src/mo_two_e_ints/mo_bi_integrals_erf.irp.f b/src/mo_two_e_ints/mo_bi_integrals_erf.irp.f
index e009b7d9..1afc1f3c 100644
--- a/src/mo_two_e_ints/mo_bi_integrals_erf.irp.f
+++ b/src/mo_two_e_ints/mo_bi_integrals_erf.irp.f
@@ -9,25 +9,27 @@ subroutine mo_two_e_integrals_erf_index(i,j,k,l,i1)
   integer(key_kind)              :: p,q,r,s,i2
   p = min(i,k)
   r = max(i,k)
-  p = p+ishft(r*r-r,-1)
+  p = p+shiftr(r*r-r,1)
   q = min(j,l)
   s = max(j,l)
-  q = q+ishft(s*s-s,-1)
+  q = q+shiftr(s*s-s,1)
   i1 = min(p,q)
   i2 = max(p,q)
-  i1 = i1+ishft(i2*i2-i2,-1)
+  i1 = i1+shiftr(i2*i2-i2,1)
 end
 
 
 BEGIN_PROVIDER [ logical, mo_two_e_integrals_erf_in_map ]
   use map_module
   implicit none
-  integer(bit_kind)              :: mask_ijkl(N_int,4)
-  integer(bit_kind)              :: mask_ijk(N_int,3)
-
   BEGIN_DOC
   ! If True, the map of MO two-electron integrals is provided
   END_DOC
+  integer(bit_kind)              :: mask_ijkl(N_int,4)
+  integer(bit_kind)              :: mask_ijk(N_int,3)
+  double precision               :: cpu_1, cpu_2, wall_1, wall_2
+
+  PROVIDE mo_class
 
   real                           :: map_mb
 
@@ -55,7 +57,7 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_erf_in_map ]
   if (write_mo_two_e_integrals_erf) then
     call ezfio_set_work_empty(.False.)
     call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
-    call ezfio_set_mo_two_e_erf_ints_io_mo_two_e_integrals_erf("Read")
+    call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals_erf("Read")
   endif
 
 END_PROVIDER
diff --git a/src/mo_two_e_ints/routines_save_integrals_erf.irp.f b/src/mo_two_e_ints/routines_save_integrals_erf.irp.f
index 52fb8f63..9915b206 100644
--- a/src/mo_two_e_ints/routines_save_integrals_erf.irp.f
+++ b/src/mo_two_e_ints/routines_save_integrals_erf.irp.f
@@ -4,7 +4,7 @@ subroutine save_erf_two_e_integrals_mo
  PROVIDE mo_two_e_integrals_erf_in_map
  call ezfio_set_work_empty(.False.)
  call map_save_to_disk(trim(ezfio_filename)//'/work/mo_ints_erf',mo_integrals_erf_map)
- call ezfio_set_mo_two_e_erf_ints_io_mo_two_e_integrals_erf('Read')
+ call ezfio_set_mo_two_e_ints_io_mo_two_e_integrals_erf('Read')
 end
 
 
diff --git a/src/mu_of_r/mu_of_r_conditions.irp.f b/src/mu_of_r/mu_of_r_conditions.irp.f
index f9c3b3b3..959950a6 100644
--- a/src/mu_of_r/mu_of_r_conditions.irp.f
+++ b/src/mu_of_r/mu_of_r_conditions.irp.f
@@ -128,7 +128,7 @@
  BEGIN_PROVIDER [double precision, mu_average_prov, (N_states)]
  implicit none
  BEGIN_DOC
- ! average value of mu(r) weighted with the total one-e density and divised by the number of electrons 
+ ! average value of mu(r) weighted with the total one-e density and divided by the number of electrons 
  !
  ! !!!!!! WARNING !!!!!! if no_core_density == .True. then all contributions from the core orbitals 
  !
diff --git a/src/non_h_ints_mu/tc_integ_num.irp.f b/src/non_h_ints_mu/tc_integ_num.irp.f
index ee34f531..5a088331 100644
--- a/src/non_h_ints_mu/tc_integ_num.irp.f
+++ b/src/non_h_ints_mu/tc_integ_num.irp.f
@@ -47,7 +47,7 @@
   call total_memory(mem)
   mem      = max(1.d0, qp_max_mem - mem)
   n_double = mem * 1.d8
-  n_blocks = min(n_double / (n_points_extra_final_grid * 4), 1.d0*n_points_final_grid)
+  n_blocks = int(min(n_double / (n_points_extra_final_grid * 4.d0), 1.d0*n_points_final_grid))
   n_rest   = int(mod(n_points_final_grid, n_blocks))
   n_pass   = int((n_points_final_grid - n_rest) / n_blocks)
 
diff --git a/src/tc_scf/11.tc_scf.bats b/src/tc_scf/11.tc_scf.bats
index 91b52540..b81c2f4b 100644
--- a/src/tc_scf/11.tc_scf.bats
+++ b/src/tc_scf/11.tc_scf.bats
@@ -8,15 +8,15 @@ function run_Ne() {
   rm -rf Ne_tc_scf
   echo Ne > Ne.xyz
   qp create_ezfio -b cc-pcvdz Ne.xyz -o Ne_tc_scf
-  qp run scf 
+  qp run scf
 
-  qp set ao_two_e_erf_ints mu_erf 0.87
+  qp set hamiltonian mu_erf 0.87
   qp set tc_keywords j1b_type 3
-  qp set tc_keywords j1b_pen [1.5] 
-  qp set tc_keywords bi_ortho True 
+  qp set tc_keywords j1b_pen [1.5]
+  qp set tc_keywords bi_ortho True
   qp set tc_keywords test_cycle_tc True
 
-  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out 
+  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
   eref=-128.552134
   energy="$(qp get tc_scf bitc_energy)"
   eq $energy $eref 1e-6
@@ -24,22 +24,22 @@ function run_Ne() {
 
 
 @test "Ne" {
- run_Ne 
+ run_Ne
 }
 
 function run_C() {
   rm -rf C_tc_scf
   echo C  > C.xyz
   qp create_ezfio -b cc-pcvdz C.xyz -o C_tc_scf -m 3
-  qp run scf 
+  qp run scf
 
-  qp set ao_two_e_erf_ints mu_erf 0.87
+  qp set hamiltonian mu_erf 0.87
   qp set tc_keywords j1b_type 3
-  qp set tc_keywords j1b_pen [1.5] 
-  qp set tc_keywords bi_ortho True 
+  qp set tc_keywords j1b_pen [1.5]
+  qp set tc_keywords bi_ortho True
   qp set tc_keywords test_cycle_tc True
 
-  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out 
+  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
   eref=-37.691254356408791
   energy="$(qp get tc_scf bitc_energy)"
   eq $energy $eref 1e-6
@@ -47,7 +47,7 @@ function run_C() {
 
 
 @test "C" {
- run_C 
+ run_C
 }
 
 
@@ -55,15 +55,15 @@ function run_O() {
   rm -rf O_tc_scf
   echo O  > O.xyz
   qp create_ezfio -b cc-pcvdz O.xyz -o O_tc_scf -m 3
-  qp run scf 
+  qp run scf
 
-  qp set ao_two_e_erf_ints mu_erf 0.87
+  qp set hamiltonian mu_erf 0.87
   qp set tc_keywords j1b_type 3
-  qp set tc_keywords j1b_pen [1.5] 
-  qp set tc_keywords bi_ortho True 
+  qp set tc_keywords j1b_pen [1.5]
+  qp set tc_keywords bi_ortho True
   qp set tc_keywords test_cycle_tc True
 
-  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out 
+  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
   eref=-74.814687229354590
   energy="$(qp get tc_scf bitc_energy)"
   eq $energy $eref 1e-6
@@ -71,7 +71,7 @@ function run_O() {
 
 
 @test "O" {
- run_O 
+ run_O
 }
 
 
@@ -79,16 +79,16 @@ function run_O() {
 function run_ch2() {
   rm -rf ch2_tc_scf
   cp ${QP_ROOT}/tests/input/ch2.xyz .
-  qp create_ezfio -b "C:cc-pcvdz|H:cc-pvdz" ch2.xyz -o ch2_tc_scf 
-  qp run scf 
+  qp create_ezfio -b "C:cc-pcvdz|H:cc-pvdz" ch2.xyz -o ch2_tc_scf
+  qp run scf
 
-  qp set ao_two_e_erf_ints mu_erf 0.87
+  qp set hamiltonian mu_erf 0.87
   qp set tc_keywords j1b_type 3
   qp set tc_keywords j1b_pen '[1.5,10000,10000]'
   qp set tc_keywords bi_ortho True
   qp set tc_keywords test_cycle_tc True
 
-  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out 
+  qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
   eref=-38.903247818077737
   energy="$(qp get tc_scf bitc_energy)"
   eq $energy $eref 1e-6
@@ -96,6 +96,6 @@ function run_ch2() {
 
 
 @test "ch2" {
- run_ch2 
+ run_ch2
 }
 
diff --git a/src/tc_scf/fock_three_hermit.irp.f b/src/tc_scf/fock_three_hermit.irp.f
index 6c132189..00d47fae 100644
--- a/src/tc_scf/fock_three_hermit.irp.f
+++ b/src/tc_scf/fock_three_hermit.irp.f
@@ -95,7 +95,12 @@ BEGIN_PROVIDER [double precision, diag_three_elem_hf]
 
   if(.not. three_body_h_tc) then
 
-    diag_three_elem_hf = 0.d0
+   if(noL_standard) then
+      PROVIDE noL_0e
+      diag_three_elem_hf = noL_0e
+    else
+      diag_three_elem_hf = 0.d0
+    endif
 
   else
 
diff --git a/src/tc_scf/rh_tcscf_diis.irp.f b/src/tc_scf/rh_tcscf_diis.irp.f
index c7f35451..12678500 100644
--- a/src/tc_scf/rh_tcscf_diis.irp.f
+++ b/src/tc_scf/rh_tcscf_diis.irp.f
@@ -71,10 +71,7 @@ subroutine rh_tcscf_diis()
   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
-  endif
+  etc_3e  = diag_three_elem_hf
   !tc_grad = grad_non_hermit
   er_DIIS = maxval(abs(FQS_SQF_mo))
   e_delta = dabs(etc_tot - e_save)
@@ -202,10 +199,7 @@ subroutine rh_tcscf_diis()
     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
-    endif
+    etc_3e  = diag_three_elem_hf
     !tc_grad  = grad_non_hermit
     er_DIIS  = maxval(abs(FQS_SQF_mo))
     e_delta  = dabs(etc_tot - e_save)