added cosgtos_ao_int

src/cosgtos_ao_int/EZFIO.cfg

@@ -0,0 +1,19 @@
+[ao_expoim_cosgtos]
+type: double precision
+doc: imag part for Exponents for each primitive of each cosGTOs |AO|
+size: (ao_basis.ao_num,ao_basis.ao_prim_num_max)
+interface: ezfio, provider
+
+[use_cosgtos]
+type: logical
+doc: If true, use cosgtos for AO integrals
+interface: ezfio,provider,ocaml
+default: False
+
+[ao_integrals_threshold]
+type: Threshold
+doc: If | (pq|rs) | < `ao_integrals_threshold` then (pq|rs) is zero
+interface: ezfio,provider,ocaml
+default: 1.e-15
+ezfio_name: threshold_ao
+

src/cosgtos_ao_int/NEED

@@ -0,0 +1 @@
+ao_basis

src/cosgtos_ao_int/README.rst

@@ -0,0 +1,4 @@
+==============
+cosgtos_ao_int
+==============
+

src/cosgtos_ao_int/aos_cosgtos.irp.f

@@ -0,0 +1,210 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
+
+  implicit none
+  integer :: i, j
+
+  do j = 1, ao_num
+    do i = 1, ao_prim_num_max
+      ao_coef_norm_ord_transp_cosgtos(i,j) = ao_coef_norm_ord_cosgtos(j,i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ complex*16, ao_expo_ord_transp_cosgtos, (ao_prim_num_max, ao_num) ]
+
+  implicit none
+  integer :: i, j
+
+  do j = 1, ao_num
+    do i = 1, ao_prim_num_max
+      ao_expo_ord_transp_cosgtos(i,j) = ao_expo_ord_cosgtos(j,i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ double precision, ao_coef_norm_cosgtos, (ao_num, ao_prim_num_max) ]
+
+  implicit none
+
+  integer          :: i, j, powA(3), nz
+  double precision :: norm
+  complex*16       :: overlap_x, overlap_y, overlap_z, C_A(3)
+  complex*16       :: integ1, integ2, expo
+
+  nz = 100
+
+  C_A(1) = (0.d0, 0.d0)
+  C_A(2) = (0.d0, 0.d0)
+  C_A(3) = (0.d0, 0.d0)
+
+  ao_coef_norm_cosgtos = 0.d0
+
+  do i = 1, ao_num
+
+    powA(1) = ao_power(i,1)
+    powA(2) = ao_power(i,2)
+    powA(3) = ao_power(i,3)
+
+    ! Normalization of the primitives
+    if(primitives_normalized) then
+
+      do j = 1, ao_prim_num(i)
+
+        expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expoim_cosgtos(i,j)
+
+        call overlap_cgaussian_xyz(C_A, C_A,        expo, expo, powA, powA, overlap_x, overlap_y, overlap_z, integ1, nz)
+        call overlap_cgaussian_xyz(C_A, C_A, conjg(expo), expo, powA, powA, overlap_x, overlap_y, overlap_z, integ2, nz)
+
+        norm = 2.d0 * real( integ1 + integ2 )
+
+        ao_coef_norm_cosgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
+      enddo
+
+    else
+
+      do j = 1, ao_prim_num(i)
+        ao_coef_norm_cosgtos(i,j) = ao_coef(i,j)
+      enddo
+
+    endif
+
+  enddo
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, ao_coef_norm_ord_cosgtos, (ao_num, ao_prim_num_max) ]
+&BEGIN_PROVIDER [ complex*16      , ao_expo_ord_cosgtos,      (ao_num, ao_prim_num_max) ]
+
+  implicit none
+  integer          :: i, j
+  integer          :: iorder(ao_prim_num_max)
+  double precision :: d(ao_prim_num_max,3)
+
+  d = 0.d0
+
+  do i = 1, ao_num
+
+    do j = 1, ao_prim_num(i)
+      iorder(j) = j
+      d(j,1) = ao_expo(i,j)
+      d(j,2) = ao_coef_norm_cosgtos(i,j)
+      d(j,3) = ao_expoim_cosgtos(i,j)
+    enddo
+
+    call dsort     (d(1,1), iorder, ao_prim_num(i))
+    call dset_order(d(1,2), iorder, ao_prim_num(i))
+    call dset_order(d(1,3), iorder, ao_prim_num(i))
+
+    do j = 1, ao_prim_num(i)
+      ao_expo_ord_cosgtos     (i,j) = d(j,1) + (0.d0, 1.d0) * d(j,3)
+      ao_coef_norm_ord_cosgtos(i,j) = d(j,2)
+    enddo
+
+  enddo
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos,   (ao_num, ao_num) ]
+&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_x, (ao_num, ao_num) ]
+&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_y, (ao_num, ao_num) ]
+&BEGIN_PROVIDER [ double precision, ao_overlap_cosgtos_z, (ao_num, ao_num) ]
+
+  implicit none
+  integer          :: i, j, n, l, dim1, power_A(3), power_B(3)
+  double precision :: c, overlap, overlap_x, overlap_y, overlap_z
+  complex*16       :: alpha, beta, A_center(3), B_center(3)
+  complex*16       :: overlap1, overlap_x1, overlap_y1, overlap_z1
+  complex*16       :: overlap2, overlap_x2, overlap_y2, overlap_z2
+
+  ao_overlap_cosgtos   = 0.d0
+  ao_overlap_cosgtos_x = 0.d0
+  ao_overlap_cosgtos_y = 0.d0
+  ao_overlap_cosgtos_z = 0.d0
+
+  dim1 = 100
+
+ !$OMP PARALLEL DO SCHEDULE(GUIDED)                                                                 &
+ !$OMP DEFAULT(NONE)                                                                                &
+ !$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, n, l, c                    &
+ !$OMP        , overlap_x , overlap_y , overlap_z , overlap                                         &
+ !$OMP        , overlap_x1, overlap_y1, overlap_z1, overlap1                                        &
+ !$OMP        , overlap_x2, overlap_y2, overlap_z2, overlap2 )                                      &
+ !$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1                             &
+ !$OMP       , ao_overlap_cosgtos_x, ao_overlap_cosgtos_y, ao_overlap_cosgtos_z, ao_overlap_cosgtos &
+ !$OMP       , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos )
+
+  do j = 1, ao_num
+
+    A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0) 
+    A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
+    A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
+    power_A(1)  = ao_power(j,1)
+    power_A(2)  = ao_power(j,2)
+    power_A(3)  = ao_power(j,3)
+
+    do i = 1, ao_num
+
+      B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0) 
+      B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
+      B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
+      power_B(1)  = ao_power(i,1)
+      power_B(2)  = ao_power(i,2)
+      power_B(3)  = ao_power(i,3)
+
+      do n = 1, ao_prim_num(j)
+        alpha = ao_expo_ord_transp_cosgtos(n,j)
+
+        do l = 1, ao_prim_num(i)
+          c    = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
+          beta = ao_expo_ord_transp_cosgtos(l,i)
+
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
+                                    , overlap_x1, overlap_y1, overlap_z1, overlap1, dim1 )
+
+          call overlap_cgaussian_xyz( A_center, B_center, conjg(alpha), beta, power_A, power_B &
+                                    , overlap_x2, overlap_y2, overlap_z2, overlap2, dim1       )
+
+          overlap_x = 2.d0 * real( overlap_x1 + overlap_x2 )
+          overlap_y = 2.d0 * real( overlap_y1 + overlap_y2 )
+          overlap_z = 2.d0 * real( overlap_z1 + overlap_z2 )
+          overlap   = 2.d0 * real( overlap1   + overlap2   )
+
+          ao_overlap_cosgtos(i,j) = ao_overlap_cosgtos(i,j) + c * overlap
+
+          if( isnan(ao_overlap_cosgtos(i,j)) ) then
+            print*,'i, j', i, j
+            print*,'l, n', l, n
+            print*,'c, overlap', c, overlap
+            print*, overlap_x, overlap_y, overlap_z
+            stop
+          endif
+
+          ao_overlap_cosgtos_x(i,j) = ao_overlap_cosgtos_x(i,j) + c * overlap_x
+          ao_overlap_cosgtos_y(i,j) = ao_overlap_cosgtos_y(i,j) + c * overlap_y
+          ao_overlap_cosgtos_z(i,j) = ao_overlap_cosgtos_z(i,j) + c * overlap_z
+
+        enddo
+      enddo
+    enddo
+  enddo
+ !$OMP END PARALLEL DO
+
+END_PROVIDER
+
+! ---
+
+
+

