Merge pull request #282 from AbdAmmar/dev-stable-tc-scf

Dev stable tc scf

config/ifort_2021_xHost.cfg

@@ -6,7 +6,7 @@
 # --align=32                 : Align all provided arrays on a 32-byte boundary
 #
 [COMMON]
-FC           : ifort -fpic
+FC           : ifort -fpic -diag-disable 5462
 LAPACK_LIB   : -mkl=parallel
 IRPF90       : irpf90
 IRPF90_FLAGS : --ninja --align=64 -DINTEL 

configure

@@ -232,7 +232,7 @@ EOF
 
             execute << EOF
               cd "\${QP_ROOT}"/external
-              tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.2.tar.gz
+              tar --gunzip --extract --file qp2-dependencies/f77-zmq-4.3.?.tar.gz
               cd f77-zmq-*
               ./configure --prefix=\$QP_ROOT
               export ZMQ_H="\$QP_ROOT"/include/zmq.h

scripts/compilation/qp_create_ninja

@@ -25,7 +25,7 @@ except ImportError:
                                       "quantum_package.rc"))
 
     print("\n".join(["", "Error:", "source %s" % f, ""]))
-    sys.exit(1)
+    raise
 
 # Compress path
 def comp_path(path):

scripts/qp_exc_energy.py

@@ -0,0 +1,186 @@
+#!/usr/bin/env python
+# Computes the error on the excitation energy of a CIPSI run.
+
+def student(p,df):
+    import scipy
+    from scipy.stats import t
+    return t.ppf(p, df)
+
+
+def chi2cdf(x,k):
+    import scipy
+    import scipy.stats
+    return scipy.stats.chi2.cdf(x,k)
+
+
+def jarque_bera(data):
+
+    n    = max(len(data), 2)
+    norm = 1./ sum( [ w for (_,w) in data ] )
+
+    mu     =   sum( [ w* x        for (x,w) in data ] ) * norm
+    sigma2 =   sum( [ w*(x-mu)**2 for (x,w) in data ] ) * norm
+    if sigma2 > 0.:
+      S      = ( sum( [ w*(x-mu)**3 for (x,w) in data ] ) * norm ) / sigma2**(3./2.)
+      K      = ( sum( [ w*(x-mu)**4 for (x,w) in data ] ) * norm ) / sigma2**2
+    else:
+      S = 0.
+      K = 0.
+
+    # Value of the Jarque-Bera test
+    JB = n/6. * (S**2 + 1./4. * (K-3.)**2)
+
+    # Probability that the data comes from a Gaussian distribution
+    p = 1. - chi2cdf(JB,2)
+
+    return JB, mu, sqrt(sigma2/(n-1)), p
+
+
+
+to_eV = 27.2107362681
+import sys, os
+import scipy
+import scipy.stats
+from math import sqrt, gamma, exp
+import json
+
+
+def read_data(filename,state):
+  """ Read energies and PT2 from input file """
+  with open(filename,'r') as f:
+      lines = json.load(f)['fci']
+
+  print(f"State: {state}")
+
+  gs = []
+  es = []
+
+  for l in lines:
+      try:
+        pt2_0 = l['states'][0]['pt2']
+        e_0   = l['states'][0]['energy']
+        pt2_1 = l['states'][state]['pt2']
+        e_1   = l['states'][state]['energy']
+        gs.append( (e_0, pt2_0) )
+        es.append( (e_1, pt2_1) )
+      except: pass
+
+  def f(p_1, p0, p1):
+      e, pt2 = p0
+      y0, x0 = p_1
+      y1, x1 = p1
+      try:
+        alpha = (y1-y0)/(x0-x1)
+      except ZeroDivisionError:
+        alpha = 1.
+      return [e, pt2, alpha]
+
+  for l in (gs, es):
+      p_1, p0, p1 = l[0], l[0], l[1]
+      l[0] = f(p_1, p0, p1)
+
+      for i in range(1,len(l)-1):
+          p_1 = (l[i-1][0], l[i-1][1])
+          p0  = l[i]
+          p1  = l[i+1]
+          l[i] = f(p_1, p0, p1)
+
+      i = len(l)-1
+      p_1 = (l[i-1][0], l[i-1][1])
+      p0  = l[i]
+      p1  = l[-1]
+      l[i] = f(p_1, p0, p1)
+
+  return [ x+y for x,y in zip(gs,es) ]
+
+
+def compute(data):
+
+    d = []
+    for e0, p0, a0, e1, p1, a1 in data:
+        x   = (e1+p1)-(e0+p0)
+        w   = 1./sqrt(p0**2 + p1**2)
+        bias = (a1-1.)*p1 - (a0-1.)*p0
+        d.append( (x,w,bias) )
+
+    x = []
+    target = (scipy.stats.norm.cdf(1.)-0.5)*2   # = 0.6827
+
+    print("| %2s | %8s | %8s | %8s | %8s | %8s |"%( "N", "DE", "+/-", "bias", "P(G)", "J"))
+    print("|----+----------+----------+----------+----------+----------|")
+    xmax = (0.,0.,0.,0.,0.,0,0.)
+    for i in range(len(data)-1):
+        jb, mu, sigma, p = jarque_bera( [ (x,w) for (x,w,bias) in d[i:] ] )
+        bias = sum ( [ w * e for (_,w,e) in d[i:] ] ) / sum ( [ w for (_,w,_) in d[i:] ] )
+        mu = (mu+0.5*bias) * to_eV
+        sigma = sigma * to_eV
+        bias = bias * to_eV
+        n = len(data[i:])
+        beta = student(0.5*(1.+target/p) ,n)
+        err = sigma * beta + 0.5*abs(bias)
+        print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |"%( n, mu, err, bias, p, jb))
+        if n < 3 :
+            continue
+        y = (err, p, mu, err, jb,n,bias)
+        if p > xmax[1]: xmax = y
+        if p < 0.8:
+            continue
+        x.append(y)
+
+    x = sorted(x)
+
+    print("|----+----------+----------+----------+----------+----------|")
+    if x != []:
+      xmax = x[0]
+    _, p, mu, err, jb, n, bias = xmax
+    beta = student(0.5*(1.+target/p),n)
+    print("| %2d | %8.3f | %8.3f | %8.3f | %8.3f | %8.3f |\n"%(n, mu, err, bias, p, jb))
+
+    return mu, err, bias, p
+
+filename = sys.argv[1]
+print(filename)
+if len(sys.argv) > 2:
+    state = int(sys.argv[2])
+else:
+    state = 1
+data = read_data(filename,state)
+mu, err, bias, _ = compute(data)
+print(" %s: %8.3f +/- %5.3f eV\n"%(filename, mu, err))
+
+import numpy as np
+A = np.array( [ [ data[-1][1], 1. ],
+                [ data[-2][1], 1. ] ] )
+B = np.array( [ [ data[-1][0] ],
+                [ data[-2][0] ] ] )
+E0 = np.linalg.solve(A,B)[1]
+A = np.array( [ [ data[-1][4], 1. ],
+                [ data[-2][4], 1. ] ] )
+B = np.array( [ [ data[-1][3] ],
+                [ data[-2][3] ] ] )
+E1 = np.linalg.solve(A,B)[1]
+average_2 = (E1-E0)*to_eV
+
+A = np.array( [ [ data[-1][1], 1. ],
+                [ data[-2][1], 1. ],
+                [ data[-3][1], 1. ] ] )
+B = np.array( [ [ data[-1][0] ],
+                [ data[-2][0] ],
+                [ data[-3][0] ] ] )
+E0 = np.linalg.lstsq(A,B,rcond=None)[0][1]
+A = np.array( [ [ data[-1][4], 1. ],
+                [ data[-2][4], 1. ],
+                [ data[-3][4], 1. ] ] )
+B = np.array( [ [ data[-1][3] ],
+                [ data[-2][3] ],
+                [ data[-3][3] ] ] )
+E1 = np.linalg.lstsq(A,B,rcond=None)[0][1]
+average_3 = (E1-E0)*to_eV
+
+exc = ((data[-1][3] + data[-1][4]) - (data[-1][0] + data[-1][1])) * to_eV
+error_2 = abs(average_2 - average_3)
+error_3 = abs(average_3 - exc)
+print(" 2-3 points: %.3f +/- %.3f "% (average_3, error_2))
+print(" largest wf: %.3f +/- %.3f  "%(average_3, error_3))
+
+

src/ao_many_one_e_ints/NEED

@@ -3,3 +3,4 @@ ao_two_e_ints
 becke_numerical_grid
 mo_one_e_ints
 dft_utils_in_r
+tc_keywords

src/ao_many_one_e_ints/grad2_jmu_modif.irp.f

@@ -1,4 +1,72 @@
 
+
+! ---
+
+BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2, (ao_num, ao_num, n_points_final_grid)]
+
+  BEGIN_DOC
+  !
+  ! -\frac{1}{4} x int dr2 phi_i(r2) phi_j(r2) [1 - erf(mu r12)]^2
+  !
+  END_DOC
+
+  implicit none
+  integer                       :: i, j, ipoint, i_fit
+  double precision              :: r(3), expo_fit, coef_fit
+  double precision              :: tmp
+  double precision              :: wall0, wall1
+
+  double precision, external    :: overlap_gauss_r12_ao
+
+  print*, ' providing int2_grad1u2_grad2u2 ...'
+  call wall_time(wall0)
+
+  provide mu_erf final_grid_points j1b_pen
+
+  int2_grad1u2_grad2u2 = 0.d0
+
+  !$OMP PARALLEL DEFAULT (NONE)                                               &
+  !$OMP PRIVATE (ipoint, i, j, i_fit, r, coef_fit, expo_fit, tmp)             & 
+  !$OMP SHARED  (n_points_final_grid, ao_num, final_grid_points, ng_fit_jast, &
+  !$OMP          expo_gauss_1_erf_x_2, coef_gauss_1_erf_x_2,int2_grad1u2_grad2u2)
+  !$OMP DO
+  do ipoint = 1, n_points_final_grid
+    r(1) = final_grid_points(1,ipoint)
+    r(2) = final_grid_points(2,ipoint)
+    r(3) = final_grid_points(3,ipoint)
+
+    do i = 1, ao_num
+      do j = i, ao_num
+
+        tmp = 0.d0
+        do i_fit = 1, ng_fit_jast
+
+          expo_fit = expo_gauss_1_erf_x_2(i_fit)
+          coef_fit = coef_gauss_1_erf_x_2(i_fit)
+
+          tmp += -0.25d0 * coef_fit * overlap_gauss_r12_ao(r, expo_fit, i, j)
+        enddo
+
+        int2_grad1u2_grad2u2(j,i,ipoint) = tmp
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  do ipoint = 1, n_points_final_grid
+    do i = 2, ao_num
+      do j = 1, i-1
+        int2_grad1u2_grad2u2(j,i,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint)
+      enddo
+    enddo
+  enddo
+
+  call wall_time(wall1)
+  print*, ' wall time for int2_grad1u2_grad2u2 =', wall1 - wall0
+
+END_PROVIDER 
+
 ! ---
 
 BEGIN_PROVIDER [ double precision, int2_grad1u2_grad2u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
@@ -96,7 +164,7 @@ END_PROVIDER
 
 ! ---
 
-BEGIN_PROVIDER [ double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
+BEGIN_PROVIDER [double precision, int2_u2_j1b2, (ao_num, ao_num, n_points_final_grid)]
 
   BEGIN_DOC
   !

src/ao_many_one_e_ints/listj1b.irp.f

@@ -1,17 +1,34 @@
 
 ! ---
 
-BEGIN_PROVIDER [ integer, List_all_comb_b2_size]
+BEGIN_PROVIDER [integer, List_all_comb_b2_size]
 
   implicit none
 
+  PROVIDE j1b_type
+
+  if(j1b_type .eq. 3) then
+
     List_all_comb_b2_size = 2**nucl_num
 
+  elseif(j1b_type .eq. 4) then
+
+    List_all_comb_b2_size = nucl_num + 1
+
+  else 
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
+  print *, ' nb of linear terms in the envelope is ', List_all_comb_b2_size
+
 END_PROVIDER
 
 ! ---
 
-BEGIN_PROVIDER [ integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)]
+BEGIN_PROVIDER [integer, List_all_comb_b2, (nucl_num, List_all_comb_b2_size)]
 
   implicit none
   integer :: i, j
@@ -50,6 +67,8 @@ END_PROVIDER
   List_all_comb_b2_expo = 0.d0
   List_all_comb_b2_cent = 0.d0
 
+  if(j1b_type .eq. 3) then
+
     do i = 1, List_all_comb_b2_size
 
       tmp_cent_x = 0.d0
@@ -102,6 +121,26 @@ END_PROVIDER
       List_all_comb_b2_coef(i) = (-1.d0)**dble(phase) * dexp(-List_all_comb_b2_coef(i))
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    List_all_comb_b2_coef(    1) = 1.d0
+    List_all_comb_b2_expo(    1) = 0.d0
+    List_all_comb_b2_cent(1:3,1) = 0.d0
+    do i = 1, nucl_num
+      List_all_comb_b2_coef(  i+1) = -1.d0
+      List_all_comb_b2_expo(  i+1) = j1b_pen(   i)
+      List_all_comb_b2_cent(1,i+1) = nucl_coord(i,1)
+      List_all_comb_b2_cent(2,i+1) = nucl_coord(i,2)
+      List_all_comb_b2_cent(3,i+1) = nucl_coord(i,3)
+    enddo
+
+  else
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
   !print *, ' coeff, expo & cent of list b2'
   !do i = 1, List_all_comb_b2_size
   !  print*, i, List_all_comb_b2_coef(i), List_all_comb_b2_expo(i)
@@ -115,14 +154,31 @@ END_PROVIDER
 BEGIN_PROVIDER [ integer, List_all_comb_b3_size]
 
   implicit none
+  double precision :: tmp
+
+  if(j1b_type .eq. 3) then
 
     List_all_comb_b3_size = 3**nucl_num
 
+  elseif(j1b_type .eq. 4) then
+
+    tmp                   = 0.5d0 * dble(nucl_num) * (dble(nucl_num) + 3.d0)
+    List_all_comb_b3_size = int(tmp) + 1
+
+  else
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
+  print *, ' nb of linear terms in the square of the envelope is ', List_all_comb_b3_size
+
 END_PROVIDER
 
 ! ---
 
-BEGIN_PROVIDER [ integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)]
+BEGIN_PROVIDER [integer, List_all_comb_b3, (nucl_num, List_all_comb_b3_size)]
 
   implicit none
   integer              :: i, j, ii, jj
@@ -162,7 +218,11 @@ END_PROVIDER
 
   implicit none
   integer          :: i, j, k, phase
+  integer          :: ii
   double precision :: tmp_alphaj, tmp_alphak, facto
+  double precision :: tmp1, tmp2, tmp3, tmp4
+  double precision :: xi, yi, zi, xj, yj, zj
+  double precision :: dx, dy, dz, r2
 
   provide j1b_pen
 
@@ -170,6 +230,8 @@ END_PROVIDER
   List_all_comb_b3_expo = 0.d0
   List_all_comb_b3_cent = 0.d0
 
+  if(j1b_type .eq. 3) then
+
     do i = 1, List_all_comb_b3_size
 
       do j = 1, nucl_num
@@ -225,6 +287,70 @@ END_PROVIDER
       List_all_comb_b3_coef(i) = (-1.d0)**dble(phase) * facto * dexp(-List_all_comb_b3_coef(i))
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    ii = 1
+    List_all_comb_b3_coef(    ii) = 1.d0
+    List_all_comb_b3_expo(    ii) = 0.d0
+    List_all_comb_b3_cent(1:3,ii) = 0.d0
+
+    do i = 1, nucl_num
+      ii = ii + 1
+      List_all_comb_b3_coef(  ii) = -2.d0
+      List_all_comb_b3_expo(  ii) = j1b_pen(   i)
+      List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
+      List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
+      List_all_comb_b3_cent(3,ii) = nucl_coord(i,3)
+    enddo
+
+    do i = 1, nucl_num
+      ii = ii + 1
+      List_all_comb_b3_coef(  ii) = 1.d0
+      List_all_comb_b3_expo(  ii) = 2.d0 * j1b_pen(i)
+      List_all_comb_b3_cent(1,ii) = nucl_coord(i,1)
+      List_all_comb_b3_cent(2,ii) = nucl_coord(i,2)
+      List_all_comb_b3_cent(3,ii) = nucl_coord(i,3)
+    enddo
+
+    do i = 1, nucl_num-1
+
+      tmp1 = j1b_pen(i)
+
+      xi = nucl_coord(i,1)
+      yi = nucl_coord(i,2)
+      zi = nucl_coord(i,3)
+
+      do j = i+1, nucl_num
+
+        tmp2 = j1b_pen(j)
+        tmp3 = tmp1 + tmp2
+        tmp4 = 1.d0 / tmp3
+
+        xj = nucl_coord(j,1)
+        yj = nucl_coord(j,2)
+        zj = nucl_coord(j,3)
+
+        dx = xi - xj
+        dy = yi - yj
+        dz = zi - zj
+        r2 = dx*dx + dy*dy + dz*dz
+        
+        ii = ii + 1
+        List_all_comb_b3_coef(  ii) = dexp(-tmp1*tmp2*tmp4*r2)
+        List_all_comb_b3_expo(  ii) = tmp3
+        List_all_comb_b3_cent(1,ii) = tmp4 * (tmp1 * xi + tmp2 * xj)
+        List_all_comb_b3_cent(2,ii) = tmp4 * (tmp1 * yi + tmp2 * yj)
+        List_all_comb_b3_cent(3,ii) = tmp4 * (tmp1 * zi + tmp2 * zj)
+      enddo
+    enddo
+
+  else
+
+    print *, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
   !print *, ' coeff, expo & cent of list b3'
   !do i = 1, List_all_comb_b3_size
   !  print*, i, List_all_comb_b3_coef(i), List_all_comb_b3_expo(i)

