Use mctc-lib Wigner-Seitz cell and cutoff-based EEQBC lattice sums

src/multicharge/CMakeLists.txt

@@ -28,7 +28,6 @@ list(
   "${dir}/output.f90"
   "${dir}/param.f90"
   "${dir}/version.f90"
-  "${dir}/wignerseitz.f90"
   "${dir}/cache.f90"
 )
 

src/multicharge/cache.f90

@@ -20,7 +20,7 @@
 module multicharge_model_cache
    use mctc_env, only: wp
    use mctc_io, only: structure_type
-   use multicharge_wignerseitz, only: wignerseitz_cell_type
+   use mctc_wignerseitz, only: wignerseitz_cell
    implicit none
    private
 
@@ -42,7 +42,7 @@ module multicharge_model_cache
       !> Ewald separation parameter
       real(wp) :: alpha
       !> Wigner-Seitz cell
-      type(wignerseitz_cell_type) :: wsc
+      type(wignerseitz_cell) :: wsc
    end type model_cache
 
 end module multicharge_model_cache

src/multicharge/meson.build

@@ -25,6 +25,5 @@ srcs += files(
   'output.f90',
   'param.f90',
   'version.f90',
-  'wignerseitz.f90',
   'cache.f90',
 )

src/multicharge/model/eeq.f90

@@ -28,12 +28,12 @@ module multicharge_model_eeq
    use mctc_io_constants, only: pi
    use mctc_io_math, only: matdet_3x3
    use mctc_ncoord, only: new_ncoord, cn_count
+   use mctc_wignerseitz, only: wignerseitz_cell, new_wignerseitz_cell
    use multicharge_blas, only: gemm
    use multicharge_ewald, only: get_alpha
    use multicharge_model_cache, only: cache_container, model_cache
    use multicharge_model_hessian, only: coulomb_kernel, radial_hessian, add_pair_hessian, add_pair_block
    use multicharge_model_type, only: mchrg_model_type, get_dir_trans, get_rec_trans
-   use multicharge_wignerseitz, only: wignerseitz_cell_type, new_wignerseitz_cell
    implicit none
    private
 
@@ -258,7 +258,7 @@ end subroutine add_model_hessian
 subroutine add_periodic_hessian(self, mol, wsc, alpha, qvec, adj, hessian, diagonal)
    class(eeq_model), intent(in) :: self
    type(structure_type), intent(in) :: mol
-   type(wignerseitz_cell_type), intent(in) :: wsc
+   type(wignerseitz_cell), intent(in) :: wsc
    real(wp), intent(in) :: alpha, qvec(:), adj(:)
    real(wp), intent(inout) :: hessian(:, :)
    real(wp), intent(inout) :: diagonal(:, :, :)
