diff --git a/CMakeLists.txt b/CMakeLists.txt
index 756e96ec..2fd224f5 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -33,6 +33,10 @@ if(NOT TARGET "OpenMP::OpenMP_Fortran" AND WITH_OpenMP)
   find_package("OpenMP" REQUIRED)
 endif()
 
+if(WITH_ILP64)
+  message(STATUS "Using LAPACK/BLAS ILP64 interface")
+endif()
+
 if(NOT TARGET "LAPACK::LAPACK")
   find_package("custom-lapack" REQUIRED)
 endif()
@@ -42,7 +46,7 @@ if(NOT TARGET "mctc-lib::mctc-lib")
   find_package("mctc-lib" REQUIRED)
 endif()
 
-if(NOT TARGET "mstore::mstore")
+if(NOT TARGET "mstore::mstore" AND WITH_TESTS)
   find_package("mstore")
 endif()
 
@@ -83,6 +87,17 @@ if(WITH_OpenMP)
     "OpenMP::OpenMP_Fortran"
   )
 endif()
+if(WITH_ILP64)
+  target_compile_definitions(
+    "${PROJECT_NAME}-lib"
+    PUBLIC -DIK=i8
+  )
+else()
+  target_compile_definitions(
+    "${PROJECT_NAME}-lib"
+    PUBLIC -DIK=i4
+  )
+endif()
 target_include_directories(
   "${PROJECT_NAME}-lib"
   PUBLIC
@@ -126,5 +141,7 @@ install(
 )
 
 # add the testsuite
-enable_testing()
-add_subdirectory("test")
+if(WITH_TESTS)
+  enable_testing()
+  add_subdirectory("test")
+endif()
diff --git a/config/CMakeLists.txt b/config/CMakeLists.txt
index 8ddc7bcc..873a5ef4 100644
--- a/config/CMakeLists.txt
+++ b/config/CMakeLists.txt
@@ -14,9 +14,12 @@
 # limitations under the License.
 
 option(BUILD_SHARED_LIBS "Whether the libraries built should be shared" FALSE)
-
 option(WITH_OpenMP "Enable support for shared memory parallelisation with OpenMP" TRUE)
-if(NOT DEFINED "${PROJECT_NAME}-dependeny-method")
+option(WITH_ILP64 "Enable support for ILP64 BLAS/LAPACK calls" FALSE)
+option(WITH_TESTS "Enable compilation of unit tests" TRUE)
+
+
+if(NOT DEFINED "${PROJECT_NAME}-dependency-method")
   set(
     "${PROJECT_NAME}-dependency-method"
     "subproject" "cmake" "pkgconf" "fetch"
diff --git a/config/cmake/Findcustom-blas.cmake b/config/cmake/Findcustom-blas.cmake
index 3dfe5c14..7f71e627 100644
--- a/config/cmake/Findcustom-blas.cmake
+++ b/config/cmake/Findcustom-blas.cmake
@@ -13,6 +13,14 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+if ((BLA_VENDOR MATCHES ^Intel) OR (DEFINED ENV{MKLROOT}))
+  enable_language("C")
+endif()
+
+if(WITH_ILP64)
+  set(BLA_SIZEOF_INTEGER 8)
+endif()
+
 if(NOT BLAS_FOUND)
   find_package("BLAS")
   if(NOT TARGET "BLAS::BLAS")
diff --git a/config/cmake/Findcustom-lapack.cmake b/config/cmake/Findcustom-lapack.cmake
index c78569d6..8cd9a6a6 100644
--- a/config/cmake/Findcustom-lapack.cmake
+++ b/config/cmake/Findcustom-lapack.cmake
@@ -13,6 +13,14 @@
 # See the License for the specific language governing permissions and
 # limitations under the License.
 
+if ((BLA_VENDOR MATCHES ^Intel) OR (DEFINED ENV{MKLROOT}))
+  enable_language("C")
+endif()
+
+if(WITH_ILP64)
+  set(BLA_SIZEOF_INTEGER 8)
+endif()
+
 if(NOT LAPACK_FOUND)
   find_package("LAPACK")
 
diff --git a/config/meson.build b/config/meson.build
index 09a5abb7..e35c7acd 100644
--- a/config/meson.build
+++ b/config/meson.build
@@ -49,6 +49,14 @@ if lapack_vendor == 'auto'
   endif
 endif
 
+if get_option('ilp64')
+  ilp64 = true
+  add_project_arguments('-DIK=i8', language:'fortran')
+else
+  add_project_arguments('-DIK=i4', language:'fortran')
+  ilp64 = false
+endif
+
 if lapack_vendor == 'mkl'
   mkl_dep = []
   cc = fc
@@ -56,12 +64,12 @@ if lapack_vendor == 'mkl'
     cc = meson.get_compiler('c')
   endif
   if fc_id == 'intel'
-    mkl_dep += cc.find_library('mkl_intel_lp64')
+    mkl_dep += cc.find_library(ilp64 ? 'mkl_intel_ilp64' : 'mkl_intel_lp64')
     if get_option('openmp')
       mkl_dep += cc.find_library('mkl_intel_thread')
     endif
   elif fc_id == 'gcc'
-    mkl_dep += cc.find_library('mkl_gf_lp64')
+    mkl_dep += cc.find_library(ilp64 ? 'mkl_gf_ilp64' : 'mkl_gf_lp64')
     if get_option('openmp')
       mkl_dep += cc.find_library('mkl_gnu_thread')
     endif
@@ -79,33 +87,50 @@ elif lapack_vendor == 'mkl-rt'
   lib_deps += mkl_dep
 
 elif lapack_vendor == 'openblas'
-  openblas_dep = dependency('openblas', required: false)
+  openblas_dep = dependency(ilp64 ? 'openblas64' : 'openblas', required: false)
   if not openblas_dep.found()
-    openblas_dep = fc.find_library('openblas')
+    openblas_dep = fc.find_library(ilp64 ? 'openblas64' : 'openblas') 
   endif
   lib_deps += openblas_dep
   if not fc.links('external dsytrs; call dsytrs(); end', dependencies: openblas_dep)
-    lapack_dep = dependency('lapack', required: false)
+    lapack_dep = dependency(ilp64 ? 'lapack64' : 'lapack', required: false)
     if not lapack_dep.found()
-      lapack_dep = fc.find_library('lapack')
+      lapack_dep = fc.find_library(ilp64 ? 'lapack64' : 'lapack')
     endif
     lib_deps += lapack_dep
   endif
 
 elif lapack_vendor == 'custom'
-  foreach lib: get_option('custom_libraries')
-    lib_deps += fc.find_library(lib)
-  endforeach
+  custom_deps = []
+  libs = get_option('custom_libraries')
+  if libs[0].startswith('-L')
+    foreach lib: libs
+      if lib != libs[0]
+        custom_deps += fc.find_library(lib, dirs: libs[0].substring(2))
+      endif
+    endforeach
+  else
+    foreach lib: libs
+      custom_deps += fc.find_library(lib)
+    endforeach
+  endif
+  if (not fc.links('external dsytrs; call dsytrs(); end', dependencies: [custom_deps,omp_dep]))
+    error('Custom LAPACK libraries do not link')
+  elif (not fc.links('external dsytrs; call dgemm(); end', dependencies: [custom_deps,omp_dep]))
+    error('Custom BLAS libraries do not link')
+  endif
+  lib_deps += custom_deps
+
 
 else
-  lapack_dep = dependency('lapack', required: false)
+  lapack_dep = dependency(ilp64 ? 'lapack64' : 'lapack', required: false)
   if not lapack_dep.found()
-    lapack_dep = fc.find_library('lapack')
+    lapack_dep = fc.find_library(ilp64 ? 'lapack64' : 'lapack') 
   endif
   lib_deps += lapack_dep
-  blas_dep = dependency('blas', required: false)
+  blas_dep = dependency(ilp64 ? 'blas64' : 'blas', required: false) 
   if not blas_dep.found()
-    blas_dep = fc.find_library('blas')
+    blas_dep = fc.find_library(ilp64 ? 'blas64' : 'blas') 
   endif
   lib_deps += blas_dep
 endif
diff --git a/meson_options.txt b/meson_options.txt
index 4d279d2d..2b3d156e 100644
--- a/meson_options.txt
+++ b/meson_options.txt
@@ -37,3 +37,11 @@ option(
   yield: true,
   description: 'use OpenMP parallelisation',
 )
+
+option(
+  'ilp64',
+  type: 'boolean',
+  value: false,
+  yield: true,
+  description: 'enable ILP64 LAPACK/BLAS',
+)
diff --git a/src/multicharge/CMakeLists.txt b/src/multicharge/CMakeLists.txt
index 5b87d544..0677d47f 100644
--- a/src/multicharge/CMakeLists.txt
+++ b/src/multicharge/CMakeLists.txt
@@ -20,12 +20,12 @@ set(dir "${CMAKE_CURRENT_SOURCE_DIR}")
 
 list(
   APPEND srcs
-  "${dir}/blas.f90"
+  "${dir}/blas.F90"
   "${dir}/cutoff.f90"
   "${dir}/data.f90"
   "${dir}/ewald.f90"
-  "${dir}/lapack.f90"
-  "${dir}/model.f90"
+  "${dir}/lapack.F90"
+  "${dir}/model.F90"
   "${dir}/ncoord.f90"
   "${dir}/output.f90"
   "${dir}/param.f90"
diff --git a/src/multicharge/blas.f90 b/src/multicharge/blas.F90
similarity index 91%
rename from src/multicharge/blas.f90
rename to src/multicharge/blas.F90
index f7475828..3e345a52 100644
--- a/src/multicharge/blas.f90
+++ b/src/multicharge/blas.F90
@@ -13,9 +13,13 @@
 ! See the License for the specific language governing permissions and
 ! limitations under the License.
 
+#ifndef IK
+#define IK i4
+#endif
+
 !> Interface to BLAS library for matrix-vector and matrix-matrix operations
 module multicharge_blas
-   use mctc_env, only : sp, dp
+   use mctc_env, only : sp, dp, ik => IK
    implicit none
    private
 
@@ -78,32 +82,32 @@ module multicharge_blas
    !> m by n matrix.
    interface blas_gemv
       pure subroutine sgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
-         import :: sp
+         import :: sp, ik
+         integer(ik), intent(in) :: lda
          real(sp), intent(in) :: a(lda, *)
          real(sp), intent(in) :: x(*)
          real(sp), intent(inout) :: y(*)
          real(sp), intent(in) :: alpha
          real(sp), intent(in) :: beta
          character(len=1), intent(in) :: trans
-         integer, intent(in) :: incx
-         integer, intent(in) :: incy
-         integer, intent(in) :: m
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(in) :: incx
+         integer(ik), intent(in) :: incy
+         integer(ik), intent(in) :: m
+         integer(ik), intent(in) :: n
       end subroutine sgemv
       pure subroutine dgemv(trans, m, n, alpha, a, lda, x, incx, beta, y, incy)
-         import :: dp
+         import :: dp, ik
+         integer(ik), intent(in) :: lda
          real(dp), intent(in) :: a(lda, *)
          real(dp), intent(in) :: x(*)
          real(dp), intent(inout) :: y(*)
          real(dp), intent(in) :: alpha
          real(dp), intent(in) :: beta
          character(len=1), intent(in) :: trans
-         integer, intent(in) :: incx
-         integer, intent(in) :: incy
-         integer, intent(in) :: m
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(in) :: incx
+         integer(ik), intent(in) :: incy
+         integer(ik), intent(in) :: m
+         integer(ik), intent(in) :: n
       end subroutine dgemv
    end interface blas_gemv
 
@@ -115,30 +119,30 @@ end subroutine dgemv
    !> A is an n by n symmetric matrix.
    interface blas_symv
       pure subroutine ssymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
-         import :: sp
+         import :: sp, ik
+         integer(ik), intent(in) :: lda
          real(sp), intent(in) :: a(lda, *)
          real(sp), intent(in) :: x(*)
          real(sp), intent(inout) :: y(*)
          character(len=1), intent(in) :: uplo
          real(sp), intent(in) :: alpha
          real(sp), intent(in) :: beta
-         integer, intent(in) :: incx
-         integer, intent(in) :: incy
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(in) :: incx
+         integer(ik), intent(in) :: incy
+         integer(ik), intent(in) :: n
       end subroutine ssymv
       pure subroutine dsymv(uplo, n, alpha, a, lda, x, incx, beta, y, incy)
-         import :: dp
+         import :: dp, ik
+         integer(ik), intent(in) :: lda
          real(dp), intent(in) :: a(lda, *)
          real(dp), intent(in) :: x(*)
          real(dp), intent(inout) :: y(*)
          character(len=1), intent(in) :: uplo
          real(dp), intent(in) :: alpha
          real(dp), intent(in) :: beta
-         integer, intent(in) :: incx
-         integer, intent(in) :: incy
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(in) :: incx
+         integer(ik), intent(in) :: incy
+         integer(ik), intent(in) :: n
       end subroutine dsymv
    end interface blas_symv
 
@@ -155,7 +159,10 @@ end subroutine dsymv
    interface blas_gemm
       pure subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, &
             & beta, c, ldc)
-         import :: sp
+         import :: sp, ik
+         integer(ik), intent(in) :: lda
+         integer(ik), intent(in) :: ldb
+         integer(ik), intent(in) :: ldc
          real(sp), intent(in) :: a(lda, *)
          real(sp), intent(in) :: b(ldb, *)
          real(sp), intent(inout) :: c(ldc, *)
@@ -163,16 +170,16 @@ pure subroutine sgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, &
          character(len=1), intent(in) :: transb
          real(sp), intent(in) :: alpha
          real(sp), intent(in) :: beta
-         integer, intent(in) :: m
-         integer, intent(in) :: n
-         integer, intent(in) :: k
-         integer, intent(in) :: lda
-         integer, intent(in) :: ldb
-         integer, intent(in) :: ldc
+         integer(ik), intent(in) :: m
+         integer(ik), intent(in) :: n
+         integer(ik), intent(in) :: k
       end subroutine sgemm
       pure subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, &
             & beta, c, ldc)
-         import :: dp
+         import :: dp, ik
+         integer(ik), intent(in) :: lda
+         integer(ik), intent(in) :: ldb
+         integer(ik), intent(in) :: ldc
          real(dp), intent(in) :: a(lda, *)
          real(dp), intent(in) :: b(ldb, *)
          real(dp), intent(inout) :: c(ldc, *)
@@ -180,12 +187,9 @@ pure subroutine dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, &
          character(len=1), intent(in) :: transb
          real(dp), intent(in) :: alpha
          real(dp), intent(in) :: beta
-         integer, intent(in) :: m
-         integer, intent(in) :: n
-         integer, intent(in) :: k
-         integer, intent(in) :: lda
-         integer, intent(in) :: ldb
-         integer, intent(in) :: ldc
+         integer(ik), intent(in) :: m
+         integer(ik), intent(in) :: n
+         integer(ik), intent(in) :: k
       end subroutine dgemm
    end interface blas_gemm
 
@@ -302,7 +306,7 @@ pure subroutine mchrg_sgemv(amat, xvec, yvec, alpha, beta, trans)
    character(len=1), intent(in), optional :: trans
    real(sp) :: a, b
    character(len=1) :: tra
-   integer :: incx, incy, m, n, lda
+   integer(ik) :: incx, incy, m, n, lda
    if (present(alpha)) then
       a = alpha
    else
@@ -318,8 +322,8 @@ pure subroutine mchrg_sgemv(amat, xvec, yvec, alpha, beta, trans)
    else
       tra = 'n'
    end if
-   incx = 1
-   incy = 1
+   incx = 1_ik
+   incy = 1_ik
    lda = max(1, size(amat, 1))
    m = size(amat, 1)
    n = size(amat, 2)
@@ -336,7 +340,7 @@ pure subroutine mchrg_dgemv(amat, xvec, yvec, alpha, beta, trans)
    character(len=1), intent(in), optional :: trans
    real(dp) :: a, b
    character(len=1) :: tra
-   integer :: incx, incy, m, n, lda
+   integer(ik) :: incx, incy, m, n, lda
    if (present(alpha)) then
       a = alpha
    else
@@ -352,8 +356,8 @@ pure subroutine mchrg_dgemv(amat, xvec, yvec, alpha, beta, trans)
    else
       tra = 'n'
    end if
-   incx = 1
-   incy = 1
+   incx = 1_ik
+   incy = 1_ik
    lda = max(1, size(amat, 1))
    m = size(amat, 1)
    n = size(amat, 2)
@@ -370,7 +374,7 @@ pure subroutine mchrg_ssymv(amat, xvec, yvec, uplo, alpha, beta)
    real(sp), intent(in), optional :: beta
    character(len=1) :: ula
    real(sp) :: a, b
-   integer :: incx, incy, n, lda
+   integer(ik) :: incx, incy, n, lda
    if (present(alpha)) then
       a = alpha
    else
@@ -386,8 +390,8 @@ pure subroutine mchrg_ssymv(amat, xvec, yvec, uplo, alpha, beta)
    else
       ula = 'u'
    end if
-   incx = 1
-   incy = 1
+   incx = 1_ik
+   incy = 1_ik
    lda = max(1, size(amat, 1))
    n = size(amat, 2)
    call blas_symv(ula, n, a, amat, lda, xvec, incx, b, yvec, incy)
@@ -403,7 +407,7 @@ pure subroutine mchrg_dsymv(amat, xvec, yvec, uplo, alpha, beta)
    real(dp), intent(in), optional :: beta
    character(len=1) :: ula
    real(dp) :: a, b
-   integer :: incx, incy, n, lda
+   integer(ik) :: incx, incy, n, lda
    if (present(alpha)) then
       a = alpha
    else
@@ -419,8 +423,8 @@ pure subroutine mchrg_dsymv(amat, xvec, yvec, uplo, alpha, beta)
    else
       ula = 'u'
    end if
-   incx = 1
-   incy = 1
+   incx = 1_ik
+   incy = 1_ik
    lda = max(1, size(amat, 1))
    n = size(amat, 2)
    call blas_symv(ula, n, a, amat, lda, xvec, incx, b, yvec, incy)
@@ -437,7 +441,7 @@ pure subroutine mchrg_sgemm(amat, bmat, cmat, transa, transb, alpha, beta)
    real(sp), intent(in), optional :: beta
    character(len=1) :: tra, trb
    real(sp) :: a, b
-   integer :: m, n, k, lda, ldb, ldc
+   integer(ik) :: m, n, k, lda, ldb, ldc
    if (present(alpha)) then
       a = alpha
    else
@@ -482,7 +486,7 @@ pure subroutine mchrg_dgemm(amat, bmat, cmat, transa, transb, alpha, beta)
    real(dp), intent(in), optional :: beta
    character(len=1) :: tra, trb
    real(dp) :: a, b
-   integer :: m, n, k, lda, ldb, ldc
+   integer(ik) :: m, n, k, lda, ldb, ldc
    if (present(alpha)) then
       a = alpha
    else
diff --git a/src/multicharge/lapack.f90 b/src/multicharge/lapack.F90
similarity index 75%
rename from src/multicharge/lapack.f90
rename to src/multicharge/lapack.F90
index 7d44516d..47e3abbf 100644
--- a/src/multicharge/lapack.f90
+++ b/src/multicharge/lapack.F90
@@ -13,8 +13,12 @@
 ! See the License for the specific language governing permissions and
 ! limitations under the License.
 
+#ifndef IK
+#define IK i4
+#endif
+
 module multicharge_lapack
-   use mctc_env, only : sp, dp
+   use mctc_env, only : sp, dp, ik => IK
    implicit none
    private
 
@@ -42,75 +46,75 @@ module multicharge_lapack
 
    interface lapack_sytrf
       pure subroutine ssytrf(uplo, n, a, lda, ipiv, work, lwork, info)
-         import :: sp
+         import :: sp, ik
+         integer(ik), intent(in) :: lda
          real(sp), intent(inout) :: a(lda, *)
          character(len=1), intent(in) :: uplo
-         integer, intent(out) :: ipiv(*)
-         integer, intent(out) :: info
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(out) :: ipiv(*)
+         integer(ik), intent(out) :: info
+         integer(ik), intent(in) :: n
          real(sp), intent(inout) :: work(*)
-         integer, intent(in) :: lwork
+         integer(ik), intent(in) :: lwork
       end subroutine ssytrf
       pure subroutine dsytrf(uplo, n, a, lda, ipiv, work, lwork, info)
-         import :: dp
+         import :: dp, ik
+         integer(ik), intent(in) :: lda
          real(dp), intent(inout) :: a(lda, *)
          character(len=1), intent(in) :: uplo
-         integer, intent(out) :: ipiv(*)
-         integer, intent(out) :: info
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(out) :: ipiv(*)
+         integer(ik), intent(out) :: info
+         integer(ik), intent(in) :: n
          real(dp), intent(inout) :: work(*)
-         integer, intent(in) :: lwork
+         integer(ik), intent(in) :: lwork
       end subroutine dsytrf
    end interface lapack_sytrf
 
    interface lapack_sytrs
       pure subroutine ssytrs(uplo, n, nrhs, a, lda, ipiv, b, ldb, info)
-         import :: sp
+         import :: sp, ik
+         integer(ik), intent(in) :: lda
+         integer(ik), intent(in) :: ldb
          real(sp), intent(in) :: a(lda, *)
          real(sp), intent(inout) :: b(ldb, *)
-         integer, intent(in) :: ipiv(*)
+         integer(ik), intent(in) :: ipiv(*)
          character(len=1), intent(in) :: uplo
-         integer, intent(out) :: info
-         integer, intent(in) :: n
-         integer, intent(in) :: nrhs
-         integer, intent(in) :: lda
-         integer, intent(in) :: ldb
+         integer(ik), intent(out) :: info
+         integer(ik), intent(in) :: n
+         integer(ik), intent(in) :: nrhs
       end subroutine ssytrs
       pure subroutine dsytrs(uplo, n, nrhs, a, lda, ipiv, b, ldb, info)
-         import :: dp
+         import :: dp, ik
+         integer(ik), intent(in) :: lda
+         integer(ik), intent(in) :: ldb
          real(dp), intent(in) :: a(lda, *)
          real(dp), intent(inout) :: b(ldb, *)
-         integer, intent(in) :: ipiv(*)
+         integer(ik), intent(in) :: ipiv(*)
          character(len=1), intent(in) :: uplo
-         integer, intent(out) :: info
-         integer, intent(in) :: n
-         integer, intent(in) :: nrhs
-         integer, intent(in) :: lda
-         integer, intent(in) :: ldb
+         integer(ik), intent(out) :: info
+         integer(ik), intent(in) :: n
+         integer(ik), intent(in) :: nrhs
       end subroutine dsytrs
    end interface lapack_sytrs
 
    interface lapack_sytri
       pure subroutine ssytri(uplo, n, a, lda, ipiv, work, info)
-         import :: sp
+         import :: sp, ik
+         integer(ik), intent(in) :: lda
          real(sp), intent(inout) :: a(lda, *)
-         integer, intent(in) :: ipiv(*)
+         integer(ik), intent(in) :: ipiv(*)
          character(len=1), intent(in) :: uplo
-         integer, intent(out) :: info
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(out) :: info
+         integer(ik), intent(in) :: n
          real(sp), intent(in) :: work(*)
       end subroutine ssytri
       pure subroutine dsytri(uplo, n, a, lda, ipiv, work, info)
-         import :: dp
+         import :: dp, ik
+         integer(ik), intent(in) :: lda
          real(dp), intent(inout) :: a(lda, *)
-         integer, intent(in) :: ipiv(*)
+         integer(ik), intent(in) :: ipiv(*)
          character(len=1), intent(in) :: uplo
-         integer, intent(out) :: info
-         integer, intent(in) :: n
-         integer, intent(in) :: lda
+         integer(ik), intent(out) :: info
+         integer(ik), intent(in) :: n
          real(dp), intent(in) :: work(*)
       end subroutine dsytri
    end interface lapack_sytri
@@ -121,11 +125,11 @@ end subroutine dsytri
 
 subroutine mchrg_ssytrf(amat, ipiv, uplo, info)
    real(sp), intent(inout) :: amat(:, :)
-   integer, intent(out) :: ipiv(:)
+   integer(ik), intent(out) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    character(len=1) :: ula
-   integer :: stat, n, lda, lwork, stat_alloc, stat_dealloc
+   integer(ik) :: stat, n, lda, lwork, stat_alloc, stat_dealloc
    real(sp), allocatable :: work(:)
    real(sp) :: test(1)
    if (present(uplo)) then
@@ -135,8 +139,8 @@ subroutine mchrg_ssytrf(amat, ipiv, uplo, info)
    end if
    lda = max(1, size(amat, 1))
    n = size(amat, 2)
-   stat_alloc = 0
-   lwork = -1
+   stat_alloc = 0_ik
+   lwork = -1_ik
    call lapack_sytrf(ula, n, amat, lda, ipiv, test, lwork, stat)
    if (stat == 0) then
       lwork = nint(test(1))
@@ -146,7 +150,7 @@ subroutine mchrg_ssytrf(amat, ipiv, uplo, info)
       if (stat_alloc==0) then
          call lapack_sytrf(ula, n, amat, lda, ipiv, work, lwork, stat)
       else
-         stat = -1000
+         stat = -1000_ik
       end if
       deallocate(work, stat=stat_dealloc)
    end if
@@ -160,11 +164,11 @@ end subroutine mchrg_ssytrf
 
 subroutine mchrg_dsytrf(amat, ipiv, uplo, info)
    real(dp), intent(inout) :: amat(:, :)
-   integer, intent(out) :: ipiv(:)
+   integer(ik), intent(out) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    character(len=1) :: ula
-   integer :: stat, n, lda, lwork, stat_alloc, stat_dealloc
+   integer(ik) :: stat, n, lda, lwork, stat_alloc, stat_dealloc
    real(dp), allocatable :: work(:)
    real(dp) :: test(1)
    if (present(uplo)) then
@@ -174,8 +178,8 @@ subroutine mchrg_dsytrf(amat, ipiv, uplo, info)
    end if
    lda = max(1, size(amat, 1))
    n = size(amat, 2)
-   stat_alloc = 0
-   lwork = -1
+   stat_alloc = 0_ik
+   lwork = -1_ik
    call lapack_sytrf(ula, n, amat, lda, ipiv, test, lwork, stat)
    if (stat == 0) then
       lwork = nint(test(1))
@@ -185,7 +189,7 @@ subroutine mchrg_dsytrf(amat, ipiv, uplo, info)
       if (stat_alloc==0) then
          call lapack_sytrf(ula, n, amat, lda, ipiv, work, lwork, stat)
       else
-         stat = -1000
+         stat = -1000_ik
       end if
       deallocate(work, stat=stat_dealloc)
    end if
@@ -200,11 +204,11 @@ end subroutine mchrg_dsytrf
 subroutine mchrg_ssytrs(amat, bmat, ipiv, uplo, info)
    real(sp), intent(in) :: amat(:, :)
    real(sp), intent(inout) :: bmat(:, :)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    character(len=1) :: ula
-   integer :: stat, n, nrhs, lda, ldb
+   integer(ik) :: stat, n, nrhs, lda, ldb
    if (present(uplo)) then
       ula = uplo
    else
@@ -226,11 +230,11 @@ end subroutine mchrg_ssytrs
 subroutine mchrg_dsytrs(amat, bmat, ipiv, uplo, info)
    real(dp), intent(in) :: amat(:, :)
    real(dp), intent(inout) :: bmat(:, :)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    character(len=1) :: ula
-   integer :: stat, n, nrhs, lda, ldb
+   integer(ik) :: stat, n, nrhs, lda, ldb
    if (present(uplo)) then
       ula = uplo
    else
@@ -252,9 +256,9 @@ end subroutine mchrg_dsytrs
 subroutine mchrg_ssytrs1(amat, bvec, ipiv, uplo, info)
    real(sp), intent(in) :: amat(:, :)
    real(sp), intent(inout), target :: bvec(:)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    real(sp), pointer :: bptr(:, :)
    bptr(1:size(bvec), 1:1) => bvec
    call sytrs(amat, bptr, ipiv, uplo, info)
@@ -264,9 +268,9 @@ end subroutine mchrg_ssytrs1
 subroutine mchrg_ssytrs3(amat, bmat, ipiv, uplo, info)
    real(sp), intent(in) :: amat(:, :)
    real(sp), intent(inout), contiguous, target :: bmat(:, :, :)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    real(sp), pointer :: bptr(:, :)
    bptr(1:size(bmat, 1), 1:size(bmat, 2)*size(bmat, 3)) => bmat
    call sytrs(amat, bptr, ipiv, uplo, info)
@@ -276,9 +280,9 @@ end subroutine mchrg_ssytrs3
 subroutine mchrg_dsytrs1(amat, bvec, ipiv, uplo, info)
    real(dp), intent(in) :: amat(:, :)
    real(dp), intent(inout), target :: bvec(:)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    real(dp), pointer :: bptr(:, :)
    bptr(1:size(bvec), 1:1) => bvec
    call sytrs(amat, bptr, ipiv, uplo, info)
@@ -288,9 +292,9 @@ end subroutine mchrg_dsytrs1
 subroutine mchrg_dsytrs3(amat, bmat, ipiv, uplo, info)
    real(dp), intent(in) :: amat(:, :)
    real(dp), intent(inout), contiguous, target :: bmat(:, :, :)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    real(dp), pointer :: bptr(:, :)
    bptr(1:size(bmat, 1), 1:size(bmat, 2)*size(bmat, 3)) => bmat
    call sytrs(amat, bptr, ipiv, uplo, info)
@@ -299,11 +303,11 @@ end subroutine mchrg_dsytrs3
 
 subroutine mchrg_ssytri(amat, ipiv, uplo, info)
    real(sp), intent(inout) :: amat(:, :)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    character(len=1) :: ula
-   integer :: stat, n, lda, stat_alloc, stat_dealloc
+   integer(ik) :: stat, n, lda, stat_alloc, stat_dealloc
    real(sp), allocatable :: work(:)
    if (present(uplo)) then
       ula = uplo
@@ -312,12 +316,12 @@ subroutine mchrg_ssytri(amat, ipiv, uplo, info)
    end if
    lda = max(1, size(amat, 1))
    n = size(amat, 2)
-   stat_alloc = 0
+   stat_alloc = 0_ik
    allocate(work(n), stat=stat_alloc)
    if (stat_alloc==0) then
       call lapack_sytri(ula, n, amat, lda, ipiv, work, stat)
    else
-      stat = -1000
+      stat = -1000_ik
    end if
    deallocate(work, stat=stat_dealloc)
    if (present(info)) then
@@ -330,11 +334,11 @@ end subroutine mchrg_ssytri
 
 subroutine mchrg_dsytri(amat, ipiv, uplo, info)
    real(dp), intent(inout) :: amat(:, :)
-   integer, intent(in) :: ipiv(:)
+   integer(ik), intent(in) :: ipiv(:)
    character(len=1), intent(in), optional :: uplo
-   integer, intent(out), optional :: info
+   integer(ik), intent(out), optional :: info
    character(len=1) :: ula
-   integer :: stat, n, lda, stat_alloc, stat_dealloc
+   integer(ik) :: stat, n, lda, stat_alloc, stat_dealloc
    real(dp), allocatable :: work(:)
    if (present(uplo)) then
       ula = uplo
@@ -343,12 +347,12 @@ subroutine mchrg_dsytri(amat, ipiv, uplo, info)
    end if
    lda = max(1, size(amat, 1))
    n = size(amat, 2)
-   stat_alloc = 0
+   stat_alloc = 0_ik
    allocate(work(n), stat=stat_alloc)
    if (stat_alloc==0) then
       call lapack_sytri(ula, n, amat, lda, ipiv, work, stat)
    else
-      stat = -1000
+      stat = -1000_ik
    end if
    deallocate(work, stat=stat_dealloc)
    if (present(info)) then
diff --git a/src/multicharge/meson.build b/src/multicharge/meson.build
index 191fa022..2fcb6a5c 100644
--- a/src/multicharge/meson.build
+++ b/src/multicharge/meson.build
@@ -17,12 +17,12 @@ subdir('data')
 subdir('param')
 
 srcs += files(
-  'blas.f90',
+  'blas.F90',
   'cutoff.f90',
   'data.f90',
   'ewald.f90',
-  'lapack.f90',
-  'model.f90',
+  'lapack.F90',
+  'model.F90',
   'ncoord.f90',
   'output.f90',
   'param.f90',
diff --git a/src/multicharge/model.f90 b/src/multicharge/model.F90
similarity index 98%
rename from src/multicharge/model.f90
rename to src/multicharge/model.F90
index a6f385a9..6d2b28e3 100644
--- a/src/multicharge/model.f90
+++ b/src/multicharge/model.F90
@@ -13,8 +13,12 @@
 ! See the License for the specific language governing permissions and
 ! limitations under the License.
 
+#ifndef IK
+#define IK i4
+#endif
+
 module multicharge_model
-   use mctc_env, only : error_type, wp
+   use mctc_env, only : error_type, wp, ik => IK
    use mctc_io, only : structure_type
    use mctc_io_constants, only : pi
    use mctc_io_math, only : matdet_3x3, matinv_3x3
@@ -431,10 +435,12 @@ subroutine solve(self, mol, cn, dcndr, dcndL, energy, gradient, sigma, qvec, dqd
    real(wp), intent(inout), contiguous, optional :: gradient(:, :)
    real(wp), intent(inout), contiguous, optional :: sigma(:, :)
 
-   integer :: ic, jc, iat, ndim, info
+   integer :: ic, jc, iat, ndim
    logical :: grad, cpq, dcn
    real(wp) :: alpha
-   integer, allocatable :: ipiv(:)
+   integer(ik) :: info
+   integer(ik), allocatable :: ipiv(:)
+
    real(wp), allocatable :: xvec(:), vrhs(:), amat(:, :), ainv(:, :)
    real(wp), allocatable :: dxdcn(:), atrace(:, :), dadr(:, :, :), dadL(:, :, :)
    type(wignerseitz_cell_type) :: wsc