src/ao_tc_eff_map/compute_ints_eff_pot.irp.f

@@ -53,13 +53,13 @@ subroutine compute_ao_tc_sym_two_e_pot_jl(j, l, n_integrals, buffer_i, buffer_va
       integral_erf = ao_two_e_integral_erf(i, k, j, l)
       integral     = integral_erf + integral_pot
 
-      if( j1b_type .eq. 1 ) then
+      !if( j1b_type .eq. 1 ) then
-        !print *, ' j1b type 1 is added'
+      !  !print *, ' j1b type 1 is added'
-        integral = integral + j1b_gauss_2e_j1(i, k, j, l)
+      !  integral = integral + j1b_gauss_2e_j1(i, k, j, l)
-      elseif( j1b_type .eq. 2 ) then
+      !elseif( j1b_type .eq. 2 ) then
-        !print *, ' j1b type 2 is added'
+      !  !print *, ' j1b type 2 is added'
-        integral = integral + j1b_gauss_2e_j2(i, k, j, l)
+      !  integral = integral + j1b_gauss_2e_j2(i, k, j, l)
-      endif
+      !endif
 
       if(abs(integral) < thr) then
         cycle

src/ao_two_e_ints/EZFIO.cfg

@@ -18,3 +18,8 @@ interface: ezfio,provider,ocaml
 default: False
 ezfio_name: direct
 
+[do_ao_cholesky]
+type: logical
+doc: Perform Cholesky decomposition of AO integrals
+interface: ezfio,provider,ocaml
+default: True

src/ao_two_e_ints/cholesky.irp.f

@@ -0,0 +1,100 @@
+BEGIN_PROVIDER [ integer, cholesky_ao_num_guess ]
+ implicit none
+ BEGIN_DOC
+ ! Number of Cholesky vectors in AO basis
+ END_DOC
+
+ integer :: i,j,k,l
+ double precision :: xnorm0, x, integral
+ double precision, external :: ao_two_e_integral
+
+ cholesky_ao_num_guess = 0
+ xnorm0 = 0.d0
+ x = 0.d0
+ do j=1,ao_num
+   do i=1,ao_num
+     integral = ao_two_e_integral(i,i,j,j)
+     if (integral > ao_integrals_threshold) then
+       cholesky_ao_num_guess += 1
+     else
+       x += integral
+     endif
+   enddo
+ enddo
+ print *, 'Cholesky decomposition of AO integrals'
+ print *, '--------------------------------------'
+ print *, ''
+ print *, 'Estimated Error: ', x
+ print *, 'Guess size: ', cholesky_ao_num_guess, '(', 100.d0*dble(cholesky_ao_num_guess)/dble(ao_num*ao_num), ' %)'
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, cholesky_ao_num ]
+&BEGIN_PROVIDER [ double precision, cholesky_ao, (ao_num, ao_num, cholesky_ao_num_guess) ]
+ use mmap_module
+ implicit none
+ BEGIN_DOC
+ ! Cholesky vectors in AO basis: (ik|a):
+ ! <ij|kl> = (ik|jl) = sum_a (ik|a).(a|jl)
+ END_DOC
+
+ type(c_ptr) :: ptr
+ integer :: fd, i,j,k,l, rank
+ double precision, pointer :: ao_integrals(:,:,:,:)
+ double precision, external :: ao_two_e_integral
+
+ ! Store AO integrals in a memory mapped file
+ call mmap(trim(ezfio_work_dir)//'ao_integrals', &
+   (/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), &
+   8, fd, .False., ptr)
+ call c_f_pointer(ptr, ao_integrals, (/ao_num, ao_num, ao_num, ao_num/))
+
+ double precision :: integral
+ logical, external :: ao_two_e_integral_zero
+ !$OMP PARALLEL DO DEFAULT(SHARED) PRIVATE(i,j,k,l, integral) SCHEDULE(dynamic)
+ do l=1,ao_num
+  do j=1,l
+   do k=1,ao_num
+    do i=1,k
+     if (ao_two_e_integral_zero(i,j,k,l)) cycle
+     integral = ao_two_e_integral(i,k,j,l)
+     ao_integrals(i,k,j,l) = integral
+     ao_integrals(k,i,j,l) = integral
+     ao_integrals(i,k,l,j) = integral
+     ao_integrals(k,i,l,j) = integral
+    enddo
+   enddo
+  enddo
+ enddo
+ !$OMP END PARALLEL DO
+
+ ! Call Lapack
+ cholesky_ao_num = cholesky_ao_num_guess
+ call pivoted_cholesky(ao_integrals, cholesky_ao_num, ao_integrals_threshold, ao_num*ao_num, cholesky_ao)
+ print *, 'Rank: ', cholesky_ao_num, '(', 100.d0*dble(cholesky_ao_num)/dble(ao_num*ao_num), ' %)'
+
+ ! Remove mmap
+ double precision, external :: getUnitAndOpen
+ call munmap( &
+   (/ int(ao_num,8), int(ao_num,8), int(ao_num,8), int(ao_num,8) /), &
+   8, fd, ptr)
+ open(unit=99,file=trim(ezfio_work_dir)//'ao_integrals')
+ close(99, status='delete')
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, cholesky_ao_transp, (cholesky_ao_num, ao_num, ao_num) ]
+ implicit none
+ BEGIN_DOC
+! Transposed of the Cholesky vectors in AO basis set
+ END_DOC
+ integer :: i,j,k
+ do j=1,ao_num
+  do i=1,ao_num
+   do k=1,ao_num
+    cholesky_ao_transp(k,i,j) = cholesky_ao(i,j,k)
+   enddo
+  enddo
+ enddo
+END_PROVIDER
+

src/ao_two_e_ints/map_integrals.irp.f

@@ -486,7 +486,7 @@ subroutine get_ao_two_e_integrals(j,k,l,sze,out_val)
   PROVIDE ao_two_e_integrals_in_map ao_integrals_map
 
   if (ao_one_e_integral_zero(j,l)) then
-    out_val = 0.d0
+    out_val(1:sze) = 0.d0
     return
   endif
 

src/bi_ort_ints/one_e_bi_ort.irp.f

@@ -8,22 +8,22 @@ BEGIN_PROVIDER [double precision, ao_one_e_integrals_tc_tot, (ao_num,ao_num)]
 
   ao_one_e_integrals_tc_tot = ao_one_e_integrals      
 
-  provide j1b_type
+  !provide j1b_type
 
-  if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
+  !if( (j1b_type .eq. 1) .or. (j1b_type .eq. 2) ) then
+  !
+  !  print *, ' do things properly !'
+  !  stop
 
-    print *, ' do things properly !'
+  !  !do i = 1, ao_num
-    stop
+  !  !  do j = 1, ao_num
+  !  !    ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI  (j,i) &
+  !  !                                      + j1b_gauss_hermII (j,i) &
+  !  !                                      + j1b_gauss_nonherm(j,i) )
+  !  !  enddo
+  !  !enddo
 
-    !do i = 1, ao_num
+  !endif
-    !  do j = 1, ao_num
-    !    ao_one_e_integrals_tc_tot(j,i) += ( j1b_gauss_hermI  (j,i) &
-    !                                      + j1b_gauss_hermII (j,i) &
-    !                                      + j1b_gauss_nonherm(j,i) )
-    !  enddo
-    !enddo
-
-  endif
 
 END_PROVIDER
 

src/cipsi/NEED

@@ -1,3 +1,4 @@
+json
 perturbation
 zmq
 mpi

src/cipsi/cipsi.irp.f

@@ -16,7 +16,6 @@ subroutine run_cipsi
   double precision, external :: memory_of_double
   PROVIDE H_apply_buffer_allocated
 
-  N_iter = 1
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators
 
@@ -76,7 +75,6 @@ subroutine run_cipsi
         )
       write(*,'(A)')  '--------------------------------------------------------------------------------'
 
-
     to_select = int(sqrt(dble(N_states))*dble(N_det)*selection_factor)
     to_select = max(N_states_diag, to_select)
     if (do_pt2) then
@@ -106,10 +104,10 @@ subroutine run_cipsi
 
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
 
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
     call print_mol_properties()
-    N_iter += 1
+    call write_cipsi_json(pt2_data,pt2_data_err)
 
     if (qp_stop()) exit
 
@@ -130,13 +128,13 @@ subroutine run_cipsi
     if (qp_stop()) exit
   enddo
 
-  if (.not.qp_stop()) then
+  ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
-    if (N_det < N_det_max) then
+  if ((.not.qp_stop()).and.                                          &
-        call diagonalize_CI
+        (N_det > N_det_max) .and.                                    &
-        call save_wavefunction
+        (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and.    &
-        call save_energy(psi_energy_with_nucl_rep, zeros)
+        (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
-    endif
+        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
-
+        ) then
     if (do_pt2) then
       call pt2_dealloc(pt2_data)
       call pt2_dealloc(pt2_data_err)
@@ -155,11 +153,13 @@ subroutine run_cipsi
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
     call print_summary(psi_energy_with_nucl_rep(1:N_states), &
       pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
     call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
   endif
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)
 
 end
+

src/cipsi/environment.irp.f

@@ -7,7 +7,9 @@ BEGIN_PROVIDER [ integer, nthreads_pt2 ]
  character*(32) :: env
  call getenv('QP_NTHREADS_PT2',env)
  if (trim(env) /= '') then
+   call lock_io()
    read(env,*) nthreads_pt2
+   call unlock_io()
    call write_int(6,nthreads_pt2,'Target number of threads for PT2')
  endif
 END_PROVIDER

src/cipsi/stochastic_cipsi.irp.f

@@ -15,7 +15,6 @@ subroutine run_stochastic_cipsi
   double precision, external :: memory_of_double
   PROVIDE H_apply_buffer_allocated distributed_davidson mo_two_e_integrals_in_map
 
-  N_iter = 1
   threshold_generators = 1.d0
   SOFT_TOUCH threshold_generators
 
@@ -96,10 +95,10 @@ subroutine run_stochastic_cipsi
 
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
 
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
     call print_mol_properties()
-    N_iter += 1
+    call write_cipsi_json(pt2_data,pt2_data_err)
 
     if (qp_stop()) exit
 
@@ -119,13 +118,13 @@ subroutine run_stochastic_cipsi
     if (qp_stop()) exit
   enddo
 
-  if (.not.qp_stop()) then
+  ! If stopped because N_det > N_det_max, do an extra iteration to compute the PT2
-    if (N_det < N_det_max) then
+  if ((.not.qp_stop()).and.                                          &
-        call diagonalize_CI
+        (N_det > N_det_max) .and.                                    &
-        call save_wavefunction
+        (maxval(abs(pt2_data % pt2(1:N_states))) > pt2_max) .and.    &
-        call save_energy(psi_energy_with_nucl_rep, zeros)
+        (maxval(abs(pt2_data % variance(1:N_states))) > variance_max) .and.&
-    endif
+        (correlation_energy_ratio <= correlation_energy_ratio_max)   &
-
+        ) then
     call pt2_dealloc(pt2_data)
     call pt2_dealloc(pt2_data_err)
     call pt2_alloc(pt2_data, N_states)
@@ -135,9 +134,10 @@ subroutine run_stochastic_cipsi
     call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
     call print_summary(psi_energy_with_nucl_rep, &
        pt2_data , pt2_data_err, N_det, N_configuration, N_states, psi_s2)
-    call save_iterations(psi_energy_with_nucl_rep(1:N_states),pt2_data % rpt2,N_det)
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
     call print_extrapolated_energy()
     call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
   endif
   call pt2_dealloc(pt2_data)
   call pt2_dealloc(pt2_data_err)

src/cipsi/write_cipsi_json.irp.f

@@ -0,0 +1,53 @@
+subroutine write_cipsi_json(pt2_data, pt2_data_err)
+    use selection_types
+    implicit none
+    BEGIN_DOC
+!  Writes JSON data for CIPSI runs
+    END_DOC
+    type(pt2_type), intent(in) :: pt2_data, pt2_data_err
+    integer :: i,j,k
+
+    call lock_io
+    character*(64), allocatable    :: fmtk(:)
+    integer :: N_states_p, N_iter_p
+    N_states_p = min(N_states,N_det)
+    N_iter_p = min(N_iter,8)
+    allocate(fmtk(0:N_iter_p))
+    fmtk(:) = '(''   '',E22.15,'','')'
+    fmtk(N_iter_p) = '(''   '',E22.15)'
+
+    write(json_unit, json_dict_uopen_fmt)
+    write(json_unit, json_int_fmt) 'n_det', N_det
+    if (s2_eig) then
+      write(json_unit, json_int_fmt) 'n_cfg', N_configuration
+      if (only_expected_s2) then
+        write(json_unit, json_int_fmt) 'n_csf', N_csf
+      endif
+    endif
+    write(json_unit, json_array_open_fmt) 'states'
+    do k=1,N_states_p
+      write(json_unit, json_dict_uopen_fmt)
+      write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k)
+      write(json_unit, json_real_fmt) 's2', psi_s2(k)
+      write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k)
+      write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k)
+      write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k)
+      write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k)
+      write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k)
+      write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k)
+      write(json_unit, json_array_open_fmt) 'ex_energy'
+      do i=2,N_iter_p
+          write(json_unit, fmtk(i)) extrapolated_energy(i,k)
+      enddo
+      write(json_unit, json_array_close_fmtx)
+      if (k < N_states_p) then
+        write(json_unit, json_dict_close_fmt)
+      else
+        write(json_unit, json_dict_close_fmtx)
+      endif
+    enddo
+    write(json_unit, json_array_close_fmtx)
+    write(json_unit, json_dict_close_fmt)
+    deallocate(fmtk)
+    call unlock_io
+end

src/cipsi_tc_bi_ortho/NEED

@@ -1,6 +1,7 @@
+json
 mpi
 perturbation
 zmq
-iterations_tc
+iterations
 csf
 tc_bi_ortho

src/cipsi_tc_bi_ortho/stochastic_cipsi.irp.f

@@ -94,7 +94,15 @@ subroutine run_stochastic_cipsi
     call ZMQ_pt2(E_denom, pt2_data, pt2_data_err, relative_error,to_select) ! Stochastic PT2 and selection
 !    stop
 
-    N_iter += 1
+    call print_summary(psi_energy_with_nucl_rep, &
+       pt2_data, pt2_data_err, N_det,N_configuration,N_states,psi_s2)
+
+    call save_energy(psi_energy_with_nucl_rep, pt2_data % pt2)
+
+    call increment_n_iter(psi_energy_with_nucl_rep, pt2_data)
+    call print_extrapolated_energy()
+!    call print_mol_properties()
+    call write_cipsi_json(pt2_data,pt2_data_err)
 
     if (qp_stop()) exit
 

src/cipsi_tc_bi_ortho/write_cipsi_json.irp.f

@@ -0,0 +1,53 @@
+subroutine write_cipsi_json(pt2_data, pt2_data_err)
+    use selection_types
+    implicit none
+    BEGIN_DOC
+!  Writes JSON data for CIPSI runs
+    END_DOC
+    type(pt2_type), intent(in) :: pt2_data, pt2_data_err
+    integer :: i,j,k
+
+    call lock_io
+    character*(64), allocatable    :: fmtk(:)
+    integer :: N_states_p, N_iter_p
+    N_states_p = min(N_states,N_det)
+    N_iter_p = min(N_iter,8)
+    allocate(fmtk(0:N_iter_p))
+    fmtk(:) = '(''   '',E22.15,'','')'
+    fmtk(N_iter_p) = '(''   '',E22.15)'
+
+    write(json_unit, json_dict_uopen_fmt)
+    write(json_unit, json_int_fmt) 'n_det', N_det
+    if (s2_eig) then
+      write(json_unit, json_int_fmt) 'n_cfg', N_configuration
+      if (only_expected_s2) then
+        write(json_unit, json_int_fmt) 'n_csf', N_csf
+      endif
+    endif
+    write(json_unit, json_array_open_fmt) 'states'
+    do k=1,N_states_p
+      write(json_unit, json_dict_uopen_fmt)
+      write(json_unit, json_real_fmt) 'energy', psi_energy_with_nucl_rep(k)
+      write(json_unit, json_real_fmt) 's2', psi_s2(k)
+      write(json_unit, json_real_fmt) 'pt2', pt2_data % pt2(k)
+      write(json_unit, json_real_fmt) 'pt2_err', pt2_data_err % pt2(k)
+      write(json_unit, json_real_fmt) 'rpt2', pt2_data % rpt2(k)
+      write(json_unit, json_real_fmt) 'rpt2_err', pt2_data_err % rpt2(k)
+      write(json_unit, json_real_fmt) 'variance', pt2_data % variance(k)
+      write(json_unit, json_real_fmt) 'variance_err', pt2_data_err % variance(k)
+      write(json_unit, json_array_open_fmt) 'ex_energy'
+      do i=2,N_iter_p
+          write(json_unit, fmtk(i)) extrapolated_energy(i,k)
+      enddo
+      write(json_unit, json_array_close_fmtx)
+      if (k < N_states_p) then
+        write(json_unit, json_dict_close_fmt)
+      else
+        write(json_unit, json_dict_close_fmtx)
+      endif
+    enddo
+    write(json_unit, json_array_close_fmtx)
+    write(json_unit, json_dict_close_fmt)
+    deallocate(fmtk)
+    call unlock_io
+end

src/davidson/davidson_parallel.irp.f

@@ -150,7 +150,9 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, N_st, sze,
       exit
     endif
     if(task_id == 0) exit
+    call lock_io()
     read (msg,*) imin, imax, ishift, istep
