starting to create cart

external/irpf90

@@ -1 +1 @@
-Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b
+Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6

plugins/local/mo_localization/localization_fboys_sub.irp.f

@@ -0,0 +1,193 @@
+
+! Criterion FB (old)
+
+! The criterion is just computed as
+
+! \begin{align*}
+! C = - \sum_i^{mo_{num}} (<i|x|i>^2 + <i|y|i>^2 + <i|z|i>^2)
+! \end{align*}
+
+! The minus sign is here in order to minimize this criterion
+
+! Output:
+! | criterion | double precision | criterion for the Foster-Boys localization |
+
+
+subroutine criterion_FB_old(criterion)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute the Foster-Boys localization criterion
+  END_DOC
+
+  double precision, intent(out) :: criterion
+  integer                       :: i
+
+  ! Criterion (= \sum_i <i|r|i>^2 )
+  criterion = 0d0
+  do i = 1, mo_num
+    criterion = criterion + mo_dipole_x(i,i)**2 + mo_dipole_y(i,i)**2 + mo_dipole_z(i,i)**2
+  enddo
+  criterion = - criterion
+
+end subroutine
+
+! Criterion FB
+
+subroutine criterion_FB(tmp_list_size, tmp_list, criterion)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute the Foster-Boys localization criterion
+  END_DOC
+
+  integer, intent(in)           :: tmp_list_size, tmp_list(tmp_list_size)
+  double precision, intent(out) :: criterion
+  integer                       :: i, tmp_i
+
+  ! Criterion (= - \sum_i <i|r|i>^2 )
+  criterion = 0d0
+  do tmp_i = 1, tmp_list_size
+    i = tmp_list(tmp_i)
+    criterion = criterion + mo_dipole_x(i,i)**2 + mo_dipole_y(i,i)**2 + mo_dipole_z(i,i)**2
+  enddo
+  criterion = - criterion
+
+end subroutine
+
+subroutine theta_FB(l, n, m_x, max_elem)
+
+  include 'pi.h'
+
+  BEGIN_DOC
+  ! Compute the angles to minimize the Foster-Boys criterion by using pairwise rotations of the MOs
+  ! Warning: you must give - the angles to build the rotation matrix...
+  END_DOC
+
+  implicit none
+
+  integer, intent(in)           :: n, l(n)
+  double precision, intent(out) :: m_x(n,n), max_elem
+
+  integer                       :: i,j, tmp_i, tmp_j
+  double precision, allocatable :: cos4theta(:,:), sin4theta(:,:)
+  double precision, allocatable :: A(:,:), B(:,:), beta(:,:), gamma(:,:)
+  integer                       :: idx_i,idx_j
+
+  allocate(cos4theta(n, n), sin4theta(n, n))
+  allocate(A(n,n), B(n,n), beta(n,n), gamma(n,n))
+
+  do tmp_j = 1, n
+    j = l(tmp_j)
+    do tmp_i = 1, n
+      i = l(tmp_i)
+      A(tmp_i,tmp_j) = mo_dipole_x(i,j)**2 - 0.25d0 * (mo_dipole_x(i,i) - mo_dipole_x(j,j))**2 &
+                     + mo_dipole_y(i,j)**2 - 0.25d0 * (mo_dipole_y(i,i) - mo_dipole_y(j,j))**2 &
+                     + mo_dipole_z(i,j)**2 - 0.25d0 * (mo_dipole_z(i,i) - mo_dipole_z(j,j))**2
+    enddo
+    A(j,j) = 0d0
+  enddo
+
+  do tmp_j = 1, n
+    j = l(tmp_j)
+    do tmp_i = 1, n
+      i = l(tmp_i)
+      B(tmp_i,tmp_j) = mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) &
+                     + mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) &
+                     + mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j))
+    enddo
+  enddo
+
+  !do tmp_j = 1, n
+  !  j = l(tmp_j)
+  !  do tmp_i = 1, n
+  !    i = l(tmp_i)
+  !    beta(tmp_i,tmp_j) =  (mo_dipole_x(i,i) - mo_dipole_x(j,j)) - 4d0 * mo_dipole_x(i,j)**2 &
+  !               + (mo_dipole_y(i,i) - mo_dipole_y(j,j)) - 4d0 * mo_dipole_y(i,j)**2 &
+  !               + (mo_dipole_z(i,i) - mo_dipole_z(j,j)) - 4d0 * mo_dipole_z(i,j)**2
+  !  enddo
+  !enddo
+
+  !do tmp_j = 1, n
+  !  j = l(tmp_j)
+  !  do tmp_i = 1, n
+  !    i = l(tmp_i)
+  !    gamma(tmp_i,tmp_j) = 4d0 * ( mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) &
+  !                       + mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) &
+  !                       + mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j)))
+  !  enddo
+  !enddo
+
+  !
+  !do j = 1, n
+  !  do i = 1, n
+  !    cos4theta(i,j) = - A(i,j) / dsqrt(A(i,j)**2 + B(i,j)**2)
+  !  enddo
+  !enddo
+
+  !do j = 1, n
+  !  do i = 1, n
+  !    sin4theta(i,j) = B(i,j) / dsqrt(A(i,j)**2 + B(i,j)**2)
+  !  enddo
+  !enddo
+
+  ! Theta
+  do j = 1, n
+    do i = 1, n
+      m_x(i,j) = 0.25d0 * atan2(B(i,j), -A(i,j))
+      !m_x(i,j) = 0.25d0 * atan2(sin4theta(i,j), cos4theta(i,j))
+    enddo
+  enddo
+
+  ! Enforce a perfect antisymmetry
+  do j = 1, n-1
+    do i = j+1, n
+      m_x(j,i) = - m_x(i,j)
+    enddo
+  enddo
+  do i = 1, n
+    m_x(i,i) = 0d0
+  enddo
+
+  ! Max
+  max_elem = 0d0
+  do j = 1, n-1
+    do i = j+1, n
+      if (dabs(m_x(i,j)) > dabs(max_elem)) then
+        max_elem = m_x(i,j)
+        !idx_i = i
+        !idx_j = j
+      endif
+    enddo
+  enddo
+
+  ! Debug
+  !print*,''
+  !print*,'sin/B'
+  !do i = 1, n
+  !  write(*,'(100F10.4)') sin4theta(i,:)
+  !  !B(i,:)
+  !enddo
+  !print*,'cos/A'
+  !do i = 1, n
+  !  write(*,'(100F10.4)') cos4theta(i,:)
+  !  !A(i,:)
+  !enddo
+  !print*,'X'
+  !!m_x = 0d0
+  !!m_x(idx_i,idx_j) = max_elem
+  !!m_x(idx_j,idx_i) = -max_elem
+  !do i = 1, n
+  !  write(*,'(100F10.4)') m_x(i,:)
+  !enddo
+  !print*,idx_i,idx_j,max_elem
+
+  max_elem = dabs(max_elem)
+  
+  deallocate(cos4theta, sin4theta)
+  deallocate(A,B,beta,gamma)
+  
+end
+

plugins/local/mo_localization/localization_pipek_sub.irp.f

