Merge branch 'QuantumPackage-dev-stable-tc-scf' into dev-stable-tc-scf

.travis.yml

@@ -1,52 +0,0 @@
-#sudo: true
-#before_script:
-#  - sudo apt-get update -q
-#  - sudo apt-get remove curl
-#  - sudo apt-get remove zlib1g-dev
-#  - sudo apt-get install autoconf
-#  - sudo rm /usr/local/bin/bats
-
-os: linux
-
-dist: bionic
-
-sudo: false
-
-compiler: gfortran
-
-addons:
-  apt:
-    packages:
-    - gfortran
-    - gcc
-    - libatlas-base-dev
-#    - liblapack-dev
-#    - libblas-dev
-    - wget
-
-env:
-  - OPAMROOT=$HOME/.opam
-
-cache:
-  directories:
-  - $HOME/.opam/
-  - $HOME/cache
-
-language: python
-python:
-    - "3.7"
-
-stages:
-  - configuration
-  - compilation
-  - testing
-
-jobs:
-  include:
-    - stage: configuration
-      script: travis/configuration.sh
-    - stage: compilation
-      script: travis/compilation.sh
-    - stage: testing
-      script: travis/testing.sh
-

RELEASE_NOTES.org

@@ -9,15 +9,23 @@
   - Configure adapted for ARM
   - Added many types of integrals
   - Accelerated four-index transformation
-
-*** TODO: take from dev
-  - [ ] Added GTOs with complex exponent
-  - Updated version of f77-zmq
   - Added transcorrelated SCF
   - Added transcorrelated CIPSI
+  - Added CCSD and CCSD(T)
+  - Added MO localization
+  - Changed coupling parameters for ROHF
+  - General Davidson algorithm
+  - Accelerated restore_symmetry
+  - Point charges in the Hamiltonian
+  - Removed cryptokit dependency in OCaml
+  - Using now standard convention in RDM
+  - Added molecular properties
+  - [ ] Added GTOs with complex exponent
+
+*** TODO: take from dev
+  - Updated version of f77-zmq
   - Started to introduce shells in AOs
   - Added ECMD UEG functional
-  - General Davidson algorithm
 
 * Version 2.2
 

scripts/Hn.py

@@ -0,0 +1,43 @@
+#!/usr/bin/env python                                                                                                     
+import sys
+from math import *
+arg = sys.argv
+#f = open('data_dft','r')
+n = int(sys.argv[1])
+r = float(sys.argv[2])
+f = open('H'+str(n)+'_'+str(r)+'.xyz','w')
+string=str(n)+"\n"
+f.write(string)
+string="\n"
+f.write(string)
+for i in range(n):
+ x = r * cos(2.* i* pi/n)
+ y = r * sin(2.* i* pi/n)
+ z = 0.
+ string="H   "+str(x)+"  "+str(y)+"  "+str(z)+"\n"
+ f.write(string)
+ 
+#lines = f.readlines()
+#cipsi_dft= []
+#
+#dissoc = []
+#dissoc.append(float(-76.0179223470363))
+#dissoc.append(float(-76.0592367866993))
+#dissoc.append(float(-76.0678739715659))
+#delta_e = []
+#
+#for line in lines:
+#  data = line.split()
+#  if(len(data)>0):
+#   dft=float(data[1])
+#   fci=float(data[2])
+#   e=fci+dft
+#   cipsi_dft.append(e)
+#
+#print(*cipsi_dft,sep="  &  ")
+#
+#for i in 0,1,2:
+# delta_e.append(1000.*(dissoc[i] - cipsi_dft[i]))
+#
+#print(*delta_e,sep="  &  ")
+#

scripts/ezfio_interface/ezfio_generate_provider.py

@@ -52,7 +52,7 @@ BEGIN_PROVIDER [ %(type)s, %(name)s %(size)s ]
     %(test_null_size)s
     call ezfio_has_%(ezfio_dir)s_%(ezfio_name)s(has)
     if (has) then
-      write(6,'(A)') '.. >>>>> [ IO READ: %(name)s ] <<<<< ..'
+!      write(6,'(A)') '.. >>>>> [ IO READ: %(name)s ] <<<<< ..'
       call ezfio_get_%(ezfio_dir)s_%(ezfio_name)s(%(name)s)
     else
       print *, '%(ezfio_dir)s/%(ezfio_name)s not found in EZFIO file'
@@ -117,7 +117,7 @@ END_PROVIDER
         output = self.output
         name = self.name
         l_write = ["",
-                   "  call write_time(%(output)s)",
+                   "!  call write_time(%(output)s)",
                    ""]
 
         self.write = "\n".join(l_write) % locals()
@@ -129,7 +129,7 @@ END_PROVIDER
                 write = self.write_correspondance[self.type]
 
                 l_write = ["",
-                          "  call write_time(%(output)s)",
+                          "!  call write_time(%(output)s)",
                           "  call %(write)s(%(output)s, %(name)s, &",
                           "     '%(name)s')",
                           ""]

scripts/get_fci_conv.sh

@@ -0,0 +1,7 @@
+file=$1
+grep "N_det             =" $1 | cut -d "=" -f 2 > N_det_tmp
+grep "E               =" $file | cut -d "=" -f 2 > E_tmp
+grep "E+PT2           =" $file | cut -d "=" -f 2 | cut -d "+" -f 1 > E+PT2_tmp
+grep "E+rPT2          =" $file | cut -d "=" -f 2 | cut -d "+" -f 1 > E+rPT2_tmp
+paste N_det_tmp E_tmp E+PT2_tmp E+rPT2_tmp | column -s ' ' -t > $file.conv_fci
+rm N_det_tmp E_tmp E+PT2_tmp E+rPT2_tmp

scripts/get_fci_tc_conv.sh

@@ -0,0 +1,2 @@
+file=$1
+grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" $file | cut -d "=" -f 2 > ${file}.conv_fci_tc

scripts/qp_extract_cipsi_data.py

