Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 851784930 : 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 851784930 : r = SQRT(r2)
51 851784930 : 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 851772120 : IF (ewald_type == do_ewald_none) THEN
66 102856 : potential_coulomb = 1.0_dp/r
67 102856 : fscalar = fscalar + qfac*potential_coulomb/r2
68 : ELSE
69 851669264 : x1 = alpha*r
70 851669264 : potential_coulomb = erfc(x1)/r
71 : fscalar = fscalar + qfac*(potential_coulomb + &
72 851669264 : two_over_sqrt_pi*EXP(-x1*x1)*alpha)/r2
73 : END IF
74 : END IF
75 :
76 851784930 : potential_coulomb = qfac*(potential_coulomb - interaction_cutoff)
77 :
78 851784930 : END FUNCTION potential_coulomb
79 :
80 : END MODULE pair_potential_coulomb
|