+    call unlock_io()
     integer :: k
     do k=imin,imax
       v_t(:,k) = 0.d0
@@ -546,6 +548,8 @@ end
 
 
 
+
+
 integer function zmq_put_N_states_diag(zmq_to_qp_run_socket,worker_id)
   use f77_zmq
   implicit none

src/davidson/diagonalization_hs2_dressed.irp.f

@@ -461,9 +461,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,s2_out,energies,dim_in,sze,N_
        integer :: lwork, info
        double precision, allocatable :: work(:)
 
-
        y = h
-       !y = h_p
+!       y = h_p   ! Doesn't work for non-singlets
        lwork = -1
        allocate(work(1))
        call dsygv(1,'V','U',shift2,y,size(y,1), &

src/ezfio_files/lock.irp.f

@@ -9,4 +9,21 @@ BEGIN_PROVIDER [ integer(omp_lock_kind), file_lock ]
   call omp_init_lock(file_lock)
 END_PROVIDER
 
+! These functions need to be called because internal read and write are not thread safe.
+subroutine lock_io()
+  implicit none
+  BEGIN_DOC
+! Needs to be called because before doing I/O because internal read and write
+! are not thread safe.
+  END_DOC
+  call omp_set_lock(file_lock)
+end subroutine lock_io()
 
+subroutine unlock_io()
+  implicit none
+  BEGIN_DOC
+! Needs to be called because afterdoing I/O because internal read and write
+! are not thread safe.
+  END_DOC
+  call omp_unset_lock(file_lock)
+end subroutine lock_io()

src/fci/fci.irp.f

@@ -39,12 +39,19 @@ program fci
   if (.not.is_zmq_slave) then
     PROVIDE psi_det psi_coef mo_two_e_integrals_in_map
 
+    write(json_unit,json_array_open_fmt) 'fci'
+
     if (do_pt2) then
       call run_stochastic_cipsi
     else
       call run_cipsi
     endif
 
+    write(json_unit,json_dict_uopen_fmt)
+    write(json_unit,json_dict_close_fmtx)
+    write(json_unit,json_array_close_fmtx)
+    call json_close
+
   else
     PROVIDE mo_two_e_integrals_in_map pt2_min_parallel_tasks
 

src/fci_tc_bi/NEED

@@ -1,3 +1,4 @@
+json
 tc_bi_ortho
 davidson_undressed
 cipsi_tc_bi_ortho

src/fci_tc_bi/fci_tc_bi_ortho.irp.f

@@ -62,6 +62,7 @@ subroutine run_cipsi_tc
       endif
     endif
     ! ---
+    write(json_unit,json_array_open_fmt) 'fci_tc'
 
     if (do_pt2) then
       call run_stochastic_cipsi
@@ -69,6 +70,11 @@ subroutine run_cipsi_tc
       call run_cipsi
     endif
 
+    write(json_unit,json_dict_uopen_fmt)
+    write(json_unit,json_dict_close_fmtx)
+    write(json_unit,json_array_close_fmtx)
+    call json_close
+
   else
     PROVIDE mo_bi_ortho_tc_one_e mo_bi_ortho_tc_two_e pt2_min_parallel_tasks
     if(elec_alpha_num+elec_beta_num.ge.3)then

src/fci_tc_bi/save_energy.irp.f

@@ -4,6 +4,6 @@ subroutine save_energy(E,pt2)
 ! Saves the energy in |EZFIO|.
   END_DOC
   double precision, intent(in) :: E(N_states), pt2(N_states)
-  call ezfio_set_fci_tc_energy(E(1:N_states))
+  call ezfio_set_fci_tc_bi_energy(E(1:N_states))
-  call ezfio_set_fci_tc_energy_pt2(E(1:N_states)+pt2(1:N_states))
+  call ezfio_set_fci_tc_bi_energy_pt2(E(1:N_states)+pt2(1:N_states))
 end

src/hartree_fock/NEED

@@ -1,3 +1,4 @@
 ao_one_e_ints
 ao_two_e_ints
 scf_utils
+json

src/hartree_fock/fock_matrix_hf.irp.f

@@ -15,16 +15,22 @@
  double precision, allocatable  :: ao_two_e_integral_alpha_tmp(:,:)
  double precision, allocatable  :: ao_two_e_integral_beta_tmp(:,:)
 
+ if (do_ao_cholesky) then   ! Use Cholesky-decomposed integrals
+   ao_two_e_integral_alpha(:,:) = ao_two_e_integral_alpha_chol(:,:)
+   ao_two_e_integral_beta (:,:) = ao_two_e_integral_beta_chol (:,:)
+
+ else ! Use integrals in AO basis set
+
    ao_two_e_integral_alpha = 0.d0
    ao_two_e_integral_beta  = 0.d0
    if (do_direct_integrals) then
 
      !$OMP PARALLEL DEFAULT(NONE)                                    &
-       !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0, &
+         !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,keys,values,p,q,r,s,i0,j0,k0,l0,&
-       !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2, &
+         !$OMP ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp, c0, c1, c2,&
-       !$OMP local_threshold)&
+         !$OMP local_threshold)                                      &
          !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
-       !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz, &
+         !$OMP ao_integrals_map,ao_integrals_threshold, ao_two_e_integral_schwartz,&
          !$OMP ao_two_e_integral_alpha, ao_two_e_integral_beta)
 
      allocate(keys(1), values(1))
@@ -105,7 +111,7 @@
      double precision, allocatable  :: values(:)
 
      !$OMP PARALLEL DEFAULT(NONE)                                    &
-       !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max, &
+         !$OMP PRIVATE(i,j,l,k1,k,integral,ii,jj,kk,ll,i8,keys,values,n_elements_max,&
          !$OMP  n_elements,ao_two_e_integral_alpha_tmp,ao_two_e_integral_beta_tmp)&
          !$OMP SHARED(ao_num,SCF_density_matrix_ao_alpha,SCF_density_matrix_ao_beta,&
          !$OMP  ao_integrals_map, ao_two_e_integral_alpha, ao_two_e_integral_beta)
@@ -150,6 +156,81 @@
      !$OMP END PARALLEL
 
    endif
+ endif
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, ao_two_e_integral_alpha_chol, (ao_num, ao_num) ]
+&BEGIN_PROVIDER [ double precision, ao_two_e_integral_beta_chol ,  (ao_num, ao_num) ]
+ use map_module
+ implicit none
+ BEGIN_DOC
+ ! Alpha and Beta Fock matrices in AO basis set
+ END_DOC
+
+ integer :: m,n,l,s,j
+ double precision :: integral
+ double precision, allocatable :: X(:), X2(:,:,:,:), X3(:,:,:,:)
+
+ allocate (X(cholesky_ao_num))
+
+
+ ! X(j) = \sum_{mn} SCF_density_matrix_ao(m,n) * cholesky_ao(m,n,j)
+ call dgemm('T','N',cholesky_ao_num,1,ao_num*ao_num,1.d0,            &
+     cholesky_ao, ao_num*ao_num,                                     &
+     SCF_density_matrix_ao, ao_num*ao_num,0.d0,                      &
+     X, cholesky_ao_num)
+!
+
+ ! ao_two_e_integral_alpha(m,n) = \sum_{j} cholesky_ao(m,n,j) * X(j)
+ call dgemm('N','N',ao_num*ao_num,1,cholesky_ao_num, 1.d0,           &
+     cholesky_ao, ao_num*ao_num,                                     &
+     X, cholesky_ao_num, 0.d0,                                       &
+     ao_two_e_integral_alpha_chol, ao_num*ao_num)
+
+ deallocate(X)
+
+ ao_two_e_integral_beta_chol = ao_two_e_integral_alpha_chol
+
+
+ allocate(X2(ao_num,ao_num,cholesky_ao_num,2))
+
+! ao_two_e_integral_alpha_chol (l,s) -= cholesky_ao(l,m,j) * SCF_density_matrix_ao_beta (m,n) * cholesky_ao(n,s,j)
+
+ call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0,      &
+     SCF_density_matrix_ao_alpha, ao_num,                            &
+     cholesky_ao, ao_num, 0.d0,                                      &
+     X2(1,1,1,1), ao_num)
+
+ call dgemm('N','N',ao_num,ao_num*cholesky_ao_num,ao_num, 1.d0,      &
+     SCF_density_matrix_ao_beta, ao_num,                             &
+     cholesky_ao, ao_num, 0.d0,                                      &
+     X2(1,1,1,2), ao_num)
+
+ allocate(X3(ao_num,cholesky_ao_num,ao_num,2))
+
+ do s=1,ao_num
+  do j=1,cholesky_ao_num
+   do m=1,ao_num
+    X3(m,j,s,1) = X2(m,s,j,1)
+    X3(m,j,s,2) = X2(m,s,j,2)
+   enddo
+  enddo
+ enddo
+
+ deallocate(X2)
+
+ call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0,     &
+     cholesky_ao, ao_num,                                            &
+     X3(1,1,1,1), ao_num*cholesky_ao_num, 1.d0,                      &
+     ao_two_e_integral_alpha_chol, ao_num)
+
+ call dgemm('N','N',ao_num,ao_num,ao_num*cholesky_ao_num, -1.d0,     &
+     cholesky_ao, ao_num,                                            &
+     X3(1,1,1,2), ao_num*cholesky_ao_num, 1.d0,                      &
+     ao_two_e_integral_beta_chol, ao_num)
+
+ deallocate(X3)
 
 END_PROVIDER
 

src/hartree_fock/scf.irp.f

@@ -80,9 +80,14 @@ subroutine run
 
   mo_label = 'Orthonormalized'
 
-  call Roothaan_Hall_SCF
+  write(json_unit,json_array_open_fmt) 'scf'
-  call ezfio_set_hartree_fock_energy(SCF_energy)
 
+  call Roothaan_Hall_SCF
+
+  write(json_unit,json_array_close_fmtx)
+  call json_close
+
+  call ezfio_set_hartree_fock_energy(SCF_energy)
 end
 
 

src/iterations/EZFIO.cfg

@@ -1,24 +0,0 @@
-[n_iter]
-interface: ezfio
-doc: Number of saved iterations
-type:integer
-default: 1
-
-[n_det_iterations]
-interface: ezfio, provider
-doc: Number of determinants at each iteration
-type: integer
-size: (100)
- 
-[energy_iterations]
-interface: ezfio, provider
-doc: The variational energy at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-
-[pt2_iterations]
-interface: ezfio, provider
-doc: The |PT2| correction at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-

src/iterations/io.irp.f

@@ -1,37 +0,0 @@
-BEGIN_PROVIDER [ integer, n_iter  ]
-  implicit none
-  BEGIN_DOC
-! number of iterations
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  if (mpi_master) then
-
-      double precision :: zeros(N_states,100)
-      integer :: izeros(100)
-      zeros = 0.d0
-      izeros = 0
-      call ezfio_set_iterations_n_iter(0)
-      call ezfio_set_iterations_energy_iterations(zeros)
-      call ezfio_set_iterations_pt2_iterations(zeros)
-      call ezfio_set_iterations_n_det_iterations(izeros)
-      n_iter = 1
-  endif
-  IRP_IF MPI_DEBUG
-    print *,  irp_here, mpi_rank
-    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
-  IRP_ENDIF
-  IRP_IF MPI
-    include 'mpif.h'
-    integer :: ierr
-    call MPI_BCAST( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    if (ierr /= MPI_SUCCESS) then
-      stop 'Unable to read n_iter with MPI'
-    endif
-  IRP_ENDIF
-
-  call write_time(6)
-
-END_PROVIDER
-

src/iterations/iterations.irp.f

@@ -1,42 +1,65 @@
-BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ]
+BEGIN_PROVIDER [ integer, N_iter  ]
   implicit none
   BEGIN_DOC
- ! Extrapolated energy, using E_var = f(PT2) where PT2=0
+! Number of CIPSI iterations
   END_DOC
- integer :: i
+
- do i=1,min(N_states,N_det)
+  N_iter = 0
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, N_iter_max ]
+ implicit none
+ BEGIN_DOC
+ ! Max number of iterations for extrapolations
+ END_DOC
+ N_iter_max = 8
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, energy_iterations , (n_states,N_iter_max) ]
+&BEGIN_PROVIDER [ double precision, pt2_iterations , (n_states,N_iter_max) ]
+&BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter_max,N_states) ]
+  implicit none
+  BEGIN_DOC
+! The energy at each iteration for the extrapolations
+  END_DOC
+
+   energy_iterations = 0.d0
+   pt2_iterations = 0.d0
+   extrapolated_energy = 0.d0
+END_PROVIDER
+
+subroutine increment_n_iter(e, pt2_data)
+  use selection_types
+  implicit none
+  BEGIN_DOC
+! Does what is necessary to increment n_iter
+  END_DOC
+  double precision, intent(in) :: e(*)
+  type(pt2_type), intent(in)   :: pt2_data
+  integer :: k, i
+
+  if (N_det < N_states) return
+
+  if (N_iter < N_iter_max) then
+    N_iter += 1
+  else
+    do k=2,N_iter
+      energy_iterations(1:N_states,k-1) = energy_iterations(1:N_states,k)
+      pt2_iterations(1:N_states,k-1) = pt2_iterations(1:N_states,k)
+    enddo
+  endif
+  energy_iterations(1:N_states,N_iter) = e(1:N_states)
+  pt2_iterations(1:N_states,N_iter) = pt2_data % rpt2(1:N_states)
+
+  if (N_iter < 2) then
+    extrapolated_energy(1,:) = energy_iterations(:,1) + pt2_iterations(:,1)
+    extrapolated_energy(2,:) = energy_iterations(:,2) + pt2_iterations(:,2)
+  else
+    do i=1,N_states
       call extrapolate_data(N_iter,                               &
           energy_iterations(i,1:N_iter),                          &
              pt2_iterations(i,1:N_iter),                          &
           extrapolated_energy(1:N_iter,i))
     enddo
-END_PROVIDER
-
-
-subroutine save_iterations(e_, pt2_,n_)
-  implicit none
-  BEGIN_DOC
-! Update the energy in the EZFIO file.
-  END_DOC
-  integer, intent(in) :: n_
-  double precision, intent(in) :: e_(N_states), pt2_(N_states)
-  integer :: i
-
-  if (N_iter == 101) then
-    do i=2,N_iter-1
-      energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter)
-      pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) 
-    enddo
-    N_iter = N_iter-1
-    TOUCH N_iter
   endif
-
-  energy_iterations(1:N_states,N_iter) = e_(1:N_states)
-     pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states)
-  n_det_iterations(N_iter) = n_
-  call ezfio_set_iterations_N_iter(N_iter)
-  call ezfio_set_iterations_energy_iterations(energy_iterations)
-  call ezfio_set_iterations_pt2_iterations(pt2_iterations)
-  call ezfio_set_iterations_n_det_iterations(n_det_iterations)
 end
-

src/iterations/print_extrapolation.irp.f

@@ -5,10 +5,14 @@ subroutine print_extrapolated_energy
   END_DOC
 
   integer :: i,k
+  integer :: N_states_p, N_iter_p
 
   if (N_iter< 2) then
     return
   endif
+  N_states_p = min(N_states,N_det)
+  N_iter_p = min(N_iter, 8)
+
   write(*,'(A)') ''
   write(*,'(A)') 'Extrapolated energies'
   write(*,'(A)') '------------------------'
@@ -20,20 +24,20 @@ subroutine print_extrapolated_energy
   write(*,*)  '=========== ', '==================='
   write(*,*)  'minimum PT2 ', 'Extrapolated energy'
   write(*,*)  '=========== ', '==================='
-  do k=2,min(N_iter,8)
+  do k=2,N_iter_p
-    write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1)
+    write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,k), extrapolated_energy(k,1)
   enddo
   write(*,*)  '=========== ', '==================='
 
-  do i=2, min(N_states,N_det)
+  do i=2, N_states_p
     print *,  ''
     print *,  'State ', i
     print *,  ''
     write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
     write(*,*)  'minimum PT2 ', 'Extrapolated energy ', '  Excitation (a.u)  ', '  Excitation (eV)  '
     write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-    do k=2,min(N_iter,8)