src/cosgtos_ao_int/cosgtos_ao_int.irp.f

@@ -0,0 +1,7 @@
+program cosgtos_ao_int
+  implicit none
+  BEGIN_DOC
+! TODO : Put the documentation of the program here
+  END_DOC
+  print *, 'Hello world'
+end

src/cosgtos_ao_int/expoim_opt.py

@@ -0,0 +1,172 @@
+#!/usr/bin/env python
+
+import sys
+import os
+import subprocess
+from datetime import datetime
+import time
+import numpy as np
+from modif_powell_imp import my_fmin_powell
+
+QP_PATH=os.environ["QP_ROOT"]
+sys.path.insert(0, QP_PATH+"external/ezfio/Python")
+from ezfio import ezfio
+
+
+
+
+
+#------------------------------------------------------------------------------
+#
+def get_expoim():
+
+    expo_im = np.array(ezfio.get_cosgtos_ao_int_ao_expoim_cosgtos()).T
+    #print(expo_im.shape)
+
+    x = []
+    for i in range(ao_num):
+        for j in range(ao_prim_num[i]):
+            x.append(expo_im[i,j])
+
+    return x
+
+# ---
+
+def set_expoim(x):
+
+    expo_im = np.zeros((ao_num, ao_prim_num_max))
+
+    ii = 0
+    for i in range(ao_num):
+        for j in range(ao_prim_num[i]):
+            expo_im[i,j] = x[ii]
+            ii = ii + 1
+
+    ezfio.set_cosgtos_ao_int_ao_expoim_cosgtos(expo_im.T)
+#
+#------------------------------------------------------------------------------
+
+
+#------------------------------------------------------------------------------
+#
+def save_res(results, file_output):
+
+    lines = results.splitlines()
+    with open(file_output, "w") as f:
+        for line in lines:
+            f.write(f"{line}\n")
+
+#
+#------------------------------------------------------------------------------
+
+
+
+#------------------------------------------------------------------------------
+#
+def get_scfenergy(results):
+
+    scf_energy = 0.0
+
+    lines = results.splitlines()
+    for line in lines:
+        if("SCF energy" in line):
+            scf_energy = float(line.split()[-1])
+
+    return scf_energy
+#
+#------------------------------------------------------------------------------
+
+
+
+
+#------------------------------------------------------------------------------
+#
+def run_scf():
+
+    return subprocess.check_output( ['qp_run', 'scf', EZFIO_file]
+                                  , encoding = "utf-8" )
+#
+#------------------------------------------------------------------------------
+
+
+
+#------------------------------------------------------------------------------
+#
+def f_scf(x):
+
+    global i_call
+    i_call += 1
+
+    #print(x)
+
+    # set expo
+    set_expoim(x)
+
+    # run scf
+    results = run_scf()
+    #save_res(results, "scf_"+str(i_call))
+    
+
+    # get scf_energy
+    scf_energy = get_scfenergy(results)
+    print( scf_energy )
+    sys.stdout.flush()
+
+    return scf_energy
+#
+#------------------------------------------------------------------------------
+
+
+
+
+if __name__ == '__main__':
+
+    t0 = time.time()
+
+    EZFIO_file = sys.argv[1]
+    ezfio.set_file(EZFIO_file)
+    print(" Today's date:", datetime.now() )
+    print(" EZFIO file = {}".format(EZFIO_file))
+
+
+    ao_num = ezfio.get_ao_basis_ao_num()
+    print(f" ao_num = {ao_num}")
+
+    ao_prim_num = ezfio.get_ao_basis_ao_prim_num()
+
+    ao_prim_num_max = np.amax(ao_prim_num)
+    print(f" ao_prim_num_max = {ao_prim_num_max}")
+
+    ezfio.set_ao_basis_ao_prim_num_max(ao_prim_num_max)
+
+    x = get_expoim()
+
+    n_par = len(x)
+    print(' nb of parameters = {}'.format(n_par))
+
+    sys.stdout.flush()
+
+    #x    = (np.random.rand(n_par) - 0.5) * 1.0
+    x     = [ (+0.00) for _ in range(n_par)]
+
+    x_min = [ (-10.0) for _ in range(n_par)]
+    x_max = [ (+10.0) for _ in range(n_par)]
+
+    i_call = 0
+    memo_val = {'fmin': 100.}
+
+    opt = my_fmin_powell( f_scf
+                        , x, x_min, x_max
+                        #, xtol        = 1e-1
+                        #, ftol        = 1e-1
+                        , maxfev      = 1e8
+                        , full_output = 1
+                        , verbose     = 1 )
+
+
+    print(" x = " + str(opt))
+
+    print(" end after {:.3f} minutes".format((time.time()-t0)/60.) )
+
+    # !!!
+# !!!

