(git:a145afa)
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-2026 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 ELSE IF (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 ELSE IF (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 ELSE IF (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 ELSE IF (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 ELSE IF (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) THEN
160 cpabort("Error evaluating generic potential energy function")
161 END IF
162 ELSE
163 lvalue = 0.0_dp
164 END IF
165 value = value + lvalue
166 END DO
167 value = value - energy_cutoff
168 END FUNCTION ener_pot
169
170! **************************************************************************************************
171!> \brief Evaluates the ZBL scattering potential, very short range
172!> Only shell-model for interactions among pairs without repulsive term
173!> \param pot ...
174!> \param r ...
175!> \return ...
176!> \author i soliti ignoti
177! **************************************************************************************************
178 FUNCTION ener_zbl(pot, r)
179
180 TYPE(pair_potential_single_type), POINTER :: pot
181 REAL(kind=dp), INTENT(IN) :: r
182 REAL(kind=dp) :: ener_zbl
183
184 REAL(kind=dp) :: au, fac, x
185
186 ener_zbl = 0.0_dp
187 IF (r <= pot%zbl_rcut(1)) THEN
188 au = 0.88534_dp*bohr/(pot%z1**0.23_dp + pot%z2**0.23_dp)
189 x = r/au
190 fac = pot%z1*pot%z2/evolt
191 ener_zbl = fac/r*(0.1818_dp*exp(-3.2_dp*x) + 0.5099_dp*exp(-0.9423_dp*x) + &
192 0.2802_dp*exp(-0.4029_dp*x) + 0.02817_dp*exp(-0.2016_dp*x))
193 ELSE IF (r > pot%zbl_rcut(1) .AND. r <= pot%zbl_rcut(2)) THEN
194 ener_zbl = pot%zbl_poly(0) + pot%zbl_poly(1)*r + pot%zbl_poly(2)*r*r + pot%zbl_poly(3)*r*r*r + &
195 pot%zbl_poly(4)*r*r*r*r + pot%zbl_poly(5)*r*r*r*r*r
196 ELSE
197 ener_zbl = 0.0_dp
198 END IF
199
200 END FUNCTION ener_zbl
201
202! **************************************************************************************************
203!> \brief Determine the polinomial coefficients used to set to zero the zbl potential
204!> at the cutoff radius, with continuity in function, first and second derivative
205!> Only shell-model for interactions among pairs without repulsive term
206!> \param pot ...
207!> \param rcov1 ...
208!> \param rcov2 ...
209!> \param z1 ...
210!> \param z2 ...
211!> \author i soliti ignoti
212! **************************************************************************************************
213 SUBROUTINE zbl_matching_polinomial(pot, rcov1, rcov2, z1, z2)
214
215 TYPE(pair_potential_single_type), POINTER :: pot
216 REAL(kind=dp), INTENT(IN) :: rcov1, rcov2, z1, z2
217
218 REAL(kind=dp) :: au, d1, d2, dd1, dd2, fac, v1, v2, x, &
219 x1, x2
220
221 pot%zbl_rcut(1) = (rcov1 + rcov2)*(1.0_dp - 0.2_dp)*bohr
222 pot%zbl_rcut(2) = (rcov1 + rcov2)*bohr
223 x1 = pot%zbl_rcut(1)
224 x2 = pot%zbl_rcut(2)
225 pot%z1 = z1
226 pot%z2 = z2
227
228 au = 0.88534_dp*bohr/(z1**0.23_dp + z2**0.23_dp)
229 x = x1/au
230 fac = z1*z2/evolt
231 v1 = fac/x1*(0.1818_dp*exp(-3.2_dp*x) + 0.5099_dp*exp(-0.9423_dp*x) + &
232 0.2802_dp*exp(-0.4029_dp*x) + 0.02817_dp*exp(-0.2016_dp*x))
233 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) &
234 - 0.4029_dp*0.2802_dp*exp(-0.4029_dp*x) - 0.2016_dp*0.02817_dp*exp(-0.2016_dp*x)) &
235 - fac/x1/x1*(0.1818_dp*exp(-3.2_dp*x) + 0.5099_dp*exp(-0.9423_dp*x) + &
236 0.2802_dp*exp(-0.4029_dp*x) + 0.02817_dp*exp(-0.2016_dp*x))
237
238 dd1 = 2.0_dp*fac/x1**3*(0.1818_dp*exp(-0.32e1_dp*x) &
239 + 0.5099_dp*exp(-0.9423_dp*x) + 0.2802_dp*exp(-0.4029_dp*x) &
240 + 0.2817e-1_dp*exp(-0.2016_dp*x)) &
241 - 0.2e1_dp*fac/x1**2/au*(-0.58176_dp*exp(-0.32e1_dp*x) - 0.48047877_dp*exp(-0.9423_dp*x) &
242 - 0.11289258_dp*exp(-0.4029_dp*x) - 0.5679072e-2_dp*exp(-0.2016_dp*x)) &
243 + fac/x1/au**2*(0.1861632e1_dp*exp(-0.32e1_dp*x) + &
244 0.4527551450_dp*exp(-0.9423_dp*x) + 0.4548442048e-1_dp*exp(-0.4029_dp*x) + &
245 0.1144900915e-2_dp*exp(-0.2016_dp*x))
246
247 v2 = 0.0_dp
248 d2 = 0.0_dp
249 dd2 = 0.0_dp
250
251 CALL compute_polinomial_5th(x1, v1, d1, dd1, x2, v2, d2, dd2, pot%zbl_poly)
252
253 END SUBROUTINE zbl_matching_polinomial
254
255! **************************************************************************************************
256!> \brief ...
257!> \param r1 ...
258!> \param v1 ...
259!> \param d1 ...
260!> \param dd1 ...
261!> \param r2 ...
262!> \param v2 ...
263!> \param d2 ...
264!> \param dd2 ...
265!> \param poly ...
266! **************************************************************************************************
267 SUBROUTINE compute_polinomial_5th(r1, v1, d1, dd1, r2, v2, d2, dd2, poly)
268
269 REAL(kind=dp) :: r1, v1, d1, dd1, r2, v2, d2, dd2, &
270 poly(0:5)
271
272 REAL(kind=dp) :: a0, a1, a2, a3, a4, a5
273
274! 5th order
275
276 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 - &
277 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 - &
278 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 + &
279 8._dp*r1**3*d1*r2**3 - 2._dp*r1**5*r2*d2 + r1**5*r2**2*dd2)/ &
280 (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)
281
282 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 + &
283 r2**4*r1**2*dd1 - 2._dp*r2**5*r1*dd1 - 10._dp*r2**4*r1*d1 + 10._dp*r1**4*r2*d2 - &
284 r1**4*r2**2*dd2 - 3._dp*r1**4*r2**2*dd1 + &
285 2._dp*r1**5*r2*dd2 - 24._dp*r2**3*r1**2*d2 - 16._dp*r2**3*r1**2*d1 + &
286 16._dp*r2**2*r1**3*d2 - 2._dp*r1**5*d2 + 2._dp*r2**5*d1 - &
287 60._dp*r1**2*r2**2*v1 + 60._dp*r1**2*r2**2*v2)/ &
288 (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)
289
290 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 - &
291 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 + &
292 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 - &
293 60._dp*r1**2*v1*r2 + r1**5*dd2 - r2**5*dd1)/ &
294 (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)
295
296 a3 = -.5_dp*(3._dp*r1**4*dd2 - r1**4*dd1 + 8.*r1**3*d1 - 4.*r1**3*r2*dd1 + &
297 12._dp*r1**3*d2 + 32._dp*r1**2*r2*d1 - 8._dp*r1**2*r2**2*dd2 - &
298 20._dp*r1**2*v1 + 8._dp*r1**2*r2**2*dd1 + 28._dp*r1**2*r2*d2 + &
299 20._dp*r1**2*v2 + 80._dp*r1*r2*v2 - 28._dp*r2**2*r1*d1 - 80._dp*r1*v1*r2 - &
300 32._dp*r2**2*r1*d2 + 4._dp*r1*r2**3*dd2 - 8._dp*r2**3*d2 - 12._dp*r2**3*d1 + &
301 r2**4*dd2 - 3._dp*r2**4*dd1 + 20._dp*r2**2*v2 - 20._dp*r2**2*v1)/ &
302 (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)
303
304 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 + &
305 16._dp*r1**2*d2 - 2._dp*r1*r2*d2 - r1*r2**2*dd2 - &
306 30._dp*r1*v1 + 30.*r1*v2 + 2._dp*r1*r2*d1 + 4._dp*r1*r2**2*dd1 - 16._dp*r2**2*d1 + &
307 2._dp*r2**3*dd2 - 14._dp*r2**2*d2 + 30._dp*r2*v2 - 30._dp*v1*r2 - &
308 3._dp*r2**3*dd1)/(10._dp*r2**2*r1**3 - 5._dp*r2*r1**4 - &
309 10._dp*r1**2*r2**3 + 5._dp*r2**4*r1 - r2**5 + r1**5)
310
311 a5 = -.5_dp*(6._dp*r1*d1 + 2._dp*r2*r1*dd1 + 6._dp*r1*d2 - 2.*r2*r1*dd2 - &
312 r2**2*dd1 - r1**2*dd1 - 12.*v1 + 12._dp*v2 + r1**2*dd2 - &
313 6._dp*r2*d1 + r2**2*dd2 - 6._dp*r2*d2)/ &
314 (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)
315
316 poly(0) = a0
317 poly(1) = a1
318 poly(2) = a2
319 poly(3) = a3
320 poly(4) = a4
321 poly(5) = a5
322
323 END SUBROUTINE compute_polinomial_5th
324
325END MODULE pair_potential_util
326
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:206
fparser::evalerrtype
integer, public evalerrtype
Definition fparser.F:33
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:36
pair_potential_types::bm_type
integer, parameter, public bm_type
Definition pair_potential_types.F:36
pair_potential_types::wl_type
integer, parameter, public wl_type
Definition pair_potential_types.F:36
pair_potential_types::ft_type
integer, parameter, public ft_type
Definition pair_potential_types.F:36
pair_potential_types::tab_type
integer, parameter, public tab_type
Definition pair_potential_types.F:36
pair_potential_types::ftd_type
integer, parameter, public ftd_type
Definition pair_potential_types.F:36
pair_potential_types::ip_type
integer, parameter, public ip_type
Definition pair_potential_types.F:36
pair_potential_types::lj_type
integer, parameter, public lj_type
Definition pair_potential_types.F:36
pair_potential_types::gp_type
integer, parameter, public gp_type
Definition pair_potential_types.F:36
pair_potential_types::gw_type
integer, parameter, public gw_type
Definition pair_potential_types.F:36
pair_potential_types::not_initialized
real(kind=dp), parameter, public not_initialized
Definition pair_potential_types.F:91
pair_potential_types::b4_type
integer, parameter, public b4_type
Definition pair_potential_types.F:36
pair_potential_types::ea_type
integer, parameter, public ea_type
Definition pair_potential_types.F:36
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:179
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:214
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:371