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   849329631 :    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   849329631 :       r = SQRT(r2)
      51   849329631 :       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   849316821 :          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   849214037 :             x1 = alpha*r
      70   849214037 :             potential_coulomb = erfc(x1)/r
      71             :             fscalar = fscalar + qfac*(potential_coulomb + &
      72   849214037 :                                       two_over_sqrt_pi*EXP(-x1*x1)*alpha)/r2
      73             :          END IF
      74             :       END IF
      75             : 
      76   849329631 :       potential_coulomb = qfac*(potential_coulomb - interaction_cutoff)
      77             : 
      78   849329631 :    END FUNCTION potential_coulomb
      79             : 
      80             : END MODULE pair_potential_coulomb