From 95e8c805b6efb27c70fc001a5093e9d50def6e22 Mon Sep 17 00:00:00 2001 From: eginer Date: Sat, 12 Apr 2025 14:30:24 +0200 Subject: [PATCH] starting to create cart --- external/irpf90 | 2 +- .../localization_fboys_sub.irp.f | 193 ++++ .../localization_pipek_sub.irp.f | 776 ++++++++++++++ .../mo_localization/localization_sub.irp.f | 970 +----------------- src/ao_cart_basis/EZFIO.cfg | 75 ++ src/ao_cart_basis/README.rst | 45 + src/ao_cart_basis/aos.irp.f | 445 ++++++++ src/ao_cart_basis/aos_in_r.irp.f | 360 +++++++ src/ao_cart_basis/aos_transp.irp.f | 68 ++ src/ao_cart_basis/cgtos.irp.f | 201 ++++ src/ao_cart_basis/dimensions_integrals.irp.f | 19 + .../spherical_to_cartesian.irp.f | 807 +++++++++++++++ src/basis/EZFIO.cfg | 7 + 13 files changed, 2998 insertions(+), 970 deletions(-) create mode 100644 plugins/local/mo_localization/localization_fboys_sub.irp.f create mode 100644 plugins/local/mo_localization/localization_pipek_sub.irp.f create mode 100644 src/ao_cart_basis/EZFIO.cfg create mode 100644 src/ao_cart_basis/README.rst create mode 100644 src/ao_cart_basis/aos.irp.f create mode 100644 src/ao_cart_basis/aos_in_r.irp.f create mode 100644 src/ao_cart_basis/aos_transp.irp.f create mode 100644 src/ao_cart_basis/cgtos.irp.f create mode 100644 src/ao_cart_basis/dimensions_integrals.irp.f create mode 100644 src/ao_cart_basis/spherical_to_cartesian.irp.f diff --git a/external/irpf90 b/external/irpf90 index 43160c60..4ab1b175 160000 --- a/external/irpf90 +++ b/external/irpf90 @@ -1 +1 @@ -Subproject commit 43160c60d88d9f61fb97cc0b35477c8eb0df862b +Subproject commit 4ab1b175fc7ed0d96c1912f13dc53579b24157a6 diff --git a/plugins/local/mo_localization/localization_fboys_sub.irp.f b/plugins/local/mo_localization/localization_fboys_sub.irp.f new file mode 100644 index 00000000..23c83171 --- /dev/null +++ b/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}} (^2 + ^2 + ^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 ^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 ^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 + diff --git a/plugins/local/mo_localization/localization_pipek_sub.irp.f b/plugins/local/mo_localization/localization_pipek_sub.irp.f new file mode 100644 index 00000000..9e3bdba1 --- /dev/null +++ b/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 + 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 \left[ - \right] +! \end{align*} + +! \begin{align*} +! = \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( ^2 - \frac{1}{4} \left[ - \right]^2 \right) +! \end{align*} + +! with +! \begin{align*} +! = \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[ \right]^2 +! \end{align*} +! with +! \begin{align*} +! = \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 + 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 + diff --git a/plugins/local/mo_localization/localization_sub.irp.f b/plugins/local/mo_localization/localization_sub.irp.f index f5afed07..7438dc42 100644 --- a/plugins/local/mo_localization/localization_sub.irp.f +++ b/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 - 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 \left[ - \right] -! \end{align*} - -! \begin{align*} -! = \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( ^2 - \frac{1}{4} \left[ - \right]^2 \right) -! \end{align*} - -! with -! \begin{align*} -! = \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[ \right]^2 -! \end{align*} -! with -! \begin{align*} -! = \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}} (^2 + ^2 + ^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 ^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 ^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 - 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 diff --git a/src/ao_cart_basis/EZFIO.cfg b/src/ao_cart_basis/EZFIO.cfg new file mode 100644 index 00000000..644ca375 --- /dev/null +++ b/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 + diff --git a/src/ao_cart_basis/README.rst b/src/ao_cart_basis/README.rst new file mode 100644 index 00000000..0b7e5c3e --- /dev/null +++ b/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``). diff --git a/src/ao_cart_basis/aos.irp.f b/src/ao_cart_basis/aos.irp.f new file mode 100644 index 00000000..e2ede308 --- /dev/null +++ b/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 + diff --git a/src/ao_cart_basis/aos_in_r.irp.f b/src/ao_cart_basis/aos_in_r.irp.f new file mode 100644 index 00000000..2a9b0712 --- /dev/null +++ b/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 + +! --- + diff --git a/src/ao_cart_basis/aos_transp.irp.f b/src/ao_cart_basis/aos_transp.irp.f new file mode 100644 index 00000000..d06b358e --- /dev/null +++ b/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 + diff --git a/src/ao_cart_basis/cgtos.irp.f b/src/ao_cart_basis/cgtos.irp.f new file mode 100644 index 00000000..abf46dcb --- /dev/null +++ b/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 + + diff --git a/src/ao_cart_basis/dimensions_integrals.irp.f b/src/ao_cart_basis/dimensions_integrals.irp.f new file mode 100644 index 00000000..97fd83e1 --- /dev/null +++ b/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 + diff --git a/src/ao_cart_basis/spherical_to_cartesian.irp.f b/src/ao_cart_basis/spherical_to_cartesian.irp.f new file mode 100644 index 00000000..53a60413 --- /dev/null +++ b/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 + + + + diff --git a/src/basis/EZFIO.cfg b/src/basis/EZFIO.cfg index a9ec2c1b..c20440dc 100644 --- a/src/basis/EZFIO.cfg +++ b/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 +