@@ -0,0 +1,54 @@
+#!/usr/bin/env python3
+
+import re
+import sys
+
+# Read output file
+with open(sys.argv[1], 'r') as file:
+  output = file.read()
+
+
+def extract_data(output):
+    lines = output.split("\n")
+    data = []
+
+    n_det = None
+    e = None
+    pt2 = None
+    err_pt2 = None
+    rpt2 = None
+    err_rpt2 = None
+    e_ex = None
+
+
+    reading = False
+    for iline, line in enumerate(lines):
+        if not reading and line.startswith(" N_det "):
+            n_det = int(re.search(r"N_det\s+=\s+(\d+)", line).group(1))
+            reading = True
+
+        if reading:
+            if line.startswith(" E "):
+                e = float(re.search(r"E\s+=\s+(-?\d+\.\d+)", line).group(1))
+            elif line.startswith(" PT2 "):
+                pt2 = float(re.search(r"PT2\s+=\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+                err_pt2 = float(re.search(r"\+/-\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+            elif line.startswith(" rPT2 "):
+                rpt2 = float(re.search(r"rPT2\s+=\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+                err_rpt2 = float(re.search(r"\+/-\s+(-?\d+\.\d+E?.\d*)", line).group(1))
+            elif "minimum PT2 Extrapolated energy" in line:
+                e_ex_line = lines[iline+2]
+                e_ex = float(e_ex_line.split()[1])
+                reading = False
+
+                data.append((n_det, e, pt2, err_pt2, rpt2, err_rpt2, e_ex))
+                n_det = e = pt2 = err_pt2 = rpt2 = err_rpt2 = e_ex = None
+
+    return data
+
+data = extract_data(output)
+
+for item in data:
+    print(" ".join(str(x) for x in item))
+
+

scripts/script_fci_tc.sh

@@ -0,0 +1,33 @@
+ source ~/qp2/quantum_package.rc  
+ alpha=1.8
+ input=O
+ basis=cc-pvdz
+ mult=3
+ output=${input}_${basis}_al_${alpha}
+ qp create_ezfio -b ${basis} ${input}.xyz -m $mult
+ qp run scf 
+ qp set perturbation pt2_max 0.0001
+ qp set_frozen_core 
+
+########## FCI CALCULATION FOR REFERENCE 
+ qp run fci | tee ${EZFIO_FILE}.fci.out
+ qp run sort_wf 
+ mv ${EZFIO_FILE}.wf_sorted ${EZFIO_FILE}_fci.wf_sorted
+########### TC SCF CALCULATION
+ qp reset -d 
+ qp set ao_two_e_erf_ints mu_erf 0.87 
+ qp set tc_keywords j1b_type 3
+ qp set tc_keywords j1b_pen "[${alpha}]"
+ qp set tc_keywords bi_ortho True 
+ qp set tc_keywords test_cycle_tc True
+ qp set tc_keywords write_tc_integ True
+ qp set tc_keywords read_tc_integ False
+ qp run tc_scf | tee ${EZFIO_FILE}.tc_scf.out
+ qp set tc_keywords write_tc_integ False
+ qp set tc_keywords read_tc_integ True
+############ TC-FCI CALCULATION
+ qp run fci_tc_bi_ortho | tee ${EZFIO_FILE}.fci_tc_bi_ortho.out
+ grep "Ndet,E,E+PT2,E+RPT2,|PT2|=" ${EZFIO_FILE}.fci_tc_bi_ortho.out | cut -d "=" -f 2 > data_al_$alpha
+ qp run sort_wf 
+ mv ${EZFIO_FILE}.wf_sorted ${EZFIO_FILE}_tc_fci.wf_sorted
+ 

src/ao_many_one_e_ints/grad2_jmu_manu.irp.f

@@ -38,7 +38,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
      !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,                                             &
      !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,                                           &
      !$OMP          List_comb_thr_b3_cent, int2_grad1u2_grad2u2_j1b2_test, ao_abs_comb_b3_j1b,              &
-     !$OMP          ao_overlap_abs,sq_pi_3_2)
+     !$OMP          ao_overlap_abs,sq_pi_3_2,thrsh_cycle_tc)
  !$OMP DO SCHEDULE(dynamic)
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -46,7 +46,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
     r(3) = final_grid_points(3,ipoint)
     do i = 1, ao_num
       do j = i, ao_num
-        if(ao_overlap_abs(j,i) .lt. 1.d-12) then
+        if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then
           cycle
         endif
 
@@ -58,7 +58,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_1_erf_x_2(i_fit)
           coef_fit = -0.25d0 *  coef_gauss_1_erf_x_2(i_fit)
-          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle
+!          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
           int_gauss = overlap_gauss_r12_ao(r, expo_fit, i, j)
           int2_grad1u2_grad2u2_j1b2_test(j,i,ipoint) += coef_fit * int_gauss 
         enddo
@@ -81,8 +81,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test, (ao_num, ao_n
             !DIR$ FORCEINLINE
             call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
             coef_fit = -0.25d0 *  coef_gauss_1_erf_x_2(i_fit) * coef
-!           if(dabs(coef_fit*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version
+!            if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
-            if(dabs(coef_fit*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
 !           call overlap_gauss_r12_ao_with1s_v(B_center, beta, final_grid_points_transp, &
 !                 expo_fit, i, j, int_fit_v, n_points_final_grid)
             int_gauss = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
@@ -145,14 +144,14 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
      !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,      &
      !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,    &
      !$OMP          List_comb_thr_b3_cent, big_array,&
-     !$OMP          ao_abs_comb_b3_j1b,ao_overlap_abs)
+     !$OMP          ao_abs_comb_b3_j1b,ao_overlap_abs,thrsh_cycle_tc)
 !
  allocate(int_fit_v(n_points_final_grid))
  !$OMP DO SCHEDULE(dynamic)
  do i = 1, ao_num
    do j = i, ao_num
 
-     if(ao_overlap_abs(j,i) .lt. 1.d-12) then
+     if(ao_overlap_abs(j,i) .lt. thrsh_cycle_tc) then
        cycle
      endif
 
@@ -161,7 +160,6 @@ BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2_test_v, (ao_num, ao
          coef        = List_comb_thr_b3_coef  (i_1s,j,i)
          beta        = List_comb_thr_b3_expo  (i_1s,j,i)
          int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-!         if(dabs(coef)*dabs(int_j1b).lt.1.d-15)cycle
          B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
          B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
          B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -243,7 +241,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
  !$OMP          final_grid_points, ng_fit_jast,                     &
  !$OMP          expo_gauss_j_mu_x_2, coef_gauss_j_mu_x_2,           &
  !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,sq_pi_3_2,       & 
- !$OMP          List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b)
+ !$OMP          List_comb_thr_b3_cent, int2_u2_j1b2_test,ao_abs_comb_b3_j1b,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -260,11 +258,11 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b3_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+        if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_j_mu_x_2(i_fit)
           coef_fit = coef_gauss_j_mu_x_2(i_fit)
-          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.1.d-10)cycle
+!          if(dabs(coef_fit*int_j1b*sq_pi_3_2*(expo_fit)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
           int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
           tmp += coef_fit * int_fit
         enddo
@@ -278,7 +276,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
           coef        = List_comb_thr_b3_coef  (i_1s,j,i)
           beta        = List_comb_thr_b3_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!         if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -288,8 +286,7 @@ BEGIN_PROVIDER [ double precision, int2_u2_j1b2_test, (ao_num, ao_num, n_points_
             coef_fit = coef_gauss_j_mu_x_2(i_fit)
             !DIR$ FORCEINLINE
             call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
-!            if(dabs(coef_fit*coef*factor_ij_1s*int_j1b).lt.1.d-10)cycle ! old version
+!            if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.thrsh_cycle_tc)cycle
-            if(dabs(coef_fit*coef*factor_ij_1s*int_j1b*sq_pi_3_2*(beta_ij)**(-1.5d0)).lt.1.d-10)cycle
             int_fit = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
             tmp += coef * coef_fit * int_fit
           enddo
@@ -350,7 +347,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
  !$OMP          final_grid_points, ng_fit_jast,                     &
  !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
  !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,       & 
- !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2)
+ !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_x_j1b2_test,ao_abs_comb_b3_j1b,sq_pi_3_2,thrsh_cycle_tc)
  !$OMP DO
 
   do ipoint = 1, n_points_final_grid
@@ -369,7 +366,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
           coef        = List_comb_thr_b3_coef  (i_1s,j,i)
           beta        = List_comb_thr_b3_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -392,8 +389,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_x_j1b2_test, (ao_num, ao_num, n
             expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist 
             coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
             sq_alpha = alpha_1s_inv * dsqrt(alpha_1s_inv)
-!            if(dabs(coef_tmp*int_j1b) .lt. 1d-10) cycle ! old version
+!            if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. thrsh_cycle_tc) cycle
-            if(dabs(coef_tmp*int_j1b*sq_pi_3_2*sq_alpha) .lt. 1d-10) cycle
             
             call NAI_pol_x_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s, 1.d+9, r, int_fit)
 
@@ -470,13 +466,13 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
  !$OMP          expo_gauss_j_mu_1_erf, coef_gauss_j_mu_1_erf,       &
  !$OMP          ao_prod_dist_grid, ao_prod_sigma, ao_overlap_abs_grid,ao_prod_center,dsqpi_3_2,   &
  !$OMP          List_comb_thr_b3_coef, List_comb_thr_b3_expo,  ao_abs_comb_b3_j1b,     &
- !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test)
+ !$OMP          List_comb_thr_b3_cent, int2_u_grad1u_j1b2_test,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     do i = 1, ao_num
       do j = i, ao_num
 
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10) cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc) cycle
 
         r(1) = final_grid_points(1,ipoint)
         r(2) = final_grid_points(2,ipoint)
@@ -489,10 +485,10 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b3_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+!        if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_j_mu_1_erf(i_fit)
-          if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.1.d-15) cycle
+!          if(dabs(int_j1b)*dsqpi_3_2*expo_fit**(-1.5d0).lt.thrsh_cycle_tc) cycle
           coef_fit = coef_gauss_j_mu_1_erf(i_fit)
           int_fit  = NAI_pol_mult_erf_ao_with1s(i, j, expo_fit, r,  1.d+9, r)
           tmp += coef_fit * int_fit
@@ -507,7 +503,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
           coef        = List_comb_thr_b3_coef  (i_1s,j,i)
           beta        = List_comb_thr_b3_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b3_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b3_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b3_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b3_cent(3,i_1s,j,i)
@@ -517,7 +513,7 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
           do i_fit = 1, ng_fit_jast
             expo_fit = expo_gauss_j_mu_1_erf(i_fit)
             call gaussian_product(expo_fit,r,beta,B_center,factor_ij_1s,beta_ij,center_ij_1s)
-            if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.1.d-15)cycle
+!            if(factor_ij_1s*dabs(coef*int_j1b)*dsqpi_3_2*beta_ij**(-1.5d0).lt.thrsh_cycle_tc)cycle
             coef_fit = coef_gauss_j_mu_1_erf(i_fit)
 
             alpha_1s     = beta + expo_fit
