! This file is part of dftd4. ! SPDX-Identifier: LGPL-3.0-or-later ! ! dftd4 is free software: you can redistribute it and/or modify it under ! the terms of the Lesser GNU General Public License as published by ! the Free Software Foundation, either version 3 of the License, or ! (at your option) any later version. ! ! dftd4 is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! Lesser GNU General Public License for more details. ! ! You should have received a copy of the Lesser GNU General Public License ! along with dftd4. If not, see <https://www.gnu.org/licenses/>. module dftd4_ncoord use mctc_env, only : wp use mctc_io, only : structure_type use mctc_io_constants, only : pi implicit none private public :: get_coordination_number, add_coordination_number_derivs !> Steepness of counting function real(wp), parameter :: kcn = 7.5_wp !> Parameter for electronegativity scaling real(wp), parameter :: k4 = 4.10451_wp !> Parameter for electronegativity scaling real(wp), parameter :: k5 = 19.08857_wp !> Parameter for electronegativity scaling real(wp), parameter :: k6 = 2*11.28174_wp**2 contains !> Geometric fractional coordination number, supports error function counting. subroutine get_coordination_number(mol, trans, cutoff, rcov, en, cn, dcndr, dcndL) !DEC$ ATTRIBUTES DLLEXPORT :: get_coordination_number !> Molecular structure data type(structure_type), intent(in) :: mol !> Lattice points real(wp), intent(in) :: trans(:, :) !> Real space cutoff real(wp), intent(in) :: cutoff !> Covalent radius real(wp), intent(in) :: rcov(:) !> Electronegativity real(wp), intent(in) :: en(:) !> Error function coordination number. real(wp), intent(out) :: cn(:) !> Derivative of the CN with respect to the Cartesian coordinates. real(wp), intent(out), optional :: dcndr(:, :, :) !> Derivative of the CN with respect to strain deformations. real(wp), intent(out), optional :: dcndL(:, :, :) if (present(dcndr) .and. present(dcndL)) then call ncoord_derf(mol, trans, cutoff, rcov, en, cn, dcndr, dcndL) else call ncoord_erf(mol, trans, cutoff, rcov, en, cn) end if end subroutine get_coordination_number subroutine ncoord_erf(mol, trans, cutoff, rcov, en, cn) !> Molecular structure data type(structure_type), intent(in) :: mol !> Lattice points real(wp), intent(in) :: trans(:, :) !> Real space cutoff real(wp), intent(in) :: cutoff !> Covalent radius real(wp), intent(in) :: rcov(:) !> Electronegativity real(wp), intent(in) :: en(:) !> Error function coordination number. real(wp), intent(out) :: cn(:) integer :: iat, jat, izp, jzp, itr real(wp) :: r2, r1, rc, rij(3), countf, cutoff2, den cn(:) = 0.0_wp cutoff2 = cutoff**2 !$omp parallel do default(none) reduction(+:cn) & !$omp shared(mol, trans, cutoff2, rcov, en) & !$omp private(jat, itr, izp, jzp, r2, rij, r1, rc, countf, den) do iat = 1, mol%nat izp = mol%id(iat) do jat = 1, iat jzp = mol%id(jat) den = k4*exp(-(abs(en(izp)-en(jzp)) + k5)**2/k6) do itr = 1, size(trans, dim=2) rij = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, itr)) r2 = sum(rij**2) if (r2 > cutoff2 .or. r2 < 1.0e-12_wp) cycle r1 = sqrt(r2) rc = rcov(izp) + rcov(jzp) countf = den*erf_count(kcn, r1, rc) cn(iat) = cn(iat) + countf if (iat /= jat) then cn(jat) = cn(jat) + countf end if end do end do end do end subroutine ncoord_erf subroutine ncoord_derf(mol, trans, cutoff, rcov, en, cn, dcndr, dcndL) !> Molecular structure data type(structure_type), intent(in) :: mol !> Lattice points real(wp), intent(in) :: trans(:, :) !> Real space cutoff real(wp), intent(in) :: cutoff !> Covalent radius real(wp), intent(in) :: rcov(:) !> Electronegativity real(wp), intent(in) :: en(:) !> Error function coordination number. real(wp), intent(out) :: cn(:) !> Derivative of the CN with respect to the Cartesian coordinates. real(wp), intent(out) :: dcndr(:, :, :) !> Derivative of the CN with respect to strain deformations. real(wp), intent(out) :: dcndL(:, :, :) integer :: iat, jat, izp, jzp, itr real(wp) :: r2, r1, rc, rij(3), countf, countd(3), sigma(3, 3), cutoff2, den cn(:) = 0.0_wp dcndr(:, :, :) = 0.0_wp dcndL(:, :, :) = 0.0_wp cutoff2 = cutoff**2 !$omp parallel do schedule(runtime) default(none) & !$omp reduction(+:cn, dcndr, dcndL) & !$omp shared(mol, trans, cutoff2, rcov, en) & !$omp private(jat, itr, izp, jzp, r2, rij, r1, rc, countf, countd, sigma, den) do iat = 1, mol%nat izp = mol%id(iat) do jat = 1, iat jzp = mol%id(jat) den = k4*exp(-(abs(en(izp)-en(jzp)) + k5)**2/k6) do itr = 1, size(trans, dim=2) rij = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, itr)) r2 = sum(rij**2) if (r2 > cutoff2 .or. r2 < 1.0e-12_wp) cycle r1 = sqrt(r2) rc = rcov(izp) + rcov(jzp) countf = den*erf_count(kcn, r1, rc) countd = den*derf_count(kcn, r1, rc) * rij/r1 cn(iat) = cn(iat) + countf if (iat /= jat) then cn(jat) = cn(jat) + countf end if dcndr(:, iat, iat) = dcndr(:, iat, iat) + countd dcndr(:, jat, jat) = dcndr(:, jat, jat) - countd dcndr(:, iat, jat) = dcndr(:, iat, jat) + countd dcndr(:, jat, iat) = dcndr(:, jat, iat) - countd sigma = spread(countd, 1, 3) * spread(rij, 2, 3) dcndL(:, :, iat) = dcndL(:, :, iat) + sigma if (iat /= jat) then dcndL(:, :, jat) = dcndL(:, :, jat) + sigma end if end do end do end do end subroutine ncoord_derf subroutine add_coordination_number_derivs(mol, trans, cutoff, rcov, en, dEdcn, gradient, sigma) !> Molecular structure data type(structure_type), intent(in) :: mol !> Lattice points real(wp), intent(in) :: trans(:, :) !> Real space cutoff real(wp), intent(in) :: cutoff !> Covalent radius real(wp), intent(in) :: rcov(:) !> Electronegativity real(wp), intent(in) :: en(:) !> Derivative of expression with respect to the coordination number real(wp), intent(in) :: dEdcn(:) !> Derivative of the CN with respect to the Cartesian coordinates real(wp), intent(inout) :: gradient(:, :) !> Derivative of the CN with respect to strain deformations real(wp), intent(inout) :: sigma(:, :) integer :: iat, jat, izp, jzp, itr real(wp) :: r2, r1, rc, rij(3), countf, countd(3), ds(3, 3), cutoff2, den cutoff2 = cutoff**2 !$omp parallel do schedule(runtime) default(none) & !$omp reduction(+:gradient, sigma) & !$omp shared(mol, trans, cutoff2, rcov, en, dEdcn) & !$omp private(iat, jat, itr, izp, jzp, r2, rij, r1, rc, countd, ds, den) do iat = 1, mol%nat izp = mol%id(iat) do jat = 1, iat jzp = mol%id(jat) den = k4*exp(-(abs(en(izp)-en(jzp)) + k5)**2/k6) do itr = 1, size(trans, dim=2) rij = mol%xyz(:, iat) - (mol%xyz(:, jat) + trans(:, itr)) r2 = sum(rij**2) if (r2 > cutoff2 .or. r2 < 1.0e-12_wp) cycle r1 = sqrt(r2) rc = rcov(izp) + rcov(jzp) countd = den*derf_count(kcn, r1, rc) * rij/r1 gradient(:, iat) = gradient(:, iat) + countd * (dEdcn(iat) + dEdcn(jat)) gradient(:, jat) = gradient(:, jat) - countd * (dEdcn(iat) + dEdcn(jat)) ds = spread(countd, 1, 3) * spread(rij, 2, 3) sigma(:, :) = sigma(:, :) & & + ds * (dEdcn(iat) + merge(dEdcn(jat), 0.0_wp, jat /= iat)) end do end do end do end subroutine add_coordination_number_derivs !> Error function counting function for coordination number contributions. pure function erf_count(k, r, r0) result(count) !> Steepness of the counting function. real(wp), intent(in) :: k !> Current distance. real(wp), intent(in) :: r !> Cutoff radius. real(wp), intent(in) :: r0 real(wp) :: count count = 0.5_wp * (1.0_wp + erf(-k*(r-r0)/r0)) end function erf_count !> Derivative of the counting function w.r.t. the distance. pure function derf_count(k, r, r0) result(count) !> Steepness of the counting function. real(wp), intent(in) :: k !> Current distance. real(wp), intent(in) :: r !> Cutoff radius. real(wp), intent(in) :: r0 real(wp), parameter :: sqrtpi = sqrt(pi) real(wp) :: count count = -k/sqrtpi/r0*exp(-k**2*(r-r0)**2/r0**2) end function derf_count end module dftd4_ncoord