@@ -267,7 +267,7 @@ subroutine add_periodic_hessian(self, mol, wsc, alpha, qvec, adj, hessian, diago
    real(wp) :: val, dr, drr, ar, arr
    real(wp), allocatable :: dtrans(:, :), rtrans(:, :), rec_weight(:)
 
-   call get_dir_trans(mol%lattice, dtrans)
+   call get_dir_trans(mol, dtrans)
    call get_rec_trans(mol%lattice, rtrans)
    call get_rec_weights(abs(matdet_3x3(mol%lattice)), alpha, rtrans, rec_weight)
 
@@ -368,15 +368,15 @@ end subroutine get_amat_0d
 subroutine get_amat_3d(self, mol, wsc, alpha, amat)
    class(eeq_model), intent(in) :: self
    type(structure_type), intent(in) :: mol
-   type(wignerseitz_cell_type), intent(in) :: wsc
+   type(wignerseitz_cell), intent(in) :: wsc
    real(wp), intent(in) :: alpha
    real(wp), intent(out) :: amat(:, :)
 
    integer :: iat, jat, izp, jzp, img, ipair, npair
    real(wp) :: vec(3), gam, wsw, dtmp, rtmp, tmp
    real(wp), allocatable :: dtrans(:, :), rtrans(:, :), rec_weight(:)
 
-   call get_dir_trans(mol%lattice, dtrans)
+   call get_dir_trans(mol, dtrans)
    call get_rec_trans(mol%lattice, rtrans)
    call get_rec_weights(abs(matdet_3x3(mol%lattice)), alpha, rtrans, rec_weight)
 
@@ -512,7 +512,7 @@ end subroutine get_damat_0d
 subroutine get_damat_3d(self, mol, wsc, alpha, qvec, dadr, dadL, atrace)
    class(eeq_model), intent(in) :: self
    type(structure_type), intent(in) :: mol
-   type(wignerseitz_cell_type), intent(in) :: wsc
+   type(wignerseitz_cell), intent(in) :: wsc
    real(wp), intent(in) :: alpha, qvec(:)
    real(wp), intent(out) :: dadr(:, :, :), dadL(:, :, :), atrace(:, :)
 
@@ -523,7 +523,7 @@ subroutine get_damat_3d(self, mol, wsc, alpha, qvec, dadr, dadL, atrace)
 
    atrace(:, :) = 0.0_wp
    dadL(:, :, :) = 0.0_wp
-   call get_dir_trans(mol%lattice, dtrans)
+   call get_dir_trans(mol, dtrans)
    call get_rec_trans(mol%lattice, rtrans)
    call get_rec_weights(abs(matdet_3x3(mol%lattice)), alpha, rtrans, rec_weight, rec_strain)
 

src/multicharge/model/eeqbc.f90

@@ -27,11 +27,11 @@ module multicharge_model_eeqbc
    use mctc_io, only: structure_type
    use mctc_io_constants, only: pi
    use mctc_ncoord, only: new_ncoord, cn_count
+   use mctc_wignerseitz, only: wignerseitz_cell, new_wignerseitz_cell
    use multicharge_blas, only: gemv, gemm
    use multicharge_model_cache, only: cache_container, model_cache
    use multicharge_model_hessian, only: coulomb_kernel, add_pair_hessian, add_chain_hessian
    use multicharge_model_type, only: mchrg_model_type, get_dir_trans
-   use multicharge_wignerseitz, only: new_wignerseitz_cell, wignerseitz_cell_type
    implicit none
    private
 
@@ -286,7 +286,7 @@ subroutine get_xvec(self, mol, cache, xvec)
    call gemv(ptr%cmat, ptr%xtmp, xvec)
 
    if (any(mol%periodic)) then
-      call get_dir_trans(mol%lattice, dtrans)
+      call get_dir_trans(mol, dtrans, self%ncoord%cutoff)
       !$omp parallel default(none) &
       !$omp shared(mol, self, ptr, xvec, dtrans) private(iat, izp, img, wsw) &
       !$omp private(capi, vec, rvdw, ctmp, xvec_local)
@@ -366,7 +366,7 @@ subroutine get_xvec_derivs(self, mol, cache, dxdr, dxdL)
    call gemm(dtmpdL, ptr%cmat, dxdL)
 
    if (any(mol%periodic)) then
-      call get_dir_trans(mol%lattice, dtrans)
+      call get_dir_trans(mol, dtrans, self%ncoord%cutoff)
       !$omp parallel default(none) &
       !$omp shared(mol, self, ptr, dxdr, dxdL, dtrans) &
       !$omp private(iat, izp, jat, jzp, img, wsw) &
@@ -491,7 +491,7 @@ subroutine add_model_hessian(self, mol, error, cache, adj, hessian)
    nat = mol%nat
    periodic = any(mol%periodic)
    if (periodic) then
-      call get_dir_trans(mol%lattice, dtrans)
+      call get_dir_trans(mol, dtrans, self%ncoord%cutoff)
    else
       allocate(dtrans(3, 1), source=0.0_wp)
    end if
@@ -666,7 +666,7 @@ end subroutine get_amat_0d
 subroutine get_amat_3d(self, mol, wsc, cn, qloc, cmat, amat)
    class(eeqbc_model), intent(in) :: self
    type(structure_type), intent(in) :: mol
-   type(wignerseitz_cell_type), intent(in) :: wsc
+   type(wignerseitz_cell), intent(in) :: wsc
    real(wp), intent(in) :: cn(:), qloc(:), cmat(:, :)
    real(wp), intent(out) :: amat(:, :)
 
@@ -677,7 +677,7 @@ subroutine get_amat_3d(self, mol, wsc, cn, qloc, cmat, amat)
    ! Thread-private array for reduction
    real(wp), allocatable :: amat_local(:, :)
 
-   call get_dir_trans(mol%lattice, dtrans)
+   call get_dir_trans(mol, dtrans, self%ncoord%cutoff)
 
    amat(:, :) = 0.0_wp
 
@@ -923,7 +923,7 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, &
    & dqlocdL, cmat, dcdr, dcdL, dadr, dadL, atrace)
    class(eeqbc_model), intent(in) :: self
    type(structure_type), intent(in) :: mol