@@ -527,9 +523,9 @@ BEGIN_PROVIDER [ double precision, int2_u_grad1u_j1b2_test, (ao_num, ao_num, n_p
             centr_1s(3)  = alpha_1s_inv * (beta * B_center(3) + expo_fit * r(3))
 
             expo_coef_1s = beta * expo_fit * alpha_1s_inv * dist
-            if(expo_coef_1s .gt. 20.d0) cycle
+!            if(expo_coef_1s .gt. 20.d0) cycle
             coef_tmp = coef * coef_fit * dexp(-expo_coef_1s)
-            if(dabs(coef_tmp) .lt. 1d-08) cycle
+!            if(dabs(coef_tmp) .lt. 1d-08) cycle
 
             int_fit = NAI_pol_mult_erf_ao_with1s(i, j, alpha_1s, centr_1s,  1.d+9, r)
 

src/ao_many_one_e_ints/grad_lapl_jmu_manu.irp.f

@@ -31,7 +31,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
  !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points, &
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,ao_abs_comb_b2_j1b,  &
  !$OMP          v_ij_erf_rk_cst_mu_j1b_test, mu_erf,                                   &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
  !$OMP DO
   !do ipoint = 1, 10
   do ipoint = 1, n_points_final_grid
@@ -41,7 +41,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp = 0.d0
         do i_1s = 1, List_comb_thr_b2_size(j,i)
@@ -49,7 +49,7 @@ BEGIN_PROVIDER [ double precision, v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_num,
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!         if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@@ -110,7 +110,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
  !$OMP SHARED  (n_points_final_grid, ao_num, List_comb_thr_b2_size, final_grid_points,&
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo, List_comb_thr_b2_cent,  &
  !$OMP          x_v_ij_erf_rk_cst_mu_j1b_test, mu_erf,ao_abs_comb_b2_j1b,         &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,thrsh_cycle_tc)
 ! !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,expo_erfc_mu_gauss)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
@@ -120,7 +120,7 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-10)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp_x = 0.d0
         tmp_y = 0.d0
@@ -130,19 +130,11 @@ BEGIN_PROVIDER [ double precision, x_v_ij_erf_rk_cst_mu_j1b_test, (ao_num, ao_nu
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+ !        if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
 
-!          if(ao_prod_center(1,j,i).ne.10000.d0)then
-!           ! approximate 1 - erf(mu r12) by a gaussian * 10
-!           !DIR$ FORCEINLINE
-!           call gaussian_product(expo_erfc_mu_gauss,r,     &
-!                ao_prod_sigma(j,i),ao_prod_center(1,j,i),  & 
-!                factor_ij_1s,beta_ij,center_ij_1s)
-!           if(dabs(coef * factor_ij_1s*int_j1b*10.d0 * dsqpi_3_2 * beta_ij**(-1.5d0)).lt.1.d-10)cycle 
-!          endif
           call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center, mu_erf, r, ints        )
           call NAI_pol_x_mult_erf_ao_with1s(i, j, beta, B_center,  1.d+9, r, ints_coulomb)
 
@@ -216,7 +208,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
  !$OMP          expo_gauss_j_mu_x, coef_gauss_j_mu_x,               &
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size, & 
  !$OMP          List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_test,ao_abs_comb_b2_j1b,   &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -225,7 +217,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp = 0.d0
 
@@ -234,11 +226,11 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b2_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+ !      if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         do i_fit = 1, ng_fit_jast
           expo_fit = expo_gauss_j_mu_x(i_fit)
           coef_fit = coef_gauss_j_mu_x(i_fit)
-          if(ao_overlap_abs_grid(j,i).lt.1.d-15) cycle
+ !        if(ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc) cycle
           int_fit = overlap_gauss_r12_ao(r, expo_fit, i, j)
           tmp += coef_fit * int_fit
         enddo
@@ -251,7 +243,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@@ -259,9 +251,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_test, (ao_num, ao_num, n_po
             expo_fit = expo_gauss_j_mu_x(i_fit)
             coef_fit = coef_gauss_j_mu_x(i_fit)
             coeftot = coef * coef_fit
-            if(dabs(coeftot).lt.1.d-15)cycle
+!            if(dabs(coeftot).lt.thrsh_cycle_tc)cycle
             call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
-            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
+!            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle
             int_fit  = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
             tmp += coef * coef_fit * int_fit
           enddo
@@ -325,7 +317,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
  !$OMP          expo_gauss_j_mu_x, coef_gauss_j_mu_x,                                  &
  !$OMP          List_comb_thr_b2_coef, List_comb_thr_b2_expo,List_comb_thr_b2_size,    & 
  !$OMP          List_comb_thr_b2_cent, v_ij_u_cst_mu_j1b_ng_1_test,ao_abs_comb_b2_j1b, &
- !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2)
+ !$OMP          ao_overlap_abs_grid,ao_prod_center,ao_prod_sigma,dsqpi_3_2,thrsh_cycle_tc)
  !$OMP DO
   do ipoint = 1, n_points_final_grid
     r(1) = final_grid_points(1,ipoint)
@@ -334,7 +326,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
 
     do i = 1, ao_num
       do j = i, ao_num
-        if(dabs(ao_overlap_abs_grid(j,i)).lt.1.d-20)cycle
+        if(dabs(ao_overlap_abs_grid(j,i)).lt.thrsh_cycle_tc)cycle
 
         tmp = 0.d0
 
@@ -343,7 +335,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
         ! --- --- ---
 
         int_j1b = ao_abs_comb_b2_j1b(1,j,i)
-        if(dabs(int_j1b).lt.1.d-10) cycle
+!        if(dabs(int_j1b).lt.thrsh_cycle_tc) cycle
         expo_fit = expo_good_j_mu_1gauss
         int_fit  = overlap_gauss_r12_ao(r, expo_fit, i, j)
         tmp += int_fit
@@ -356,7 +348,7 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
           coef        = List_comb_thr_b2_coef  (i_1s,j,i)
           beta        = List_comb_thr_b2_expo  (i_1s,j,i)
           int_j1b = ao_abs_comb_b2_j1b(i_1s,j,i)
-          if(dabs(coef)*dabs(int_j1b).lt.1.d-10)cycle
+!          if(dabs(coef)*dabs(int_j1b).lt.thrsh_cycle_tc)cycle
           B_center(1) = List_comb_thr_b2_cent(1,i_1s,j,i)
           B_center(2) = List_comb_thr_b2_cent(2,i_1s,j,i)
           B_center(3) = List_comb_thr_b2_cent(3,i_1s,j,i)
@@ -364,9 +356,9 @@ BEGIN_PROVIDER [ double precision, v_ij_u_cst_mu_j1b_ng_1_test, (ao_num, ao_num,
             expo_fit = expo_good_j_mu_1gauss
             coef_fit = 1.d0
             coeftot = coef * coef_fit
-            if(dabs(coeftot).lt.1.d-15)cycle
+            if(dabs(coeftot).lt.thrsh_cycle_tc)cycle
             call gaussian_product(beta,B_center,expo_fit,r,factor_ij_1s_u,beta_ij_u,center_ij_1s_u)
-            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.1.d-15)cycle
+            if(factor_ij_1s_u*ao_overlap_abs_grid(j,i).lt.thrsh_cycle_tc)cycle
             int_fit  = overlap_gauss_r12_ao_with1s(B_center, beta, r, expo_fit, i, j)
             tmp += coef * coef_fit * int_fit
 !          enddo

src/ao_many_one_e_ints/listj1b_sorted.irp.f

@@ -3,15 +3,16 @@
 &BEGIN_PROVIDER [ integer, max_List_comb_thr_b2_size]
  implicit none
  integer :: i_1s,i,j,ipoint
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  List_comb_thr_b2_size = 0
+ print*,'List_all_comb_b2_size = ',List_all_comb_b2_size
+! pause
  do i = 1, ao_num
   do j = i, ao_num
    do i_1s = 1, List_all_comb_b2_size
      coef        = List_all_comb_b2_coef  (i_1s)
-     if(dabs(coef).lt.1.d-15)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      beta        = List_all_comb_b2_expo  (i_1s)
      beta = max(beta,1.d-12)
      center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
@@ -24,7 +25,7 @@
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       List_comb_thr_b2_size(j,i) += 1
      endif
    enddo
@@ -40,6 +41,7 @@
   list(i) = maxval(List_comb_thr_b2_size(:,i))
  enddo
  max_List_comb_thr_b2_size = maxval(list) 
+ print*,'max_List_comb_thr_b2_size = ',max_List_comb_thr_b2_size
  
 END_PROVIDER 
 
@@ -49,16 +51,15 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, ao_abs_comb_b2_j1b, ( max_List_comb_thr_b2_size ,ao_num, ao_num)]
  implicit none
  integer :: i_1s,i,j,ipoint,icount
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  ao_abs_comb_b2_j1b = 10000000.d0
  do i = 1, ao_num
   do j = i, ao_num
    icount = 0
    do i_1s = 1, List_all_comb_b2_size
      coef        = List_all_comb_b2_coef  (i_1s)
-     if(dabs(coef).lt.1.d-12)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      beta        = List_all_comb_b2_expo  (i_1s)
      center(1:3) = List_all_comb_b2_cent(1:3,i_1s)
      int_j1b = 0.d0
@@ -70,7 +71,7 @@ END_PROVIDER
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       icount += 1
       List_comb_thr_b2_coef(icount,j,i) = coef
       List_comb_thr_b2_expo(icount,j,i) = beta
@@ -98,17 +99,17 @@ END_PROVIDER
 &BEGIN_PROVIDER [ integer, max_List_comb_thr_b3_size]
  implicit none
  integer :: i_1s,i,j,ipoint
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  List_comb_thr_b3_size = 0
+ print*,'List_all_comb_b3_size = ',List_all_comb_b3_size
  do i = 1, ao_num
   do j = 1, ao_num
    do i_1s = 1, List_all_comb_b3_size
      coef        = List_all_comb_b3_coef  (i_1s)
      beta        = List_all_comb_b3_expo  (i_1s)
      center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
-     if(dabs(coef).lt.thr)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      int_j1b = 0.d0
      do ipoint = 1, n_points_extra_final_grid
       r(1:3) = final_grid_points_extra(1:3,ipoint)
@@ -118,7 +119,7 @@ END_PROVIDER
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       List_comb_thr_b3_size(j,i) += 1
      endif
    enddo
@@ -144,9 +145,8 @@ END_PROVIDER
 &BEGIN_PROVIDER [ double precision, ao_abs_comb_b3_j1b, ( max_List_comb_thr_b3_size ,ao_num, ao_num)]
  implicit none
  integer :: i_1s,i,j,ipoint,icount
- double precision :: coef,beta,center(3),int_j1b,thr
+ double precision :: coef,beta,center(3),int_j1b
  double precision :: r(3),weight,dist
- thr = 1.d-15
  ao_abs_comb_b3_j1b = 10000000.d0
  do i = 1, ao_num
   do j = 1, ao_num
@@ -156,7 +156,7 @@ END_PROVIDER
      beta        = List_all_comb_b3_expo  (i_1s)
      beta = max(beta,1.d-12)
      center(1:3) = List_all_comb_b3_cent(1:3,i_1s)
-     if(dabs(coef).lt.thr)cycle
+     if(dabs(coef).lt.thrsh_cycle_tc)cycle
      int_j1b = 0.d0
      do ipoint = 1, n_points_extra_final_grid
       r(1:3) = final_grid_points_extra(1:3,ipoint)
@@ -166,7 +166,7 @@ END_PROVIDER
       dist += ( center(3) - r(3) )*( center(3) - r(3) )
       int_j1b += dabs(aos_in_r_array_extra_transp(ipoint,i) * aos_in_r_array_extra_transp(ipoint,j))*dexp(-beta*dist) * weight
      enddo
-     if(dabs(coef)*dabs(int_j1b).gt.thr)then
+     if(dabs(coef)*dabs(int_j1b).gt.thrsh_cycle_tc)then
       icount += 1
       List_comb_thr_b3_coef(icount,j,i) = coef
       List_comb_thr_b3_expo(icount,j,i) = beta
@@ -177,15 +177,5 @@ END_PROVIDER
   enddo 
  enddo
 
-! do i = 1, ao_num
-!  do j = 1, i-1
-!    do icount = 1, List_comb_thr_b3_size(j,i)
-!     List_comb_thr_b3_coef(icount,j,i) = List_comb_thr_b3_coef(icount,i,j)
-!     List_comb_thr_b3_expo(icount,j,i) = List_comb_thr_b3_expo(icount,i,j)
-!     List_comb_thr_b3_cent(1:3,icount,j,i) = List_comb_thr_b3_cent(1:3,icount,i,j)
-!    enddo
-!  enddo
-! enddo
- 
 END_PROVIDER 
 

src/bi_ort_ints/three_body_ijmk.irp.f

@@ -27,7 +27,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_direct_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_direct_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -74,7 +74,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_1_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_cycle_1_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -121,7 +121,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_cycle_2_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_cycle_2_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -168,7 +168,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch23_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_exch23_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -214,7 +214,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch13_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_exch13_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -261,7 +261,7 @@ BEGIN_PROVIDER [ double precision, three_e_4_idx_exch12_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)             &
  !$OMP PRIVATE (i,j,k,m,integral) & 
  !$OMP SHARED (mo_num,three_e_4_idx_exch12_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num

src/bi_ort_ints/three_body_ijmkl.irp.f

@@ -26,7 +26,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_direct_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_direct_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -75,7 +75,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_1_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_cycle_1_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -124,7 +124,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_cycle_2_bi_ort, (mo_num, mo_num
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_cycle_2_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -173,7 +173,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch23_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_exch23_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -222,7 +222,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch13_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_exch13_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num
@@ -271,7 +271,7 @@ BEGIN_PROVIDER [ double precision, three_e_5_idx_exch12_bi_ort, (mo_num, mo_num,
  !$OMP DEFAULT (NONE)               &
  !$OMP PRIVATE (i,j,k,m,l,integral) & 
  !$OMP SHARED (mo_num,three_e_5_idx_exch12_bi_ort)
- !$OMP DO SCHEDULE (dynamic)
+ !$OMP DO SCHEDULE (dynamic) COLLAPSE(2)
   do i = 1, mo_num
     do k = 1, mo_num
       do j = 1, mo_num

src/bi_ort_ints/total_twoe_pot.irp.f

@@ -57,6 +57,7 @@ BEGIN_PROVIDER [double precision, ao_two_e_tc_tot, (ao_num, ao_num, ao_num, ao_n
 
     PROVIDE ao_tc_sym_two_e_pot_in_map
 
+    !!! TODO :: OPENMP
     do j = 1, ao_num
       do l = 1, ao_num
         do i = 1, ao_num

src/davidson_keywords/input.irp.f

@@ -33,11 +33,6 @@ BEGIN_PROVIDER [ integer, n_states_diag  ]
     endif
   IRP_ENDIF
 
-  call write_time(6)
-  if (mpi_master) then
-    write(6, *) 'Read  n_states_diag'
-  endif
-
 END_PROVIDER
 
 ! ---

src/mo_basis/mo_class.irp.f

@@ -1,8 +1,3 @@
-! DO NOT MODIFY BY HAND
-! Created by $QP_ROOT/scripts/ezfio_interface/ei_handler.py
-! from file /home/eginer/programs/qp2/src/mo_basis/EZFIO.cfg
-
-
 BEGIN_PROVIDER [ character*(32), mo_class , (mo_num) ]
   implicit none
   BEGIN_DOC
@@ -35,6 +30,4 @@ BEGIN_PROVIDER [ character*(32), mo_class , (mo_num) ]
     endif
   IRP_ENDIF
 
-  call write_time(6)
-
 END_PROVIDER

src/non_h_ints_mu/grad_squared.irp.f

@@ -2,7 +2,7 @@
 ! ---
 
 !       TODO : strong optmization : write the loops in a different way
-!            : for each couple of AO, the gaussian product are done once for all 
+!            : for each couple of AO, the gaussian product are done once for all
 
 BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_points_final_grid) ]
 
@@ -20,14 +20,14 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
   ! gradu_squared_u_ij_mu = -0.50 x \int r2 \phi_i(2) \phi_j(2) [ v1^2  v2^2 ((grad_1 u12)^2 + (grad_2 u12^2)]) + u12^2  v2^2 (grad_1 v1)^2 + 2 u12 v1 v2^2 (grad_1 u12) . (grad_1 v1) ]
   !                       = -0.25 x        v1^2      \int r2 \phi_i(2) \phi_j(2) [1 - erf(mu r12)]^2 v2^2
   !                       + -0.50 x    (grad_1 v1)^2 \int r2 \phi_i(2) \phi_j(2)      u12^2          v2^2
-  !                       + -1.00 x v1 (grad_1 v1)   \int r2 \phi_i(2) \phi_j(2)   (grad_1 u12)      v2^2 
+  !                       + -1.00 x v1 (grad_1 v1)   \int r2 \phi_i(2) \phi_j(2)   (grad_1 u12)      v2^2
   !                       =                 v1^2        x   int2_grad1u2_grad2u2_j1b2
   !                       + -0.5 x      (grad_1 v1)^2   x   int2_u2_j1b2
   !                       + -1.0 X V1 x (grad_1 v1)   \cdot [ int2_u_grad1u_j1b2 x r - int2_u_grad1u_x_j1b ]
   !
   !
   END_DOC
- 
+
   implicit none
   integer                    :: ipoint, i, j, m, igauss
   double precision           :: x, y, z, r(3), delta, coef
@@ -100,7 +100,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
   call wall_time(time1)
   print*, ' Wall time for gradu_squared_u_ij_mu = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
@@ -151,7 +151,7 @@ END_PROVIDER
 !
 !  deallocate(ac_mat)
 !
-!END_PROVIDER 
+!END_PROVIDER
 
 ! ---
 
@@ -214,12 +214,12 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_loop, (ao_num, ao_num, ao_nu
   call wall_time(time1)
   print*, ' Wall time for tc_grad_square_ao_loop = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
 BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_grid) ]
- 
+
   implicit none
   integer                    :: ipoint, i, j, m, igauss
   double precision           :: r(3), delta, coef
@@ -267,7 +267,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12, (ao_num, ao_num, n_points_final_g
   call wall_time(time1)
   print*, ' Wall time for grad12_j12 = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
@@ -297,12 +297,12 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_
   call wall_time(time1)
   print*, ' Wall time for u12sq_j1bsq = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
 BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num, n_points_final_grid) ]
- 
+
   implicit none
   integer                    :: ipoint, i, j, m, igauss
   double precision           :: x, y, z
@@ -347,7 +347,7 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b, (ao_num, ao_num,
   call wall_time(time1)
   print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
@@ -370,26 +370,18 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   if(read_tc_integ) then
 
-    open(unit=11, form="unformatted", file='tc_grad_square_ao', action="read")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao', action="read")
-      do i = 1, ao_num
+    read(11) tc_grad_square_ao
-        do j = 1, ao_num
-          do k = 1, ao_num
-            do l = 1, ao_num
-              read(11) tc_grad_square_ao(l,k,j,i)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
 
   else
 
     allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
-  
+
     b_mat = 0.d0
    !$OMP PARALLEL               &
    !$OMP DEFAULT (NONE)         &
-   !$OMP PRIVATE (i, k, ipoint) & 
+   !$OMP PRIVATE (i, k, ipoint) &
    !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
    !$OMP DO SCHEDULE (static)
     do i = 1, ao_num
@@ -401,11 +393,11 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
     enddo
    !$OMP END DO
    !$OMP END PARALLEL
-  
+
     tmp = 0.d0
    !$OMP PARALLEL               &
    !$OMP DEFAULT (NONE)         &
-   !$OMP PRIVATE (j, l, ipoint) & 
+   !$OMP PRIVATE (j, l, ipoint) &
    !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
    !$OMP DO SCHEDULE (static)
     do ipoint = 1, n_points_final_grid
@@ -417,25 +409,25 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
     enddo
    !$OMP END DO
    !$OMP END PARALLEL
-  
+
     tc_grad_square_ao = 0.d0
     call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
               , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid      &
               , 1.d0, tc_grad_square_ao, ao_num*ao_num)
     deallocate(tmp, b_mat)
-  
+
     call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)
-  
+
    !!$OMP PARALLEL             &
    !!$OMP DEFAULT (NONE)       &
-   !!$OMP PRIVATE (i, j, k, l) & 
+   !!$OMP PRIVATE (i, j, k, l) &
    !!$OMP SHARED (ac_mat, tc_grad_square_ao, ao_num)
    !!$OMP DO SCHEDULE (static)
    ! do j = 1, ao_num
    !   do l = 1, ao_num
    !     do i = 1, ao_num
    !       do k = 1, ao_num
-   !         tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i) 
+   !         tc_grad_square_ao(k,i,l,j) = ac_mat(k,i,l,j) + ac_mat(l,j,k,i)
    !       enddo
    !     enddo
    !   enddo
