(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
14MODULE 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
37CONTAINS
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
324END MODULE pair_potential_util
325
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
mathconstants::fac
real(kind=dp), dimension(0:maxfac), parameter, public fac
Definition mathconstants.F:37
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
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:379