Periodic EEQ-BC implementation complete

app/main.f90

@@ -17,7 +17,7 @@ program main
    use, intrinsic :: iso_fortran_env, only: output_unit, error_unit, input_unit
    use mctc_env, only: error_type, fatal_error, get_argument, wp
    use mctc_io, only: structure_type, read_structure, filetype, get_filetype
-   use mctc_cutoff, only : get_lattice_points
+   use mctc_cutoff, only: get_lattice_points
    use multicharge, only: mchrg_model_type, mchargeModel, new_eeq2019_model, &
       & new_eeqbc2025_model, get_multicharge_version, &
       & write_ascii_model, write_ascii_properties, write_ascii_results

src/multicharge/cache.f90

@@ -20,8 +20,7 @@
 module multicharge_model_cache
    use mctc_env, only: wp
    use mctc_io, only: structure_type
-   use multicharge_wignerseitz, only: wignerseitz_cell_type, new_wignerseitz_cell
-   use multicharge_ewald, only: get_alpha
+   use multicharge_wignerseitz, only: wignerseitz_cell_type
    implicit none
    private
 
@@ -40,23 +39,8 @@ module multicharge_model_cache
       real(wp), allocatable :: dcndL(:, :, :)
       !> Ewald separation parameter
       real(wp) :: alpha
+      !> Wigner-Seitz cell
       type(wignerseitz_cell_type) :: wsc
-   contains
-      !> Create WSC
-      procedure :: update
    end type model_cache
 
-contains
-   subroutine update(self, mol)
-      class(model_cache), intent(inout) :: self
-      type(structure_type), intent(in) :: mol
-
-      ! Create WSC
-      if (any(mol%periodic)) then
-         call new_wignerseitz_cell(self%wsc, mol)
-         call get_alpha(mol%lattice, self%alpha)
-      end if
-
-   end subroutine update
-
 end module multicharge_model_cache

src/multicharge/model/eeq.f90

@@ -17,7 +17,7 @@
 !> Provides implementation of the electronegativity equilibration model (EEQ)
 
 !> Electronegativity equlibration charge model published in