@@ -0,0 +1,776 @@
+! Gradient v1
+subroutine grad_pipek(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute gradient for the Pipek-Mezey localization
+  END_DOC
+
+  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
+  double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
+  double precision, allocatable :: m_grad(:,:), tmp_int(:,:)
+  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho 
+
+  ! Allocation
+  allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size))
+
+  ! Initialization 
+  m_grad = 0d0
+
+  do a = 1, nucl_num ! loop over the nuclei
+     tmp_int = 0d0 ! Initialization for each nuclei
+
+     ! Loop over the MOs of the a given mo_class to compute <i|P_a|j>
+     do tmp_j = 1, tmp_list_size
+        j = tmp_list(tmp_j) 
+        do tmp_i = 1, tmp_list_size
+           i = tmp_list(tmp_i)
+           do rho = 1, ao_num ! loop over all the AOs
+              do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+                 mu = nucl_aos(a,b) ! AO centered on atom a 
+
+                 tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+                      + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
+
+              enddo
+           enddo
+        enddo
+     enddo
+
+     ! Gradient
+     do tmp_j = 1, tmp_list_size
+        do tmp_i = 1, tmp_list_size
+
+           m_grad(tmp_i,tmp_j) = m_grad(tmp_i,tmp_j) +  4d0 * tmp_int(tmp_i,tmp_j) * (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))
+
+        enddo
+     enddo
+
+  enddo
+
+  ! 2D -> 1D
+  do tmp_k = 1, tmp_n
+     call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
+     v_grad(tmp_k) = m_grad(tmp_i,tmp_j) 
+  enddo
+
+  ! Maximum element in the gradient
+  max_elem = 0d0
+  do tmp_k = 1, tmp_n
+     if (ABS(v_grad(tmp_k)) > max_elem) then
+        max_elem = ABS(v_grad(tmp_k))
+     endif
+  enddo
+
+  ! Norm of the gradient
+  norm_grad = 0d0
+  do tmp_k = 1, tmp_n
+     norm_grad = norm_grad + v_grad(tmp_k)**2
+  enddo
+  norm_grad = dsqrt(norm_grad)
+
+  print*, 'Maximal element in the gradient:', max_elem
+  print*, 'Norm of the gradient:', norm_grad
+
+  ! Deallocation
+  deallocate(m_grad,tmp_int)
+
+end subroutine grad_pipek
+
+! Gradient
+
+! The gradient is
+
+! \begin{align*}
+! \left. \frac{\partial \mathcal{P} (\theta)}{\partial \theta} \right|_{\theta=0}= \gamma^{PM}
+! \end{align*}
+! with
+! \begin{align*}
+! \gamma_{st}^{PM} = \sum_{A=1}^N <s|P_A|t> \left[ <s| P_A |s> - <t|P_A|t> \right]
+! \end{align*}
+
+! \begin{align*}
+! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right]
+! \end{align*}
+! $\sum_{\rho}$ -> sum over all the AOs
+! $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A
+! $c^t$ -> expansion coefficient of orbital |t>
+
+! Input: 
+! | tmp_n                   | integer          | Number of parameters in the MO subspace           |
+! | tmp_list_size           | integer          | Number of MOs in the mo_class we want to localize |
+! | tmp_list(tmp_list_size) | integer          | MOs in the mo_class                               |
+
+! Output:
+! | v_grad(tmp_n)           | double precision | Gradient in the subspace                          |
+! | max_elem                | double precision | Maximal element in the gradient                   |
+! | norm_grad               | double precision | Norm of the gradient                              |
+
+! Internal:
+! | m_grad(tmp_list_size,tmp_list_size)       | double precision | Gradient in a 2D array                                   |
+! | tmp_int(tmp_list_size,tmp_list_size)      |                  | Temporary array to store the integrals                   |
+! | tmp_accu(tmp_list_size,tmp_list_size)     |                  | Temporary array to store a matrix                        |
+! |                                           |                  | product and compute tmp_int                              |
+! | CS(tmp_list_size,ao_num)                  |                  | Array to store the result of mo_coef * ao_overlap        |
+! | tmp_mo_coef(ao_num,tmp_list_size)         |                  | Array to store just the useful MO coefficients           |
+! |                                           |                  | depending of the mo_class                                |
+! | tmp_mo_coef2(nucl_n_aos(a),tmp_list_size) |                  | Array to store just the useful MO coefficients           |
+! |                                           |                  | depending of the nuclei                                  |
+! | tmp_CS(tmp_list_size,nucl_n_aos(a))       |                  | Array to store just the useful mo_coef * ao_overlap      |
+! |                                           |                  | values depending of the nuclei                           |
+! | a                                         |                  | index to loop over the nuclei                            |
+! | b                                         |                  | index to loop over the AOs which belongs to the nuclei a |
+! | mu                                        |                  | index to refer to an AO which belongs to the nuclei a    |
+! | rho                                       |                  | index to loop over all the AOs                           |
+
+
+subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute gradient for the Pipek-Mezey localization
+  END_DOC
+  
+  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
+  double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
+  double precision, allocatable :: m_grad(:,:), tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
+  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho
+  double precision              :: t1,t2,t3
+
+  print*,''
+  print*,'---gradient_PM---'
+
+  call wall_time(t1)
+
+  ! Allocation
+  allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size),tmp_accu(tmp_list_size, tmp_list_size))
+  allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size))
+
+
+  ! submatrix of the mo_coef
+  do tmp_i = 1, tmp_list_size
+    i = tmp_list(tmp_i)
+    do j = 1, ao_num
+
+      tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
+ 
+    enddo
+  enddo
+
+  call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
+ 
+  m_grad = 0d0
+
+  do a = 1, nucl_num ! loop over the nuclei
+    tmp_int = 0d0
+
+    !do tmp_j = 1, tmp_list_size
+    !  do tmp_i = 1, tmp_list_size
+    !    do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+    !      mu = nucl_aos(a,b)
+
+    !      tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu))
+
+    !                             !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+    !                             !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
+
+    !    enddo
+    !  enddo
+    !enddo
+
+    allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
+
+    do tmp_i = 1, tmp_list_size
+      do b = 1, nucl_n_aos(a)
+        mu = nucl_aos(a,b)
+
+        tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i)
+
+      enddo
+    enddo
+    
+    do b = 1, nucl_n_aos(a)
+      mu = nucl_aos(a,b)
+      do tmp_i = 1, tmp_list_size
+
+        tmp_CS(tmp_i,b) = CS(tmp_i,mu)
+
+      enddo
+    enddo   
+
+    call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
+ 
+    do tmp_j = 1, tmp_list_size
+      do tmp_i = 1, tmp_list_size
+
+        tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))   
+
+      enddo
+    enddo
+
+    deallocate(tmp_mo_coef2,tmp_CS)
+
+    do tmp_j = 1, tmp_list_size
+      do tmp_i = 1, tmp_list_size
+
+        m_grad(tmp_i,tmp_j) = m_grad(tmp_i,tmp_j) +  4d0 * tmp_int(tmp_i,tmp_j) * (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))
+
+      enddo
+    enddo
+
+  enddo
+
+  ! 2D -> 1D
+  do tmp_k = 1, tmp_n
+    call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
+    v_grad(tmp_k) = m_grad(tmp_i,tmp_j) 
+  enddo
+
+  ! Maximum element in the gradient
+  max_elem = 0d0
+  do tmp_k = 1, tmp_n
+    if (ABS(v_grad(tmp_k)) > max_elem) then
+      max_elem = ABS(v_grad(tmp_k))
+    endif
+  enddo 
+
+  ! Norm of the gradient
+  norm_grad = 0d0
+  do tmp_k = 1, tmp_n
+    norm_grad = norm_grad + v_grad(tmp_k)**2
+  enddo
+  norm_grad = dsqrt(norm_grad)
+
+  print*, 'Maximal element in the gradient:', max_elem
+  print*, 'Norm of the gradient:', norm_grad
+
+  ! Deallocation
+  deallocate(m_grad,tmp_int,CS,tmp_mo_coef)
+
+  call wall_time(t2)
+  t3 = t2 - t1
+  print*,'Time in gradient_PM:', t3
+
+  print*,'---End gradient_PM---'
+
+end
+
+! Hessian v1
+
+subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H)
+
+  implicit none
+
+  BEGIN_DOC
+  ! Compute diagonal hessian for the Pipek-Mezey localization
+  END_DOC
+
+  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
+  double precision, intent(out) :: H(tmp_n)
+  double precision, allocatable :: beta(:,:),tmp_int(:,:)
+  integer                       :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu
+  double precision              :: max_elem
+    
+  ! Allocation
+  allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size))
+
+  beta = 0d0
+
+  do a = 1, nucl_num
+    tmp_int = 0d0
+
+    do tmp_j = 1, tmp_list_size
+      j = tmp_list(tmp_j)
+      do tmp_i = 1, tmp_list_size
+        i = tmp_list(tmp_i)
+        do rho = 1, ao_num
+          do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+            mu = nucl_aos(a,b)
+
+            tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+                                   + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
+
+          enddo
+        enddo  
+      enddo
+    enddo
+
+    ! Calculation
+    do tmp_j = 1, tmp_list_size
+      do tmp_i = 1, tmp_list_size
+
+        beta(tmp_i,tmp_j) = beta(tmp_i, tmp_j) +  (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))**2 - 4d0 * tmp_int(tmp_i,tmp_j)**2
+
+      enddo
+    enddo
+  
+  enddo
+
+  H = 0d0
+  do tmp_k = 1, tmp_n
+    call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
+    H(tmp_k) = 4d0 * beta(tmp_i, tmp_j)
+  enddo
+  
+  ! Deallocation
+  deallocate(beta,tmp_int)
+
+end
+
+! Hessian
+
+! The hessian is
+! \begin{align*}
+! \left. \frac{\partial^2 \mathcal{P} (\theta)}{\partial \theta^2}\right|_{\theta=0} = 4 \beta^{PM}
+! \end{align*}
+! \begin{align*}
+! \beta_{st}^{PM} = \sum_{A=1}^N \left( <s|P_A|t>^2 - \frac{1}{4} \left[<s|P_A|s> - <t|P_A|t> \right]^2 \right)
+! \end{align*}
+
+! with
+! \begin{align*}
+! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right]
+! \end{align*}
+! $\sum_{\rho}$ -> sum over all the AOs
+! $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A
+! $c^t$ -> expansion coefficient of orbital |t>
+
+
+subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
+
+  implicit none
+  
+  BEGIN_DOC
+  ! Compute diagonal hessian for the Pipek-Mezey localization
+  END_DOC
+
+  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
+  double precision, intent(out) :: H(tmp_n)
+  double precision, allocatable :: beta(:,:),tmp_int(:,:),CS(:,:),tmp_mo_coef(:,:),tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
+  integer                       :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu
+  double precision              :: max_elem, t1,t2,t3
+    
+  print*,''
+  print*,'---hessian_PM---'
+
+  call wall_time(t1)
+
+  ! Allocation
+  allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size),tmp_accu(tmp_list_size,tmp_list_size))
+  allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size))
+
+  beta = 0d0
+
+  do tmp_i = 1, tmp_list_size
+    i = tmp_list(tmp_i)
+    do j = 1, ao_num
+
+      tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
+ 
+    enddo
+  enddo
+
+  call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
+ 
+  do a = 1, nucl_num ! loop over the nuclei
+    tmp_int = 0d0
+
+    !do tmp_j = 1, tmp_list_size
+    !  do tmp_i = 1, tmp_list_size
+    !    do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+    !      mu = nucl_aos(a,b)
+
+    !      tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu))
+
+    !                             !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+    !                             !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
+
+    !    enddo
+    !  enddo
+    !enddo
+ 
+    allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
+
+    do tmp_i = 1, tmp_list_size
+      do b = 1, nucl_n_aos(a)
+        mu = nucl_aos(a,b)
+
+        tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i)
+
+      enddo
+    enddo
+    
+    do b = 1, nucl_n_aos(a)
+      mu = nucl_aos(a,b)
+      do tmp_i = 1, tmp_list_size
+
+        tmp_CS(tmp_i,b) = CS(tmp_i,mu)
+
+      enddo
+    enddo   
+
+    call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
+ 
+    do tmp_j = 1, tmp_list_size
+      do tmp_i = 1, tmp_list_size
+
+        tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))   
+
+      enddo
+    enddo
+
+    deallocate(tmp_mo_coef2,tmp_CS)
+
+    ! Calculation
+    do tmp_j = 1, tmp_list_size
+      do tmp_i = 1, tmp_list_size
+
+        beta(tmp_i,tmp_j) = beta(tmp_i, tmp_j) +  (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))**2 - 4d0 * tmp_int(tmp_i,tmp_j)**2
+
+      enddo
+    enddo
+  
+  enddo
+
+  H = 0d0
+  do tmp_k = 1, tmp_n
+    call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
+    H(tmp_k) = 4d0 * beta(tmp_i, tmp_j)
+  enddo
+  
+  ! Deallocation
+  deallocate(beta,tmp_int)
+
+  call wall_time(t2)
+  t3 = t2 - t1
+  print*,'Time in hessian_PM:', t3
+
+  print*,'---End hessian_PM---'
+
+end
+
+! Criterion PM (old)
+
+subroutine compute_crit_pipek(criterion)
+
+  implicit none
+  
+  BEGIN_DOC
+  ! Compute the Pipek-Mezey localization criterion
+  END_DOC
+
+  double precision, intent(out) :: criterion
+  double precision, allocatable :: tmp_int(:,:)
+  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho 
+
+  ! Allocation
+  allocate(tmp_int(mo_num, mo_num))
+ 
+  criterion = 0d0
+
+  do a = 1, nucl_num ! loop over the nuclei
+    tmp_int = 0d0
+
+    do i = 1, mo_num
+      do rho = 1, ao_num ! loop over all the AOs
+        do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+          mu = nucl_aos(a,b)
+
+          tmp_int(i,i) = tmp_int(i,i) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,i) &
+                                 + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,i))
+
+        enddo
+      enddo  
+    enddo
+
+    do i = 1, mo_num 
+      criterion = criterion + tmp_int(i,i)**2
+    enddo
+
+  enddo
+  
+  criterion = - criterion 
+
+  deallocate(tmp_int)
+
+end
+
+! Criterion PM
+
+! The criterion is computed as
+! \begin{align*}
+! \mathcal{P} = \sum_{i=1}^n \sum_{A=1}^N \left[ <i|P_A|i> \right]^2
+! \end{align*}
+! with 
+! \begin{align*}
+! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right]
+! \end{align*}
+
+
+subroutine criterion_PM(tmp_list_size,tmp_list,criterion)
+
+  implicit none
+  
+  BEGIN_DOC
+  ! Compute the Pipek-Mezey localization criterion
+  END_DOC
+
+  integer, intent(in)           :: tmp_list_size, tmp_list(tmp_list_size)
+  double precision, intent(out) :: criterion
+  double precision, allocatable :: tmp_int(:,:),CS(:,:)
+  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho
+  
+  print*,''
+  print*,'---criterion_PM---'
+  
+  ! Allocation
+  allocate(tmp_int(tmp_list_size, tmp_list_size),CS(mo_num,ao_num))
+  
+  ! Initialization
+  criterion = 0d0
+
+  call dgemm('T','N',mo_num,ao_num,ao_num,1d0,mo_coef,size(mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
+ 
+  do a = 1, nucl_num ! loop over the nuclei
+    tmp_int = 0d0
+
+      do tmp_i = 1, tmp_list_size
+        i = tmp_list(tmp_i)
+        do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+          mu = nucl_aos(a,b)
+          
+          tmp_int(tmp_i,tmp_i) = tmp_int(tmp_i,tmp_i) + 0.5d0 * (CS(i,mu) * mo_coef(mu,i) + mo_coef(mu,i) * CS(i,mu))
+
+                                 !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+                                 !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
+
+      enddo
+    enddo
+
+    do tmp_i = 1, tmp_list_size
+      criterion = criterion + tmp_int(tmp_i,tmp_i)**2
+    enddo
+
+  enddo
+  
+  criterion = - criterion 
+
+  deallocate(tmp_int,CS)
+
+  print*,'---End criterion_PM---'
+
+end
+
+! Criterion PM v3
+
+subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
+
+  implicit none
+  
+  BEGIN_DOC
+  ! Compute the Pipek-Mezey localization criterion
+  END_DOC
+  
+  integer, intent(in)           :: tmp_list_size, tmp_list(tmp_list_size)
+  double precision, intent(out) :: criterion
+  double precision, allocatable :: tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
+  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho,nu,c
+  double precision              :: t1,t2,t3
+
+  print*,''
+  print*,'---criterion_PM_v3---'
+  
+  call wall_time(t1)
+
+  ! Allocation
+  allocate(tmp_int(tmp_list_size, tmp_list_size),tmp_accu(tmp_list_size, tmp_list_size))
+  allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size))
+
+  criterion = 0d0
+
+  ! submatrix of the mo_coef
+  do tmp_i = 1, tmp_list_size
+    i = tmp_list(tmp_i)
+    do j = 1, ao_num
+
+      tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
+ 
+    enddo
+  enddo
+
+  ! ao_overlap(ao_num,ao_num)
+  ! mo_coef(ao_num,mo_num)
+  call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
+ 
+  do a = 1, nucl_num ! loop over the nuclei
+  
+    do j = 1, tmp_list_size
+      do i = 1, tmp_list_size
+        tmp_int(i,j) = 0d0
+      enddo
+    enddo
+
+    !do tmp_j = 1, tmp_list_size
+    !  do tmp_i = 1, tmp_list_size
+    !    do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+    !      mu = nucl_aos(a,b)
+
+    !      tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu))
+
+    !                             !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+    !                             !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
+
+    !    enddo
+    !  enddo
+    !enddo
+
+    allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
+ 
+    do tmp_i = 1, tmp_list_size
+      do b = 1, nucl_n_aos(a)
+        mu = nucl_aos(a,b)
+
+        tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i)
+
+      enddo
+    enddo
+    
+    do b = 1, nucl_n_aos(a)
+      mu = nucl_aos(a,b)
+      do tmp_i = 1, tmp_list_size
+
+         tmp_CS(tmp_i,b) = CS(tmp_i,mu)
+
+      enddo
+    enddo   
+
+    call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
+ 
+    ! Integrals
+    do tmp_j = 1, tmp_list_size
+      do tmp_i = 1, tmp_list_size
+
+        tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))   
+
+      enddo
+    enddo
+
+    deallocate(tmp_mo_coef2,tmp_CS)
+
+    ! Criterion
+    do tmp_i = 1, tmp_list_size
+      criterion = criterion + tmp_int(tmp_i,tmp_i)**2
+    enddo
+
+  enddo
+
+  criterion = - criterion 
+
+  deallocate(tmp_int,CS,tmp_accu,tmp_mo_coef)
+
+  call wall_time(t2)
+  t3 = t2 - t1
+  print*,'Time in criterion_PM_v3:', t3
+
+  print*,'---End criterion_PM_v3---'
+
+end
+
+subroutine theta_PM(l, n, m_x, max_elem)
+  
+  include 'pi.h'
+
+  BEGIN_DOC
+  ! Compute the angles to minimize the Pipek-Mezey criterion by using pairwise rotations of the MOs
+  ! Warning: you must give - the angles to build the rotation matrix...
+  END_DOC
+
+  implicit none
+
+  integer, intent(in)           :: n, l(n)
+  double precision, intent(out) :: m_x(n,n), max_elem
+
+  integer                       :: a,b,i,j,tmp_i,tmp_j,rho,mu,nu,idx_i,idx_j
+  double precision, allocatable :: Aij(:,:), Bij(:,:), Pa(:,:)
+
+  allocate(Aij(n,n), Bij(n,n), Pa(n,n))
+
+  do a = 1, nucl_num ! loop over the nuclei
+    Pa = 0d0 ! Initialization for each nuclei
+
+    ! Loop over the MOs of the a given mo_class to compute <i|P_a|j>
+    do tmp_j = 1, n
+      j = l(tmp_j) 
+      do tmp_i = 1, n
+         i = l(tmp_i)
+        do rho = 1, ao_num ! loop over all the AOs
+          do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
+            mu = nucl_aos(a,b) ! AO centered on atom a 
+
+            Pa(tmp_i,tmp_j) = Pa(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
+                                   + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
+
+          enddo
+        enddo  
+      enddo
+    enddo
+
+    ! A
+    do j = 1, n
+      do i = 1, n
+        Aij(i,j) = Aij(i,j) + Pa(i,j)**2 - 0.25d0 * (Pa(i,i) - Pa(j,j))**2
+      enddo
+    enddo
+    
+    ! B
+    do j = 1, n
+      do i = 1, n
+        Bij(i,j) = Bij(i,j) + Pa(i,j) * (Pa(i,i) - Pa(j,j))
+      enddo
+    enddo
+
+  enddo
+
+  ! Theta
+  do j = 1, n
+    do i = 1, n
+      m_x(i,j) = 0.25d0 * atan2(Bij(i,j), -Aij(i,j))
+    enddo
+  enddo
+
+  ! Enforce a perfect antisymmetry
+  do j = 1, n-1
+    do i = j+1, n
+      m_x(j,i) = - m_x(i,j)
+    enddo
+  enddo
+  do i = 1, n
+    m_x(i,i) = 0d0
+  enddo
+
+  ! Max
+  max_elem = 0d0
+  do j = 1, n-1
+    do i = j+1, n
+      if (dabs(m_x(i,j)) > dabs(max_elem)) then
+        max_elem = m_x(i,j)
+        idx_i = i
+        idx_j = j
+      endif
+    enddo
+  enddo
+
+  ! Debug
+  !do i = 1, n
+  !  write(*,'(100F10.4)') m_x(i,:)
+  !enddo
+  !print*,'Max',idx_i,idx_j,max_elem
+
+  max_elem = dabs(max_elem)
+
+  deallocate(Aij,Bij,Pa)
+
+end
+

plugins/local/mo_localization/localization_sub.irp.f

@@ -486,975 +486,6 @@ subroutine hessian_FB_omp(tmp_n, tmp_list_size, tmp_list, H)
 
 end subroutine
 
-! Gradient v1
-
-subroutine grad_pipek(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
-
-  implicit none
-
-  BEGIN_DOC
-  ! Compute gradient for the Pipek-Mezey localization
-  END_DOC
-
-  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
-  double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
-  double precision, allocatable :: m_grad(:,:), tmp_int(:,:)
-  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho 
-
-  ! Allocation
-  allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size))
-
-  ! Initialization 
-  m_grad = 0d0
-
-  do a = 1, nucl_num ! loop over the nuclei
-     tmp_int = 0d0 ! Initialization for each nuclei
-
-     ! Loop over the MOs of the a given mo_class to compute <i|P_a|j>
-     do tmp_j = 1, tmp_list_size
-        j = tmp_list(tmp_j) 
-        do tmp_i = 1, tmp_list_size
-           i = tmp_list(tmp_i)
-           do rho = 1, ao_num ! loop over all the AOs
-              do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-                 mu = nucl_aos(a,b) ! AO centered on atom a 
-
-                 tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
-                      + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
-
-              enddo
-           enddo
-        enddo
-     enddo
-
-     ! Gradient
-     do tmp_j = 1, tmp_list_size
-        do tmp_i = 1, tmp_list_size
-
-           m_grad(tmp_i,tmp_j) = m_grad(tmp_i,tmp_j) +  4d0 * tmp_int(tmp_i,tmp_j) * (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))
-
-        enddo
-     enddo
-
-  enddo
-
-  ! 2D -> 1D
-  do tmp_k = 1, tmp_n
-     call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
-     v_grad(tmp_k) = m_grad(tmp_i,tmp_j) 
-  enddo
-
-  ! Maximum element in the gradient
-  max_elem = 0d0
-  do tmp_k = 1, tmp_n
-     if (ABS(v_grad(tmp_k)) > max_elem) then
-        max_elem = ABS(v_grad(tmp_k))
-     endif
-  enddo
-
-  ! Norm of the gradient
-  norm_grad = 0d0
-  do tmp_k = 1, tmp_n
-     norm_grad = norm_grad + v_grad(tmp_k)**2
-  enddo
-  norm_grad = dsqrt(norm_grad)
-
-  print*, 'Maximal element in the gradient:', max_elem
-  print*, 'Norm of the gradient:', norm_grad
-
-  ! Deallocation
-  deallocate(m_grad,tmp_int)
-
-end subroutine grad_pipek
-
-! Gradient
-
-! The gradient is
-
-! \begin{align*}
-! \left. \frac{\partial \mathcal{P} (\theta)}{\partial \theta} \right|_{\theta=0}= \gamma^{PM}
-! \end{align*}
-! with
-! \begin{align*}
-! \gamma_{st}^{PM} = \sum_{A=1}^N <s|P_A|t> \left[ <s| P_A |s> - <t|P_A|t> \right]
-! \end{align*}
-
-! \begin{align*}
-! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right]
-! \end{align*}
-! $\sum_{\rho}$ -> sum over all the AOs
-! $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A
-! $c^t$ -> expansion coefficient of orbital |t>
-
-! Input: 
-! | tmp_n                   | integer          | Number of parameters in the MO subspace           |
-! | tmp_list_size           | integer          | Number of MOs in the mo_class we want to localize |
-! | tmp_list(tmp_list_size) | integer          | MOs in the mo_class                               |
-
-! Output:
-! | v_grad(tmp_n)           | double precision | Gradient in the subspace                          |
-! | max_elem                | double precision | Maximal element in the gradient                   |
-! | norm_grad               | double precision | Norm of the gradient                              |
-
-! Internal:
-! | m_grad(tmp_list_size,tmp_list_size)       | double precision | Gradient in a 2D array                                   |
-! | tmp_int(tmp_list_size,tmp_list_size)      |                  | Temporary array to store the integrals                   |
-! | tmp_accu(tmp_list_size,tmp_list_size)     |                  | Temporary array to store a matrix                        |
-! |                                           |                  | product and compute tmp_int                              |
-! | CS(tmp_list_size,ao_num)                  |                  | Array to store the result of mo_coef * ao_overlap        |
-! | tmp_mo_coef(ao_num,tmp_list_size)         |                  | Array to store just the useful MO coefficients           |
-! |                                           |                  | depending of the mo_class                                |
-! | tmp_mo_coef2(nucl_n_aos(a),tmp_list_size) |                  | Array to store just the useful MO coefficients           |
-! |                                           |                  | depending of the nuclei                                  |
-! | tmp_CS(tmp_list_size,nucl_n_aos(a))       |                  | Array to store just the useful mo_coef * ao_overlap      |
-! |                                           |                  | values depending of the nuclei                           |
-! | a                                         |                  | index to loop over the nuclei                            |
-! | b                                         |                  | index to loop over the AOs which belongs to the nuclei a |
-! | mu                                        |                  | index to refer to an AO which belongs to the nuclei a    |
-! | rho                                       |                  | index to loop over all the AOs                           |
-
-
-subroutine gradient_PM(tmp_n, tmp_list_size, tmp_list, v_grad, max_elem, norm_grad)
-
-  implicit none
-
-  BEGIN_DOC
-  ! Compute gradient for the Pipek-Mezey localization
-  END_DOC
-  
-  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
-  double precision, intent(out) :: v_grad(tmp_n), max_elem, norm_grad
-  double precision, allocatable :: m_grad(:,:), tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
-  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho
-  double precision              :: t1,t2,t3
-
-  print*,''
-  print*,'---gradient_PM---'
-
-  call wall_time(t1)
-
-  ! Allocation
-  allocate(m_grad(tmp_list_size, tmp_list_size), tmp_int(tmp_list_size, tmp_list_size),tmp_accu(tmp_list_size, tmp_list_size))
-  allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size))
-
-
-  ! submatrix of the mo_coef
-  do tmp_i = 1, tmp_list_size
-    i = tmp_list(tmp_i)
-    do j = 1, ao_num
-
-      tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
- 
-    enddo
-  enddo
-
-  call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
- 
-  m_grad = 0d0
-
-  do a = 1, nucl_num ! loop over the nuclei
-    tmp_int = 0d0
-
-    !do tmp_j = 1, tmp_list_size
-    !  do tmp_i = 1, tmp_list_size
-    !    do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-    !      mu = nucl_aos(a,b)
-
-    !      tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu))
-
-    !                             !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
-    !                             !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
-
-    !    enddo
-    !  enddo
-    !enddo
-
-    allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
-
-    do tmp_i = 1, tmp_list_size
-      do b = 1, nucl_n_aos(a)
-        mu = nucl_aos(a,b)
-
-        tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i)
-
-      enddo
-    enddo
-    
-    do b = 1, nucl_n_aos(a)
-      mu = nucl_aos(a,b)
-      do tmp_i = 1, tmp_list_size
-
-        tmp_CS(tmp_i,b) = CS(tmp_i,mu)
-
-      enddo
-    enddo   
-
-    call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
- 
-    do tmp_j = 1, tmp_list_size
-      do tmp_i = 1, tmp_list_size
-
-        tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))   
-
-      enddo
-    enddo
-
-    deallocate(tmp_mo_coef2,tmp_CS)
-
-    do tmp_j = 1, tmp_list_size
-      do tmp_i = 1, tmp_list_size
-
-        m_grad(tmp_i,tmp_j) = m_grad(tmp_i,tmp_j) +  4d0 * tmp_int(tmp_i,tmp_j) * (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))
-
-      enddo
-    enddo
-
-  enddo
-
-  ! 2D -> 1D
-  do tmp_k = 1, tmp_n
-    call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
-    v_grad(tmp_k) = m_grad(tmp_i,tmp_j) 
-  enddo
-
-  ! Maximum element in the gradient
-  max_elem = 0d0
-  do tmp_k = 1, tmp_n
-    if (ABS(v_grad(tmp_k)) > max_elem) then
-      max_elem = ABS(v_grad(tmp_k))
-    endif
-  enddo 
-
-  ! Norm of the gradient
-  norm_grad = 0d0
-  do tmp_k = 1, tmp_n
-    norm_grad = norm_grad + v_grad(tmp_k)**2
-  enddo
-  norm_grad = dsqrt(norm_grad)
-
-  print*, 'Maximal element in the gradient:', max_elem
-  print*, 'Norm of the gradient:', norm_grad
-
-  ! Deallocation
-  deallocate(m_grad,tmp_int,CS,tmp_mo_coef)
-
-  call wall_time(t2)
-  t3 = t2 - t1
-  print*,'Time in gradient_PM:', t3
-
-  print*,'---End gradient_PM---'
-
-end
-
-! Hessian v1
-
-subroutine hess_pipek(tmp_n, tmp_list_size, tmp_list, H)
-
-  implicit none
-
-  BEGIN_DOC
-  ! Compute diagonal hessian for the Pipek-Mezey localization
-  END_DOC
-
-  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
-  double precision, intent(out) :: H(tmp_n)
-  double precision, allocatable :: beta(:,:),tmp_int(:,:)
-  integer                       :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu
-  double precision              :: max_elem
-    
-  ! Allocation
-  allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size))
-
-  beta = 0d0
-
-  do a = 1, nucl_num
-    tmp_int = 0d0
-
-    do tmp_j = 1, tmp_list_size
-      j = tmp_list(tmp_j)
-      do tmp_i = 1, tmp_list_size
-        i = tmp_list(tmp_i)
-        do rho = 1, ao_num
-          do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-            mu = nucl_aos(a,b)
-
-            tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
-                                   + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
-
-          enddo
-        enddo  
-      enddo
-    enddo
-
-    ! Calculation
-    do tmp_j = 1, tmp_list_size
-      do tmp_i = 1, tmp_list_size
-
-        beta(tmp_i,tmp_j) = beta(tmp_i, tmp_j) +  (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))**2 - 4d0 * tmp_int(tmp_i,tmp_j)**2
-
-      enddo
-    enddo
-  
-  enddo
-
-  H = 0d0
-  do tmp_k = 1, tmp_n
-    call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
-    H(tmp_k) = 4d0 * beta(tmp_i, tmp_j)
-  enddo
-  
-  ! Deallocation
-  deallocate(beta,tmp_int)
-
-end
-
-! Hessian
-
-! The hessian is
-! \begin{align*}
-! \left. \frac{\partial^2 \mathcal{P} (\theta)}{\partial \theta^2}\right|_{\theta=0} = 4 \beta^{PM}
-! \end{align*}
-! \begin{align*}
-! \beta_{st}^{PM} = \sum_{A=1}^N \left( <s|P_A|t>^2 - \frac{1}{4} \left[<s|P_A|s> - <t|P_A|t> \right]^2 \right)
-! \end{align*}
-
-! with
-! \begin{align*}
-! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right]
-! \end{align*}
-! $\sum_{\rho}$ -> sum over all the AOs
-! $\sum_{\mu \in A}$ -> sum over the AOs which belongs to atom A
-! $c^t$ -> expansion coefficient of orbital |t>
-
-
-subroutine hessian_PM(tmp_n, tmp_list_size, tmp_list, H)
-
-  implicit none
-  
-  BEGIN_DOC
-  ! Compute diagonal hessian for the Pipek-Mezey localization
-  END_DOC
-
-  integer, intent(in)           :: tmp_n, tmp_list_size, tmp_list(tmp_list_size)
-  double precision, intent(out) :: H(tmp_n)
-  double precision, allocatable :: beta(:,:),tmp_int(:,:),CS(:,:),tmp_mo_coef(:,:),tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
-  integer                       :: i,j,tmp_k,tmp_i, tmp_j, a,b,rho,mu
-  double precision              :: max_elem, t1,t2,t3
-    
-  print*,''
-  print*,'---hessian_PM---'
-
-  call wall_time(t1)
-
-  ! Allocation
-  allocate(beta(tmp_list_size,tmp_list_size),tmp_int(tmp_list_size,tmp_list_size),tmp_accu(tmp_list_size,tmp_list_size))
-  allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size))
-
-  beta = 0d0
-
-  do tmp_i = 1, tmp_list_size
-    i = tmp_list(tmp_i)
-    do j = 1, ao_num
-
-      tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
- 
-    enddo
-  enddo
-
-  call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
- 
-  do a = 1, nucl_num ! loop over the nuclei
-    tmp_int = 0d0
-
-    !do tmp_j = 1, tmp_list_size
-    !  do tmp_i = 1, tmp_list_size
-    !    do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-    !      mu = nucl_aos(a,b)
-
-    !      tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu))
-
-    !                             !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
-    !                             !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
-
-    !    enddo
-    !  enddo
-    !enddo
- 
-    allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
-
-    do tmp_i = 1, tmp_list_size
-      do b = 1, nucl_n_aos(a)
-        mu = nucl_aos(a,b)
-
-        tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i)
-
-      enddo
-    enddo
-    
-    do b = 1, nucl_n_aos(a)
-      mu = nucl_aos(a,b)
-      do tmp_i = 1, tmp_list_size
-
-        tmp_CS(tmp_i,b) = CS(tmp_i,mu)
-
-      enddo
-    enddo   
-
-    call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
- 
-    do tmp_j = 1, tmp_list_size
-      do tmp_i = 1, tmp_list_size
-
-        tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))   
-
-      enddo
-    enddo
-
-    deallocate(tmp_mo_coef2,tmp_CS)
-
-    ! Calculation
-    do tmp_j = 1, tmp_list_size
-      do tmp_i = 1, tmp_list_size
-
-        beta(tmp_i,tmp_j) = beta(tmp_i, tmp_j) +  (tmp_int(tmp_i,tmp_i) - tmp_int(tmp_j,tmp_j))**2 - 4d0 * tmp_int(tmp_i,tmp_j)**2
-
-      enddo
-    enddo
-  
-  enddo
-
-  H = 0d0
-  do tmp_k = 1, tmp_n
-    call vec_to_mat_index(tmp_k,tmp_i,tmp_j)
-    H(tmp_k) = 4d0 * beta(tmp_i, tmp_j)
-  enddo
-  
-  ! Deallocation
-  deallocate(beta,tmp_int)
-
-  call wall_time(t2)
-  t3 = t2 - t1
-  print*,'Time in hessian_PM:', t3
-
-  print*,'---End hessian_PM---'
-
-end
-
-! Criterion PM (old)
-
-subroutine compute_crit_pipek(criterion)
-
-  implicit none
-  
-  BEGIN_DOC
-  ! Compute the Pipek-Mezey localization criterion
-  END_DOC
-
-  double precision, intent(out) :: criterion
-  double precision, allocatable :: tmp_int(:,:)
-  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho 
-
-  ! Allocation
-  allocate(tmp_int(mo_num, mo_num))
- 
-  criterion = 0d0
-
-  do a = 1, nucl_num ! loop over the nuclei
-    tmp_int = 0d0
-
-    do i = 1, mo_num
-      do rho = 1, ao_num ! loop over all the AOs
-        do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-          mu = nucl_aos(a,b)
-
-          tmp_int(i,i) = tmp_int(i,i) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,i) &
-                                 + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,i))
-
-        enddo
-      enddo  
-    enddo
-
-    do i = 1, mo_num 
-      criterion = criterion + tmp_int(i,i)**2
-    enddo
-
-  enddo
-  
-  criterion = - criterion 
-
-  deallocate(tmp_int)
-
-end
-
-! Criterion PM
-
-! The criterion is computed as
-! \begin{align*}
-! \mathcal{P} = \sum_{i=1}^n \sum_{A=1}^N \left[ <i|P_A|i> \right]^2
-! \end{align*}
-! with 
-! \begin{align*}
-! <s|P_A|t> = \frac{1}{2} \sum_{\rho} \sum_{\mu \in A} \left[ c_{\rho}^{s*} S_{\rho \nu} c_{\mu}^{t} +c_{\mu}^{s*} S_{\mu \rho} c_{\rho}^t \right]
-! \end{align*}
-
-
-subroutine criterion_PM(tmp_list_size,tmp_list,criterion)
-
-  implicit none
-  
-  BEGIN_DOC
-  ! Compute the Pipek-Mezey localization criterion
-  END_DOC
-
-  integer, intent(in)           :: tmp_list_size, tmp_list(tmp_list_size)
-  double precision, intent(out) :: criterion
-  double precision, allocatable :: tmp_int(:,:),CS(:,:)
-  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho
-  
-  print*,''
-  print*,'---criterion_PM---'
-  
-  ! Allocation
-  allocate(tmp_int(tmp_list_size, tmp_list_size),CS(mo_num,ao_num))
-  
-  ! Initialization
-  criterion = 0d0
-
-  call dgemm('T','N',mo_num,ao_num,ao_num,1d0,mo_coef,size(mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
- 
-  do a = 1, nucl_num ! loop over the nuclei
-    tmp_int = 0d0
-
-      do tmp_i = 1, tmp_list_size
-        i = tmp_list(tmp_i)
-        do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-          mu = nucl_aos(a,b)
-          
-          tmp_int(tmp_i,tmp_i) = tmp_int(tmp_i,tmp_i) + 0.5d0 * (CS(i,mu) * mo_coef(mu,i) + mo_coef(mu,i) * CS(i,mu))
-
-                                 !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
-                                 !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
-
-      enddo
-    enddo
-
-    do tmp_i = 1, tmp_list_size
-      criterion = criterion + tmp_int(tmp_i,tmp_i)**2
-    enddo
-
-  enddo
-  
-  criterion = - criterion 
-
-  deallocate(tmp_int,CS)
-
-  print*,'---End criterion_PM---'
-
-end
-
-! Criterion PM v3
-
-subroutine criterion_PM_v3(tmp_list_size,tmp_list,criterion)
-
-  implicit none
-  
-  BEGIN_DOC
-  ! Compute the Pipek-Mezey localization criterion
-  END_DOC
-  
-  integer, intent(in)           :: tmp_list_size, tmp_list(tmp_list_size)
-  double precision, intent(out) :: criterion
-  double precision, allocatable :: tmp_int(:,:), CS(:,:), tmp_mo_coef(:,:), tmp_mo_coef2(:,:),tmp_accu(:,:),tmp_CS(:,:)
-  integer                       :: i,j,k,tmp_i,tmp_j,tmp_k, a, b, mu ,rho,nu,c
-  double precision              :: t1,t2,t3
-
-  print*,''
-  print*,'---criterion_PM_v3---'
-  
-  call wall_time(t1)
-
-  ! Allocation
-  allocate(tmp_int(tmp_list_size, tmp_list_size),tmp_accu(tmp_list_size, tmp_list_size))
-  allocate(CS(tmp_list_size,ao_num),tmp_mo_coef(ao_num,tmp_list_size))
-
-  criterion = 0d0
-
-  ! submatrix of the mo_coef
-  do tmp_i = 1, tmp_list_size
-    i = tmp_list(tmp_i)
-    do j = 1, ao_num
-
-      tmp_mo_coef(j,tmp_i) = mo_coef(j,i)
- 
-    enddo
-  enddo
-
-  ! ao_overlap(ao_num,ao_num)
-  ! mo_coef(ao_num,mo_num)
-  call dgemm('T','N',tmp_list_size,ao_num,ao_num,1d0,tmp_mo_coef,size(tmp_mo_coef,1),ao_overlap,size(ao_overlap,1),0d0,CS,size(CS,1))
- 
-  do a = 1, nucl_num ! loop over the nuclei
-  
-    do j = 1, tmp_list_size
-      do i = 1, tmp_list_size
-        tmp_int(i,j) = 0d0
-      enddo
-    enddo
-
-    !do tmp_j = 1, tmp_list_size
-    !  do tmp_i = 1, tmp_list_size
-    !    do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-    !      mu = nucl_aos(a,b)
-
-    !      tmp_int(tmp_i,tmp_j) = tmp_int(tmp_i,tmp_j) + 0.5d0 * (CS(tmp_i,mu) * tmp_mo_coef(mu,tmp_j) + tmp_mo_coef(mu,tmp_i) * CS(tmp_j,mu))
-
-    !                             !  (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
-    !                             !+ mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
-
-    !    enddo
-    !  enddo
-    !enddo
-
-    allocate(tmp_mo_coef2(nucl_n_aos(a),tmp_list_size),tmp_CS(tmp_list_size,nucl_n_aos(a)))
- 
-    do tmp_i = 1, tmp_list_size
-      do b = 1, nucl_n_aos(a)
-        mu = nucl_aos(a,b)
-
-        tmp_mo_coef2(b,tmp_i) = tmp_mo_coef(mu,tmp_i)
-
-      enddo
-    enddo
-    
-    do b = 1, nucl_n_aos(a)
-      mu = nucl_aos(a,b)
-      do tmp_i = 1, tmp_list_size
-
-         tmp_CS(tmp_i,b) = CS(tmp_i,mu)
-
-      enddo
-    enddo   
-
-    call dgemm('N','N',tmp_list_size,tmp_list_size,nucl_n_aos(a),1d0,tmp_CS,size(tmp_CS,1),tmp_mo_coef2,size(tmp_mo_coef2,1),0d0,tmp_accu,size(tmp_accu,1))
- 
-    ! Integrals
-    do tmp_j = 1, tmp_list_size
-      do tmp_i = 1, tmp_list_size
-
-        tmp_int(tmp_i,tmp_j) = 0.5d0 * (tmp_accu(tmp_i,tmp_j) + tmp_accu(tmp_j,tmp_i))   
-
-      enddo
-    enddo
-
-    deallocate(tmp_mo_coef2,tmp_CS)
-
-    ! Criterion
-    do tmp_i = 1, tmp_list_size
-      criterion = criterion + tmp_int(tmp_i,tmp_i)**2
-    enddo
-
-  enddo
-
-  criterion = - criterion 
-
-  deallocate(tmp_int,CS,tmp_accu,tmp_mo_coef)
-
-  call wall_time(t2)
-  t3 = t2 - t1
-  print*,'Time in criterion_PM_v3:', t3
-
-  print*,'---End criterion_PM_v3---'
-
-end
-
-! Criterion FB (old)
-
-! The criterion is just computed as
-
-! \begin{align*}
-! C = - \sum_i^{mo_{num}} (<i|x|i>^2 + <i|y|i>^2 + <i|z|i>^2)
-! \end{align*}
-
-! The minus sign is here in order to minimize this criterion
-
-! Output:
-! | criterion | double precision | criterion for the Foster-Boys localization |
-
-
-subroutine criterion_FB_old(criterion)
-
-  implicit none
-
-  BEGIN_DOC
-  ! Compute the Foster-Boys localization criterion
-  END_DOC
-
-  double precision, intent(out) :: criterion
-  integer                       :: i
-
-  ! Criterion (= \sum_i <i|r|i>^2 )
-  criterion = 0d0
-  do i = 1, mo_num
-    criterion = criterion + mo_dipole_x(i,i)**2 + mo_dipole_y(i,i)**2 + mo_dipole_z(i,i)**2
-  enddo
-  criterion = - criterion
-
-end subroutine
-
-! Criterion FB
-
-subroutine criterion_FB(tmp_list_size, tmp_list, criterion)
-
-  implicit none
-
-  BEGIN_DOC
-  ! Compute the Foster-Boys localization criterion
-  END_DOC
-
-  integer, intent(in)           :: tmp_list_size, tmp_list(tmp_list_size)
-  double precision, intent(out) :: criterion
-  integer                       :: i, tmp_i
-
-  ! Criterion (= - \sum_i <i|r|i>^2 )
-  criterion = 0d0
-  do tmp_i = 1, tmp_list_size
-    i = tmp_list(tmp_i)
-    criterion = criterion + mo_dipole_x(i,i)**2 + mo_dipole_y(i,i)**2 + mo_dipole_z(i,i)**2
-  enddo
-  criterion = - criterion
-
-end subroutine
-
-subroutine theta_FB(l, n, m_x, max_elem)
-
-  include 'pi.h'
-
-  BEGIN_DOC
-  ! Compute the angles to minimize the Foster-Boys criterion by using pairwise rotations of the MOs
-  ! Warning: you must give - the angles to build the rotation matrix...
-  END_DOC
-
-  implicit none
-
-  integer, intent(in)           :: n, l(n)
-  double precision, intent(out) :: m_x(n,n), max_elem
-
-  integer                       :: i,j, tmp_i, tmp_j
-  double precision, allocatable :: cos4theta(:,:), sin4theta(:,:)
-  double precision, allocatable :: A(:,:), B(:,:), beta(:,:), gamma(:,:)
-  integer                       :: idx_i,idx_j
-
-  allocate(cos4theta(n, n), sin4theta(n, n))
-  allocate(A(n,n), B(n,n), beta(n,n), gamma(n,n))
-
-  do tmp_j = 1, n
-    j = l(tmp_j)
-    do tmp_i = 1, n
-      i = l(tmp_i)
-      A(tmp_i,tmp_j) = mo_dipole_x(i,j)**2 - 0.25d0 * (mo_dipole_x(i,i) - mo_dipole_x(j,j))**2 &
-                     + mo_dipole_y(i,j)**2 - 0.25d0 * (mo_dipole_y(i,i) - mo_dipole_y(j,j))**2 &
-                     + mo_dipole_z(i,j)**2 - 0.25d0 * (mo_dipole_z(i,i) - mo_dipole_z(j,j))**2
-    enddo
-    A(j,j) = 0d0
-  enddo
-
-  do tmp_j = 1, n
-    j = l(tmp_j)
-    do tmp_i = 1, n
-      i = l(tmp_i)
-      B(tmp_i,tmp_j) = mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) &
-                     + mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) &
-                     + mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j))
-    enddo
-  enddo
-
-  !do tmp_j = 1, n
-  !  j = l(tmp_j)
-  !  do tmp_i = 1, n
-  !    i = l(tmp_i)
-  !    beta(tmp_i,tmp_j) =  (mo_dipole_x(i,i) - mo_dipole_x(j,j)) - 4d0 * mo_dipole_x(i,j)**2 &
-  !               + (mo_dipole_y(i,i) - mo_dipole_y(j,j)) - 4d0 * mo_dipole_y(i,j)**2 &
-  !               + (mo_dipole_z(i,i) - mo_dipole_z(j,j)) - 4d0 * mo_dipole_z(i,j)**2
-  !  enddo
-  !enddo
-
-  !do tmp_j = 1, n
-  !  j = l(tmp_j)
-  !  do tmp_i = 1, n
-  !    i = l(tmp_i)
-  !    gamma(tmp_i,tmp_j) = 4d0 * ( mo_dipole_x(i,j) * (mo_dipole_x(i,i) - mo_dipole_x(j,j)) &
-  !                       + mo_dipole_y(i,j) * (mo_dipole_y(i,i) - mo_dipole_y(j,j)) &
-  !                       + mo_dipole_z(i,j) * (mo_dipole_z(i,i) - mo_dipole_z(j,j)))
-  !  enddo
-  !enddo
-
-  !
-  !do j = 1, n
-  !  do i = 1, n
-  !    cos4theta(i,j) = - A(i,j) / dsqrt(A(i,j)**2 + B(i,j)**2)
-  !  enddo
-  !enddo
-
-  !do j = 1, n
-  !  do i = 1, n
-  !    sin4theta(i,j) = B(i,j) / dsqrt(A(i,j)**2 + B(i,j)**2)
-  !  enddo
-  !enddo
-
-  ! Theta
-  do j = 1, n
-    do i = 1, n
-      m_x(i,j) = 0.25d0 * atan2(B(i,j), -A(i,j))
-      !m_x(i,j) = 0.25d0 * atan2(sin4theta(i,j), cos4theta(i,j))
-    enddo
-  enddo
-
-  ! Enforce a perfect antisymmetry
-  do j = 1, n-1
-    do i = j+1, n
-      m_x(j,i) = - m_x(i,j)
-    enddo
-  enddo
-  do i = 1, n
-    m_x(i,i) = 0d0
-  enddo
-
-  ! Max
-  max_elem = 0d0
-  do j = 1, n-1
-    do i = j+1, n
-      if (dabs(m_x(i,j)) > dabs(max_elem)) then
-        max_elem = m_x(i,j)
-        !idx_i = i
-        !idx_j = j
-      endif
-    enddo
-  enddo
-
-  ! Debug
-  !print*,''
-  !print*,'sin/B'
-  !do i = 1, n
-  !  write(*,'(100F10.4)') sin4theta(i,:)
-  !  !B(i,:)
-  !enddo
-  !print*,'cos/A'
-  !do i = 1, n
-  !  write(*,'(100F10.4)') cos4theta(i,:)
-  !  !A(i,:)
-  !enddo
-  !print*,'X'
-  !!m_x = 0d0
-  !!m_x(idx_i,idx_j) = max_elem
-  !!m_x(idx_j,idx_i) = -max_elem
-  !do i = 1, n
-  !  write(*,'(100F10.4)') m_x(i,:)
-  !enddo
-  !print*,idx_i,idx_j,max_elem
-
-  max_elem = dabs(max_elem)
-  
-  deallocate(cos4theta, sin4theta)
-  deallocate(A,B,beta,gamma)
-  
-end
-
-subroutine theta_PM(l, n, m_x, max_elem)
-  
-  include 'pi.h'
-
-  BEGIN_DOC
-  ! Compute the angles to minimize the Foster-Boys criterion by using pairwise rotations of the MOs
-  ! Warning: you must give - the angles to build the rotation matrix...
-  END_DOC
-
-  implicit none
-
-  integer, intent(in)           :: n, l(n)
-  double precision, intent(out) :: m_x(n,n), max_elem
-
-  integer                       :: a,b,i,j,tmp_i,tmp_j,rho,mu,nu,idx_i,idx_j
-  double precision, allocatable :: Aij(:,:), Bij(:,:), Pa(:,:)
-
-  allocate(Aij(n,n), Bij(n,n), Pa(n,n))
-
-  do a = 1, nucl_num ! loop over the nuclei
-    Pa = 0d0 ! Initialization for each nuclei
-
-    ! Loop over the MOs of the a given mo_class to compute <i|P_a|j>
-    do tmp_j = 1, n
-      j = l(tmp_j) 
-      do tmp_i = 1, n
-         i = l(tmp_i)
-        do rho = 1, ao_num ! loop over all the AOs
-          do b = 1, nucl_n_aos(a) ! loop over the number of AOs which belongs to the nuclei a
-            mu = nucl_aos(a,b) ! AO centered on atom a 
-
-            Pa(tmp_i,tmp_j) = Pa(tmp_i,tmp_j) + 0.5d0 * (mo_coef(rho,i) * ao_overlap(rho,mu) * mo_coef(mu,j) &
-                                   + mo_coef(mu,i) * ao_overlap(mu,rho) * mo_coef(rho,j))
-
-          enddo
-        enddo  
-      enddo
-    enddo
-
-    ! A
-    do j = 1, n
-      do i = 1, n
-        Aij(i,j) = Aij(i,j) + Pa(i,j)**2 - 0.25d0 * (Pa(i,i) - Pa(j,j))**2
-      enddo
-    enddo
-    
-    ! B
-    do j = 1, n
-      do i = 1, n
-        Bij(i,j) = Bij(i,j) + Pa(i,j) * (Pa(i,i) - Pa(j,j))
-      enddo
-    enddo
-
-  enddo
-
-  ! Theta
-  do j = 1, n
-    do i = 1, n
-      m_x(i,j) = 0.25d0 * atan2(Bij(i,j), -Aij(i,j))
-    enddo
-  enddo
-
-  ! Enforce a perfect antisymmetry
-  do j = 1, n-1
-    do i = j+1, n
-      m_x(j,i) = - m_x(i,j)
-    enddo
-  enddo
-  do i = 1, n
-    m_x(i,i) = 0d0
-  enddo
-
-  ! Max
-  max_elem = 0d0
-  do j = 1, n-1
-    do i = j+1, n
-      if (dabs(m_x(i,j)) > dabs(max_elem)) then
-        max_elem = m_x(i,j)
-        idx_i = i
-        idx_j = j
-      endif
-    enddo
-  enddo
-
-  ! Debug
-  !do i = 1, n
-  !  write(*,'(100F10.4)') m_x(i,:)
-  !enddo
-  !print*,'Max',idx_i,idx_j,max_elem
-
-  max_elem = dabs(max_elem)
-
-  deallocate(Aij,Bij,Pa)
-
-end
-
 ! Spatial extent
 
 ! The spatial extent of an orbital $i$ is computed as 
