(git:b279b6b)
pair_potential_util.F
Go to the documentation of this file.
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 ! **************************************************************************************************
9 !> \par History
10 !> September 2005 - Introduced the Born-Mayer-Huggins-Fumi-Tosi Potential (BMHTF)
11 !> 2006 - Major rewriting of the routines.. Linear scaling setup of splines
12 !> \author CJM
13 ! **************************************************************************************************
14 MODULE pair_potential_util
15 
16  USE fparser, ONLY: evalerrtype,&
17  evalf
18  USE kinds, ONLY: dp
19  USE mathconstants, ONLY: ifac
20  USE pair_potential_types, ONLY: &
21  b4_type, bm_type, ea_type, ft_type, ftd_type, gp_type, gw_type, ip_type, lj_charmm_type, &
22  lj_type, not_initialized, pair_potential_single_type, tab_type, wl_type
23  USE physcon, ONLY: bohr,&
24  evolt
25 #include "./base/base_uses.f90"
26 
27  IMPLICIT NONE
28 
29  PRIVATE
30 
31  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pair_potential_util'
32  REAL(KIND=dp), PARAMETER, PRIVATE :: min_hicut_value = 1.0e-15_dp, &
33  default_hicut_value = 1.0e3_dp
34 
35  PUBLIC :: ener_pot, ener_zbl, zbl_matching_polinomial
36 
37 CONTAINS
38 
39 ! **************************************************************************************************
40 !> \brief Evaluates the nonbond potential energy for the implemented FF kinds
41 !> \param pot ...
42 !> \param r ...
43 !> \param energy_cutoff ...
44 !> \return ...
45 ! **************************************************************************************************
46  FUNCTION ener_pot(pot, r, energy_cutoff) RESULT(value)
47  TYPE(pair_potential_single_type), POINTER :: pot
48  REAL(kind=dp), INTENT(IN) :: r, energy_cutoff
49  REAL(kind=dp) :: value
50 
51  INTEGER :: i, index, index1, index2, j, n
52  REAL(kind=dp) :: bd6, bd8, dampsum, f6, f8, lvalue, pp, &
53  qq, scale, xf
54 
55  value = 0.0_dp
56  DO j = 1, SIZE(pot%type)
57  ! A lower boundary for the potential definition was defined
58  IF ((pot%set(j)%rmin /= not_initialized) .AND. (r < pot%set(j)%rmin)) cycle
59  ! An upper boundary for the potential definition was defined
60  IF ((pot%set(j)%rmax /= not_initialized) .AND. (r >= pot%set(j)%rmax)) cycle
61  ! If within limits let's compute the potential...
62  IF (pot%type(j) == lj_charmm_type) THEN
63  lvalue = &
64  4.0_dp*pot%set(j)%lj%epsilon*(pot%set(j)%lj%sigma12*r**(-12) - pot%set(j)%lj% &
65  sigma6*r**(-6))
66  ELSE IF (pot%type(j) == lj_type) THEN
67  lvalue = pot%set(j)%lj%epsilon* &
68  (pot%set(j)%lj%sigma12*r**(-12) - pot%set(j)%lj%sigma6*r**(-6))
69  ELSE IF (pot%type(j) == ip_type) THEN
70  lvalue = 0._dp
71  IF (r > pot%set(j)%ipbv%rcore) THEN
72  DO i = 2, 15
73  lvalue = lvalue + pot%set(j)%ipbv%a(i)/(r**(i - 1)*real(i - 1, dp))
74  END DO
75  ELSE
76  ! use a linear potential
77  lvalue = pot%set(j)%ipbv%m*r + pot%set(j)%ipbv%b
78  END IF
79  lvalue = lvalue
80  ELSE IF (pot%type(j) == wl_type) THEN
81  lvalue = pot%set(j)%willis%a*exp(-pot%set(j)%willis%b*r) - pot%set(j)%willis%c/r**6
82  ELSE IF (pot%type(j) == gw_type) THEN
83  scale = exp(pot%set(j)%goodwin%m*(-(r/pot%set(j)%goodwin%dc)**pot%set(j)%goodwin%mc + &
84  (pot%set(j)%goodwin%d/pot%set(j)%goodwin%dc)**pot%set(j)%goodwin%mc))
85  lvalue = scale*pot%set(j)%goodwin%vr0*(pot%set(j)%goodwin%d/r)**pot%set(j)%goodwin%m
86  ELSE IF (pot%type(j) == ft_type) THEN
87  lvalue = pot%set(j)%ft%a*exp(-pot%set(j)%ft%b*r) - pot%set(j)%ft%c/r**6 - pot%set(j)%ft%d/r**8
88  ELSE IF (pot%type(j) == ftd_type) THEN
89  ! Compute 6th order dispersion correction term
90  bd6 = pot%set(j)%ftd%bd(1)
91  dampsum = 1.0_dp
92  xf = 1.0_dp
93  DO i = 1, 6
94  xf = xf*bd6*r
95  dampsum = dampsum + xf*ifac(i)
96  END DO
97  f6 = 1.0_dp - exp(-bd6*r)*dampsum
98  ! Compute 8th order dispersion correction term
99  bd8 = pot%set(j)%ftd%bd(2)
100  dampsum = 1.0_dp
101  xf = 1.0_dp
102  DO i = 1, 8
103  xf = xf*bd8*r
104  dampsum = dampsum + xf*ifac(i)
105  END DO
106  f8 = 1.0_dp - exp(-bd8*r)*dampsum
107  lvalue = pot%set(j)%ftd%a*exp(-pot%set(j)%ftd%b*r) - f6*pot%set(j)%ftd%c/r**6 - f8*pot%set(j)%ftd%d/r**8
108  ELSE IF (pot%type(j) == ea_type) THEN
109  index = int(r/pot%set(j)%eam%drar) + 1
110  IF (index > pot%set(j)%eam%npoints) THEN
111  index = pot%set(j)%eam%npoints
112  ELSEIF (index < 1) THEN
113  index = 1
114  END IF
115  qq = r - pot%set(j)%eam%rval(index)
116  pp = pot%set(j)%eam%phi(index) + &
117  qq*pot%set(j)%eam%phip(index)
118  lvalue = pp
119  ELSE IF (pot%type(j) == b4_type) THEN
120  IF (r <= pot%set(j)%buck4r%r1) THEN
121  pp = pot%set(j)%buck4r%a*exp(-pot%set(j)%buck4r%b*r)
122  ELSEIF (r > pot%set(j)%buck4r%r1 .AND. r <= pot%set(j)%buck4r%r2) THEN
123  pp = 0.0_dp
124  DO n = 0, pot%set(j)%buck4r%npoly1
125  pp = pp + pot%set(j)%buck4r%poly1(n)*r**n
126  END DO
127  ELSEIF (r > pot%set(j)%buck4r%r2 .AND. r <= pot%set(j)%buck4r%r3) THEN
128  pp = 0.0_dp
129  DO n = 0, pot%set(j)%buck4r%npoly2
130  pp = pp + pot%set(j)%buck4r%poly2(n)*r**n
131  END DO
132  ELSEIF (r > pot%set(j)%buck4r%r3) THEN
133  pp = -pot%set(j)%buck4r%c/r**6
134  END IF
135  lvalue = pp
136  ELSE IF (pot%type(j) == tab_type) THEN
137  index1 = floor((r - pot%set(j)%tab%r(1))/pot%set(j)%tab%dr) + 1
138  index2 = index1 + 1
139  IF (index2 > pot%set(j)%tab%npoints) THEN
140  index2 = pot%set(j)%tab%npoints
141  index1 = index2 - 1
142  ELSEIF (index1 < 1) THEN
143  index1 = 1
144  index2 = 2
145  END IF
146  pp = pot%set(j)%tab%e(index1) + (r - pot%set(j)%tab%r(index1))* &
147  (pot%set(j)%tab%e(index2) - pot%set(j)%tab%e(index1))/ &
148  (pot%set(j)%tab%r(index2) - pot%set(j)%tab%r(index1))
149  lvalue = pp
150  ELSE IF (pot%type(j) == bm_type) THEN
151  lvalue = pot%set(j)%buckmo%f0*(pot%set(j)%buckmo%b1 + pot%set(j)%buckmo%b2)* &
152  exp((pot%set(j)%buckmo%a1 + pot%set(j)%buckmo%a2 - r)/(pot%set(j)%buckmo%b1 + pot%set(j)%buckmo%b2)) &
153  - pot%set(j)%buckmo%c/r**6 &
154  + pot%set(j)%buckmo%d*(exp(-2._dp*pot%set(j)%buckmo%beta*(r - pot%set(j)%buckmo%r0)) - &
155  2.0_dp*exp(-pot%set(j)%buckmo%beta*(r - pot%set(j)%buckmo%r0)))
156  ELSE IF (pot%type(j) == gp_type) THEN
157  pot%set(j)%gp%values(1) = r
158  lvalue = evalf(pot%set(j)%gp%myid, pot%set(j)%gp%values)
159  IF (evalerrtype > 0) &
160  cpabort("Error evaluating generic potential energy function")
161  ELSE
162  lvalue = 0.0_dp
163  END IF
164  value = value + lvalue
165  END DO
166  value = value - energy_cutoff
167  END FUNCTION ener_pot
168 
169 ! **************************************************************************************************
170 !> \brief Evaluates the ZBL scattering potential, very short range
171 !> Only shell-model for interactions among pairs without repulsive term
172 !> \param pot ...
173 !> \param r ...
174 !> \return ...
175 !> \author i soliti ignoti
176 ! **************************************************************************************************
177  FUNCTION ener_zbl(pot, r)
178 
179  TYPE(pair_potential_single_type), POINTER :: pot
180  REAL(kind=dp), INTENT(IN) :: r
181  REAL(kind=dp) :: ener_zbl
182 
183  REAL(kind=dp) :: au, fac, x
184 
185  ener_zbl = 0.0_dp
186  IF (r <= pot%zbl_rcut(1)) THEN
187  au = 0.88534_dp*bohr/(pot%z1**0.23_dp + pot%z2**0.23_dp)
188  x = r/au
189  fac = pot%z1*pot%z2/evolt
190  ener_zbl = fac/r*(0.1818_dp*exp(-3.2_dp*x) + 0.5099_dp*exp(-0.9423_dp*x) + &
191  0.2802_dp*exp(-0.4029_dp*x) + 0.02817_dp*exp(-0.2016_dp*x))
192  ELSEIF (r > pot%zbl_rcut(1) .AND. r <= pot%zbl_rcut(2)) THEN
193  ener_zbl = pot%zbl_poly(0) + pot%zbl_poly(1)*r + pot%zbl_poly(2)*r*r + pot%zbl_poly(3)*r*r*r + &
194  pot%zbl_poly(4)*r*r*r*r + pot%zbl_poly(5)*r*r*r*r*r
195  ELSE
196  ener_zbl = 0.0_dp
197  END IF
198 
199  END FUNCTION ener_zbl
200 
201 ! **************************************************************************************************
202 !> \brief Determine the polinomial coefficients used to set to zero the zbl potential
203 !> at the cutoff radius, with continuity in function, first and second derivative
204 !> Only shell-model for interactions among pairs without repulsive term
205 !> \param pot ...
206 !> \param rcov1 ...
207 !> \param rcov2 ...
208 !> \param z1 ...
209 !> \param z2 ...
210 !> \author i soliti ignoti
211 ! **************************************************************************************************
212  SUBROUTINE zbl_matching_polinomial(pot, rcov1, rcov2, z1, z2)
213 
214  TYPE(pair_potential_single_type), POINTER :: pot
215  REAL(kind=dp), INTENT(IN) :: rcov1, rcov2, z1, z2
216 
217  REAL(kind=dp) :: au, d1, d2, dd1, dd2, fac, v1, v2, x, &
218  x1, x2
219 
220  pot%zbl_rcut(1) = (rcov1 + rcov2)*(1.0_dp - 0.2_dp)*bohr
221  pot%zbl_rcut(2) = (rcov1 + rcov2)*bohr
222  x1 = pot%zbl_rcut(1)
223  x2 = pot%zbl_rcut(2)
224  pot%z1 = z1
225  pot%z2 = z2
226 
227  au = 0.88534_dp*bohr/(z1**0.23_dp + z2**0.23_dp)
228  x = x1/au
229  fac = z1*z2/evolt
230  v1 = fac/x1*(0.1818_dp*exp(-3.2_dp*x) + 0.5099_dp*exp(-0.9423_dp*x) + &
231  0.2802_dp*exp(-0.4029_dp*x) + 0.02817_dp*exp(-0.2016_dp*x))
232  d1 = fac/x1/au*(-3.2_dp*0.1818_dp*exp(-3.2_dp*x) - 0.9423_dp*0.5099_dp*exp(-0.9423_dp*x) &
233  - 0.4029_dp*0.2802_dp*exp(-0.4029_dp*x) - 0.2016_dp*0.02817_dp*exp(-0.2016_dp*x)) &
234  - fac/x1/x1*(0.1818_dp*exp(-3.2_dp*x) + 0.5099_dp*exp(-0.9423_dp*x) + &
235  0.2802_dp*exp(-0.4029_dp*x) + 0.02817_dp*exp(-0.2016_dp*x))
236 
237  dd1 = 2.0_dp*fac/x1**3*(0.1818_dp*exp(-0.32e1_dp*x) &
238  + 0.5099_dp*exp(-0.9423_dp*x) + 0.2802_dp*exp(-0.4029_dp*x) &
239  + 0.2817e-1_dp*exp(-0.2016_dp*x)) &
240  - 0.2e1_dp*fac/x1**2/au*(-0.58176_dp*exp(-0.32e1_dp*x) - 0.48047877_dp*exp(-0.9423_dp*x) &
241  - 0.11289258_dp*exp(-0.4029_dp*x) - 0.5679072e-2_dp*exp(-0.2016_dp*x)) &
242  + fac/x1/au**2*(0.1861632e1_dp*exp(-0.32e1_dp*x) + &
243  0.4527551450_dp*exp(-0.9423_dp*x) + 0.4548442048e-1_dp*exp(-0.4029_dp*x) + &
244  0.1144900915e-2_dp*exp(-0.2016_dp*x))
245 
246  v2 = 0.0_dp
247  d2 = 0.0_dp
248  dd2 = 0.0_dp
249 
250  CALL compute_polinomial_5th(x1, v1, d1, dd1, x2, v2, d2, dd2, pot%zbl_poly)
251 
252  END SUBROUTINE zbl_matching_polinomial
253 
254 ! **************************************************************************************************
255 !> \brief ...
256 !> \param r1 ...
257 !> \param v1 ...
258 !> \param d1 ...
259 !> \param dd1 ...
260 !> \param r2 ...
261 !> \param v2 ...
262 !> \param d2 ...
263 !> \param dd2 ...
264 !> \param poly ...
265 ! **************************************************************************************************
266  SUBROUTINE compute_polinomial_5th(r1, v1, d1, dd1, r2, v2, d2, dd2, poly)
267 
268  REAL(kind=dp) :: r1, v1, d1, dd1, r2, v2, d2, dd2, &
269  poly(0:5)
270 
271  REAL(kind=dp) :: a0, a1, a2, a3, a4, a5
272 
273 ! 5th order
274 
275  a0 = .5_dp*(2._dp*r1**5*v2 - 2._dp*v1*r2**5 + 10._dp*v1*r2**4*r1 - 20._dp*v1*r1**2*r2**3 - r1**2*dd1*r2**5 - &
276  r1**4*r2**3*dd1 + 20._dp*r1**3*r2**2*v2 + 2._dp*r1**3*r2**4*dd1 + r1**3*r2**4*dd2 - 8._dp*r1**3*r2**3*d2 - &
277  2._dp*r1**4*r2**3*dd2 + 10._dp*r1**4*r2**2*d2 - 10._dp*r1**4*r2*v2 + 2._dp*r1*d1*r2**5 - 10._dp*r1**2*d1*r2**4 + &
278  8._dp*r1**3*d1*r2**3 - 2._dp*r1**5*r2*d2 + r1**5*r2**2*dd2)/ &
279  (10.*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
280 
281  a1 = -.5_dp*(-4._dp*r2**3*r1**3*dd2 + 24._dp*r2**2*r1**3*d1 + 4._dp*r2**3*r1**3*dd1 + 3._dp*r2**4*r1**2*dd2 + &
282  r2**4*r1**2*dd1 - 2._dp*r2**5*r1*dd1 - 10._dp*r2**4*r1*d1 + 10._dp*r1**4*r2*d2 - &
283  r1**4*r2**2*dd2 - 3._dp*r1**4*r2**2*dd1 + &
284  2._dp*r1**5*r2*dd2 - 24._dp*r2**3*r1**2*d2 - 16._dp*r2**3*r1**2*d1 + &
285  16._dp*r2**2*r1**3*d2 - 2._dp*r1**5*d2 + 2._dp*r2**5*d1 - &
286  60._dp*r1**2*r2**2*v1 + 60._dp*r1**2*r2**2*v2)/ &
287  (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
288 
289  a2 = .5_dp*(60._dp*r1**2*r2*v2 - 60._dp*v1*r1*r2**2 - 12._dp*r1**2*r2**2*d2 - 36._dp*r1*d1*r2**3 + 3._dp*r2**4*r1*dd2 - &
290  24._dp*r2**3*r1*d2 - 4._dp*r2**4*r1*dd1 + 12._dp*r1**2*r2**2*d1 - 8._dp*r1**3*r2**2*dd2 + 24._dp*r1**3*r2*d1 + &
291  4._dp*r1**4*r2*dd2 + 36._dp*r1**3*r2*d2 - 3._dp*r1**4*r2*dd1 + 8._dp*r2**3*r1**2*dd1 + 60._dp*r2**2*r1*v2 - &
292  60._dp*r1**2*v1*r2 + r1**5*dd2 - r2**5*dd1)/ &
293  (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
294 
295  a3 = -.5_dp*(3._dp*r1**4*dd2 - r1**4*dd1 + 8.*r1**3*d1 - 4.*r1**3*r2*dd1 + &
296  12._dp*r1**3*d2 + 32._dp*r1**2*r2*d1 - 8._dp*r1**2*r2**2*dd2 - &
297  20._dp*r1**2*v1 + 8._dp*r1**2*r2**2*dd1 + 28._dp*r1**2*r2*d2 + &
298  20._dp*r1**2*v2 + 80._dp*r1*r2*v2 - 28._dp*r2**2*r1*d1 - 80._dp*r1*v1*r2 - &
299  32._dp*r2**2*r1*d2 + 4._dp*r1*r2**3*dd2 - 8._dp*r2**3*d2 - 12._dp*r2**3*d1 + &
300  r2**4*dd2 - 3._dp*r2**4*dd1 + 20._dp*r2**2*v2 - 20._dp*r2**2*v1)/ &
301  (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
302 
303  a4 = .5_dp*(3._dp*r1**3*dd2 - 2._dp*r1**3*dd1 + r1**2*r2*dd1 + 14.*r1**2*d1 - 4._dp*r1**2*r2*dd2 + &
304  16._dp*r1**2*d2 - 2._dp*r1*r2*d2 - r1*r2**2*dd2 - &
305  30._dp*r1*v1 + 30.*r1*v2 + 2._dp*r1*r2*d1 + 4._dp*r1*r2**2*dd1 - 16._dp*r2**2*d1 + &
306  2._dp*r2**3*dd2 - 14._dp*r2**2*d2 + 30._dp*r2*v2 - 30._dp*v1*r2 - &
307  3._dp*r2**3*dd1)/(10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - &
308  10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
309 
310  a5 = -.5_dp*(6._dp*r1*d1 + 2._dp*r2*r1*dd1 + 6._dp*r1*d2 - 2.*r2*r1*dd2 - &
311  r2**2*dd1 - r1**2*dd1 - 12.*v1 + 12._dp*v2 + r1**2*dd2 - &
312  6._dp*r2*d1 + r2**2*dd2 - 6._dp*r2*d2)/ &
313  (10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
314 
315  poly(0) = a0
316  poly(1) = a1
317  poly(2) = a2
318  poly(3) = a3
319  poly(4) = a4
320  poly(5) = a5
321 
322  END SUBROUTINE compute_polinomial_5th
323 
324 END MODULE pair_potential_util
325 
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
fparser
This public domain function parser module is intended for applications where a set of mathematical ex...
Definition: fparser.F:17
fparser::evalf
real(rn) function, public evalf(i, Val)
...
Definition: fparser.F:180
fparser::evalerrtype
integer, public evalerrtype
Definition: fparser.F:32
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
mathconstants::ifac
real(kind=dp), dimension(0:maxfac), parameter, public ifac
Definition: mathconstants.F:49
pair_potential_types
Definition: pair_potential_types.F:14
pair_potential_types::lj_charmm_type
integer, parameter, public lj_charmm_type
Definition: pair_potential_types.F:35
pair_potential_types::bm_type
integer, parameter, public bm_type
Definition: pair_potential_types.F:35
pair_potential_types::wl_type
integer, parameter, public wl_type
Definition: pair_potential_types.F:35
pair_potential_types::ft_type
integer, parameter, public ft_type
Definition: pair_potential_types.F:35
pair_potential_types::tab_type
integer, parameter, public tab_type
Definition: pair_potential_types.F:35
pair_potential_types::ftd_type
integer, parameter, public ftd_type
Definition: pair_potential_types.F:35
pair_potential_types::ip_type
integer, parameter, public ip_type
Definition: pair_potential_types.F:35
pair_potential_types::lj_type
integer, parameter, public lj_type
Definition: pair_potential_types.F:35
pair_potential_types::gp_type
integer, parameter, public gp_type
Definition: pair_potential_types.F:35
pair_potential_types::gw_type
integer, parameter, public gw_type
Definition: pair_potential_types.F:35
pair_potential_types::not_initialized
real(kind=dp), parameter, public not_initialized
Definition: pair_potential_types.F:88
pair_potential_types::b4_type
integer, parameter, public b4_type
Definition: pair_potential_types.F:35
pair_potential_types::ea_type
integer, parameter, public ea_type
Definition: pair_potential_types.F:35
pair_potential_util
Definition: pair_potential_util.F:14
pair_potential_util::ener_zbl
real(kind=dp) function, public ener_zbl(pot, r)
Evaluates the ZBL scattering potential, very short range Only shell-model for interactions among pair...
Definition: pair_potential_util.F:178
pair_potential_util::ener_pot
real(kind=dp) function, public ener_pot(pot, r, energy_cutoff)
Evaluates the nonbond potential energy for the implemented FF kinds.
Definition: pair_potential_util.F:47
pair_potential_util::zbl_matching_polinomial
subroutine, public zbl_matching_polinomial(pot, rcov1, rcov2, z1, z2)
Determine the polinomial coefficients used to set to zero the zbl potential at the cutoff radius,...
Definition: pair_potential_util.F:213
physcon
Definition of physical constants:
Definition: physcon.F:68
physcon::evolt
real(kind=dp), parameter, public evolt
Definition: physcon.F:183
physcon::bohr
real(kind=dp), parameter, public bohr
Definition: physcon.F:147