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 : ! **************************************************************************************************
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 779729834 : 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 779729834 : value = 0.0_dp
56 1559639036 : DO j = 1, SIZE(pot%type)
57 : ! A lower boundary for the potential definition was defined
58 779909202 : 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 779908544 : IF ((pot%set(j)%rmax /= not_initialized) .AND. (r >= pot%set(j)%rmax)) CYCLE
61 : ! If within limits let's compute the potential...
62 779906444 : 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 273074548 : sigma6*r**(-6))
66 : ELSE IF (pot%type(j) == lj_type) THEN
67 : lvalue = pot%set(j)%lj%epsilon* &
68 604304 : (pot%set(j)%lj%sigma12*r**(-12) - pot%set(j)%lj%sigma6*r**(-6))
69 : ELSE IF (pot%type(j) == ip_type) THEN
70 6057566 : lvalue = 0._dp
71 6057566 : IF (r > pot%set(j)%ipbv%rcore) THEN
72 66181110 : DO i = 2, 15
73 66181110 : 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 1645492 : 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 132970017 : 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 0 : (pot%set(j)%goodwin%d/pot%set(j)%goodwin%dc)**pot%set(j)%goodwin%mc))
85 0 : 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 25762 : 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 8994460 : bd6 = pot%set(j)%ftd%bd(1)
91 8994460 : dampsum = 1.0_dp
92 8994460 : xf = 1.0_dp
93 62961220 : DO i = 1, 6
94 53966760 : xf = xf*bd6*r
95 62961220 : dampsum = dampsum + xf*ifac(i)
96 : END DO
97 8994460 : f6 = 1.0_dp - EXP(-bd6*r)*dampsum
98 : ! Compute 8th order dispersion correction term
99 8994460 : bd8 = pot%set(j)%ftd%bd(2)
100 8994460 : dampsum = 1.0_dp
101 8994460 : xf = 1.0_dp
102 80950140 : DO i = 1, 8
103 71955680 : xf = xf*bd8*r
104 80950140 : dampsum = dampsum + xf*ifac(i)
105 : END DO
106 8994460 : f8 = 1.0_dp - EXP(-bd8*r)*dampsum
107 8994460 : 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 3854 : index = INT(r/pot%set(j)%eam%drar) + 1
110 3854 : IF (index > pot%set(j)%eam%npoints) THEN
111 : index = pot%set(j)%eam%npoints
112 3702 : ELSEIF (index < 1) THEN
113 : index = 1
114 : END IF
115 3854 : qq = r - pot%set(j)%eam%rval(index)
116 : pp = pot%set(j)%eam%phi(index) + &
117 3854 : qq*pot%set(j)%eam%phip(index)
118 3854 : lvalue = pp
119 : ELSE IF (pot%type(j) == b4_type) THEN
120 285055946 : IF (r <= pot%set(j)%buck4r%r1) THEN
121 73446226 : pp = pot%set(j)%buck4r%a*EXP(-pot%set(j)%buck4r%b*r)
122 211609720 : ELSEIF (r > pot%set(j)%buck4r%r1 .AND. r <= pot%set(j)%buck4r%r2) THEN
123 17904118 : pp = 0.0_dp
124 125328826 : DO n = 0, pot%set(j)%buck4r%npoly1
125 125328826 : pp = pp + pot%set(j)%buck4r%poly1(n)*r**n
126 : END DO
127 193705602 : ELSEIF (r > pot%set(j)%buck4r%r2 .AND. r <= pot%set(j)%buck4r%r3) THEN
128 19525932 : pp = 0.0_dp
129 97629660 : DO n = 0, pot%set(j)%buck4r%npoly2
130 97629660 : pp = pp + pot%set(j)%buck4r%poly2(n)*r**n
131 : END DO
132 174179670 : ELSEIF (r > pot%set(j)%buck4r%r3) THEN
133 174179670 : pp = -pot%set(j)%buck4r%c/r**6
134 : END IF
135 : lvalue = pp
136 : ELSE IF (pot%type(j) == tab_type) THEN
137 5647200 : index1 = FLOOR((r - pot%set(j)%tab%r(1))/pot%set(j)%tab%dr) + 1
138 5647200 : index2 = index1 + 1
139 5647200 : IF (index2 > pot%set(j)%tab%npoints) THEN
140 48 : index2 = pot%set(j)%tab%npoints
141 48 : index1 = index2 - 1
142 5647152 : ELSEIF (index1 < 1) THEN
143 882796 : index1 = 1
144 882796 : 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 5647200 : (pot%set(j)%tab%r(index2) - pot%set(j)%tab%r(index1))
149 5647200 : 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 1038716 : 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 50359896 : pot%set(j)%gp%values(1) = r
158 50359896 : lvalue = evalf(pot%set(j)%gp%myid, pot%set(j)%gp%values)
159 50359896 : IF (EvalErrType > 0) &
160 0 : CPABORT("Error evaluating generic potential energy function")
161 : ELSE
162 : lvalue = 0.0_dp
163 : END IF
164 1559639036 : value = value + lvalue
165 : END DO
166 779729834 : value = value - energy_cutoff
167 779729834 : 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 40236728 : 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 40236728 : ener_zbl = 0.0_dp
186 40236728 : IF (r <= pot%zbl_rcut(1)) THEN
187 13579482 : au = 0.88534_dp*bohr/(pot%z1**0.23_dp + pot%z2**0.23_dp)
188 13579482 : x = r/au
189 13579482 : 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 13579482 : 0.2802_dp*EXP(-0.4029_dp*x) + 0.02817_dp*EXP(-0.2016_dp*x))
192 26657246 : 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 1872516 : 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 40236728 : 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 48 : 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 48 : pot%zbl_rcut(1) = (rcov1 + rcov2)*(1.0_dp - 0.2_dp)*bohr
221 48 : pot%zbl_rcut(2) = (rcov1 + rcov2)*bohr
222 48 : x1 = pot%zbl_rcut(1)
223 48 : x2 = pot%zbl_rcut(2)
224 48 : pot%z1 = z1
225 48 : pot%z2 = z2
226 :
227 48 : au = 0.88534_dp*bohr/(z1**0.23_dp + z2**0.23_dp)
228 48 : x = x1/au
229 48 : 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 48 : 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 48 : 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 48 : 0.1144900915E-2_dp*EXP(-0.2016_dp*x))
245 :
246 48 : v2 = 0.0_dp
247 48 : d2 = 0.0_dp
248 48 : dd2 = 0.0_dp
249 :
250 48 : CALL compute_polinomial_5th(x1, v1, d1, dd1, x2, v2, d2, dd2, pot%zbl_poly)
251 :
252 48 : 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 48 : 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 48 : (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 48 : (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 48 : (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 48 : (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 48 : 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 48 : (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 48 : poly(0) = a0
316 48 : poly(1) = a1
317 48 : poly(2) = a2
318 48 : poly(3) = a3
319 48 : poly(4) = a4
320 48 : poly(5) = a5
321 :
322 48 : END SUBROUTINE compute_polinomial_5th
323 :
324 : END MODULE pair_potential_util
325 :
|