+    do k=2,N_iter_p
-      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
+      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,k), extrapolated_energy(k,i), &
           extrapolated_energy(k,i) - extrapolated_energy(k,1), &
           (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
     enddo

src/iterations_tc/EZFIO.cfg

@@ -1,24 +0,0 @@
-[n_iter]
-interface: ezfio
-doc: Number of saved iterations
-type:integer
-default: 1
-
-[n_det_iterations]
-interface: ezfio, provider
-doc: Number of determinants at each iteration
-type: integer
-size: (100)
- 
-[energy_iterations]
-interface: ezfio, provider
-doc: The variational energy at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-
-[pt2_iterations]
-interface: ezfio, provider
-doc: The |PT2| correction at each iteration
-type: double precision 
-size: (determinants.n_states,100)
-

src/iterations_tc/io.irp.f

@@ -1,37 +0,0 @@
-BEGIN_PROVIDER [ integer, n_iter  ]
-  implicit none
-  BEGIN_DOC
-! number of iterations
-  END_DOC
-
-  logical                        :: has
-  PROVIDE ezfio_filename
-  if (mpi_master) then
-
-      double precision :: zeros(N_states,100)
-      integer :: izeros(100)
-      zeros = 0.d0
-      izeros = 0
-      call ezfio_set_iterations_n_iter(0)
-      call ezfio_set_iterations_energy_iterations(zeros)
-      call ezfio_set_iterations_pt2_iterations(zeros)
-      call ezfio_set_iterations_n_det_iterations(izeros)
-      n_iter = 1
-  endif
-  IRP_IF MPI_DEBUG
-    print *,  irp_here, mpi_rank
-    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
-  IRP_ENDIF
-  IRP_IF MPI
-    include 'mpif.h'
-    integer :: ierr
-    call MPI_BCAST( n_iter, 1, MPI_INTEGER, 0, MPI_COMM_WORLD, ierr)
-    if (ierr /= MPI_SUCCESS) then
-      stop 'Unable to read n_iter with MPI'
-    endif
-  IRP_ENDIF
-
-  call write_time(6)
-
-END_PROVIDER
-

src/iterations_tc/iterations.irp.f

@@ -1,43 +0,0 @@
-BEGIN_PROVIDER [ double precision, extrapolated_energy, (N_iter,N_states) ]
- implicit none
- BEGIN_DOC
- ! Extrapolated energy, using E_var = f(PT2) where PT2=0
- END_DOC
-! integer :: i
- extrapolated_energy = 0.D0
-END_PROVIDER 
-
- subroutine get_extrapolated_energy(Niter,ept2,pt1,extrap_energy)
- implicit none
- integer, intent(in)  :: Niter
- double precision, intent(in) :: ept2(Niter),pt1(Niter),extrap_energy(Niter)
- call extrapolate_data(Niter,ept2,pt1,extrap_energy)
- end
-
-subroutine save_iterations(e_, pt2_,n_)
-  implicit none
-  BEGIN_DOC
-! Update the energy in the EZFIO file.
-  END_DOC
-  integer, intent(in) :: n_
-  double precision, intent(in) :: e_(N_states), pt2_(N_states)
-  integer :: i
-
-  if (N_iter == 101) then
-    do i=2,N_iter-1
-      energy_iterations(1:N_states,N_iter-1) = energy_iterations(1:N_states,N_iter)
-      pt2_iterations(1:N_states,N_iter-1) = pt2_iterations(1:N_states,N_iter) 
-    enddo
-    N_iter = N_iter-1
-    TOUCH N_iter
-  endif
-
-  energy_iterations(1:N_states,N_iter) = e_(1:N_states)
-     pt2_iterations(1:N_states,N_iter) = pt2_(1:N_states)
-  n_det_iterations(N_iter) = n_
-  call ezfio_set_iterations_N_iter(N_iter)
-  call ezfio_set_iterations_energy_iterations(energy_iterations)
-  call ezfio_set_iterations_pt2_iterations(pt2_iterations)
-  call ezfio_set_iterations_n_det_iterations(n_det_iterations)
-end
-

src/iterations_tc/print_extrapolation.irp.f

@@ -1,46 +0,0 @@
-subroutine print_extrapolated_energy
-  implicit none
-  BEGIN_DOC
-! Print the extrapolated energy in the output
-  END_DOC
-
-  integer :: i,k
-
-  if (N_iter< 2) then
-    return
-  endif
-  write(*,'(A)') ''
-  write(*,'(A)') 'Extrapolated energies'
-  write(*,'(A)') '------------------------'
-  write(*,'(A)') ''
-
-  print *,  ''
-  print *,  'State ', 1
-  print *,  ''
-  write(*,*)  '=========== ', '==================='
-  write(*,*)  'minimum PT2 ', 'Extrapolated energy'
-  write(*,*)  '=========== ', '==================='
-  do k=2,min(N_iter,8)
-    write(*,'(F11.4,2X,F18.8)') pt2_iterations(1,N_iter+1-k), extrapolated_energy(k,1)
-  enddo
-  write(*,*)  '=========== ', '==================='
-
-  do i=2, min(N_states,N_det)
-    print *,  ''
-    print *,  'State ', i
-    print *,  ''
-    write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-    write(*,*)  'minimum PT2 ', 'Extrapolated energy ', '  Excitation (a.u)  ', '  Excitation (eV)  '
-    write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-    do k=2,min(N_iter,8)
-      write(*,'(F11.4,X,3(X,F18.8))') pt2_iterations(i,N_iter+1-k), extrapolated_energy(k,i), &
-          extrapolated_energy(k,i) - extrapolated_energy(k,1), &
-          (extrapolated_energy(k,i) - extrapolated_energy(k,1) ) * 27.211396641308d0
-    enddo
-    write(*,*)  '=========== ', '=================== ', '=================== ', '==================='
-  enddo
-
-  print *,  ''
-
-end subroutine
-

src/iterations_tc/print_summary.irp.f

@@ -1,104 +0,0 @@
-subroutine print_summary(e_,pt2_data,pt2_data_err,n_det_,n_configuration_,n_st,s2_)
-  use selection_types
-  implicit none
-  BEGIN_DOC
-! Print the extrapolated energy in the output
-  END_DOC
-
-  integer, intent(in)            :: n_det_, n_configuration_, n_st
-  double precision, intent(in)   :: e_(n_st), s2_(n_st)
-  type(pt2_type)  , intent(in)   :: pt2_data, pt2_data_err
-  integer                        :: i, k
-  integer                        :: N_states_p
-  character*(9)                  :: pt2_string
-  character*(512)                :: fmt
-
-  if (do_pt2) then
-    pt2_string = '        '
-  else
-    pt2_string = '(approx)'
-  endif
-
-  N_states_p = min(N_det_,n_st)
-
-  print *, ''
-  print '(A,I12)',  'Summary at N_det = ', N_det_
-  print '(A)',      '-----------------------------------'
-  print *, ''
-
-  write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
-  write(*,fmt)
-  write(fmt,*) '(13X,', N_states_p, '(6X,A7,1X,I6,10X))'
-  write(*,fmt) ('State',k, k=1,N_states_p)
-  write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
-  write(*,fmt)
-  write(fmt,*) '(A13,', N_states_p, '(1X,F14.8,15X))'
-  write(*,fmt) '# E          ', e_(1:N_states_p)
-  if (N_states_p > 1) then
-    write(*,fmt) '# Excit. (au)', e_(1:N_states_p)-e_(1)
-    write(*,fmt) '# Excit. (eV)', (e_(1:N_states_p)-e_(1))*27.211396641308d0
-  endif
-  write(fmt,*) '(A13,', 2*N_states_p, '(1X,F14.8))'
-  write(*,fmt) '# PT2 '//pt2_string, (pt2_data % pt2(k), pt2_data_err % pt2(k), k=1,N_states_p)
-  write(*,fmt) '# rPT2'//pt2_string, (pt2_data % rpt2(k), pt2_data_err % rpt2(k), k=1,N_states_p)
-  write(*,'(A)') '#'
-  write(*,fmt) '# E+PT2      ', (e_(k)+pt2_data % pt2(k),pt2_data_err % pt2(k), k=1,N_states_p)
-  write(*,fmt) '# E+rPT2     ', (e_(k)+pt2_data % rpt2(k),pt2_data_err % rpt2(k), k=1,N_states_p)
-  if (N_states_p > 1) then
-    write(*,fmt) '# Excit. (au)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1)), &
-      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1)), k=1,N_states_p)
-    write(*,fmt) '# Excit. (eV)', ( (e_(k)+pt2_data % pt2(k)-e_(1)-pt2_data % pt2(1))*27.211396641308d0, &
-      dsqrt(pt2_data_err % pt2(k)*pt2_data_err % pt2(k)+pt2_data_err % pt2(1)*pt2_data_err % pt2(1))*27.211396641308d0, k=1,N_states_p)
-  endif
-  write(fmt,*) '(''# ============'',', N_states_p, '(1X,''=============================''))'
-  write(*,fmt)
-  print *,  ''
-
-  print *,  'N_det             = ', N_det_
-  print *,  'N_states          = ', n_st
-  if (s2_eig) then
-    print *,  'N_cfg             = ', N_configuration_
-    if (only_expected_s2) then
-      print *,  'N_csf             = ', N_csf
-    endif
-  endif
-  print *,  ''
-
-  do k=1, N_states_p
-    print*,'* State ',k
-    print *,  '< S^2 >         = ', s2_(k)
-    print *,  'E               = ', e_(k)
-    print *,  'Variance        = ', pt2_data % variance(k), ' +/- ', pt2_data_err % variance(k)
-    print *,  'PT norm         = ', dsqrt(pt2_data % overlap(k,k)), ' +/- ', 0.5d0*dsqrt(pt2_data % overlap(k,k)) * pt2_data_err % overlap(k,k) / (pt2_data % overlap(k,k))
-    print *,  'PT2             = ', pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
-    print *,  'rPT2            = ', pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k)
-    print *,  'E+PT2 '//pt2_string//' = ', e_(k)+pt2_data % pt2(k), ' +/- ', pt2_data_err % pt2(k)
-    print *,  'E+rPT2'//pt2_string//' = ', e_(k)+pt2_data % rpt2(k), ' +/- ', pt2_data_err % rpt2(k)
-    print *,  ''
-  enddo
-
-  print *,  '-----'
-  if(n_st.gt.1)then
-    print *, 'Variational Energy difference (au | eV)'
-    do i=2, N_states_p
-      print*,'Delta E = ', (e_(i) - e_(1)), &
-        (e_(i) - e_(1)) * 27.211396641308d0
-    enddo
-    print *,  '-----'
-    print*, 'Variational + perturbative Energy difference (au | eV)'
-    do i=2, N_states_p
-      print*,'Delta E = ', (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))), &
-        (e_(i)+ pt2_data % pt2(i) - (e_(1) + pt2_data % pt2(1))) * 27.211396641308d0
-    enddo
-    print *,  '-----'
-    print*, 'Variational + renormalized perturbative Energy difference (au | eV)'
-    do i=2, N_states_p
-      print*,'Delta E = ', (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))), &
-        (e_(i)+ pt2_data % rpt2(i) - (e_(1) + pt2_data % rpt2(1))) * 27.211396641308d0
-    enddo
-  endif
-
-!  call print_energy_components()
-
-end subroutine
-

src/json/EZFIO.cfg

@@ -0,0 +1,5 @@
+[empty]
+type: logical
+doc: Needed to create the json directory
+interface: ezfio
+

src/json/NEED

@@ -0,0 +1 @@
+ezfio_files

src/json/README.rst

@@ -0,0 +1,5 @@
+====
+json
+====
+
+JSON files to simplify getting output information from QP.

src/json/json.irp.f

@@ -0,0 +1,45 @@
+BEGIN_PROVIDER [ character*(128), json_filename ]
+ implicit none
+ BEGIN_DOC
+ ! Fortran unit of the JSON file
+ END_DOC
+ integer, external            :: getUnitAndOpen
+ integer                      :: counter
+ character*(128)              :: prefix
+ logical                      :: exists
+
+ prefix = trim(ezfio_filename)//'/json/'
+
+ call lock_io
+ exists = .True.
+ counter = 0
+ do while (exists)
+   counter += 1
+   write(json_filename, '(A,I5.5,A)') trim(prefix), counter, '.json'
+   INQUIRE(FILE=trim(json_filename), EXIST=exists)
+ enddo
+ call unlock_io
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, json_unit]
+ implicit none
+ BEGIN_DOC
+ ! Unit file for JSON output
+ END_DOC
+ integer, external :: getUnitAndOpen
+ call ezfio_set_json_empty(.False.)
+ call lock_io
+ json_unit = getUnitAndOpen(json_filename, 'w')
+ write(json_unit, '(A)') '{'
+ call unlock_io
+END_PROVIDER
+
+subroutine json_close
+ call lock_io
+ write(json_unit, '(A)') '}'
+ close(json_unit)
+ call unlock_io
+ FREE json_unit
+end
+

src/json/json_formats.irp.f

@@ -0,0 +1,46 @@
+ BEGIN_PROVIDER [ character*(64), json_int_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_int_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_real_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_real_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_str_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_str_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_true_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_true_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_false_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_false_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_array_open_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_close_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_close_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_array_close_fmtx ]
+&BEGIN_PROVIDER [ character*(64), json_dict_open_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_close_uopen_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_close_fmt ]
+&BEGIN_PROVIDER [ character*(64), json_dict_close_fmtx ]
+ implicit none
+ BEGIN_DOC
+ ! Formats for JSON output.
+ ! x: used to mark the last write (no comma)
+ END_DOC
+ json_int_fmt    = '(''   "'',A,''": '',I10,'','')'
+ json_int_fmtx   = '(''   "'',A,''": '',I10)'
+ json_real_fmt   = '(''   "'',A,''": '',E22.15,'','')'
+ json_real_fmtx  = '(''   "'',A,''": '',E22.15)'
+ json_str_fmt    = '(''   "'',A,''": "'',A,''",'')'
+ json_str_fmtx   = '(''   "'',A,''": "'',A,''"'')'
+ json_true_fmt   = '(''   "'',A,''": true,'')'
+ json_true_fmtx  = '(''   "'',A,''": true'')'
+ json_false_fmt  = '(''   "'',A,''": false,'')'
+ json_false_fmtx = '(''   "'',A,''": false'')'
+ json_array_open_fmt   = '(''   "'',A,''": ['')'
+ json_array_uopen_fmt  = '(''   ['')'
+ json_array_close_fmt  = '(''   ],'')'
+ json_array_close_uopen_fmt  = '(''   ], ['')'
+ json_array_close_fmtx = '(''   ]'')'
+ json_dict_open_fmt    = '(''   "'',A,''": {'')'
+ json_dict_uopen_fmt   = '(''   {'')'
+ json_dict_close_fmt   = '(''   },'')'
+ json_dict_close_uopen_fmt   = '(''   }, {'')'
+ json_dict_close_fmtx  = '(''   }'')'
+END_PROVIDER

src/kohn_sham/ks_scf.irp.f

@@ -90,7 +90,11 @@ subroutine run
 
 ! Choose SCF algorithm
 
+  write(json_unit,*) '"scf" : ['
   call Roothaan_Hall_SCF
+  write(json_unit,*) ']'
+
+  call json_close
 
 end
 

src/kohn_sham_rs/rs_ks_scf.irp.f

@@ -93,7 +93,10 @@ subroutine run
 
   level_shift += 1.d0
   touch level_shift
+  write(json_unit,*) '"scf" : ['
   call Roothaan_Hall_SCF
+  write(json_unit,*) ']'
+  call json_close
   call ezfio_set_kohn_sham_rs_energy(SCF_energy)
 
  write(*, '(A22,X,F16.10)') 'one_e_energy = ',one_e_energy

src/mo_basis/EZFIO.cfg

@@ -9,6 +9,12 @@ doc: Coefficient of the i-th |AO| on the j-th |MO|
 interface: ezfio
 size: (ao_basis.ao_num,mo_basis.mo_num)
 
+[mo_coef_aux]
+type: double precision
+doc: AUX Coefficient of the i-th |AO| on the j-th |MO|
+interface: ezfio
+size: (ao_basis.ao_num,mo_basis.mo_num)
+
 [mo_coef_imag]
 type: double precision
 doc: Imaginary part of the MO coefficient of the i-th |AO| on the j-th |MO|

src/mo_basis/mos_aux.irp.f