@@ -444,23 +436,17 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
    !!$OMP END PARALLEL
   endif
 
-  if(write_tc_integ) then
+  if(write_tc_integ.and.mpi_master) then
-    open(unit=11, form="unformatted", file='tc_grad_square_ao', action="write")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao', action="write")
-      do i = 1, ao_num
+    call ezfio_set_work_empty(.False.)
-        do j = 1, ao_num
+    write(11) tc_grad_square_ao
-          do k = 1, ao_num
-            do l = 1, ao_num
-              write(11) tc_grad_square_ao(l,k,j,i)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
+    call ezfio_set_tc_keywords_io_tc_integ('Read')
   endif
 
   call wall_time(time1)
   print*, ' Wall time for tc_grad_square_ao = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---

src/non_h_ints_mu/grad_squared_manu.irp.f

@@ -17,29 +17,21 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
   call wall_time(time0)
 
   if(read_tc_integ) then
-  
+
-    open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="read")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao_test', action="read")
-      do i = 1, ao_num
+    read(11) tc_grad_square_ao_test
-        do j = 1, ao_num
-          do k = 1, ao_num
-            do l = 1, ao_num
-              read(11) tc_grad_square_ao_test(l,k,j,i)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
 
   else
 
     provide u12sq_j1bsq_test u12_grad1_u12_j1b_grad1_j1b_test grad12_j12_test
