diff --git a/app/main.f90 b/app/main.f90 index 0bd4542b..013ef744 100644 --- a/app/main.f90 +++ b/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 diff --git a/src/multicharge/cache.f90 b/src/multicharge/cache.f90 index 6de222e7..a73fbb79 100644 --- a/src/multicharge/cache.f90 +++ b/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 diff --git a/src/multicharge/model/eeq.f90 b/src/multicharge/model/eeq.f90 index 5283f8b6..0736a92f 100644 --- a/src/multicharge/model/eeq.f90 +++ b/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,7 +108,6 @@ 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 @@ -123,6 +115,13 @@ subroutine update(self, mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) 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,11 +456,11 @@ 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) @@ -469,31 +468,31 @@ subroutine get_damat_3d(self, mol, wsc, alpha, qvec, dadr, dadL, atrace) 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 diff --git a/src/multicharge/model/eeqbc.f90 b/src/multicharge/model/eeqbc.f90 index f1de07bb..dc838536 100644 --- a/src/multicharge/model/eeqbc.f90 +++ b/src/multicharge/model/eeqbc.f90 @@ -17,22 +17,17 @@ !> Provides implementation of the bond capacitor electronegativity equilibration model (EEQ_BC) !> Bond capacitor electronegativity equilibration charge model published in -!> -!> Thomas Froitzheim, Marcel Müller, Andreas Hansen, and Stefan Grimme, +!> +!> Thomas Froitzheim, Marcel Müller, Andreas Hansen, and Stefan Grimme, !> *J. Chem. Phys.*, **2025**, 162, 214109. !> DOI: [10.1063/5.0268978](https://dx.doi.org/10.1063/5.0268978) module multicharge_model_eeqbc - - use iso_fortran_env, only: output_unit - use mctc_env, only: error_type, wp use mctc_io, only: structure_type use mctc_io_constants, only: pi - use mctc_io_convert, only: autoaa - use mctc_io_math, only: matdet_3x3 use mctc_ncoord, only: new_ncoord, cn_count - use multicharge_wignerseitz, only: wignerseitz_cell_type - use multicharge_model_type, only: mchrg_model_type, get_dir_trans, get_rec_trans + use multicharge_wignerseitz, only: new_wignerseitz_cell, wignerseitz_cell_type + use multicharge_model_type, only: mchrg_model_type, get_dir_trans use multicharge_blas, only: gemv, gemm use multicharge_model_cache, only: cache_container, model_cache implicit none @@ -48,10 +43,8 @@ module multicharge_model_eeqbc real(wp), allocatable :: dqlocdr(:, :, :) !> Local charge dL derivative real(wp), allocatable :: dqlocdL(:, :, :) - !> Full Maxwell capacitance matrix for 0d case + !> Full Maxwell capacitance matrix real(wp), allocatable :: cmat(:, :) - !> Diagonal elements of Maxwell capacitance matrix for every WSC image - real(wp), allocatable :: cdiag(:, :) !> Derivative of Maxwell capacitance matrix w.r.t positions real(wp), allocatable :: dcdr(:, :, :) !> Derivative of Maxwell capacitance matrix w.r.t lattice vectors @@ -82,21 +75,14 @@ module multicharge_model_eeqbc procedure :: get_xvec !> Calculate derivatives of EN vector 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 - !> Calculate constraint matrix (molecular case) + !> Calculate constraint matrix (molecular) procedure :: get_cmat_0d - !> Calculate diagonal contributions (periodic case) - procedure :: get_cdiag_3d + !> Calculate full constraint matrix (periodic) + procedure :: get_cmat_3d !> Calculate constraint matrix derivatives (molecular) procedure :: get_dcmat_0d - ! procedure :: get_dcmat_3d + !> Calculate constraint matrix derivatives (periodic) + procedure :: get_dcmat_3d end type eeqbc_model real(wp), parameter :: sqrtpi = sqrt(pi) @@ -111,7 +97,7 @@ module multicharge_model_eeqbc contains subroutine new_eeqbc_model(self, mol, error, chi, rad, & - & eta, kcnchi, kqchi, kqeta, kcnrad, cap, avg_cn, & + & eta, kcnchi, kqchi, kqeta, kcnrad, cap, avg_cn, & & kbc, cutoff, cn_exp, rcov, en, cn_max, norm_exp, rvdw) !> Bond capacitor electronegativity equilibration model type(eeqbc_model), intent(out) :: self @@ -182,7 +168,7 @@ subroutine new_eeqbc_model(self, mol, error, chi, rad, & & cutoff=cutoff, kcn=cn_exp, rcov=rcov, cut=cn_max, & & norm_exp=self%norm_exp) ! Electronegativity weighted coordination number for local charge - call new_ncoord(self%ncoord_en, mol, cn_count%erf_en, error, & + call new_ncoord(self%ncoord_en, mol, cn_count%erf_en, error, & & cutoff=cutoff, kcn=cn_exp, rcov=rcov, en=en, cut=cn_max, & & norm_exp=self%norm_exp) @@ -204,7 +190,6 @@ subroutine update(self, mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) type(eeqbc_cache), pointer :: ptr call taint(cache, ptr) - call ptr%update(mol) grad = present(dcndr) .and. present(dcndL) .and. present(dqlocdr) .and. present(dqlocdL) @@ -228,29 +213,32 @@ subroutine update(self, mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) allocate (ptr%xtmp(mol%nat + 1)) end if + ! Allocate cmat + if (.not. allocated(ptr%cmat)) then + allocate (ptr%cmat(mol%nat + 1, mol%nat + 1)) + end if + if (any(mol%periodic)) then - ! Allocate cmat diagonal WSC image contributions - if (.not. allocated(ptr%cdiag)) then - allocate (ptr%cdiag(mol%nat, ptr%wsc%nimg_max)) + ! Create WSC + call new_wignerseitz_cell(ptr%wsc, mol) + + ! Get full cmat sum over all WSC images (for get_xvec and xvec_derivs) + call get_cmat_3d(self, mol, ptr%wsc, ptr%cmat) + if (grad) then + if (.not. allocated(ptr%dcdr) .and. .not. allocated(ptr%dcdL)) then + allocate (ptr%dcdr(3, mol%nat, mol%nat + 1), ptr%dcdL(3, 3, mol%nat + 1)) + end if + call get_dcmat_3d(self, mol, ptr%wsc, ptr%dcdr, ptr%dcdL) end if - ! Get cmat diagonal contributions for all WSC images - call self%get_cdiag_3d(mol, ptr%wsc, ptr%cdiag) - ! if (grad) then - ! call self%get_dcmat_3d() - ! end if else - ! Allocate cmat - if (.not. allocated(ptr%cmat)) then - allocate (ptr%cmat(mol%nat + 1, mol%nat + 1)) - end if - call self%get_cmat_0d(mol, ptr%cmat) + call get_cmat_0d(self, mol, ptr%cmat) ! cmat gradients if (grad) then if (.not. allocated(ptr%dcdr) .and. .not. allocated(ptr%dcdL)) then allocate (ptr%dcdr(3, mol%nat, mol%nat + 1), ptr%dcdL(3, 3, mol%nat + 1)) end if - call self%get_dcmat_0d(mol, ptr%dcdr, ptr%dcdL) + call get_dcmat_0d(self, mol, ptr%dcdr, ptr%dcdL) end if end if @@ -264,37 +252,74 @@ subroutine get_xvec(self, mol, cache, xvec) type(eeqbc_cache), pointer :: ptr - integer :: iat, izp + integer :: iat, izp, img + real(wp) :: ctmp, vec(3), rvdw, capi, wsw + real(wp), allocatable :: dtrans(:, :) + + ! Thread-private array for reduction + real(wp), allocatable :: xvec_local(:) call view(cache, ptr) + xvec(:) = 0.0_wp !$omp parallel do default(none) schedule(runtime) & - !$omp shared(mol, self, ptr) private(iat, izp) + !$omp shared(mol, self, ptr, xvec) & + !$omp private(iat, izp) do iat = 1, mol%nat izp = mol%id(iat) ptr%xtmp(iat) = -self%chi(izp) + self%kcnchi(izp)*ptr%cn(iat) & & + self%kqchi(izp)*ptr%qloc(iat) end do ptr%xtmp(mol%nat + 1) = mol%charge + call gemv(ptr%cmat, ptr%xtmp, xvec) + if (any(mol%periodic)) then + call get_dir_trans(mol%lattice, dtrans) + !$omp parallel default(none) & + !$omp shared(mol, self, ptr, xvec, dtrans) private(iat, izp, img, wsw) & + !$omp private(capi, vec, rvdw, ctmp, xvec_local) + allocate (xvec_local, mold=xvec) + xvec_local(:) = 0.0_wp + !$omp do schedule(runtime) + do iat = 1, mol%nat + izp = mol%id(iat) + capi = self%cap(izp) + ! eliminate self-interaction (quasi off-diagonal) + rvdw = self%rvdw(iat, iat) + wsw = 1.0_wp/real(ptr%wsc%nimg(iat, iat), wp) + do img = 1, ptr%wsc%nimg(iat, iat) + vec = ptr%wsc%trans(:, ptr%wsc%tridx(img, iat, iat)) + + call get_cpair_dir(self%kbc, vec, dtrans, rvdw, capi, capi, ctmp) + xvec_local(iat) = xvec_local(iat) - wsw*ctmp*ptr%xtmp(iat) + end do + end do + !$omp end do + !$omp critical (get_xvec_) + xvec(:) = xvec + xvec_local + !$omp end critical (get_xvec_) + deallocate (xvec_local) + !$omp end parallel + end if end subroutine get_xvec subroutine get_xvec_derivs(self, mol, cache, dxdr, dxdL) class(eeqbc_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(:, :, :) type(eeqbc_cache), pointer :: ptr - integer :: iat, izp, jat - real(wp) :: tmp(3) + integer :: iat, izp, jat, jzp, img + real(wp) :: capi, capj, wsw, vec(3), ctmp, rvdw, dG(3), dS(3, 3) real(wp), allocatable :: dtmpdr(:, :, :), dtmpdL(:, :, :) + real(wp), allocatable :: dtrans(:, :) ! Thread-private arrays for reduction - real(wp), allocatable :: dxdr_local(:, :, :), dxdL_local(:, :, :) + real(wp), allocatable :: dxdr_local(:, :, :), dxdL_local(:, :, :), dtmpdr_local(:, :, :), dtmpdL_local(:, :, :) call view(cache, ptr) allocate (dtmpdr(3, mol%nat, mol%nat + 1), dtmpdL(3, 3, mol%nat + 1)) @@ -304,35 +329,123 @@ subroutine get_xvec_derivs(self, mol, cache, dxdr, dxdL) dtmpdr(:, :, :) = 0.0_wp dtmpdL(:, :, :) = 0.0_wp - !$omp parallel do default(none) schedule(runtime) & + !$omp parallel default(none) & !$omp shared(mol, self, ptr, dtmpdr, dtmpdL) & - !$omp private(iat, izp) + !$omp private(iat, izp, dtmpdr_local, dtmpdL_local) + allocate (dtmpdr_local, source=dtmpdr) + allocate (dtmpdL_local, source=dtmpdL) + !$omp do schedule(runtime) do iat = 1, mol%nat izp = mol%id(iat) ! CN and effective charge derivative - dtmpdr(:, :, iat) = self%kcnchi(izp)*ptr%dcndr(:, :, iat) + dtmpdr(:, :, iat) - dtmpdL(:, :, iat) = self%kcnchi(izp)*ptr%dcndL(:, :, iat) + dtmpdL(:, :, iat) - dtmpdr(:, :, iat) = self%kqchi(izp)*ptr%dqlocdr(:, :, iat) + dtmpdr(:, :, iat) - dtmpdL(:, :, iat) = self%kqchi(izp)*ptr%dqlocdL(:, :, iat) + dtmpdL(:, :, iat) + dtmpdr_local(:, :, iat) = self%kcnchi(izp)*ptr%dcndr(:, :, iat) + dtmpdr_local(:, :, iat) + dtmpdL_local(:, :, iat) = self%kcnchi(izp)*ptr%dcndL(:, :, iat) + dtmpdL_local(:, :, iat) + dtmpdr_local(:, :, iat) = self%kqchi(izp)*ptr%dqlocdr(:, :, iat) + dtmpdr_local(:, :, iat) + dtmpdL_local(:, :, iat) = self%kqchi(izp)*ptr%dqlocdL(:, :, iat) + dtmpdL_local(:, :, iat) end do + !$omp end do + !$omp critical (get_xvec_derivs_) + dtmpdr(:, :, :) = dtmpdr + dtmpdr_local + dtmpdL(:, :, :) = dtmpdL + dtmpdL_local + !$omp end critical (get_xvec_derivs_) + deallocate (dtmpdL_local, dtmpdr_local) + !$omp end parallel call gemm(dtmpdr, ptr%cmat, dxdr) call gemm(dtmpdL, ptr%cmat, dxdL) - !$omp parallel do default(none) schedule(runtime) & - !$omp shared(mol, self, ptr, dxdr) & - !$omp private(iat, izp, tmp) - do iat = 1, mol%nat - tmp = 0.0_wp - do jat = 1, mol%nat - ! Diagonal elements - tmp(:) = tmp(:) + ptr%xtmp(jat)*ptr%dcdr(:, iat, jat) - ! Derivative of capacitance matrix - dxdr(:, iat, jat) = (ptr%xtmp(iat) - ptr%xtmp(jat))*ptr%dcdr(:, iat, jat) & - & + dxdr(:, iat, jat) + if (any(mol%periodic)) then + call get_dir_trans(mol%lattice, dtrans) + !$omp parallel default(none) & + !$omp shared(mol, self, ptr, dxdr, dxdL, dtrans) & + !$omp private(iat, izp, jat, jzp, img, wsw) & + !$omp private(capi, capj, vec, rvdw, ctmp, dG, dS) & + !$omp private(dxdr_local, dxdL_local) + allocate (dxdr_local, mold=dxdr) + allocate (dxdL_local, mold=dxdL) + dxdr_local(:, :, :) = 0.0_wp + dxdL_local(:, :, :) = 0.0_wp + !$omp do schedule(runtime) + do iat = 1, mol%nat + izp = mol%id(iat) + capi = self%cap(izp) + do jat = 1, mol%nat + rvdw = self%rvdw(iat, jat) + jzp = mol%id(jat) + capj = self%cap(jzp) + + ! Diagonal elements + dxdr_local(:, iat, iat) = dxdr_local(:, iat, iat) + ptr%xtmp(jat)*ptr%dcdr(:, iat, jat) + + ! Derivative of capacitance matrix + dxdr_local(:, iat, jat) = dxdr_local(:, iat, jat) & + & + (ptr%xtmp(iat) - ptr%xtmp(jat))*ptr%dcdr(:, iat, jat) + + wsw = 1.0_wp/real(ptr%wsc%nimg(iat, jat), wp) + do img = 1, ptr%wsc%nimg(iat, jat) + vec = mol%xyz(:, jat) - mol%xyz(:, iat) + ptr%wsc%trans(:, ptr%wsc%tridx(img, jat, iat)) + call get_dcpair_dir(self%kbc, vec, dtrans, rvdw, capi, capj, dG, dS) + dxdL_local(:, :, iat) = dxdL_local(:, :, iat) + wsw*dS*ptr%xtmp(jat) + end do + end do + dxdL_local(:, :, iat) = dxdL_local(:, :, iat) + ptr%xtmp(iat)*ptr%dcdL(:, :, iat) + + ! Capacitance terms for i = j, T != 0 + rvdw = self%rvdw(iat, iat) + wsw = 1.0_wp/real(ptr%wsc%nimg(iat, iat), wp) + do img = 1, ptr%wsc%nimg(iat, iat) + vec = ptr%wsc%trans(:, ptr%wsc%tridx(img, iat, iat)) + + call get_cpair_dir(self%kbc, vec, dtrans, rvdw, capi, capi, ctmp) + ctmp = ctmp*wsw + ! EN derivative + dxdr_local(:, :, iat) = dxdr_local(:, :, iat) - ctmp*self%kcnchi(izp)*ptr%dcndr(:, :, iat) + dxdL_local(:, :, iat) = dxdL_local(:, :, iat) - ctmp*self%kcnchi(izp)*ptr%dcndL(:, :, iat) + dxdr_local(:, :, iat) = dxdr_local(:, :, iat) - ctmp*self%kqchi(izp)*ptr%dqlocdr(:, :, iat) + dxdL_local(:, :, iat) = dxdL_local(:, :, iat) - ctmp*self%kqchi(izp)*ptr%dqlocdL(:, :, iat) + end do end do - dxdr(:, iat, iat) = dxdr(:, iat, iat) + tmp(:) - end do + !$omp end do + !$omp critical (get_xvec_derivs_update) + dxdr(:, :, :) = dxdr + dxdr_local + dxdL(:, :, :) = dxdL + dxdL_local + !$omp end critical (get_xvec_derivs_update) + deallocate (dxdL_local, dxdr_local) + !$omp end parallel + else + !$omp parallel default(none) & + !$omp shared(mol, self, ptr, dxdr, dxdL) & + !$omp private(iat, izp, jat, jzp, vec, dxdr_local, dxdL_local) + allocate (dxdr_local, mold=dxdr) + allocate (dxdL_local, mold=dxdL) + dxdr_local(:, :, :) = 0.0_wp + dxdL_local(:, :, :) = 0.0_wp + !$omp do schedule(runtime) + do iat = 1, mol%nat + do jat = 1, iat - 1 + ! Diagonal elements + dxdr_local(:, iat, iat) = dxdr_local(:, iat, iat) + ptr%xtmp(jat)*ptr%dcdr(:, iat, jat) + dxdr_local(:, jat, jat) = dxdr_local(:, jat, jat) + ptr%xtmp(iat)*ptr%dcdr(:, jat, iat) + + ! Derivative of capacitance matrix + dxdr_local(:, iat, jat) = (ptr%xtmp(iat) - ptr%xtmp(jat))*ptr%dcdr(:, iat, jat) + dxdr_local(:, iat, jat) + dxdr_local(:, jat, iat) = (ptr%xtmp(jat) - ptr%xtmp(iat))*ptr%dcdr(:, jat, iat) + dxdr_local(:, jat, iat) + + vec = mol%xyz(:, iat) - mol%xyz(:, jat) + dxdL_local(:, :, iat) = dxdL_local(:, :, iat) + ptr%xtmp(jat)*spread(ptr%dcdr(:, iat, jat), 1, 3)*spread(vec, 2, 3) + dxdL_local(:, :, jat) = dxdL_local(:, :, jat) + ptr%xtmp(iat)*spread(ptr%dcdr(:, jat, iat), 1, 3)*spread(-vec, 2, 3) + end do + dxdr_local(:, iat, iat) = dxdr_local(:, iat, iat) + ptr%xtmp(iat)*ptr%dcdr(:, iat, iat) + dxdL_local(:, :, iat) = dxdL_local(:, :, iat) + ptr%xtmp(iat)*ptr%dcdL(:, :, iat) + end do + !$omp end do + !$omp critical (get_xvec_derivs_) + dxdr(:, :, :) = dxdr + dxdr_local + dxdL(:, :, :) = dxdL + dxdL_local + !$omp end critical (get_xvec_derivs_) + deallocate (dxdL_local, dxdr_local) + !$omp end parallel + end if end subroutine get_xvec_derivs @@ -346,10 +459,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, ptr%cn, & - & ptr%qloc, amat, ptr%cdiag) + call get_amat_3d(self, mol, ptr%wsc, ptr%cn, ptr%qloc, ptr%cmat, amat) else - call self%get_amat_0d(mol, ptr%cn, ptr%qloc, ptr%cmat, amat) + call get_amat_0d(self, mol, ptr%cn, ptr%qloc, ptr%cmat, amat) end if end subroutine get_coulomb_matrix @@ -366,42 +478,42 @@ subroutine get_amat_0d(self, mol, cn, qloc, cmat, amat) ! Thread-private array for reduction real(wp), allocatable :: amat_local(:, :) - + amat(:, :) = 0.0_wp - + !$omp parallel default(none) & !$omp shared(amat, mol, self, cn, qloc, cmat) & !$omp private(iat, izp, jat, jzp, gam2, vec, r2, tmp) & !$omp private(norm_cn, radi, radj, 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) ! Effective charge width of i - norm_cn = 1.0_wp / self%avg_cn(izp)**self%norm_exp - radi = self%rad(izp) * (1.0_wp - self%kcnrad*cn(iat)*norm_cn) + norm_cn = 1.0_wp/self%avg_cn(izp)**self%norm_exp + radi = self%rad(izp)*(1.0_wp - self%kcnrad*cn(iat)*norm_cn) 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 ! Effective charge width of j - norm_cn = cn(jat) / self%avg_cn(jzp)**self%norm_exp - radj = self%rad(jzp) * (1.0_wp - self%kcnrad*norm_cn) + norm_cn = cn(jat)/self%avg_cn(jzp)**self%norm_exp + radj = self%rad(jzp)*(1.0_wp - self%kcnrad*norm_cn) ! Coulomb interaction of Gaussian charges - gam2 = 1.0_wp / (radi**2 + radj**2) - tmp = erf(sqrt(r2*gam2)) / sqrt(r2) * cmat(jat, iat) - amat_local(jat, iat) = amat_local(jat, iat) + tmp - amat_local(iat, jat) = amat_local(iat, jat) + tmp + gam2 = 1.0_wp/(radi**2 + radj**2) + tmp = erf(sqrt(r2*gam2))/sqrt(r2)*cmat(jat, iat) + amat_local(jat, iat) = tmp + amat_local(iat, jat) = tmp end do ! Effective hardness - tmp = self%eta(izp) + self%kqeta(izp) * qloc(iat) + sqrt2pi / radi + tmp = self%eta(izp) + self%kqeta(izp)*qloc(iat) + sqrt2pi/radi amat_local(iat, iat) = amat_local(iat, iat) + tmp*cmat(iat, iat) + 1.0_wp end do !$omp end do !$omp critical (get_amat_0d_) - amat(:, :) = amat(:, :) + amat_local(:, :) + 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 @@ -410,82 +522,74 @@ subroutine get_amat_0d(self, mol, cn, qloc, cmat, amat) end subroutine get_amat_0d - subroutine get_amat_3d(self, mol, wsc, alpha, cn, qloc, amat, cdiag) + 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 - real(wp), intent(in) :: alpha - real(wp), intent(in) :: cn(:), qloc(:) + real(wp), intent(in) :: cn(:), qloc(:), cmat(:, :) real(wp), intent(out) :: amat(:, :) - real(wp), intent(out) :: cdiag(:, :) - integer :: iat, jat, isp, jsp, izp, jzp, img - real(wp) :: vec(3), gam, wsw, dtmp, rtmp, vol, ctmp, capi, capj, radi, radj, norm_cn, rvdw - real(wp), allocatable :: dtrans(:, :), rtrans(:, :) + integer :: iat, jat, izp, jzp, img + real(wp) :: vec(3), r1, gam, dtmp, ctmp, capi, capj, radi, radj, norm_cn, rvdw, wsw + real(wp), allocatable :: dtrans(:, :) ! Thread-private array for reduction real(wp), allocatable :: amat_local(:, :) - amat(:, :) = 0.0_wp - - vol = abs(matdet_3x3(mol%lattice)) call get_dir_trans(mol%lattice, dtrans) - call get_rec_trans(mol%lattice, rtrans) + + amat(:, :) = 0.0_wp !$omp parallel default(none) & - !$omp shared(amat, mol, cn, qloc, self, wsc, dtrans, rtrans, alpha, vol, cdiag) & - !$omp private(iat, izp, jat, jzp, gam, wsw, vec, dtmp, rtmp, ctmp, norm_cn) & - !$omp private(isp, jsp, radi, radj, capi, capj, rvdw, amat_local) - allocate(amat_local, source=amat) - !$omp do schedule(runtime) + !$omp shared(amat, cmat, mol, cn, qloc, self, wsc, dtrans) & + !$omp private(iat, izp, jat, jzp, gam, vec, dtmp, ctmp, norm_cn) & + !$omp private(radi, radj, capi, capj, rvdw, r1, wsw, amat_local) + allocate (amat_local, source=amat) + !$omp do schedule(runtime) do iat = 1, mol%nat izp = mol%id(iat) - isp = mol%num(izp) ! Effective charge width of i - norm_cn = 1.0_wp / self%avg_cn(izp)**self%norm_exp - radi = self%rad(izp) * (1.0_wp - self%kcnrad*cn(iat)*norm_cn) - capi = self%cap(isp) + norm_cn = 1.0_wp/self%avg_cn(izp)**self%norm_exp + radi = self%rad(izp)*(1.0_wp - self%kcnrad*cn(iat)*norm_cn) + capi = self%cap(izp) do jat = 1, iat - 1 jzp = mol%id(jat) - jsp = mol%num(jzp) ! vdw distance in Angstrom (approximate factor 2) rvdw = self%rvdw(iat, jat) ! Effective charge width of j - norm_cn = cn(jat) / self%avg_cn(jzp)**self%norm_exp - radj = self%rad(jzp) * (1.0_wp - self%kcnrad*norm_cn) - capj = self%cap(jsp) + norm_cn = cn(jat)/self%avg_cn(jzp)**self%norm_exp + radj = self%rad(jzp)*(1.0_wp - self%kcnrad*norm_cn) + capj = self%cap(jzp) ! Coulomb interaction of Gaussian charges - gam = 1.0_wp / sqrt(radi**2 + radj**2) - wsw = 1.0_wp / real(wsc%nimg(jat, iat), wp) + gam = 1.0_wp/sqrt(radi**2 + radj**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)) - call get_cpair(self%kbc, ctmp, vec, rvdw, capi, capj) - 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) + ctmp*(dtmp + rtmp)*wsw - amat_local(iat, jat) = amat_local(iat, jat) + ctmp*(dtmp + rtmp)*wsw + vec = mol%xyz(:, jat) - mol%xyz(:, iat) + wsc%trans(:, wsc%tridx(img, jat, iat)) + call get_amat_dir_3d(vec, gam, dtrans, self%kbc, rvdw, capi, capj, dtmp) + amat_local(jat, iat) = amat_local(jat, iat) + dtmp*wsw + amat_local(iat, jat) = amat_local(iat, jat) + dtmp*wsw end do end do - ! WSC image contributions - gam = 1.0_wp / sqrt(2.0_wp * self%rad(izp)**2) - wsw = 1.0_wp / real(wsc%nimg(iat, iat), wp) + ! diagonal Coulomb interaction terms + gam = 1.0_wp/sqrt(2.0_wp*radi**2) + rvdw = self%rvdw(iat, iat) + 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)) - ctmp = cdiag(iat, img) - 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) + ctmp*(dtmp + rtmp)*wsw + call get_amat_dir_3d(vec, gam, dtrans, self%kbc, rvdw, capi, capi, dtmp) + amat_local(iat, iat) = amat_local(iat, iat) + dtmp*wsw end do + ! Effective hardness - dtmp = self%eta(izp) + self%kqeta(izp) * qloc(iat) + sqrt2pi / radi - amat_local(iat, iat) = amat_local(iat, iat) + cdiag(iat, 1) * dtmp + 1.0_wp + dtmp = self%eta(izp) + self%kqeta(izp)*qloc(iat) + sqrt2pi/radi + amat_local(iat, iat) = amat_local(iat, iat) + cmat(iat, iat)*dtmp + 1.0_wp end do !$omp end do !$omp critical (get_amat_3d_) - amat(:, :) = amat(:, :) + amat_local(:, :) + 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 @@ -494,15 +598,18 @@ subroutine get_amat_3d(self, mol, wsc, alpha, cn, qloc, amat, cdiag) end subroutine get_amat_3d - subroutine get_amat_dir_3d(rij, gam, alp, trans, amat) + subroutine get_amat_dir_3d(rij, gam, trans, kbc, rvdw, capi, capj, amat) real(wp), intent(in) :: rij(3) real(wp), intent(in) :: gam - real(wp), intent(in) :: alp + real(wp), intent(in) :: kbc + real(wp), intent(in) :: rvdw + real(wp), intent(in) :: capi + real(wp), intent(in) :: capj real(wp), intent(in) :: trans(:, :) real(wp), intent(out) :: amat integer :: itr - real(wp) :: vec(3), r1, tmp + real(wp) :: vec(3), r1, tmp, ctmp amat = 0.0_wp @@ -510,35 +617,13 @@ subroutine get_amat_dir_3d(rij, gam, alp, trans, amat) vec(:) = rij + trans(:, itr) r1 = norm2(vec) if (r1 < eps) cycle - tmp = erf(gam*r1)/r1 - erf(alp*r1)/r1 + call get_cpair(kbc, ctmp, r1, rvdw, capi, capj) + tmp = -ctmp*erf(gam*r1)/r1 amat = amat + tmp end do end subroutine get_amat_dir_3d - subroutine get_amat_rec_3d(rij, vol, alp, trans, amat) - real(wp), intent(in) :: rij(3) - real(wp), intent(in) :: vol - real(wp), intent(in) :: alp - real(wp), intent(in) :: trans(:, :) - real(wp), intent(out) :: amat - - integer :: itr - real(wp) :: fac, vec(3), g2, tmp - - amat = 0.0_wp - fac = 4*pi/vol - - do itr = 1, size(trans, 2) - vec(:) = trans(:, itr) - g2 = dot_product(vec, vec) - if (g2 < eps) cycle - tmp = cos(dot_product(rij, vec))*fac*exp(-0.25_wp*g2/(alp*alp))/g2 - amat = amat + tmp - end do - - end subroutine get_amat_rec_3d - subroutine get_coulomb_derivs(self, mol, cache, qvec, dadr, dadL, atrace) class(eeqbc_model), intent(in) :: self type(structure_type), intent(in) :: mol @@ -550,9 +635,12 @@ 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%cn, & + & ptr%qloc, qvec, ptr%dcndr, ptr%dcndL, ptr%dqlocdr, & + & ptr%dqlocdL, ptr%cmat, ptr%dcdr, ptr%dcdL, dadr, dadL, atrace) + else - call self%get_damat_0d(mol, ptr%cn, & + call get_damat_0d(self, mol, ptr%cn, & & ptr%qloc, qvec, ptr%dcndr, ptr%dcndL, ptr%dqlocdr, & & ptr%dqlocdL, ptr%cmat, ptr%dcdr, ptr%dcdL, dadr, dadL, atrace) end if @@ -597,9 +685,9 @@ subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & !$omp private(iat, izp, jat, jzp, gam, vec, r2, dtmp, norm_cn, arg) & !$omp private(radi, radj, dradi, dradj, dgamdr, dgamdL, dG, dS) & !$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) @@ -627,14 +715,14 @@ subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & arg = gam*gam*r2 dtmp = 2.0_wp*gam*exp(-arg)/(sqrtpi*r2) & & - erf(sqrt(arg))/(r2*sqrt(r2)) - dG(:) = -dtmp*vec ! questionable sign + dG(:) = dtmp*vec dS(:, :) = spread(dG, 1, 3)*spread(vec, 2, 3) - atrace_local(:, iat) = +dG*qvec(jat)*cmat(jat, iat) + atrace_local(:, iat) - atrace_local(:, jat) = -dG*qvec(iat)*cmat(iat, jat) + atrace_local(:, jat) - dadr_local(:, iat, jat) = +dG*qvec(iat)*cmat(iat, jat) + dadr_local(:, iat, jat) - dadr_local(:, jat, iat) = -dG*qvec(jat)*cmat(jat, iat) + dadr_local(:, jat, iat) - dadL_local(:, :, jat) = +dS*qvec(iat)*cmat(iat, jat) + dadL_local(:, :, jat) + atrace_local(:, iat) = -dG*qvec(jat)*cmat(jat, iat) + atrace_local(:, iat) + atrace_local(:, jat) = +dG*qvec(iat)*cmat(iat, jat) + atrace_local(:, jat) + dadr_local(:, iat, jat) = -dG*qvec(iat)*cmat(iat, jat) + dadr_local(:, iat, jat) + dadr_local(:, jat, iat) = +dG*qvec(jat)*cmat(jat, iat) + dadr_local(:, jat, iat) dadL_local(:, :, iat) = +dS*qvec(jat)*cmat(jat, iat) + dadL_local(:, :, iat) + dadL_local(:, :, jat) = +dS*qvec(iat)*cmat(iat, jat) + dadL_local(:, :, jat) ! Effective charge width derivative dtmp = 2.0_wp*exp(-arg)/(sqrtpi) @@ -642,229 +730,312 @@ subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & atrace_local(:, jat) = -dtmp*qvec(iat)*dgamdr(:, iat)*cmat(iat, jat) + atrace_local(:, jat) dadr_local(:, iat, jat) = +dtmp*qvec(iat)*dgamdr(:, iat)*cmat(iat, jat) + dadr_local(:, iat, jat) dadr_local(:, jat, iat) = +dtmp*qvec(jat)*dgamdr(:, jat)*cmat(jat, iat) + dadr_local(:, jat, iat) - dadL_local(:, :, jat) = +dtmp*qvec(iat)*dgamdL(:, :)*cmat(iat, jat) + dadL_local(:, :, jat) dadL_local(:, :, iat) = +dtmp*qvec(jat)*dgamdL(:, :)*cmat(jat, iat) + dadL_local(:, :, iat) + dadL_local(:, :, jat) = +dtmp*qvec(iat)*dgamdL(:, :)*cmat(iat, jat) + dadL_local(:, :, jat) ! Capacitance derivative off-diagonal dtmp = erf(sqrt(r2)*gam)/(sqrt(r2)) - ! potentially switch indices for dcdr atrace_local(:, iat) = -dtmp*qvec(jat)*dcdr(:, jat, iat) + atrace_local(:, iat) atrace_local(:, jat) = -dtmp*qvec(iat)*dcdr(:, iat, jat) + atrace_local(:, jat) dadr_local(:, iat, jat) = +dtmp*qvec(iat)*dcdr(:, iat, jat) + dadr_local(:, iat, jat) dadr_local(:, jat, iat) = +dtmp*qvec(jat)*dcdr(:, jat, iat) + dadr_local(:, jat, iat) - dadL_local(:, :, jat) = +dtmp*qvec(iat)*dcdL(:, :, iat) + dadL_local(:, :, jat) - dadL_local(:, :, iat) = +dtmp*qvec(jat)*dcdL(:, :, jat) + dadL_local(:, :, iat) + dadL_local(:, :, iat) = -dtmp*qvec(jat)*spread(dcdr(:, iat, jat), 2, 3)*spread(vec, 1, 3) & + & + dadL_local(:, :, iat) + dadL_local(:, :, jat) = -dtmp*qvec(iat)*spread(dcdr(:, iat, jat), 2, 3)*spread(vec, 1, 3) & + & + dadL_local(:, :, jat) ! Capacitance derivative diagonal dtmp = (self%eta(izp) + self%kqeta(izp)*qloc(iat) + sqrt2pi/radi)*qvec(iat) dadr_local(:, jat, iat) = -dtmp*dcdr(:, jat, iat) + dadr_local(:, jat, iat) + dtmp = (self%eta(jzp) + self%kqeta(jzp)*qloc(jat) + sqrt2pi/radj)*qvec(jat) dadr_local(:, iat, jat) = -dtmp*dcdr(:, iat, jat) + dadr_local(:, iat, jat) end do ! Hardness derivative dtmp = self%kqeta(izp)*qvec(iat)*cmat(iat, iat) - !atrace_local(:, iat) = +dtmp*dqlocdr(:, iat, iat) + atrace_local(:, iat) dadr_local(:, :, iat) = +dtmp*dqlocdr(:, :, iat) + dadr_local(:, :, iat) dadL_local(:, :, iat) = +dtmp*dqlocdL(:, :, iat) + dadL_local(:, :, iat) ! Effective charge width derivative dtmp = -sqrt2pi*dradi/(radi**2)*qvec(iat)*cmat(iat, iat) - !atrace_local(:, iat) = -dtmp*dcndr(:, iat, iat) + atrace_local(:, iat) dadr_local(:, :, iat) = +dtmp*dcndr(:, :, iat) + dadr_local(:, :, iat) dadL_local(:, :, iat) = +dtmp*dcndL(:, :, iat) + dadL_local(:, :, iat) ! Capacitance derivative dtmp = (self%eta(izp) + self%kqeta(izp)*qloc(iat) + sqrt2pi/radi)*qvec(iat) - !atrace_local(:, iat) = -dtmp*dcdr(:, iat, iat) + atrace_local(:, iat) dadr_local(:, iat, iat) = +dtmp*dcdr(:, iat, iat) + dadr_local(:, iat, iat) dadL_local(:, :, iat) = +dtmp*dcdL(:, :, iat) + dadL_local(:, :, iat) end do !$omp end do !$omp critical (get_damat_0d_) - atrace(:, :) = atrace(:, :) + atrace_local(:, :) - dadr(:, :, :) = dadr(:, :, :) + dadr_local(:, :, :) - dadL(:, :, :) = dadL(:, :, :) + dadL_local(:, :, :) + atrace(:, :) = atrace + atrace_local + 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) + 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 - real(wp), intent(in) :: alpha + real(wp), intent(in) :: cn(:) + real(wp), intent(in) :: qloc(:) real(wp), intent(in) :: qvec(:) + real(wp), intent(in) :: dcndr(:, :, :) + real(wp), intent(in) :: dcndL(:, :, :) + real(wp), intent(in) :: dqlocdr(:, :, :) + real(wp), intent(in) :: dqlocdL(:, :, :) + real(wp), intent(in) :: cmat(:, :) + real(wp), intent(in) :: dcdr(:, :, :) + real(wp), intent(in) :: dcdL(:, :, :) real(wp), intent(out) :: dadr(:, :, :) real(wp), intent(out) :: dadL(:, :, :) real(wp), intent(out) :: atrace(:, :) integer :: iat, jat, izp, jzp, img - real(wp) :: vol, gam, wsw, vec(3), dG(3), dS(3, 3) - real(wp) :: dGd(3), dSd(3, 3), dGr(3), dSr(3, 3) - real(wp), allocatable :: dtrans(:, :), rtrans(:, :) + real(wp) :: vec(3), r2, gam, arg, dtmp, norm_cn, rvdw, wsw, dgam + real(wp) :: radi, radj, dradi, dradj, dG(3), dS(3, 3) + real(wp) :: dgamdL(3, 3), capi, capj + real(wp), allocatable :: dgamdr(:, :), dtrans(:, :) ! Thread-private arrays for reduction real(wp), allocatable :: atrace_local(:, :) real(wp), allocatable :: dadr_local(:, :, :), dadL_local(:, :, :) + call get_dir_trans(mol%lattice, dtrans) + + allocate (dgamdr(3, mol%nat)) + atrace(:, :) = 0.0_wp dadr(:, :, :) = 0.0_wp dadL(:, :, :) = 0.0_wp - vol = abs(matdet_3x3(mol%lattice)) - call get_dir_trans(mol%lattice, dtrans) - call get_rec_trans(mol%lattice, rtrans) - !$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(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) + !$omp shared(self, mol, cn, qloc, qvec, wsc, dadr, dadL, atrace) & + !$omp shared (cmat, dcdr, dcdL, dcndr, dcndL, dqlocdr, dqlocdL, dtrans) & + !$omp private(iat, izp, jat, jzp, img, gam, vec, r2, dtmp, norm_cn, arg, rvdw) & + !$omp private(radi, radj, dradi, dradj, capi, capj, dgamdr, dgamdL, dG, dS, wsw) & + !$omp private(dgam, dadr_local, dadL_local, atrace_local) + 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) + ! Effective charge width of i + norm_cn = 1.0_wp/self%avg_cn(izp)**self%norm_exp + radi = self%rad(izp)*(1.0_wp - self%kcnrad*cn(iat)*norm_cn) + dradi = -self%rad(izp)*self%kcnrad*norm_cn + capi = self%cap(izp) 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) + capj = self%cap(jzp) + rvdw = self%rvdw(iat, jat) + + ! Effective charge width of j + norm_cn = 1.0_wp/self%avg_cn(jzp)**self%norm_exp + radj = self%rad(jzp)*(1.0_wp - self%kcnrad*cn(jat)*norm_cn) + dradj = -self%rad(jzp)*self%kcnrad*norm_cn + + ! Coulomb interaction of Gaussian charges + gam = 1.0_wp/sqrt(radi**2 + radj**2) + dgamdr(:, :) = -(radi*dradi*dcndr(:, :, iat) + radj*dradj*dcndr(:, :, jat)) & + & *gam**3.0_wp + dgamdL(:, :) = -(radi*dradi*dcndL(:, :, iat) + radj*dradj*dcndL(:, :, jat)) & + & *gam**3.0_wp + + 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)) - 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 + vec = mol%xyz(:, jat) - mol%xyz(:, iat) + wsc%trans(:, wsc%tridx(img, jat, iat)) + + call get_damat_dir(vec, dtrans, capi, capj, rvdw, self%kbc, gam, dG, dS, dgam) + dG = dG*wsw + dS = dS*wsw + dgam = dgam*wsw + + ! Explicit derivative + 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) + + ! Effective charge width derivative + atrace_local(:, iat) = +dgam*qvec(jat)*dgamdr(:, jat) + atrace_local(:, iat) + atrace_local(:, jat) = +dgam*qvec(iat)*dgamdr(:, iat) + atrace_local(:, jat) + dadr_local(:, iat, jat) = -dgam*qvec(iat)*dgamdr(:, iat) + dadr_local(:, iat, jat) + dadr_local(:, jat, iat) = -dgam*qvec(jat)*dgamdr(:, jat) + dadr_local(:, jat, iat) + dadL_local(:, :, iat) = -dgam*qvec(jat)*dgamdL(:, :) + dadL_local(:, :, iat) + dadL_local(:, :, jat) = -dgam*qvec(iat)*dgamdL(:, :) + dadL_local(:, :, jat) + + call get_damat_dc_dir(vec, dtrans, capi, capj, rvdw, self%kbc, gam, dG, dS) + dG = dG*wsw + dS = dS*wsw + + ! Capacitance derivative off-diagonal + atrace_local(:, iat) = -qvec(jat)*dG(:) + atrace_local(:, iat) + atrace_local(:, jat) = +qvec(iat)*dG(:) + atrace_local(:, jat) + dadr_local(:, jat, iat) = +qvec(jat)*dG(:) + dadr_local(:, jat, iat) + dadr_local(:, iat, jat) = -qvec(iat)*dG(:) + dadr_local(:, iat, jat) + dadL_local(:, :, jat) = +qvec(iat)*dS(:, :) + dadL_local(:, :, jat) + dadL_local(:, :, iat) = +qvec(jat)*dS(:, :) + dadL_local(:, :, iat) + + call get_dcpair_dir(self%kbc, vec, dtrans, rvdw, capi, capj, dG, dS) + dG = dG*wsw + + ! Capacitance derivative diagonal + dtmp = (self%eta(izp) + self%kqeta(izp)*qloc(iat) + sqrt2pi/radi)*qvec(iat) + dadr_local(:, jat, iat) = -dtmp*dG(:) + dadr_local(:, jat, iat) + dtmp = (self%eta(jzp) + self%kqeta(jzp)*qloc(jat) + sqrt2pi/radj)*qvec(jat) + dadr_local(:, iat, jat) = +dtmp*dG(:) + dadr_local(:, iat, jat) 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) - 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) + ! diagonal explicit, charge width, and capacitance derivative terms + gam = 1.0_wp/sqrt(2.0_wp*radi**2) + dtmp = -sqrt2pi*dradi/(radi**2)*qvec(iat) + rvdw = self%rvdw(iat, iat) + 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 + call get_damat_dir(vec, dtrans, capi, capi, rvdw, self%kbc, gam, dG, dS, dgam) + dgam = dgam*wsw + + ! Explicit derivative + dadL_local(:, :, iat) = +dS*wsw*qvec(iat) + dadL_local(:, :, iat) + + ! Effective charge width derivative + atrace_local(:, iat) = +dtmp*dcndr(:, iat, iat)*dgam + atrace_local(:, iat) + dadr_local(:, iat, iat) = -dtmp*dcndr(:, iat, iat)*dgam + dadr_local(:, iat, iat) + dadL_local(:, :, iat) = -dtmp*dcndL(:, :, iat)*dgam + dadL_local(:, :, iat) + + ! Capacitance derivative + call get_damat_dc_dir(vec, dtrans, capi, capi, rvdw, self%kbc, gam, dG, dS) + dadL_local(:, :, iat) = +qvec(iat)*dS*wsw + dadL_local(:, :, iat) end do - dadL_local(:, :, iat) = +dS*qvec(iat) + dadL_local(:, :, iat) + + ! Hardness derivative + dtmp = self%kqeta(izp)*qvec(iat)*cmat(iat, iat) + dadr_local(:, :, iat) = +dtmp*dqlocdr(:, :, iat) + dadr_local(:, :, iat) + dadL_local(:, :, iat) = +dtmp*dqlocdL(:, :, iat) + dadL_local(:, :, iat) + + ! Effective charge width derivative + dtmp = -sqrt2pi*dradi/(radi**2)*qvec(iat)*cmat(iat, iat) + dadr_local(:, :, iat) = +dtmp*dcndr(:, :, iat) + dadr_local(:, :, iat) + dadL_local(:, :, iat) = +dtmp*dcndL(:, :, iat) + dadL_local(:, :, iat) + + dtmp = (self%eta(izp) + self%kqeta(izp)*qloc(iat) + sqrt2pi/radi)*qvec(iat) + dadr_local(:, iat, iat) = +dtmp*dcdr(:, iat, iat) + dadr_local(:, iat, iat) + dadL_local(:, :, iat) = +dtmp*dcdL(:, :, iat) + dadL_local(:, :, iat) + end do !$omp end do !$omp critical (get_damat_3d_) - atrace(:, :) = atrace(:, :) + atrace_local(:, :) - dadr(:, :, :) = dadr(:, :, :) + dadr_local(:, :, :) - dadL(:, :, :) = dadL(:, :, :) + dadL_local(:, :, :) + atrace(:, :) = atrace + atrace_local + 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 - subroutine get_damat_dir_3d(rij, gam, alp, trans, dg, ds) + subroutine get_damat_dir(rij, trans, capi, capj, rvdw, kbc, gam, dG, dS, dgam) real(wp), intent(in) :: rij(3) - real(wp), intent(in) :: gam - real(wp), intent(in) :: alp real(wp), intent(in) :: trans(:, :) - real(wp), intent(out) :: dg(3) - real(wp), intent(out) :: ds(3, 3) + real(wp), intent(in) :: gam + real(wp), intent(in) :: capi, capj, rvdw, kbc + real(wp), intent(out) :: dG(3) + real(wp), intent(out) :: dS(3, 3) + real(wp), intent(out) :: dgam integer :: itr - real(wp) :: vec(3), r1, r2, gtmp, atmp, gam2, alp2 + real(wp) :: vec(3), r1, r2, gtmp, gam2, cmat - dg(:) = 0.0_wp - ds(:, :) = 0.0_wp + dG(:) = 0.0_wp + dS(:, :) = 0.0_wp + dgam = 0.0_wp gam2 = gam*gam - alp2 = alp*alp do itr = 1, size(trans, 2) - vec(:) = rij + trans(:, itr) + vec(:) = rij(:) + trans(:, itr) r1 = norm2(vec) if (r1 < eps) cycle r2 = r1*r1 - gtmp = +2*gam*exp(-r2*gam2)/(sqrtpi*r2) - erf(r1*gam)/(r2*r1) - atmp = -2*alp*exp(-r2*alp2)/(sqrtpi*r2) + erf(r1*alp)/(r2*r1) - dg(:) = dg + (gtmp + atmp)*vec - ds(:, :) = ds + (gtmp + atmp)*spread(vec, 1, 3)*spread(vec, 2, 3) + call get_cpair(kbc, cmat, r1, rvdw, capi, capj) + gtmp = 2.0_wp*gam*exp(-r2*gam2)/(sqrtpi*r2) - erf(r1*gam)/(r2*r1) + dG(:) = dG - cmat*gtmp*vec + dS(:, :) = dS - cmat*gtmp*spread(vec, 1, 3)*spread(vec, 2, 3) + dgam = dgam + cmat*2.0_wp*exp(-gam2*r2)/sqrtpi end do - end subroutine get_damat_dir_3d + end subroutine get_damat_dir - subroutine get_damat_rec_3d(rij, vol, alp, trans, dg, ds) + subroutine get_damat_dc_dir(rij, trans, capi, capj, rvdw, kbc, gam, dG, dS) real(wp), intent(in) :: rij(3) - real(wp), intent(in) :: vol - real(wp), intent(in) :: alp real(wp), intent(in) :: trans(:, :) - real(wp), intent(out) :: dg(3) - real(wp), intent(out) :: ds(3, 3) + real(wp), intent(in) :: gam + real(wp), intent(in) :: capi, capj, rvdw, kbc + real(wp), intent(out) :: dG(3) + real(wp), intent(out) :: dS(3, 3) integer :: itr - real(wp) :: fac, vec(3), g2, gv, etmp, dtmp, alp2 - real(wp), parameter :: unity(3, 3) = reshape(& - & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp) :: vec(3), r1, gtmp(3), stmp(3, 3), tmp - dg(:) = 0.0_wp - ds(:, :) = 0.0_wp - fac = 4*pi/vol - alp2 = alp*alp + dG(:) = 0.0_wp + dS(:, :) = 0.0_wp do itr = 1, size(trans, 2) - vec(:) = trans(:, itr) - g2 = dot_product(vec, vec) - if (g2 < eps) cycle - gv = dot_product(rij, vec) - etmp = fac*exp(-0.25_wp*g2/alp2)/g2 - dtmp = -sin(gv)*etmp - dg(:) = dg + dtmp*vec - ds(:, :) = ds + etmp*cos(gv) & - & *((2.0_wp/g2 + 0.5_wp/alp2)*spread(vec, 1, 3)*spread(vec, 2, 3) - unity) + vec(:) = rij(:) + trans(:, itr) + r1 = norm2(vec) + if (r1 < eps) cycle + call get_dcpair(kbc, vec, rvdw, capi, capj, gtmp, stmp) + tmp = erf(gam*r1)/r1 + dG(:) = dG(:) + tmp*gtmp + dS(:, :) = dS(:, :) + tmp*stmp end do - end subroutine get_damat_rec_3d + end subroutine get_damat_dc_dir subroutine get_cmat_0d(self, mol, cmat) class(eeqbc_model), intent(in) :: self type(structure_type), intent(in) :: mol real(wp), intent(out) :: cmat(:, :) - integer :: iat, jat, izp, jzp, isp, jsp - real(wp) :: vec(3), rvdw, tmp, capi, capj + integer :: iat, jat, izp, jzp + real(wp) :: vec(3), rvdw, tmp, capi, capj, r1 ! Thread-private array for reduction real(wp), allocatable :: cmat_local(:, :) - + cmat(:, :) = 0.0_wp - + !$omp parallel default(none) & !$omp shared(cmat, mol, self) & - !$omp private(iat, izp, isp, jat, jzp, jsp) & - !$omp private(vec, rvdw, tmp, capi, capj, cmat_local) - allocate(cmat_local, source=cmat) - !$omp do schedule(runtime) + !$omp private(iat, izp, jat, jzp) & + !$omp private(vec, r1, rvdw, tmp, capi, capj, cmat_local) + allocate (cmat_local, source=cmat) + !$omp do schedule(runtime) do iat = 1, mol%nat izp = mol%id(iat) - isp = mol%num(izp) capi = self%cap(izp) do jat = 1, iat - 1 jzp = mol%id(jat) - jsp = mol%num(jzp) vec = mol%xyz(:, jat) - mol%xyz(:, iat) + r1 = norm2(vec) rvdw = self%rvdw(iat, jat) capj = self%cap(jzp) - call get_cpair(self%kbc, tmp, vec, rvdw, capi, capj) + + call get_cpair(self%kbc, tmp, r1, rvdw, capi, capj) + ! Off-diagonal elements cmat_local(jat, iat) = -tmp cmat_local(iat, jat) = -tmp @@ -875,72 +1046,138 @@ subroutine get_cmat_0d(self, mol, cmat) end do !$omp end do !$omp critical (get_cmat_0d_) - cmat(:, :) = cmat(:, :) + cmat_local(:, :) + cmat(:, :) = cmat + cmat_local !$omp end critical (get_cmat_0d_) - deallocate(cmat_local) + deallocate (cmat_local) !$omp end parallel cmat(mol%nat + 1, mol%nat + 1) = 1.0_wp end subroutine get_cmat_0d - subroutine get_cpair(kbc, cpair, vec, rvdw, capi, capj) - real(wp), intent(in) :: vec(3), capi, capj, rvdw, kbc - real(wp), intent(out) :: cpair - - real(wp) :: r2, arg - - r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 - ! Capacitance of bond between atom i and j - arg = -kbc*(sqrt(r2) - rvdw)/rvdw - cpair = sqrt(capi*capj)*0.5_wp*(1.0_wp + erf(arg)) - end subroutine get_cpair - - subroutine get_cdiag_3d(self, mol, wsc, cdiag) + 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 - real(wp), intent(out) :: cdiag(:, :) + real(wp), intent(out) :: cmat(:, :) - integer :: iat, jat, izp, jzp, isp, jsp, img - real(wp) :: vec(3), rvdw, capi, capj, tmp + integer :: iat, jat, izp, jzp, img + real(wp) :: vec(3), rvdw, tmp, capi, capj, wsw + real(wp), allocatable :: dtrans(:, :) ! Thread-private array for reduction - real(wp), allocatable :: cdiag_local(:, :) + real(wp), allocatable :: cmat_local(:, :) + + call get_dir_trans(mol%lattice, dtrans) - cdiag(:, :) = 0.0_wp + cmat(:, :) = 0.0_wp !$omp parallel default(none) & - !$omp shared(cdiag, mol, self, wsc) & - !$omp private(iat, izp, isp, jat, jzp, jsp, img) & - !$omp private(vec, rvdw, tmp, capi, capj, cdiag_local) - allocate(cdiag_local, source=cdiag) - !$omp do schedule(runtime) + !$omp shared(cmat, mol, self, wsc, dtrans) & + !$omp private(iat, izp, jat, jzp, img) & + !$omp private(vec, rvdw, tmp, capi, capj, wsw, cmat_local) + allocate (cmat_local, source=cmat) + !$omp do schedule(runtime) do iat = 1, mol%nat izp = mol%id(iat) - isp = mol%num(izp) - capi = self%cap(isp) + capi = self%cap(izp) do jat = 1, iat - 1 jzp = mol%id(jat) - jsp = mol%num(jzp) rvdw = self%rvdw(iat, jat) - capj = self%cap(jsp) + capj = self%cap(jzp) + 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)) - call get_cpair(self%kbc, tmp, vec, rvdw, capi, capj) - cdiag_local(iat, img) = cdiag_local(iat, img) + tmp - cdiag_local(jat, img) = cdiag_local(jat, img) + tmp + + call get_cpair_dir(self%kbc, vec, dtrans, rvdw, capi, capj, tmp) + + ! Off-diagonal elements + cmat_local(jat, iat) = cmat_local(jat, iat) - tmp*wsw + cmat_local(iat, jat) = cmat_local(iat, jat) - tmp*wsw + ! Diagonal elements + cmat_local(iat, iat) = cmat_local(iat, iat) + tmp*wsw + cmat_local(jat, jat) = cmat_local(jat, jat) + tmp*wsw end do end do + + ! diagonal capacitance (interaction with images) + rvdw = self%rvdw(iat, iat) + 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_cpair_dir(self%kbc, vec, dtrans, rvdw, capi, capi, tmp) + cmat_local(iat, iat) = cmat_local(iat, iat) + tmp*wsw + end do end do !$omp end do - !$omp critical (get_cdiag_3d_) - cdiag(:, :) = cdiag(:, :) + cdiag_local(:, :) - !$omp end critical (get_cdiag_3d_) - deallocate(cdiag_local) + !$omp critical (get_cmat_3d_) + cmat(:, :) = cmat + cmat_local + !$omp end critical (get_cmat_3d_) + deallocate (cmat_local) !$omp end parallel + ! + cmat(mol%nat + 1, mol%nat + 1) = 1.0_wp + + end subroutine get_cmat_3d + + subroutine get_cpair(kbc, cpair, r1, rvdw, capi, capj) + real(wp), intent(in) :: kbc + real(wp), intent(in) :: r1 + real(wp), intent(in) :: capi + real(wp), intent(in) :: capj + real(wp), intent(in) :: rvdw + real(wp), intent(out) :: cpair + + real(wp) :: arg + + ! Capacitance of bond between atom i and j + arg = -kbc*(r1 - rvdw)/rvdw + cpair = sqrt(capi*capj)*0.5_wp*(1.0_wp + erf(arg)) + end subroutine get_cpair + + subroutine get_cpair_dir(kbc, rij, trans, rvdw, capi, capj, cpair) + real(wp), intent(in) :: kbc + real(wp), intent(in) :: rij(3) + real(wp), intent(in) :: trans(:, :) + real(wp), intent(in) :: rvdw + real(wp), intent(in) :: capi + real(wp), intent(in) :: capj + real(wp), intent(out) :: cpair + + integer :: itr + real(wp) :: vec(3), r1, tmp + + cpair = 0.0_wp + do itr = 1, size(trans, 2) + vec(:) = rij + trans(:, itr) + r1 = norm2(vec) + if (r1 < eps) cycle + call get_cpair(kbc, tmp, r1, rvdw, capi, capj) + cpair = cpair + tmp + end do + end subroutine get_cpair_dir + + subroutine get_dcpair(kbc, vec, rvdw, capi, capj, dgpair, dspair) + real(wp), intent(in) :: kbc + real(wp), intent(in) :: vec(3) + real(wp), intent(in) :: rvdw + real(wp), intent(in) :: capi + real(wp), intent(in) :: capj + real(wp), intent(out) :: dgpair(3) + real(wp), intent(out) :: dspair(3, 3) - end subroutine get_cdiag_3d + real(wp) :: r1, arg, dtmp + + dgpair(:) = 0.0_wp + dspair(:, :) = 0.0_wp + + r1 = norm2(vec) + ! Capacitance of bond between atom i and j + arg = -(kbc*(r1 - rvdw)/rvdw)**2 + dtmp = sqrt(capi*capj)*kbc*exp(arg)/(sqrtpi*rvdw) + dgpair = dtmp*vec/r1 + dspair = spread(dgpair, 1, 3)*spread(vec, 2, 3) + end subroutine get_dcpair subroutine get_dcmat_0d(self, mol, dcdr, dcdL) class(eeqbc_model), intent(in) :: self @@ -962,8 +1199,8 @@ subroutine get_dcmat_0d(self, mol, dcdr, dcdL) !$omp private(iat, izp, jat, jzp, r2, vec, rvdw) & !$omp private(dG, dS, dtmp, arg, capi, capj) & !$omp private(dcdr_local, dcdL_local) - allocate(dcdr_local, source=dcdr) - allocate(dcdL_local, source=dcdL) + allocate (dcdr_local, source=dcdr) + allocate (dcdL_local, source=dcdL) !$omp do schedule(runtime) do iat = 1, mol%nat izp = mol%id(iat) @@ -971,16 +1208,10 @@ subroutine get_dcmat_0d(self, mol, dcdr, dcdL) do jat = 1, iat - 1 jzp = mol%id(jat) capj = self%cap(jzp) - vec = mol%xyz(:, jat) - mol%xyz(:, iat) - r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 rvdw = self%rvdw(iat, jat) + vec = mol%xyz(:, jat) - mol%xyz(:, iat) - ! Capacitance of bond between atom i and j - arg = -(self%kbc*(sqrt(r2) - rvdw)/rvdw)**2 - dtmp = sqrt(capi*capj)* & - & self%kbc*exp(arg)/(sqrtpi*rvdw) - dG = dtmp*vec/sqrt(r2) - dS = spread(dG, 1, 3)*spread(vec, 2, 3) + call get_dcpair(self%kbc, vec, rvdw, capi, capj, dG, dS) ! Negative off-diagonal elements dcdr_local(:, iat, jat) = -dG @@ -988,63 +1219,111 @@ subroutine get_dcmat_0d(self, mol, dcdr, dcdL) ! Positive diagonal elements dcdr_local(:, iat, iat) = +dG + dcdr_local(:, iat, iat) dcdr_local(:, jat, jat) = -dG + dcdr_local(:, jat, jat) - dcdL_local(:, :, jat) = +dS + dcdL_local(:, :, jat) - dcdL_local(:, :, iat) = +dS + dcdL_local(:, :, iat) + dcdL_local(:, :, iat) = -dS + dcdL_local(:, :, iat) + dcdL_local(:, :, jat) = -dS + dcdL_local(:, :, jat) end do end do !$omp end do !$omp critical (get_dcmat_0d_) - dcdr(:, :, :) = dcdr(:, :, :) + dcdr_local(:, :, :) - dcdL(:, :, :) = dcdL(:, :, :) + dcdL_local(:, :, :) + dcdr(:, :, :) = dcdr + dcdr_local + dcdL(:, :, :) = dcdL + dcdL_local !$omp end critical (get_dcmat_0d_) - deallocate(dcdL_local, dcdr_local) + deallocate (dcdL_local, dcdr_local) !$omp end parallel end subroutine get_dcmat_0d - subroutine write_2d_matrix(matrix, name, unit, step) - implicit none - real(wp), intent(in) :: matrix(:, :) - character(len=*), intent(in), optional :: name - integer, intent(in), optional :: unit - integer, intent(in), optional :: step - integer :: d1, d2 - integer :: i, j, k, l, istep, iunit + 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 + real(wp), intent(out) :: dcdr(:, :, :) + real(wp), intent(out) :: dcdL(:, :, :) - d1 = size(matrix, dim=1) - d2 = size(matrix, dim=2) + integer :: iat, jat, izp, jzp, img + real(wp) :: vec(3), r2, rvdw, dtmp, arg, dG(3), dS(3, 3), capi, capj, wsw + real(wp), allocatable :: dtrans(:, :) - if (present(unit)) then - iunit = unit - else - iunit = output_unit - end if + ! Thread-private arrays for reduction + real(wp), allocatable :: dcdr_local(:, :, :), dcdL_local(:, :, :) - if (present(step)) then - istep = step - else - istep = 6 - end if + call get_dir_trans(mol%lattice, dtrans) - if (present(name)) write (iunit, '(/,"matrix printed:",1x,a)') name + dcdr(:, :, :) = 0.0_wp + dcdL(:, :, :) = 0.0_wp - do i = 1, d2, istep - l = min(i + istep - 1, d2) - write (iunit, '(/,6x)', advance='no') - do k = i, l - write (iunit, '(6x,i7,3x)', advance='no') k - end do - write (iunit, '(a)') - do j = 1, d1 - write (iunit, '(i6)', advance='no') j - do k = i, l - write (iunit, '(1x,f15.8)', advance='no') matrix(j, k) + !$omp parallel default(none) & + !$omp shared(dcdr, dcdL, mol, self, dtrans, wsc) & + !$omp private(iat, izp, jat, jzp, r2, vec, rvdw) & + !$omp private(dG, dS, dtmp, arg, capi, capj, wsw) & + !$omp private(dcdr_local, dcdL_local) + allocate (dcdr_local, source=dcdr) + allocate (dcdL_local, source=dcdL) + !$omp do schedule(runtime) + do iat = 1, mol%nat + izp = mol%id(iat) + capi = self%cap(izp) + do jat = 1, iat - 1 + jzp = mol%id(jat) + capj = self%cap(jzp) + rvdw = self%rvdw(iat, jat) + wsw = 1/real(wsc%nimg(jat, iat), wp) + do img = 1, wsc%nimg(jat, iat) + vec = mol%xyz(:, jat) - mol%xyz(:, iat) + wsc%trans(:, wsc%tridx(img, jat, iat)) + + call get_dcpair_dir(self%kbc, vec, dtrans, rvdw, capi, capj, dG, dS) + + ! Negative off-diagonal elements + dcdr_local(:, iat, jat) = -dG*wsw + dcdr_local(:, iat, jat) + dcdr_local(:, jat, iat) = +dG*wsw + dcdr_local(:, jat, iat) + ! Positive diagonal elements + dcdr_local(:, iat, iat) = +dG*wsw + dcdr_local(:, iat, iat) + dcdr_local(:, jat, jat) = -dG*wsw + dcdr_local(:, jat, jat) + dcdL_local(:, :, jat) = -dS*wsw + dcdL_local(:, :, jat) + dcdL_local(:, :, iat) = -dS*wsw + dcdL_local(:, :, iat) end do - write (iunit, '(a)') + end do + + rvdw = self%rvdw(iat, iat) + wsw = 1/real(wsc%nimg(iat, iat), wp) + do img = 1, wsc%nimg(iat, iat) + vec = wsc%trans(:, wsc%tridx(img, iat, iat)) + + call get_dcpair_dir(self%kbc, vec, dtrans, rvdw, capi, capi, dG, dS) + + ! Positive diagonal elements + dcdL_local(:, :, iat) = -dS*wsw + dcdL_local(:, :, iat) end do end do + !$omp end do + !$omp critical (get_dcmat_3d_) + dcdr(:, :, :) = dcdr + dcdr_local + dcdL(:, :, :) = dcdL + dcdL_local + !$omp end critical (get_dcmat_3d_) + deallocate (dcdL_local, dcdr_local) + !$omp end parallel + + end subroutine get_dcmat_3d + + subroutine get_dcpair_dir(kbc, rij, trans, rvdw, capi, capj, dgpair, dspair) + real(wp), intent(in) :: rij(3), capi, capj, rvdw, kbc, trans(:, :) + real(wp), intent(out) :: dgpair(3) + real(wp), intent(out) :: dspair(3, 3) - end subroutine write_2d_matrix + integer :: itr + real(wp) :: r1, arg, dtmp, dgtmp(3), dstmp(3, 3), vec(3) + + dgpair(:) = 0.0_wp + dspair(:, :) = 0.0_wp + do itr = 1, size(trans, 2) + vec(:) = rij + trans(:, itr) + r1 = norm2(vec) + if (r1 < eps) cycle + call get_dcpair(kbc, vec, rvdw, capi, capj, dgtmp, dstmp) + dgpair(:) = dgpair + dgtmp + dspair(:, :) = dspair + dstmp + end do + end subroutine get_dcpair_dir ! NOTE: the following is basically identical to tblite versions of this pattern diff --git a/src/multicharge/model/type.F90 b/src/multicharge/model/type.F90 index 1df6173d..77623719 100644 --- a/src/multicharge/model/type.F90 +++ b/src/multicharge/model/type.F90 @@ -23,8 +23,6 @@ !> General charge model module multicharge_model_type - use iso_fortran_env, only: output_unit - use mctc_env, only: error_type, fatal_error, wp, ik => IK use mctc_io, only: structure_type use mctc_io_constants, only: pi @@ -122,61 +120,9 @@ subroutine get_xvec_derivs(self, mol, cache, dxdr, dxdL) class(mchrg_model_type), 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(:, :, :) end subroutine get_xvec_derivs - - !subroutine get_amat_0d(self, mol, amat, cn, qloc, cmat) - ! import :: mchrg_model_type, cache_container, structure_type, wp - ! class(mchrg_model_type), intent(in) :: self - ! type(structure_type), intent(in) :: mol - ! real(wp), intent(out) :: amat(:, :) - ! real(wp), intent(in), optional :: cn(:) - ! real(wp), intent(in), optional :: qloc(:) - ! real(wp), intent(in), optional :: cmat(:, :) - !end subroutine get_amat_0d - - !subroutine get_amat_3d(self, mol, cache, wsc, alpha, amat) - ! import :: mchrg_model_type, cache_container, structure_type, & - ! & wignerseitz_cell_type, wp - ! class(mchrg_model_type), intent(in) :: self - ! type(structure_type), intent(in) :: mol - ! type(cache_container), intent(inout) :: cache - ! type(wignerseitz_cell_type), intent(in) :: wsc - ! real(wp), intent(in) :: alpha - ! real(wp), intent(out) :: amat(:, :) - !end subroutine get_amat_3d - - !subroutine get_damat_0d(self, mol, cn, qloc, qvec, dcndr, dcndL, & - ! & dqlocdr, dqlocdL, dadr, dadL, atrace) - ! import :: mchrg_model_type, structure_type, wp - ! class(mchrg_model_type), intent(in) :: self - ! type(structure_type), intent(in) :: mol - ! real(wp), intent(in) :: qvec(:) - ! real(wp), intent(out) :: dadr(:, :, :) - ! real(wp), intent(out) :: dadL(:, :, :) - ! real(wp), intent(out) :: atrace(:, :) - ! real(wp), intent(in), optional :: cn(:) - ! real(wp), intent(in), optional :: qloc(:) - ! real(wp), intent(in), optional :: dcndr(:, :, :) - ! real(wp), intent(in), optional :: dcndL(:, :, :) - ! real(wp), intent(in), optional :: dqlocdr(:, :, :) - ! real(wp), intent(in), optional :: dqlocdL(:, :, :) - !end subroutine get_damat_0d - - !subroutine get_damat_3d(self, mol, wsc, alpha, qvec, dadr, dadL, atrace) - ! import :: mchrg_model_type, structure_type, & - ! & wignerseitz_cell_type, wp - ! class(mchrg_model_type), intent(in) :: self - ! type(structure_type), intent(in) :: mol - ! type(wignerseitz_cell_type), intent(in) :: wsc - ! real(wp), intent(in) :: alpha - ! real(wp), intent(in) :: qvec(:) - ! real(wp), intent(out) :: dadr(:, :, :) - ! real(wp), intent(out) :: dadL(:, :, :) - ! real(wp), intent(out) :: atrace(:, :) - !end subroutine get_damat_3d - end interface real(wp), parameter :: twopi = 2*pi @@ -315,7 +261,7 @@ subroutine solve(self, mol, error, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL, & if (present(energy)) then ! Extract only the Coulomb matrix without the constraints - allocate(jmat(mol%nat, mol%nat)) + allocate (jmat(mol%nat, mol%nat)) jmat = amat(:mol%nat, :mol%nat) call symv(jmat, vrhs(:mol%nat), xvec(:mol%nat), & & alpha=0.5_wp, beta=-1.0_wp, uplo='l') @@ -330,7 +276,7 @@ subroutine solve(self, mol, error, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL, & call self%get_coulomb_derivs(mol, cache, vrhs, dadr, dadL, atrace) do iat = 1, mol%nat dadr(:, iat, iat) = atrace(:, iat) + dadr(:, iat, iat) - end do + end do end if if (grad) then @@ -382,48 +328,4 @@ subroutine local_charge(self, mol, trans, qloc, dqlocdr, dqlocdL) end subroutine local_charge - subroutine write_2d_matrix(matrix, name, unit, step) - implicit none - real(wp), intent(in) :: matrix(:, :) - character(len=*), intent(in), optional :: name - integer, intent(in), optional :: unit - integer, intent(in), optional :: step - integer :: d1, d2 - integer :: i, j, k, l, istep, iunit - - d1 = size(matrix, dim=1) - d2 = size(matrix, dim=2) - - if (present(unit)) then - iunit = unit - else - iunit = output_unit - end if - - if (present(step)) then - istep = step - else - istep = 6 - end if - - if (present(name)) write (iunit, '(/,"matrix printed:",1x,a)') name - - do i = 1, d2, istep - l = min(i + istep - 1, d2) - write (iunit, '(/,6x)', advance='no') - do k = i, l - write (iunit, '(6x,i7,3x)', advance='no') k - end do - write (iunit, '(a)') - do j = 1, d1 - write (iunit, '(i6)', advance='no') j - do k = i, l - write (iunit, '(1x,f15.8)', advance='no') matrix(j, k) - end do - write (iunit, '(a)') - end do - end do - - end subroutine write_2d_matrix - end module multicharge_model_type diff --git a/src/multicharge/param.f90 b/src/multicharge/param.f90 index 80300c7f..3cc0da86 100644 --- a/src/multicharge/param.f90 +++ b/src/multicharge/param.f90 @@ -112,10 +112,10 @@ subroutine new_eeqbc2025_model(mol, model, error) &.or. mol%num == 97 .or. mol%num == 103) en = en/3.98_wp rvdw = get_vdw_rad(spread(mol%num(mol%id), 2, mol%nat), & - & spread(mol%num(mol%id), 1, mol%nat)) * autoaa + & spread(mol%num(mol%id), 1, mol%nat))*autoaa allocate (eeqbc) - call new_eeqbc_model(eeqbc, mol=mol, error=error, chi=chi, & + call new_eeqbc_model(eeqbc, mol=mol, error=error, chi=chi, & & rad=rad, eta=eta, kcnchi=kcnchi, kqchi=kqchi, kqeta=kqeta, & & kcnrad=kcnrad, cap=cap, avg_cn=avg_cn, kbc=kbc, & & cutoff=cutoff, cn_exp=cn_exp, rcov=rcov, en=en, & diff --git a/src/multicharge/wignerseitz.f90 b/src/multicharge/wignerseitz.f90 index 797b7c19..ac7ba5f7 100644 --- a/src/multicharge/wignerseitz.f90 +++ b/src/multicharge/wignerseitz.f90 @@ -14,9 +14,9 @@ ! limitations under the License. module multicharge_wignerseitz - use mctc_env, only : wp - use mctc_io, only : structure_type - use mctc_cutoff, only : get_lattice_points + use mctc_env, only: wp + use mctc_io, only: structure_type + use mctc_cutoff, only: get_lattice_points implicit none private @@ -29,96 +29,91 @@ module multicharge_wignerseitz real(wp), allocatable :: trans(:, :) end type wignerseitz_cell_type - !> Small cutoff threshold to create only closest cells real(wp), parameter :: thr = sqrt(epsilon(0.0_wp)) !> Tolerance to consider equivalent images real(wp), parameter :: tol = 0.01_wp - contains + subroutine new_wignerseitz_cell(self, mol) -subroutine new_wignerseitz_cell(self, mol) - - !> Wigner-Seitz cell instance - type(wignerseitz_cell_type), intent(out) :: self + !> Wigner-Seitz cell instance + type(wignerseitz_cell_type), intent(out) :: self - !> Molecular structure data - type(structure_type), intent(in) :: mol + !> Molecular structure data + type(structure_type), intent(in) :: mol - integer :: iat, jat, ntr, nimg - integer, allocatable :: tridx(:) - real(wp) :: vec(3) - real(wp), allocatable :: trans(:, :) - - call get_lattice_points(mol%periodic, mol%lattice, thr, trans) - ntr = size(trans, 2) - allocate(self%nimg(mol%nat, mol%nat), self%tridx(ntr, mol%nat, mol%nat), & - & tridx(ntr)) + integer :: iat, jat, ntr, nimg + integer, allocatable :: tridx(:) + real(wp) :: vec(3) + real(wp), allocatable :: trans(:, :) - self%nimg_max = 0 - !$omp parallel do default(none) schedule(runtime) collapse(2) & - !$omp shared(mol, trans, self) private(iat, jat, vec, nimg, tridx) - do iat = 1, mol%nat - do jat = 1, mol%nat - vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) - call get_pairs(nimg, trans, vec, tridx) - self%nimg(jat, iat) = nimg - self%tridx(:, jat, iat) = tridx - self%nimg_max = max(nimg, self%nimg_max) + call get_lattice_points(mol%periodic, mol%lattice, thr, trans) + ntr = size(trans, 2) + allocate (self%nimg(mol%nat, mol%nat), self%tridx(ntr, mol%nat, mol%nat), & + & tridx(ntr)) + + self%nimg_max = 0 + !$omp parallel do default(none) schedule(runtime) collapse(2) & + !$omp shared(mol, trans, self) private(iat, jat, vec, nimg, tridx) + do iat = 1, mol%nat + do jat = 1, mol%nat + vec(:) = mol%xyz(:, iat) - mol%xyz(:, jat) + call get_pairs(nimg, trans, vec, tridx) + self%nimg(jat, iat) = nimg + self%tridx(:, jat, iat) = tridx + self%nimg_max = max(nimg, self%nimg_max) + end do end do - end do - - call move_alloc(trans, self%trans) - -end subroutine new_wignerseitz_cell + call move_alloc(trans, self%trans) -subroutine get_pairs(iws, trans, rij, list) - integer, intent(out) :: iws - real(wp), intent(in) :: rij(3) - real(wp), intent(in) :: trans(:, :) - integer, intent(out) :: list(:) + end subroutine new_wignerseitz_cell - logical :: mask(size(list)) - real(wp) :: dist(size(list)), vec(3), r2 - integer :: itr, img, pos + subroutine get_pairs(iws, trans, rij, list) + integer, intent(out) :: iws + real(wp), intent(in) :: rij(3) + real(wp), intent(in) :: trans(:, :) + integer, intent(out) :: list(:) - iws = 0 - img = 0 - list(:) = 0 - mask(:) = .true. + logical :: mask(size(list)) + real(wp) :: dist(size(list)), vec(3), r2 + integer :: itr, img, pos - do itr = 1, size(trans, 2) - vec(:) = rij - trans(:, itr) - r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 - if (r2 < thr) cycle - img = img + 1 - dist(img) = r2 - end do + iws = 0 + img = 0 + list(:) = 0 + mask(:) = .true. - if (img == 0) return - - pos = minloc(dist(:img), dim=1) + do itr = 1, size(trans, 2) + vec(:) = rij - trans(:, itr) + r2 = vec(1)**2 + vec(2)**2 + vec(3)**2 + if (r2 < thr) cycle + img = img + 1 + dist(img) = r2 + end do - r2 = dist(pos) - mask(pos) = .false. + if (img == 0) return - iws = 1 - list(iws) = pos - if (img <= iws) return + pos = minloc(dist(:img), dim=1) - do - pos = minloc(dist(:img), dim=1, mask=mask(:img)) - if (abs(dist(pos) - r2) > tol) exit + r2 = dist(pos) mask(pos) = .false. - iws = iws + 1 - list(iws) = pos - end do -end subroutine get_pairs + iws = 1 + list(iws) = pos + if (img <= iws) return + + do + pos = minloc(dist(:img), dim=1, mask=mask(:img)) + if (abs(dist(pos) - r2) > tol) exit + mask(pos) = .false. + iws = iws + 1 + list(iws) = pos + end do + end subroutine get_pairs end module multicharge_wignerseitz diff --git a/test/unit/test_model.f90 b/test/unit/test_model.f90 index d7128601..56b652d6 100644 --- a/test/unit/test_model.f90 +++ b/test/unit/test_model.f90 @@ -22,6 +22,7 @@ module test_model use mctc_io_structure, only: structure_type, new use mstore, only: get_structure use multicharge_model, only: mchrg_model_type + use multicharge_model_eeqbc, only: eeqbc_model use multicharge_param, only: new_eeq2019_model, new_eeqbc2025_model use multicharge_model_cache, only: cache_container use multicharge_blas, only: gemv @@ -44,8 +45,9 @@ subroutine collect_model(testsuite) testsuite = [ & & new_unittest("eeq-dadr-mb01", test_eeq_dadr_mb01), & - ! !& new_unittest("eeq-dadL-mb01", test_eeq_dadL_mb01), & + & new_unittest("eeq-dadL-mb01", test_eeq_dadL_mb01), & & new_unittest("eeq-dbdr-mb01", test_eeq_dbdr_mb01), & + & new_unittest("eeq-dbdL-mb01", test_eeq_dbdL_mb01), & & new_unittest("eeq-charges-mb01", test_eeq_q_mb01), & & new_unittest("eeq-charges-mb02", test_eeq_q_mb02), & & new_unittest("eeq-charges-actinides", test_eeq_q_actinides), & @@ -64,10 +66,12 @@ subroutine collect_model(testsuite) & new_unittest("eeq-dbdr-znooh", test_eeq_dbdr_znooh), & & new_unittest("gradient-znooh", test_g_znooh), & & new_unittest("dqdr-znooh", test_dqdr_znooh), & - !& new_unittest("eeqbc-dadr-mb01", test_eeqbc_dadr_mb01), & - !& new_unittest("eeqbc-dadL-mb01", test_eeqbc_dadL_mb01), & + & new_unittest("eeqbc-dadr-mb01", test_eeqbc_dadr_mb01), & ! fails randomly due to numerical noise? + & new_unittest("eeqbc-dadL-mb01", test_eeqbc_dadL_mb01), & & new_unittest("eeqbc-dbdr-mb01", test_eeqbc_dbdr_mb01), & + & new_unittest("eeqbc-dbdL-mb01", test_eeqbc_dbdL_mb01), & & new_unittest("eeqbc-dadr-mb05", test_eeqbc_dadr_mb05), & + & new_unittest("eeqbc-dadL-mb05", test_eeqbc_dadL_mb05), & & new_unittest("eeqbc-dbdr-mb05", test_eeqbc_dbdr_mb05), & & new_unittest("eeqbc-charges-mb01", test_eeqbc_q_mb01), & & new_unittest("eeqbc-charges-mb02", test_eeqbc_q_mb02), & @@ -76,12 +80,12 @@ subroutine collect_model(testsuite) & new_unittest("eeqbc-energy-mb04", test_eeqbc_e_mb04), & & new_unittest("eeqbc-gradient-mb05", test_eeqbc_g_mb05), & & new_unittest("eeqbc-gradient-mb06", test_eeqbc_g_mb06), & - !& new_unittest("eeqbc-sigma-mb07", test_eeqbc_s_mb07), & - !& new_unittest("eeqbc-sigma-mb08", test_eeqbc_s_mb08), & + & new_unittest("eeqbc-sigma-mb07", test_eeqbc_s_mb07), & + & new_unittest("eeqbc-sigma-mb08", test_eeqbc_s_mb08), & & new_unittest("eeqbc-dqdr-mb09", test_eeqbc_dqdr_mb09), & - & new_unittest("eeqbc-dqdr-mb10", test_eeqbc_dqdr_mb10) & - !& new_unittest("eeqbc-dqdL-mb11", test_eeqbc_dqdL_mb11), & - !& new_unittest("eeqbc-dqdL-mb12", test_eeqbc_dqdL_mb12) & + & new_unittest("eeqbc-dqdr-mb10", test_eeqbc_dqdr_mb10), & + & new_unittest("eeqbc-dqdL-mb11", test_eeqbc_dqdL_mb11), & + & new_unittest("eeqbc-dqdL-mb12", test_eeqbc_dqdL_mb12) & & ] end subroutine collect_model @@ -98,56 +102,89 @@ subroutine test_dadr(error, mol, model) type(error_type), allocatable, intent(out) :: error integer :: iat, ic, jat, kat + real(wp) :: thr2_local real(wp), parameter :: trans(3, 1) = 0.0_wp real(wp), parameter :: step = 1.0e-6_wp real(wp), allocatable :: cn(:) real(wp), allocatable :: qloc(:) real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :), dqlocdr(:, :, :), dqlocdL(:, :, :) real(wp), allocatable :: dadr(:, :, :), dadL(:, :, :), atrace(:, :) - real(wp), allocatable :: qvec(:), numgrad(:, :, :), amatr(:, :), amatl(:, :), numtrace(:, :) + real(wp), allocatable :: qvec(:), numgrad(:, :, :), amatr1(:, :), amatr2(:, :), amatl1(:, :), amatl2(:, :), numtrace(:, :) type(cache_container), allocatable :: cache allocate (cache) - allocate (cn(mol%nat), qloc(mol%nat), amatr(mol%nat + 1, mol%nat + 1), amatl(mol%nat + 1, mol%nat + 1), & + allocate (cn(mol%nat), qloc(mol%nat), amatr1(mol%nat + 1, mol%nat + 1), amatl1(mol%nat + 1, mol%nat + 1), & + & amatr2(mol%nat + 1, mol%nat + 1), amatl2(mol%nat + 1, mol%nat + 1), & & dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), dqlocdr(3, mol%nat, mol%nat), & & dqlocdL(3, 3, mol%nat), dadr(3, mol%nat, mol%nat + 1), dadL(3, 3, mol%nat + 1), & & atrace(3, mol%nat), numtrace(3, mol%nat), numgrad(3, mol%nat, mol%nat + 1), qvec(mol%nat)) + ! Set tolerance higher if testing eeqbc model + select type (model) + type is (eeqbc_model) + thr2_local = 3.0_wp*thr2 + class default + thr2_local = thr2 + end select + ! Obtain the vector of charges call model%ncoord%get_coordination_number(mol, trans, cn) call model%local_charge(mol, trans, qloc) call model%solve(mol, error, cn, qloc, qvec=qvec) - if(allocated(error)) return + if (allocated(error)) return numgrad = 0.0_wp lp: do iat = 1, mol%nat do ic = 1, 3 - ! Right-hand side - amatr(:, :) = 0.0_wp + amatr1(:, :) = 0.0_wp + amatr2(:, :) = 0.0_wp + amatl1(:, :) = 0.0_wp + amatl2(:, :) = 0.0_wp + + ! First right-hand side (x+h) mol%xyz(ic, iat) = mol%xyz(ic, iat) + step call model%ncoord%get_coordination_number(mol, trans, cn) call model%local_charge(mol, trans, qloc) call model%update(mol, cache, cn, qloc) - call model%get_coulomb_matrix(mol, cache, amatr) + call model%get_coulomb_matrix(mol, cache, amatr1) - ! Left-hand side - amatl(:, :) = 0.0_wp + ! Second right-hand side (x+2h) + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_coulomb_matrix(mol, cache, amatr2) + + ! Return to original position before calculating left sides mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + + ! First left-hand side (x-h) + mol%xyz(ic, iat) = mol%xyz(ic, iat) - step call model%ncoord%get_coordination_number(mol, trans, cn) call model%local_charge(mol, trans, qloc) call model%update(mol, cache, cn, qloc) - call model%get_coulomb_matrix(mol, cache, amatl) + call model%get_coulomb_matrix(mol, cache, amatl1) + + ! Second left-hand side (x-2h) + mol%xyz(ic, iat) = mol%xyz(ic, iat) - step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_coulomb_matrix(mol, cache, amatl2) + + ! Return to original position + mol%xyz(ic, iat) = mol%xyz(ic, iat) + 2*step - mol%xyz(ic, iat) = mol%xyz(ic, iat) + step - do kat = 1, mol%nat do jat = 1, mol%nat - ! Numerical gradient of the A matrix - numgrad(ic, iat, kat) = 0.5_wp*qvec(jat)*(amatr(kat, jat) - amatl(kat, jat))/step & - & + numgrad(ic, iat, kat) - end do - end do + ! Numerical gradient using 4-step central difference formula + ! f'(x) ≈ [f(x-2h) - 8f(x-h) + 8f(x+h) - f(x+2h)]/(12h) + numgrad(ic, iat, kat) = numgrad(ic, iat, kat) + & + & qvec(jat)*(amatl2(kat, jat) - 8.0_wp*amatl1(kat, jat) + & + & 8.0_wp*amatr1(kat, jat) - amatr2(kat, jat))/(12.0_wp*step) + end do + end do end do end do lp @@ -160,37 +197,18 @@ subroutine test_dadr(error, mol, model) ! Add trace of the A matrix do iat = 1, mol%nat dadr(:, iat, iat) = atrace(:, iat) + dadr(:, iat, iat) - end do + end do - if (any(abs(dadr(:, :, :) - numgrad(:, :, :)) > thr2)) then + if (any(abs(dadr(:, :, :) - numgrad(:, :, :)) > thr2_local)) then call test_failed(error, "Derivative of the A matrix does not match") print'(a)', "dadr:" - print'(3es21.14)', dadr + print'(3es21.12)', dadr print'(a)', "numgrad:" - print'(3es21.14)', numgrad + print'(3es21.12)', numgrad print'(a)', "diff:" - print'(3es21.14)', dadr - numgrad + print'(3es21.12)', dadr - numgrad end if - ! numtrace(:, :) = 0.0_wp - ! do iat = 1, mol%nat - ! do jat = 1, iat - 1 - ! ! Numerical trace of the a matrix - ! numtrace(:, iat) = - numgrad(:, jat, iat) + numtrace(:, iat) - ! numtrace(:, jat) = - numgrad(:, iat, jat) + numtrace(:, jat) - ! end do - ! end do - - ! if (any(abs(atrace(:, :) - numtrace(:, :)) > thr2)) then - ! call test_failed(error, "Derivative of the A matrix trace does not match") - ! print'(a)', "atrace:" - ! print'(3es21.14)', atrace - ! print'(a)', "numtrace:" - ! print'(3es21.14)', numtrace - ! print'(a)', "diff:" - ! print'(3es21.14)', atrace - numtrace - ! end if - end subroutine test_dadr subroutine test_dadL(error, mol, model) @@ -211,7 +229,7 @@ subroutine test_dadL(error, mol, model) real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :) real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) real(wp), allocatable :: dadr(:, :, :), dadL(:, :, :), atrace(:, :) - real(wp), allocatable :: lattr(:, :), xyz(:, :) + real(wp), allocatable :: xyz(:, :) real(wp), allocatable :: qvec(:), numsigma(:, :, :), amatr(:, :), amatl(:, :) real(wp) :: eps(3, 3) type(cache_container), allocatable :: cache @@ -226,19 +244,17 @@ subroutine test_dadL(error, mol, model) call model%ncoord%get_coordination_number(mol, trans, cn) call model%local_charge(mol, trans, qloc) call model%solve(mol, error, cn, qloc, qvec=qvec) - if(allocated(error)) return + if (allocated(error)) return - qvec = 1.0_wp + numsigma = 0.0_wp eps(:, :) = unity xyz(:, :) = mol%xyz - lattr = trans lp: do ic = 1, 3 do jc = 1, 3 amatr(:, :) = 0.0_wp eps(jc, ic) = eps(jc, ic) + step mol%xyz(:, :) = matmul(eps, xyz) - lattr(:, :) = matmul(eps, trans) call model%ncoord%get_coordination_number(mol, trans, cn) call model%local_charge(mol, trans, qloc) call model%update(mol, cache, cn, qloc) @@ -248,7 +264,6 @@ subroutine test_dadL(error, mol, model) amatl(:, :) = 0.0_wp eps(jc, ic) = eps(jc, ic) - 2*step mol%xyz(:, :) = matmul(eps, xyz) - lattr(:, :) = matmul(eps, trans) call model%ncoord%get_coordination_number(mol, trans, cn) call model%local_charge(mol, trans, qloc) call model%update(mol, cache, cn, qloc) @@ -257,7 +272,6 @@ subroutine test_dadL(error, mol, model) eps(jc, ic) = eps(jc, ic) + step mol%xyz(:, :) = xyz - lattr(:, :) = trans do iat = 1, mol%nat ! Numerical sigma of the a matrix numsigma(jc, ic, :) = 0.5_wp*qvec(iat)*(amatr(iat, :) - amatl(iat, :))/step + numsigma(jc, ic, :) @@ -270,34 +284,17 @@ subroutine test_dadL(error, mol, model) call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) call model%update(mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) - ! dcndr(:, :, :) = 0.0_wp - ! dcndL(:, :, :) = 0.0_wp - ! dqlocdr(:, :, :) = 0.0_wp - ! dqlocdL(:, :, :) = 0.0_wp - call model%get_coulomb_derivs(mol, cache, qvec, dadr, dadL, atrace) if (allocated(error)) return - ! do iat = 1, mol%nat - ! write(*,*) "iat", iat - ! call write_2d_matrix(dadL(:, :, iat), "dadL", unit=output_unit) - ! call write_2d_matrix(numsigma(:, :, iat), "numsigma", unit=output_unit) - ! call write_2d_matrix(dadL(:, :, iat) - numsigma(:, :, iat), "diff", unit=output_unit) - ! end do - - ! do ic = 1, 3 - ! do jc = 1, 3 - ! write(*,*) "ic, jc", ic, jc - ! write(*,*) dadL(ic, jc, :) - numsigma(ic, jc, :) - ! end do - ! end do - if (any(abs(dadL(:, :, :) - numsigma(:, :, :)) > thr2)) then call test_failed(error, "Derivative of the A matrix does not match") - !print'(a)', "dadr:" - !print'(3es21.14)', dadr - !print'(a)', "diff:" - !print'(3es21.14)', dadr - numsigma + print'(a)', "dadL:" + print'(3es21.12)', dadL + print'(a)', "numsigma:" + print'(3es21.12)', numsigma + print'(a)', "diff:" + print'(3es21.12)', dadL - numsigma end if end subroutine test_dadL @@ -370,49 +367,85 @@ subroutine test_dbdr(error, mol, model) end subroutine test_dbdr - subroutine write_2d_matrix(matrix, name, unit, step) - implicit none - real(wp), intent(in) :: matrix(:, :) - character(len=*), intent(in), optional :: name - integer, intent(in), optional :: unit - integer, intent(in), optional :: step - integer :: d1, d2 - integer :: i, j, k, l, istep, iunit - - d1 = size(matrix, dim=1) - d2 = size(matrix, dim=2) - - if (present(unit)) then - iunit = unit - else - iunit = output_unit - end if + subroutine test_dbdL(error, mol, model) - if (present(step)) then - istep = step - else - istep = 6 - end if + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model - if (present(name)) write (iunit, '(/,"matrix printed:",1x,a)') name + !> Error handling + type(error_type), allocatable, intent(out) :: error - do i = 1, d2, istep - l = min(i + istep - 1, d2) - write (iunit, '(/,6x)', advance='no') - do k = i, l - write (iunit, '(6x,i7,3x)', advance='no') k - end do - write (iunit, '(a)') - do j = 1, d1 - write (iunit, '(i6)', advance='no') j - do k = i, l - write (iunit, '(1x,f15.8)', advance='no') matrix(j, k) + integer :: iat, ic, jc + real(wp), parameter :: trans(3, 1) = 0.0_wp + real(wp), parameter :: step = 1.0e-6_wp, unity(3, 3) = reshape(& + & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: dbdr(:, :, :), dbdL(:, :, :) + real(wp), allocatable :: numsigma(:, :, :), xvecr(:), xvecl(:) + real(wp), allocatable :: xyz(:, :) + real(wp) :: eps(3, 3) + type(cache_container), allocatable :: cache + allocate (cache) + + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & xvecr(mol%nat + 1), xvecl(mol%nat + 1), numsigma(3, 3, mol%nat + 1), & + & dbdr(3, mol%nat, mol%nat + 1), dbdL(3, 3, mol%nat + 1), xyz(3, mol%nat)) + + numsigma = 0.0_wp + + eps(:, :) = unity + xyz(:, :) = mol%xyz + lp: do ic = 1, 3 + do jc = 1, 3 + ! Right-hand side + xvecr(:) = 0.0_wp + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_xvec(mol, cache, xvecr) + + ! Left-hand side + xvecl(:) = 0.0_wp + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_xvec(mol, cache, xvecl) + + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + do iat = 1, mol%nat + numsigma(jc, ic, iat) = 0.5_wp*(xvecr(iat) - xvecl(iat))/step end do - write (iunit, '(a)') end do - end do + end do lp - end subroutine write_2d_matrix + ! Analytical gradient + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + call model%update(mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) + call model%get_xvec(mol, cache, xvecl) ! need to call this for xtmp in cache (eeqbc) + call model%get_xvec_derivs(mol, cache, dbdr, dbdL) + + if (any(abs(dbdL(:, :, :) - numsigma(:, :, :)) > thr2)) then + call test_failed(error, "Derivative of the b vector does not match") + print'(a)', "dbdL:" + print'(3es21.14)', dbdL + print'(a)', "numsigma:" + print'(3es21.14)', numsigma + print'(a)', "diff:" + print'(3es21.14)', dbdL - numsigma + end if + + end subroutine test_dbdL subroutine gen_test(error, mol, model, qref, eref) @@ -816,6 +849,21 @@ subroutine test_eeq_dbdr_mb01(error) end subroutine test_eeq_dbdr_mb01 + subroutine test_eeq_dbdL_mb01(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "MB16-43", "01") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_dbdL(error, mol, model) + + end subroutine test_eeq_dbdL_mb01 + subroutine test_eeq_q_mb01(error) !> Error handling @@ -840,7 +888,7 @@ subroutine test_eeq_q_mb01(error) if (allocated(error)) return ! Check wrapper functions - allocate (qvec(mol%nat), source = 0.0_wp) + allocate (qvec(mol%nat), source=0.0_wp) call get_charges(model, mol, error, qvec) if (allocated(error)) return @@ -853,7 +901,7 @@ subroutine test_eeq_q_mb01(error) end if if (allocated(error)) return - qvec = 0.0_wp + qvec = 0.0_wp call get_eeq_charges(mol, error, qvec) if (allocated(error)) return @@ -1281,6 +1329,21 @@ subroutine test_eeqbc_dbdr_mb01(error) end subroutine test_eeqbc_dbdr_mb01 + subroutine test_eeqbc_dbdL_mb01(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "MB16-43", "01") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_dbdL(error, mol, model) + + end subroutine test_eeqbc_dbdL_mb01 + subroutine test_eeqbc_dadr_mb05(error) !> Error handling @@ -1296,6 +1359,21 @@ subroutine test_eeqbc_dadr_mb05(error) end subroutine test_eeqbc_dadr_mb05 + subroutine test_eeqbc_dadL_mb05(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "MB16-43", "05") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_dadL(error, mol, model) + + end subroutine test_eeqbc_dadL_mb05 + subroutine test_eeqbc_dbdr_mb05(error) !> Error handling @@ -1334,7 +1412,7 @@ subroutine test_eeqbc_q_mb01(error) call gen_test(error, mol, model, qref=ref) ! Check wrapper functions - allocate (qvec(mol%nat), source = 0.0_wp) + allocate (qvec(mol%nat), source=0.0_wp) call get_charges(model, mol, error, qvec) if (allocated(error)) return @@ -1347,7 +1425,7 @@ subroutine test_eeqbc_q_mb01(error) end if if (allocated(error)) return - qvec = 0.0_wp + qvec = 0.0_wp call get_eeqbc_charges(mol, error, qvec) if (allocated(error)) return diff --git a/test/unit/test_pbc.f90 b/test/unit/test_pbc.f90 index 6f62323f..dc92828e 100644 --- a/test/unit/test_pbc.f90 +++ b/test/unit/test_pbc.f90 @@ -14,14 +14,16 @@ ! limitations under the License. module test_pbc - use mctc_env, only : wp - use mctc_env_testing, only : new_unittest, unittest_type, error_type, & + use mctc_env, only: wp + use mctc_env_testing, only: new_unittest, unittest_type, error_type, & & test_failed - use mctc_io_structure, only : structure_type - use mctc_cutoff, only : get_lattice_points - use mstore, only : get_structure - use multicharge_model, only : mchrg_model_type - use multicharge_param, only : new_eeq2019_model + use mctc_io_structure, only: structure_type + use mctc_cutoff, only: get_lattice_points + use mstore, only: get_structure + use multicharge_model, only: mchrg_model_type + use multicharge_model_eeqbc, only: eeqbc_model + use multicharge_param, only: new_eeq2019_model, new_eeqbc2025_model + use multicharge_model_cache, only: cache_container implicit none private @@ -30,498 +32,1069 @@ module test_pbc real(wp), parameter :: thr = 1000*epsilon(1.0_wp) real(wp), parameter :: thr2 = sqrt(epsilon(1.0_wp)) - contains - !> Collect all exported unit tests -subroutine collect_pbc(testsuite) + subroutine collect_pbc(testsuite) - !> Collection of tests - type(unittest_type), allocatable, intent(out) :: testsuite(:) + !> Collection of tests + type(unittest_type), allocatable, intent(out) :: testsuite(:) - testsuite = [ & - & new_unittest("charges-cyanamide", test_q_cyanamide), & - & new_unittest("energy-formamide", test_e_formamide), & - & new_unittest("gradient-co2", test_g_co2), & - & new_unittest("sigma-ice", test_s_ice), & - & new_unittest("dqdr-urea", test_dqdr_urea), & - & new_unittest("dqdL-oxacb", test_dqdL_oxacb) & - & ] + testsuite = [ & + & new_unittest("eeq-charges-cyanamide", test_eeq_q_cyanamide), & + & new_unittest("eeq-energy-formamide", test_eeq_e_formamide), & + & new_unittest("eeq-dbdr-co2", test_eeq_dbdr_co2), & + & new_unittest("eeq-dbdL-co2", test_eeq_dbdL_co2), & + & new_unittest("eeq-dadr-ice", test_eeq_dadr_ice), & + & new_unittest("eeq-dadL-ice", test_eeq_dadL_ice), & + & new_unittest("eeq-gradient-co2", test_eeq_g_co2), & + & new_unittest("eeq-sigma-ice", test_eeq_s_ice), & + & new_unittest("eeq-dqdr-urea", test_eeq_dqdr_urea), & + & new_unittest("eeq-dqdL-oxacb", test_eeq_dqdL_oxacb), & + & new_unittest("eeqbc-dbdr-co2", test_eeqbc_dbdr_co2), & + & new_unittest("eeqbc-dbdL-co2", test_eeqbc_dbdL_co2), & + & new_unittest("eeqbc-dadr-ice", test_eeqbc_dadr_ice), & + & new_unittest("eeqbc-dadL-ice", test_eeqbc_dadL_ice), & + & new_unittest("eeqbc-gradient-co2", test_eeqbc_g_co2), & + & new_unittest("eeqbc-sigma-ice", test_eeqbc_s_ice), & + & new_unittest("eeqbc-dqdr-urea", test_eeqbc_dqdr_urea), & + & new_unittest("eeqbc-dqdL-oxacb", test_eeqbc_dqdL_oxacb) & + & ] -end subroutine collect_pbc + end subroutine collect_pbc + subroutine gen_test(error, mol, model, qref, eref) -subroutine gen_test(error, mol, model, qref, eref) + !> Molecular structure data + type(structure_type), intent(in) :: mol - !> Molecular structure data - type(structure_type), intent(in) :: mol + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model - !> Electronegativity equilibration model - class(mchrg_model_type), intent(in) :: model + !> Reference charges + real(wp), intent(in), optional :: qref(:) - !> Reference charges - real(wp), intent(in), optional :: qref(:) + !> Reference energies + real(wp), intent(in), optional :: eref(:) - !> Reference energies - real(wp), intent(in), optional :: eref(:) + !> Error handling + type(error_type), allocatable, intent(out) :: error - !> Error handling - type(error_type), allocatable, intent(out) :: error + real(wp), parameter :: cutoff = 25.0_wp + real(wp), allocatable :: cn(:), qloc(:), trans(:, :) + real(wp), allocatable :: energy(:) + real(wp), allocatable :: qvec(:) - real(wp), parameter :: cutoff = 25.0_wp - real(wp), allocatable :: cn(:), qloc(:), trans(:, :) - real(wp), allocatable :: energy(:) - real(wp), allocatable :: qvec(:) + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) - call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + allocate (cn(mol%nat), qloc(mol%nat)) - allocate(cn(mol%nat), qloc(mol%nat)) + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) - call model%ncoord%get_coordination_number(mol, trans, cn) - call model%local_charge(mol, trans, qloc) + if (present(eref)) then + allocate (energy(mol%nat)) + energy(:) = 0.0_wp + end if + if (present(qref)) then + allocate (qvec(mol%nat)) + end if - if (present(eref)) then - allocate(energy(mol%nat)) - energy(:) = 0.0_wp - end if - if (present(qref)) then - allocate(qvec(mol%nat)) - end if - - call model%solve(mol, error, cn, qloc, energy=energy, qvec=qvec) - if (allocated(error)) return - - if (present(qref)) then - if (any(abs(qvec - qref) > thr)) then - call test_failed(error, "Partial charges do not match") - print'(a)', "Charges:" - print'(3es21.14)', qvec - print'("---")' - print'(3es21.14)', qref - print'("---")' - print'(3es21.14)', qvec - qref + call model%solve(mol, error, cn, qloc, energy=energy, qvec=qvec) + if (allocated(error)) return + + if (present(qref)) then + if (any(abs(qvec - qref) > thr)) then + call test_failed(error, "Partial charges do not match") + print'(a)', "Charges:" + print'(3es21.14)', qvec + print'("---")' + print'(3es21.14)', qref + print'("---")' + print'(3es21.14)', qvec - qref + end if end if - end if - if (allocated(error)) return - - if (present(eref)) then - if (any(abs(energy - eref) > thr)) then - call test_failed(error, "Energies do not match") - print'(a)', "Energy:" - print'(3es21.14)', energy - print'("---")' - print'(3es21.14)', eref - print'("---")' - print'(3es21.14)', energy - eref + if (allocated(error)) return + + if (present(eref)) then + if (any(abs(energy - eref) > thr)) then + call test_failed(error, "Energies do not match") + print'(a)', "Energy:" + print'(3es21.14)', energy + print'("---")' + print'(3es21.14)', eref + print'("---")' + print'(3es21.14)', energy - eref + end if end if - end if -end subroutine gen_test + end subroutine gen_test + subroutine test_numgrad(error, mol, model) -subroutine test_numgrad(error, mol, model) + !> Molecular structure data + type(structure_type), intent(inout) :: mol - !> Molecular structure data - type(structure_type), intent(inout) :: mol + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model - !> Electronegativity equilibration model - class(mchrg_model_type), intent(in) :: model + !> Error handling + type(error_type), allocatable, intent(out) :: error - !> Error handling - type(error_type), allocatable, intent(out) :: error + integer :: iat, ic + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: energy(:), gradient(:, :), sigma(:, :) + real(wp), allocatable :: numgrad(:, :) + real(wp) :: er, el - integer :: iat, ic - real(wp), parameter :: cutoff = 25.0_wp - real(wp), parameter :: step = 1.0e-6_wp - real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) - real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) - real(wp), allocatable :: energy(:), gradient(:, :), sigma(:, :) - real(wp), allocatable :: numgrad(:, :) - real(wp) :: er, el + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) - call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & energy(mol%nat), gradient(3, mol%nat), sigma(3, 3), numgrad(3, mol%nat)) + energy(:) = 0.0_wp + gradient(:, :) = 0.0_wp + sigma(:, :) = 0.0_wp + + lp: do iat = 1, mol%nat + do ic = 1, 3 + energy(:) = 0.0_wp + er = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%solve(mol, error, cn, qloc, energy=energy) + if (allocated(error)) exit lp + er = sum(energy) + + energy(:) = 0.0_wp + el = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%solve(mol, error, cn, qloc, energy=energy) + if (allocated(error)) exit lp + el = sum(energy) + + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + numgrad(ic, iat) = 0.5_wp*(er - el)/step + end do + end do lp + if (allocated(error)) return + + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) - allocate(cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & - & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & - & energy(mol%nat), gradient(3, mol%nat), sigma(3, 3), numgrad(3, mol%nat)) - energy(:) = 0.0_wp - gradient(:, :) = 0.0_wp - sigma(:, :) = 0.0_wp + energy(:) = 0.0_wp + call model%solve(mol, error, cn, qloc, dcndr, dcndL, & + & dqlocdr, dqlocdL, energy, gradient, sigma) + if (allocated(error)) return + + if (any(abs(gradient(:, :) - numgrad(:, :)) > thr2)) then + call test_failed(error, "Derivative of energy does not match") + print'(a)', "gradient:" + print'(3es21.14)', gradient + print'(a)', "numgrad:" + print'(3es21.14)', numgrad + print'(a)', "diff:" + print'(3es21.14)', gradient - numgrad + end if - lp: do iat = 1, mol%nat - do ic = 1, 3 - energy(:) = 0.0_wp - mol%xyz(ic, iat) = mol%xyz(ic, iat) + step - call model%ncoord%get_coordination_number(mol, trans, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, energy=energy) - if (allocated(error)) exit lp - er = sum(energy) + end subroutine test_numgrad - energy(:) = 0.0_wp - mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step - call model%ncoord%get_coordination_number(mol, trans, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, energy=energy) - if (allocated(error)) exit lp - el = sum(energy) - - mol%xyz(ic, iat) = mol%xyz(ic, iat) + step - numgrad(ic, iat) = 0.5_wp*(er - el)/step - end do - end do lp - if (allocated(error)) return + subroutine test_numsigma(error, mol, model) + + !> Molecular structure data + type(structure_type), intent(inout) :: mol - call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) - call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model - energy(:) = 0.0_wp - call model%solve(mol, error, cn, qloc, dcndr, dcndL, & - & dqlocdr, dqlocdL, energy, gradient, sigma) - if (allocated(error)) return + !> Error handling + type(error_type), allocatable, intent(out) :: error - if (any(abs(gradient(:, :) - numgrad(:, :)) > thr2)) then - call test_failed(error, "Derivative of energy does not match") - end if + integer :: ic, jc + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp, unity(3, 3) = reshape(& + & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: energy(:), gradient(:, :) + real(wp), allocatable :: lattr(:, :), xyz(:, :) + real(wp) :: er, el, eps(3, 3), numsigma(3, 3), sigma(3, 3), lattice(3, 3) -end subroutine test_numgrad + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & energy(mol%nat), gradient(3, mol%nat), xyz(3, mol%nat)) + energy(:) = 0.0_wp + gradient(:, :) = 0.0_wp + sigma(:, :) = 0.0_wp + + eps(:, :) = unity + xyz(:, :) = mol%xyz + lattice(:, :) = mol%lattice + lattr = trans + lp: do ic = 1, 3 + do jc = 1, 3 + energy(:) = 0.0_wp + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%solve(mol, error, cn, qloc, energy=energy) + if (allocated(error)) exit lp + er = sum(energy) + + energy(:) = 0.0_wp + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%solve(mol, error, cn, qloc, energy=energy) + if (allocated(error)) exit lp + el = sum(energy) + + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + mol%lattice(:, :) = lattice + lattr(:, :) = trans + numsigma(jc, ic) = 0.5_wp*(er - el)/step + end do + end do lp + if (allocated(error)) return + + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) -subroutine test_numsigma(error, mol, model) + energy(:) = 0.0_wp + call model%solve(mol, error, cn, qloc, dcndr, dcndL, & + & dqlocdr, dqlocdL, energy, gradient, sigma) + if (allocated(error)) return + + if (any(abs(sigma(:, :) - numsigma(:, :)) > thr2)) then + call test_failed(error, "Derivative of energy does not match") + print'(a)', "sigma:" + print'(3es21.14)', sigma + print'(a)', "numgrad:" + print'(3es21.14)', numsigma + print'(a)', "diff:" + print'(3es21.14)', sigma(:, :) - numsigma(:, :) + end if - !> Molecular structure data - type(structure_type), intent(inout) :: mol + end subroutine test_numsigma + + subroutine test_dbdr(error, mol, model) + + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: iat, ic + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: dbdr(:, :, :), dbdL(:, :, :) + real(wp), allocatable :: numgrad(:, :, :), xvecr(:), xvecl(:) + type(cache_container), allocatable :: cache + allocate (cache) + + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & xvecr(mol%nat + 1), xvecl(mol%nat + 1), numgrad(3, mol%nat, mol%nat + 1), & + & dbdr(3, mol%nat, mol%nat + 1), dbdL(3, 3, mol%nat + 1)) + + numgrad = 0.0_wp + + lp: do iat = 1, mol%nat + do ic = 1, 3 + ! Right-hand side + xvecr(:) = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_xvec(mol, cache, xvecr) + + ! Left-hand side + xvecl(:) = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_xvec(mol, cache, xvecl) + + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + numgrad(ic, iat, :) = 0.5_wp*(xvecr(:) - xvecl(:))/step + end do + end do lp + + ! Analytical gradient + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + call model%update(mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) + call model%get_xvec(mol, cache, xvecl) ! need to call this for xtmp in cache (eeqbc) + call model%get_xvec_derivs(mol, cache, dbdr, dbdL) + + if (any(abs(dbdr(:, :, :) - numgrad(:, :, :)) > thr2)) then + call test_failed(error, "Derivative of the b vector does not match") + print'(a)', "dbdr:" + print'(3es21.14)', dbdr + print'(a)', "numgrad:" + print'(3es21.14)', numgrad + print'(a)', "diff:" + print'(3es21.14)', dbdr - numgrad + end if - !> Electronegativity equilibration model - class(mchrg_model_type), intent(in) :: model + end subroutine test_dbdr - !> Error handling - type(error_type), allocatable, intent(out) :: error + subroutine test_dbdL(error, mol, model) - integer :: ic, jc - real(wp), parameter :: cutoff = 25.0_wp - real(wp), parameter :: step = 1.0e-6_wp, unity(3, 3) = reshape(& + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: iat, ic, jc + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp, unity(3, 3) = reshape(& & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) - real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) - real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) - real(wp), allocatable :: energy(:), gradient(:, :) - real(wp), allocatable :: lattr(:, :), xyz(:, :) - real(wp) :: er, el, eps(3, 3), numsigma(3, 3), sigma(3, 3), lattice(3, 3) - - call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) - - allocate(cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & - & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & - & energy(mol%nat), gradient(3, mol%nat), xyz(3, mol%nat)) - energy(:) = 0.0_wp - gradient(:, :) = 0.0_wp - sigma(:, :) = 0.0_wp - - eps(:, :) = unity - xyz(:, :) = mol%xyz - lattice(:, :) = mol%lattice - lattr = trans - lp: do ic = 1, 3 - do jc = 1, 3 - energy(:) = 0.0_wp - eps(jc, ic) = eps(jc, ic) + step - mol%xyz(:, :) = matmul(eps, xyz) - mol%lattice(:, :) = matmul(eps, lattice) - lattr(:, :) = matmul(eps, trans) - call model%ncoord%get_coordination_number(mol, lattr, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, energy=energy) - if (allocated(error)) exit lp - er = sum(energy) + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: dbdr(:, :, :), dbdL(:, :, :) + real(wp), allocatable :: numsigma(:, :, :), xvecr(:), xvecl(:) + real(wp), allocatable :: xyz(:, :), lattr(:, :), trans(:, :) + real(wp) :: lattice(3, 3) + real(wp) :: eps(3, 3) + type(cache_container), allocatable :: cache + allocate (cache) + + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & xvecr(mol%nat + 1), xvecl(mol%nat + 1), numsigma(3, 3, mol%nat + 1), & + & dbdr(3, mol%nat, mol%nat + 1), dbdL(3, 3, mol%nat + 1), xyz(3, mol%nat)) + + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + + numsigma = 0.0_wp + + eps(:, :) = unity + xyz(:, :) = mol%xyz + lattice(:, :) = mol%lattice + lattr = trans + lp: do ic = 1, 3 + do jc = 1, 3 + ! Right-hand side + xvecr(:) = 0.0_wp + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_xvec(mol, cache, xvecr) + + ! Left-hand side + xvecl(:) = 0.0_wp + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_xvec(mol, cache, xvecl) + + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + mol%lattice(:, :) = lattice + lattr(:, :) = trans + do iat = 1, mol%nat + numsigma(jc, ic, iat) = 0.5_wp*(xvecr(iat) - xvecl(iat))/step + end do + end do + end do lp + + ! Analytical gradient + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + call model%update(mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) + call model%get_xvec(mol, cache, xvecl) ! need to call this for xtmp in cache (eeqbc) + call model%get_xvec_derivs(mol, cache, dbdr, dbdL) + + if (any(abs(dbdL(:, :, :) - numsigma(:, :, :)) > thr2)) then + call test_failed(error, "Derivative of the b vector does not match") + print'(a)', "dbdL:" + print'(3es21.14)', dbdL + print'(a)', "numsigma:" + print'(3es21.14)', numsigma + print'(a)', "diff:" + print'(3es21.14)', dbdL - numsigma + end if - energy(:) = 0.0_wp - eps(jc, ic) = eps(jc, ic) - 2*step - mol%xyz(:, :) = matmul(eps, xyz) - mol%lattice(:, :) = matmul(eps, lattice) - lattr(:, :) = matmul(eps, trans) - call model%ncoord%get_coordination_number(mol, lattr, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, energy=energy) - if (allocated(error)) exit lp - el = sum(energy) - - eps(jc, ic) = eps(jc, ic) + step - mol%xyz(:, :) = xyz - mol%lattice(:, :) = lattice - lattr(:, :) = trans - numsigma(jc, ic) = 0.5_wp*(er - el)/step + end subroutine test_dbdL + + subroutine test_dadr(error, mol, model) + + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: iat, ic, jat, kat + real(wp) :: thr2_local + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp + real(wp), allocatable :: cn(:) + real(wp), allocatable :: qloc(:) + real(wp), allocatable :: trans(:, :) + real(wp), allocatable :: dcndr(:, :, :), dcndL(:, :, :), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: dadr(:, :, :), dadL(:, :, :), atrace(:, :) + real(wp), allocatable :: qvec(:), numgrad(:, :, :), amatr(:, :), amatl(:, :), numtrace(:, :) + type(cache_container), allocatable :: cache + allocate (cache) + + allocate (cn(mol%nat), qloc(mol%nat), amatr(mol%nat + 1, mol%nat + 1), amatl(mol%nat + 1, mol%nat + 1), & + & dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), dqlocdr(3, mol%nat, mol%nat), & + & dqlocdL(3, 3, mol%nat), dadr(3, mol%nat, mol%nat + 1), dadL(3, 3, mol%nat + 1), & + & atrace(3, mol%nat), numtrace(3, mol%nat), numgrad(3, mol%nat, mol%nat + 1), qvec(mol%nat)) + + ! Set tolerance higher if testing eeqbc model + select type (model) + type is (eeqbc_model) + thr2_local = 3.0_wp*thr2 + class default + thr2_local = thr2 + end select + + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + + ! Obtain the vector of charges + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%solve(mol, error, cn, qloc, qvec=qvec) + if (allocated(error)) return + + numgrad = 0.0_wp + + lp: do iat = 1, mol%nat + do ic = 1, 3 + ! Right-hand side + amatr(:, :) = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_coulomb_matrix(mol, cache, amatr) + + ! Left-hand side + amatl(:, :) = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_coulomb_matrix(mol, cache, amatl) + + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + + do kat = 1, mol%nat + do jat = 1, mol%nat + ! Numerical gradient of the A matrix + numgrad(ic, iat, kat) = 0.5_wp*qvec(jat)*(amatr(kat, jat) - amatl(kat, jat))/step & + & + numgrad(ic, iat, kat) + end do + end do + end do + end do lp + + ! Analytical gradient + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + call model%update(mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) + call model%get_coulomb_derivs(mol, cache, qvec, dadr, dadL, atrace) + + ! Add trace of the A matrix + do iat = 1, mol%nat + dadr(:, iat, iat) = atrace(:, iat) + dadr(:, iat, iat) end do - end do lp - if (allocated(error)) return - call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) - call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + ! higher tolerance for numerical gradient + if (any(abs(dadr(:, :, :) - numgrad(:, :, :)) > thr2_local)) then + call test_failed(error, "Derivative of the A matrix does not match") + print'(a)', "dadr:" + print'(3es21.14)', dadr + print'(a)', "numgrad:" + print'(3es21.14)', numgrad + print'(a)', "diff:" + print'(3es21.14)', dadr - numgrad + end if + + end subroutine test_dadr - energy(:) = 0.0_wp - call model%solve(mol, error, cn, qloc, dcndr, dcndL, & - & dqlocdr, dqlocdL, energy, gradient, sigma) - if (allocated(error)) return + subroutine test_dadL(error, mol, model) - if (any(abs(sigma(:, :) - numsigma(:, :)) > thr2)) then - call test_failed(error, "Derivative of energy does not match") - end if + !> Molecular structure data + type(structure_type), intent(inout) :: mol -end subroutine test_numsigma + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model + !> Error handling + type(error_type), allocatable, intent(out) :: error -subroutine test_numdqdr(error, mol, model) + integer :: ic, jc, iat + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp, unity(3, 3) = reshape(& + & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: dadr(:, :, :), dadL(:, :, :), atrace(:, :) + real(wp), allocatable :: xyz(:, :), lattr(:, :), trans(:, :) + real(wp), allocatable :: qvec(:), numsigma(:, :, :), amatr(:, :), amatl(:, :) + real(wp) :: lattice(3, 3) + real(wp) :: eps(3, 3) + type(cache_container), allocatable :: cache + allocate (cache) + + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & amatr(mol%nat + 1, mol%nat + 1), amatl(mol%nat + 1, mol%nat + 1), & + & dadr(3, mol%nat, mol%nat + 1), dadL(3, 3, mol%nat + 1), atrace(3, mol%nat), & + & numsigma(3, 3, mol%nat + 1), qvec(mol%nat), xyz(3, mol%nat)) + + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%solve(mol, error, cn, qloc, qvec=qvec) + if (allocated(error)) return + + numsigma = 0.0_wp + + eps(:, :) = unity + xyz(:, :) = mol%xyz + lattice(:, :) = mol%lattice + lattr = trans + lp: do ic = 1, 3 + do jc = 1, 3 + amatr(:, :) = 0.0_wp + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_coulomb_matrix(mol, cache, amatr) + if (allocated(error)) exit lp + + amatl(:, :) = 0.0_wp + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%update(mol, cache, cn, qloc) + call model%get_coulomb_matrix(mol, cache, amatl) + if (allocated(error)) exit lp + + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + mol%lattice(:, :) = lattice + lattr(:, :) = trans + do iat = 1, mol%nat + ! Numerical sigma of the a matrix + numsigma(jc, ic, :) = 0.5_wp*qvec(iat)*(amatr(iat, :) - amatl(iat, :))/step + numsigma(jc, ic, :) + end do + end do + end do lp + if (allocated(error)) return + + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + call model%update(mol, cache, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL) + + call model%get_coulomb_derivs(mol, cache, qvec, dadr, dadL, atrace) + if (allocated(error)) return + + if (any(abs(dadL(:, :, :) - numsigma(:, :, :)) > thr2)) then + call test_failed(error, "Derivative of the A matrix does not match") + print'(a)', "dadL:" + print'(3es21.14)', dadL + print'(a)', "numsigma:" + print'(3es21.14)', numsigma + print'(a)', "diff:" + print'(3es21.14)', dadL - numsigma + end if + + end subroutine test_dadL + + subroutine test_numdqdr(error, mol, model) + + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: iat, ic + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: ql(:), qr(:), dqdr(:, :, :), dqdL(:, :, :) + real(wp), allocatable :: numdr(:, :, :) + + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & ql(mol%nat), qr(mol%nat), dqdr(3, mol%nat, mol%nat), dqdL(3, 3, mol%nat), & + & numdr(3, mol%nat, mol%nat)) + + lp: do iat = 1, mol%nat + do ic = 1, 3 + qr = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%solve(mol, error, cn, qloc, qvec=qr) + if (allocated(error)) exit lp + + ql = 0.0_wp + mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step + call model%ncoord%get_coordination_number(mol, trans, cn) + call model%local_charge(mol, trans, qloc) + call model%solve(mol, error, cn, qloc, qvec=ql) + if (allocated(error)) exit lp + + mol%xyz(ic, iat) = mol%xyz(ic, iat) + step + numdr(ic, iat, :) = 0.5_wp*(qr - ql)/step + end do + end do lp + if (allocated(error)) return + + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + + call model%solve(mol, error, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL, dqdr=dqdr, dqdL=dqdL) + if (allocated(error)) return + + if (any(abs(dqdr(:, :, :) - numdr(:, :, :)) > thr2)) then + call test_failed(error, "Derivative of charges does not match") + print'(a)', "dqdr:" + print'(3es21.14)', dqdr + print'(a)', "numdr:" + print'(3es21.14)', numdr + print'(a)', "diff:" + print'(3es21.14)', dqdr - numdr + end if - !> Molecular structure data - type(structure_type), intent(inout) :: mol + end subroutine test_numdqdr + + subroutine test_numdqdL(error, mol, model) + + !> Molecular structure data + type(structure_type), intent(inout) :: mol + + !> Electronegativity equilibration model + class(mchrg_model_type), intent(in) :: model + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + integer :: ic, jc + real(wp), parameter :: cutoff = 25.0_wp + real(wp), parameter :: step = 1.0e-6_wp, unity(3, 3) = reshape(& + & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) + real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) + real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) + real(wp), allocatable :: qr(:), ql(:), dqdr(:, :, :), dqdL(:, :, :) + real(wp), allocatable :: lattr(:, :), xyz(:, :), numdL(:, :, :) + real(wp) :: eps(3, 3), lattice(3, 3) + + call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + + allocate (cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & + & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & + & qr(mol%nat), ql(mol%nat), dqdr(3, mol%nat, mol%nat), dqdL(3, 3, mol%nat), & + & xyz(3, mol%nat), numdL(3, 3, mol%nat)) + + eps(:, :) = unity + xyz(:, :) = mol%xyz + lattice(:, :) = mol%lattice + lattr = trans + lp: do ic = 1, 3 + do jc = 1, 3 + qr = 0.0_wp + ql = 0.0_wp + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%solve(mol, error, cn, qloc, qvec=qr) + if (allocated(error)) exit lp + + eps(jc, ic) = eps(jc, ic) - 2*step + mol%xyz(:, :) = matmul(eps, xyz) + mol%lattice(:, :) = matmul(eps, lattice) + lattr(:, :) = matmul(eps, trans) + call model%ncoord%get_coordination_number(mol, lattr, cn) + call model%local_charge(mol, lattr, qloc) + call model%solve(mol, error, cn, qloc, qvec=ql) + if (allocated(error)) exit lp + + eps(jc, ic) = eps(jc, ic) + step + mol%xyz(:, :) = xyz + mol%lattice(:, :) = lattice + lattr(:, :) = trans + numdL(jc, ic, :) = 0.5_wp*(qr - ql)/step + end do + end do lp + if (allocated(error)) return + + call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) + call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + + call model%solve(mol, error, cn, qloc, dcndr, dcndL, & + & dqlocdr, dqlocdL, dqdr=dqdr, dqdL=dqdL) + if (allocated(error)) return + + if (any(abs(dqdL(:, :, :) - numdL(:, :, :)) > thr2)) then + call test_failed(error, "Derivative of charges does not match") + print'(a)', "dqdL:" + print'(3es21.14)', dqdL + print'(a)', "numdL:" + print'(3es21.14)', numdL + print'(a)', "diff:" + print'(3es21.14)', dqdL - numdL + end if - !> Electronegativity equilibration model - class(mchrg_model_type), intent(in) :: model + end subroutine test_numdqdL - !> Error handling - type(error_type), allocatable, intent(out) :: error + subroutine test_eeq_q_cyanamide(error) - integer :: iat, ic - real(wp), parameter :: cutoff = 25.0_wp - real(wp), parameter :: step = 1.0e-6_wp - real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) - real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) - real(wp), allocatable :: ql(:), qr(:), dqdr(:, :, :), dqdL(:, :, :) - real(wp), allocatable :: numdr(:, :, :) + !> Error handling + type(error_type), allocatable, intent(out) :: error - call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + real(wp), parameter :: ref(40) = [& + & 3.47274007973765E-1_wp, 3.47273017259661E-1_wp, 3.47276806231462E-1_wp, & + & 3.47267235121650E-1_wp, 3.47288751410308E-1_wp, 3.47296840650318E-1_wp, & + & 3.47299868243747E-1_wp, 3.47299359754750E-1_wp, 3.45885710186856E-1_wp, & + & 3.45874246032015E-1_wp, 3.45888242047875E-1_wp, 3.45877451600398E-1_wp, & + & 3.45902365123333E-1_wp, 3.45902162041418E-1_wp, 3.45900336974539E-1_wp, & + & 3.45903770121699E-1_wp, 3.58996736955583E-1_wp, 3.58987991603151E-1_wp, & + & 3.59000859720791E-1_wp, 3.58990960169615E-1_wp, 3.58990422147650E-1_wp, & + & 3.58997842028106E-1_wp, 3.58997124263551E-1_wp, 3.59001728635141E-1_wp, & + &-5.86637659146611E-1_wp, -5.86601835343945E-1_wp, -5.86625078062709E-1_wp, & + &-5.86579999312624E-1_wp, -5.86654394920376E-1_wp, -5.86669360667441E-1_wp, & + &-5.86667481028248E-1_wp, -5.86674873278027E-1_wp, -4.65529092171338E-1_wp, & + &-4.65527713392387E-1_wp, -4.65546450764424E-1_wp, -4.65540905286785E-1_wp, & + &-4.65529860758785E-1_wp, -4.65533757413217E-1_wp, -4.65527683275370E-1_wp, & + &-4.65527691475100E-1_wp] - allocate(cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & - & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & - & ql(mol%nat), qr(mol%nat), dqdr(3, mol%nat, mol%nat), dqdL(3, 3, mol%nat), & - & numdr(3, mol%nat, mol%nat)) + call get_structure(mol, "X23", "cyanamide") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call gen_test(error, mol, model, qref=ref) - lp: do iat = 1, mol%nat - do ic = 1, 3 - mol%xyz(ic, iat) = mol%xyz(ic, iat) + step - call model%ncoord%get_coordination_number(mol, trans, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, qvec=qr) - if (allocated(error)) exit lp + end subroutine test_eeq_q_cyanamide - mol%xyz(ic, iat) = mol%xyz(ic, iat) - 2*step - call model%ncoord%get_coordination_number(mol, trans, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, qvec=ql) - if (allocated(error)) exit lp + subroutine test_eeq_e_formamide(error) - mol%xyz(ic, iat) = mol%xyz(ic, iat) + step - numdr(ic, iat, :) = 0.5_wp*(qr - ql)/step - end do - end do lp - if (allocated(error)) return + !> Error handling + type(error_type), allocatable, intent(out) :: error - call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) - call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + real(wp), parameter :: ref(24) = [& + & 4.01878116759118E-1_wp, 4.01884779867146E-1_wp, 4.01836430415694E-1_wp, & + & 4.01859412981015E-1_wp, 4.30132679242368E-1_wp, 4.30105708698129E-1_wp, & + & 4.30085668404029E-1_wp, 4.30108945050621E-1_wp, 1.90400193163051E-1_wp, & + & 1.90401822781550E-1_wp, 1.90404130981992E-1_wp, 1.90407447306916E-1_wp, & + & 3.33868335196572E-1_wp, 3.33848682655264E-1_wp, 3.33875977331594E-1_wp, & + & 3.33893030511429E-1_wp, -8.35635584733599E-1_wp, -8.35614926719694E-1_wp, & + &-8.35542684637248E-1_wp, -8.35589488474056E-1_wp, -6.31079122091240E-1_wp, & + &-6.31059677948463E-1_wp, -6.31085206912995E-1_wp, -6.31081747027041E-1_wp] - call model%solve(mol, error, cn, qloc, dcndr, dcndL, dqlocdr, dqlocdL, dqdr=dqdr, dqdL=dqdL) - if (allocated(error)) return + call get_structure(mol, "X23", "formamide") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call gen_test(error, mol, model, eref=ref) - if (any(abs(dqdr(:, :, :) - numdr(:, :, :)) > thr2)) then - call test_failed(error, "Derivative of charges does not match") - end if + end subroutine test_eeq_e_formamide -end subroutine test_numdqdr + subroutine test_eeq_dbdr_co2(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error -subroutine test_numdqdL(error, mol, model) + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model - !> Molecular structure data - type(structure_type), intent(inout) :: mol + call get_structure(mol, "X23", "CO2") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_dbdr(error, mol, model) - !> Electronegativity equilibration model - class(mchrg_model_type), intent(in) :: model + end subroutine test_eeq_dbdr_co2 - !> Error handling - type(error_type), allocatable, intent(out) :: error + subroutine test_eeq_dbdL_co2(error) - integer :: ic, jc - real(wp), parameter :: cutoff = 25.0_wp - real(wp), parameter :: step = 1.0e-6_wp, unity(3, 3) = reshape(& - & [1, 0, 0, 0, 1, 0, 0, 0, 1], shape(unity)) - real(wp), allocatable :: cn(:), dcndr(:, :, :), dcndL(:, :, :), trans(:, :) - real(wp), allocatable :: qloc(:), dqlocdr(:, :, :), dqlocdL(:, :, :) - real(wp), allocatable :: qr(:), ql(:), dqdr(:, :, :), dqdL(:, :, :) - real(wp), allocatable :: lattr(:, :), xyz(:, :), numdL(:, :, :) - real(wp) :: eps(3, 3), lattice(3, 3) - - call get_lattice_points(mol%periodic, mol%lattice, cutoff, trans) - - allocate(cn(mol%nat), dcndr(3, mol%nat, mol%nat), dcndL(3, 3, mol%nat), & - & qloc(mol%nat), dqlocdr(3, mol%nat, mol%nat), dqlocdL(3, 3, mol%nat), & - & qr(mol%nat), ql(mol%nat), dqdr(3, mol%nat, mol%nat), dqdL(3, 3, mol%nat), & - & xyz(3, mol%nat), numdL(3, 3, mol%nat)) - - eps(:, :) = unity - xyz(:, :) = mol%xyz - lattice(:, :) = mol%lattice - lattr = trans - lp: do ic = 1, 3 - do jc = 1, 3 - eps(jc, ic) = eps(jc, ic) + step - mol%xyz(:, :) = matmul(eps, xyz) - mol%lattice(:, :) = matmul(eps, lattice) - lattr(:, :) = matmul(eps, trans) - call model%ncoord%get_coordination_number(mol, lattr, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, qvec=qr) - if (allocated(error)) exit lp - - eps(jc, ic) = eps(jc, ic) - 2*step - mol%xyz(:, :) = matmul(eps, xyz) - mol%lattice(:, :) = matmul(eps, lattice) - lattr(:, :) = matmul(eps, trans) - call model%ncoord%get_coordination_number(mol, lattr, cn) - call model%local_charge(mol, trans, qloc) - call model%solve(mol, error, cn, qloc, qvec=ql) - if (allocated(error)) exit lp - - eps(jc, ic) = eps(jc, ic) + step - mol%xyz(:, :) = xyz - mol%lattice(:, :) = lattice - lattr(:, :) = trans - numdL(jc, ic, :) = 0.5_wp*(qr - ql)/step - end do - end do lp - if (allocated(error)) return + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "X23", "CO2") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_dbdL(error, mol, model) + + end subroutine test_eeq_dbdL_co2 + + subroutine test_eeq_dadr_ice(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "ICE10", "vi") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_dadr(error, mol, model) + + end subroutine test_eeq_dadr_ice + + subroutine test_eeq_dadL_ice(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "ICE10", "vi") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_dadL(error, mol, model) + + end subroutine test_eeq_dadL_ice + + subroutine test_eeq_g_co2(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "X23", "CO2") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_numgrad(error, mol, model) + + end subroutine test_eeq_g_co2 + + subroutine test_eeq_s_ice(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "ICE10", "vi") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_numsigma(error, mol, model) + + end subroutine test_eeq_s_ice + + subroutine test_eeq_dqdr_urea(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "X23", "urea") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_numdqdr(error, mol, model) + + end subroutine test_eeq_dqdr_urea + + subroutine test_eeq_dqdL_oxacb(error) + + !> Error handling + type(error_type), allocatable, intent(out) :: error + + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + + call get_structure(mol, "X23", "oxacb") + call new_eeq2019_model(mol, model, error) + if (allocated(error)) return + call test_numdqdL(error, mol, model) - call model%ncoord%get_coordination_number(mol, trans, cn, dcndr, dcndL) - call model%local_charge(mol, trans, qloc, dqlocdr, dqlocdL) + end subroutine test_eeq_dqdL_oxacb - call model%solve(mol, error, cn, qloc, dcndr, dcndL, & - & dqlocdr, dqlocdL, dqdr=dqdr, dqdL=dqdL) - if (allocated(error)) return + subroutine test_eeqbc_dbdr_co2(error) - if (any(abs(dqdL(:, :, :) - numdL(:, :, :)) > thr2)) then - call test_failed(error, "Derivative of charges does not match") - end if + !> Error handling + type(error_type), allocatable, intent(out) :: error -end subroutine test_numdqdL + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + call get_structure(mol, "X23", "CO2") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_dbdr(error, mol, model) -subroutine test_q_cyanamide(error) + end subroutine test_eeqbc_dbdr_co2 - !> Error handling - type(error_type), allocatable, intent(out) :: error + subroutine test_eeqbc_dbdL_co2(error) - type(structure_type) :: mol - class(mchrg_model_type), allocatable :: model - real(wp), parameter :: ref(40) = [& - & 3.47274007973765E-1_wp, 3.47273017259661E-1_wp, 3.47276806231462E-1_wp, & - & 3.47267235121650E-1_wp, 3.47288751410308E-1_wp, 3.47296840650318E-1_wp, & - & 3.47299868243747E-1_wp, 3.47299359754750E-1_wp, 3.45885710186856E-1_wp, & - & 3.45874246032015E-1_wp, 3.45888242047875E-1_wp, 3.45877451600398E-1_wp, & - & 3.45902365123333E-1_wp, 3.45902162041418E-1_wp, 3.45900336974539E-1_wp, & - & 3.45903770121699E-1_wp, 3.58996736955583E-1_wp, 3.58987991603151E-1_wp, & - & 3.59000859720791E-1_wp, 3.58990960169615E-1_wp, 3.58990422147650E-1_wp, & - & 3.58997842028106E-1_wp, 3.58997124263551E-1_wp, 3.59001728635141E-1_wp, & - &-5.86637659146611E-1_wp,-5.86601835343945E-1_wp,-5.86625078062709E-1_wp, & - &-5.86579999312624E-1_wp,-5.86654394920376E-1_wp,-5.86669360667441E-1_wp, & - &-5.86667481028248E-1_wp,-5.86674873278027E-1_wp,-4.65529092171338E-1_wp, & - &-4.65527713392387E-1_wp,-4.65546450764424E-1_wp,-4.65540905286785E-1_wp, & - &-4.65529860758785E-1_wp,-4.65533757413217E-1_wp,-4.65527683275370E-1_wp, & - &-4.65527691475100E-1_wp] + !> Error handling + type(error_type), allocatable, intent(out) :: error - call get_structure(mol, "X23", "cyanamide") - call new_eeq2019_model(mol, model, error) - if (allocated(error)) return - call gen_test(error, mol, model, qref=ref) + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model -end subroutine test_q_cyanamide + call get_structure(mol, "X23", "CO2") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_dbdL(error, mol, model) + end subroutine test_eeqbc_dbdL_co2 -subroutine test_e_formamide(error) + subroutine test_eeqbc_dadr_ice(error) - !> Error handling - type(error_type), allocatable, intent(out) :: error + !> Error handling + type(error_type), allocatable, intent(out) :: error - type(structure_type) :: mol - class(mchrg_model_type), allocatable :: model - real(wp), parameter :: ref(24) = [& - & 4.01878116759118E-1_wp, 4.01884779867146E-1_wp, 4.01836430415694E-1_wp, & - & 4.01859412981015E-1_wp, 4.30132679242368E-1_wp, 4.30105708698129E-1_wp, & - & 4.30085668404029E-1_wp, 4.30108945050621E-1_wp, 1.90400193163051E-1_wp, & - & 1.90401822781550E-1_wp, 1.90404130981992E-1_wp, 1.90407447306916E-1_wp, & - & 3.33868335196572E-1_wp, 3.33848682655264E-1_wp, 3.33875977331594E-1_wp, & - & 3.33893030511429E-1_wp,-8.35635584733599E-1_wp,-8.35614926719694E-1_wp, & - &-8.35542684637248E-1_wp,-8.35589488474056E-1_wp,-6.31079122091240E-1_wp, & - &-6.31059677948463E-1_wp,-6.31085206912995E-1_wp,-6.31081747027041E-1_wp] + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model - call get_structure(mol, "X23", "formamide") - call new_eeq2019_model(mol, model, error) - if (allocated(error)) return - call gen_test(error, mol, model, eref=ref) + call get_structure(mol, "ICE10", "vi") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_dadr(error, mol, model) -end subroutine test_e_formamide + end subroutine test_eeqbc_dadr_ice + subroutine test_eeqbc_dadL_ice(error) -subroutine test_g_co2(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error - !> Error handling - type(error_type), allocatable, intent(out) :: error + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model - type(structure_type) :: mol - class(mchrg_model_type), allocatable :: model + call get_structure(mol, "ICE10", "vi") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_dadL(error, mol, model) - call get_structure(mol, "X23", "CO2") - call new_eeq2019_model(mol, model, error) - if (allocated(error)) return - call test_numgrad(error, mol, model) + end subroutine test_eeqbc_dadL_ice -end subroutine test_g_co2 + subroutine test_eeqbc_g_co2(error) + !> Error handling + type(error_type), allocatable, intent(out) :: error -subroutine test_s_ice(error) + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model - !> Error handling - type(error_type), allocatable, intent(out) :: error + call get_structure(mol, "X23", "CO2") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_numgrad(error, mol, model) - type(structure_type) :: mol - class(mchrg_model_type), allocatable :: model + end subroutine test_eeqbc_g_co2 - call get_structure(mol, "ICE10", "vi") - call new_eeq2019_model(mol, model, error) - if (allocated(error)) return - call test_numsigma(error, mol, model) + subroutine test_eeqbc_s_ice(error) -end subroutine test_s_ice + !> Error handling + type(error_type), allocatable, intent(out) :: error + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model -subroutine test_dqdr_urea(error) + call get_structure(mol, "ICE10", "vi") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_numsigma(error, mol, model) - !> Error handling - type(error_type), allocatable, intent(out) :: error + end subroutine test_eeqbc_s_ice - type(structure_type) :: mol - class(mchrg_model_type), allocatable :: model + subroutine test_eeqbc_dqdr_urea(error) - call get_structure(mol, "X23", "urea") - call new_eeq2019_model(mol, model, error) - if (allocated(error)) return - call test_numdqdr(error, mol, model) + !> Error handling + type(error_type), allocatable, intent(out) :: error -end subroutine test_dqdr_urea + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model + call get_structure(mol, "X23", "urea") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_numdqdr(error, mol, model) -subroutine test_dqdL_oxacb(error) + end subroutine test_eeqbc_dqdr_urea - !> Error handling - type(error_type), allocatable, intent(out) :: error + subroutine test_eeqbc_dqdL_oxacb(error) - type(structure_type) :: mol - class(mchrg_model_type), allocatable :: model + !> Error handling + type(error_type), allocatable, intent(out) :: error - call get_structure(mol, "X23", "oxacb") - call new_eeq2019_model(mol, model, error) - if (allocated(error)) return - call test_numdqdL(error, mol, model) + type(structure_type) :: mol + class(mchrg_model_type), allocatable :: model -end subroutine test_dqdL_oxacb + call get_structure(mol, "X23", "oxacb") + call new_eeqbc2025_model(mol, model, error) + if (allocated(error)) return + call test_numdqdL(error, mol, model) + end subroutine test_eeqbc_dqdL_oxacb end module test_pbc