! ! MiMiC: A Framework for Multiscale Modeling in Computational Chemistry ! Copyright (C) 2015-2022 The MiMiC Contributors (see CONTRIBUTORS file for details). ! ! This file is part of MiMiC. ! ! MiMiC is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License as ! published by the Free Software Foundation, either version 3 of ! the License, or (at your option) any later version. ! ! MiMiC is distributed in the hope that it will be useful, but ! WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public License ! along with this program. If not, see <http://www.gnu.org/licenses/>. ! !> author: Jógvan Magnus Haugaard Olsen !> This module provides procedures for computing interactions tensors in !> packed form (one-dimensional arrays) and includes several operations !> involving tensors. module mimic_tensors use mimic_precision implicit none private public :: initialize_tensors public :: terminate_tensors public :: contract_tensors public :: interaction_tensor public :: folded_interaction_tensor public :: tensor_element public :: detrace_tensor public :: convert_index public :: tensor_trace public :: tensor_rank public :: tensor_size public :: polytensor_size !> maximum order of derivative integer, parameter :: MAX_ORDER = 42 !> needed for open-ended computation of derivatives of \(\frac{1}{|\mathbf{r}_{a}-\mathbf{r}_{b}|}\) real(dp), dimension(0:MAX_ORDER,0:MAX_ORDER,1:2*MAX_ORDER+1), save :: tensor_coefficient = huge(0.0_dp) !> factorials computed and stored during initialization real(dp), dimension(0:MAX_ORDER), public, save :: factorial = huge(0.0_dp) !> double factorials computed and stored during initialization real(dp), dimension(-1:MAX_ORDER), public, save :: double_factorial = huge(0.0_dp) !> binomial coefficients computed and stores during initialization real(dp), dimension(0:3*MAX_ORDER,0:MAX_ORDER), public, save :: binomial = huge(0.0_dp) !> trinomial coefficients computed and stores during initialization real(dp), dimension(0:MAX_ORDER,0:MAX_ORDER,0:MAX_ORDER), public, save :: trinomial = huge(0.0_dp) contains !> initialize tensor module subroutine initialize_tensors() call compute_tensor_coefficients() call compute_factorials() call compute_double_factorials() call compute_binomial_coefficients() call compute_trinomial_coefficients() end subroutine initialize_tensors subroutine terminate_tensors() end subroutine terminate_tensors !> compute tensor coefficients needed to compute open-ended derivatives of \(\frac{1}{|\mathbf{r}_{a}-\mathbf{r}_{b}|}\) !> (from C. E. Dykstra, J. Comput. Chem., 9 (1988), 476, where C^(n)_ij are stored as tensor_coefficient(j,i,n)) subroutine compute_tensor_coefficients() integer :: i, j, k, n tensor_coefficient(:,:,:) = 0.0_dp tensor_coefficient(0,0,:) = 1.0_dp do n = 1, 2*MAX_ORDER + 1 if (mod(n, 2) == 0) cycle do i = 1, MAX_ORDER k = merge(i, i - 1, (mod(i, 2) == 0)) do j = 0, i if (mod(i + j, 2) /= 0) cycle if (j == 0) then tensor_coefficient(j,i,n) = tensor_coefficient(j+1,i-1,n) else if (j /= i) then tensor_coefficient(j,i,n) = (real(j, dp) + 1) * tensor_coefficient(j+1,i-1,n) tensor_coefficient(j,i,n) = tensor_coefficient(j,i,n) & - (real(n, dp) + real(k, dp)) * tensor_coefficient(j-1,i-1,n) k = k + 2 else if (j == i) then tensor_coefficient(j,i,n) = - (real(n, dp) + real(k, dp)) * tensor_coefficient(j-1,i-1,n) end if end do end do end do end subroutine compute_tensor_coefficients !> compute and store factorials upto 12! subroutine compute_factorials() integer :: i factorial(0) = 1.0_dp do i = 1, MAX_ORDER factorial(i) = factorial(i-1) * real(i, dp) end do end subroutine compute_factorials !> compute and store double factorials upto 12!! subroutine compute_double_factorials() integer :: i double_factorial(-1) = 1.0_dp double_factorial(0) = 1.0_dp double_factorial(1) = 1.0_dp do i = 2, MAX_ORDER double_factorial(i) = double_factorial(i-2) * real(i, dp) end do end subroutine compute_double_factorials !> Calculate binomial coefficients using recursive algorithm by !> Sten Rettrup and Ruben Pauncz, Int. J. Quantum Chem., 60: 91–98 (1996) !> DOI: 10.1002/(SICI)1097-461X(1996)60:1<91::AID-QUA10>3.0.CO;2-A subroutine compute_binomial_coefficients() integer :: i, j integer :: n, k, nk real(dp), dimension(0:3*MAX_ORDER) :: weights do n = 0, 3*MAX_ORDER do k = 0, MAX_ORDER if (k == n .or. k == 0) then binomial(n,k) = 1.0_dp end if weights(:) = 0.0_dp nk = n - k do i = 0, nk weights(i) = 1.0_dp end do do i = 1, k binomial(n,k) = 0.0_dp do j = 0, nk binomial(n,k) = binomial(n,k) + weights(j) weights(j) = binomial(n,k) end do end do end do end do end subroutine compute_binomial_coefficients !> compute trinomial coefficients (\(\frac{(i+j+k)!}{i!+j!+k!}\)) subroutine compute_trinomial_coefficients() integer :: i, j, k do i = 0, MAX_ORDER do j = 0, MAX_ORDER do k = 0, MAX_ORDER trinomial(i,j,k) = binomial(i,i) * binomial(i+j,j) * binomial(i+j+k,k) end do end do end do end subroutine compute_trinomial_coefficients !> contract tensors A and B according to rank of C and return result in C pure subroutine contract_tensors(A, B, C) real(dp), dimension(:), intent(in) :: A real(dp), dimension(:), intent(in) :: B real(dp), dimension(:), intent(inout) :: C integer :: x_a, y_a, z_a integer :: x_b, y_b, z_b integer :: x_c, y_c, z_c integer :: idx_a, idx_b, idx_c integer :: rank_a, rank_b, rank_c integer :: error_code rank_a = tensor_rank(A) rank_b = tensor_rank(B) rank_c = tensor_rank(C) error_code = 0 if (rank_a < rank_b) then ! TODO: return control to caller with error code ! 'MiMiC ERROR: rank(A) must be equal to or larger than rank(B)' error_code = error_code + 1 end if if (rank_c /= rank_a - rank_b) then ! TODO: return control to caller with error code ! 'MiMiC ERROR: rank(C) must be equal to rank(A) - rank(B)' error_code = error_code + 1 end if do x_c = rank_c, 0, -1 do y_c = rank_c - x_c, 0, -1 z_c = rank_c - x_c - y_c idx_c = convert_index(x_c, y_c, z_c) do x_a = rank_a, x_c, -1 do y_a = rank_a - x_a, y_c, -1 z_a = rank_a - x_a - y_a if (z_a < z_c) cycle idx_a = convert_index(x_a, y_a, z_a) x_b = x_a - x_c y_b = y_a - y_c z_b = z_a - z_c idx_b = convert_index(x_b, y_b, z_b) C(idx_c) = C(idx_c) + trinomial(x_b, y_b, z_b) * A(idx_a) * B(idx_b) end do end do end do end do end subroutine contract_tensors !> Compute multipole interaction tensor containing derivatives of \(\frac{1}{r_{ab}} = \frac{1}{|\mathbf{r}_{b}-\mathbf{r}_{a}|}\), !> i.e. the da'th derivative wrt \(\mathbf{r}_{a}\) and db'th derivative wrt \(\mathbf{r}_{b}\). !> Tensor is returned in a one-dimensional array with elements in anticanonical ordering, !> e.g. xx, xy, xz, yy, yz, zz. pure subroutine interaction_tensor(r_ab, T, da, db) !> coordinate \(\mathbf{r}_{ab} = \mathbf{r}_{b}-\mathbf{r}_{a}\) real(dp), dimension(:), intent(in) :: r_ab !> the interaction tensor real(dp), dimension(:), intent(out) :: T !> order of derivative wrt \(\mathbf{r}_{a}\) integer, intent(in) :: da !> order of derivative wrt \(\mathbf{r}_{b}\) integer, intent(in) :: db integer :: i, j, k, l integer :: error_code real(dp) :: sgn error_code = 0 if (size(T) /= (da + db + 1) * (da + db + 2) / 2) then ! 'MiMiC ERROR in tensor: wrong size of T_ba' error_code = error_code + 1 end if ! tensor_elements gives derivatives of 1 / |r_b - r_a| wrt r_b, so therefore the sign ! alternates for each order derivative, and thus we here change sign for equal order ! derivatives with respect to r_a, because it should not change sign sgn = merge(-1.0_dp, 1.0_dp, (mod(da, 2) /= 0)) l = 1 do i = da + db, 0, - 1 do j = da + db - i, 0, - 1 k = da + db - i - j T(l) = sgn * tensor_element(i, j, k, r_ab) l = l + 1 end do end do end subroutine interaction_tensor !> Compute folded multipole interaction tensor containing derivatives of !> \(\frac{1}{r_{ab}} = \frac{1}{|\mathbf{r}_{b}-\mathbf{r}_{a}|}\), !> i.e. the da'th derivative wrt \(\mathbf{r}_{a}\) and db'th derivative wrt \(\mathbf{r}_{b}\). !> Tensor is returned in a one-dimensional array with elements in anticanonical ordering, !> e.g. xx, xy, xz, yy, yz, zz. pure subroutine folded_interaction_tensor(r_ab, T, da, db) !> coordinate \(\mathbf{r}_{ab} = \mathbf{r}_{b}-\mathbf{r}_{a}\) real(dp), dimension(:), intent(in) :: r_ab !> the interaction tensor real(dp), dimension(:), intent(out) :: T !> order of derivative wrt \(\mathbf{r}_{a}\) integer, intent(in) :: da !> order of derivative wrt \(\mathbf{r}_{b}\) integer, intent(in) :: db integer :: i, j, k, l integer :: error_code real(dp) :: sgn error_code = 0 if (size(T) /= (da + db + 1) * (da + db + 2) / 2) then ! 'MiMiC ERROR in tensor: wrong size of T_ba' error_code = error_code + 1 end if ! tensor_elements gives derivatives of 1 / |r_b - r_a| wrt r_b, so therefore the sign ! alternates for each order derivative, and thus we here change sign for equal order ! derivatives with respect to r_a, because it should not change sign sgn = merge(-1.0_dp, 1.0_dp, (mod(da, 2) /= 0)) l = 1 do i = da + db, 0, - 1 do j = da + db - i, 0, - 1 k = da + db - i - j T(l) = sgn * trinomial(i, j, k) * tensor_element(i, j, k, r_ab) l = l + 1 end do end do end subroutine folded_interaction_tensor !> compute element of multipole interaction tensor (C. E. Dykstra, J. Comput. Chem., 9 (1988), 476) !> which in multi-index notation this is !> \(\mathbf{T}^{(i,j,k)} !> = \frac{\mathrm{d}^{i}\mathrm{d}^{j}\mathrm{d}^{k}}{\mathrm{d}x_{a}^{i}\mathrm{d}y_{a}^{j}\mathrm{d}z_{a}^{k}} !> \frac{1}{|\mathbf{r}_{a} - \mathbf{r}_{b}|}\) real(dp) pure function tensor_element(i, j, k, r) !> order of the derivative wrt \(x\)-component of \(\mathbf{r}_{a}\) integer, intent(in) :: i !> order of the derivative wrt \(y\)-component of \(\mathbf{r}_{a}\) integer, intent(in) :: j !> order of the derivative wrt \(z\)-component of \(\mathbf{r}_{a}\) integer, intent(in) :: k !> coordinate vector \(\mathbf{r} = \mathbf{r}_{a}\) - \(\mathbf{r}_{b}\) real(dp), dimension(:), intent(in) :: r integer :: l, m, n, o, p real(dp) :: cl, cm, cn real(dp) :: norm tensor_element = 0.0_dp norm = norm2(r) do l = 0, i cl = tensor_coefficient(l,i,1) * (r(1) / norm)**l o = l + i + 1 do m = 0, j cm = cl * tensor_coefficient(m,j,o) * (r(2) / norm)**m p = o + j + m do n = 0, k cn = cm * tensor_coefficient(n,k,p) * (r(3) / norm)**n tensor_element = tensor_element + cn end do end do end do tensor_element = tensor_element / norm**(i + j + k + 1) end function tensor_element !> convert tensor to traceless tensor pure subroutine detrace_tensor(A) real(dp), dimension(:), intent(inout) :: A integer :: rank integer :: norm integer :: ti integer :: ax, ay, az integer :: bx, by, bz real(dp) :: xn, yn, zn real(dp) :: sgn real(dp) :: divisor real(dp), dimension(:), allocatable :: trace_tensor rank = tensor_rank(A) divisor = double_factorial(2 * rank - 1) allocate(trace_tensor(size(A))) trace_tensor = 0.0_dp do ax = rank, 0, -1 do ay = rank - ax, 0, -1 az = rank - ax - ay ti = convert_index(ax, ay, az) do bx = 0, floor(real(ax)/2.0) xn = some_nomial(ax, bx) do by = 0, floor(real(ay)/2.0) yn = some_nomial(ay, by) do bz = 0, floor(real(az)/2.0) zn = some_nomial(az, bz) norm = bx + by + bz if (norm /= 1) then cycle end if sgn = (-1.0_dp)**(norm) trace_tensor(ti) = trace_tensor(ti) & + sgn * double_factorial(2*rank-2*norm-1) * xn * yn * zn & * tensor_trace(A, ax-2*bx, ay-2*by, az-2*bz) / divisor end do end do end do end do end do A = A + trace_tensor deallocate(trace_tensor) contains real(dp) pure function some_nomial(ai, bi) integer, intent(in) :: ai integer, intent(in) :: bi some_nomial = factorial(ai) / (2.0_dp**(bi) * factorial(bi) * factorial(ai-2*bi)) end function some_nomial end subroutine detrace_tensor !> compute trace of a packed tensor real(dp) pure function tensor_trace(A, ax, ay, az, order) real(dp), dimension(:), intent(in) :: A integer, intent(in) :: ax integer, intent(in) :: ay integer, intent(in) :: az integer, intent(in), optional :: order integer :: ti integer :: bx, by, bz integer :: trace_order if (present(order)) then trace_order = order else trace_order = 1 end if tensor_trace = 0.0_dp do bx = trace_order, 0, -1 do by = trace_order - bx, 0, -1 bz = trace_order - bx - by ti = convert_index(ax+2*bx, ay+2*by, az+2*bz) tensor_trace = tensor_trace + trinomial(bx, by, bz) * A(ti) end do end do end function tensor_trace !> convert from multi-index to packed tensor index integer pure function convert_index(ax, ay, az) !> first (\(x\)) component of multi-index integer, intent(in) :: ax !> second (\(y\)) component of multi-index integer, intent(in) :: ay !> third (\(z\)) component of multi-index integer, intent(in) :: az integer :: norm integer :: bx integer :: by integer :: bz norm = ax + ay + az convert_index = 1 do bx = norm, 0, -1 do by = norm - bx, 0, -1 bz = norm - bx - by if (bx /= ax .or. by /= ay .or. bz /= az) then convert_index = convert_index + 1 else return end if end do end do end function convert_index !> compute rank of tensor integer pure function tensor_rank(A) real(dp), dimension(:), intent(in) :: A tensor_rank = int(0.5_dp * (sqrt(8.0_dp * real(size(A), dp) + 1.0_dp) - 3.0_dp)) end function tensor_rank !> compute size of tensor of given rank integer pure function tensor_size(rank) integer, intent(in) :: rank tensor_size = (rank + 2) * (rank + 1) / 2 end function tensor_size !> compute size of first degree polytensor with tensors up to maxmimum rank integer pure function polytensor_size(max_rank) integer, intent(in) :: max_rank polytensor_size = (max_rank + 3) * (max_rank + 2) * (max_rank + 1) / 6 end function polytensor_size end module mimic_tensors