src/cosgtos_ao_int/gauss_legendre.irp.f

@@ -0,0 +1,57 @@
+ BEGIN_PROVIDER [ double precision, gauleg_t2, (n_pt_max_integrals,n_pt_max_integrals/2) ]
+&BEGIN_PROVIDER [ double precision, gauleg_w, (n_pt_max_integrals,n_pt_max_integrals/2) ]
+   implicit none
+   BEGIN_DOC
+   ! t_w(i,1,k) = w(i)
+   ! t_w(i,2,k) = t(i)
+   END_DOC
+   integer                        :: i,j,l
+   l=0
+   do i = 2,n_pt_max_integrals,2
+     l = l+1
+     call gauleg(0.d0,1.d0,gauleg_t2(1,l),gauleg_w(1,l),i)
+     do j=1,i
+       gauleg_t2(j,l) *= gauleg_t2(j,l)
+     enddo
+   enddo
+
+END_PROVIDER
+
+subroutine gauleg(x1,x2,x,w,n)
+   implicit none
+   BEGIN_DOC
+   ! Gauss-Legendre
+   END_DOC
+   integer, intent(in)            :: n
+   double precision, intent(in)   :: x1, x2
+   double precision, intent (out) :: x(n),w(n)
+   double precision, parameter    :: eps=3.d-14
+
+   integer                        :: m,i,j
+   double precision               :: xm, xl, z, z1, p1, p2, p3, pp, dn
+   m=(n+1)/2
+   xm=0.5d0*(x2+x1)
+   xl=0.5d0*(x2-x1)
+   dn = dble(n)
+   do i=1,m
+     z=dcos(3.141592654d0*(dble(i)-.25d0)/(dble(n)+.5d0))
+     z1 = z+1.d0
+     do while (dabs(z-z1) > eps)
+       p1=1.d0
+       p2=0.d0
+       do j=1,n
+         p3=p2
+         p2=p1
+         p1=(dble(j+j-1)*z*p2-dble(j-1)*p3)/j
+       enddo
+       pp=dn*(z*p1-p2)/(z*z-1.d0)
+       z1=z
+       z=z1-p1/pp
+     end do
+     x(i)=xm-xl*z
+     x(n+1-i)=xm+xl*z
+     w(i)=(xl+xl)/((1.d0-z*z)*pp*pp)
+     w(n+1-i)=w(i)
+   enddo
+end
+

src/cosgtos_ao_int/one_e_Coul_integrals.irp.f