-  
+
     allocate(b_mat(n_points_final_grid,ao_num,ao_num), tmp(ao_num,ao_num,n_points_final_grid))
-  
+
     b_mat = 0.d0
    !$OMP PARALLEL               &
    !$OMP DEFAULT (NONE)         &
-   !$OMP PRIVATE (i, k, ipoint) & 
+   !$OMP PRIVATE (i, k, ipoint) &
    !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
    !$OMP DO SCHEDULE (static)
     do i = 1, ao_num
@@ -51,11 +43,11 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
     enddo
    !$OMP END DO
    !$OMP END PARALLEL
-  
+
     tmp = 0.d0
    !$OMP PARALLEL               &
    !$OMP DEFAULT (NONE)         &
-   !$OMP PRIVATE (j, l, ipoint) & 
+   !$OMP PRIVATE (j, l, ipoint) &
    !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test)
    !$OMP DO SCHEDULE (static)
     do ipoint = 1, n_points_final_grid
@@ -67,23 +59,23 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
     enddo
    !$OMP END DO
    !$OMP END PARALLEL
-  
+
     tc_grad_square_ao_test = 0.d0
     call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0 &
               , tmp(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid      &
               , 1.d0, tc_grad_square_ao_test, ao_num*ao_num)
     deallocate(tmp, b_mat)
-  
+
     call sum_A_At(tc_grad_square_ao_test(1,1,1,1), ao_num*ao_num)
     !do i = 1, ao_num
     !  do j = 1, ao_num
     !    do k = i, ao_num
-  
+
     !      do l = max(j,k), ao_num
     !        tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
     !        tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
     !      end do
-  
+
     !      !if (j.eq.k) then
     !      !  do l = j+1, ao_num
     !      !    tc_grad_square_ao_test(i,j,k,l) = 0.5d0 * (tc_grad_square_ao_test(i,j,k,l) + tc_grad_square_ao_test(k,l,i,j))
@@ -95,14 +87,14 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
     !      !    tc_grad_square_ao_test(k,l,i,j) = tc_grad_square_ao_test(i,j,k,l)
     !      !  enddo
     !      !endif
-  
+
     !    enddo
     !  enddo
     !enddo
     !tc_grad_square_ao_test = 2.d0 * tc_grad_square_ao_test
   ! !$OMP PARALLEL             &
   ! !$OMP DEFAULT (NONE)       &
-  ! !$OMP PRIVATE (i, j, k, l) & 
+  ! !$OMP PRIVATE (i, j, k, l) &
   ! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
   ! !$OMP DO SCHEDULE (static)
   !  integer :: ii
@@ -121,10 +113,10 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
   !  print *, ' ii =', ii
   ! !$OMP END DO
   ! !$OMP END PARALLEL
-  
+
   ! !$OMP PARALLEL             &
   ! !$OMP DEFAULT (NONE)       &
-  ! !$OMP PRIVATE (i, j, k, l) & 
+  ! !$OMP PRIVATE (i, j, k, l) &
   ! !$OMP SHARED (tc_grad_square_ao_test, ao_num)
   ! !$OMP DO SCHEDULE (static)
   !  do j = 1, ao_num
@@ -144,24 +136,18 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test, (ao_num, ao_num, ao_nu
 
   endif
 
-  if(write_tc_integ) then
+  if(write_tc_integ.and.mpi_master) then
-    open(unit=11, form="unformatted", file='tc_grad_square_ao_test', action="write")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_square_ao_test', action="write")
-      do i = 1, ao_num
+    call ezfio_set_work_empty(.False.)
-        do j = 1, ao_num
+    write(11) tc_grad_square_ao_test
-          do k = 1, ao_num
-            do l = 1, ao_num
-              write(11) tc_grad_square_ao_test(l,k,j,i)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
+    call ezfio_set_tc_keywords_io_tc_integ('Read')
   endif
 
   call wall_time(time1)
   print*, ' Wall time for tc_grad_square_ao_test = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
@@ -189,7 +175,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
   b_mat = 0.d0
  !$OMP PARALLEL               &
  !$OMP DEFAULT (NONE)         &
- !$OMP PRIVATE (i, k, ipoint) & 
+ !$OMP PRIVATE (i, k, ipoint) &
  !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector)
  !$OMP DO SCHEDULE (static)
   do i = 1, ao_num
@@ -205,7 +191,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
   tmp = 0.d0
  !$OMP PARALLEL               &
  !$OMP DEFAULT (NONE)         &
- !$OMP PRIVATE (j, l, ipoint) & 
+ !$OMP PRIVATE (j, l, ipoint) &
  !$OMP SHARED (tmp, ao_num, n_points_final_grid, u12sq_j1bsq_test, u12_grad1_u12_j1b_grad1_j1b_test, grad12_j12_test)
  !$OMP DO SCHEDULE (static)
   do ipoint = 1, n_points_final_grid
@@ -226,7 +212,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
 
  !$OMP PARALLEL             &
  !$OMP DEFAULT (NONE)       &