@@ -0,0 +1,53 @@
+
+! ---
+
+BEGIN_PROVIDER [double precision, mo_coef_aux, (ao_num,mo_num)]
+
+  implicit none
+  integer                        :: i, j
+  logical                        :: exists
+  double precision, allocatable  :: buffer(:,:)
+
+  PROVIDE ezfio_filename
+
+  if (mpi_master) then
+    ! Coefs
+    call ezfio_has_mo_basis_mo_coef_aux(exists)
+  endif
+
+  IRP_IF MPI_DEBUG
+    print *,  irp_here, mpi_rank
+    call MPI_BARRIER(MPI_COMM_WORLD, ierr)
+  IRP_ENDIF
+
+  IRP_IF MPI
+    include 'mpif.h'
+    integer :: ierr
+    call MPI_BCAST(exists, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read mo_coef_aux with MPI'
+    endif
+  IRP_ENDIF
+
+  if (exists) then
+    if (mpi_master) then
+      call ezfio_get_mo_basis_mo_coef_aux(mo_coef_aux)
+      write(*,*) 'Read  mo_coef_aux'
+    endif
+    IRP_IF MPI
+      call MPI_BCAST(mo_coef_aux, mo_num*ao_num, MPI_DOUBLE_PRECISION, 0, MPI_COMM_WORLD, ierr)
+      if (ierr /= MPI_SUCCESS) then
+        stop 'Unable to read mo_coef_aux with MPI'
+      endif
+    IRP_ENDIF
+  else
+    ! Orthonormalized AO basis
+    do i = 1, mo_num
+      do j = 1, ao_num
+        mo_coef_aux(j,i) = ao_ortho_canonical_coef(j,i)
+      enddo
+    enddo
+  endif
+
+END_PROVIDER
+

src/mo_two_e_ints/cholesky.irp.f

@@ -0,0 +1,16 @@
+BEGIN_PROVIDER [ double precision, cholesky_mo, (mo_num, mo_num, cholesky_ao_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Cholesky vectors in MO basis
+ END_DOC
+
+ integer :: k
+
+ !$OMP PARALLEL DO PRIVATE(k)
+ do k=1,cholesky_ao_num
+  call ao_to_mo(cholesky_ao(1,1,k),ao_num,cholesky_mo(1,1,k),mo_num)
+ enddo
+ !$OMP END PARALLEL DO
+
+END_PROVIDER
+

src/mo_two_e_ints/map_integrals.irp.f

@@ -377,6 +377,7 @@ integer function load_mo_integrals(filename)
   integer*8                      :: n, j
   load_mo_integrals = 1
   open(unit=66,file=filename,FORM='unformatted',STATUS='UNKNOWN')
+  call lock_io()
   read(66,err=98,end=98) iknd, kknd
   if (iknd /= integral_kind) then
     print *,  'Wrong integrals kind in file :', iknd
@@ -399,6 +400,7 @@ integer function load_mo_integrals(filename)
     n = mo_integrals_map%map(i)%n_elements
     read(66,err=99,end=99) (key(j), j=1,n), (val(j), j=1,n)
   enddo
+  call unlock_io()
   call map_sort(mo_integrals_map)
   load_mo_integrals = 0
   return

src/mo_two_e_ints/mo_bi_integrals.irp.f

@@ -50,8 +50,10 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
   call cpu_time(cpu_1)
 
   if(no_vvvv_integrals)then
-!    call four_idx_novvvv
     call four_idx_novvvv_old
+  else
+    if (do_ao_cholesky) then
+      call add_integrals_to_map_cholesky
     else
       if (dble(ao_num)**4 * 32.d-9 < dble(qp_max_mem)) then
         call four_idx_dgemm
@@ -59,6 +61,7 @@ BEGIN_PROVIDER [ logical, mo_two_e_integrals_in_map ]
         call add_integrals_to_map(full_ijkl_bitmask_4)
       endif
     endif
+  endif
 
   call wall_time(wall_2)
   call cpu_time(cpu_2)
@@ -175,7 +178,7 @@ subroutine add_integrals_to_map(mask_ijkl)
   implicit none
 
   BEGIN_DOC
-  ! Adds integrals to tha MO map according to some bitmask
+  ! Adds integrals to the MO map according to some bitmask
   END_DOC
 
   integer(bit_kind), intent(in)  :: mask_ijkl(N_int,4)
@@ -450,13 +453,72 @@ subroutine add_integrals_to_map(mask_ijkl)
 
 end
 
+subroutine add_integrals_to_map_cholesky
+  use bitmasks
+  implicit none
+
+  BEGIN_DOC
+  ! Adds integrals to the MO map using Cholesky vectors
+  END_DOC
+
+  integer :: i,j,k,l,m
+  integer :: size_buffer, n_integrals
+  size_buffer = min(mo_num*mo_num*mo_num,16000000)
+
+  double precision, allocatable :: Vtmp(:,:,:,:)
+  integer(key_kind)  , allocatable :: buffer_i(:)
+  real(integral_kind), allocatable :: buffer_value(:)
+
+  if (.True.) then
+    ! In-memory transformation
+
+    allocate (Vtmp(mo_num,mo_num,mo_num,mo_num))
+
+    call dgemm('N','T',mo_num*mo_num,mo_num*mo_num,cholesky_ao_num,1.d0, &
+       cholesky_mo, mo_num*mo_num, &
+       cholesky_mo, mo_num*mo_num, 0.d0, &
+       Vtmp, mo_num*mo_num)
+
+    !$OMP PARALLEL PRIVATE(i,j,k,l,n_integrals,buffer_value, buffer_i)
+    allocate (buffer_i(size_buffer), buffer_value(size_buffer))
+    n_integrals = 0
+    !$OMP DO
+    do l=1,mo_num
+      do k=1,l
+        do j=1,mo_num
+          do i=1,j
+            if (abs(Vtmp(i,j,k,l)) > mo_integrals_threshold) then
+              n_integrals += 1
+              buffer_value(n_integrals) = Vtmp(i,j,k,l)
+              !DIR$ FORCEINLINE
+              call mo_two_e_integrals_index(i,k,j,l,buffer_i(n_integrals))
+              if (n_integrals == size_buffer) then
+                call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
+                n_integrals = 0
+              endif
+            endif
+          enddo
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+    call map_append(mo_integrals_map, buffer_i, buffer_value, n_integrals)
+    deallocate(buffer_i, buffer_value)
+    !$OMP END PARALLEL
+
+    deallocate(Vtmp)
+    call map_unique(mo_integrals_map)
+
+  endif
+
+end
 
 subroutine add_integrals_to_map_three_indices(mask_ijk)
   use bitmasks
   implicit none
 
   BEGIN_DOC
-  ! Adds integrals to tha MO map according to some bitmask
+  ! Adds integrals to the MO map according to some bitmask
   END_DOC
 
   integer(bit_kind), intent(in)  :: mask_ijk(N_int,3)

src/non_h_ints_mu/debug_fit.irp.f

@@ -18,13 +18,13 @@ program debug_fit
   PROVIDE mu_erf j1b_pen
 
   !call test_j1b_nucl()
-  call test_grad_j1b_nucl()
+  !call test_grad_j1b_nucl()
   !call test_lapl_j1b_nucl()
 
   !call test_list_b2()
-  !call test_list_b3()
+  call test_list_b3()
 
-  call test_fit_u()
+  !call test_fit_u()
   !call test_fit_u2()
   !call test_fit_ugradu()
 
@@ -236,16 +236,25 @@ subroutine test_list_b3()
   integer                    :: ipoint
   double precision           :: acc_ij, acc_tot, eps_ij, i_exc, i_tmp, i_num, normalz
   double precision           :: r(3)
-  double precision, external :: j1b_nucl
+  double precision           :: grad_num(3), eps_der, eps_lap, tmp_der, tmp_lap, i0, ip, im
+  double precision, external :: j1b_nucl_square
 
   print*, ' test_list_b3 ...'
 
+  eps_ij  = 1d-7
+
+  eps_der = 1d-5
+  tmp_der = 0.5d0 / eps_der
+
+  eps_lap = 1d-4
+  tmp_lap = 1.d0 / (eps_lap*eps_lap)
+
+  ! ---
+
   PROVIDE v_1b_list_b3
 
-  eps_ij  = 1d-7
   acc_tot = 0.d0
   normalz = 0.d0
-
   do ipoint = 1, n_points_final_grid
 
     r(1) = final_grid_points(1,ipoint)
@@ -253,8 +262,7 @@ subroutine test_list_b3()
     r(3) = final_grid_points(3,ipoint)
 
     i_exc  = v_1b_list_b3(ipoint) 
-    i_tmp  = j1b_nucl(r)
+    i_num  = j1b_nucl_square(r)
-    i_num  = i_tmp * i_tmp
     acc_ij = dabs(i_exc - i_num)
     if(acc_ij .gt. eps_ij) then
       print *, ' problem in list_b3 on', ipoint
@@ -267,8 +275,136 @@ subroutine test_list_b3()
     normalz += dabs(i_num)
   enddo
 
-  print*, ' acc_tot = ', acc_tot
+  print*, ' acc_tot on val = ', acc_tot
-  print*, ' normalz = ', normalz
+  print*, ' normalz on val = ', normalz
+
+  ! ---
+
+  PROVIDE v_1b_square_grad
+
+  acc_tot = 0.d0
+  normalz = 0.d0
+  do ipoint = 1, n_points_final_grid
+
+    r(1) = final_grid_points(1,ipoint)
+    r(2) = final_grid_points(2,ipoint)
+    r(3) = final_grid_points(3,ipoint)
+
+    i_exc  = v_1b_square_grad(ipoint,1)
+    r(1)   = r(1) + eps_der
+    ip     = j1b_nucl_square(r)
+    r(1)   = r(1) - 2.d0 * eps_der
+    im     = j1b_nucl_square(r)
+    r(1)   = r(1) + eps_der
+    i_num  = tmp_der * (ip - im)
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in grad_x list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+
+    i_exc  = v_1b_square_grad(ipoint,2)
+    r(2)   = r(2) + eps_der
+    ip     = j1b_nucl_square(r)
+    r(2)   = r(2) - 2.d0 * eps_der
+    im     = j1b_nucl_square(r)
+    r(2)   = r(2) + eps_der
+    i_num  = tmp_der * (ip - im)
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in grad_y list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+
+    i_exc  = v_1b_square_grad(ipoint,3)
+    r(3)   = r(3) + eps_der
+    ip     = j1b_nucl_square(r)
+    r(3)   = r(3) - 2.d0 * eps_der
+    im     = j1b_nucl_square(r)
+    r(3)   = r(3) + eps_der
+    i_num  = tmp_der * (ip - im)
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in grad_z list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+  enddo
+
+  print*, ' acc_tot on grad = ', acc_tot
+  print*, ' normalz on grad = ', normalz
+
+  ! ---
+
+  PROVIDE v_1b_square_lapl
+
+  acc_tot = 0.d0
+  normalz = 0.d0
+  do ipoint = 1, n_points_final_grid
+
+    r(1) = final_grid_points(1,ipoint)
+    r(2) = final_grid_points(2,ipoint)
+    r(3) = final_grid_points(3,ipoint)
+    i0   = j1b_nucl_square(r)
+
+    i_exc = v_1b_square_lapl(ipoint)
+
+    r(1)  = r(1) + eps_lap
+    ip    = j1b_nucl_square(r)
+    r(1)  = r(1) - 2.d0 * eps_lap
+    im    = j1b_nucl_square(r)
+    r(1)  = r(1) + eps_lap
+    i_num = tmp_lap * (ip - 2.d0 * i0 + im)
+
+    r(2)  = r(2) + eps_lap
+    ip    = j1b_nucl_square(r)
+    r(2)  = r(2) - 2.d0 * eps_lap
+    im    = j1b_nucl_square(r)
+    r(2)  = r(2) + eps_lap
+    i_num = i_num + tmp_lap * (ip - 2.d0 * i0 + im)
+
+    r(3)  = r(3) + eps_lap
+    ip    = j1b_nucl_square(r)
+    r(3)  = r(3) - 2.d0 * eps_lap
+    im    = j1b_nucl_square(r)
+    r(3)  = r(3) + eps_lap
+    i_num = i_num + tmp_lap * (ip - 2.d0 * i0 + im)
+
+    acc_ij = dabs(i_exc - i_num)
+    if(acc_ij .gt. eps_ij) then
+      print *, ' problem in lapl list_b3 on', ipoint
+      print *, ' r      = ', r
+      print *, ' r2     = ', r(1)*r(1) + r(2)*r(2) + r(3)*r(3)
+      print *, ' analyt = ', i_exc
+      print *, ' numeri = ', i_num
+      print *, ' diff   = ', acc_ij
+    endif
+
+    acc_tot += acc_ij
+    normalz += dabs(i_num)
+  enddo
+
+  print*, ' acc_tot on lapl = ', acc_tot
+  print*, ' normalz on lapl = ', normalz
+
+  ! ---
 
   return
 end subroutine test_list_b3

src/non_h_ints_mu/grad_squared.irp.f

@@ -17,7 +17,7 @@ BEGIN_PROVIDER [ double precision, gradu_squared_u_ij_mu, (ao_num, ao_num, n_poi
   !
   ! if J(r1,r2) = u12 x v1 x v2
   !
-  ! 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) ]
+  ! 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
@@ -267,7 +267,7 @@ END_PROVIDER
 
 ! ---
 
-BEGIN_PROVIDER [ double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid) ]
+BEGIN_PROVIDER [double precision, u12sq_j1bsq, (ao_num, ao_num, n_points_final_grid)]
 
   implicit none
   integer                    :: ipoint, i, j
@@ -358,7 +358,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   implicit none
   integer                       :: ipoint, i, j, k, l
-  double precision              :: weight1, ao_ik_r, ao_i_r
+  double precision              :: weight1, ao_k_r, ao_i_r
+  double precision              :: der_envsq_x, der_envsq_y, der_envsq_z, lap_envsq
   double precision              :: time0, time1
   double precision, allocatable :: b_mat(:,:,:), tmp(:,:,:)
 
@@ -373,6 +374,8 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
 
   else
 
+    ! ---
+
     PROVIDE int2_grad1_u12_square_ao
 
     allocate(b_mat(n_points_final_grid,ao_num,ao_num))
@@ -397,6 +400,50 @@ BEGIN_PROVIDER [double precision, tc_grad_square_ao, (ao_num, ao_num, ao_num, ao
     call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0                 &
               , int2_grad1_u12_square_ao(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
               , 0.d0, tc_grad_square_ao, ao_num*ao_num)
+
+    ! ---
+
+    if((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
+
+      ! an additional term is added here directly instead of 
+      ! being added in int2_grad1_u12_square_ao for performance
+      ! note that the factor
+
+      PROVIDE int2_u2_j1b2
+
+      b_mat = 0.d0
+      !$OMP PARALLEL                                                                                     &
+      !$OMP DEFAULT (NONE)                                                                               &
+      !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r)                                              &
+      !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector, &
+      !$OMP         v_1b_square_grad, v_1b_square_lapl, aos_grad_in_r_array_transp_bis)
+      !$OMP DO SCHEDULE (static)
+      do i = 1, ao_num
+        do k = 1, ao_num
+          do ipoint = 1, n_points_final_grid
+
+            weight1 = 0.25d0 * 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) = weight1 * ( ao_k_r * ao_i_r * v_1b_square_lapl(ipoint)                                                                                   &
+                              + (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)) * v_1b_square_grad(ipoint,1) &
+                              + (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)) * v_1b_square_grad(ipoint,2) &
+                              + (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)) * v_1b_square_grad(ipoint,3) )
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+      call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0     &
+                , int2_u2_j1b2(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)
+    endif
+
+    ! ---
+
     deallocate(b_mat)
 
     call sum_A_At(tc_grad_square_ao(1,1,1,1), ao_num*ao_num)

src/non_h_ints_mu/j12_nucl_utils.irp.f

@@ -8,6 +8,10 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
   double precision :: x, y, z, dx, dy, dz
   double precision :: a, d, e, fact_r
 
+  if(j1b_type .eq. 3) then
+
+    ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)]
+
     do ipoint = 1, n_points_final_grid
 
       x = final_grid_points(1,ipoint)
@@ -29,19 +33,56 @@ BEGIN_PROVIDER [ double precision, v_1b, (n_points_final_grid)]
       v_1b(ipoint) = fact_r
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2)
+
+    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)
+
+      fact_r = 1.d0
+      do j = 1, nucl_num
+        a  = j1b_pen(j)
+        dx = x - nucl_coord(j,1)
+        dy = y - nucl_coord(j,2)
+        dz = z - nucl_coord(j,3)
+        d  = dx*dx + dy*dy + dz*dz
+
+        fact_r = fact_r - dexp(-a*d)
+      enddo
+
+      v_1b(ipoint) = fact_r
+    enddo
+
+  else
+
+    print*, 'j1b_type = ', j1b_pen, 'is not implemented for v_1b'
+    stop
+
+  endif
+
 END_PROVIDER
 
 ! ---
 
-BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)]
+BEGIN_PROVIDER [double precision, v_1b_grad, (3, n_points_final_grid)]
 
   implicit none
   integer          :: ipoint, i, j, phase
-  double precision :: x, y, z, dx, dy, dz
+  double precision :: x, y, z, dx, dy, dz, r2
   double precision :: a, d, e
   double precision :: fact_x, fact_y, fact_z
   double precision :: ax_der, ay_der, az_der, a_expo
 
+  PROVIDE j1b_type
+
+  if(j1b_type .eq. 3) then
+
+    ! v(r) = \Pi_{a} [1 - \exp(-\alpha_a (r - r_a)^2)]
+
     do ipoint = 1, n_points_final_grid
 
       x = final_grid_points(1,ipoint)
@@ -82,6 +123,46 @@ BEGIN_PROVIDER [ double precision, v_1b_grad, (3, n_points_final_grid)]
       v_1b_grad(3,ipoint) = fact_z
     enddo
 
+  elseif(j1b_type .eq. 4) then
+
+    ! v(r) = 1 - \sum_{a} \exp(-\alpha_a (r - r_a)^2)
+
+    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)
+
+      ax_der = 0.d0 
+      ay_der = 0.d0 
+      az_der = 0.d0 
+      do j = 1, nucl_num
+
+        dx = x - nucl_coord(j,1)
+        dy = y - nucl_coord(j,2)
+        dz = z - nucl_coord(j,3)
+        r2 = dx*dx + dy*dy + dz*dz
+
+        a = j1b_pen(j)
+        e = a * dexp(-a * r2)
+
+        ax_der += e * dx
+        ay_der += e * dy
+        az_der += e * dz
+      enddo
+
+      v_1b_grad(1,ipoint) = 2.d0 * ax_der
+      v_1b_grad(2,ipoint) = 2.d0 * ay_der
+      v_1b_grad(3,ipoint) = 2.d0 * az_der
+    enddo
+
+  else
+
+    print*, 'j1b_type = ', j1b_pen, 'is not implemented'
+    stop
+
+  endif
+
 END_PROVIDER
 
 ! ---
@@ -91,7 +172,7 @@ BEGIN_PROVIDER [ double precision, v_1b_lapl, (n_points_final_grid)]
   implicit none
   integer          :: ipoint, i, j, phase
   double precision :: x, y, z, dx, dy, dz
-  double precision :: a, d, e, b
+  double precision :: a, e, b
   double precision :: fact_r
   double precision :: ax_der, ay_der, az_der, a_expo
 
@@ -202,6 +283,56 @@ BEGIN_PROVIDER [ double precision, v_1b_list_b3, (n_points_final_grid)]
 
 END_PROVIDER
 
+! ---
+
+ BEGIN_PROVIDER [double precision, v_1b_square_grad, (n_points_final_grid,3)]
+&BEGIN_PROVIDER [double precision, v_1b_square_lapl, (n_points_final_grid)  ]
+
+  implicit none
+  integer          :: ipoint, i
+  double precision :: x, y, z, dx, dy, dz, r2
+  double precision :: coef, expo, a_expo, tmp
+  double precision :: fact_x, fact_y, fact_z, fact_r
+
+  PROVIDE List_all_comb_b3_coef List_all_comb_b3_expo List_all_comb_b3_cent
+
+  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)
+
+    fact_x = 0.d0
+    fact_y = 0.d0
+    fact_z = 0.d0
+    fact_r = 0.d0
+    do i = 1, List_all_comb_b3_size
+
+      coef = List_all_comb_b3_coef(i)
+      expo = List_all_comb_b3_expo(i)
+
+      dx = x - List_all_comb_b3_cent(1,i)
+      dy = y - List_all_comb_b3_cent(2,i)
+      dz = z - List_all_comb_b3_cent(3,i)
+      r2 = dx * dx + dy * dy + dz * dz
+
+      a_expo = expo * r2
+      tmp    = coef * expo * dexp(-a_expo)
+
+      fact_x += tmp * dx
+      fact_y += tmp * dy
+      fact_z += tmp * dz
+      fact_r += tmp * (3.d0 - 2.d0 * a_expo)
+    enddo
+
+    v_1b_square_grad(ipoint,1) = -2.d0 * fact_x
+    v_1b_square_grad(ipoint,2) = -2.d0 * fact_y
+    v_1b_square_grad(ipoint,3) = -2.d0 * fact_z
+    v_1b_square_lapl(ipoint)   = -2.d0 * fact_r
+  enddo
+
+END_PROVIDER
+
 ! ---
 
 double precision function j12_mu_r12(r12)