-   type(wignerseitz_cell_type), intent(in) :: wsc
+   type(wignerseitz_cell), intent(in) :: wsc
    real(wp), intent(in) :: cn(:)
    real(wp), intent(in) :: qloc(:)
    real(wp), intent(in) :: qvec(:)
@@ -948,7 +948,7 @@ subroutine get_damat_3d(self, mol, wsc, cn, qloc, qvec, dcndr, dcndL, dqlocdr, &
    real(wp), allocatable :: atrace_local(:, :)
    real(wp), allocatable :: dadr_local(:, :, :), dadL_local(:, :, :)
 
-   call get_dir_trans(mol%lattice, dtrans)
+   call get_dir_trans(mol, dtrans, self%ncoord%cutoff)
 
    allocate(dgamdr(3, mol%nat))
 
@@ -1197,7 +1197,7 @@ end subroutine get_cmat_0d
 subroutine get_cmat_3d(self, mol, wsc, cmat)
    class(eeqbc_model), intent(in) :: self
    type(structure_type), intent(in) :: mol
-   type(wignerseitz_cell_type), intent(in) :: wsc
+   type(wignerseitz_cell), intent(in) :: wsc
    real(wp), intent(out) :: cmat(:, :)
 
    integer :: iat, jat, izp, jzp, img
@@ -1207,7 +1207,7 @@ subroutine get_cmat_3d(self, mol, wsc, cmat)
    ! Thread-private array for reduction
    real(wp), allocatable :: cmat_local(:, :)
 
-   call get_dir_trans(mol%lattice, dtrans)
+   call get_dir_trans(mol, dtrans, self%ncoord%cutoff)
 
    cmat(:, :) = 0.0_wp
 
@@ -1375,7 +1375,7 @@ end subroutine get_dcmat_0d
 subroutine get_dcmat_3d(self, mol, wsc, dcdr, dcdL)
    class(eeqbc_model), intent(in) :: self
    type(structure_type), intent(in) :: mol
-   type(wignerseitz_cell_type), intent(in) :: wsc
+   type(wignerseitz_cell), intent(in) :: wsc
    real(wp), intent(out) :: dcdr(:, :, :)
    real(wp), intent(out) :: dcdL(:, :, :)
 
@@ -1386,7 +1386,7 @@ subroutine get_dcmat_3d(self, mol, wsc, dcdr, dcdL)
    ! Thread-private arrays for reduction
    real(wp), allocatable :: dcdr_local(:, :, :), dcdL_local(:, :, :)
 
-   call get_dir_trans(mol%lattice, dtrans)
+   call get_dir_trans(mol, dtrans, self%ncoord%cutoff)
 
    dcdr(:, :, :) = 0.0_wp
    dcdL(:, :, :) = 0.0_wp

src/multicharge/model/type.F90

@@ -32,7 +32,6 @@ module multicharge_model_type
    use multicharge_blas, only: gemv, symv, gemm
    use multicharge_lapack, only: sytrf, sytrs, sytri
    use multicharge_model_cache, only: model_cache, cache_container
-   use multicharge_wignerseitz, only: wignerseitz_cell_type, new_wignerseitz_cell
    implicit none
    private
 
@@ -132,12 +131,24 @@ end subroutine get_xvec_derivs
 
 contains
 
-subroutine get_dir_trans(lattice, trans)
-   real(wp), intent(in) :: lattice(:, :)
+!> Generate direct lattice translation vectors for a periodic structure
+subroutine get_dir_trans(mol, trans, cutoff)
+   !> Molecular structure data
+   type(structure_type), intent(in) :: mol
+
+   !> Translation vectors
    real(wp), allocatable, intent(out) :: trans(:, :)
+
+   !> Optional lattice-vector cutoff
+   real(wp), intent(in), optional :: cutoff
+
    integer, parameter :: rep(3) = [2, 2, 2]
 
-   call get_lattice_points(lattice, rep, .true., trans)
+   if (present(cutoff)) then
+      call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans)
+   else
+      call get_lattice_points(mol%lattice, rep, .true., trans)
+   end if
 
 end subroutine get_dir_trans
 

subprojects/.gitignore

@@ -3,3 +3,4 @@ json-fortran-8.2.5.wrap
 jonquil.wrap
 toml-f.wrap
 test-drive.wrap
+.wraplock