- !$OMP PRIVATE (i, j, k, l) & 
+ !$OMP PRIVATE (i, j, k, l) &
  !$OMP SHARED (ac_mat, tc_grad_square_ao_test_ref, ao_num)
  !$OMP DO SCHEDULE (static)
   do j = 1, ao_num
@@ -246,7 +232,7 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao_test_ref, (ao_num, ao_num, a
   call wall_time(time1)
   print*, ' Wall time for tc_grad_square_ao_test_ref = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
@@ -276,12 +262,12 @@ BEGIN_PROVIDER [ double precision, u12sq_j1bsq_test, (ao_num, ao_num, n_points_f
   call wall_time(time1)
   print*, ' Wall time for u12sq_j1bsq_test = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
 BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b_test, (ao_num, ao_num, n_points_final_grid) ]
- 
+
   implicit none
   integer                    :: ipoint, i, j, m, igauss
   double precision           :: x, y, z
@@ -328,12 +314,12 @@ BEGIN_PROVIDER [ double precision, u12_grad1_u12_j1b_grad1_j1b_test, (ao_num, ao
   call wall_time(time1)
   print*, ' Wall time for u12_grad1_u12_j1b_grad1_j1b_test = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
 BEGIN_PROVIDER [ double precision, grad12_j12_test, (ao_num, ao_num, n_points_final_grid) ]
- 
+
   implicit none
   integer                    :: ipoint, i, j, m, igauss
   double precision           :: r(3), delta, coef
@@ -381,7 +367,7 @@ BEGIN_PROVIDER [ double precision, grad12_j12_test, (ao_num, ao_num, n_points_fi
   call wall_time(time1)
   print*, ' Wall time for grad12_j12_test = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 

src/non_h_ints_mu/new_grad_tc.irp.f

@@ -36,16 +36,8 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_
 
   if(read_tc_integ) then
 
-    open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="read")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read")
-      do m = 1, 3
+    read(11) int2_grad1_u12_ao
-        do ipoint = 1, n_points_final_grid
-          do j = 1, ao_num
-            do i = 1, ao_num
-              read(11) int2_grad1_u12_ao(i,j,ipoint,m) 
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
 
   else
@@ -89,18 +81,12 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_
 
   endif
 
-  if(write_tc_integ) then
+  if(write_tc_integ.and.mpi_master) then
-    open(unit=11, form="unformatted", file='int2_grad1_u12_ao', action="write")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write")
-      do m = 1, 3
+    call ezfio_set_work_empty(.False.)
-        do ipoint = 1, n_points_final_grid
+    write(11) int2_grad1_u12_ao
-          do j = 1, ao_num
-            do i = 1, ao_num
-              write(11) int2_grad1_u12_ao(i,j,ipoint,m) 
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
+    call ezfio_set_tc_keywords_io_tc_integ('Read')
   endif
 
   call wall_time(time1)

src/non_h_ints_mu/new_grad_tc_manu.irp.f

@@ -3,7 +3,7 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
 
   BEGIN_DOC
   !
-  ! int2_grad1_u12_ao_test(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2) 
+  ! int2_grad1_u12_ao_test(i,j,ipoint,:) = \int dr2 [-1 * \grad_r1 J(r1,r2)] \phi_i(r2) \phi_j(r2)
   !
   ! where r1 = r(ipoint)
   !
@@ -15,9 +15,9 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
   ! if J(r1,r2) = u12 x v1 x v2
   !
   ! int2_grad1_u12_ao_test(i,j,ipoint,:) =      v1    x [ 0.5 x \int dr2 [(r1 - r2) (erf(mu * r12)-1)r_12] v2 \phi_i(r2) \phi_j(r2) ]
-  !                                 - \grad_1 v1 x [       \int dr2                  u12              v2 \phi_i(r2) \phi_j(r2) ] 
+  !                                 - \grad_1 v1 x [       \int dr2                  u12              v2 \phi_i(r2) \phi_j(r2) ]
-  !                                 =    0.5 v_1b(ipoint) * v_ij_erf_rk_cst_mu_j1b(i,j,ipoint) * r(:) 
+  !                                 =    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,:) 
+  !                                 -    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)
   !
   !
@@ -35,25 +35,18 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
 
   if(read_tc_integ) then
 
-    open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="read")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao_test', action="read")
-      do m = 1, 3
+    read(11) int2_grad1_u12_ao_test
-        do ipoint = 1, n_points_final_grid
-          do j = 1, ao_num
-            do i = 1, ao_num
-              read(11) int2_grad1_u12_ao_test(i,j,ipoint,m)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
 
+
   else
-  
+
     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) 
+        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)
@@ -87,24 +80,18 @@ BEGIN_PROVIDER [ double precision, int2_grad1_u12_ao_test, (ao_num, ao_num, n_po
 
   endif
 
-  if(write_tc_integ) then
+  if(write_tc_integ.and.mpi_master) then
-    open(unit=11, form="unformatted", file='int2_grad1_u12_ao_test', action="write")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao_test', action="write")
-      do m = 1, 3
+    call ezfio_set_work_empty(.False.)
-        do ipoint = 1, n_points_final_grid
+    write(11) int2_grad1_u12_ao_test
-          do j = 1, ao_num
-            do i = 1, ao_num
-              write(11) int2_grad1_u12_ao_test(i,j,ipoint,m)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
+    call ezfio_set_tc_keywords_io_tc_integ('Read')
   endif
 
   call wall_time(time1)
   print*, ' Wall time for int2_grad1_u12_ao_test = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 
@@ -114,9 +101,9 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
   !
   ! tc_grad_and_lapl_ao_test(k,i,l,j) = < k l | -1/2 \Delta_1 u(r1,r2) - \grad_1 u(r1,r2) | ij >
   !
-  ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2) 
+  ! = 1/2 \int dr1 (phi_k(r1) \grad_r1 phi_i(r1) - phi_i(r1) \grad_r1 phi_k(r1)) . \int dr2 \grad_r1 u(r1,r2) \phi_l(r2) \phi_j(r2)
   !
-  ! This is obtained by integration by parts. 
+  ! This is obtained by integration by parts.
   !
   END_DOC
 
@@ -131,40 +118,32 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
   call wall_time(time0)
 
   if(read_tc_integ) then
-  
+
-    open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="read")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao_test', action="read")
-      do i = 1, ao_num
+    read(11) tc_grad_and_lapl_ao_test
-        do j = 1, ao_num
-          do k = 1, ao_num
-            do l = 1, ao_num
-              read(11) tc_grad_and_lapl_ao_test(l,k,j,i)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
 
   else
 
-    provide int2_grad1_u12_ao_test 
+    provide int2_grad1_u12_ao_test
-   
+
     allocate(b_mat(n_points_final_grid,ao_num,ao_num,3), ac_mat(ao_num,ao_num,ao_num,ao_num))
-  
+
     b_mat = 0.d0
    !$OMP PARALLEL                                                              &
    !$OMP DEFAULT (NONE)                                                        &
-   !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r)                       & 
+   !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r)                       &
-   !$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, & 
+   !$OMP SHARED (aos_in_r_array_transp, aos_grad_in_r_array_transp_bis, b_mat, &
    !$OMP         ao_num, n_points_final_grid, final_weight_at_r_vector)
    !$OMP DO SCHEDULE (static)
     do i = 1, ao_num
       do k = 1, ao_num
         do ipoint = 1, n_points_final_grid
-  
+
           weight1 = 0.5d0 * final_weight_at_r_vector(ipoint)
           ao_i_r  = aos_in_r_array_transp(ipoint,i)
           ao_k_r  = aos_in_r_array_transp(ipoint,k)
-  
+
           b_mat(ipoint,k,i,1) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,1) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,1))
           b_mat(ipoint,k,i,2) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,2) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,2))
           b_mat(ipoint,k,i,3) = weight1 * (ao_k_r * aos_grad_in_r_array_transp_bis(ipoint,i,3) - ao_i_r * aos_grad_in_r_array_transp_bis(ipoint,k,3))
@@ -173,19 +152,19 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
     enddo
    !$OMP END DO
    !$OMP END PARALLEL
-  
+
     ac_mat = 0.d0
     do m = 1, 3
       call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0                   &
                 , int2_grad1_u12_ao_test(1,1,1,m), ao_num*ao_num, b_mat(1,1,1,m), n_points_final_grid &
                 , 1.d0, ac_mat, ao_num*ao_num)
-  
+
     enddo
     deallocate(b_mat)
-  
+
    !$OMP PARALLEL             &
    !$OMP DEFAULT (NONE)       &
-   !$OMP PRIVATE (i, j, k, l) & 
+   !$OMP PRIVATE (i, j, k, l) &
    !$OMP SHARED (ac_mat, tc_grad_and_lapl_ao_test, ao_num)
    !$OMP DO SCHEDULE (static)
     do j = 1, ao_num