@@ -1465,6 +496,7 @@ end
 ! From that we can also compute the average and the standard deviation
 
 
+
 subroutine compute_spatial_extent(spatial_extent)
 
   implicit none

src/ao_cart_basis/EZFIO.cfg

@@ -0,0 +1,75 @@
+[ao_cart_basis]
+type: character*(256)
+doc: Name of the |AO| basis set
+interface: ezfio
+
+[ao_cart_num]
+type: integer
+doc: Number of |AOs|
+interface: ezfio, provider
+
+[ao_cart_prim_num]
+type: integer
+doc: Number of primitives per |AO|
+size: (ao_cart_basis.ao_cart_num)
+interface: ezfio, provider
+
+[ao_cart_prim_num_max]
+type: integer
+doc: Maximum number of primitives
+default: =maxval(ao_cart_basis.ao_cart_prim_num)
+interface: ezfio
+
+[ao_cart_nucl]
+type: integer
+doc: Index of the nucleus on which the |AO| is centered
+size: (ao_cart_basis.ao_cart_num)
+interface: ezfio, provider
+
+[ao_cart_power]
+type: integer
+doc: Powers of x, y and z for each |AO|
+size: (ao_cart_basis.ao_cart_num,3)
+interface: ezfio, provider
+
+[ao_cart_coef]
+type: double precision
+doc: Primitive coefficients, read from input. Those should not be used directly, as the MOs are expressed on the basis of **normalized** AOs.
+size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_prim_num_max)
+interface: ezfio, provider
+
+[ao_cart_expo]
+type: double precision
+doc: Exponents for each primitive of each |AO|
+size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_prim_num_max)
+interface: ezfio, provider
+
+[ao_cart_md5]
+type: character*(32)
+doc: MD5 key, specific of the |AO| basis
+interface: ezfio, provider
+
+[ao_cart_cartesian]
+type: logical
+doc: If |true|, use |AOs| in Cartesian coordinates (6d,10f,...)
+interface: ezfio, provider
+default: false
+
+[ao_cart_expo_im]
+type: double precision
+doc: imag part for Exponents for each primitive of each cGTOs |AO|
+size: (ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_prim_num_max)
+interface: ezfio, provider
+
+[ao_cart_expo_pw]
+type: double precision
+doc: plane wave part for each primitive GTOs |AO|
+size: (3,ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_prim_num_max)
+interface: ezfio, provider
+
+[ao_cart_expo_phase]
+type: double precision
+doc: phase shift for each primitive GTOs |AO|
+size: (3,ao_cart_basis.ao_cart_num,ao_cart_basis.ao_cart_prim_num_max)
+interface: ezfio, provider
+