src/non_h_ints_mu/jast_deriv.irp.f

@@ -21,6 +21,10 @@
   implicit none
   integer                    :: ipoint, jpoint
   double precision           :: r1(3), r2(3)
+  double precision           :: v1b_r1, v1b_r2, u2b_r12
+  double precision           :: grad1_v1b(3), grad1_u2b(3)
+  double precision           :: dx, dy, dz
+  double precision, external :: j12_mu, j1b_nucl
 
   PROVIDE j1b_type
   PROVIDE final_grid_points_extra
@@ -28,12 +32,43 @@
   grad1_u12_num         = 0.d0
   grad1_u12_squared_num = 0.d0
 
-  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+  if(j1b_type .eq. 100) then
 
-    double precision           :: v1b_r1, v1b_r2, u2b_r12
+    !$OMP PARALLEL                                                                                    &
-    double precision           :: grad1_v1b(3), grad1_u2b(3)
+    !$OMP DEFAULT (NONE)                                                                              &
-    double precision           :: dx, dy, dz
+    !$OMP PRIVATE (ipoint, jpoint, r1, r2, v1b_r1, v1b_r2, u2b_r12, grad1_v1b, grad1_u2b, dx, dy, dz) &
-    double precision, external :: j12_mu, j1b_nucl
+    !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points,                  &
+    !$OMP         final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num)
+    !$OMP DO SCHEDULE (static)
+    do ipoint = 1, n_points_final_grid  ! r1
+
+      r1(1) = final_grid_points(1,ipoint)
+      r1(2) = final_grid_points(2,ipoint)
+      r1(3) = final_grid_points(3,ipoint)
+
+      do jpoint = 1, n_points_extra_final_grid ! r2 
+
+        r2(1) = final_grid_points_extra(1,jpoint)
+        r2(2) = final_grid_points_extra(2,jpoint)
+        r2(3) = final_grid_points_extra(3,jpoint)
+
+        call grad1_j12_mu(r1, r2, grad1_u2b)
+
+        dx = grad1_u2b(1)
+        dy = grad1_u2b(2)
+        dz = grad1_u2b(3)
+
+        grad1_u12_num(jpoint,ipoint,1) = dx 
+        grad1_u12_num(jpoint,ipoint,2) = dy 
+        grad1_u12_num(jpoint,ipoint,3) = dz 
+
+        grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+  elseif((j1b_type .gt. 100) .and. (j1b_type .lt. 200)) then
 
     !$OMP PARALLEL                                                                                    &
     !$OMP DEFAULT (NONE)                                                                              &
@@ -74,6 +109,42 @@
     !$OMP END DO
     !$OMP END PARALLEL
 
+  elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then
+
+    !$OMP PARALLEL                                                                   &
+    !$OMP DEFAULT (NONE)                                                             &
+    !$OMP PRIVATE (ipoint, jpoint, r1, r2, grad1_u2b, dx, dy, dz)                    &
+    !$OMP SHARED (n_points_final_grid, n_points_extra_final_grid, final_grid_points, &
+    !$OMP         final_grid_points_extra, grad1_u12_num, grad1_u12_squared_num)
+    !$OMP DO SCHEDULE (static)
+    do ipoint = 1, n_points_final_grid  ! r1
+
+      r1(1) = final_grid_points(1,ipoint)
+      r1(2) = final_grid_points(2,ipoint)
+      r1(3) = final_grid_points(3,ipoint)
+
+      do jpoint = 1, n_points_extra_final_grid ! r2 
+
+        r2(1) = final_grid_points_extra(1,jpoint)
+        r2(2) = final_grid_points_extra(2,jpoint)
+        r2(3) = final_grid_points_extra(3,jpoint)
+
+        call grad1_j12_mu(r1, r2, grad1_u2b)
+
+        dx = grad1_u2b(1)
+        dy = grad1_u2b(2)
+        dz = grad1_u2b(3)
+
+        grad1_u12_num(jpoint,ipoint,1) = dx 
+        grad1_u12_num(jpoint,ipoint,2) = dy 
+        grad1_u12_num(jpoint,ipoint,3) = dz 
+
+        grad1_u12_squared_num(jpoint,ipoint) = dx*dx + dy*dy + dz*dz
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
   else
 
     print *, ' j1b_type = ', j1b_type, 'not implemented yet'
@@ -91,20 +162,20 @@ double precision function j12_mu(r1, r2)
 
   implicit none
   double precision, intent(in) :: r1(3), r2(3)
-  double precision             :: mu_r12, r12
+  double precision             :: mu_tmp, r12
 
-  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+  if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then
 
     r12 = dsqrt( (r1(1) - r2(1)) * (r1(1) - r2(1)) &
                + (r1(2) - r2(2)) * (r1(2) - r2(2)) &
                + (r1(3) - r2(3)) * (r1(3) - r2(3)) )
-    mu_r12 = mu_erf * r12
+    mu_tmp = mu_erf * r12
 
-    j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_r12)) - inv_sq_pi_2 * dexp(-mu_r12*mu_r12) / mu_erf
+    j12_mu = 0.5d0 * r12 * (1.d0 - derf(mu_tmp)) - inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / mu_erf
 
   else
 
-    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    print *, ' j1b_type = ', j1b_type, 'not implemented for j12_mu'
     stop
 
   endif
@@ -116,6 +187,8 @@ end function j12_mu
 
 subroutine grad1_j12_mu(r1, r2, grad)
 
+  include 'constants.include.F'
+
   implicit none
   double precision, intent(in)  :: r1(3), r2(3)
   double precision, intent(out) :: grad(3)
@@ -123,13 +196,13 @@ subroutine grad1_j12_mu(r1, r2, grad)
 
   grad = 0.d0
 
-  if((j1b_type .ge. 100) .and. (j1b_type .lt. 200)) then
+  if((j1b_type .ge. 0) .and. (j1b_type .lt. 200)) then
 
     dx = r1(1) - r2(1)
     dy = r1(2) - r2(2)
     dz = r1(3) - r2(3)
 
-    r12 = dsqrt( dx * dx + dy * dy + dz * dz )
+    r12 = dsqrt(dx * dx + dy * dy + dz * dz)
     if(r12 .lt. 1d-10) return
 
     tmp = 0.5d0 * (1.d0 - derf(mu_erf * r12)) / r12
@@ -138,6 +211,28 @@ subroutine grad1_j12_mu(r1, r2, grad)
     grad(2) = tmp * dy
     grad(3) = tmp * dz
 
+  elseif((j1b_type .ge. 200) .and. (j1b_type .lt. 300)) then
+
+    double precision :: mu_val, mu_tmp, mu_der(3)
+
+    dx  = r1(1) - r2(1)
+    dy  = r1(2) - r2(2)
+    dz  = r1(3) - r2(3)
+    r12 = dsqrt(dx * dx + dy * dy + dz * dz)
+
+    call mu_r_val_and_grad(r1, r2, mu_val, mu_der)
+    mu_tmp  = mu_val * r12
+    tmp     = inv_sq_pi_2 * dexp(-mu_tmp*mu_tmp) / (mu_val * mu_val)
+    grad(1) = tmp * mu_der(1)
+    grad(2) = tmp * mu_der(2)
+    grad(3) = tmp * mu_der(3)
+
+    if(r12 .lt. 1d-10) return
+    tmp     = 0.5d0 * (1.d0 - derf(mu_tmp)) / r12
+    grad(1) = grad(1) + tmp * dx
+    grad(2) = grad(2) + tmp * dy
+    grad(3) = grad(3) + tmp * dz
+
   else
 
     print *, ' j1b_type = ', j1b_type, 'not implemented yet'
@@ -157,7 +252,18 @@ double precision function j1b_nucl(r)
   integer                      :: i
   double precision             :: a, d, e, x, y, z
 
-  if(j1b_type .eq. 103) then
+  if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then
+
+    j1b_nucl = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      j1b_nucl = j1b_nucl - dexp(-a*dsqrt(d))
+    enddo
+
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then
 
     j1b_nucl = 1.d0
     do i = 1, nucl_num
@@ -169,7 +275,7 @@ double precision function j1b_nucl(r)
       j1b_nucl = j1b_nucl * e
     enddo
 
-  elseif(j1b_type .eq. 104) then
+  elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then
 
     j1b_nucl = 1.d0
     do i = 1, nucl_num
@@ -180,7 +286,7 @@ double precision function j1b_nucl(r)
       j1b_nucl = j1b_nucl - dexp(-a*d)
     enddo
 
-  elseif(j1b_type .eq. 105) then
+  elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
 
     j1b_nucl = 1.d0
     do i = 1, nucl_num
@@ -194,7 +300,7 @@ double precision function j1b_nucl(r)
 
   else
 
-    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl'
     stop
 
   endif
@@ -204,6 +310,75 @@ end function j1b_nucl
 
 ! ---
 
+double precision function j1b_nucl_square(r)
+
+  implicit none
+  double precision, intent(in) :: r(3)
+  integer                      :: i
+  double precision             :: a, d, e, x, y, z
+
+  if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      j1b_nucl_square = j1b_nucl_square - dexp(-a*dsqrt(d))
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      e = 1.d0 - dexp(-a*d)
+      j1b_nucl_square = j1b_nucl_square * e
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      d = ( (r(1) - nucl_coord(i,1)) * (r(1) - nucl_coord(i,1)) &
+          + (r(2) - nucl_coord(i,2)) * (r(2) - nucl_coord(i,2)) &
+          + (r(3) - nucl_coord(i,3)) * (r(3) - nucl_coord(i,3)) )
+      j1b_nucl_square = j1b_nucl_square - dexp(-a*d)
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
+
+    j1b_nucl_square = 1.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      x = r(1) - nucl_coord(i,1)
+      y = r(2) - nucl_coord(i,2)
+      z = r(3) - nucl_coord(i,3)
+      d = x*x + y*y + z*z
+      j1b_nucl_square = j1b_nucl_square - dexp(-a*d*d)
+    enddo
+    j1b_nucl_square = j1b_nucl_square * j1b_nucl_square
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented for j1b_nucl_square'
+    stop
+
+  endif
+
+  return
+end function j1b_nucl_square
+
+! ---
+
 subroutine grad1_j1b_nucl(r, grad)
 
   implicit none
@@ -215,7 +390,29 @@ subroutine grad1_j1b_nucl(r, grad)
   double precision              :: fact_x, fact_y, fact_z
   double precision              :: ax_der, ay_der, az_der, a_expo
 
-  if(j1b_type .eq. 103) then
+  if((j1b_type .eq. 2) .or. (j1b_type .eq. 102)) then
+
+    fact_x = 0.d0
+    fact_y = 0.d0
+    fact_z = 0.d0
+    do i = 1, nucl_num
+      a = j1b_pen(i)
+      x = r(1) - nucl_coord(i,1)
+      y = r(2) - nucl_coord(i,2)
+      z = r(3) - nucl_coord(i,3)
+      d = dsqrt(x*x + y*y + z*z)
+      e = a * dexp(-a*d) / d
+
+      fact_x += e * x
+      fact_y += e * y
+      fact_z += e * z
+    enddo
+
+    grad(1) = fact_x
+    grad(2) = fact_y
+    grad(3) = fact_z
+
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 103)) then
 
     x = r(1)
     y = r(2)
@@ -254,7 +451,7 @@ subroutine grad1_j1b_nucl(r, grad)
     grad(2) = fact_y
     grad(3) = fact_z
 
-  else if(j1b_type .eq. 104) then
+  elseif((j1b_type .eq. 4) .or. (j1b_type .eq. 104)) then
 
     fact_x = 0.d0
     fact_y = 0.d0
@@ -276,7 +473,7 @@ subroutine grad1_j1b_nucl(r, grad)
     grad(2) = 2.d0 * fact_y
     grad(3) = 2.d0 * fact_z
 
-  else if(j1b_type .eq. 105) then
+  elseif((j1b_type .eq. 5) .or. (j1b_type .eq. 105)) then
 
     fact_x = 0.d0
     fact_y = 0.d0
@@ -300,7 +497,7 @@ subroutine grad1_j1b_nucl(r, grad)
 
   else
 
-    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    print *, ' j1b_type = ', j1b_type, 'not implemented for grad1_j1b_nucl'
     stop
 
   endif
@@ -310,3 +507,180 @@ end subroutine grad1_j1b_nucl
 
 ! ---
 
+subroutine mu_r_val_and_grad(r1, r2, mu_val, mu_der)
+
+  implicit none
+  double precision, intent(in)  :: r1(3), r2(3)
+  double precision, intent(out) :: mu_val, mu_der(3)
+  double precision              :: r(3)
+  double precision              :: dm_a(1), dm_b(1), grad_dm_a(3,1), grad_dm_b(3,1)
+  double precision              :: dm_tot, tmp1, tmp2, tmp3
+
+  if(j1b_type .eq. 200) then
+
+    !
+    ! r = 0.5 (r1 + r2)
+    !
+    ! mu[rho(r)] = alpha sqrt(rho) + mu0 exp(-rho)
+    !
+    ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
+    ! d mu[rho(r)] / dx  = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
+    !
+
+    PROVIDE mu_r_ct mu_erf
+
+    r(1) = 0.5d0 * (r1(1) + r2(1))
+    r(2) = 0.5d0 * (r1(2) + r2(2))
+    r(3) = 0.5d0 * (r1(3) + r2(3))
+
+    call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)
+
+    dm_tot = dm_a(1) + dm_b(1)
+    tmp1   = dsqrt(dm_tot)
+    tmp2   = mu_erf * dexp(-dm_tot)
+
+    mu_val = mu_r_ct * tmp1 + tmp2
+
+    mu_der = 0.d0
+    if(dm_tot .lt. 1d-7) return
+
+    tmp3      = 0.25d0 * mu_r_ct / tmp1 - 0.5d0 * tmp2
+    mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
+    mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
+    mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))
+
+  elseif(j1b_type .eq. 201) then
+
+    !
+    ! r = 0.5 (r1 + r2)
+    !
+    ! mu[rho(r)] = alpha rho + mu0 exp(-rho)
+    !
+    ! d mu[rho(r)] / dx1 = 0.5 d mu[rho(r)] / dx
+    ! d mu[rho(r)] / dx  = [0.5 alpha / sqrt(rho) - mu0 exp(-rho)] (d rho(r) / dx)
+    !
+
+    PROVIDE mu_r_ct mu_erf
+
+    r(1) = 0.5d0 * (r1(1) + r2(1))
+    r(2) = 0.5d0 * (r1(2) + r2(2))
+    r(3) = 0.5d0 * (r1(3) + r2(3))
+
+    call density_and_grad_alpha_beta(r, dm_a, dm_b, grad_dm_a, grad_dm_b)
+
+    dm_tot = dm_a(1) + dm_b(1)
+    tmp2   = mu_erf * dexp(-dm_tot)
+
+    mu_val = mu_r_ct * dm_tot + tmp2
+
+    tmp3      = 0.5d0 * (mu_r_ct - tmp2)
+    mu_der(1) = tmp3 * (grad_dm_a(1,1) + grad_dm_b(1,1))
+    mu_der(2) = tmp3 * (grad_dm_a(2,1) + grad_dm_b(2,1))
+    mu_der(3) = tmp3 * (grad_dm_a(3,1) + grad_dm_b(3,1))
+
+  else
+
+    print *, ' j1b_type = ', j1b_type, 'not implemented yet'
+    stop
+
+  endif
+
+  return
+end subroutine mu_r_val_and_grad
+
+! ---
+
+subroutine grad1_j1b_nucl_square_num(r1, grad)
+
+  implicit none
+  double precision, intent(in)  :: r1(3)
+  double precision, intent(out) :: grad(3)
+  double precision              :: r(3), eps, tmp_eps, vp, vm
+  double precision, external    :: j1b_nucl_square
+
+  eps     = 1d-5
+  tmp_eps = 0.5d0 / eps
+
+  r(1:3) = r1(1:3)
+
+  r(1) = r(1) + eps
+  vp   = j1b_nucl_square(r)
+  r(1) = r(1) - 2.d0 * eps
+  vm   = j1b_nucl_square(r)
+  r(1) = r(1) + eps
+  grad(1) = tmp_eps * (vp - vm)
+
+  r(2) = r(2) + eps
+  vp   = j1b_nucl_square(r)
+  r(2) = r(2) - 2.d0 * eps
+  vm   = j1b_nucl_square(r)
+  r(2) = r(2) + eps
+  grad(2) = tmp_eps * (vp - vm)
+
+  r(3) = r(3) + eps
+  vp   = j1b_nucl_square(r)
+  r(3) = r(3) - 2.d0 * eps
+  vm   = j1b_nucl_square(r)
+  r(3) = r(3) + eps
+  grad(3) = tmp_eps * (vp - vm)
+  
+  return
+end subroutine grad1_j1b_nucl_square_num
+
+! ---
+
+subroutine grad1_j12_mu_square_num(r1, r2, grad)
+
+  include 'constants.include.F'
+
+  implicit none
+  double precision, intent(in)  :: r1(3), r2(3)
+  double precision, intent(out) :: grad(3)
+  double precision              :: r(3)
+  double precision              :: eps, tmp_eps, vp, vm
+  double precision, external    :: j12_mu_square
+
+  eps     = 1d-5
+  tmp_eps = 0.5d0 / eps
+
+  r(1:3) = r1(1:3)
+
+  r(1)    = r(1) + eps
+  vp      = j12_mu_square(r, r2)
+  r(1)    = r(1) - 2.d0 * eps
+  vm      = j12_mu_square(r, r2)
+  r(1)    = r(1) + eps
+  grad(1) = tmp_eps * (vp - vm)
+
+  r(2)    = r(2) + eps
+  vp      = j12_mu_square(r, r2)
+  r(2)    = r(2) - 2.d0 * eps
+  vm      = j12_mu_square(r, r2)
+  r(2)    = r(2) + eps
+  grad(2) = tmp_eps * (vp - vm)
+
+  r(3)    = r(3) + eps
+  vp      = j12_mu_square(r, r2)
+  r(3)    = r(3) - 2.d0 * eps
+  vm      = j12_mu_square(r, r2)
+  r(3)    = r(3) + eps
+  grad(3) = tmp_eps * (vp - vm)
+
+  return
+end subroutine grad1_j12_mu_square_num
+
+! ---
+
+double precision function j12_mu_square(r1, r2)
+
+  implicit none
+  double precision, intent(in) :: r1(3), r2(3)
+  double precision, external   :: j12_mu
+
+  j12_mu_square = j12_mu(r1, r2) * j12_mu(r1, r2)
+
+  return
+end function j12_mu_square
+
+! ---
+