@@ -0,0 +1,535 @@
+
+! ---
+
+BEGIN_PROVIDER [ double precision, ao_integrals_n_e_cosgtos, (ao_num, ao_num)]
+
+  BEGIN_DOC
+  !
+  !  Nucleus-electron interaction, in the cosgtos |AO| basis set.
+  !
+  !  :math:`\langle \chi_i | -\sum_A \frac{1}{|r-R_A|} | \chi_j \rangle`
+  !
+  END_DOC
+
+  implicit none
+  integer          :: num_A, num_B, power_A(3), power_B(3)
+  integer          :: i, j, k, l, n_pt_in, m
+  double precision :: c, Z, A_center(3), B_center(3), C_center(3)
+  complex*16       :: alpha, beta, c1, c2
+
+  complex*16       :: NAI_pol_mult_cosgtos
+
+  ao_integrals_n_e_cosgtos = 0.d0
+
+ !$OMP PARALLEL                                                                           &
+ !$OMP DEFAULT (NONE)                                                                     &
+ !$OMP PRIVATE ( i, j, k, l, m, alpha, beta, A_center, B_center, C_center                 &
+ !$OMP         , power_A, power_B, num_A, num_B, Z, c, c1, c2, n_pt_in )                  &
+ !$OMP SHARED ( ao_num, ao_prim_num, ao_nucl, nucl_coord, ao_power, nucl_num, nucl_charge &
+ !$OMP        , ao_expo_ord_transp_cosgtos, ao_coef_norm_ord_transp_cosgtos               &
+ !$OMP        , n_pt_max_integrals, ao_integrals_n_e_cosgtos )
+
+ n_pt_in = n_pt_max_integrals
+
+ !$OMP DO SCHEDULE (dynamic)
+
+  do j = 1, ao_num
+    num_A         = ao_nucl(j)
+    power_A(1:3)  = ao_power(j,1:3)
+    A_center(1:3) = nucl_coord(num_A,1:3)
+
+    do i = 1, ao_num
+      num_B         = ao_nucl(i)
+      power_B(1:3)  = ao_power(i,1:3)
+      B_center(1:3) = nucl_coord(num_B,1:3)
+
+      do l = 1, ao_prim_num(j)
+        alpha = ao_expo_ord_transp_cosgtos(l,j)
+
+        do m = 1, ao_prim_num(i)
+          beta = ao_expo_ord_transp_cosgtos(m,i)
+
+          c = 0.d0
+          do k = 1, nucl_num
+
+            Z = nucl_charge(k)
+
+            C_center(1:3) = nucl_coord(k,1:3)
+
+            !print *, ' '
+            !print *, A_center, B_center, C_center, power_A, power_B
+            !print *, real(alpha), real(beta)
+
+            c1 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B &
+                                     , alpha, beta, C_center, n_pt_in )
+
+            !c2 = c1
+            c2 = NAI_pol_mult_cosgtos( A_center, B_center, power_A, power_B  &
+                                     , conjg(alpha), beta, C_center, n_pt_in )
+
+            !print *, ' c1 = ', real(c1)
+            !print *, ' c2 = ', real(c2)
+
+            c = c - Z * 2.d0 * real(c1 + c2)
+
+          enddo
+          ao_integrals_n_e_cosgtos(i,j) = ao_integrals_n_e_cosgtos(i,j) &
+                      + ao_coef_norm_ord_transp_cosgtos(l,j)            &
+                      * ao_coef_norm_ord_transp_cosgtos(m,i) * c
+        enddo
+      enddo
+    enddo
+  enddo
+
+ !$OMP END DO
+ !$OMP END PARALLEL
+
+END_PROVIDER
+
+! ---
+
+complex*16 function NAI_pol_mult_cosgtos(A_center, B_center, power_A, power_B, alpha, beta, C_center, n_pt_in)
+
+  BEGIN_DOC
+  !
+  ! Computes the electron-nucleus attraction with two primitves cosgtos.
+  !
+  ! :math:`\langle g_i | \frac{1}{|r-R_c|} | g_j \rangle`
+  !
+  END_DOC
+
+  implicit none
+  include 'utils/constants.include.F'
+
+  integer,          intent(in) :: n_pt_in, power_A(3), power_B(3)
+  double precision, intent(in) :: C_center(3), A_center(3), B_center(3)
+  complex*16,       intent(in) :: alpha, beta
+
+  integer                      :: i, n_pt, n_pt_out
+  double precision             :: dist, const_mod
+  complex*16                   :: p, p_inv, rho, dist_integral, const, const_factor, coeff, factor
+  complex*16                   :: accu, P_center(3)
+  complex*16                   :: d(0:n_pt_in)
+
+  complex*16                   :: V_n_e_cosgtos
+  complex*16                   :: crint
+
+  if ( (A_center(1)/=B_center(1)) .or. (A_center(2)/=B_center(2)) .or. (A_center(3)/=B_center(3)) .or. &
+       (A_center(1)/=C_center(1)) .or. (A_center(2)/=C_center(2)) .or. (A_center(3)/=C_center(3)) ) then
+
+    continue
+
+  else
+
+    NAI_pol_mult_cosgtos = V_n_e_cosgtos( power_A(1), power_A(2), power_A(3) &
+                                        , power_B(1), power_B(2), power_B(3) &
+                                        , alpha, beta )
+    return
+
+  endif
+
+  p     = alpha + beta
+  p_inv = (1.d0, 0.d0) / p
+  rho   = alpha * beta * p_inv
+
+  dist          = 0.d0
+  dist_integral = (0.d0, 0.d0)
+  do i = 1, 3
+    P_center(i)    = (alpha * A_center(i) + beta * B_center(i)) * p_inv
+    dist          += (A_center(i) - B_center(i)) * (A_center(i) - B_center(i))
+    dist_integral += (P_center(i) - C_center(i)) * (P_center(i) - C_center(i))
+  enddo
+
+  const_factor = dist * rho
+  const        = p * dist_integral
+
+  const_mod = dsqrt(real(const_factor)*real(const_factor) + aimag(const_factor)*aimag(const_factor))
+  if(const_mod > 80.d0) then
+    NAI_pol_mult_cosgtos = (0.d0, 0.d0)
+    return
+  endif
+
+  factor = zexp(-const_factor)
+  coeff  = dtwo_pi * factor * p_inv
+
+  do i = 0, n_pt_in
+    d(i) = (0.d0, 0.d0)
+  enddo
+
+  n_pt = 2 * ( (power_A(1) + power_B(1)) + (power_A(2) + power_B(2)) + (power_A(3) + power_B(3)) )
+  if(n_pt == 0) then
+    NAI_pol_mult_cosgtos = coeff * crint(0, const)
+    return
+  endif
+
+  call give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
+                                         , power_A, power_B, C_center, n_pt_in, d, n_pt_out)
+
+  if(n_pt_out < 0) then
+    NAI_pol_mult_cosgtos = (0.d0, 0.d0)
+    return
+  endif
+
+  accu = (0.d0, 0.d0)
+  do i = 0, n_pt_out, 2
+    accu += crint(shiftr(i, 1), const) * d(i)
+
+!    print *, shiftr(i, 1), real(const), real(d(i)), real(crint(shiftr(i, 1), const))
+  enddo
+  NAI_pol_mult_cosgtos = accu * coeff
+
+end function NAI_pol_mult_cosgtos
+
+! ---
+
+subroutine give_cpolynomial_mult_center_one_e( A_center, B_center, alpha, beta &
+                                             , power_A, power_B, C_center, n_pt_in, d, n_pt_out)
+
+  BEGIN_DOC
+  ! Returns the explicit polynomial in terms of the "t" variable of the following
+  !
+  ! $I_{x1}(a_x, d_x,p,q) \times I_{x1}(a_y, d_y,p,q) \times I_{x1}(a_z, d_z,p,q)$.
+  END_DOC
+
+  implicit none
+
+  integer,          intent(in) :: n_pt_in, power_A(3), power_B(3)
+  double precision, intent(in) :: A_center(3), B_center(3), C_center(3)
+  complex*16,       intent(in) :: alpha, beta
+  integer,         intent(out) :: n_pt_out
+  complex*16,      intent(out) :: d(0:n_pt_in)
+
+  integer                      :: a_x, b_x, a_y, b_y, a_z, b_z
+  integer                      :: n_pt1, n_pt2, n_pt3, dim, i, n_pt_tmp
+  complex*16                   :: p, P_center(3), rho, p_inv, p_inv_2
+  complex*16                   :: R1x(0:2), B01(0:2), R1xp(0:2),R2x(0:2)
+  complex*16                   :: d1(0:n_pt_in), d2(0:n_pt_in), d3(0:n_pt_in)
+
+  ASSERT (n_pt_in > 1)
+
+  p       = alpha + beta
+  p_inv   = (1.d0, 0.d0) / p
+  p_inv_2 = 0.5d0 * p_inv
+
+  do i = 1, 3
+    P_center(i) = (alpha * A_center(i) + beta * B_center(i)) * p_inv
+  enddo
+
+  do i = 0, n_pt_in
+    d(i)  = (0.d0, 0.d0)
+    d1(i) = (0.d0, 0.d0)
+    d2(i) = (0.d0, 0.d0)
+    d3(i) = (0.d0, 0.d0)
+  enddo
+
+  ! ---
+
+  n_pt1 = n_pt_in
+
+  R1x(0)  =  (P_center(1) - A_center(1))
+  R1x(1)  =  (0.d0, 0.d0)
+  R1x(2)  = -(P_center(1) - C_center(1))
+
+  R1xp(0) =  (P_center(1) - B_center(1))
+  R1xp(1) =  (0.d0, 0.d0)
+  R1xp(2) = -(P_center(1) - C_center(1))
+
+  R2x(0)  =  p_inv_2
+  R2x(1)  = (0.d0, 0.d0)
+  R2x(2)  = -p_inv_2
+
+  a_x = power_A(1)
+  b_x = power_B(1)
+  call I_x1_pol_mult_one_e_cosgtos(a_x, b_x, R1x, R1xp, R2x, d1, n_pt1, n_pt_in)
+
+  if(n_pt1 < 0) then
+    n_pt_out = -1
+    do i = 0, n_pt_in
+      d(i) = (0.d0, 0.d0)
+    enddo
+    return
+  endif
+
+  ! ---
+
+  n_pt2 = n_pt_in
+
+  R1x(0)  =  (P_center(2) - A_center(2))
+  R1x(1)  =  (0.d0, 0.d0)
+  R1x(2)  = -(P_center(2) - C_center(2))
+
+  R1xp(0) =  (P_center(2) - B_center(2))
+  R1xp(1) =  (0.d0, 0.d0)
+  R1xp(2) = -(P_center(2) - C_center(2))
+
+  a_y = power_A(2)
+  b_y = power_B(2)
+  call I_x1_pol_mult_one_e_cosgtos(a_y, b_y, R1x, R1xp, R2x, d2, n_pt2, n_pt_in)
+
+  if(n_pt2 < 0) then
+    n_pt_out = -1
+    do i = 0, n_pt_in
+      d(i) = (0.d0, 0.d0)
+    enddo
+    return
+  endif
+
+  ! ---
+
+  n_pt3 = n_pt_in
+
+  R1x(0)  =  (P_center(3) - A_center(3))
+  R1x(1)  =  (0.d0, 0.d0)
+  R1x(2)  = -(P_center(3) - C_center(3))
+
+  R1xp(0) =  (P_center(3) - B_center(3))
+  R1xp(1) =  (0.d0, 0.d0)
+  R1xp(2) = -(P_center(3) - C_center(3))
+
+  a_z = power_A(3)
+  b_z = power_B(3)
+  call I_x1_pol_mult_one_e_cosgtos(a_z, b_z, R1x, R1xp, R2x, d3, n_pt3, n_pt_in)
+
+  if(n_pt3 < 0) then
+    n_pt_out = -1
+    do i = 0, n_pt_in
+      d(i) = (0.d0, 0.d0)
+    enddo
+    return
+  endif
+
+  ! ---
+
+  n_pt_tmp = 0
+  call multiply_cpoly(d1, n_pt1, d2, n_pt2, d, n_pt_tmp)
+  do i = 0, n_pt_tmp
+    d1(i) = (0.d0, 0.d0)
+  enddo
+
+  n_pt_out = 0
+  call multiply_cpoly(d, n_pt_tmp, d3, n_pt3, d1, n_pt_out)
+  do i = 0, n_pt_out
+    d(i) = d1(i)
+  enddo
+
+end subroutine give_cpolynomial_mult_center_one_e
+
+! ---
+
+recursive subroutine I_x1_pol_mult_one_e_cosgtos(a, c, R1x, R1xp, R2x, d, nd, n_pt_in)
+
+  BEGIN_DOC
+  !  Recursive routine involved in the electron-nucleus potential
+  END_DOC
+
+  implicit none
+  include 'utils/constants.include.F'
+
+  integer,    intent(in)    :: a, c, n_pt_in
+  complex*16, intent(in)    :: R1x(0:2), R1xp(0:2), R2x(0:2)
+  integer,    intent(inout) :: nd
+  complex*16, intent(inout) :: d(0:n_pt_in)
+
+  integer                   :: nx, ix, dim, iy, ny
+  complex*16                :: X(0:max_dim)
+  complex*16                :: Y(0:max_dim)
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
+
+  dim = n_pt_in
+
+  if( (a==0) .and. (c==0)) then
+
+    nd   = 0
+    d(0) = (1.d0, 0.d0)
+    return
+
+  elseif( (c < 0) .or. (nd < 0) ) then
+
+    nd = -1
+    return
+
+  elseif((a == 0) .and. (c .ne. 0)) then
+
+    call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, n_pt_in)
+
+  elseif(a == 1) then
+
+    nx = nd
+    do ix = 0, n_pt_in
+      X(ix) = (0.d0, 0.d0)
+      Y(ix) = (0.d0, 0.d0)
+    enddo
+
+    call I_x2_pol_mult_one_e_cosgtos(c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
+
+    do ix = 0, nx
+      X(ix) *= dble(c)
+    enddo
+
+    call multiply_cpoly(X, nx, R2x, 2, d, nd)
+
+    ny = 0
+    call I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, Y, ny, n_pt_in)
+    call multiply_cpoly(Y, ny, R1x, 2, d, nd)
+
+  else
+
+    nx = 0
+    do ix = 0, n_pt_in
+      X(ix) = (0.d0, 0.d0)
+      Y(ix) = (0.d0, 0.d0)
+    enddo
+
+    call I_x1_pol_mult_one_e_cosgtos(a-2, c, R1x, R1xp, R2x, X, nx, n_pt_in)
+
+    do ix = 0, nx
+      X(ix) *= dble(a-1)
+    enddo
+    call multiply_cpoly(X, nx, R2x, 2, d, nd)
+
+    nx = nd
+    do ix = 0, n_pt_in
+      X(ix) = (0.d0, 0.d0)
+    enddo
+
+    call I_x1_pol_mult_one_e_cosgtos(a-1, c-1, R1x, R1xp, R2x, X, nx, n_pt_in)
+    do ix = 0, nx
+      X(ix) *= dble(c)
+    enddo
+
+    call multiply_cpoly(X, nx, R2x, 2, d, nd)
+
+    ny = 0
+    call I_x1_pol_mult_one_e_cosgtos(a-1, c, R1x, R1xp, R2x, Y, ny, n_pt_in)
+    call multiply_cpoly(Y, ny, R1x, 2, d, nd)
+
+  endif
+
+end subroutine I_x1_pol_mult_one_e_cosgtos
+
+! ---
+
+recursive subroutine I_x2_pol_mult_one_e_cosgtos(c, R1x, R1xp, R2x, d, nd, dim)
+
+  BEGIN_DOC
+  !  Recursive routine involved in the electron-nucleus potential
+  END_DOC
+
+  implicit none
+  include 'utils/constants.include.F'
+
+  integer,    intent(in)    :: dim, c
+  complex*16, intent(in)    :: R1x(0:2), R1xp(0:2), R2x(0:2)
+  integer,    intent(inout) :: nd
+  complex*16, intent(out)   :: d(0:max_dim)
+
+  integer                   :: i, nx, ix, ny
+  complex*16                :: X(0:max_dim), Y(0:max_dim)
+  !DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: X, Y
+
+  if(c == 0) then
+
+    nd = 0
+    d(0) = (1.d0, 0.d0)
+    return
+
+  elseif((nd < 0) .or. (c < 0)) then
+
+    nd = -1
+    return
+
+  else
+
+    nx = 0
+    do ix = 0, dim
+      X(ix) = (0.d0, 0.d0)
+      Y(ix) = (0.d0, 0.d0)
+    enddo
+
+    call I_x1_pol_mult_one_e_cosgtos(0, c-2, R1x, R1xp, R2x, X, nx, dim)
+
+    do ix = 0, nx
+      X(ix) *= dble(c-1)
+    enddo
+
+    call multiply_cpoly(X, nx, R2x, 2, d, nd)
+
+    ny = 0
+    do ix = 0, dim
+      Y(ix) = (0.d0, 0.d0)
+    enddo
+
+    call I_x1_pol_mult_one_e_cosgtos(0, c-1, R1x, R1xp, R2x, Y, ny, dim)
+
+    if(ny .ge. 0) then
+      call multiply_cpoly(Y, ny, R1xp, 2, d, nd)
+    endif
+
+  endif
+
+end subroutine I_x2_pol_mult_one_e_cosgtos
+
+! ---
+
+complex*16 function V_n_e_cosgtos(a_x, a_y, a_z, b_x, b_y, b_z, alpha, beta)
+
+  BEGIN_DOC
+  ! Primitve nuclear attraction between the two primitves centered on the same atom.
+  !
+  ! $p_1 = x^{a_x} y^{a_y} z^{a_z} \exp(-\alpha r^2)$
+  !
+  ! $p_2 = x^{b_x} y^{b_y} z^{b_z} \exp(-\beta  r^2)$
+  END_DOC
+
+  implicit none
+
+  integer,    intent(in) :: a_x, a_y, a_z, b_x, b_y, b_z
+  complex*16, intent(in) :: alpha, beta
+
+  double precision       :: V_phi, V_theta
+  complex*16             :: V_r_cosgtos
+
+  if( (iand(a_x + b_x, 1) == 1) .or. &
+      (iand(a_y + b_y, 1) == 1) .or. &
+      (iand(a_z + b_z, 1) == 1) ) then
+
+    V_n_e_cosgtos = (0.d0, 0.d0)
+
+  else
+
+    V_n_e_cosgtos = V_r_cosgtos(a_x + b_x + a_y + b_y + a_z + b_z + 1, alpha + beta) &
+                  * V_phi(a_x + b_x, a_y + b_y)                                      &
+                  * V_theta(a_z + b_z, a_x + b_x + a_y + b_y + 1)
+  endif
+
+end function V_n_e_cosgtos
+
+! ---
+
+complex*16 function V_r_cosgtos(n, alpha)
+
+  BEGIN_DOC
+  ! Computes the radial part of the nuclear attraction integral:
+  !
+  ! $\int_{0}^{\infty} r^n  \exp(-\alpha  r^2)  dr$
+  !
+  END_DOC
+
+  implicit none
+  include 'utils/constants.include.F'
+
+  integer   , intent(in) :: n
+  complex*16, intent(in) :: alpha
+
+  double precision       :: fact
+
+  if(iand(n, 1) .eq. 1) then
+    V_r_cosgtos = 0.5d0 * fact(shiftr(n, 1)) / (alpha**(shiftr(n, 1) + 1))
+  else
+    V_r_cosgtos = sqpi * fact(n) / fact(shiftr(n, 1)) * (0.5d0/zsqrt(alpha))**(n+1)
+  endif
+
+end function V_r_cosgtos
+
+! ---
+