src/ao_cart_basis/README.rst

@@ -0,0 +1,45 @@
+========
+ao_basis
+========
+
+This module describes the atomic orbitals basis set.
+
+An |AO| :math:`\chi` centered on nucleus A is represented as:
+
+.. math::
+
+   \chi_i({\bf r}) = (x-X_A)^a (y-Y_A)^b (z-Z_A)^c \sum_k c_{ki} e^{-\gamma_{ki} |{\bf r} - {\bf R}_A|^2}
+
+
+The |AO| coefficients are normalized as:
+
+.. math::
+
+  {\tilde c}_{ki} = \frac{c_{ki}}{ \int \left( (x-X_A)^a (y-Y_A)^b (z-Z_A)^c  e^{-\gamma_{ki} |{\bf r} - {\bf R}_A|^2} \right)^2 dr}
+
+
+.. warning::
+
+  `ao_coef` contains the |AO| coefficients given in input. These do not
+  include the normalization constant of the |AO|. The `ao_coef_normalized`
+  provider includes this normalization factor.
+
+
+The |AOs| are also sorted by increasing exponent to accelerate the calculation of
+the two electron integrals.
+
+
+
+Complex Gaussian-Type Orbitals (cGTOs)
+=====================================
+
+Complex Gaussian-Type Orbitals (cGTOs) are also supported:
+
+.. math::
+
+   \chi_i(\mathbf{r}) = x^a y^b z^c \sum_k c_{ki} \left( e^{-\alpha_{ki} \mathbf{r}^2 - \imath \mathbf{k}_{ki} \cdot \mathbf{r} - \imath \phi_{ki}} + \text{C.C.} \right)
+
+where:
+   - :math:`\alpha \in \mathbb{C}` and :math:`\Re(\alpha) > 0` (specified by ``ao_expo`` and ``ao_expo_im``),
+   - :math:`\mathbf{k} = (k_x, k_y, k_z) \in \mathbb{R}^3` (specified by ``ao_expo_pw``),
+   - :math:`\phi = \phi_x + \phi_y + \phi_z \in \mathbb{R}` (specified by ``ao_expo_phase``).

src/ao_cart_basis/aos.irp.f