src/non_h_ints_mu/tc_integ.irp.f

@@ -35,13 +35,40 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
   PROVIDE j1b_type
 
   if(read_tc_integ) then
+
     open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="read")
     read(11) int2_grad1_u12_ao
-  endif
 
-  if(j1b_type .eq. 3) then
+  else
 
-    if(.not.read_tc_integ) then
+    if(j1b_type .eq. 0) then
+
+      PROVIDE v_ij_erf_rk_cst_mu x_v_ij_erf_rk_cst_mu
+
+      int2_grad1_u12_ao = 0.d0
+      !$OMP PARALLEL                                                &
+      !$OMP DEFAULT (NONE)                                          &
+      !$OMP PRIVATE (ipoint, i, j, x, y, z, tmp1)                   &
+      !$OMP SHARED ( ao_num, n_points_final_grid, final_grid_points &
+      !$OMP        , v_ij_erf_rk_cst_mu, x_v_ij_erf_rk_cst_mu, int2_grad1_u12_ao)
+      !$OMP DO SCHEDULE (static)
+      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)
+        do j = 1, ao_num
+          do i = 1, ao_num
+            tmp1 = v_ij_erf_rk_cst_mu(i,j,ipoint)
+            int2_grad1_u12_ao(i,j,ipoint,1) = 0.5d0 * (tmp1 * x - x_v_ij_erf_rk_cst_mu(i,j,ipoint,1))
+            int2_grad1_u12_ao(i,j,ipoint,2) = 0.5d0 * (tmp1 * y - x_v_ij_erf_rk_cst_mu(i,j,ipoint,2))
+            int2_grad1_u12_ao(i,j,ipoint,3) = 0.5d0 * (tmp1 * z - x_v_ij_erf_rk_cst_mu(i,j,ipoint,3))
+          enddo
+        enddo
+      enddo
+      !$OMP END DO
+      !$OMP END PARALLEL
+
+    elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
 
       PROVIDE v_1b_grad v_ij_erf_rk_cst_mu_j1b v_ij_u_cst_mu_j1b x_v_ij_erf_rk_cst_mu_j1b
 
@@ -73,8 +100,6 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
       !$OMP END DO
       !$OMP END PARALLEL
 
-    endif
-
     elseif(j1b_type .ge. 100) then
 
       PROVIDE final_weight_at_r_vector_extra aos_in_r_array_extra
@@ -98,7 +123,6 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
       !$OMP END DO
       !$OMP END PARALLEL
 
-    if(.not.read_tc_integ) then
       int2_grad1_u12_ao = 0.d0
       do m = 1, 3
         !call dgemm( "T", "N", ao_num*ao_num, n_points_final_grid, n_points_extra_final_grid, +1.d0         &
@@ -108,7 +132,7 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
                   , 0.d0, int2_grad1_u12_ao(1,1,1,m), ao_num*ao_num) 
       enddo
 
-      !! these dgemm are equivalen to
+      !! these dgemm are equivalent to
       !!$OMP PARALLEL                                                           &
       !!$OMP DEFAULT (NONE)                                                     &
       !!$OMP PRIVATE (j, i, ipoint, jpoint, w)                                  &
@@ -132,16 +156,15 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_ao, (ao_num, ao_num, n_points_f
       !enddo
       !!$OMP END DO
       !!$OMP END PARALLEL
-    endif
 
       deallocate(tmp)
-  
     else
 
       print *, ' j1b_type = ', j1b_type, 'not implemented yet'
       stop
 
     endif
+  endif
 
   if(write_tc_integ.and.mpi_master) then
     open(unit=11, form="unformatted", file=trim(ezfio_filename)//'/work/int2_grad1_u12_ao', action="write")
@@ -176,20 +199,43 @@ BEGIN_PROVIDER [double precision, int2_grad1_u12_square_ao, (ao_num, ao_num, n_p
 
   PROVIDE j1b_type
 
-  if(j1b_type .eq. 3) then
+  if(j1b_type .eq. 0) then
 
-    PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+    PROVIDE int2_grad1u2_grad2u2
 
     int2_grad1_u12_square_ao = 0.d0
     !$OMP PARALLEL               &
     !$OMP DEFAULT (NONE)         &
     !$OMP PRIVATE (i, j, ipoint) &
-    !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, u12_grad1_u12_j1b_grad1_j1b, grad12_j12)
+    !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, int2_grad1u2_grad2u2)
     !$OMP DO SCHEDULE (static)
     do ipoint = 1, n_points_final_grid
       do j = 1, ao_num
         do i = 1, ao_num
-          int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+          int2_grad1_u12_square_ao(i,j,ipoint) = int2_grad1u2_grad2u2(i,j,ipoint)
+        enddo
+      enddo
+    enddo
+    !$OMP END DO
+    !$OMP END PARALLEL
+
+  elseif((j1b_type .eq. 3) .or. (j1b_type .eq. 4)) then
+
+    ! the term u12_grad1_u12_j1b_grad1_j1b is added directly for performance
+    !PROVIDE u12sq_j1bsq u12_grad1_u12_j1b_grad1_j1b grad12_j12
+    PROVIDE u12sq_j1bsq grad12_j12
+
+    int2_grad1_u12_square_ao = 0.d0
+    !$OMP PARALLEL               &
+    !$OMP DEFAULT (NONE)         &
+    !$OMP PRIVATE (i, j, ipoint) &
+    !$OMP SHARED (int2_grad1_u12_square_ao, ao_num, n_points_final_grid, u12sq_j1bsq, grad12_j12)
+    !$OMP DO SCHEDULE (static)
+    do ipoint = 1, n_points_final_grid
+      do j = 1, ao_num
+        do i = 1, ao_num
+          !int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + u12_grad1_u12_j1b_grad1_j1b(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
+          int2_grad1_u12_square_ao(i,j,ipoint) = u12sq_j1bsq(i,j,ipoint) + 0.5d0 * grad12_j12(i,j,ipoint)
         enddo
       enddo
     enddo

src/non_h_ints_mu/test_non_h_ints.irp.f

@@ -1,15 +1,25 @@
 program test_non_h
+
   implicit none
+
   my_grid_becke  = .True.
   my_n_pt_r_grid = 50
   my_n_pt_a_grid = 74
+  !my_n_pt_r_grid = 400
+  !my_n_pt_a_grid = 974
+
 !  my_n_pt_r_grid = 10 ! small grid for quick debug
 !  my_n_pt_a_grid = 26 ! small grid for quick debug
   touch my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
-!call routine_grad_squared
+
- call routine_fit
+  !call routine_grad_squared
+  !call routine_fit
+  
+  call test_ipp()
 end
 
+! ---
+
 subroutine routine_lapl_grad
  implicit none
  integer :: i,j,k,l
@@ -100,3 +110,445 @@ subroutine routine_fit
  enddo
 
 end
+
+
+subroutine test_ipp()
+
+  implicit none
+  integer                       :: i, j, k, l, ipoint
+  double precision              :: accu, norm, diff, old, new, eps, int_num
+  double precision              :: weight1, ao_i_r, ao_k_r
+  double precision, allocatable :: b_mat(:,:,:), I1(:,:,:,:), I2(:,:,:,:)
+
+  eps = 1d-7
+
+  allocate(b_mat(n_points_final_grid,ao_num,ao_num))
+  b_mat = 0.d0
+
+  ! ---
+
+  ! first way
+
+  allocate(I1(ao_num,ao_num,ao_num,ao_num))
+  I1 = 0.d0
+
+  PROVIDE u12_grad1_u12_j1b_grad1_j1b
+
+  !$OMP PARALLEL               &
+  !$OMP DEFAULT (NONE)         &
+  !$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
+    do k = 1, ao_num
+      do ipoint = 1, n_points_final_grid
+        b_mat(ipoint,k,i) = final_weight_at_r_vector(ipoint) * aos_in_r_array_transp(ipoint,i) * aos_in_r_array_transp(ipoint,k)
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0                    &
+            , u12_grad1_u12_j1b_grad1_j1b(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
+            , 0.d0, I1, ao_num*ao_num)
+
+  ! ---
+
+  ! 2nd way
+
+  allocate(I2(ao_num,ao_num,ao_num,ao_num))
+  I2 = 0.d0
+
+  PROVIDE int2_u2_j1b2
+
+  b_mat = 0.d0
+  !$OMP PARALLEL                                                                                     &
+  !$OMP DEFAULT (NONE)                                                                               &
+  !$OMP PRIVATE (i, k, ipoint, weight1, ao_i_r, ao_k_r)                                              &
+  !$OMP SHARED (aos_in_r_array_transp, b_mat, ao_num, n_points_final_grid, final_weight_at_r_vector, &
+  !$OMP         v_1b_square_grad, v_1b_square_lapl, aos_grad_in_r_array_transp_bis)
+  !$OMP DO SCHEDULE (static)
+  do i = 1, ao_num
+    do k = 1, ao_num
+      do ipoint = 1, n_points_final_grid
+
+        weight1 = 0.25d0 * 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) = weight1 * ( ao_k_r * ao_i_r * v_1b_square_lapl(ipoint)                                                                                   &
+                          + (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)) * v_1b_square_grad(ipoint,1) &
+                          + (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)) * v_1b_square_grad(ipoint,2) &
+                          + (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)) * v_1b_square_grad(ipoint,3) )
+      enddo
+    enddo
+  enddo
+  !$OMP END DO
+  !$OMP END PARALLEL
+
+  call dgemm( "N", "N", ao_num*ao_num, ao_num*ao_num, n_points_final_grid, 1.d0     &
+            , int2_u2_j1b2(1,1,1), ao_num*ao_num, b_mat(1,1,1), n_points_final_grid &
+            , 0.d0, I2, ao_num*ao_num)
+ 
+  ! ---
+
+  deallocate(b_mat)
+
+  accu = 0.d0
+  norm = 0.d0
+  do i = 1, ao_num
+    do k = 1, ao_num
+      do l = 1, ao_num
+        do j = 1, ao_num
+
+          old = I1(j,l,k,i)
+          new = I2(j,l,k,i)
+
+          !print*, l, k, j, i
+          !print*, old, new
+          
+          diff = new - old
+          if(dabs(diff) .gt. eps) then
+            print*, ' problem on :', j, l, k, i
+            print*, ' diff      = ', diff
+            print*, ' old value = ', old
+            print*, ' new value = ', new
+            call I_grade_gradu_naive1(i, j, k, l, int_num)
+            print*, ' full num1 = ', int_num
+            call I_grade_gradu_naive2(i, j, k, l, int_num)
+            print*, ' full num2 = ', int_num
+            call I_grade_gradu_naive3(i, j, k, l, int_num)
+            print*, ' full num3 = ', int_num
+            call I_grade_gradu_naive4(i, j, k, l, int_num)
+            print*, ' full num4 = ', int_num
+            call I_grade_gradu_seminaive(i, j, k, l, int_num)
+            print*, ' semi num  = ', int_num
+          endif
+
+          accu += dabs(diff)
+          norm += dabs(old)
+        enddo
+      enddo
+    enddo
+  enddo
+
+  deallocate(I1, I2)
+
+  print*, ' accu = ', accu
+  print*, ' norm = ', norm
+
+  return
+end subroutine test_ipp
+
+! ---
+
+subroutine I_grade_gradu_naive1(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1_x, weight1_y, weight1_z
+  double precision              :: weight2_x, weight2_y, weight2_z
+  double precision              :: aor_i, aor_j, aor_k, aor_l
+  double precision              :: e1_val, e2_val, e1_der(3), u12_val, u12_der(3)
+  double precision, external    :: j1b_nucl, j12_mu
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    e1_val = j1b_nucl(r1)
+    call grad1_j1b_nucl(r1, e1_der)
+
+    weight1_x = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(1)
+    weight1_y = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(2)
+    weight1_z = aor_i * aor_k * e1_val * final_weight_at_r_vector(ipoint) * e1_der(3)
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val = j1b_nucl(r2)
+
+      u12_val = j12_mu(r1, r2)
+      call grad1_j12_mu(r1, r2, u12_der)
+
+      weight2_x = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(1)
+      weight2_y = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(2)
+      weight2_z = aor_j * aor_l * e2_val * e2_val * u12_val * final_weight_at_r_vector_extra(jpoint) * u12_der(3)
+
+      int = int - (weight1_x * weight2_x  + weight1_y * weight2_y + weight1_z * weight2_z)
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive1
+
+! ---
+
+subroutine I_grade_gradu_naive2(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1_x, weight1_y, weight1_z
+  double precision              :: weight2_x, weight2_y, weight2_z
+  double precision              :: aor_i, aor_j, aor_k, aor_l
+  double precision              :: e1_square_der(3), e2_val, u12_square_der(3)
+  double precision, external    :: j1b_nucl
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    call grad1_j1b_nucl_square_num(r1, e1_square_der)
+
+    weight1_x = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(1)
+    weight1_y = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(2)
+    weight1_z = aor_i * aor_k * final_weight_at_r_vector(ipoint) * e1_square_der(3)
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val = j1b_nucl(r2)
+      call grad1_j12_mu_square_num(r1, r2, u12_square_der)
+
+      weight2_x = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(1)
+      weight2_y = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(2)
+      weight2_z = aor_j * aor_l * e2_val * e2_val * final_weight_at_r_vector_extra(jpoint) * u12_square_der(3)
+
+      int = int - 0.25d0 * (weight1_x * weight2_x  + weight1_y * weight2_y + weight1_z * weight2_z)
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive2
+
+! ---
+
+subroutine I_grade_gradu_naive3(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1, weight2
+  double precision              :: aor_j, aor_l
+  double precision              :: grad(3), e2_val, u12_val
+  double precision, external    :: j1b_nucl, j12_mu
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    call grad1_aos_ik_grad1_esquare(i, k, r1, grad)
+
+    weight1 = final_weight_at_r_vector(ipoint) * (grad(1) + grad(2) + grad(3))
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val  = j1b_nucl(r2)
+      u12_val = j12_mu(r1, r2)
+
+      weight2 = aor_j * aor_l * e2_val * e2_val * u12_val * u12_val * final_weight_at_r_vector_extra(jpoint)
+
+      int = int + 0.25d0 * weight1 * weight2
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive3
+
+! ---
+
+subroutine I_grade_gradu_naive4(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint, jpoint
+  double precision              :: r1(3), r2(3)
+  double precision              :: weight1, weight2
+  double precision              :: aor_j, aor_l, aor_k, aor_i
+  double precision              :: grad(3), e2_val, u12_val
+  double precision, external    :: j1b_nucl, j12_mu
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    r1(1) = final_grid_points(1,ipoint)
+    r1(2) = final_grid_points(2,ipoint)
+    r1(3) = final_grid_points(3,ipoint)
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    weight1 = final_weight_at_r_vector(ipoint) * ( aor_k * aor_i * v_1b_square_lapl(ipoint)                                                          &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,1) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,2) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,3) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) )
+
+    do jpoint = 1, n_points_extra_final_grid ! r2 
+
+      r2(1) = final_grid_points_extra(1,jpoint)
+      r2(2) = final_grid_points_extra(2,jpoint)
+      r2(3) = final_grid_points_extra(3,jpoint)
+
+      aor_j = aos_in_r_array_extra_transp(jpoint,j)
+      aor_l = aos_in_r_array_extra_transp(jpoint,l)
+
+      e2_val  = j1b_nucl(r2)
+      u12_val = j12_mu(r1, r2)
+
+      weight2 = aor_j * aor_l * e2_val * e2_val * u12_val * u12_val * final_weight_at_r_vector_extra(jpoint)
+
+      int = int + 0.25d0 * weight1 * weight2
+    enddo
+  enddo
+
+  return
+end subroutine I_grade_gradu_naive4
+
+! ---
+
+subroutine I_grade_gradu_seminaive(i, j, k, l, int)
+
+  implicit none
+  integer,          intent(in)  :: i, j, k, l
+  double precision, intent(out) :: int
+  integer                       :: ipoint
+  double precision              :: r1(3)
+  double precision              :: weight1
+  double precision              :: aor_i, aor_k
+
+  int = 0.d0
+
+  do ipoint = 1, n_points_final_grid ! r1
+
+    aor_i = aos_in_r_array_transp(ipoint,i)
+    aor_k = aos_in_r_array_transp(ipoint,k)
+
+    weight1 = 0.25d0 * final_weight_at_r_vector(ipoint) * ( aor_k * aor_i * v_1b_square_lapl(ipoint)                                                 &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,1) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,1)) * v_1b_square_grad(ipoint,1) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,2) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,2)) * v_1b_square_grad(ipoint,2) &
+            + (aor_k * aos_grad_in_r_array_transp_bis(ipoint,i,3) + aor_i * aos_grad_in_r_array_transp_bis(ipoint,k,3)) * v_1b_square_grad(ipoint,3) )
+
+    int = int + weight1 * int2_u2_j1b2(j,l,ipoint)
+  enddo
+
+  return
+end subroutine I_grade_gradu_seminaive
+
+! ---
+
+subroutine aos_ik_grad1_esquare(i, k, r1, val)
+
+  implicit none
+  integer,          intent(in)  :: i, k
+  double precision, intent(in)  :: r1(3)
+  double precision, intent(out) :: val(3)
+  double precision              :: tmp
+  double precision              :: der(3), aos_array(ao_num), aos_grad_array(3,ao_num)
+
+  call give_all_aos_and_grad_at_r(r1, aos_array, aos_grad_array)
+  call grad1_j1b_nucl_square_num(r1, der)
+
+  tmp    = aos_array(i) * aos_array(k)
+  val(1) = tmp * der(1)
+  val(2) = tmp * der(2)
+  val(3) = tmp * der(3)
+
+  return
+end subroutine phi_ik_grad1_esquare
+
+! ---
+
+subroutine grad1_aos_ik_grad1_esquare(i, k, r1, grad)
+
+  implicit none
+  integer,          intent(in)  :: i, k
+  double precision, intent(in)  :: r1(3)
+  double precision, intent(out) :: grad(3)
+  double precision              :: r(3), eps, tmp_eps, val_p(3), val_m(3)
+
+  eps     = 1d-5
+  tmp_eps = 0.5d0 / eps
+
+  r(1:3) = r1(1:3)
+
+  r(1) = r(1) + eps
+  call aos_ik_grad1_esquare(i, k, r, val_p)
+  r(1) = r(1) - 2.d0 * eps
+  call aos_ik_grad1_esquare(i, k, r, val_m)
+  r(1) = r(1) + eps
+  grad(1) = tmp_eps * (val_p(1) - val_m(1))
+
+  r(2) = r(2) + eps
+  call aos_ik_grad1_esquare(i, k, r, val_p)
+  r(2) = r(2) - 2.d0 * eps
+  call aos_ik_grad1_esquare(i, k, r, val_m)
+  r(2) = r(2) + eps
+  grad(2) = tmp_eps * (val_p(2) - val_m(2))
+
+  r(3) = r(3) + eps
+  call aos_ik_grad1_esquare(i, k, r, val_p)
+  r(3) = r(3) - 2.d0 * eps
+  call aos_ik_grad1_esquare(i, k, r, val_m)
+  r(3) = r(3) + eps
+  grad(3) = tmp_eps * (val_p(3) - val_m(3))
+
+  return
+end subroutine grad1_aos_ik_grad1_esquare
+
+! ---
+
+
+
+
+
+
+

