Inter-particle distances

1. Squared distance
- 1.1. qmckl_distance_sq
  - 1.1.1. Performance
2. Distance
- 2.1. qmckl_distance
3. Rescaled Distance
- 3.1. qmckl_distance_rescaled
4. Rescaled Distance Derivatives
- 4.1. qmckl_distance_rescaled_gl

1 Squared distance

1.1 `qmckl_distance_sq`

qmckl_distance_sq computes the matrix of the squared distances between all pairs of points in two sets, one point within each set:

\[ C_{ij} = \sum_{k=1}^3 (A_{k,i}-B_{k,j})^2 \]

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`transa`	`char`	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
`transb`	`char`	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
`m`	`int64_t`	in	Number of points in the first set
`n`	`int64_t`	in	Number of points in the second set
`A`	`double[][lda]`	in	Array containing the \(m \times 3\) matrix \(A\)
`lda`	`int64_t`	in	Leading dimension of array `A`
`B`	`double[][ldb]`	in	Array containing the \(n \times 3\) matrix \(B\)
`ldb`	`int64_t`	in	Leading dimension of array `B`
`C`	`double[n][ldc]`	out	Array containing the \(m \times n\) matrix \(C\)
`ldc`	`int64_t`	in	Leading dimension of array `C`

Requirements:

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least \(3 \times m \times 8\) bytes
B is allocated with at least \(3 \times n \times 8\) bytes
C is allocated with at least \(m \times n \times 8\) bytes

qmckl_exit_code qmckl_distance_sq (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc );

integer function qmckl_distance_sq_f(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  character  , intent(in)  :: transa, transb
  integer*8  , intent(in)  :: m, n
  integer*8  , intent(in)  :: lda
  real*8     , intent(in)  :: A(lda,*)
  integer*8  , intent(in)  :: ldb
  real*8     , intent(in)  :: B(ldb,*)
  integer*8  , intent(in)  :: ldc
  real*8     , intent(out) :: C(ldc,*)

  integer*8 :: i,j
  real*8    :: x, y, z
  integer   :: transab

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (m <= 0_8) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (n <= 0_8) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  if (transa == 'N' .or. transa == 'n') then
     transab = 0
  else if (transa == 'T' .or. transa == 't') then
     transab = 1
  else
     transab = -100
  endif

  if (transb == 'N' .or. transb == 'n') then
     continue
  else if (transb == 'T' .or. transb == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  if (transab < 0) then
     info = QMCKL_INVALID_ARG_1
     return
  endif

  if (iand(transab,1) == 0 .and. LDA < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (iand(transab,1) == 1 .and. LDA < m) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (iand(transab,2) == 2 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_7
     return
  endif


  select case (transab)

  case(0)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(1,j)
           y = A(2,i) - B(2,j)
           z = A(3,i) - B(3,j)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  case(1)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(1,j)
           y = A(i,2) - B(2,j)
           z = A(i,3) - B(3,j)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  case(2)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(j,1)
           y = A(2,i) - B(j,2)
           z = A(3,i) - B(j,3)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  case(3)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(j,1)
           y = A(i,2) - B(j,2)
           z = A(i,3) - B(j,3)
           C(i,j) = x*x + y*y + z*z
        end do
     end do

  end select

end function qmckl_distance_sq_f

1.1.1 Performance

This function is more efficient when A and B are transposed.

2 Distance

2.1 `qmckl_distance`

qmckl_distance computes the matrix of the distances between all pairs of points in two sets, one point within each set:

\[ C_{ij} = \sqrt{\sum_{k=1}^3 (A_{k,i}-B_{k,j})^2} \]

If the input array is normal ('N'), the xyz coordinates are in the leading dimension: [n][3] in C and (3,n) in Fortran.

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`transa`	`char`	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
`transb`	`char`	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
`m`	`int64_t`	in	Number of points in the first set
`n`	`int64_t`	in	Number of points in the second set
`A`	`double[][lda]`	in	Array containing the \(m \times 3\) matrix \(A\)
`lda`	`int64_t`	in	Leading dimension of array `A`
`B`	`double[][ldb]`	in	Array containing the \(n \times 3\) matrix \(B\)
`ldb`	`int64_t`	in	Leading dimension of array `B`
`C`	`double[n][ldc]`	out	Array containing the \(m \times n\) matrix \(C\)
`ldc`	`int64_t`	in	Leading dimension of array `C`

2.1.1 Requirements

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least \(3 \times m \times 8\) bytes
B is allocated with at least \(3 \times n \times 8\) bytes
C is allocated with at least \(m \times n \times 8\) bytes

2.1.2 C header

qmckl_exit_code qmckl_distance (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc );

2.1.3 Source

integer function qmckl_distance_f(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  character  , intent(in)  :: transa, transb
  integer*8  , intent(in)  :: m, n
  integer*8  , intent(in)  :: lda
  real*8     , intent(in)  :: A(lda,*)
  integer*8  , intent(in)  :: ldb
  real*8     , intent(in)  :: B(ldb,*)
  integer*8  , intent(in)  :: ldc
  real*8     , intent(out) :: C(ldc,*)

  integer*8 :: i,j
  real*8    :: x, y, z
  integer   :: transab

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (m <= 0_8) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (n <= 0_8) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  if (transa == 'N' .or. transa == 'n') then
     transab = 0
  else if (transa == 'T' .or. transa == 't') then
     transab = 1
  else
     transab = -100
  endif

  if (transb == 'N' .or. transb == 'n') then
     continue
  else if (transb == 'T' .or. transb == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  if (transab < 0) then
     info = QMCKL_INVALID_ARG_1
     return
  endif

  ! check for LDA
  if (iand(transab,1) == 0 .and. LDA < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (iand(transab,1) == 1 .and. LDA < m) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  ! check for LDB
  if (iand(transab,1) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,1) == 1 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  ! check for LDC
  if (LDC < m) then
     info = QMCKL_INVALID_ARG_11
     return
  endif


  select case (transab)

  case(0)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(1,j)
           y = A(2,i) - B(2,j)
           z = A(3,i) - B(3,j)
           C(i,j) = x*x + y*y + z*z
        end do
        C(:,j) = dsqrt(C(:,j))
     end do

  case(1)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(1,j)
           y = A(i,2) - B(2,j)
           z = A(i,3) - B(3,j)
           C(i,j) = x*x + y*y + z*z
        end do
        C(:,j) = dsqrt(C(:,j))
     end do

  case(2)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(j,1)
           y = A(2,i) - B(j,2)
           z = A(3,i) - B(j,3)
           C(i,j) = x*x + y*y + z*z
        end do
        C(:,j) = dsqrt(C(:,j))
     end do

  case(3)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(j,1)
           y = A(i,2) - B(j,2)
           z = A(i,3) - B(j,3)
           C(i,j) = x*x + y*y + z*z
        end do
        C(:,j) = dsqrt(C(:,j))
     end do

  end select

end function qmckl_distance_f

2.1.4 Performance

This function is more efficient when A and B are transposed.

3 Rescaled Distance

3.1 `qmckl_distance_rescaled`

qmckl_distance_rescaled computes the matrix of the rescaled distances between all pairs of points in two sets, one point within each set:

\[ C_{ij} = \frac{ 1 - e^{-\kappa r_{ij}}}{\kappa} \]

If the input array is normal ('N'), the xyz coordinates are in the leading dimension: [n][3] in C and (3,n) in Fortran.

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`transa`	`char`	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
`transb`	`char`	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
`m`	`int64_t`	in	Number of points in the first set
`n`	`int64_t`	in	Number of points in the second set
`A`	`double[][lda]`	in	Array containing the \(m \times 3\) matrix \(A\)
`lda`	`int64_t`	in	Leading dimension of array `A`
`B`	`double[][ldb]`	in	Array containing the \(n \times 3\) matrix \(B\)
`ldb`	`int64_t`	in	Leading dimension of array `B`
`C`	`double[n][ldc]`	out	Array containing the \(m \times n\) matrix \(C\)
`ldc`	`int64_t`	in	Leading dimension of array `C`
`rescale_factor_kappa`	`double`	in	Factor for calculating rescaled distances

3.1.1 Requirements

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least \(3 \times m \times 8\) bytes
B is allocated with at least \(3 \times n \times 8\) bytes
C is allocated with at least \(m \times n \times 8\) bytes

3.1.2 C header

qmckl_exit_code qmckl_distance_rescaled (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc,
      const double rescale_factor_kappa );

3.1.3 Source

integer function qmckl_distance_rescaled_f(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  character  , intent(in)  :: transa, transb
  integer*8  , intent(in)  :: m, n
  integer*8  , intent(in)  :: lda
  real*8     , intent(in)  :: A(lda,*)
  integer*8  , intent(in)  :: ldb
  real*8     , intent(in)  :: B(ldb,*)
  integer*8  , intent(in)  :: ldc
  real*8     , intent(out) :: C(ldc,*)
  real*8     , intent(in)  :: rescale_factor_kappa

  integer*8 :: i,j
  real*8    :: x, y, z, dist, rescale_factor_kappa_inv
  integer   :: transab

  rescale_factor_kappa_inv = 1.0d0/rescale_factor_kappa;

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (m <= 0_8) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (n <= 0_8) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  if (transa == 'N' .or. transa == 'n') then
     transab = 0
  else if (transa == 'T' .or. transa == 't') then
     transab = 1
  else
     transab = -100
  endif

  if (transb == 'N' .or. transb == 'n') then
     continue
  else if (transb == 'T' .or. transb == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  ! check for LDA
  if (transab < 0) then
     info = QMCKL_INVALID_ARG_1
     return
  endif

  if (iand(transab,1) == 0 .and. LDA < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (iand(transab,1) == 1 .and. LDA < m) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  ! check for LDB
  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 2 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  ! check for LDC
  if (LDC < m) then
     info = QMCKL_INVALID_ARG_11
     return
  endif


  select case (transab)

  case(0)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(1,j)
           y = A(2,i) - B(2,j)
           z = A(3,i) - B(3,j)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        end do
     end do

  case(1)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(1,j)
           y = A(i,2) - B(2,j)
           z = A(i,3) - B(3,j)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        end do
     end do

  case(2)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(j,1)
           y = A(2,i) - B(j,2)
           z = A(3,i) - B(j,3)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        end do
     end do

  case(3)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(j,1)
           y = A(i,2) - B(j,2)
           z = A(i,3) - B(j,3)
           dist = dsqrt(x*x + y*y + z*z)
           C(i,j) = (1.0d0 - dexp(-rescale_factor_kappa * dist)) * rescale_factor_kappa_inv
        end do
     end do

  end select

end function qmckl_distance_rescaled_f

3.1.4 Performance

This function is more efficient when A and B are transposed.

4 Rescaled Distance Derivatives

4.1 `qmckl_distance_rescaled_gl`

qmckl_distance_rescaled_gl computes the matrix of the gradient and Laplacian of the rescaled distance with respect to the electron coordinates. The derivative is a rank 3 tensor. The first dimension has a dimension of 4 of which the first three coordinates contains the gradient vector and the last index is the Laplacian.

\[ C(r_{ij}) = \left( 1 - \exp(-\kappa\, r_{ij})\right)/\kappa \]

Here the gradient is defined as follows:

\[ \nabla_i C(r_{ij}) = \left(\frac{\partial C(r_{ij})}{\partial x_i},\frac{\partial C(r_{ij})}{\partial y_i},\frac{\partial C(r_{ij})}{\partial z_i} \right) \] and the Laplacian is defined as follows:

\[ \Delta_i C(r_{ij}) = \frac{\partial^2}{\partial x_i^2} + \frac{\partial^2}{\partial y_i^2} + \frac{\partial^2}{\partial z_i^2} \]

Using the above three formulas, the expression for the gradient and Laplacian is:

\[ \frac{\partial C(r_{ij})}{\partial x_i} = \frac{|(x_i - x_j)|}{r_{ij}} \exp (- \kappa \, r_{ij}) \]

\[ \Delta C_{ij}(r_{ij}) = \left[ \frac{2}{r_{ij}} - \kappa \right] \exp (- \kappa \, r_{ij}) \]

If the input array is normal ('N'), the xyz coordinates are in the leading dimension: [n][3] in C and (3,n) in Fortran.

Variable	Type	In/Out	Description
`context`	`qmckl_context`	in	Global state
`transa`	`char`	in	Array `A` is `'N'`: Normal, `'T'`: Transposed
`transb`	`char`	in	Array `B` is `'N'`: Normal, `'T'`: Transposed
`m`	`int64_t`	in	Number of points in the first set
`n`	`int64_t`	in	Number of points in the second set
`A`	`double[][lda]`	in	Array containing the \(m \times 3\) matrix \(A\)
`lda`	`int64_t`	in	Leading dimension of array `A`
`B`	`double[][ldb]`	in	Array containing the \(n \times 3\) matrix \(B\)
`ldb`	`int64_t`	in	Leading dimension of array `B`
`C`	`double[4][n][ldc]`	out	Array containing the \(4 \times m \times n\) matrix \(C\)
`ldc`	`int64_t`	in	Leading dimension of array `C`
`rescale_factor_kappa`	`double`	in	Factor for calculating rescaled distances derivatives

Requirements:

context is not QMCKL_NULL_CONTEXT
m > 0
n > 0
lda >= 3 if transa == 'N'
lda >= m if transa == 'T'
ldb >= 3 if transb == 'N'
ldb >= n if transb == 'T'
ldc >= m
A is allocated with at least \(3 \times m \times 8\) bytes
B is allocated with at least \(3 \times n \times 8\) bytes
C is allocated with at least \(4 \times m \times n \times 8\) bytes

qmckl_exit_code qmckl_distance_rescaled_gl (
      const qmckl_context context,
      const char transa,
      const char transb,
      const int64_t m,
      const int64_t n,
      const double* A,
      const int64_t lda,
      const double* B,
      const int64_t ldb,
      double* const C,
      const int64_t ldc,
      const double rescale_factor_kappa );

integer function qmckl_distance_rescaled_gl_f(context, transa, transb, m, n, &
     A, LDA, B, LDB, C, LDC, rescale_factor_kappa) &
     result(info)
  use qmckl
  implicit none
  integer(qmckl_context)  , intent(in)  :: context
  character  , intent(in)  :: transa, transb
  integer*8  , intent(in)  :: m, n
  integer*8  , intent(in)  :: lda
  real*8     , intent(in)  :: A(lda,*)
  integer*8  , intent(in)  :: ldb
  real*8     , intent(in)  :: B(ldb,*)
  integer*8  , intent(in)  :: ldc
  real*8     , intent(out) :: C(4,ldc,*)
  real*8     , intent(in)  :: rescale_factor_kappa

  integer*8 :: i,j
  real*8    :: x, y, z, dist, dist_inv
  real*8    :: rij
  integer   :: transab

  info = QMCKL_SUCCESS

  if (context == QMCKL_NULL_CONTEXT) then
     info = QMCKL_INVALID_CONTEXT
     return
  endif

  if (m <= 0_8) then
     info = QMCKL_INVALID_ARG_4
     return
  endif

  if (n <= 0_8) then
     info = QMCKL_INVALID_ARG_5
     return
  endif

  if (transa == 'N' .or. transa == 'n') then
     transab = 0
  else if (transa == 'T' .or. transa == 't') then
     transab = 1
  else
     transab = -100
  endif

  if (transb == 'N' .or. transb == 'n') then
     continue
  else if (transb == 'T' .or. transb == 't') then
     transab = transab + 2
  else
     transab = -100
  endif

  ! check for LDA
  if (transab < 0) then
     info = QMCKL_INVALID_ARG_1
     return
  endif

  if (iand(transab,1) == 0 .and. LDA < 3) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  if (iand(transab,1) == 1 .and. LDA < m) then
     info = QMCKL_INVALID_ARG_7
     return
  endif

  ! check for LDB
  if (iand(transab,2) == 0 .and. LDB < 3) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  if (iand(transab,2) == 2 .and. LDB < n) then
     info = QMCKL_INVALID_ARG_9
     return
  endif

  ! check for LDC
  if (LDC < m) then
     info = QMCKL_INVALID_ARG_11
     return
  endif


  select case (transab)

  case(0)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(1,j)
           y = A(2,i) - B(2,j)
           z = A(3,i) - B(3,j)
           dist = dsqrt(x*x + y*y + z*z)
           ! Avoid floating-point exception
           if (dist == 0.d0) then
             dist = 69.d0/rescale_factor_kappa
           endif
           dist_inv = 1.0d0/dist
           rij = dexp(-rescale_factor_kappa * dist)
           C(1,i,j) = x * dist_inv * rij
           C(2,i,j) = y * dist_inv * rij
           C(3,i,j) = z * dist_inv * rij
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

  case(1)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(1,j)
           y = A(i,2) - B(2,j)
           z = A(i,3) - B(3,j)
           dist = dsqrt(x*x + y*y + z*z)
           ! Avoid floating-point exception
           if (dist == 0.d0) then
             dist = 69.d0/rescale_factor_kappa
           endif
           dist_inv = 1.0d0/dist
           rij = dexp(-rescale_factor_kappa * dist)
           C(1,i,j) = x * dist_inv * rij
           C(2,i,j) = y * dist_inv * rij
           C(3,i,j) = z * dist_inv * rij
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

  case(2)

     do j=1,n
        do i=1,m
           x = A(1,i) - B(j,1)
           y = A(2,i) - B(j,2)
           z = A(3,i) - B(j,3)
           dist = dsqrt(x*x + y*y + z*z)
           ! Avoid floating-point exception
           if (dist == 0.d0) then
             dist = 69.d0/rescale_factor_kappa
           endif
           dist_inv = 1.0d0/dist
           rij = dexp(-rescale_factor_kappa * dist)
           C(1,i,j) = x * dist_inv * rij
           C(2,i,j) = y * dist_inv * rij
           C(3,i,j) = z * dist_inv * rij
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

  case(3)

     do j=1,n
        do i=1,m
           x = A(i,1) - B(j,1)
           y = A(i,2) - B(j,2)
           z = A(i,3) - B(j,3)
           dist = dsqrt(x*x + y*y + z*z)
           ! Avoid floating-point exception
           if (dist == 0.d0) then
             dist = 69.d0/rescale_factor_kappa
           endif
           dist_inv = 1.0d0/dist
           rij = dexp(-rescale_factor_kappa * dist)
           C(1,i,j) = x * dist_inv * rij
           C(2,i,j) = y * dist_inv * rij
           C(3,i,j) = z * dist_inv * rij
           C(4,i,j) = (2.0d0 * dist_inv - rescale_factor_kappa) * rij
        end do
     end do

  end select

end function qmckl_distance_rescaled_gl_f

This function is more efficient when A and B are transposed.

Inter-particle distances

Table of Contents

1 Squared distance

1.1 qmckl_distance_sq

1.1.1 Performance

2 Distance

2.1 qmckl_distance

2.1.1 Requirements

2.1.2 C header

2.1.3 Source

2.1.4 Performance

3 Rescaled Distance

3.1 qmckl_distance_rescaled

3.1.1 Requirements

3.1.2 C header

3.1.3 Source

3.1.4 Performance

4 Rescaled Distance Derivatives

4.1 qmckl_distance_rescaled_gl

1.1 `qmckl_distance_sq`

2.1 `qmckl_distance`

3.1 `qmckl_distance_rescaled`

4.1 `qmckl_distance_rescaled_gl`