@@ -0,0 +1,445 @@
+BEGIN_PROVIDER [ integer, ao_cart_prim_num_max ]
+ implicit none
+ BEGIN_DOC
+ ! Max number of primitives.
+ END_DOC
+ ao_cart_prim_num_max = maxval(ao_cart_prim_num)
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_cart_shell, (ao_cart_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the shell to which the AO corresponds
+ END_DOC
+ integer :: i, j, k, n
+ k=0
+ do i=1,shell_num
+   n = shell_ang_mom(i)+1
+   do j=1,(n*(n+1))/2
+     k = k+1
+     ao_cart_shell(k) = i
+   enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_cart_sphe_num ]
+ implicit none
+ BEGIN_DOC
+ ! Number of spherical AOs
+ END_DOC
+ integer :: n, i
+ if (ao_cart_cartesian) then
+   ao_cart_sphe_num = ao_cart_num
+ else
+   ao_cart_sphe_num=0
+   do i=1,shell_num
+     n = shell_ang_mom(i)
+     ao_cart_sphe_num += 2*n+1
+   enddo
+ endif
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_cart_sphe_shell, (ao_cart_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the shell to which the AO corresponds
+ END_DOC
+ integer :: i, j, k, n
+ k=0
+ do i=1,shell_num
+   n = shell_ang_mom(i)
+   do j=-n,n
+     k = k+1
+     ao_cart_sphe_shell(k) = i
+   enddo
+ enddo
+END_PROVIDER
+
+BEGIN_PROVIDER [ integer, ao_cart_first_of_shell, (shell_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Index of the shell to which the AO corresponds
+ END_DOC
+ integer :: i, j, k, n
+ k=1
+ do i=1,shell_num
+   ao_cart_first_of_shell(i) = k
+   n = shell_ang_mom(i)+1
+   k = k+(n*(n+1))/2
+ enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, ao_cart_coef_normalized, (ao_cart_num,ao_cart_prim_num_max) ]
+&BEGIN_PROVIDER [ double precision, ao_cart_coef_normalization_factor, (ao_cart_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Coefficients including the |AO| normalization
+  END_DOC
+  double precision               :: norm,overlap_x,overlap_y,overlap_z,C_A(3), c
+  integer                        :: l, powA(3), nz
+  integer                        :: i,j,k
+  nz=100
+  C_A(1) = 0.d0
+  C_A(2) = 0.d0
+  C_A(3) = 0.d0
+  ao_cart_coef_normalized = 0.d0
+
+  if (primitives_normalized) then
+
+    if (ezfio_convention >= 20250211) then
+      ! Same primitive normalization factors for all AOs of the same shell, or read from trexio file
+
+      do i=1,ao_cart_num
+        k=1
+        do while (k<=prim_num .and. shell_index(k) /= ao_cart_shell(i))
+          k = k+1
+        end do
+        do j=1,ao_cart_prim_num(i)
+          ao_cart_coef_normalized(i,j) = ao_cart_coef(i,j)*prim_normalization_factor(k+j-1)
+        enddo
+      enddo
+
+    else
+        ! GAMESS convention for primitive factors
+
+      do i=1,ao_cart_num
+        powA(1) = ao_cart_power(i,1)
+        powA(2) = ao_cart_power(i,2)
+        powA(3) = ao_cart_power(i,3)
+
+        do j=1,ao_cart_prim_num(i)
+          call overlap_gaussian_xyz(C_A,C_A,ao_cart_expo(i,j),ao_cart_expo(i,j), &
+             powA,powA,overlap_x,overlap_y,overlap_z,norm,nz)
+          ao_cart_coef_normalized(i,j) = ao_cart_coef(i,j)/dsqrt(norm)
+        enddo
+      enddo
+
+    endif
+
+  else
+
+    do i=1,ao_cart_num
+      do j=1,ao_cart_prim_num(i)
+        ao_cart_coef_normalized(i,j) = ao_cart_coef(i,j)
+      enddo
+    enddo
+
+  endif
+
+  double precision, allocatable :: self_overlap(:)
+  allocate(self_overlap(ao_cart_num))
+
+  do i=1,ao_cart_num
+    powA(1) = ao_cart_power(i,1)
+    powA(2) = ao_cart_power(i,2)
+    powA(3) = ao_cart_power(i,3)
+    self_overlap(i) = 0.d0
+    do j=1,ao_cart_prim_num(i)
+      do k=1,j-1
+        call overlap_gaussian_xyz(C_A,C_A,ao_cart_expo(i,j),ao_cart_expo(i,k),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+        self_overlap(i) = self_overlap(i) + 2.d0*c*ao_cart_coef_normalized(i,j)*ao_cart_coef_normalized(i,k)
+      enddo
+      call overlap_gaussian_xyz(C_A,C_A,ao_cart_expo(i,j),ao_cart_expo(i,j),powA,powA,overlap_x,overlap_y,overlap_z,c,nz)
+      self_overlap(i) = self_overlap(i) +c*ao_cart_coef_normalized(i,j)*ao_cart_coef_normalized(i,j)
+    enddo
+  enddo
+
+  if (ao_cart_normalized) then
+
+    do i=1,ao_cart_num
+      ao_cart_coef_normalization_factor(i) = 1.d0/dsqrt(self_overlap(i))
+    enddo
+
+  else
+
+    do i=1,ao_cart_num
+      ao_cart_coef_normalization_factor(i) = 1.d0
+    enddo
+
+  endif
+
+  do i=1,ao_cart_num
+    do j=1,ao_cart_prim_num(i)
+      ao_cart_coef_normalized(i,j) = ao_cart_coef_normalized(i,j) * ao_cart_coef_normalization_factor(i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ao_cart_sphe_coef_normalization_factor, (ao_cart_sphe_num) ]
+ implicit none
+ BEGIN_DOC
+ ! Normalization factor in spherical AO basis
+ END_DOC
+
+ ao_cart_sphe_coef_normalization_factor(:) = 1.d0
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ double precision, ao_cart_coef_normalized_ordered, (ao_cart_num,ao_cart_prim_num_max) ]
+&BEGIN_PROVIDER [ double precision, ao_cart_expo_ordered, (ao_cart_num,ao_cart_prim_num_max) ]
+  implicit none
+  BEGIN_DOC
+  ! Sorted primitives to accelerate 4 index |MO| transformation
+  END_DOC
+
+  integer                        :: iorder(ao_cart_prim_num_max)
+  double precision               :: d(ao_cart_prim_num_max,2)
+  integer                        :: i,j
+  do i=1,ao_cart_num
+    do j=1,ao_cart_prim_num(i)
+      iorder(j) = j
+      d(j,1) = ao_cart_expo(i,j)
+      d(j,2) = ao_cart_coef_normalized(i,j)
+    enddo
+    call dsort(d(1,1),iorder,ao_cart_prim_num(i))
+    call dset_order(d(1,2),iorder,ao_cart_prim_num(i))
+    do j=1,ao_cart_prim_num(i)
+      ao_cart_expo_ordered(i,j) = d(j,1)
+      ao_cart_coef_normalized_ordered(i,j) = d(j,2)
+    enddo
+  enddo
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ double precision, ao_cart_coef_normalized_ordered_transp, (ao_cart_prim_num_max,ao_cart_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Transposed :c:data:`ao_cart_coef_normalized_ordered`
+  END_DOC
+  integer                        :: i,j
+  do j=1, ao_cart_num
+    do i=1, ao_cart_prim_num_max
+      ao_cart_coef_normalized_ordered_transp(i,j) = ao_cart_coef_normalized_ordered(j,i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ao_cart_expo_ordered_transp, (ao_cart_prim_num_max,ao_cart_num) ]
+  implicit none
+  BEGIN_DOC
+  ! Transposed :c:data:`ao_cart_expo_ordered`
+  END_DOC
+  integer                        :: i,j
+  do j=1, ao_cart_num
+    do i=1, ao_cart_prim_num_max
+      ao_cart_expo_ordered_transp(i,j) = ao_cart_expo_ordered(j,i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, ao_cart_l, (ao_cart_num) ]
+&BEGIN_PROVIDER [ integer, ao_cart_l_max  ]
+&BEGIN_PROVIDER [ character*(128), ao_cart_l_char, (ao_cart_num) ]
+ implicit none
+ BEGIN_DOC
+! :math:`l` value of the |AO|: :math`a+b+c` in :math:`x^a y^b z^c`
+ END_DOC
+ integer :: i
+ do i=1,ao_cart_num
+   ao_cart_l(i) = ao_cart_power(i,1) + ao_cart_power(i,2) + ao_cart_power(i,3)
+   ao_cart_l_char(i) = l_to_character(ao_cart_l(i))
+ enddo
+ ao_cart_l_max = maxval(ao_cart_l)
+END_PROVIDER
+
+integer function ao_cart_power_index(nx,ny,nz)
+  implicit none
+  integer, intent(in)            :: nx, ny, nz
+  BEGIN_DOC
+  ! Unique index given to a triplet of powers:
+  !
+  ! :math:`\frac{1}{2} (l-n_x) (l-n_x+1) + n_z + 1`
+  END_DOC
+  integer                        :: l
+  l = nx + ny + nz
+  ao_cart_power_index = ((l-nx)*(l-nx+1))/2 + nz + 1
+end
+
+
+BEGIN_PROVIDER [ character*(128), l_to_character, (0:7)]
+ BEGIN_DOC
+ ! Character corresponding to the "l" value of an |AO|
+ END_DOC
+ implicit none
+ l_to_character(0)='s'
+ l_to_character(1)='p'
+ l_to_character(2)='d'
+ l_to_character(3)='f'
+ l_to_character(4)='g'
+ l_to_character(5)='h'
+ l_to_character(6)='i'
+ l_to_character(7)='j'
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, Nucl_N_Aos, (nucl_num)]
+&BEGIN_PROVIDER [ integer, N_AOs_max ]
+ implicit none
+ BEGIN_DOC
+ ! Number of |AOs| per atom
+ END_DOC
+ integer :: i
+ Nucl_N_Aos = 0
+ do i = 1, ao_cart_num
+  Nucl_N_Aos(ao_cart_nucl(i)) +=1
+ enddo
+ N_AOs_max = maxval(Nucl_N_Aos)
+END_PROVIDER
+
+ BEGIN_PROVIDER [ integer, nucl_aos, (nucl_num,N_AOs_max)]
+ implicit none
+ BEGIN_DOC
+ ! List of |AOs| centered on each atom
+ END_DOC
+ integer :: i
+ integer, allocatable :: nucl_tmp(:)
+ allocate(nucl_tmp(nucl_num))
+ nucl_tmp = 0
+ Nucl_Aos = 0
+ do i = 1, ao_cart_num
+  nucl_tmp(ao_cart_nucl(i))+=1
+  Nucl_Aos(ao_cart_nucl(i),nucl_tmp(ao_cart_nucl(i))) = i
+ enddo
+ deallocate(nucl_tmp)
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ integer, Nucl_list_shell_Aos, (nucl_num,N_AOs_max)]
+&BEGIN_PROVIDER [ integer, Nucl_num_shell_Aos, (nucl_num)]
+ implicit none
+ integer :: i,j,k
+ BEGIN_DOC
+ ! Index of the shell type |AOs| and of the corresponding |AOs|
+ ! By convention, for p,d,f and g |AOs|, we take the index
+ ! of the |AO| with the the corresponding power in the x axis
+ END_DOC
+ do i = 1, nucl_num
+  Nucl_num_shell_Aos(i) = 0
+  do j = 1, Nucl_N_Aos(i)
+    if (ao_cart_power(Nucl_Aos(i,j),1) == ao_cart_l(Nucl_Aos(i,j))) then
+     Nucl_num_shell_Aos(i)+=1
+     Nucl_list_shell_Aos(i,Nucl_num_shell_Aos(i))=Nucl_Aos(i,j)
+    endif
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ character*(4), ao_cart_l_char_space, (ao_cart_num) ]
+ implicit none
+ BEGIN_DOC
+! Converts an l value to a string
+ END_DOC
+ integer :: i
+ character*(4) :: give_ao_cart_character_space
+ do i=1,ao_cart_num
+
+  if(ao_cart_l(i)==0)then
+  ! S type AO
+   give_ao_cart_character_space = 'S   '
+  elseif(ao_cart_l(i) == 1)then
+  ! P type AO
+   if(ao_cart_power(i,1)==1)then
+    give_ao_cart_character_space = 'X   '
+   elseif(ao_cart_power(i,2) == 1)then
+    give_ao_cart_character_space = 'Y   '
+   else
+    give_ao_cart_character_space = 'Z   '
+   endif
+  elseif(ao_cart_l(i) == 2)then
+  ! D type AO
+   if(ao_cart_power(i,1)==2)then
+    give_ao_cart_character_space = 'XX  '
+   elseif(ao_cart_power(i,2) == 2)then
+    give_ao_cart_character_space = 'YY  '
+   elseif(ao_cart_power(i,3) == 2)then
+    give_ao_cart_character_space = 'ZZ  '
+   elseif(ao_cart_power(i,1) == 1 .and. ao_cart_power(i,2) == 1)then
+    give_ao_cart_character_space = 'XY  '
+   elseif(ao_cart_power(i,1) == 1 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'XZ  '
+   else
+    give_ao_cart_character_space = 'YZ  '
+   endif
+  elseif(ao_cart_l(i) == 3)then
+  ! F type AO
+   if(ao_cart_power(i,1)==3)then
+    give_ao_cart_character_space = 'XXX '
+   elseif(ao_cart_power(i,2) == 3)then
+    give_ao_cart_character_space = 'YYY '
+   elseif(ao_cart_power(i,3) == 3)then
+    give_ao_cart_character_space = 'ZZZ '
+   elseif(ao_cart_power(i,1) == 2 .and. ao_cart_power(i,2) == 1)then
+    give_ao_cart_character_space = 'XXY '
+   elseif(ao_cart_power(i,1) == 2 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'XXZ '
+   elseif(ao_cart_power(i,2) == 2 .and. ao_cart_power(i,1) == 1)then
+    give_ao_cart_character_space = 'YYX '
+   elseif(ao_cart_power(i,2) == 2 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'YYZ '
+   elseif(ao_cart_power(i,3) == 2 .and. ao_cart_power(i,1) == 1)then
+    give_ao_cart_character_space = 'ZZX '
+   elseif(ao_cart_power(i,3) == 2 .and. ao_cart_power(i,2) == 1)then
+    give_ao_cart_character_space = 'ZZY '
+   elseif(ao_cart_power(i,3) == 1 .and. ao_cart_power(i,2) == 1 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'XYZ '
+   endif
+  elseif(ao_cart_l(i) == 4)then
+  ! G type AO
+   if(ao_cart_power(i,1)==4)then
+    give_ao_cart_character_space = 'XXXX'
+   elseif(ao_cart_power(i,2) == 4)then
+    give_ao_cart_character_space = 'YYYY'
+   elseif(ao_cart_power(i,3) == 4)then
+    give_ao_cart_character_space = 'ZZZZ'
+   elseif(ao_cart_power(i,1) == 3 .and. ao_cart_power(i,2) == 1)then
+    give_ao_cart_character_space = 'XXXY'
+   elseif(ao_cart_power(i,1) == 3 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'XXXZ'
+   elseif(ao_cart_power(i,2) == 3 .and. ao_cart_power(i,1) == 1)then
+    give_ao_cart_character_space = 'YYYX'
+   elseif(ao_cart_power(i,2) == 3 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'YYYZ'
+   elseif(ao_cart_power(i,3) == 3 .and. ao_cart_power(i,1) == 1)then
+    give_ao_cart_character_space = 'ZZZX'
+   elseif(ao_cart_power(i,3) == 3 .and. ao_cart_power(i,2) == 1)then
+    give_ao_cart_character_space = 'ZZZY'
+   elseif(ao_cart_power(i,1) == 2 .and. ao_cart_power(i,2) == 2)then
+    give_ao_cart_character_space = 'XXYY'
+   elseif(ao_cart_power(i,2) == 2 .and. ao_cart_power(i,3) == 2)then
+    give_ao_cart_character_space = 'YYZZ'
+   elseif(ao_cart_power(i,1) == 2 .and. ao_cart_power(i,2) == 1 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'XXYZ'
+   elseif(ao_cart_power(i,2) == 2 .and. ao_cart_power(i,1) == 1 .and. ao_cart_power(i,3) == 1)then
+    give_ao_cart_character_space = 'YYXZ'
+   elseif(ao_cart_power(i,3) == 2 .and. ao_cart_power(i,1) == 1 .and. ao_cart_power(i,2) == 1)then
+    give_ao_cart_character_space = 'ZZXY'
+   endif
+  endif
+   ao_cart_l_char_space(i) = give_ao_cart_character_space
+ enddo
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [ logical, use_pw ]
+
+  implicit none
+
+  logical :: exist
+
+  use_pw = .false.
+
+  call ezfio_has_ao_cart_basis_ao_cart_expo_pw(exist)
+  if(exist) then
+    PROVIDE ao_cart_expo_pw_ord_transp
+    if(maxval(dabs(ao_cart_expo_pw_ord_transp(4,:,:))) .gt. 1d-15) use_pw = .true.
+  endif
+
+END_PROVIDER
+

src/ao_cart_basis/aos_in_r.irp.f

@@ -0,0 +1,360 @@
+
+! ---
+
+double precision function ao_cart_value(i, r)
+
+  BEGIN_DOC
+  ! Returns the value of the i-th ao at point $\textbf{r}$
+  END_DOC
+
+  implicit none
+  integer,          intent(in) :: i
+  double precision, intent(in) :: r(3)
+
+  integer          :: m, num_ao
+  integer          :: power_ao(3)
+  double precision :: center_ao(3)
+  double precision :: beta
+  double precision :: accu, dx, dy, dz, r2
+
+  num_ao = ao_cart_nucl(i)
+  power_ao(1:3) = ao_cart_power(i,1:3)
+  center_ao(1:3) = nucl_coord(num_ao,1:3)
+  dx = r(1) - center_ao(1)
+  dy = r(2) - center_ao(2)
+  dz = r(3) - center_ao(3)
+  r2 = dx*dx + dy*dy + dz*dz
+  dx = dx**power_ao(1)
+  dy = dy**power_ao(2)
+  dz = dz**power_ao(3)
+ 
+  accu = 0.d0
+  do m = 1, ao_cart_prim_num(i)
+    beta = ao_cart_expo_ordered_transp(m,i)
+    accu += ao_cart_coef_normalized_ordered_transp(m,i) * dexp(-beta*r2)
+  enddo
+  ao_cart_value = accu * dx * dy * dz
+
+end
+
+
+double precision function primitive_value(i, j, r)
+
+  BEGIN_DOC
+  ! Returns the value of the j-th primitive of the i-th |AO| at point $\textbf{r}
+  ! **without the coefficient**
+  END_DOC
+
+  implicit none
+  integer,          intent(in) :: i, j
+  double precision, intent(in) :: r(3)
+ 
+  integer          :: m, num_ao
+  integer          :: power_ao(3)
+  double precision :: center_ao(3)
+  double precision :: beta
+  double precision :: accu, dx, dy, dz, r2
+
+  num_ao = ao_cart_nucl(i)
+  power_ao(1:3)= ao_cart_power(i,1:3)
+  center_ao(1:3) = nucl_coord(num_ao,1:3)
+  dx = r(1) - center_ao(1)
+  dy = r(2) - center_ao(2)
+  dz = r(3) - center_ao(3)
+  r2 = dx*dx + dy*dy + dz*dz
+  dx = dx**power_ao(1)
+  dy = dy**power_ao(2)
+  dz = dz**power_ao(3)
+ 
+  accu = 0.d0
+  m = j
+  beta = ao_cart_expo_ordered_transp(m,i)
+  accu += dexp(-beta*r2)
+  primitive_value = accu * dx * dy * dz
+
+end
+
+! ---
+
+subroutine give_all_aos_at_r(r, tmp_array)
+
+  BEGIN_dOC
+  !
+  ! input  : r == r(1) = x and so on
+  !
+  ! output : tmp_array(i) = aos(i) evaluated in $\textbf{r}$
+  !
+  END_DOC
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: tmp_array(ao_cart_num)
+  integer                       :: p_ao(3)
+  integer                       :: i, j, k, l, m
+  double precision              :: dx, dy, dz, r2
+  double precision              :: dx2, dy2, dz2
+  double precision              :: c_ao(3)
+  double precision              :: beta
+
+  do i = 1, nucl_num
+
+    c_ao(1:3) = nucl_coord(i,1:3)
+    dx = r(1) - c_ao(1)
+    dy = r(2) - c_ao(2)
+    dz = r(3) - c_ao(3)
+    r2 = dx*dx + dy*dy + dz*dz
+
+    do j = 1, Nucl_N_Aos(i)
+
+      k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
+      p_ao(1:3) = ao_cart_power_ordered_transp_per_nucl(1:3,j,i)
+      dx2 = dx**p_ao(1)
+      dy2 = dy**p_ao(2)
+      dz2 = dz**p_ao(3)
+
+      tmp_array(k) = 0.d0
+      do l = 1, ao_cart_prim_num(k)
+        beta = ao_cart_expo_ordered_transp_per_nucl(l,j,i)
+        if(beta*r2.gt.50.d0) cycle
+
+        tmp_array(k) += ao_cart_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
+      enddo
+
+      tmp_array(k) = tmp_array(k) * dx2 * dy2 * dz2
+    enddo
+  enddo
+
+  return
+end
+
+! ---
+
+subroutine give_all_aos_and_grad_at_r(r, aos_array, aos_grad_array)
+
+  BEGIN_DOC
+  !
+  ! input : r(1) ==> r(1) = x, r(2) = y, r(3) = z
+  !
+  ! output : 
+  !
+  ! * aos_array(i) = ao(i) evaluated at ro
+  ! * aos_grad_array(1,i) = gradient X of the ao(i) evaluated at $\textbf{r}$
+  !
+  END_DOC
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: aos_array(ao_cart_num)
+  double precision, intent(out) :: aos_grad_array(3,ao_cart_num)
+
+  integer                       :: power_ao(3)
+  integer                       :: i, j, k, l, m
+  double precision              :: dx, dy, dz, r2
+  double precision              :: dx1, dy1, dz1
+  double precision              :: dx2, dy2, dz2
+  double precision              :: center_ao(3)
+  double precision              :: beta, accu_1, accu_2, contrib
+
+  do i = 1, nucl_num
+
+    center_ao(1:3) = nucl_coord(i,1:3)
+
+    dx = r(1) - center_ao(1)
+    dy = r(2) - center_ao(2)
+    dz = r(3) - center_ao(3)
+    r2 = dx*dx + dy*dy + dz*dz
+
+    do j = 1, Nucl_N_Aos(i)
+
+      k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
+
+      aos_array(k) = 0.d0
+      aos_grad_array(1,k) = 0.d0
+      aos_grad_array(2,k) = 0.d0
+      aos_grad_array(3,k) = 0.d0
+
+      power_ao(1:3) = ao_cart_power_ordered_transp_per_nucl(1:3,j,i)
+      dx2 = dx**power_ao(1)
+      dy2 = dy**power_ao(2)
+      dz2 = dz**power_ao(3)
+
+      dx1 = 0.d0
+      if(power_ao(1) .ne. 0) then
+        dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
+      endif
+
+      dy1 = 0.d0
+      if(power_ao(2) .ne. 0) then
+        dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
+      endif
+
+      dz1 = 0.d0
+      if(power_ao(3) .ne. 0) then
+        dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
+      endif
+
+      accu_1 = 0.d0
+      accu_2 = 0.d0
+      do l = 1, ao_cart_prim_num(k)
+        beta = ao_cart_expo_ordered_transp_per_nucl(l,j,i)
+        if(beta*r2.gt.50.d0) cycle
+        contrib = ao_cart_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
+        accu_1 += contrib
+        accu_2 += contrib * beta
+      enddo
+
+      aos_array(k) = accu_1 * dx2 * dy2 * dz2
+      aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx  * dy2 * dz2 * accu_2
+      aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy  * dz2 * accu_2
+      aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz  * accu_2
+    enddo
+  enddo
+
+end
+
+! ---
+
+subroutine give_all_aos_and_grad_and_lapl_at_r(r, aos_array, aos_grad_array, aos_lapl_array)
+
+  BEGIN_DOC
+  !
+  ! input  : r(1) ==> r(1) = x, r(2) = y, r(3) = z
+  !
+  ! output :
+  !
+  ! * aos_array(i) = ao(i) evaluated at $\textbf{r}$
+  ! * aos_grad_array(1,i) = $\nabla_x$ of the ao(i) evaluated at $\textbf{r}$
+  !
+  END_DOC
+
+  implicit none
+  double precision, intent(in)  :: r(3)
+  double precision, intent(out) :: aos_array(ao_cart_num)
+  double precision, intent(out) :: aos_grad_array(3,ao_cart_num)
+  double precision, intent(out) :: aos_lapl_array(3,ao_cart_num)
+
+  integer                       :: power_ao(3)
+  integer                       :: i, j, k, l, m
+  double precision              :: dx, dy, dz, r2
+  double precision              :: dx1, dy1, dz1
+  double precision              :: dx2, dy2, dz2
+  double precision              :: dx3, dy3, dz3
+  double precision              :: dx4, dy4, dz4
+  double precision              :: dx5, dy5, dz5
+  double precision              :: center_ao(3)
+  double precision              :: beta, accu_1, accu_2, accu_3, contrib
+
+  do i = 1, nucl_num
+
+    center_ao(1:3) = nucl_coord(i,1:3)
+
+    dx = r(1) - center_ao(1)
+    dy = r(2) - center_ao(2)
+    dz = r(3) - center_ao(3)
+    r2 = dx*dx + dy*dy + dz*dz
+    
+    do j = 1, Nucl_N_Aos(i)
+
+      k = Nucl_Aos_transposed(j,i) ! index of the ao in the ordered format
+
+      aos_array(k) = 0.d0
+      aos_grad_array(1,k) = 0.d0
+      aos_grad_array(2,k) = 0.d0
+      aos_grad_array(3,k) = 0.d0      
+      aos_lapl_array(1,k) = 0.d0
+      aos_lapl_array(2,k) = 0.d0
+      aos_lapl_array(3,k) = 0.d0
+      
+      power_ao(1:3)= ao_cart_power_ordered_transp_per_nucl(1:3,j,i)
+      dx2 = dx**power_ao(1)
+      dy2 = dy**power_ao(2)
+      dz2 = dz**power_ao(3)
+
+      ! ---
+
+      dx1 = 0.d0
+      if(power_ao(1) .ne. 0) then
+        dx1 = dble(power_ao(1)) * dx**(power_ao(1)-1)
+      endif
+
+      dx3 = 0.d0
+      if(power_ao(1) .ge. 2) then
+        dx3 = dble(power_ao(1)) * dble((power_ao(1)-1)) * dx**(power_ao(1)-2)
+      endif
+
+      if(power_ao(1) .ge. 1) then
+        dx4 = dble((2 * power_ao(1) + 1)) * dx**(power_ao(1)) 
+      else
+        dx4 = dble((power_ao(1) + 1)) * dx**(power_ao(1)) 
+      endif
+      
+      dx5 = dx**(power_ao(1)+2)
+  
+      ! ---
+      
+      dy1 = 0.d0
+      if(power_ao(2) .ne. 0) then
+        dy1 = dble(power_ao(2)) * dy**(power_ao(2)-1)
+      endif
+
+      dy3 = 0.d0
+      if(power_ao(2) .ge. 2) then
+        dy3 = dble(power_ao(2)) * dble((power_ao(2)-1))  * dy**(power_ao(2)-2)
+      endif
+      
+      if(power_ao(2) .ge. 1) then
+        dy4 = dble((2 * power_ao(2) + 1)) * dy**(power_ao(2)) 
+      else
+        dy4 = dble((power_ao(2) + 1)) * dy**(power_ao(2)) 
+      endif
+      
+      dy5 = dy**(power_ao(2)+2)
+
+      ! ---
+      
+      dz1 = 0.d0
+      if(power_ao(3) .ne. 0) then
+        dz1 = dble(power_ao(3)) * dz**(power_ao(3)-1)
+      endif
+
+      dz3 = 0.d0
+      if(power_ao(3) .ge. 2) then
+        dz3 = dble(power_ao(3)) * dble((power_ao(3)-1)) * dz**(power_ao(3)-2)
+      endif
+      
+      if(power_ao(3) .ge. 1) then
+        dz4 = dble((2 * power_ao(3) + 1)) * dz**(power_ao(3)) 
+      else
+        dz4 = dble((power_ao(3) + 1)) * dz**(power_ao(3)) 
+      endif
+      
+      dz5 = dz**(power_ao(3)+2)
+      
+      ! ---
+      
+      accu_1 = 0.d0
+      accu_2 = 0.d0
+      accu_3 = 0.d0
+      do l = 1,ao_cart_prim_num(k)
+        beta = ao_cart_expo_ordered_transp_per_nucl(l,j,i)
+        if(beta*r2.gt.50.d0) cycle
+        contrib = ao_cart_coef_normalized_ordered_transp_per_nucl(l,j,i) * dexp(-beta*r2)
+        accu_1 += contrib
+        accu_2 += contrib * beta
+        accu_3 += contrib * beta**2
+      enddo
+
+      aos_array(k) = accu_1 * dx2 * dy2 * dz2
+      aos_grad_array(1,k) = accu_1 * dx1 * dy2 * dz2 - 2.d0 * dx2 * dx  * dy2 * dz2 * accu_2
+      aos_grad_array(2,k) = accu_1 * dx2 * dy1 * dz2 - 2.d0 * dx2 * dy2 * dy  * dz2 * accu_2
+      aos_grad_array(3,k) = accu_1 * dx2 * dy2 * dz1 - 2.d0 * dx2 * dy2 * dz2 * dz  * accu_2
+      aos_lapl_array(1,k) = accu_1 * dx3 * dy2 * dz2 - 2.d0 * dx4 * dy2 * dz2 * accu_2 + 4.d0 * dx5 * dy2 * dz2 * accu_3
+      aos_lapl_array(2,k) = accu_1 * dx2 * dy3 * dz2 - 2.d0 * dx2 * dy4 * dz2 * accu_2 + 4.d0 * dx2 * dy5 * dz2 * accu_3
+      aos_lapl_array(3,k) = accu_1 * dx2 * dy2 * dz3 - 2.d0 * dx2 * dy2 * dz4 * accu_2 + 4.d0 * dx2 * dy2 * dz5 * accu_3
+    enddo
+  enddo
+
+end
+
+! ---
+

src/ao_cart_basis/aos_transp.irp.f

@@ -0,0 +1,68 @@
+
+! ---
+
+BEGIN_PROVIDER [ integer, nucl_aos_transposed, (n_AOs_max,nucl_num)]
+
+  BEGIN_DOC
+  ! List of AOs attached on each atom
+  END_DOC
+ 
+  implicit none
+  integer              :: i
+  integer, allocatable :: nucl_tmp(:)
+
+  allocate(nucl_tmp(nucl_num))
+  nucl_tmp = 0
+  do i = 1, ao_cart_num
+    nucl_tmp(ao_cart_nucl(i)) += 1
+    Nucl_Aos_transposed(nucl_tmp(ao_cart_nucl(i)),ao_cart_nucl(i)) = i
+  enddo
+  deallocate(nucl_tmp)
+
+END_PROVIDER
+
+! ---
+
+BEGIN_PROVIDER [double precision, ao_cart_expo_ordered_transp_per_nucl, (ao_cart_prim_num_max,N_AOs_max,nucl_num) ]
+ implicit none
+ integer :: i,j,k,l
+ do i = 1, nucl_num
+  do j = 1,Nucl_N_Aos(i)
+   k = Nucl_Aos_transposed(j,i)
+   do l = 1, ao_cart_prim_num(k)
+    ao_cart_expo_ordered_transp_per_nucl(l,j,i) = ao_cart_expo_ordered_transp(l,k)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+
+BEGIN_PROVIDER [ integer, ao_cart_power_ordered_transp_per_nucl, (3,N_AOs_max,nucl_num) ]
+ implicit none
+ integer :: i,j,k,l
+ do i = 1, nucl_num
+  do j = 1,Nucl_N_Aos(i)
+   k = Nucl_Aos_transposed(j,i)
+   do l = 1, 3
+    ao_cart_power_ordered_transp_per_nucl(l,j,i) = ao_cart_power(k,l)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+
+BEGIN_PROVIDER [ double precision, ao_cart_coef_normalized_ordered_transp_per_nucl, (ao_cart_prim_num_max,N_AOs_max,nucl_num) ]
+ implicit none
+ integer :: i,j,k,l
+ do i = 1, nucl_num
+  do j = 1,Nucl_N_Aos(i)
+   k = Nucl_Aos_transposed(j,i)
+   do l = 1, ao_cart_prim_num(k)
+    ao_cart_coef_normalized_ordered_transp_per_nucl(l,j,i) = ao_cart_coef_normalized_ordered_transp(l,k)
+   enddo
+  enddo
+ enddo
+
+END_PROVIDER
+

src/ao_cart_basis/cgtos.irp.f

@@ -0,0 +1,201 @@
+
+BEGIN_PROVIDER [logical, use_cgtos]
+
+  implicit none
+
+  BEGIN_DOC
+  ! If true, use cgtos for AO integrals
+  END_DOC
+
+  logical                        :: has
+  PROVIDE ezfio_filename
+  use_cgtos = .False.
+  if (mpi_master) then
+    call ezfio_has_ao_basis_use_cgtos(has)
+    if (has) then
+!      write(6,'(A)') '.. >>>>> [ IO READ: use_cgtos ] <<<<< ..'
+      call ezfio_get_ao_basis_use_cgtos(use_cgtos)
+    else
+      call ezfio_set_ao_basis_use_cgtos(use_cgtos)
+    endif
+  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( use_cgtos, 1, MPI_LOGICAL, 0, MPI_COMM_WORLD, ierr)
+    if (ierr /= MPI_SUCCESS) then
+      stop 'Unable to read use_cgtos with MPI'
+    endif
+  IRP_ENDIF
+
+!  call write_time(6)
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [complex*16, ao_expo_cgtos_ord_transp, (ao_prim_num_max, ao_num)]
+&BEGIN_PROVIDER [double precision, ao_expo_pw_ord_transp, (4, ao_prim_num_max, ao_num)]
+&BEGIN_PROVIDER [double precision, ao_expo_phase_ord_transp, (4, ao_prim_num_max, ao_num)]
+
+  implicit none
+
+  integer :: i, j, m
+
+  do j = 1, ao_num
+    do i = 1, ao_prim_num_max
+
+      ao_expo_cgtos_ord_transp(i,j) = ao_expo_cgtos_ord(j,i)
+
+      do m = 1, 4
+        ao_expo_pw_ord_transp(m,i,j) = ao_expo_pw_ord(m,j,i)
+        ao_expo_phase_ord_transp(m,i,j) = ao_expo_phase_ord(m,j,i)
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER
+
+! ---
+
+ BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos_ord, (ao_num, ao_prim_num_max)]
+&BEGIN_PROVIDER [complex*16      , ao_expo_cgtos_ord, (ao_num, ao_prim_num_max)]
+&BEGIN_PROVIDER [double precision, ao_expo_pw_ord, (4, ao_num, ao_prim_num_max)]
+&BEGIN_PROVIDER [double precision, ao_expo_phase_ord, (4, ao_num, ao_prim_num_max)]
+
+  implicit none
+
+  integer          :: i, j, m
+  integer          :: iorder(ao_prim_num_max)
+  double precision :: d(ao_prim_num_max,11)
+
+  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_cgtos(i,j)
+      d(j,3) = ao_expo_im(i,j)
+
+      do m = 1, 3
+        d(j,3+m) = ao_expo_pw(m,i,j)
+      enddo
+      d(j,7) = d(j,4) * d(j,4) + d(j,5) * d(j,5) + d(j,6) * d(j,6)
+
+      do m = 1, 3
+        d(j,7+m) = ao_expo_phase(m,i,j)
+      enddo
+      d(j,11) = d(j,8) + d(j,9) + d(j,10)
+    enddo
+
+    call dsort(d(1,1), iorder, ao_prim_num(i))
+    do j = 2, 11
+      call dset_order(d(1,j), iorder, ao_prim_num(i))
+    enddo
+
+    do j = 1, ao_prim_num(i)
+      ao_expo_cgtos_ord     (i,j) = d(j,1) + (0.d0, 1.d0) * d(j,3)
+      ao_coef_norm_cgtos_ord(i,j) = d(j,2)
+
+      do m = 1, 4
+        ao_expo_pw_ord(m,i,j) = d(j,3+m)
+        ao_expo_phase_ord(m,i,j) = d(j,7+m)
+      enddo
+    enddo
+  enddo
+
+END_PROVIDER
+
+
+
+! ---
+
+BEGIN_PROVIDER [double precision, ao_coef_cgtos_norm_ord_transp, (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_cgtos_norm_ord_transp(i,j) = ao_coef_norm_cgtos_ord(j,i)
+    enddo
+  enddo
+
+END_PROVIDER
+
+
+! ---
+
+BEGIN_PROVIDER [double precision, ao_coef_norm_cgtos, (ao_num, ao_prim_num_max)]
+
+  implicit none
+
+  integer          :: i, j, ii, m, powA(3), nz
+  double precision :: norm
+  double precision :: kA2, phiA
+  complex*16       :: expo, expo_inv, C_Ae(3), C_Ap(3)
+  complex*16       :: overlap_x, overlap_y, overlap_z
+  complex*16       :: integ1, integ2, C1, C2
+
+  nz = 100
+
+  ao_coef_norm_cgtos = 0.d0
+
+  do i = 1, ao_num
+
+    ii = ao_nucl(i)
+    powA(1) = ao_power(i,1)
+    powA(2) = ao_power(i,2)
+    powA(3) = ao_power(i,3)
+ 
+    if(primitives_normalized) then
+
+      ! Normalization of the primitives
+      do j = 1, ao_prim_num(i)
+
+        expo = ao_expo(i,j) + (0.d0, 1.d0) * ao_expo_im(i,j)
+        expo_inv = (1.d0, 0.d0) / expo
+        do m = 1, 3
+          C_Ap(m) = nucl_coord(ii,m)
+          C_Ae(m) = nucl_coord(ii,m) - (0.d0, 0.5d0) * expo_inv * ao_expo_pw(m,i,j)
+        enddo
+        phiA = ao_expo_phase(1,i,j) + ao_expo_phase(2,i,j) + ao_expo_phase(3,i,j)
+        KA2 = ao_expo_pw(1,i,j) * ao_expo_pw(1,i,j) &
+            + ao_expo_pw(2,i,j) * ao_expo_pw(2,i,j) &
+            + ao_expo_pw(3,i,j) * ao_expo_pw(3,i,j)
+
+        C1 = zexp(-(0.d0, 2.d0) * phiA - 0.5d0 * expo_inv * KA2)
+        C2 = zexp(-(0.5d0, 0.d0) * real(expo_inv) * KA2)
+
+        call overlap_cgaussian_xyz(C_Ae, C_Ae, expo, expo, powA, powA, &
+                                   C_Ap, C_Ap, overlap_x, overlap_y, overlap_z, integ1, nz)
+
+        call overlap_cgaussian_xyz(conjg(C_Ae), C_Ae, conjg(expo), expo, powA, powA, &
+                                   conjg(C_Ap), C_Ap, overlap_x, overlap_y, overlap_z, integ2, nz)
+
+        norm = 2.d0 * real(C1 * integ1 + C2 * integ2)
+
+        !ao_coef_norm_cgtos(i,j) = 1.d0 / dsqrt(norm)
+        ao_coef_norm_cgtos(i,j) = ao_coef(i,j) / dsqrt(norm)
+      enddo
+
+    else
+
+      do j = 1, ao_prim_num(i)
+        ao_coef_norm_cgtos(i,j) = ao_coef(i,j)
+      enddo
+
+    endif ! primitives_normalized
+
+  enddo
+
+END_PROVIDER
+
+

src/ao_cart_basis/dimensions_integrals.irp.f

@@ -0,0 +1,19 @@
+ BEGIN_PROVIDER [ integer, n_pt_max_integrals ]
+&BEGIN_PROVIDER [ integer, n_pt_max_i_x]
+ implicit none
+ BEGIN_DOC
+! Number of points used in the numerical integrations.
+ END_DOC
+ integer :: n_pt_sup
+ integer :: prim_power_l_max
+ include 'utils/constants.include.F'
+ prim_power_l_max = maxval(ao_power)
+ n_pt_max_integrals = 24 * prim_power_l_max + 4
+ n_pt_max_i_x = 8 * prim_power_l_max
+ ASSERT (n_pt_max_i_x-1 <= max_dim)
+ if (n_pt_max_i_x-1 > max_dim) then
+   print *, 'Increase max_dim in utils/constants.include.F to ', n_pt_max_i_x-1
+   stop 1
+ endif
+END_PROVIDER
+

src/ao_cart_basis/spherical_to_cartesian.irp.f

@@ -0,0 +1,807 @@
+! Spherical to cartesian transformation matrix obtained with
+! Horton (http://theochem.github.com/horton/, 2015)
+
+! First index is the index of the cartesian AO, obtained by ao_power_index
+! Second index is the index of the spherical AO
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_0, (1,1) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_0, (1) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=0
+  END_DOC
+  cart_to_sphe_0 = 0.d0
+
+  cart_to_sphe_0 ( 1, 1) = 1.0d0
+  cart_to_sphe_norm_0 (1) = 1.d0
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_1, (3,3) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_1, (3) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=1
+  END_DOC
+  cart_to_sphe_1 = 0.d0
+
+  cart_to_sphe_1 ( 3, 1) = 1.0d0
+  cart_to_sphe_1 ( 1, 2) = 1.0d0
+  cart_to_sphe_1 ( 2, 3) = 1.0d0
+  cart_to_sphe_norm_1 (1) = 1.d0
+  cart_to_sphe_norm_1 (2) = 1.d0
+  cart_to_sphe_norm_1 (3) = 1.d0
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_2, (6,5) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_2, (6) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=2
+  END_DOC
+  cart_to_sphe_2 = 0.d0
+
+  cart_to_sphe_2 ( 1, 1) = -0.5d0
+  cart_to_sphe_2 ( 4, 1) = -0.5d0
+  cart_to_sphe_2 ( 6, 1) = 1.0d0
+  cart_to_sphe_2 ( 3, 2) = 1.0d0
+  cart_to_sphe_2 ( 5, 3) = 1.0d0
+  cart_to_sphe_2 ( 1, 4) = 0.86602540378443864676d0
+  cart_to_sphe_2 ( 4, 4) = -0.86602540378443864676d0
+  cart_to_sphe_2 ( 2, 5) = 1.0d0
+
+  cart_to_sphe_norm_2 = (/ 1.0d0, 1.7320508075688772d0, 1.7320508075688772d0, 1.0d0, &
+                           1.7320508075688772d0, 1.0d0 /)
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_3, (10,7) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_3, (10) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=3
+  END_DOC
+  cart_to_sphe_3 = 0.d0
+
+  cart_to_sphe_3 ( 3, 1) = -0.67082039324993690892d0
+  cart_to_sphe_3 ( 8, 1) = -0.67082039324993690892d0
+  cart_to_sphe_3 (10, 1) = 1.0d0
+  cart_to_sphe_3 ( 1, 2) = -0.61237243569579452455d0
+  cart_to_sphe_3 ( 4, 2) = -0.27386127875258305673d0
+  cart_to_sphe_3 ( 6, 2) = 1.0954451150103322269d0
+  cart_to_sphe_3 ( 2, 3) = -0.27386127875258305673d0
+  cart_to_sphe_3 ( 7, 3) = -0.61237243569579452455d0
+  cart_to_sphe_3 ( 9, 3) = 1.0954451150103322269d0
+  cart_to_sphe_3 ( 3, 4) = 0.86602540378443864676d0
+  cart_to_sphe_3 ( 8, 4) = -0.86602540378443864676d0
+  cart_to_sphe_3 ( 5, 5) = 1.0d0
+  cart_to_sphe_3 ( 1, 6) = 0.790569415042094833d0
+  cart_to_sphe_3 ( 4, 6) = -1.0606601717798212866d0
+  cart_to_sphe_3 ( 2, 7) = 1.0606601717798212866d0
+  cart_to_sphe_3 ( 7, 7) = -0.790569415042094833d0
+
+  cart_to_sphe_norm_3 = (/ 1.0d0, 2.23606797749979d0, 2.23606797749979d0, &
+  2.23606797749979d0, 3.872983346207417d0, 2.23606797749979d0, 1.0d0, 2.23606797749979d0, &
+  2.23606797749979d0, 1.d00 /)
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_4, (15,9) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_4, (15) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=4
+  END_DOC
+  cart_to_sphe_4 = 0.d0
+
+  cart_to_sphe_4 ( 1, 1) = 0.375d0
+  cart_to_sphe_4 ( 4, 1) = 0.21957751641341996535d0
+  cart_to_sphe_4 ( 6, 1) = -0.87831006565367986142d0
+  cart_to_sphe_4 (11, 1) = 0.375d0
+  cart_to_sphe_4 (13, 1) = -0.87831006565367986142d0
+  cart_to_sphe_4 (15, 1) = 1.0d0
+  cart_to_sphe_4 ( 3, 2) = -0.89642145700079522998d0
+  cart_to_sphe_4 ( 8, 2) = -0.40089186286863657703d0
+  cart_to_sphe_4 (10, 2) = 1.19522860933439364d0
+  cart_to_sphe_4 ( 5, 3) = -0.40089186286863657703d0
+  cart_to_sphe_4 (12, 3) = -0.89642145700079522998d0
+  cart_to_sphe_4 (14, 3) = 1.19522860933439364d0
+  cart_to_sphe_4 ( 1, 4) = -0.5590169943749474241d0
+  cart_to_sphe_4 ( 6, 4) = 0.9819805060619657157d0
+  cart_to_sphe_4 (11, 4) = 0.5590169943749474241d0
+  cart_to_sphe_4 (13, 4) = -0.9819805060619657157d0
+  cart_to_sphe_4 ( 2, 5) = -0.42257712736425828875d0
+  cart_to_sphe_4 ( 7, 5) = -0.42257712736425828875d0
+  cart_to_sphe_4 ( 9, 5) = 1.1338934190276816816d0
+  cart_to_sphe_4 ( 3, 6) = 0.790569415042094833d0
+  cart_to_sphe_4 ( 8, 6) = -1.0606601717798212866d0
+  cart_to_sphe_4 ( 5, 7) = 1.0606601717798212866d0
+  cart_to_sphe_4 (12, 7) = -0.790569415042094833d0
+  cart_to_sphe_4 ( 1, 8) = 0.73950997288745200532d0
+  cart_to_sphe_4 ( 4, 8) = -1.2990381056766579701d0
+  cart_to_sphe_4 (11, 8) = 0.73950997288745200532d0
+  cart_to_sphe_4 ( 2, 9) = 1.1180339887498948482d0
+  cart_to_sphe_4 ( 7, 9) = -1.1180339887498948482d0
+
+  cart_to_sphe_norm_4 = (/ 1.0d0, 2.6457513110645907d0, 2.6457513110645907d0, &
+  3.4156502553198664d0, 5.916079783099616d0, 3.415650255319866d0, &
+  2.6457513110645907d0, 5.916079783099616d0, 5.916079783099616d0, &
+  2.6457513110645907d0, 1.0d0, 2.6457513110645907d0, 3.415650255319866d0, &
+  2.6457513110645907d0, 1.d00 /)
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_5, (21,11) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_5, (21) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=5
+  END_DOC
+  cart_to_sphe_5 = 0.d0
+
+  cart_to_sphe_5 ( 3, 1) = 0.625d0
+  cart_to_sphe_5 ( 8, 1) = 0.36596252735569994226d0
+  cart_to_sphe_5 (10, 1) = -1.0910894511799619063d0
+  cart_to_sphe_5 (17, 1) = 0.625d0
+  cart_to_sphe_5 (19, 1) = -1.0910894511799619063d0
+  cart_to_sphe_5 (21, 1) = 1.0d0
+  cart_to_sphe_5 ( 1, 2) = 0.48412291827592711065d0
+  cart_to_sphe_5 ( 4, 2) = 0.21128856368212914438d0
+  cart_to_sphe_5 ( 6, 2) = -1.2677313820927748663d0
+  cart_to_sphe_5 (11, 2) = 0.16137430609197570355d0
+  cart_to_sphe_5 (13, 2) = -0.56694670951384084082d0
+  cart_to_sphe_5 (15, 2) = 1.2909944487358056284d0
+  cart_to_sphe_5 ( 2, 3) = 0.16137430609197570355d0
+  cart_to_sphe_5 ( 7, 3) = 0.21128856368212914438d0
+  cart_to_sphe_5 ( 9, 3) = -0.56694670951384084082d0
+  cart_to_sphe_5 (16, 3) = 0.48412291827592711065d0
+  cart_to_sphe_5 (18, 3) = -1.2677313820927748663d0
+  cart_to_sphe_5 (20, 3) = 1.2909944487358056284d0
+  cart_to_sphe_5 ( 3, 4) = -0.85391256382996653194d0
+  cart_to_sphe_5 (10, 4) = 1.1180339887498948482d0
+  cart_to_sphe_5 (17, 4) = 0.85391256382996653194d0
+  cart_to_sphe_5 (19, 4) = -1.1180339887498948482d0
+  cart_to_sphe_5 ( 5, 5) = -0.6454972243679028142d0
+  cart_to_sphe_5 (12, 5) = -0.6454972243679028142d0
+  cart_to_sphe_5 (14, 5) = 1.2909944487358056284d0
+  cart_to_sphe_5 ( 1, 6) = -0.52291251658379721749d0
+  cart_to_sphe_5 ( 4, 6) = 0.22821773229381921394d0
+  cart_to_sphe_5 ( 6, 6) = 0.91287092917527685576d0
+  cart_to_sphe_5 (11, 6) = 0.52291251658379721749d0
+  cart_to_sphe_5 (13, 6) = -1.2247448713915890491d0
+  cart_to_sphe_5 ( 2, 7) = -0.52291251658379721749d0
+  cart_to_sphe_5 ( 7, 7) = -0.22821773229381921394d0
+  cart_to_sphe_5 ( 9, 7) = 1.2247448713915890491d0
+  cart_to_sphe_5 (16, 7) = 0.52291251658379721749d0
+  cart_to_sphe_5 (18, 7) = -0.91287092917527685576d0
+  cart_to_sphe_5 ( 3, 8) = 0.73950997288745200532d0
+  cart_to_sphe_5 ( 8, 8) = -1.2990381056766579701d0
+  cart_to_sphe_5 (17, 8) = 0.73950997288745200532d0
+  cart_to_sphe_5 ( 5, 9) = 1.1180339887498948482d0
+  cart_to_sphe_5 (12, 9) = -1.1180339887498948482d0
+  cart_to_sphe_5 ( 1,10) = 0.7015607600201140098d0
+  cart_to_sphe_5 ( 4,10) = -1.5309310892394863114d0
+  cart_to_sphe_5 (11,10) = 1.169267933366856683d0
+  cart_to_sphe_5 ( 2,11) = 1.169267933366856683d0
+  cart_to_sphe_5 ( 7,11) = -1.5309310892394863114d0
+  cart_to_sphe_5 (16,11) = 0.7015607600201140098d0
+
+  cart_to_sphe_norm_5 = (/ 1.0d0, 3.0d0, 3.0d0, 4.58257569495584d0, &
+  7.937253933193773d0, 4.58257569495584d0, 4.58257569495584d0, &
+  10.246950765959598d0, 10.246950765959598d0, 4.582575694955841d0, 3.0d0, &
+  7.937253933193773d0, 10.246950765959598d0, 7.937253933193773d0, 3.0d0, 1.0d0, &
+  3.0d0, 4.58257569495584d0, 4.582575694955841d0, 3.0d0, 1.d00 /)
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_6, (28,13) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_6, (28) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=6
+  END_DOC
+  cart_to_sphe_6 = 0.d0
+
+  cart_to_sphe_6 ( 1, 1) = -0.3125d0
+  cart_to_sphe_6 ( 4, 1) = -0.16319780245846672329d0
+  cart_to_sphe_6 ( 6, 1) = 0.97918681475080033975d0
+  cart_to_sphe_6 (11, 1) = -0.16319780245846672329d0
+  cart_to_sphe_6 (13, 1) = 0.57335309036732873772d0
+  cart_to_sphe_6 (15, 1) = -1.3055824196677337863d0
+  cart_to_sphe_6 (22, 1) = -0.3125d0
+  cart_to_sphe_6 (24, 1) = 0.97918681475080033975d0
+  cart_to_sphe_6 (26, 1) = -1.3055824196677337863d0
+  cart_to_sphe_6 (28, 1) = 1.0d0
+  cart_to_sphe_6 ( 3, 2) = 0.86356159963469679725d0
+  cart_to_sphe_6 ( 8, 2) = 0.37688918072220452831d0
+  cart_to_sphe_6 (10, 2) = -1.6854996561581052156d0
+  cart_to_sphe_6 (17, 2) = 0.28785386654489893242d0
+  cart_to_sphe_6 (19, 2) = -0.75377836144440905662d0
+  cart_to_sphe_6 (21, 2) = 1.3816985594155148756d0
+  cart_to_sphe_6 ( 5, 3) = 0.28785386654489893242d0
+  cart_to_sphe_6 (12, 3) = 0.37688918072220452831d0
+  cart_to_sphe_6 (14, 3) = -0.75377836144440905662d0
+  cart_to_sphe_6 (23, 3) = 0.86356159963469679725d0
+  cart_to_sphe_6 (25, 3) = -1.6854996561581052156d0
+  cart_to_sphe_6 (27, 3) = 1.3816985594155148756d0
+  cart_to_sphe_6 ( 1, 4) = 0.45285552331841995543d0
+  cart_to_sphe_6 ( 4, 4) = 0.078832027985861408788d0
+  cart_to_sphe_6 ( 6, 4) = -1.2613124477737825406d0
+  cart_to_sphe_6 (11, 4) = -0.078832027985861408788d0
+  cart_to_sphe_6 (15, 4) = 1.2613124477737825406d0
+  cart_to_sphe_6 (22, 4) = -0.45285552331841995543d0
+  cart_to_sphe_6 (24, 4) = 1.2613124477737825406d0
+  cart_to_sphe_6 (26, 4) = -1.2613124477737825406d0
+  cart_to_sphe_6 ( 2, 5) = 0.27308215547040717681d0
+  cart_to_sphe_6 ( 7, 5) = 0.26650089544451304287d0
+  cart_to_sphe_6 ( 9, 5) = -0.95346258924559231545d0
+  cart_to_sphe_6 (16, 5) = 0.27308215547040717681d0
+  cart_to_sphe_6 (18, 5) = -0.95346258924559231545d0
+  cart_to_sphe_6 (20, 5) = 1.4564381625088382763d0
+  cart_to_sphe_6 ( 3, 6) = -0.81924646641122153043d0
+  cart_to_sphe_6 ( 8, 6) = 0.35754847096709711829d0
+  cart_to_sphe_6 (10, 6) = 1.0660035817780521715d0
+  cart_to_sphe_6 (17, 6) = 0.81924646641122153043d0
+  cart_to_sphe_6 (19, 6) = -1.4301938838683884732d0
+  cart_to_sphe_6 ( 5, 7) = -0.81924646641122153043d0
+  cart_to_sphe_6 (12, 7) = -0.35754847096709711829d0
+  cart_to_sphe_6 (14, 7) = 1.4301938838683884732d0
+  cart_to_sphe_6 (23, 7) = 0.81924646641122153043d0
+  cart_to_sphe_6 (25, 7) = -1.0660035817780521715d0
+  cart_to_sphe_6 ( 1, 8) = -0.49607837082461073572d0
+  cart_to_sphe_6 ( 4, 8) = 0.43178079981734839863d0
+  cart_to_sphe_6 ( 6, 8) = 0.86356159963469679725d0
+  cart_to_sphe_6 (11, 8) = 0.43178079981734839863d0
+  cart_to_sphe_6 (13, 8) = -1.5169496905422946941d0
+  cart_to_sphe_6 (22, 8) = -0.49607837082461073572d0
+  cart_to_sphe_6 (24, 8) = 0.86356159963469679725d0
+  cart_to_sphe_6 ( 2, 9) = -0.59829302641309923139d0
+  cart_to_sphe_6 ( 9, 9) = 1.3055824196677337863d0
+  cart_to_sphe_6 (16, 9) = 0.59829302641309923139d0
+  cart_to_sphe_6 (18, 9) = -1.3055824196677337863d0
+  cart_to_sphe_6 ( 3,10) = 0.7015607600201140098d0
+  cart_to_sphe_6 ( 8,10) = -1.5309310892394863114d0
+  cart_to_sphe_6 (17,10) = 1.169267933366856683d0
+  cart_to_sphe_6 ( 5,11) = 1.169267933366856683d0
+  cart_to_sphe_6 (12,11) = -1.5309310892394863114d0
+  cart_to_sphe_6 (23,11) = 0.7015607600201140098d0
+  cart_to_sphe_6 ( 1,12) = 0.67169328938139615748d0
+  cart_to_sphe_6 ( 4,12) = -1.7539019000502850245d0
+  cart_to_sphe_6 (11,12) = 1.7539019000502850245d0
+  cart_to_sphe_6 (22,12) = -0.67169328938139615748d0
+  cart_to_sphe_6 ( 2,13) = 1.2151388809514737933d0
+  cart_to_sphe_6 ( 7,13) = -1.9764235376052370825d0
+  cart_to_sphe_6 (16,13) = 1.2151388809514737933d0
+
+  cart_to_sphe_norm_6 = (/ 1.0d0, 3.3166247903554003d0, 3.3166247903554003d0, &
+  5.744562646538029d0, 9.949874371066201d0, 5.744562646538029d0, &
+  6.797058187186571d0, 15.198684153570666d0, 15.198684153570664d0, &
+  6.797058187186572d0, 5.744562646538029d0, 15.198684153570666d0, &
+  19.621416870348583d0, 15.198684153570666d0, 5.744562646538029d0, &
+  3.3166247903554003d0, 9.949874371066201d0, 15.198684153570664d0, &
+  15.198684153570666d0, 9.9498743710662d0, 3.3166247903554003d0, 1.0d0, &
+  3.3166247903554003d0, 5.744562646538029d0, 6.797058187186572d0, &
+  5.744562646538029d0, 3.3166247903554003d0, 1.d00 /)
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_7, (36,15) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_7, (36) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=7
+  END_DOC
+  cart_to_sphe_7 = 0.d0
+
+  cart_to_sphe_7 ( 3, 1) = -0.60670333962134435221d0
+  cart_to_sphe_7 ( 8, 1) = -0.31684048566533184861d0
+  cart_to_sphe_7 (10, 1) = 1.4169537279434593918d0
+  cart_to_sphe_7 (17, 1) = -0.31684048566533184861d0
+  cart_to_sphe_7 (19, 1) = 0.82968314787883083417d0
+  cart_to_sphe_7 (21, 1) = -1.5208343311935928733d0
+  cart_to_sphe_7 (30, 1) = -0.60670333962134435221d0
+  cart_to_sphe_7 (32, 1) = 1.4169537279434593918d0
+  cart_to_sphe_7 (34, 1) = -1.5208343311935928733d0
+  cart_to_sphe_7 (36, 1) = 1.0d0
+  cart_to_sphe_7 ( 1, 2) = -0.41339864235384227977d0
+  cart_to_sphe_7 ( 4, 2) = -0.17963167078872714852d0
+  cart_to_sphe_7 ( 6, 2) = 1.4370533663098171882d0
+  cart_to_sphe_7 (11, 2) = -0.1338895422651523892d0
+  cart_to_sphe_7 (13, 2) = 0.62718150750531807803d0
+  cart_to_sphe_7 (15, 2) = -2.1422326762424382273d0
+  cart_to_sphe_7 (22, 2) = -0.1146561540164598136d0
+  cart_to_sphe_7 (24, 2) = 0.47901778876993906273d0
+  cart_to_sphe_7 (26, 2) = -0.95803557753987812546d0
+  cart_to_sphe_7 (28, 2) = 1.4675987714106856141d0
+  cart_to_sphe_7 ( 2, 3) = -0.1146561540164598136d0
+  cart_to_sphe_7 ( 7, 3) = -0.1338895422651523892d0
+  cart_to_sphe_7 ( 9, 3) = 0.47901778876993906273d0
+  cart_to_sphe_7 (16, 3) = -0.17963167078872714852d0
+  cart_to_sphe_7 (18, 3) = 0.62718150750531807803d0
+  cart_to_sphe_7 (20, 3) = -0.95803557753987812546d0
+  cart_to_sphe_7 (29, 3) = -0.41339864235384227977d0
+  cart_to_sphe_7 (31, 3) = 1.4370533663098171882d0
+  cart_to_sphe_7 (33, 3) = -2.1422326762424382273d0
+  cart_to_sphe_7 (35, 3) = 1.4675987714106856141d0
+  cart_to_sphe_7 ( 3, 4) = 0.84254721963085980365d0
+  cart_to_sphe_7 ( 8, 4) = 0.14666864502533059662d0
+  cart_to_sphe_7 (10, 4) = -1.7491256557036030854d0
+  cart_to_sphe_7 (17, 4) = -0.14666864502533059662d0
+  cart_to_sphe_7 (21, 4) = 1.4080189922431737275d0
+  cart_to_sphe_7 (30, 4) = -0.84254721963085980365d0
+  cart_to_sphe_7 (32, 4) = 1.7491256557036030854d0
+  cart_to_sphe_7 (34, 4) = -1.4080189922431737275d0
+  cart_to_sphe_7 ( 5, 5) = 0.50807509012231371428d0
+  cart_to_sphe_7 (12, 5) = 0.49583051751369852316d0
+  cart_to_sphe_7 (14, 5) = -1.3222147133698627284d0
+  cart_to_sphe_7 (23, 5) = 0.50807509012231371428d0
+  cart_to_sphe_7 (25, 5) = -1.3222147133698627284d0
+  cart_to_sphe_7 (27, 5) = 1.6258402883914038857d0
+  cart_to_sphe_7 ( 1, 6) = 0.42961647140211000062d0
+  cart_to_sphe_7 ( 4, 6) = -0.062226236090912312563d0
+  cart_to_sphe_7 ( 6, 6) = -1.2445247218182462513d0
+  cart_to_sphe_7 (11, 6) = -0.23190348980538452414d0
+  cart_to_sphe_7 (13, 6) = 0.54315511828342602619d0
+  cart_to_sphe_7 (15, 6) = 1.2368186122953841287d0
+  cart_to_sphe_7 (22, 6) = -0.35746251148251142922d0
+  cart_to_sphe_7 (24, 6) = 1.2445247218182462513d0
+  cart_to_sphe_7 (26, 6) = -1.6593662957576616683d0
+  cart_to_sphe_7 ( 2, 7) = 0.35746251148251142922d0
+  cart_to_sphe_7 ( 7, 7) = 0.23190348980538452414d0
+  cart_to_sphe_7 ( 9, 7) = -1.2445247218182462513d0
+  cart_to_sphe_7 (16, 7) = 0.062226236090912312563d0
+  cart_to_sphe_7 (18, 7) = -0.54315511828342602619d0
+  cart_to_sphe_7 (20, 7) = 1.6593662957576616683d0
+  cart_to_sphe_7 (29, 7) = -0.42961647140211000062d0
+  cart_to_sphe_7 (31, 7) = 1.2445247218182462513d0
+  cart_to_sphe_7 (33, 7) = -1.2368186122953841287d0
+  cart_to_sphe_7 ( 3, 8) = -0.79037935147039945351d0
+  cart_to_sphe_7 ( 8, 8) = 0.6879369240987588816d0
+  cart_to_sphe_7 (10, 8) = 1.025515817677958738d0
+  cart_to_sphe_7 (17, 8) = 0.6879369240987588816d0
+  cart_to_sphe_7 (19, 8) = -1.8014417303072302517d0
+  cart_to_sphe_7 (30, 8) = -0.79037935147039945351d0
+  cart_to_sphe_7 (32, 8) = 1.025515817677958738d0
+  cart_to_sphe_7 ( 5, 9) = -0.95323336395336381126d0
+  cart_to_sphe_7 (14, 9) = 1.5504341823651057024d0
+  cart_to_sphe_7 (23, 9) = 0.95323336395336381126d0
+  cart_to_sphe_7 (25, 9) = -1.5504341823651057024d0
+  cart_to_sphe_7 ( 1,10) = -0.47495887979908323849d0
+  cart_to_sphe_7 ( 4,10) = 0.61914323168888299344d0
+  cart_to_sphe_7 ( 6,10) = 0.82552430891851065792d0
+  cart_to_sphe_7 (11,10) = 0.25637895441948968451d0
+  cart_to_sphe_7 (13,10) = -1.8014417303072302517d0
+  cart_to_sphe_7 (22,10) = -0.65864945955866621126d0
+  cart_to_sphe_7 (24,10) = 1.3758738481975177632d0
+  cart_to_sphe_7 ( 2,11) = -0.65864945955866621126d0
+  cart_to_sphe_7 ( 7,11) = 0.25637895441948968451d0
+  cart_to_sphe_7 ( 9,11) = 1.3758738481975177632d0
+  cart_to_sphe_7 (16,11) = 0.61914323168888299344d0
+  cart_to_sphe_7 (18,11) = -1.8014417303072302517d0
+  cart_to_sphe_7 (29,11) = -0.47495887979908323849d0
+  cart_to_sphe_7 (31,11) = 0.82552430891851065792d0
+  cart_to_sphe_7 ( 3,12) = 0.67169328938139615748d0
+  cart_to_sphe_7 ( 8,12) = -1.7539019000502850245d0
+  cart_to_sphe_7 (17,12) = 1.7539019000502850245d0
+  cart_to_sphe_7 (30,12) = -0.67169328938139615748d0
+  cart_to_sphe_7 ( 5,13) = 1.2151388809514737933d0
+  cart_to_sphe_7 (12,13) = -1.9764235376052370825d0
+  cart_to_sphe_7 (23,13) = 1.2151388809514737933d0
+  cart_to_sphe_7 ( 1,14) = 0.64725984928774934788d0
+  cart_to_sphe_7 ( 4,14) = -1.96875d0
+  cart_to_sphe_7 (11,14) = 2.4456993503903949804d0
+  cart_to_sphe_7 (22,14) = -1.2566230789301937693d0
+  cart_to_sphe_7 ( 2,15) = 1.2566230789301937693d0
+  cart_to_sphe_7 ( 7,15) = -2.4456993503903949804d0
+  cart_to_sphe_7 (16,15) = 1.96875d0
+  cart_to_sphe_7 (29,15) = -0.64725984928774934788d0
+
+  cart_to_sphe_norm_7 = (/ 1.0d0, 3.6055512754639896d0, 3.605551275463989d0, &
+  6.904105059069327d0, 11.958260743101398d0, 6.904105059069326d0, &
+  9.26282894152753d0, 20.712315177207984d0, 20.71231517720798d0, &
+  9.26282894152753d0, 9.26282894152753d0, 24.507141816213494d0, &
+  31.63858403911275d0, 24.507141816213494d0, 9.262828941527529d0, &
+  6.904105059069327d0, 20.712315177207984d0, 31.63858403911275d0, &
+  31.63858403911275d0, 20.71231517720798d0, 6.904105059069327d0, &
+  3.6055512754639896d0, 11.958260743101398d0, 20.71231517720798d0, &
+  24.507141816213494d0, 20.71231517720798d0, 11.958260743101398d0, &
+  3.6055512754639896d0, 1.0d0, 3.605551275463989d0, 6.904105059069326d0, &
+  9.26282894152753d0, 9.262828941527529d0, 6.904105059069327d0, &
+  3.6055512754639896d0, 1.d00 /)
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_8, (45,17) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_8, (45) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=8
+  END_DOC
+  cart_to_sphe_8 = 0.d0
+
+  cart_to_sphe_8 ( 1, 1) = 0.2734375d0
+  cart_to_sphe_8 ( 4, 1) = 0.13566299095694674896d0
+  cart_to_sphe_8 ( 6, 1) = -1.0853039276555739917d0
+  cart_to_sphe_8 (11, 1) = 0.12099545906566282998d0
+  cart_to_sphe_8 (13, 1) = -0.56678149117738375672d0
+  cart_to_sphe_8 (15, 1) = 1.9359273450506052797d0
+  cart_to_sphe_8 (22, 1) = 0.13566299095694674896d0
+  cart_to_sphe_8 (24, 1) = -0.56678149117738375672d0
+  cart_to_sphe_8 (26, 1) = 1.1335629823547675134d0
+  cart_to_sphe_8 (28, 1) = -1.7364862842489183867d0
+  cart_to_sphe_8 (37, 1) = 0.2734375d0
+  cart_to_sphe_8 (39, 1) = -1.0853039276555739917d0
+  cart_to_sphe_8 (41, 1) = 1.9359273450506052797d0
+  cart_to_sphe_8 (43, 1) = -1.7364862842489183867d0
+  cart_to_sphe_8 (45, 1) = 1.0d0
+  cart_to_sphe_8 ( 3, 2) = -0.84721510698287244363d0
+  cart_to_sphe_8 ( 8, 2) = -0.36813537731583001376d0
+  cart_to_sphe_8 (10, 2) = 2.1951352762686132731d0
+  cart_to_sphe_8 (17, 2) = -0.27439190953357665914d0
+  cart_to_sphe_8 (19, 2) = 0.95803557753987812546d0
+  cart_to_sphe_8 (21, 2) = -2.6341623315223359277d0
+  cart_to_sphe_8 (30, 2) = -0.23497519304418891392d0
+  cart_to_sphe_8 (32, 2) = 0.73171175875620442437d0
+  cart_to_sphe_8 (34, 2) = -1.178033207410656044d0
+  cart_to_sphe_8 (36, 2) = 1.5491933384829667541d0
+  cart_to_sphe_8 ( 5, 3) = -0.23497519304418891392d0
+  cart_to_sphe_8 (12, 3) = -0.27439190953357665914d0
+  cart_to_sphe_8 (14, 3) = 0.73171175875620442437d0
+  cart_to_sphe_8 (23, 3) = -0.36813537731583001376d0
+  cart_to_sphe_8 (25, 3) = 0.95803557753987812546d0
+  cart_to_sphe_8 (27, 3) = -1.178033207410656044d0
+  cart_to_sphe_8 (38, 3) = -0.84721510698287244363d0
+  cart_to_sphe_8 (40, 3) = 2.1951352762686132731d0
+  cart_to_sphe_8 (42, 3) = -2.6341623315223359277d0
+  cart_to_sphe_8 (44, 3) = 1.5491933384829667541d0
+  cart_to_sphe_8 ( 1, 4) = -0.39218438743784791311d0
+  cart_to_sphe_8 ( 4, 4) = -0.0972889728117695298d0
+  cart_to_sphe_8 ( 6, 4) = 1.459334592176542947d0
+  cart_to_sphe_8 (13, 4) = 0.25403754506115685714d0
+  cart_to_sphe_8 (15, 4) = -2.3138757483972597747d0
+  cart_to_sphe_8 (22, 4) = 0.0972889728117695298d0
+  cart_to_sphe_8 (24, 4) = -0.25403754506115685714d0
+  cart_to_sphe_8 (28, 4) = 1.5566235649883124768d0
+  cart_to_sphe_8 (37, 4) = 0.39218438743784791311d0
+  cart_to_sphe_8 (39, 4) = -1.459334592176542947d0
+  cart_to_sphe_8 (41, 4) = 2.3138757483972597747d0
+  cart_to_sphe_8 (43, 4) = -1.5566235649883124768d0
+  cart_to_sphe_8 ( 2, 5) = -0.20252314682524563222d0
+  cart_to_sphe_8 ( 7, 5) = -0.1967766362666553471d0
+  cart_to_sphe_8 ( 9, 5) = 0.8800118701519835797d0
+  cart_to_sphe_8 (16, 5) = -0.1967766362666553471d0
+  cart_to_sphe_8 (18, 5) = 0.85880364827689588344d0
+  cart_to_sphe_8 (20, 5) = -1.7491256557036030854d0
+  cart_to_sphe_8 (29, 5) = -0.20252314682524563222d0
+  cart_to_sphe_8 (31, 5) = 0.8800118701519835797d0
+  cart_to_sphe_8 (33, 5) = -1.7491256557036030854d0
+  cart_to_sphe_8 (35, 5) = 1.7974340685458342478d0
+  cart_to_sphe_8 ( 3, 6) = 0.82265291131801144316d0
+  cart_to_sphe_8 ( 8, 6) = -0.11915417049417047641d0
+  cart_to_sphe_8 (10, 6) = -1.7762455001837659611d0
+  cart_to_sphe_8 (17, 6) = -0.44406137504594149028d0
+  cart_to_sphe_8 (19, 6) = 0.77521709118255285119d0
+  cart_to_sphe_8 (21, 6) = 1.4209964001470127689d0
+  cart_to_sphe_8 (30, 6) = -0.68448859700003543819d0
+  cart_to_sphe_8 (32, 6) = 1.7762455001837659611d0
+  cart_to_sphe_8 (34, 6) = -1.9064667279067276225d0
+  cart_to_sphe_8 ( 5, 7) = 0.68448859700003543819d0
+  cart_to_sphe_8 (12, 7) = 0.44406137504594149028d0
+  cart_to_sphe_8 (14, 7) = -1.7762455001837659611d0
+  cart_to_sphe_8 (23, 7) = 0.11915417049417047641d0
+  cart_to_sphe_8 (25, 7) = -0.77521709118255285119d0
+  cart_to_sphe_8 (27, 7) = 1.9064667279067276225d0
+  cart_to_sphe_8 (38, 7) = -0.82265291131801144316d0
+  cart_to_sphe_8 (40, 7) = 1.7762455001837659611d0
+  cart_to_sphe_8 (42, 7) = -1.4209964001470127689d0
+  cart_to_sphe_8 ( 1, 8) = 0.41132645565900572158d0
+  cart_to_sphe_8 ( 4, 8) = -0.20407507102873838124d0
+  cart_to_sphe_8 ( 6, 8) = -1.2244504261724302874d0
+  cart_to_sphe_8 (11, 8) = -0.3033516698106721761d0
+  cart_to_sphe_8 (13, 8) = 1.0657473001102595767d0
+  cart_to_sphe_8 (15, 8) = 1.2134066792426887044d0
+  cart_to_sphe_8 (22, 8) = -0.20407507102873838124d0
+  cart_to_sphe_8 (24, 8) = 1.0657473001102595767d0
+  cart_to_sphe_8 (26, 8) = -2.1314946002205191534d0
+  cart_to_sphe_8 (37, 8) = 0.41132645565900572158d0
+  cart_to_sphe_8 (39, 8) = -1.2244504261724302874d0
+  cart_to_sphe_8 (41, 8) = 1.2134066792426887044d0
+  cart_to_sphe_8 ( 2, 9) = 0.42481613669916071115d0
+  cart_to_sphe_8 ( 7, 9) = 0.13758738481975177632d0
+  cart_to_sphe_8 ( 9, 9) = -1.4767427774562605828d0
+  cart_to_sphe_8 (16, 9) = -0.13758738481975177632d0
+  cart_to_sphe_8 (20, 9) = 1.8344984642633570176d0
+  cart_to_sphe_8 (29, 9) = -0.42481613669916071115d0
+  cart_to_sphe_8 (31, 9) = 1.4767427774562605828d0
+  cart_to_sphe_8 (33, 9) = -1.8344984642633570176d0
+  cart_to_sphe_8 ( 3,10) = -0.76584818175667166625d0
+  cart_to_sphe_8 ( 8,10) = 0.99833846339806020718d0
+  cart_to_sphe_8 (10,10) = 0.99215674164922147144d0
+  cart_to_sphe_8 (17,10) = 0.41339864235384227977d0
+  cart_to_sphe_8 (19,10) = -2.1650635094610966169d0
+  cart_to_sphe_8 (30,10) = -1.0620403417479017779d0
+  cart_to_sphe_8 (32,10) = 1.6535945694153691191d0
+  cart_to_sphe_8 ( 5,11) = -1.0620403417479017779d0
+  cart_to_sphe_8 (12,11) = 0.41339864235384227977d0
+  cart_to_sphe_8 (14,11) = 1.6535945694153691191d0
+  cart_to_sphe_8 (23,11) = 0.99833846339806020718d0
+  cart_to_sphe_8 (25,11) = -2.1650635094610966169d0
+  cart_to_sphe_8 (38,11) = -0.76584818175667166625d0
+  cart_to_sphe_8 (40,11) = 0.99215674164922147144d0
+  cart_to_sphe_8 ( 1,12) = -0.45768182862115030664d0
+  cart_to_sphe_8 ( 4,12) = 0.79475821795059156217d0
+  cart_to_sphe_8 ( 6,12) = 0.79475821795059156217d0
+  cart_to_sphe_8 (13,12) = -2.0752447144854989366d0
+  cart_to_sphe_8 (22,12) = -0.79475821795059156217d0
+  cart_to_sphe_8 (24,12) = 2.0752447144854989366d0
+  cart_to_sphe_8 (37,12) = 0.45768182862115030664d0
+  cart_to_sphe_8 (39,12) = -0.79475821795059156217d0
+  cart_to_sphe_8 ( 2,13) = -0.70903764004458888811d0
+  cart_to_sphe_8 ( 7,13) = 0.53582588123382020898d0
+  cart_to_sphe_8 ( 9,13) = 1.4377717134510610478d0
+  cart_to_sphe_8 (16,13) = 0.53582588123382020898d0
+  cart_to_sphe_8 (18,13) = -2.338535866733713366d0
+  cart_to_sphe_8 (29,13) = -0.70903764004458888811d0
+  cart_to_sphe_8 (31,13) = 1.4377717134510610478d0
+  cart_to_sphe_8 ( 3,14) = 0.64725984928774934788d0
+  cart_to_sphe_8 ( 8,14) = -1.96875d0
+  cart_to_sphe_8 (17,14) = 2.4456993503903949804d0
+  cart_to_sphe_8 (30,14) = -1.2566230789301937693d0
+  cart_to_sphe_8 ( 5,15) = 1.2566230789301937693d0
+  cart_to_sphe_8 (12,15) = -2.4456993503903949804d0
+  cart_to_sphe_8 (23,15) = 1.96875d0
+  cart_to_sphe_8 (38,15) = -0.64725984928774934788d0
+  cart_to_sphe_8 ( 1,16) = 0.626706654240043952d0
+  cart_to_sphe_8 ( 4,16) = -2.176535018670731151d0
+  cart_to_sphe_8 (11,16) = 3.2353561313826025233d0
+  cart_to_sphe_8 (22,16) = -2.176535018670731151d0
+  cart_to_sphe_8 (37,16) = 0.626706654240043952d0
+  cart_to_sphe_8 ( 2,17) = 1.2945196985754986958d0
+  cart_to_sphe_8 ( 7,17) = -2.9348392204684739765d0
+  cart_to_sphe_8 (16,17) = 2.9348392204684739765d0
+  cart_to_sphe_8 (29,17) = -1.2945196985754986958d0
+
+  cart_to_sphe_norm_8 = (/ 1.0d0, 3.872983346207417d0, 3.872983346207417d0, &
+  8.062257748298551d0, 13.964240043768942d0, 8.06225774829855d0, &
+  11.958260743101398d0, 26.739483914241877d0, 26.739483914241877d0, &
+  11.958260743101398d0, 13.55939315961975d0, 35.874782229304195d0, &
+  46.31414470763765d0, 35.874782229304195d0, 13.55939315961975d0, &
+  11.958260743101398d0, 35.874782229304195d0, 54.79963503528103d0, &
+  54.79963503528103d0, 35.874782229304195d0, 11.958260743101398d0, &
+  8.062257748298551d0, 26.739483914241877d0, 46.31414470763765d0, &
+  54.79963503528103d0, 46.314144707637645d0, 26.739483914241877d0, &
+  8.06225774829855d0, 3.872983346207417d0, 13.964240043768942d0, &
+  26.739483914241877d0, 35.874782229304195d0, 35.874782229304195d0, &
+  26.739483914241877d0, 13.96424004376894d0, 3.8729833462074166d0, 1.0d0, &
+  3.872983346207417d0, 8.06225774829855d0, 11.958260743101398d0, &
+  13.55939315961975d0, 11.958260743101398d0, 8.06225774829855d0, &
+  3.8729833462074166d0, 1.d0 /)
+
+END_PROVIDER
+
+
+ BEGIN_PROVIDER [ double precision, cart_to_sphe_9, (55,19) ]
+&BEGIN_PROVIDER [ double precision, cart_to_sphe_norm_9, (55) ]
+  implicit none
+  BEGIN_DOC
+! Spherical -> Cartesian Transformation matrix for l=9
+  END_DOC
+  cart_to_sphe_9 = 0.d0
+
+  cart_to_sphe_9 ( 3, 1) = 0.59686501473785067702d0
+  cart_to_sphe_9 ( 8, 1) = 0.29612797475437320937d0
+  cart_to_sphe_9 (10, 1) = -1.7657660842403202261d0
+  cart_to_sphe_9 (17, 1) = 0.26411138361943717788d0
+  cart_to_sphe_9 (19, 1) = -0.92214126273187869253d0
+  cart_to_sphe_9 (21, 1) = 2.5354692827465969076d0
+  cart_to_sphe_9 (30, 1) = 0.29612797475437320937d0
+  cart_to_sphe_9 (32, 1) = -0.92214126273187869253d0
+  cart_to_sphe_9 (34, 1) = 1.4846187947947014119d0
+  cart_to_sphe_9 (36, 1) = -1.952374120367905548d0
+  cart_to_sphe_9 (47, 1) = 0.59686501473785067702d0
+  cart_to_sphe_9 (49, 1) = -1.7657660842403202261d0
+  cart_to_sphe_9 (51, 1) = 2.5354692827465969076d0
+  cart_to_sphe_9 (53, 1) = -1.952374120367905548d0
+  cart_to_sphe_9 (55, 1) = 1.0d0
+  cart_to_sphe_9 ( 1, 2) = 0.36685490255855924707d0
+  cart_to_sphe_9 ( 4, 2) = 0.15916400393009351387d0
+  cart_to_sphe_9 ( 6, 2) = -1.5916400393009351387d0
+  cart_to_sphe_9 (11, 2) = 0.11811420148091719529d0
+  cart_to_sphe_9 (13, 2) = -0.6916059470489090194d0
+  cart_to_sphe_9 (15, 2) = 3.1497120394911252077d0
+  cart_to_sphe_9 (22, 2) = 0.098709324918124403125d0
+  cart_to_sphe_9 (24, 2) = -0.51549263708149354579d0
+  cart_to_sphe_9 (26, 2) = 1.3746470322173161221d0
+  cart_to_sphe_9 (28, 2) = -3.1586983973799809d0
+  cart_to_sphe_9 (37, 2) = 0.088975383089683195547d0
+  cart_to_sphe_9 (39, 2) = -0.44144152106008005653d0
+  cart_to_sphe_9 (41, 2) = 1.0499040131637084026d0
+  cart_to_sphe_9 (43, 2) = -1.4126128673922561809d0
+  cart_to_sphe_9 (45, 2) = 1.62697843363992129d0
+  cart_to_sphe_9 ( 2, 3) = 0.088975383089683195547d0
+  cart_to_sphe_9 ( 7, 3) = 0.098709324918124403125d0
+  cart_to_sphe_9 ( 9, 3) = -0.44144152106008005653d0
+  cart_to_sphe_9 (16, 3) = 0.11811420148091719529d0
+  cart_to_sphe_9 (18, 3) = -0.51549263708149354579d0
+  cart_to_sphe_9 (20, 3) = 1.0499040131637084026d0
+  cart_to_sphe_9 (29, 3) = 0.15916400393009351387d0
+  cart_to_sphe_9 (31, 3) = -0.6916059470489090194d0
+  cart_to_sphe_9 (33, 3) = 1.3746470322173161221d0
+  cart_to_sphe_9 (35, 3) = -1.4126128673922561809d0
+  cart_to_sphe_9 (46, 3) = 0.36685490255855924707d0
+  cart_to_sphe_9 (48, 3) = -1.5916400393009351387d0
+  cart_to_sphe_9 (50, 3) = 3.1497120394911252077d0
+  cart_to_sphe_9 (52, 3) = -3.1586983973799809d0
+  cart_to_sphe_9 (54, 3) = 1.62697843363992129d0
+  cart_to_sphe_9 ( 3, 4) = -0.83466307816035426155d0
+  cart_to_sphe_9 ( 8, 4) = -0.2070544267420625878d0
+  cart_to_sphe_9 (10, 4) = 2.3149388661875113029d0
+  cart_to_sphe_9 (19, 4) = 0.40297913150666282783d0
+  cart_to_sphe_9 (21, 4) = -2.9546917977869539993d0
+  cart_to_sphe_9 (30, 4) = 0.2070544267420625878d0
+  cart_to_sphe_9 (32, 4) = -0.40297913150666282783d0
+  cart_to_sphe_9 (36, 4) = 1.7063893769835631924d0
+  cart_to_sphe_9 (47, 4) = 0.83466307816035426155d0
+  cart_to_sphe_9 (49, 4) = -2.3149388661875113029d0
+  cart_to_sphe_9 (51, 4) = 2.9546917977869539993d0
+  cart_to_sphe_9 (53, 4) = -1.7063893769835631924d0
+  cart_to_sphe_9 ( 5, 5) = -0.43101816018790287844d0
+  cart_to_sphe_9 (12, 5) = -0.4187881980957120927d0
+  cart_to_sphe_9 (14, 5) = 1.395960660319040309d0
+  cart_to_sphe_9 (23, 5) = -0.4187881980957120927d0
+  cart_to_sphe_9 (25, 5) = 1.3623181102386339839d0
+  cart_to_sphe_9 (27, 5) = -2.2335370565104644944d0
+  cart_to_sphe_9 (38, 5) = -0.43101816018790287844d0
+  cart_to_sphe_9 (40, 5) = 1.395960660319040309d0
+  cart_to_sphe_9 (42, 5) = -2.2335370565104644944d0
+  cart_to_sphe_9 (44, 5) = 1.9703687322875560157d0
+  cart_to_sphe_9 ( 1, 6) = -0.37548796377180986812d0
+  cart_to_sphe_9 ( 6, 6) = 1.4661859659554465543d0
+  cart_to_sphe_9 (11, 6) = 0.12089373945199884835d0
+  cart_to_sphe_9 (13, 6) = -0.21236437647040795145d0
+  cart_to_sphe_9 (15, 6) = -2.417874789039976967d0
+  cart_to_sphe_9 (22, 6) = 0.20206443016189559856d0
+  cart_to_sphe_9 (24, 6) = -0.79143530297864839268d0
+  cart_to_sphe_9 (26, 6) = 1.0552470706381978569d0
+  cart_to_sphe_9 (28, 6) = 1.6165154412951647885d0
+  cart_to_sphe_9 (37, 6) = 0.27320762396104757397d0
+  cart_to_sphe_9 (39, 6) = -1.2199404645272449631d0
+  cart_to_sphe_9 (41, 6) = 2.417874789039976967d0
+  cart_to_sphe_9 (43, 6) = -2.16878304804843549d0
+  cart_to_sphe_9 ( 2, 7) = -0.27320762396104757397d0
+  cart_to_sphe_9 ( 7, 7) = -0.20206443016189559856d0
+  cart_to_sphe_9 ( 9, 7) = 1.2199404645272449631d0
+  cart_to_sphe_9 (16, 7) = -0.12089373945199884835d0
+  cart_to_sphe_9 (18, 7) = 0.79143530297864839268d0
+  cart_to_sphe_9 (20, 7) = -2.417874789039976967d0
+  cart_to_sphe_9 (31, 7) = 0.21236437647040795145d0
+  cart_to_sphe_9 (33, 7) = -1.0552470706381978569d0
+  cart_to_sphe_9 (35, 7) = 2.16878304804843549d0
+  cart_to_sphe_9 (46, 7) = 0.37548796377180986812d0
+  cart_to_sphe_9 (48, 7) = -1.4661859659554465543d0
+  cart_to_sphe_9 (50, 7) = 2.417874789039976967d0
+  cart_to_sphe_9 (52, 7) = -1.6165154412951647885d0
+  cart_to_sphe_9 ( 3, 8) = 0.80430146722719804411d0
+  cart_to_sphe_9 ( 8, 8) = -0.39904527606894581113d0
+  cart_to_sphe_9 (10, 8) = -1.7845847267806657796d0
+  cart_to_sphe_9 (17, 8) = -0.59316922059788797031d0
+  cart_to_sphe_9 (19, 8) = 1.5532816304615888184d0
+  cart_to_sphe_9 (21, 8) = 1.4236061294349311288d0
+  cart_to_sphe_9 (30, 8) = -0.39904527606894581113d0
+  cart_to_sphe_9 (32, 8) = 1.5532816304615888184d0
+  cart_to_sphe_9 (34, 8) = -2.5007351860179508607d0
+  cart_to_sphe_9 (47, 8) = 0.80430146722719804411d0
+  cart_to_sphe_9 (49, 8) = -1.7845847267806657796d0
+  cart_to_sphe_9 (51, 8) = 1.4236061294349311288d0
+  cart_to_sphe_9 ( 5, 9) = 0.83067898344030094085d0
+  cart_to_sphe_9 (12, 9) = 0.26903627024228973454d0
+  cart_to_sphe_9 (14, 9) = -2.1522901619383178764d0
+  cart_to_sphe_9 (23, 9) = -0.26903627024228973454d0
+  cart_to_sphe_9 (27, 9) = 2.1522901619383178764d0
+  cart_to_sphe_9 (38, 9) = -0.83067898344030094085d0
+  cart_to_sphe_9 (40, 9) = 2.1522901619383178764d0
+  cart_to_sphe_9 (42, 9) = -2.1522901619383178764d0
+  cart_to_sphe_9 ( 1,10) = 0.39636409043643194293d0
+  cart_to_sphe_9 ( 4,10) = -0.34393377440500167929d0
+  cart_to_sphe_9 ( 6,10) = -1.2037682104175058775d0
+  cart_to_sphe_9 (11,10) = -0.29776858550677551679d0
+  cart_to_sphe_9 (13,10) = 1.5691988753163563388d0
+  cart_to_sphe_9 (15,10) = 1.1910743420271020672d0
+  cart_to_sphe_9 (24,10) = 0.64978432507844251538d0
+  cart_to_sphe_9 (26,10) = -2.5991373003137700615d0
+  cart_to_sphe_9 (37,10) = 0.48066206207978815025d0
+  cart_to_sphe_9 (39,10) = -1.6693261563207085231d0
+  cart_to_sphe_9 (41,10) = 1.9851239033785034453d0
+  cart_to_sphe_9 ( 2,11) = 0.48066206207978815025d0
+  cart_to_sphe_9 ( 9,11) = -1.6693261563207085231d0
+  cart_to_sphe_9 (16,11) = -0.29776858550677551679d0
+  cart_to_sphe_9 (18,11) = 0.64978432507844251538d0
+  cart_to_sphe_9 (20,11) = 1.9851239033785034453d0
+  cart_to_sphe_9 (29,11) = -0.34393377440500167929d0
+  cart_to_sphe_9 (31,11) = 1.5691988753163563388d0
+  cart_to_sphe_9 (33,11) = -2.5991373003137700615d0
+  cart_to_sphe_9 (46,11) = 0.39636409043643194293d0
+  cart_to_sphe_9 (48,11) = -1.2037682104175058775d0
+  cart_to_sphe_9 (50,11) = 1.1910743420271020672d0
+  cart_to_sphe_9 ( 3,12) = -0.74463846463549402274d0
+  cart_to_sphe_9 ( 8,12) = 1.2930544805637086353d0
+  cart_to_sphe_9 (10,12) = 0.96378590571704436469d0
+  cart_to_sphe_9 (19,12) = -2.5166038696554342464d0
+  cart_to_sphe_9 (30,12) = -1.2930544805637086353d0
+  cart_to_sphe_9 (32,12) = 2.5166038696554342464d0
+  cart_to_sphe_9 (47,12) = 0.74463846463549402274d0
+  cart_to_sphe_9 (49,12) = -0.96378590571704436469d0
+  cart_to_sphe_9 ( 5,13) = -1.1535889489914915606d0
+  cart_to_sphe_9 (12,13) = 0.87177715295353129935d0
+  cart_to_sphe_9 (14,13) = 1.7435543059070625987d0
+  cart_to_sphe_9 (23,13) = 0.87177715295353129935d0
+  cart_to_sphe_9 (25,13) = -2.8358912905407192076d0
+  cart_to_sphe_9 (38,13) = -1.1535889489914915606d0
+  cart_to_sphe_9 (40,13) = 1.7435543059070625987d0
+  cart_to_sphe_9 ( 1,14) = -0.44314852502786805507d0
+  cart_to_sphe_9 ( 4,14) = 0.96132412415957630049d0
+  cart_to_sphe_9 ( 6,14) = 0.76905929932766104039d0
+  cart_to_sphe_9 (11,14) = -0.33291539937855436029d0
+  cart_to_sphe_9 (13,14) = -2.3392235702823930554d0
+  cart_to_sphe_9 (22,14) = -0.83466307816035426155d0
+  cart_to_sphe_9 (24,14) = 2.9059238431784376645d0
+  cart_to_sphe_9 (37,14) = 0.75235513151094117345d0
+  cart_to_sphe_9 (39,14) = -1.4930907048606177933d0
+  cart_to_sphe_9 ( 2,15) = -0.75235513151094117345d0
+  cart_to_sphe_9 ( 7,15) = 0.83466307816035426155d0
+  cart_to_sphe_9 ( 9,15) = 1.4930907048606177933d0
+  cart_to_sphe_9 (16,15) = 0.33291539937855436029d0
+  cart_to_sphe_9 (18,15) = -2.9059238431784376645d0
+  cart_to_sphe_9 (29,15) = -0.96132412415957630049d0
+  cart_to_sphe_9 (31,15) = 2.3392235702823930554d0
+  cart_to_sphe_9 (46,15) = 0.44314852502786805507d0
+  cart_to_sphe_9 (48,15) = -0.76905929932766104039d0
+  cart_to_sphe_9 ( 3,16) = 0.626706654240043952d0
+  cart_to_sphe_9 ( 8,16) = -2.176535018670731151d0
+  cart_to_sphe_9 (17,16) = 3.2353561313826025233d0
+  cart_to_sphe_9 (30,16) = -2.176535018670731151d0
+  cart_to_sphe_9 (47,16) = 0.626706654240043952d0
+  cart_to_sphe_9 ( 5,17) = 1.2945196985754986958d0
+  cart_to_sphe_9 (12,17) = -2.9348392204684739765d0
+  cart_to_sphe_9 (23,17) = 2.9348392204684739765d0
+  cart_to_sphe_9 (38,17) = -1.2945196985754986958d0
+  cart_to_sphe_9 ( 1,18) = 0.60904939217552380708d0
+  cart_to_sphe_9 ( 4,18) = -2.3781845426185916576d0
+  cart_to_sphe_9 (11,18) = 4.1179360680974030877d0
+  cart_to_sphe_9 (22,18) = -3.4414040330583097636d0
+  cart_to_sphe_9 (37,18) = 1.3294455750836041652d0
+  cart_to_sphe_9 ( 2,19) = 1.3294455750836041652d0
+  cart_to_sphe_9 ( 7,19) = -3.4414040330583097636d0
+  cart_to_sphe_9 (16,19) = 4.1179360680974030877d0
+  cart_to_sphe_9 (29,19) = -2.3781845426185916576d0
+  cart_to_sphe_9 (46,19) = 0.60904939217552380708d0
+
+  cart_to_sphe_norm_9 = (/ 1.0d0, 4.1231056256176615d0, 4.1231056256176615d0, &
+  9.219544457292889d0, 15.968719422671313d0, 9.219544457292889d0, &
+  14.86606874731851d0, 33.24154027718933d0, 33.24154027718933d0, &
+  14.866068747318508d0, 18.635603405463275d0, 49.30517214248421d0, &
+  63.652703529910404d0, 49.30517214248421d0, 18.635603405463275d0, &
+  18.635603405463275d0, 55.90681021638982d0, 85.39906322671229d0, &
+  85.39906322671229d0, 55.90681021638983d0, 18.635603405463275d0, &
+  14.86606874731851d0, 49.30517214248421d0, 85.39906322671229d0, &
+  101.04553429023969d0, 85.3990632267123d0, 49.30517214248421d0, &
+  14.866068747318508d0, 9.219544457292889d0, 33.24154027718933d0, &
+  63.652703529910404d0, 85.39906322671229d0, 85.3990632267123d0, &
+  63.65270352991039d0, 33.24154027718933d0, 9.219544457292887d0, &
+  4.1231056256176615d0, 15.968719422671313d0, 33.24154027718933d0, &
+  49.30517214248421d0, 55.90681021638983d0, 49.30517214248421d0, &
+  33.24154027718933d0, 15.968719422671313d0, 4.1231056256176615d0, 1.0d0, &
+  4.1231056256176615d0, 9.219544457292889d0, 14.866068747318508d0, &
+  18.635603405463275d0, 18.635603405463275d0, 14.866068747318508d0, &
+  9.219544457292887d0, 4.1231056256176615d0, 1.d0 /)
+
+END_PROVIDER
+
+
+
+

src/basis/EZFIO.cfg

@@ -84,3 +84,10 @@ type: logical
 doc: If true, normalize the basis functions
 interface: ezfio, provider, ocaml
 default: false
+
+[use_cgtos]
+type: logical
+doc: If true, use cgtos for AO integrals
+interface: ezfio
+default: False
+