src/non_h_ints_mu/total_tc_int.irp.f

@@ -11,6 +11,13 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num,
   call wall_time(wall0)
   
   if(test_cycle_tc) then
+
+    PROVIDE j1b_type
+    if(j1b_type .ne. 3) then
+      print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type
+      stop
+    endif
+
     do j = 1, ao_num
       do l = 1, ao_num
         do i = 1, ao_num
@@ -20,7 +27,9 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num,
         enddo
       enddo
     enddo
+
   else
+
     do j = 1, ao_num
       do l = 1, ao_num
         do i = 1, ao_num
@@ -30,6 +39,7 @@ BEGIN_PROVIDER [double precision, ao_vartc_int_chemist, (ao_num, ao_num, ao_num,
         enddo
       enddo
     enddo
+
   endif
 
   call wall_time(wall1)
@@ -50,6 +60,12 @@ BEGIN_PROVIDER [double precision, ao_tc_int_chemist, (ao_num, ao_num, ao_num, ao
   
   if(test_cycle_tc) then
 
+    PROVIDE j1b_type
+    if(j1b_type .ne. 3) then
+      print*, ' TC integrals with cycle can not be used for j1b_type =', j1b_type
+      stop
+    endif
+
     ao_tc_int_chemist = ao_tc_int_chemist_test
 
   else

src/scf_utils/NEED

@@ -1,2 +1,3 @@
 mo_guess 
 bitmask 
+json

src/scf_utils/roothaan_hall_scf.irp.f

@@ -12,6 +12,7 @@ END_DOC
 
   integer                        :: iteration_SCF,dim_DIIS,index_dim_DIIS
 
+  logical                        :: converged
   integer                        :: i,j
   logical, external              :: qp_stop
   double precision, allocatable :: mo_coef_save(:,:)
@@ -50,10 +51,8 @@ END_DOC
 !
   PROVIDE FPS_SPF_matrix_AO Fock_matrix_AO 
 
-  do while ( &
+  converged = .False.
-    ( (max_error_DIIS > threshold_DIIS_nonzero) .or. &
+  do while ( .not.converged .and. (iteration_SCF < n_it_SCF_max) )
-      (dabs(Delta_energy_SCF) > thresh_SCF) &
-    ) .and. (iteration_SCF < n_it_SCF_max) )
 
 ! Increment cycle number
 
@@ -144,17 +143,49 @@ END_DOC
     SOFT_TOUCH level_shift
     energy_SCF_previous = energy_SCF
 
+    converged = ( (max_error_DIIS <= threshold_DIIS_nonzero) .and. &
+                  (dabs(Delta_energy_SCF) <= thresh_SCF) )
+
 !   Print results at the end of each iteration
 
     write(6,'(I4, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, F16.10, 1X, I3)')  &
       iteration_SCF, energy_SCF, Delta_energy_SCF, max_error_DIIS, level_shift, dim_DIIS
 
+!   Write data in JSON file
+
+    call lock_io
+    if (iteration_SCF == 1) then
+      write(json_unit, json_dict_uopen_fmt)
+    else
+      write(json_unit, json_dict_close_uopen_fmt)
+    endif
+    write(json_unit, json_int_fmt) 'iteration', iteration_SCF
+    write(json_unit, json_real_fmt) 'energy', energy_SCF
+    write(json_unit, json_real_fmt) 'delta_energy_SCF', Delta_energy_SCF
+    write(json_unit, json_real_fmt) 'max_error_DIIS', max_error_DIIS
+    write(json_unit, json_real_fmt) 'level_shift', level_shift
+    write(json_unit, json_int_fmt) 'dim_DIIS', dim_DIIS
+    call unlock_io
+
     if (Delta_energy_SCF < 0.d0) then
       call save_mos
+      write(json_unit, json_true_fmt) 'saved'
+    else
+      write(json_unit, json_false_fmt) 'saved'
     endif
+
+    call lock_io
+    if (converged) then
+      write(json_unit, json_true_fmtx) 'converged'
+    else
+      write(json_unit, json_false_fmtx) 'converged'
+    endif
+    call unlock_io
+
     if (qp_stop()) exit
 
   enddo
+  write(json_unit, json_dict_close_fmtx)
 
  if (iteration_SCF < n_it_SCF_max) then
    mo_label = 'Canonical'
@@ -166,6 +197,10 @@ END_DOC
   write(6,'(A4, 1X, A16, 1X, A16, 1X, A16, 1X, A16)')  &
     '====','================','================','================','================'
   write(6,*)
+  if (converged) then
+    write(6,*) 'SCF converged'
+  endif
+
 
   if(.not.frozen_orb_scf)then
    call mo_as_eigvectors_of_mo_matrix(Fock_matrix_mo,size(Fock_matrix_mo,1), &

src/tc_bi_ortho/print_tc_energy.irp.f

@@ -9,6 +9,10 @@ program print_tc_energy
   my_n_pt_a_grid = 50
   read_wf = .True.
   touch read_wf
+
+  PROVIDE j1b_type
+  print*, 'j1b_type = ', j1b_type
+
   touch  my_grid_becke my_n_pt_r_grid my_n_pt_a_grid
  call write_tc_energy
 end

src/tc_bi_ortho/print_tc_wf.irp.f

@@ -26,7 +26,8 @@ subroutine write_l_r_wf
  integer :: i
  print*,'Writing the left-right wf'
  do i = 1, N_det
-  write(i_unit_output,*)i,psi_l_coef_sorted_bi_ortho_left(i),psi_r_coef_sorted_bi_ortho_right(i)
+  write(i_unit_output,*)i, psi_l_coef_sorted_bi_ortho_left(i)/psi_l_coef_sorted_bi_ortho_left(1) &
+                         , psi_r_coef_sorted_bi_ortho_right(i)/psi_r_coef_sorted_bi_ortho_right(1)
  enddo
 
 

src/tc_keywords/EZFIO.cfg

@@ -124,6 +124,12 @@ doc: type of 1-body Jastrow
 interface: ezfio, provider, ocaml
 default: 0
 
+[mu_r_ct]
+type: double precision
+doc: a parameter used to define mu(r)
+interface: ezfio, provider, ocaml
+default: 6.203504908994001e-1
+
 [thr_degen_tc]
 type: Threshold
 doc: Threshold to determine if two orbitals are degenerate in TCSCF in order to avoid random quasi orthogonality between the right- and left-eigenvector for the same eigenvalue
@@ -170,7 +176,7 @@ default: 1.e-7
 type: logical
 doc: If |true|, the integrals of the three-body jastrow are computed with cycles
 interface: ezfio,provider,ocaml
-default: True
+default: False
 
 [thresh_biorthog_diag]
 type: Threshold

src/tc_scf/minimize_tc_angles.irp.f

@@ -8,6 +8,10 @@ program print_angles
   my_n_pt_a_grid = 50
   touch my_n_pt_r_grid my_n_pt_a_grid
 !  call sort_by_tc_fock
+
+  ! TODO
+  ! check if rotations of orbitals affect the TC energy
+  ! and refuse the step
   call minimize_tc_orb_angles
 end
 

src/tc_scf/tc_scf.irp.f

@@ -53,7 +53,7 @@ program tc_scf
       stop
     endif
 
-    call minimize_tc_orb_angles()
+    !call minimize_tc_orb_angles()
 
   endif
 

src/tc_scf/tc_scf_energy.irp.f

@@ -1,5 +1,5 @@
 
- BEGIN_PROVIDER [ double precision, TC_HF_energy]
+ BEGIN_PROVIDER [ double precision, TC_HF_energy      ]
 &BEGIN_PROVIDER [ double precision, TC_HF_one_e_energy]
 &BEGIN_PROVIDER [ double precision, TC_HF_two_e_energy]
 

src/tools/print_e_conv.irp.f

@@ -1,80 +0,0 @@
-program print_e_conv
- implicit none
- BEGIN_DOC
-! program that prints in a human readable format the convergence of the CIPSI algorithm.
-!
-! for all istate, this program produces
-!
-! * a file "EZFIO.istate.conv" containing the variational and var+PT2 energies as a function of N_det
-!
-! * for istate > 1, a file EZFIO.istate.delta_e.conv containing the energy difference (both var and var+PT2) with the ground state as a function of N_det
- END_DOC
-
-  provide ezfio_filename
-  call routine_e_conv
- end
-
-subroutine routine_e_conv
- implicit none
- BEGIN_DOC
-! routine called by :c:func:`print_e_conv`
- END_DOC
- integer :: N_iter_tmp
- integer :: i,istate
- character*(128) :: output
- integer :: i_unit_output,getUnitAndOpen
- character*(128) :: filename
-
- integer, allocatable :: n_det_tmp(:)
- call ezfio_get_iterations_N_iter(N_iter_tmp)
- print*,'N_iter_tmp = ',N_iter_tmp
- double precision, allocatable :: e(:,:),pt2(:,:)
- allocate(e(N_states, 100),pt2(N_states, 100),n_det_tmp(100))
- call ezfio_get_iterations_energy_iterations(e)
- call ezfio_get_iterations_pt2_iterations(pt2)
- call ezfio_get_iterations_n_det_iterations(n_det_tmp)
-
-
-  do istate = 1, N_states
-   if (istate.lt.10)then
-    write (filename, "(I1)")istate
-   else
-    write (filename, "(I2)")istate
-   endif
-   print*,filename
-   output=trim(ezfio_filename)//'.'//trim(filename)//'.conv'
-   output=trim(output)
-   print*,'output = ',trim(output)
-   i_unit_output = getUnitAndOpen(output,'w')
-   write(i_unit_output,*)'# N_det    E_var        E_var + PT2'
-   do i = 1, N_iter_tmp
-    write(i_unit_output,'(I9,X,3(F16.10,X))')n_det_tmp(i),e(istate,i),e(istate,i) + pt2(istate,i)
-   enddo
-  enddo
-
-  if(N_states.gt.1)then
-   double precision, allocatable :: deltae(:,:),deltae_pt2(:,:)
-   allocate(deltae(N_states,100),deltae_pt2(N_states,100))
-   do i = 1, N_iter_tmp
-    do istate = 1, N_states
-     deltae(istate,i) = e(istate,i) - e(1,i)
-     deltae_pt2(istate,i) = e(istate,i) + pt2(istate,i) - (e(1,i) + pt2(1,i))
-    enddo
-   enddo
-   do istate = 2, N_states
-     if (istate.lt.10)then
-      write (filename, "(I1)")istate
-     else
-      write (filename, "(I2)")istate
-     endif
-     output=trim(ezfio_filename)//'.'//trim(filename)//'.delta_e.conv'
-     print*,'output = ',trim(output)
-     i_unit_output = getUnitAndOpen(output,'w')
-     write(i_unit_output,*)'# N_det    Delta E_var   Delta (E_var + PT2)'
-     do i = 1, N_iter_tmp
-      write(i_unit_output,'(I9,X,100(F16.10,X))')n_det_tmp(i),deltae(istate,i),deltae_pt2(istate,i)
-     enddo
-   enddo
-  endif
-
-end

src/two_body_rdm/io_two_rdm.irp.f

@@ -9,7 +9,9 @@ subroutine write_array_two_rdm(n_orb,nstates,array_tmp,name_file)
  PROVIDE ezfio_filename
  output=trim(ezfio_filename)//'/work/'//trim(name_file)
  i_unit_output = getUnitAndOpen(output,'W')
+ call lock_io()
  write(i_unit_output)array_tmp
+ call unlock_io()
  close(unit=i_unit_output)
 end
 
@@ -23,7 +25,9 @@ subroutine read_array_two_rdm(n_orb,nstates,array_tmp,name_file)
  PROVIDE ezfio_filename
  output=trim(ezfio_filename)//'/work/'//trim(name_file)
  i_unit_output = getUnitAndOpen(output,'R')
+ call lock_io()
  read(i_unit_output)array_tmp
+ call unlock_io()
  close(unit=i_unit_output)
 end
 

src/utils/format_w_error.irp.f

@@ -14,7 +14,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err
 
   ! out
   ! | format_value  | character | string FX.Y for the format |
-  ! | str_error     | character | string of the error |
+  ! | str_error     | character | string of the error        |
 
   ! internal
   ! | str_size      | character | size in string format                             |
@@ -33,6 +33,7 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err
   character(len=20)              :: str_size, str_nb_digits, str_exp
   integer                        :: nb_digits
 
+  call lock_io()
   ! max_nb_digit: Y max
   ! size_nb = Size of the double: X (FX.Y)
   write(str_size,'(I3)') size_nb
@@ -67,5 +68,6 @@ subroutine format_w_error(value,error,size_nb,max_nb_digits,format_value,str_err
 
   ! FX.Y just for the value
   format_value = 'F'//trim(adjustl(str_size))//'.'//trim(adjustl(str_nb_digits))
+  call unlock_io()
 
 end

src/utils/linear_algebra.irp.f

@@ -1854,7 +1854,7 @@ do k = 1, N
 end do
 ! TODO: It should be possible to use only one vector of size (1:rank) as a buffer
 ! to do the swapping in-place
-U = 0.00D+0
+U(:,:) = 0.00D+0
 do k = 1, N
   l = piv(k)
   U(l, :) = A(1:rank, k)

src/utils/memory.irp.f

@@ -8,7 +8,9 @@ BEGIN_PROVIDER [ integer, qp_max_mem ]
  qp_max_mem = 2000
  call getenv('QP_MAXMEM',env)
  if (trim(env) /= '') then
+    call lock_io()
     read(env,*) qp_max_mem
+    call unlock_io()
  endif
  call write_int(6,qp_max_mem,'Target maximum memory (GB)')
 
@@ -25,7 +27,7 @@ subroutine resident_memory(value)
   character*(32) :: key
   double precision, intent(out) :: value
 
-  call omp_set_lock(file_lock)
+  call lock_io()
   call usleep(10)
 
   value = 0.d0
@@ -40,7 +42,7 @@ subroutine resident_memory(value)
   20 continue
   close(iunit)
   value = value / (1024.d0*1024.d0)
-  call omp_unset_lock(file_lock)
+  call unlock_io()
 end function
 
 subroutine total_memory(value)
@@ -53,6 +55,7 @@ subroutine total_memory(value)
   character*(32) :: key
   double precision, intent(out) :: value
 
+  call lock_io()
   iunit = getUnitAndOpen('/proc/self/status','r')
   do
     read(iunit,*,err=10,end=20) key, value
@@ -64,6 +67,7 @@ subroutine total_memory(value)
   20 continue
   close(iunit)
   value = value / (1024.d0*1024.d0)
+  call unlock_io()
 end function
 
 double precision function memory_of_double(n)