-!> 
+!>
 !> E. Caldeweyher, S. Ehlert, A. Hansen, H. Neugebauer, S. Spicher, C. Bannwarth
 !> and S. Grimme, *J. Chem. Phys.*, **2019**, 150, 154122.
 !> DOI: [10.1063/1.5090222](https://dx.doi.org/10.1063/1.5090222)
@@ -27,7 +27,8 @@ 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 multicharge_wignerseitz, only: wignerseitz_cell_type
+   use multicharge_wignerseitz, only: wignerseitz_cell_type, new_wignerseitz_cell
+   use multicharge_ewald, only: get_alpha
    use multicharge_model_type, only: mchrg_model_type, get_dir_trans, get_rec_trans
    use multicharge_model_cache, only: cache_container, model_cache
    implicit none
@@ -50,14 +51,6 @@ module multicharge_model_eeq
       procedure :: get_xvec
       !> Calculate EN vector derivatives
       procedure :: get_xvec_derivs
-      !> Calculate Coulomb matrix
-      procedure :: get_amat_0d
-      !> Calculate Coulomb matrix periodic
-      procedure :: get_amat_3d
-      !> Calculate Coulomb matrix derivative
-      procedure :: get_damat_0d
-      !> Calculate Coulomb matrix derivative periodic
-      procedure :: get_damat_3d
    end type eeq_model
 
    real(wp), parameter :: sqrtpi = sqrt(pi)
@@ -115,14 +108,20 @@ subroutine update(self, mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL)
       type(eeq_cache), pointer :: ptr
 
       call taint(cache, ptr)
-      call ptr%update(mol)
 
       ! Refer CN arrays in cache
       ptr%cn = cn
       if (present(dcndr) .and. present(dcndL)) then
          ptr%dcndr = dcndr
          ptr%dcndL = dcndL
       end if
+
+      if (any(mol%periodic)) then
+         ! Create WSC
+         call new_wignerseitz_cell(ptr%wsc, mol)
+         call get_alpha(mol%lattice, ptr%alpha)
+      end if
+
    end subroutine update
 
    subroutine get_xvec(self, mol, cache, xvec)
@@ -154,8 +153,8 @@ subroutine get_xvec_derivs(self, mol, cache, dxdr, dxdL)
       class(eeq_model), intent(in) :: self
       type(structure_type), intent(in) :: mol
       type(cache_container), intent(inout) :: cache
-      real(wp), intent(out) :: dxdr(:, :, :)
-      real(wp), intent(out) :: dxdL(:, :, :)
+      real(wp), intent(out), contiguous :: dxdr(:, :, :)
+      real(wp), intent(out), contiguous :: dxdL(:, :, :)
       real(wp), parameter :: reg = 1.0e-14_wp
 
       integer :: iat, izp
@@ -191,9 +190,9 @@ subroutine get_coulomb_matrix(self, mol, cache, amat)
       call view(cache, ptr)
 
       if (any(mol%periodic)) then
-         call self%get_amat_3d(mol, ptr%wsc, ptr%alpha, amat)
+         call get_amat_3d(self, mol, ptr%wsc, ptr%alpha, amat)
       else
-         call self%get_amat_0d(mol, amat)
+         call get_amat_0d(self, mol, amat)
       end if
    end subroutine get_coulomb_matrix
 
@@ -207,33 +206,33 @@ subroutine get_amat_0d(self, mol, amat)
 
       ! Thread-private array for reduction
       real(wp), allocatable :: amat_local(:, :)
-   
+
       amat(:, :) = 0.0_wp
-   
+
       !$omp parallel default(none) &
       !$omp shared(amat, mol, self) &
       !$omp private(iat, izp, jat, jzp, gam, vec, r2, tmp, amat_local)
-      allocate(amat_local, source=amat)
-      !$omp do schedule(runtime) 
+      allocate (amat_local, source=amat)
+      !$omp do schedule(runtime)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          do jat = 1, iat - 1
             jzp = mol%id(jat)
             vec = mol%xyz(:, jat) - mol%xyz(:, iat)
             r2 = vec(1)**2 + vec(2)**2 + vec(3)**2
-            gam = 1.0_wp / (self%rad(izp)**2 + self%rad(jzp)**2)
-            tmp = erf(sqrt(r2*gam)) / sqrt(r2)
+            gam = 1.0_wp/(self%rad(izp)**2 + self%rad(jzp)**2)
+            tmp = erf(sqrt(r2*gam))/sqrt(r2)
             amat_local(jat, iat) = amat_local(jat, iat) + tmp
             amat_local(iat, jat) = amat_local(iat, jat) + tmp
          end do
-         tmp = self%eta(izp) + sqrt2pi / self%rad(izp)
+         tmp = self%eta(izp) + sqrt2pi/self%rad(izp)
          amat_local(iat, iat) = amat_local(iat, iat) + tmp
       end do
       !$omp end do
       !$omp critical (get_amat_0d_)
       amat(:, :) = amat(:, :) + amat_local(:, :)
       !$omp end critical (get_amat_0d_)
-      deallocate(amat_local)
+      deallocate (amat_local)
       !$omp end parallel
 
       amat(mol%nat + 1, 1:mol%nat + 1) = 1.0_wp
@@ -265,40 +264,40 @@ subroutine get_amat_3d(self, mol, wsc, alpha, amat)
       !$omp parallel default(none) &
       !$omp shared(amat, mol, self, wsc, dtrans, rtrans, alpha, vol) &
       !$omp private(iat, izp, jat, jzp, gam, wsw, vec, dtmp, rtmp, amat_local)
-      allocate(amat_local, source=amat)
-      !$omp do schedule(runtime) 
+      allocate (amat_local, source=amat)
+      !$omp do schedule(runtime)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          do jat = 1, iat - 1
             jzp = mol%id(jat)
-            gam = 1.0_wp / sqrt(self%rad(izp)**2 + self%rad(jzp)**2)
-            wsw = 1.0_wp / real(wsc%nimg(jat, iat), wp)
+            gam = 1.0_wp/sqrt(self%rad(izp)**2 + self%rad(jzp)**2)
+            wsw = 1.0_wp/real(wsc%nimg(jat, iat), wp)
             do img = 1, wsc%nimg(jat, iat)
-               vec = mol%xyz(:, iat) - mol%xyz(:, jat) - wsc%trans(:, wsc%tridx(img, jat, iat))
+               vec = mol%xyz(:, jat) - mol%xyz(:, iat) + wsc%trans(:, wsc%tridx(img, jat, iat))
                call get_amat_dir_3d(vec, gam, alpha, dtrans, dtmp)
                call get_amat_rec_3d(vec, vol, alpha, rtrans, rtmp)
-               amat_local(jat, iat) = amat_local(jat, iat) + (dtmp + rtmp) * wsw
-               amat_local(iat, jat) = amat_local(iat, jat) + (dtmp + rtmp) * wsw
+               amat_local(jat, iat) = amat_local(jat, iat) + (dtmp + rtmp)*wsw
+               amat_local(iat, jat) = amat_local(iat, jat) + (dtmp + rtmp)*wsw
             end do
          end do
 
-         gam = 1.0_wp / sqrt(2.0_wp * self%rad(izp)**2)
-         wsw = 1.0_wp / real(wsc%nimg(iat, iat), wp)
+         gam = 1.0_wp/sqrt(2.0_wp*self%rad(izp)**2)
+         wsw = 1.0_wp/real(wsc%nimg(iat, iat), wp)
          do img = 1, wsc%nimg(iat, iat)
             vec = wsc%trans(:, wsc%tridx(img, iat, iat))
             call get_amat_dir_3d(vec, gam, alpha, dtrans, dtmp)
             call get_amat_rec_3d(vec, vol, alpha, rtrans, rtmp)
-            amat_local(iat, iat) = amat_local(iat, iat) + (dtmp + rtmp) * wsw
+            amat_local(iat, iat) = amat_local(iat, iat) + (dtmp + rtmp)*wsw
          end do
 
-         dtmp = self%eta(izp) + sqrt2pi / self%rad(izp) - 2 * alpha / sqrtpi
+         dtmp = self%eta(izp) + sqrt2pi/self%rad(izp) - 2*alpha/sqrtpi
          amat_local(iat, iat) = amat_local(iat, iat) + dtmp
       end do
       !$omp end do
       !$omp critical (get_amat_3d_)
       amat(:, :) = amat(:, :) + amat_local(:, :)
       !$omp end critical (get_amat_3d_)
-      deallocate(amat_local)
+      deallocate (amat_local)
       !$omp end parallel
 
       amat(mol%nat + 1, 1:mol%nat + 1) = 1.0_wp
@@ -364,9 +363,9 @@ subroutine get_coulomb_derivs(self, mol, cache, qvec, dadr, dadL, atrace)
       call view(cache, ptr)
 
       if (any(mol%periodic)) then
-         call self%get_damat_3d(mol, ptr%wsc, ptr%alpha, qvec, dadr, dadL, atrace)
+         call get_damat_3d(self, mol, ptr%wsc, ptr%alpha, qvec, dadr, dadL, atrace)
       else
-         call self%get_damat_0d(mol, qvec, dadr, dadL, atrace)
+         call get_damat_0d(self, mol, qvec, dadr, dadL, atrace)
       end if
    end subroutine get_coulomb_derivs
 
@@ -393,25 +392,25 @@ subroutine get_damat_0d(self, mol, qvec, dadr, dadL, atrace)
       !$omp shared(atrace, dadr, dadL, mol, self, qvec) &
       !$omp private(iat, izp, jat, jzp, gam, r2, vec, dG, dS, dtmp, arg) &
       !$omp private(atrace_local, dadr_local, dadL_local)
-      allocate(atrace_local, source=atrace)
-      allocate(dadr_local, source=dadr)
-      allocate(dadL_local, source=dadL)
+      allocate (atrace_local, source=atrace)
+      allocate (dadr_local, source=dadr)
+      allocate (dadL_local, source=dadL)
       !$omp do schedule(runtime)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          do jat = 1, iat - 1
             jzp = mol%id(jat)
-            vec = mol%xyz(:, iat) - mol%xyz(:, jat)
+            vec = mol%xyz(:, jat) - mol%xyz(:, iat)
             r2 = vec(1)**2 + vec(2)**2 + vec(3)**2
             gam = 1.0_wp/sqrt(self%rad(izp)**2 + self%rad(jzp)**2)
             arg = gam*gam*r2
-            dtmp = 2.0_wp*gam*exp(-arg)/(sqrtpi*r2)-erf(sqrt(arg))/(r2*sqrt(r2))
+            dtmp = 2.0_wp*gam*exp(-arg)/(sqrtpi*r2) - erf(sqrt(arg))/(r2*sqrt(r2))
             dG = dtmp*vec
-            dS = spread(dG, 1, 3) * spread(vec, 2, 3)
-            atrace_local(:, iat) = +dG*qvec(jat) + atrace_local(:, iat)
-            atrace_local(:, jat) = -dG*qvec(iat) + atrace_local(:, jat)
-            dadr_local(:, iat, jat) = +dG*qvec(iat)
-            dadr_local(:, jat, iat) = -dG*qvec(jat)
+            dS = spread(dG, 1, 3)*spread(vec, 2, 3)
+            atrace_local(:, iat) = -dG*qvec(jat) + atrace_local(:, iat)
+            atrace_local(:, jat) = +dG*qvec(iat) + atrace_local(:, jat)
+            dadr_local(:, iat, jat) = -dG*qvec(iat)
+            dadr_local(:, jat, iat) = +dG*qvec(jat)
             dadL_local(:, :, jat) = +dS*qvec(iat) + dadL_local(:, :, jat)
             dadL_local(:, :, iat) = +dS*qvec(jat) + dadL_local(:, :, iat)
          end do
@@ -422,9 +421,9 @@ subroutine get_damat_0d(self, mol, qvec, dadr, dadL, atrace)
       dadr(:, :, :) = dadr(:, :, :) + dadr_local(:, :, :)
       dadL(:, :, :) = dadL(:, :, :) + dadL_local(:, :, :)
       !$omp end critical (get_damat_0d_)
-      deallocate(dadL_local, dadr_local, atrace_local)
+      deallocate (dadL_local, dadr_local, atrace_local)
       !$omp end parallel
-   
+
    end subroutine get_damat_0d
 
    subroutine get_damat_3d(self, mol, wsc, alpha, qvec, dadr, dadL, atrace)
@@ -457,43 +456,43 @@ subroutine get_damat_3d(self, mol, wsc, alpha, qvec, dadr, dadL, atrace)
       !$omp parallel default(none) &
       !$omp shared(mol, self, wsc, alpha, vol, dtrans, rtrans, qvec) &
       !$omp shared(atrace, dadr, dadL) &
-      !$omp private(iat, izp, jat, jzp, img, gam, wsw, vec, dG, dS) & 
+      !$omp private(iat, izp, jat, jzp, img, gam, wsw, vec, dG, dS) &
       !$omp private(dGr, dSr, dGd, dSd, atrace_local, dadr_local, dadL_local)
-      allocate(atrace_local, source=atrace)
-      allocate(dadr_local, source=dadr)
-      allocate(dadL_local, source=dadL)
+      allocate (atrace_local, source=atrace)
+      allocate (dadr_local, source=dadr)
+      allocate (dadL_local, source=dadL)
       !$omp do schedule(runtime)
       do iat = 1, mol%nat
          izp = mol%id(iat)
          do jat = 1, iat - 1
             jzp = mol%id(jat)
             dG(:) = 0.0_wp
             dS(:, :) = 0.0_wp
-            gam = 1.0_wp / sqrt(self%rad(izp)**2 + self%rad(jzp)**2)
-            wsw = 1.0_wp / real(wsc%nimg(jat, iat), wp)
+            gam = 1.0_wp/sqrt(self%rad(izp)**2 + self%rad(jzp)**2)
+            wsw = 1.0_wp/real(wsc%nimg(jat, iat), wp)
             do img = 1, wsc%nimg(jat, iat)
-               vec = mol%xyz(:, iat) - mol%xyz(:, jat) - wsc%trans(:, wsc%tridx(img, jat, iat))
+               vec = mol%xyz(:, jat) - mol%xyz(:, iat) + wsc%trans(:, wsc%tridx(img, jat, iat))
                call get_damat_dir_3d(vec, gam, alpha, dtrans, dGd, dSd)
                call get_damat_rec_3d(vec, vol, alpha, rtrans, dGr, dSr)
-               dG = dG + (dGd + dGr) * wsw
-               dS = dS + (dSd + dSr) * wsw
+               dG = dG + (dGd + dGr)*wsw
+               dS = dS + (dSd + dSr)*wsw
             end do
-            atrace_local(:, iat) = +dG*qvec(jat) + atrace_local(:, iat)
-            atrace_local(:, jat) = -dG*qvec(iat) + atrace_local(:, jat)
-            dadr_local(:, iat, jat) = +dG*qvec(iat) + dadr_local(:, iat, jat)
-            dadr_local(:, jat, iat) = -dG*qvec(jat) + dadr_local(:, jat, iat)
+            atrace_local(:, iat) = -dG*qvec(jat) + atrace_local(:, iat)
+            atrace_local(:, jat) = +dG*qvec(iat) + atrace_local(:, jat)
+            dadr_local(:, iat, jat) = -dG*qvec(iat) + dadr_local(:, iat, jat)
+            dadr_local(:, jat, iat) = +dG*qvec(jat) + dadr_local(:, jat, iat)
             dadL_local(:, :, jat) = +dS*qvec(iat) + dadL_local(:, :, jat)
             dadL_local(:, :, iat) = +dS*qvec(jat) + dadL_local(:, :, iat)
          end do
 
          dS(:, :) = 0.0_wp
-         gam = 1.0_wp / sqrt(2.0_wp * self%rad(izp)**2)
-         wsw = 1.0_wp / real(wsc%nimg(iat, iat), wp)
+         gam = 1.0_wp/sqrt(2.0_wp*self%rad(izp)**2)
+         wsw = 1.0_wp/real(wsc%nimg(iat, iat), wp)
          do img = 1, wsc%nimg(iat, iat)
             vec = wsc%trans(:, wsc%tridx(img, iat, iat))
             call get_damat_dir_3d(vec, gam, alpha, dtrans, dGd, dSd)
             call get_damat_rec_3d(vec, vol, alpha, rtrans, dGr, dSr)
-            dS = dS + (dSd + dSr) * wsw
+            dS = dS + (dSd + dSr)*wsw
          end do
          dadL_local(:, :, iat) = +dS*qvec(iat) + dadL_local(:, :, iat)
       end do
@@ -503,7 +502,7 @@ subroutine get_damat_3d(self, mol, wsc, alpha, qvec, dadr, dadL, atrace)
       dadr(:, :, :) = dadr(:, :, :) + dadr_local(:, :, :)
       dadL(:, :, :) = dadL(:, :, :) + dadL_local(:, :, :)
       !$omp end critical (get_damat_3d_)
-      deallocate(dadL_local, dadr_local, atrace_local)
+      deallocate (dadL_local, dadr_local, atrace_local)
       !$omp end parallel
 
    end subroutine get_damat_3d