src/cosgtos_ao_int/one_e_kin_integrals.irp.f

@@ -0,0 +1,223 @@
+
+! ---
+
+ BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_x, (ao_num, ao_num) ]
+&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_y, (ao_num, ao_num) ]
+&BEGIN_PROVIDER [ double precision, ao_deriv2_cosgtos_z, (ao_num, ao_num) ]
+
+  implicit none
+  integer          :: i, j, n, l, dim1, power_A(3), power_B(3)
+  double precision :: c, deriv_tmp
+  complex*16       :: alpha, beta, A_center(3), B_center(3)
+  complex*16       :: overlap_x, overlap_y, overlap_z, overlap
+  complex*16       :: overlap_x0_1, overlap_y0_1, overlap_z0_1
+  complex*16       :: overlap_x0_2, overlap_y0_2, overlap_z0_2 
+  complex*16       :: overlap_m2_1, overlap_p2_1
+  complex*16       :: overlap_m2_2, overlap_p2_2
+  complex*16       :: deriv_tmp_1, deriv_tmp_2
+
+
+  dim1 = 100
+
+  ! -- Dummy call to provide everything
+
+  A_center(:) = (0.0d0, 0.d0)
+  B_center(:) = (1.0d0, 0.d0)
+  alpha       = (1.0d0, 0.d0)
+  beta        = (0.1d0, 0.d0)
+  power_A     = 1
+  power_B     = 0
+  call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B       &
+                            , overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
+
+  ! ---
+
+ !$OMP PARALLEL DO SCHEDULE(GUIDED)                                                                  &
+ !$OMP DEFAULT(NONE)                                                                                 &
+ !$OMP PRIVATE( A_center, B_center, power_A, power_B, alpha, beta, i, j, l, n, c                     &
+ !$OMP        , deriv_tmp, deriv_tmp_1, deriv_tmp_2                                                  &
+ !$OMP        , overlap_x, overlap_y, overlap_z, overlap                                             &
+ !$OMP        , overlap_m2_1, overlap_p2_1, overlap_m2_2, overlap_p2_2                               &
+ !$OMP        , overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap_x0_2, overlap_y0_2, overlap_z0_2 ) &
+ !$OMP SHARED( nucl_coord, ao_power, ao_prim_num, ao_num, ao_nucl, dim1                              &
+ !$OMP       , ao_coef_norm_ord_transp_cosgtos, ao_expo_ord_transp_cosgtos                           & 
+ !$OMP       , ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z ) 
+
+  do j = 1, ao_num
+    A_center(1) = nucl_coord(ao_nucl(j),1) * (1.d0, 0.d0)
+    A_center(2) = nucl_coord(ao_nucl(j),2) * (1.d0, 0.d0)
+    A_center(3) = nucl_coord(ao_nucl(j),3) * (1.d0, 0.d0)
+    power_A(1)  = ao_power(j,1)
+    power_A(2)  = ao_power(j,2)
+    power_A(3)  = ao_power(j,3)
+
+    do i = 1, ao_num
+      B_center(1) = nucl_coord(ao_nucl(i),1) * (1.d0, 0.d0) 
+      B_center(2) = nucl_coord(ao_nucl(i),2) * (1.d0, 0.d0)
+      B_center(3) = nucl_coord(ao_nucl(i),3) * (1.d0, 0.d0)
+      power_B(1)  = ao_power(i,1)
+      power_B(2)  = ao_power(i,2)
+      power_B(3)  = ao_power(i,3)
+
+      ao_deriv2_cosgtos_x(i,j) = 0.d0
+      ao_deriv2_cosgtos_y(i,j) = 0.d0
+      ao_deriv2_cosgtos_z(i,j) = 0.d0
+
+      do n = 1, ao_prim_num(j)
+        alpha = ao_expo_ord_transp_cosgtos(n,j)
+
+        do l = 1, ao_prim_num(i)
+          c    = ao_coef_norm_ord_transp_cosgtos(n,j) * ao_coef_norm_ord_transp_cosgtos(l,i)
+          beta = ao_expo_ord_transp_cosgtos(l,i)
+
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B       &
+                                    , overlap_x0_1, overlap_y0_1, overlap_z0_1, overlap, dim1 )
+
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
+                                    , overlap_x0_2, overlap_y0_2, overlap_z0_2, overlap, dim1  )
+
+          ! ---
+
+          power_A(1) = power_A(1) - 2
+          if(power_A(1) > -1) then
+            call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
+                                      , overlap_m2_1, overlap_y, overlap_z, overlap, dim1 )
+
+            call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
+                                      , overlap_m2_2, overlap_y, overlap_z, overlap, dim1        )
+          else
+            overlap_m2_1 = (0.d0, 0.d0)
+            overlap_m2_2 = (0.d0, 0.d0)
+          endif
+
+          power_A(1) = power_A(1) + 4
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
+                                    , overlap_p2_1, overlap_y, overlap_z, overlap, dim1 )
+
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
+                                    , overlap_p2_2, overlap_y, overlap_z, overlap, dim1        )
+
+          power_A(1) = power_A(1) - 2
+
+          deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_1 &
+                        + power_A(1) * (power_A(1) - 1.d0) * overlap_m2_1           &
+                        + 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_y0_1 * overlap_z0_1
+
+          deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(1) + 1.d0) * overlap_x0_2 &
+                        + power_A(1) * (power_A(1) - 1.d0) * overlap_m2_2           &
+                        + 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_y0_2 * overlap_z0_2
+
+          deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
+
+          ao_deriv2_cosgtos_x(i,j) += c * deriv_tmp
+
+          ! ---
+
+          power_A(2) = power_A(2) - 2
+          if(power_A(2) > -1) then
+            call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
+                                      , overlap_x, overlap_m2_1, overlap_y, overlap, dim1 )
+
+            call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
+                                      , overlap_x, overlap_m2_2, overlap_y, overlap, dim1        )
+          else
+            overlap_m2_1 = (0.d0, 0.d0)
+            overlap_m2_2 = (0.d0, 0.d0)
+          endif
+
+          power_A(2) = power_A(2) + 4
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
+                                    , overlap_x, overlap_p2_1, overlap_y, overlap, dim1 )
+
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
+                                    , overlap_x, overlap_p2_2, overlap_y, overlap, dim1        )
+
+          power_A(2) = power_A(2) - 2
+
+          deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_1 &
+                        + power_A(2) * (power_A(2) - 1.d0) * overlap_m2_1           &
+                        + 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_z0_1
+
+          deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(2) + 1.d0) * overlap_y0_2 &
+                        + power_A(2) * (power_A(2) - 1.d0) * overlap_m2_2           &
+                        + 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_z0_2
+
+          deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
+
+          ao_deriv2_cosgtos_y(i,j) += c * deriv_tmp
+
+          ! ---
+
+          power_A(3) = power_A(3) - 2
+          if(power_A(3) > -1) then
+            call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
+                                      , overlap_x, overlap_y, overlap_m2_1, overlap, dim1 )
+
+            call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
+                                      , overlap_x, overlap_y, overlap_m2_2, overlap, dim1        )
+          else
+            overlap_m2_1 = (0.d0, 0.d0)
+            overlap_m2_2 = (0.d0, 0.d0)
+          endif
+
+          power_A(3) = power_A(3) + 4
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, beta, power_A, power_B &
+                                    , overlap_x, overlap_y, overlap_p2_1, overlap, dim1 )
+
+          call overlap_cgaussian_xyz( A_center, B_center, alpha, conjg(beta), power_A, power_B &
+                                    , overlap_x, overlap_y, overlap_p2_2, overlap, dim1        )
+
+          power_A(3) = power_A(3) - 2
+          
+          deriv_tmp_1 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_1 &
+                        + power_A(3) * (power_A(3) - 1.d0) * overlap_m2_1           &
+                        + 4.d0 * alpha * alpha * overlap_p2_1 ) * overlap_x0_1 * overlap_y0_1
+
+          deriv_tmp_2 = ( -2.d0 * alpha * (2.d0 * power_A(3) + 1.d0) * overlap_z0_2 &
+                        + power_A(3) * (power_A(3) - 1.d0) * overlap_m2_2           &
+                        + 4.d0 * alpha * alpha * overlap_p2_2 ) * overlap_x0_2 * overlap_y0_2
+
+          deriv_tmp = 2.d0 * real(deriv_tmp_1 + deriv_tmp_2)
+
+          ao_deriv2_cosgtos_z(i,j) += c * deriv_tmp
+
+          ! ---
+
+        enddo
+      enddo
+    enddo
+  enddo
+ !$OMP END PARALLEL DO
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, ao_kinetic_integrals_cosgtos, (ao_num, ao_num)]
+
+  BEGIN_DOC
+  ! 
+  ! Kinetic energy integrals in the cosgtos |AO| basis.
+  !
+  ! $\langle \chi_i |\hat{T}| \chi_j \rangle$
+  !
+  END_DOC
+
+  implicit none
+  integer :: i, j
+
+ !$OMP PARALLEL DO DEFAULT(NONE) &
+ !$OMP  PRIVATE(i, j)            &
+ !$OMP  SHARED(ao_num, ao_kinetic_integrals_cosgtos, ao_deriv2_cosgtos_x, ao_deriv2_cosgtos_y, ao_deriv2_cosgtos_z)
+  do j = 1, ao_num
+    do i = 1, ao_num
+      ao_kinetic_integrals_cosgtos(i,j) = -0.5d0 * ( ao_deriv2_cosgtos_x(i,j) &
+                                                   + ao_deriv2_cosgtos_y(i,j) &
+                                                   + ao_deriv2_cosgtos_z(i,j) )
+    enddo
+  enddo
+ !$OMP END PARALLEL DO
+
+END_PROVIDER
+
+! ---