@@ -199,29 +178,23 @@ BEGIN_PROVIDER [double precision, tc_grad_and_lapl_ao_test, (ao_num, ao_num, ao_
     enddo
    !$OMP END DO
    !$OMP END PARALLEL
-  
+
     deallocate(ac_mat)
 
   endif
 
-  if(write_tc_integ) then
+  if(write_tc_integ.and.mpi_master) then
-    open(unit=11, form="unformatted", file='tc_grad_and_lapl_ao_test', action="write")
+    open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/tc_grad_and_lapl_ao_test', action="write")
-      do i = 1, ao_num
+    call ezfio_set_work_empty(.False.)
-        do j = 1, ao_num
+    write(11) tc_grad_and_lapl_ao_test
-          do k = 1, ao_num
-            do l = 1, ao_num
-              write(11) tc_grad_and_lapl_ao_test(l,k,j,i)
-            enddo
-          enddo
-        enddo
-      enddo
     close(11)
+    call ezfio_set_tc_keywords_io_tc_integ('Read')
   endif
 
   call wall_time(time1)
   print*, ' Wall time for tc_grad_and_lapl_ao_test = ', time1 - time0
 
-END_PROVIDER 
+END_PROVIDER
 
 ! ---
 

src/non_h_ints_mu/total_tc_int.irp.f

@@ -68,7 +68,26 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
 
 END_PROVIDER 
 
+BEGIN_PROVIDER [double precision, ao_tc_int_chemist_no_cycle, (ao_num, ao_num, ao_num, ao_num)]
 ! ---
+  implicit none
+  integer          :: i, j, k, l
+  double precision :: wall1, wall0
+  print *, ' providing ao_tc_int_chemist_no_cycle ...'
+  call wall_time(wall0)
+    do j = 1, ao_num
+      do l = 1, ao_num
+        do i = 1, ao_num
+          do k = 1, ao_num
+            ao_tc_int_chemist_no_cycle(k,i,l,j) = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
+!            ao_tc_int_chemist(k,i,l,j) = ao_two_e_coul(k,i,l,j)
+          enddo
+        enddo
+      enddo
+    enddo
+  call wall_time(wall1)
+  print *, ' wall time for ao_tc_int_chemist_no_cycle ', wall1 - wall0
+END_PROVIDER 
 
 BEGIN_PROVIDER [double precision, ao_tc_int_chemist_test, (ao_num, ao_num, ao_num, ao_num)]
 

src/tc_keywords/EZFIO.cfg

@@ -6,7 +6,7 @@ default: False
 
 [comp_left_eigv]
 type: logical
-doc: If |true|, computes also the left-eigenvector 
+doc: If |true|, computes also the left-eigenvector
 interface: ezfio,provider,ocaml
 default: False
 
@@ -14,7 +14,7 @@ default: False
 type: logical
 doc: If |true|, three-body terms are included
 interface: ezfio,provider,ocaml
-default: True 
+default: True
 
 [pure_three_body_h_tc]
 type: logical
@@ -30,13 +30,13 @@ default: False
 
 [core_tc_op]
 type: logical
-doc: If |true|, takes the usual Hamiltonian for core orbitals (assumed to be doubly occupied) 
+doc: If |true|, takes the usual Hamiltonian for core orbitals (assumed to be doubly occupied)
 interface: ezfio,provider,ocaml
 default: False
 
 [full_tc_h_solver]
 type: logical
-doc: If |true|, you diagonalize the full TC H matrix 
+doc: If |true|, you diagonalize the full TC H matrix
 interface: ezfio,provider,ocaml
 default: False
 
@@ -46,6 +46,12 @@ doc: Thresholds on the energy for iterative Davidson used in TC
 interface: ezfio,provider,ocaml
 default: 1.e-5
 
+[thrsh_cycle_tc]
+type: Threshold
+doc: Thresholds to cycle the integrals with the envelop
+interface: ezfio,provider,ocaml
+default: 1.e-10
+
 [max_it_dav]
 type: integer
 doc: nb max of iteration in Davidson used in TC
@@ -60,11 +66,11 @@ default: 0.000005
 
 [thresh_psi_r_norm]
 type: logical
-doc: If |true|, you prune the WF to compute the PT1 coef based on the norm. If False, the pruning is done through the amplitude on the right-coefficient. 
+doc: If |true|, you prune the WF to compute the PT1 coef based on the norm. If False, the pruning is done through the amplitude on the right-coefficient.
 interface: ezfio,provider,ocaml
 default: False
 
-[state_following_tc] 
+[state_following_tc]
 type: logical
 doc: If |true|, the states are re-ordered to match the input states
 default: False
@@ -78,7 +84,7 @@ default: True
 
 [symetric_fock_tc]
 type: logical
-doc: If |true|, using F+F^t as Fock TC 
+doc: If |true|, using F+F^t as Fock TC
 interface: ezfio,provider,ocaml
 default: False
 
@@ -126,7 +132,7 @@ default: 1.e-6
 
 [maxovl_tc]
 type: logical
-doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors 
+doc: If |true|, maximize the overlap between orthogonalized left- and right eigenvectors
 interface: ezfio,provider,ocaml
 default: False
 
@@ -152,7 +158,7 @@ default: 0.
 type: character*(32)
 doc: Type of TCSCF algorithm used. Possible choices are [Simple | DIIS]
 interface: ezfio,provider,ocaml
-default: Simple 
+default: DIIS
 
 [im_thresh_tcscf]
 type: Threshold
@@ -180,21 +186,15 @@ default: 1.e-6
 
 [var_tc]
 type: logical
-doc: If |true|, use VAR-TC 
+doc: If |true|, use VAR-TC
 interface: ezfio,provider,ocaml
 default: False
 
-[read_tc_integ]
+[io_tc_integ]
-type: logical
+type: Disk_access
-doc: If |true|, read integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao
+doc: Read/Write integrals int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao from/to disk [ Write | Read | None ]
 interface: ezfio,provider,ocaml
-default: False
+default: None
-
-[write_tc_integ]
-type: logical
-doc: If |true|, write integrals: int2_grad1_u12_ao, tc_grad_square_ao and tc_grad_and_lapl_ao
-interface: ezfio,provider,ocaml
-default: False
 
 [debug_tc_pt2]
 type: integer

src/tc_scf/tc_scf.irp.f

@@ -11,8 +11,8 @@ program tc_scf
   print *, ' starting ...'
 
   my_grid_becke  = .True.
-  my_n_pt_r_grid = 30
+  my_n_pt_r_grid = 60
-  my_n_pt_a_grid = 50
+  my_n_pt_a_grid = 110
 !  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/test_int.irp.f

@@ -21,25 +21,22 @@ program test_ints
   touch my_extra_grid_becke my_n_pt_r_extra_grid my_n_pt_a_extra_grid
 
 !! OK 
-!call routine_int2_u_grad1u_j1b2 
+! call routine_int2_u_grad1u_j1b2 
-!! OK
+! OK
-!call routine_v_ij_erf_rk_cst_mu_j1b
+! call routine_v_ij_erf_rk_cst_mu_j1b
-!! OK 
+! OK 
 ! call routine_x_v_ij_erf_rk_cst_mu_j1b
-!! OK
+! OK
-! call routine_v_ij_u_cst_mu_j1b
+! call routine_int2_u2_j1b2
-
+! OK
-!! OK
+! call routine_int2_u_grad1u_x_j1b2
-!call routine_int2_u2_j1b2
+! OK 
-
-!! OK
-!call routine_int2_u_grad1u_x_j1b2
-
-!! OK 
 ! call routine_int2_grad1u2_grad2u2_j1b2
 ! call routine_int2_u_grad1u_j1b2
 ! call test_total_grad_lapl
 ! call test_total_grad_square
+! call test_int2_grad1_u12_ao_test
+! call routine_v_ij_u_cst_mu_j1b_test
 ! call test_ao_tc_int_chemist
 ! call test_grid_points_ao
 ! call test_tc_scf
@@ -53,12 +50,12 @@ program test_ints
 
   !call test_two_e_tc_non_hermit_integral()
 
-  call test_tc_grad_square_ao_test()
+!  call test_tc_grad_square_ao_test()
-
-  PROVIDE TC_HF_energy VARTC_HF_energy
-  print *, '    TC_HF_energy = ',    TC_HF_energy
-  print *, ' VARTC_HF_energy = ', VARTC_HF_energy
 
+!!PROVIDE TC_HF_energy VARTC_HF_energy
+!!print *, '    TC_HF_energy = ',    TC_HF_energy
+!!print *, ' VARTC_HF_energy = ', VARTC_HF_energy
+ call test_old_ints
 end
 
 ! ---
@@ -157,6 +154,9 @@ subroutine routine_int2_u_grad1u_j1b2
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,'routine_int2_u_grad1u_j1b2'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -169,20 +169,6 @@ subroutine routine_v_ij_erf_rk_cst_mu_j1b
  integer :: i,j,ipoint,k,l
  double precision :: weight,accu_relat, accu_abs, contrib
  double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-! print*,'ao_overlap_abs = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
-! enddo
-! print*,'center = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
-! enddo
-! print*,'sigma = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
-! enddo
-
-
  allocate(array(ao_num, ao_num, ao_num, ao_num))
  array = 0.d0
  allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
@@ -215,6 +201,9 @@ subroutine routine_v_ij_erf_rk_cst_mu_j1b
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,'routine_v_ij_erf_rk_cst_mu_j1b'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -228,20 +217,6 @@ subroutine routine_x_v_ij_erf_rk_cst_mu_j1b
  integer :: i,j,ipoint,k,l,m
  double precision :: weight,accu_relat, accu_abs, contrib
  double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-! print*,'ao_overlap_abs = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
-! enddo
-! print*,'center = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
-! enddo
-! print*,'sigma = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
-! enddo
-
-
  allocate(array(ao_num, ao_num, ao_num, ao_num))
  array = 0.d0
  allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
@@ -276,6 +251,10 @@ subroutine routine_x_v_ij_erf_rk_cst_mu_j1b
     enddo
    enddo
   enddo
+
+ print*,'******'
+ print*,'******'
+ print*,'routine_x_v_ij_erf_rk_cst_mu_j1b'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -290,20 +269,6 @@ subroutine routine_v_ij_u_cst_mu_j1b_test
  integer :: i,j,ipoint,k,l
  double precision :: weight,accu_relat, accu_abs, contrib
  double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-! print*,'ao_overlap_abs = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
-! enddo
-! print*,'center = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
-! enddo
-! print*,'sigma = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
-! enddo
-
-
  allocate(array(ao_num, ao_num, ao_num, ao_num))
  array = 0.d0
  allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
@@ -336,6 +301,9 @@ subroutine routine_v_ij_u_cst_mu_j1b_test
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,'routine_v_ij_u_cst_mu_j1b_test'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -427,19 +395,6 @@ subroutine routine_int2_u2_j1b2
  integer :: i,j,ipoint,k,l
  double precision :: weight,accu_relat, accu_abs, contrib
  double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-! print*,'ao_overlap_abs = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
-! enddo
-! print*,'center = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
-! enddo
-! print*,'sigma = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
-! enddo
-
 
  allocate(array(ao_num, ao_num, ao_num, ao_num))
  array = 0.d0
@@ -473,6 +428,9 @@ subroutine routine_int2_u2_j1b2
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,'routine_int2_u2_j1b2'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -486,19 +444,6 @@ subroutine routine_int2_u_grad1u_x_j1b2
  integer :: i,j,ipoint,k,l,m
  double precision :: weight,accu_relat, accu_abs, contrib
  double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
-! print*,'ao_overlap_abs = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_overlap_abs(i,:)
-! enddo
-! print*,'center = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_center(2,i,:)
-! enddo
-! print*,'sigma = '
-! do i = 1, ao_num
-!   write(*,'(100(F10.5,X))')ao_prod_sigma(i,:)
-! enddo
-
 
  allocate(array(ao_num, ao_num, ao_num, ao_num))
  array = 0.d0
@@ -534,6 +479,9 @@ subroutine routine_int2_u_grad1u_x_j1b2
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,'routine_int2_u_grad1u_x_j1b2'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -579,6 +527,9 @@ subroutine routine_v_ij_u_cst_mu_j1b
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,'routine_v_ij_u_cst_mu_j1b'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -742,6 +693,9 @@ subroutine test_total_grad_lapl
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,' test_total_grad_lapl'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -767,6 +721,9 @@ subroutine test_total_grad_square
     enddo
    enddo
   enddo
+ print*,'******'
+ print*,'******'
+ print*,'test_total_grad_square'
  print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
  print*,'accu_relat = ',accu_relat/dble(ao_num)**4
 
@@ -1057,3 +1014,86 @@ end
 ! ---
 
 
+
+subroutine test_old_ints
+ implicit none
+ integer :: i,j,k,l
+ double precision :: old, new, contrib, get_ao_tc_sym_two_e_pot
+ double precision :: integral_sym , integral_nsym,accu
+    PROVIDE ao_tc_sym_two_e_pot_in_map
+ accu = 0.d0
+ do j = 1, ao_num
+  do l= 1, ao_num
+   do i = 1, ao_num
+    do k = 1, ao_num
+!      integral_sym  = get_ao_tc_sym_two_e_pot(i, j, k, l, ao_tc_sym_two_e_pot_map)
+      ! ao_non_hermit_term_chemist(k,i,l,j) = < k l | [erf( mu r12) - 1] d/d_r12 | i j > on the AO basis
+!      integral_nsym = ao_non_hermit_term_chemist(k,i,l,j)
+!      old = integral_sym + integral_nsym 
+!      old = tc_grad_square_ao(k,i,l,j) + tc_grad_and_lapl_ao(k,i,l,j) + ao_two_e_coul(k,i,l,j)
+      new = ao_tc_int_chemist_test(k,i,l,j)
+      old = ao_tc_int_chemist_no_cycle(k,i,l,j)
+      contrib = dabs(old - new)
+      if(contrib.gt.1.d-6)then
+       print*,'problem !!'
+       print*,i,j,k,l
+       print*,old, new, contrib
+      endif
+       accu += contrib
+    enddo
+   enddo
+  enddo
+ enddo
+ print*,'******'
+ print*,'******'
+ print*,'in test_old_ints'
+ print*,'accu = ',accu/dble(ao_num**4)
+
+end
+
+subroutine test_int2_grad1_u12_ao_test
+ implicit none
+ integer :: i,j,ipoint,m,k,l
+ double precision :: weight,accu_relat, accu_abs, contrib
+ double precision, allocatable :: array(:,:,:,:), array_ref(:,:,:,:)
+ allocate(array(ao_num, ao_num, ao_num, ao_num))
+ array = 0.d0
+ allocate(array_ref(ao_num, ao_num, ao_num, ao_num))
+ array_ref = 0.d0
+ do m = 1, 3
+  do ipoint = 1, n_points_final_grid
+   weight = final_weight_at_r_vector(ipoint)
+   do k = 1, ao_num
+    do l = 1, ao_num
+     do i = 1, ao_num
+      do j = 1, ao_num
+       array(j,i,l,k)     += int2_grad1_u12_ao_test(j,i,ipoint,m) * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
+       array_ref(j,i,l,k) += int2_grad1_u12_ao(j,i,ipoint,m)      * aos_grad_in_r_array_transp(m,k,ipoint) * aos_in_r_array(l,ipoint) * weight
+      enddo
+     enddo
+    enddo
+   enddo
+  enddo
+ enddo
+
+ accu_relat = 0.d0
+ accu_abs   = 0.d0
+  do k = 1, ao_num
+   do l = 1, ao_num
+    do i = 1, ao_num
+     do j = 1, ao_num
+      contrib = dabs(array(j,i,l,k) - array_ref(j,i,l,k))
+      accu_abs += contrib
+      if(dabs(array_ref(j,i,l,k)).gt.1.d-10)then
+       accu_relat += contrib/dabs(array_ref(j,i,l,k))
+      endif
+     enddo
+    enddo
+   enddo
+  enddo
+ print*,'******'
+ print*,'******'
+ print*,'test_int2_grad1_u12_ao_test'
+ print*,'accu_abs   = ',accu_abs/dble(ao_num)**4
+ print*,'accu_relat = ',accu_relat/dble(ao_num)**4
+end

src/tools/print_sorted_wf_coef.irp.f

@@ -0,0 +1,19 @@
+program sort_wf
+ implicit none
+ read_wf =.True.
+ call routine
+
+end
+
+subroutine routine
+ implicit none
+ integer :: i
+ character*(128) :: output
+ integer :: i_unit_output,getUnitAndOpen
+ output=trim(ezfio_filename)//'.wf_sorted'
+ i_unit_output = getUnitAndOpen(output,'w')
+ do i= 1, N_det
+  write(i_unit_output,*)i,dabs(psi_coef_sorted(i,1))
+ enddo
+
+end

src/tools/save_natorb.irp.f

@@ -12,6 +12,7 @@ program save_natorb
 ! matrices of each state with the corresponding
 ! :option:`determinants state_average_weight`
   END_DOC
+  PROVIDE nucl_coord
   read_wf = .True.
   touch read_wf
   call save_natural_mos

src/tools/sort_wf.irp.f

@@ -0,0 +1,20 @@
+program sort_wf
+ implicit none
+ read_wf = .true.
+ touch read_wf 
+ call routine
+
+end
+
+subroutine routine
+ implicit none
+ integer :: i
+ character*(128) :: output
+ integer :: i_unit_output,getUnitAndOpen
+ output=trim(ezfio_filename)//'.wf_sorted'
+ i_unit_output = getUnitAndOpen(output,'w')
+ do i = 1, N_det
+  write(i_unit_output, *)i,dabs(psi_coef_sorted(i,1))/dabs(psi_coef_sorted(1,1))
+ enddo
+
+end