Line data Source code
1 : !--------------------------------------------------------------------------------------------------! 2 : ! CP2K: A general program to perform molecular dynamics simulations ! 3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> ! 4 : ! ! 5 : ! SPDX-License-Identifier: GPL-2.0-or-later ! 6 : !--------------------------------------------------------------------------------------------------! 7 : 8 : MODULE pair_potential_coulomb 9 : 10 : USE kinds, ONLY: dp 11 : USE mathconstants, ONLY: oorootpi 12 : USE pw_poisson_types, ONLY: do_ewald_none 13 : #include "./base/base_uses.f90" 14 : 15 : IMPLICIT NONE 16 : 17 : PRIVATE 18 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pair_potential_coulomb' 19 : 20 : PUBLIC :: potential_coulomb 21 : 22 : CONTAINS 23 : 24 : ! ************************************************************************************************** 25 : !> \brief Evaluates the electrostatic energy and force 26 : !> \param r2 ... 27 : !> \param fscalar ... 28 : !> \param qfac ... 29 : !> \param ewald_type ... 30 : !> \param alpha ... 31 : !> \param beta ... 32 : !> \param interaction_cutoff ... 33 : !> \return ... 34 : !> \author Toon.Verstraelen@gmail.com 35 : ! ************************************************************************************************** 36 849329242 : FUNCTION potential_coulomb(r2, fscalar, qfac, ewald_type, alpha, beta, & 37 : interaction_cutoff) 38 : 39 : REAL(KIND=dp), INTENT(IN) :: r2 40 : REAL(KIND=dp), INTENT(INOUT) :: fscalar 41 : REAL(KIND=dp), INTENT(IN) :: qfac 42 : INTEGER, INTENT(IN) :: ewald_type 43 : REAL(KIND=dp), INTENT(IN) :: alpha, beta, interaction_cutoff 44 : REAL(KIND=dp) :: potential_coulomb 45 : 46 : REAL(KIND=dp), PARAMETER :: two_over_sqrt_pi = 2.0_dp*oorootpi 47 : 48 : REAL(KIND=dp) :: r, x1, x2 49 : 50 849329242 : r = SQRT(r2) 51 849329242 : IF (beta > 0.0_dp) THEN 52 12810 : IF (ewald_type == do_ewald_none) THEN 53 12676 : x2 = r*beta 54 12676 : potential_coulomb = erf(x2)/r 55 : fscalar = fscalar + qfac*(potential_coulomb - & 56 12676 : two_over_sqrt_pi*EXP(-x2*x2)*beta)/r2 57 : ELSE 58 134 : x1 = alpha*r 59 134 : x2 = r*beta 60 134 : potential_coulomb = (erf(x2) - erf(x1))/r 61 : fscalar = fscalar + qfac*(potential_coulomb + & 62 134 : two_over_sqrt_pi*(EXP(-x1*x1)*alpha - EXP(-x2*x2)*beta))/r2 63 : END IF 64 : ELSE 65 849316432 : IF (ewald_type == do_ewald_none) THEN 66 102784 : potential_coulomb = 1.0_dp/r 67 102784 : fscalar = fscalar + qfac*potential_coulomb/r2 68 : ELSE 69 849213648 : x1 = alpha*r 70 849213648 : potential_coulomb = erfc(x1)/r 71 : fscalar = fscalar + qfac*(potential_coulomb + & 72 849213648 : two_over_sqrt_pi*EXP(-x1*x1)*alpha)/r2 73 : END IF 74 : END IF 75 : 76 849329242 : potential_coulomb = qfac*(potential_coulomb - interaction_cutoff) 77 : 78 849329242 : END FUNCTION potential_coulomb 79 : 80 : END MODULE pair_potential_coulomb