Line data Source code
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 : !> Efficient tersoff implementation
11 : !> \author CJM, I-Feng W. Kuo, Teodoro Laino
12 : ! **************************************************************************************************
13 : MODULE manybody_tersoff
14 :
15 : USE cell_types, ONLY: cell_type
16 : USE fist_neighbor_list_types, ONLY: fist_neighbor_type,&
17 : neighbor_kind_pairs_type
18 : USE fist_nonbond_env_types, ONLY: pos_type
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: pi
21 : USE pair_potential_types, ONLY: pair_potential_pp_type,&
22 : pair_potential_single_type,&
23 : tersoff_pot_type,&
24 : tersoff_type
25 : USE util, ONLY: sort
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 :
30 : PRIVATE
31 : PUBLIC :: setup_tersoff_arrays, destroy_tersoff_arrays, &
32 : tersoff_forces, tersoff_energy
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_tersoff'
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief ...
39 : !> \param pot_loc ...
40 : !> \param tersoff ...
41 : !> \param r_last_update_pbc ...
42 : !> \param atom_a ...
43 : !> \param atom_b ...
44 : !> \param nloc_size ...
45 : !> \param full_loc_list ...
46 : !> \param loc_cell_v ...
47 : !> \param cell_v ...
48 : !> \param drij ...
49 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
50 : ! **************************************************************************************************
51 249248 : SUBROUTINE tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
52 249248 : full_loc_list, loc_cell_v, cell_v, drij)
53 :
54 : REAL(KIND=dp), INTENT(OUT) :: pot_loc
55 : TYPE(tersoff_pot_type), POINTER :: tersoff
56 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
57 : INTEGER, INTENT(IN) :: atom_a, atom_b, nloc_size
58 : INTEGER, DIMENSION(2, 1:nloc_size) :: full_loc_list
59 : REAL(KIND=dp), DIMENSION(3, 1:nloc_size) :: loc_cell_v
60 : REAL(KIND=dp), DIMENSION(3) :: cell_v
61 : REAL(KIND=dp) :: drij
62 :
63 : REAL(KIND=dp) :: b_ij, f_A, f_C, f_R
64 :
65 : b_ij = ter_b_ij(tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
66 249248 : full_loc_list, loc_cell_v, cell_v, tersoff%rcutsq)
67 249248 : f_C = ter_f_C(tersoff, drij)
68 249248 : f_A = ter_f_A(tersoff, drij)
69 249248 : f_R = ter_f_R(tersoff, drij)
70 249248 : pot_loc = f_C*(f_R + b_ij*f_A)
71 :
72 249248 : END SUBROUTINE tersoff_energy
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param tersoff ...
77 : !> \param r ...
78 : !> \return ...
79 : !> \author I-Feng W. Kuo
80 : ! **************************************************************************************************
81 4577776 : FUNCTION ter_f_C(tersoff, r)
82 : TYPE(tersoff_pot_type), POINTER :: tersoff
83 : REAL(KIND=dp), INTENT(IN) :: r
84 : REAL(KIND=dp) :: ter_f_C
85 :
86 : REAL(KIND=dp) :: bigD, bigR, RmD, RpD
87 :
88 4577776 : bigR = tersoff%bigR
89 4577776 : bigD = tersoff%bigD
90 4577776 : RmD = tersoff%bigR - tersoff%bigD
91 4577776 : RpD = tersoff%bigR + tersoff%bigD
92 4577776 : ter_f_C = 0.0_dp
93 4577776 : IF (r < RmD) ter_f_C = 1.0_dp
94 4577776 : IF (r > RpD) ter_f_C = 0.0_dp
95 4577776 : IF ((r < RpD) .AND. (r > RmD)) THEN
96 1762788 : ter_f_C = 0.5_dp*(1.0_dp - SIN(0.5_dp*PI*(r - bigR)/(bigD)))
97 : END IF
98 4577776 : END FUNCTION ter_f_C
99 :
100 : ! **************************************************************************************************
101 : !> \brief ...
102 : !> \param tersoff ...
103 : !> \param r ...
104 : !> \return ...
105 : !> \author I-Feng W. Kuo
106 : ! **************************************************************************************************
107 1269068 : FUNCTION ter_f_C_d(tersoff, r)
108 : TYPE(tersoff_pot_type), POINTER :: tersoff
109 : REAL(KIND=dp), INTENT(IN) :: r
110 : REAL(KIND=dp) :: ter_f_C_d
111 :
112 : REAL(KIND=dp) :: bigD, bigR, RmD, RpD
113 :
114 1269068 : bigR = tersoff%bigR
115 1269068 : bigD = tersoff%bigD
116 1269068 : RmD = tersoff%bigR - tersoff%bigD
117 1269068 : RpD = tersoff%bigR + tersoff%bigD
118 : ter_f_C_d = 0.0_dp
119 : IF (r < RmD) ter_f_C_d = 0.0_dp
120 : IF (r > RpD) ter_f_C_d = 0.0_dp
121 1269068 : IF ((r < RpD) .AND. (r > RmD)) THEN
122 466269 : ter_f_C_d = (0.25_dp*PI/bigD)*COS(0.5_dp*PI*(r - bigR)/(bigD))/r
123 : END IF
124 :
125 1269068 : END FUNCTION ter_f_C_d
126 :
127 : ! **************************************************************************************************
128 : !> \brief ...
129 : !> \param tersoff ...
130 : !> \param r ...
131 : !> \return ...
132 : !> \author I-Feng W. Kuo
133 : ! **************************************************************************************************
134 498496 : FUNCTION ter_f_R(tersoff, r)
135 : TYPE(tersoff_pot_type), POINTER :: tersoff
136 : REAL(KIND=dp), INTENT(IN) :: r
137 : REAL(KIND=dp) :: ter_f_R
138 :
139 : REAL(KIND=dp) :: A, lambda1
140 :
141 498496 : A = tersoff%A
142 498496 : lambda1 = tersoff%lambda1
143 498496 : ter_f_R = 0.0_dp
144 498496 : ter_f_R = A*EXP(-lambda1*r)
145 :
146 498496 : END FUNCTION ter_f_R
147 :
148 : ! **************************************************************************************************
149 : !> \brief ...
150 : !> \param tersoff ...
151 : !> \param r ...
152 : !> \return ...
153 : !> \author I-Feng W. Kuo
154 : ! **************************************************************************************************
155 249248 : FUNCTION ter_f_R_d(tersoff, r)
156 : TYPE(tersoff_pot_type), POINTER :: tersoff
157 : REAL(KIND=dp), INTENT(IN) :: r
158 : REAL(KIND=dp) :: ter_f_R_d
159 :
160 : REAL(KIND=dp) :: A, f_R, lambda1
161 :
162 249248 : A = tersoff%A
163 249248 : lambda1 = tersoff%lambda1
164 249248 : f_R = A*EXP(-lambda1*r)
165 249248 : ter_f_R_d = 0.0_dp
166 249248 : ter_f_R_d = lambda1*f_R/r
167 :
168 249248 : END FUNCTION ter_f_R_d
169 :
170 : ! **************************************************************************************************
171 : !> \brief ...
172 : !> \param tersoff ...
173 : !> \param r ...
174 : !> \return ...
175 : !> \author I-Feng W. Kuo
176 : ! **************************************************************************************************
177 498496 : FUNCTION ter_f_A(tersoff, r)
178 : TYPE(tersoff_pot_type), POINTER :: tersoff
179 : REAL(KIND=dp), INTENT(IN) :: r
180 : REAL(KIND=dp) :: ter_f_A
181 :
182 : REAL(KIND=dp) :: B, lambda2
183 :
184 498496 : B = tersoff%B
185 498496 : lambda2 = tersoff%lambda2
186 498496 : ter_f_A = 0.0_dp
187 498496 : ter_f_A = -B*EXP(-lambda2*r)
188 :
189 498496 : END FUNCTION ter_f_A
190 :
191 : ! **************************************************************************************************
192 : !> \brief ...
193 : !> \param tersoff ...
194 : !> \param r ...
195 : !> \return ...
196 : !> \author I-Feng W. Kuo
197 : ! **************************************************************************************************
198 249248 : FUNCTION ter_f_A_d(tersoff, r)
199 : TYPE(tersoff_pot_type), POINTER :: tersoff
200 : REAL(KIND=dp), INTENT(IN) :: r
201 : REAL(KIND=dp) :: ter_f_A_d
202 :
203 : REAL(KIND=dp) :: B, lambda2
204 :
205 249248 : B = tersoff%B
206 249248 : lambda2 = tersoff%lambda2
207 249248 : ter_f_A_d = 0.0_dp
208 249248 : ter_f_A_d = -B*lambda2*EXP(-lambda2*r)/r
209 :
210 249248 : END FUNCTION ter_f_A_d
211 :
212 : ! **************************************************************************************************
213 : !> \brief ...
214 : !> \param tersoff ...
215 : !> \return ...
216 : !> \author I-Feng W. Kuo
217 : ! **************************************************************************************************
218 0 : FUNCTION ter_a_ij(tersoff)
219 : TYPE(tersoff_pot_type), POINTER :: tersoff
220 : REAL(KIND=dp) :: ter_a_ij
221 :
222 : REAL(KIND=dp) :: alpha, n
223 :
224 0 : n = tersoff%n
225 0 : alpha = tersoff%alpha
226 : ter_a_ij = 0.0_dp
227 : !Note alpha = 0.0_dp for the parameters in the paper so using simplified term
228 : !ter_a_ij = (1.0_dp+(alpha*ter_n_ij(tersoff,iparticle,jparticle,r))**n)**(-0.5_dp/n)
229 0 : ter_a_ij = 1.0_dp
230 :
231 0 : END FUNCTION ter_a_ij
232 :
233 : ! **************************************************************************************************
234 : !> \brief ...
235 : !> \param tersoff ...
236 : !> \param r_last_update_pbc ...
237 : !> \param iparticle ...
238 : !> \param jparticle ...
239 : !> \param n_loc_size ...
240 : !> \param full_loc_list ...
241 : !> \param loc_cell_v ...
242 : !> \param cell_v ...
243 : !> \param rcutsq ...
244 : !> \return ...
245 : !> \author I-Feng W. Kuo, Teodoro Laino
246 : ! **************************************************************************************************
247 498496 : FUNCTION ter_b_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
248 498496 : full_loc_list, loc_cell_v, cell_v, rcutsq)
249 : TYPE(tersoff_pot_type), POINTER :: tersoff
250 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
251 : INTEGER, INTENT(IN) :: iparticle, jparticle, n_loc_size
252 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
253 : REAL(KIND=dp), DIMENSION(3, 1:n_loc_size) :: loc_cell_v
254 : REAL(KIND=dp), DIMENSION(3) :: cell_v
255 : REAL(KIND=dp), INTENT(IN) :: rcutsq
256 : REAL(KIND=dp) :: ter_b_ij
257 :
258 : REAL(KIND=dp) :: beta, n, zeta_ij
259 :
260 498496 : n = tersoff%n
261 498496 : beta = tersoff%beta
262 498496 : ter_b_ij = 0.0_dp
263 : zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, &
264 498496 : n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq)
265 498496 : ter_b_ij = (1.0_dp + (beta*zeta_ij)**n)**(-0.5_dp/n)
266 :
267 498496 : END FUNCTION ter_b_ij
268 :
269 : ! **************************************************************************************************
270 : !> \brief ...
271 : !> \param tersoff ...
272 : !> \param r_last_update_pbc ...
273 : !> \param iparticle ...
274 : !> \param jparticle ...
275 : !> \param n_loc_size ...
276 : !> \param full_loc_list ...
277 : !> \param loc_cell_v ...
278 : !> \param cell_v ...
279 : !> \param rcutsq ...
280 : !> \return ...
281 : !> \author I-Feng W. Kuo, Teodoro Laino
282 : ! **************************************************************************************************
283 249248 : FUNCTION ter_b_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
284 249248 : full_loc_list, loc_cell_v, cell_v, rcutsq)
285 : TYPE(tersoff_pot_type), POINTER :: tersoff
286 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
287 : INTEGER, INTENT(IN) :: iparticle, jparticle, n_loc_size
288 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
289 : REAL(KIND=dp), DIMENSION(3, 1:n_loc_size) :: loc_cell_v
290 : REAL(KIND=dp), DIMENSION(3) :: cell_v
291 : REAL(KIND=dp), INTENT(IN) :: rcutsq
292 : REAL(KIND=dp) :: ter_b_ij_d
293 :
294 : REAL(KIND=dp) :: beta, beta_n, n, zeta_ij, zeta_ij_n, &
295 : zeta_ij_nm1
296 :
297 249248 : n = tersoff%n
298 249248 : beta = tersoff%beta
299 249248 : beta_n = beta**n
300 : zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
301 249248 : full_loc_list, loc_cell_v, cell_v, rcutsq)
302 249248 : zeta_ij_nm1 = 0.0_dp
303 249248 : IF (zeta_ij > 0.0_dp) zeta_ij_nm1 = zeta_ij**(n - 1.0_dp)
304 249248 : zeta_ij_n = zeta_ij**(n)
305 :
306 249248 : ter_b_ij_d = 0.0_dp
307 : ter_b_ij_d = -0.5_dp*beta_n*zeta_ij_nm1* &
308 249248 : ((1.0_dp + beta_n*zeta_ij_n)**((-0.5_dp/n) - 1.0_dp))
309 :
310 249248 : END FUNCTION ter_b_ij_d
311 :
312 : ! **************************************************************************************************
313 : !> \brief ...
314 : !> \param tersoff ...
315 : !> \param r_last_update_pbc ...
316 : !> \param iparticle ...
317 : !> \param jparticle ...
318 : !> \param n_loc_size ...
319 : !> \param full_loc_list ...
320 : !> \param loc_cell_v ...
321 : !> \param cell_v ...
322 : !> \param rcutsq ...
323 : !> \return ...
324 : !> \par History
325 : !> Using a local list of neighbors - [tlaino] 2007
326 : !> \author I-Feng W. Kuo, Teodoro Laino
327 : ! **************************************************************************************************
328 747744 : FUNCTION ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
329 747744 : full_loc_list, loc_cell_v, cell_v, rcutsq)
330 : TYPE(tersoff_pot_type), POINTER :: tersoff
331 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
332 : INTEGER, INTENT(IN) :: iparticle, jparticle, n_loc_size
333 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
334 : REAL(KIND=dp), DIMENSION(3, 1:n_loc_size) :: loc_cell_v
335 : REAL(KIND=dp), DIMENSION(3) :: cell_v
336 : REAL(KIND=dp), INTENT(IN) :: rcutsq
337 : REAL(KIND=dp) :: ter_zeta_ij
338 :
339 : INTEGER :: ilist, kparticle
340 : REAL(KIND=dp) :: cell_v_2(3), costheta, drij, drik, &
341 : expterm, f_C, gterm, lambda3, n, &
342 : rab2_max, rij(3), rik(3)
343 :
344 747744 : ter_zeta_ij = 0.0_dp
345 747744 : n = tersoff%n
346 747744 : lambda3 = tersoff%lambda3
347 747744 : rab2_max = rcutsq
348 2990976 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
349 2990976 : drij = NORM2(rij)
350 747744 : ter_zeta_ij = 0.0_dp
351 80395608 : DO ilist = 1, n_loc_size
352 79647864 : kparticle = full_loc_list(2, ilist)
353 79647864 : IF (kparticle == jparticle) CYCLE
354 312317016 : cell_v_2 = loc_cell_v(:, ilist)
355 312317016 : rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
356 312317016 : drik = DOT_PRODUCT(rik, rik)
357 78079254 : IF (drik > rab2_max) CYCLE
358 3059460 : drik = SQRT(drik)
359 12237840 : costheta = DOT_PRODUCT(rij, rik)/(drij*drik)
360 3059460 : IF (costheta < -1.0_dp) costheta = -1.0_dp
361 3059460 : IF (costheta > +1.0_dp) costheta = +1.0_dp
362 3059460 : f_C = ter_f_C(tersoff, drik)
363 3059460 : gterm = ter_g(tersoff, costheta)
364 3059460 : expterm = EXP((lambda3*(drij - drik))**3)
365 80395608 : ter_zeta_ij = ter_zeta_ij + f_C*gterm*expterm
366 : END DO
367 :
368 747744 : END FUNCTION ter_zeta_ij
369 :
370 : ! **************************************************************************************************
371 : !> \brief ...
372 : !> \param tersoff ...
373 : !> \param r_last_update_pbc ...
374 : !> \param iparticle ...
375 : !> \param jparticle ...
376 : !> \param f_nonbond ...
377 : !> \param pv_nonbond ...
378 : !> \param prefactor ...
379 : !> \param n_loc_size ...
380 : !> \param full_loc_list ...
381 : !> \param loc_cell_v ...
382 : !> \param cell_v ...
383 : !> \param rcutsq ...
384 : !> \param use_virial ...
385 : !> \par History
386 : !> Using a local list of neighbors - [tlaino] 2007
387 : !> \author I-Feng W. Kuo, Teodoro Laino
388 : ! **************************************************************************************************
389 249248 : SUBROUTINE ter_zeta_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, prefactor, &
390 249248 : n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq, use_virial)
391 : TYPE(tersoff_pot_type), POINTER :: tersoff
392 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
393 : INTEGER, INTENT(IN) :: iparticle, jparticle
394 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
395 : REAL(KIND=dp), INTENT(IN) :: prefactor
396 : INTEGER, INTENT(IN) :: n_loc_size
397 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
398 : REAL(KIND=dp), DIMENSION(3, 1:n_loc_size) :: loc_cell_v
399 : REAL(KIND=dp), DIMENSION(3) :: cell_v
400 : REAL(KIND=dp), INTENT(IN) :: rcutsq
401 : LOGICAL, INTENT(IN) :: use_virial
402 :
403 : INTEGER :: ilist, kparticle, nparticle
404 : REAL(KIND=dp) :: costheta, drij, drik, expterm, &
405 : expterm_d, f_C, f_C_d, gterm, gterm_d, &
406 : lambda3, n, rab2_max
407 : REAL(KIND=dp), DIMENSION(3) :: cell_v_2, dcosdri, dcosdrj, dcosdrk, &
408 : dri, drj, drk, rij, rij_hat, rik, &
409 : rik_hat
410 :
411 249248 : n = tersoff%n
412 249248 : lambda3 = tersoff%lambda3
413 249248 : rab2_max = rcutsq
414 :
415 996992 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
416 996992 : drij = NORM2(rij)
417 996992 : rij_hat(:) = rij(:)/drij
418 :
419 26798536 : nparticle = SIZE(r_last_update_pbc)
420 26798536 : DO ilist = 1, n_loc_size
421 26549288 : kparticle = full_loc_list(2, ilist)
422 26549288 : IF (kparticle == jparticle) CYCLE
423 104105672 : cell_v_2 = loc_cell_v(:, ilist)
424 104105672 : rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
425 104105672 : drik = DOT_PRODUCT(rik, rik)
426 :
427 26026418 : IF (drik > rab2_max) CYCLE
428 1019820 : drik = SQRT(drik)
429 4079280 : rik_hat(:) = rik(:)/drik
430 4079280 : costheta = DOT_PRODUCT(rij, rik)/(drij*drik)
431 1019820 : IF (costheta < -1.0_dp) costheta = -1.0_dp
432 1019820 : IF (costheta > +1.0_dp) costheta = +1.0_dp
433 :
434 4079280 : dcosdrj(:) = (1.0_dp/(drij))*(rik_hat(:) - costheta*rij_hat(:))
435 4079280 : dcosdrk(:) = (1.0_dp/(drik))*(rij_hat(:) - costheta*rik_hat(:))
436 4079280 : dcosdri(:) = -(dcosdrj(:) + dcosdrk(:))
437 :
438 1019820 : f_C = ter_f_C(tersoff, drik)
439 1019820 : f_C_d = ter_f_C_d(tersoff, drik)
440 1019820 : gterm = ter_g(tersoff, costheta)
441 1019820 : gterm_d = ter_g_d(tersoff, costheta) !still need d(costheta)/dR term
442 1019820 : expterm = EXP((lambda3*(drij - drik))**3)
443 1019820 : expterm_d = (3.0_dp)*(lambda3**3)*((drij - drik)**2)*expterm
444 :
445 : dri = f_C_d*gterm*expterm*(rik) &
446 : + f_C*gterm_d*expterm*(dcosdri) &
447 4079280 : + f_C*gterm*expterm_d*(-rij_hat + rik_hat)
448 :
449 : !No f_C_d component for Rj
450 : drj = f_C*gterm_d*expterm*(dcosdrj) &
451 4079280 : + f_C*gterm*expterm_d*(rij_hat)
452 :
453 : drk = f_C_d*gterm*expterm*(-rik) &
454 : + f_C*gterm_d*expterm*(dcosdrk) &
455 4079280 : + f_C*gterm*expterm_d*(-rik_hat)
456 :
457 1019820 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
458 1019820 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
459 1019820 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
460 :
461 1019820 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
462 1019820 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
463 1019820 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
464 :
465 1019820 : f_nonbond(1, kparticle) = f_nonbond(1, kparticle) + prefactor*drk(1)
466 1019820 : f_nonbond(2, kparticle) = f_nonbond(2, kparticle) + prefactor*drk(2)
467 1019820 : f_nonbond(3, kparticle) = f_nonbond(3, kparticle) + prefactor*drk(3)
468 :
469 1269068 : IF (use_virial) THEN
470 279274 : pv_nonbond(1, 1) = pv_nonbond(1, 1) + prefactor*(rij(1)*drj(1) + rik(1)*drk(1))
471 279274 : pv_nonbond(1, 2) = pv_nonbond(1, 2) + prefactor*(rij(1)*drj(2) + rik(1)*drk(2))
472 279274 : pv_nonbond(1, 3) = pv_nonbond(1, 3) + prefactor*(rij(1)*drj(3) + rik(1)*drk(3))
473 :
474 279274 : pv_nonbond(2, 1) = pv_nonbond(2, 1) + prefactor*(rij(2)*drj(1) + rik(2)*drk(1))
475 279274 : pv_nonbond(2, 2) = pv_nonbond(2, 2) + prefactor*(rij(2)*drj(2) + rik(2)*drk(2))
476 279274 : pv_nonbond(2, 3) = pv_nonbond(2, 3) + prefactor*(rij(2)*drj(3) + rik(2)*drk(3))
477 :
478 279274 : pv_nonbond(3, 1) = pv_nonbond(3, 1) + prefactor*(rij(3)*drj(1) + rik(3)*drk(1))
479 279274 : pv_nonbond(3, 2) = pv_nonbond(3, 2) + prefactor*(rij(3)*drj(2) + rik(3)*drk(2))
480 279274 : pv_nonbond(3, 3) = pv_nonbond(3, 3) + prefactor*(rij(3)*drj(3) + rik(3)*drk(3))
481 : END IF
482 : END DO
483 249248 : END SUBROUTINE ter_zeta_ij_d
484 :
485 : ! **************************************************************************************************
486 : !> \brief ...
487 : !> \param tersoff ...
488 : !> \param costheta ...
489 : !> \return ...
490 : !> \author I-Feng W. Kuo
491 : ! **************************************************************************************************
492 4079280 : FUNCTION ter_g(tersoff, costheta)
493 : TYPE(tersoff_pot_type), POINTER :: tersoff
494 : REAL(KIND=dp), INTENT(IN) :: costheta
495 : REAL(KIND=dp) :: ter_g
496 :
497 : REAL(KIND=dp) :: c, c2, d, d2, h
498 :
499 4079280 : c = tersoff%c
500 4079280 : d = tersoff%d
501 4079280 : h = tersoff%h
502 4079280 : c2 = c*c
503 4079280 : d2 = d*d
504 4079280 : ter_g = 0.0_dp
505 4079280 : ter_g = 1.0_dp + (c2/d2) - (c2)/(d2 + (h - costheta)**2)
506 :
507 4079280 : END FUNCTION ter_g
508 :
509 : ! **************************************************************************************************
510 : !> \brief ...
511 : !> \param tersoff ...
512 : !> \param costheta ...
513 : !> \return ...
514 : !> \author I-Feng W. Kuo
515 : ! **************************************************************************************************
516 1019820 : FUNCTION ter_g_d(tersoff, costheta)
517 : TYPE(tersoff_pot_type), POINTER :: tersoff
518 : REAL(KIND=dp), INTENT(IN) :: costheta
519 : REAL(KIND=dp) :: ter_g_d
520 :
521 : REAL(KIND=dp) :: c, c2, d, d2, h, hc
522 :
523 1019820 : c = tersoff%c
524 1019820 : d = tersoff%d
525 1019820 : h = tersoff%h
526 1019820 : c2 = c*c
527 1019820 : d2 = d*d
528 1019820 : hc = h - costheta
529 :
530 1019820 : ter_g_d = 0.0_dp
531 : ! Still need d(costheta)/dR
532 1019820 : ter_g_d = (-2.0_dp*c2*hc)/(d2 + hc**2)**2
533 1019820 : END FUNCTION ter_g_d
534 :
535 : ! **************************************************************************************************
536 : !> \brief ...
537 : !> \param tersoff ...
538 : !> \param r_last_update_pbc ...
539 : !> \param cell_v ...
540 : !> \param n_loc_size ...
541 : !> \param full_loc_list ...
542 : !> \param loc_cell_v ...
543 : !> \param iparticle ...
544 : !> \param jparticle ...
545 : !> \param f_nonbond ...
546 : !> \param pv_nonbond ...
547 : !> \param use_virial ...
548 : !> \param rcutsq ...
549 : !> \par History
550 : !> Using a local list of neighbors - [tlaino] 2007
551 : !> \author I-Feng W. Kuo, Teodoro Laino
552 : ! **************************************************************************************************
553 498496 : SUBROUTINE tersoff_forces(tersoff, r_last_update_pbc, cell_v, n_loc_size, &
554 249248 : full_loc_list, loc_cell_v, iparticle, jparticle, f_nonbond, pv_nonbond, &
555 : use_virial, rcutsq)
556 : TYPE(tersoff_pot_type), POINTER :: tersoff
557 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
558 : REAL(KIND=dp), DIMENSION(3) :: cell_v
559 : INTEGER, INTENT(IN) :: n_loc_size
560 : INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
561 : REAL(KIND=dp), DIMENSION(3, 1:n_loc_size) :: loc_cell_v
562 : INTEGER, INTENT(IN) :: iparticle, jparticle
563 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
564 : LOGICAL, INTENT(IN) :: use_virial
565 : REAL(KIND=dp), INTENT(IN) :: rcutsq
566 :
567 : CHARACTER(LEN=*), PARAMETER :: routineN = 'tersoff_forces'
568 :
569 : INTEGER :: handle
570 : REAL(KIND=dp) :: b_ij, b_ij_d, drij, f_A, f_A1, f_A2, &
571 : f_A_d, f_C, f_C_d, f_R, f_R1, f_R2, &
572 : f_R_d, fac, prefactor, rij(3), &
573 : rij_hat(3)
574 :
575 249248 : CALL timeset(routineN, handle)
576 996992 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
577 996992 : drij = NORM2(rij)
578 996992 : rij_hat(:) = rij(:)/drij
579 :
580 249248 : fac = -0.5_dp
581 249248 : b_ij = ter_b_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq)
582 249248 : b_ij_d = ter_b_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq)
583 249248 : f_A = ter_f_A(tersoff, drij)
584 249248 : f_A_d = ter_f_A_d(tersoff, drij)
585 249248 : f_C = ter_f_C(tersoff, drij)
586 249248 : f_C_d = ter_f_C_d(tersoff, drij)
587 249248 : f_R = ter_f_R(tersoff, drij)
588 249248 : f_R_d = ter_f_R_d(tersoff, drij)
589 :
590 : ! Lets do the easy one first, the repulsive term
591 : ! Note a_ij = 1.0_dp so just going to ignore it...
592 249248 : f_R1 = f_C_d*f_R*fac
593 249248 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_R1*rij(1)
594 249248 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_R1*rij(2)
595 249248 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_R1*rij(3)
596 249248 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_R1*rij(1)
597 249248 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_R1*rij(2)
598 249248 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_R1*rij(3)
599 :
600 249248 : IF (use_virial) THEN
601 73450 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_R1*rij(1)*rij(1)
602 73450 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_R1*rij(1)*rij(2)
603 73450 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_R1*rij(1)*rij(3)
604 73450 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_R1*rij(2)*rij(1)
605 73450 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_R1*rij(2)*rij(2)
606 73450 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_R1*rij(2)*rij(3)
607 73450 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_R1*rij(3)*rij(1)
608 73450 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_R1*rij(3)*rij(2)
609 73450 : pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_R1*rij(3)*rij(3)
610 : END IF
611 :
612 249248 : f_R2 = f_C*f_R_d*fac
613 249248 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_R2*rij(1)
614 249248 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_R2*rij(2)
615 249248 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_R2*rij(3)
616 249248 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_R2*rij(1)
617 249248 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_R2*rij(2)
618 249248 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_R2*rij(3)
619 :
620 249248 : IF (use_virial) THEN
621 73450 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_R2*rij(1)*rij(1)
622 73450 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_R2*rij(1)*rij(2)
623 73450 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_R2*rij(1)*rij(3)
624 73450 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_R2*rij(2)*rij(1)
625 73450 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_R2*rij(2)*rij(2)
626 73450 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_R2*rij(2)*rij(3)
627 73450 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_R2*rij(3)*rij(1)
628 73450 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_R2*rij(3)*rij(2)
629 73450 : pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_R2*rij(3)*rij(3)
630 : END IF
631 :
632 : ! Lets do the f_A1 piece derivative of F_C
633 249248 : f_A1 = f_C_d*b_ij*f_A*fac
634 249248 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A1*rij(1)
635 249248 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A1*rij(2)
636 249248 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A1*rij(3)
637 249248 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A1*rij(1)
638 249248 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A1*rij(2)
639 249248 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A1*rij(3)
640 :
641 249248 : IF (use_virial) THEN
642 73450 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_A1*rij(1)*rij(1)
643 73450 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_A1*rij(1)*rij(2)
644 73450 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_A1*rij(1)*rij(3)
645 73450 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_A1*rij(2)*rij(1)
646 73450 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_A1*rij(2)*rij(2)
647 73450 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_A1*rij(2)*rij(3)
648 73450 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_A1*rij(3)*rij(1)
649 73450 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_A1*rij(3)*rij(2)
650 73450 : pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_A1*rij(3)*rij(3)
651 : END IF
652 :
653 : ! Lets do the f_A2 piece derivative of F_A
654 249248 : f_A2 = f_C*b_ij*f_A_d*fac
655 249248 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A2*rij(1)
656 249248 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A2*rij(2)
657 249248 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A2*rij(3)
658 249248 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A2*rij(1)
659 249248 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A2*rij(2)
660 249248 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A2*rij(3)
661 :
662 249248 : IF (use_virial) THEN
663 73450 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_A2*rij(1)*rij(1)
664 73450 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_A2*rij(1)*rij(2)
665 73450 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_A2*rij(1)*rij(3)
666 73450 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_A2*rij(2)*rij(1)
667 73450 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_A2*rij(2)*rij(2)
668 73450 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_A2*rij(2)*rij(3)
669 73450 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_A2*rij(3)*rij(1)
670 73450 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_A2*rij(3)*rij(2)
671 73450 : pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_A2*rij(3)*rij(3)
672 : END IF
673 :
674 : ! Lets do the f_A3 piece derivative of b_ij
675 249248 : prefactor = f_C*b_ij_d*f_A*fac ! Note need to do d(Zeta_ij)/dR
676 : CALL ter_zeta_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, prefactor, &
677 249248 : n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq, use_virial)
678 249248 : CALL timestop(handle)
679 249248 : END SUBROUTINE tersoff_forces
680 :
681 : ! **************************************************************************************************
682 : !> \brief ...
683 : !> \param nonbonded ...
684 : !> \param potparm ...
685 : !> \param glob_loc_list ...
686 : !> \param glob_cell_v ...
687 : !> \param glob_loc_list_a ...
688 : !> \param cell ...
689 : !> \par History
690 : !> Fast implementation of the tersoff potential - [tlaino] 2007
691 : !> \author Teodoro Laino - University of Zurich
692 : ! **************************************************************************************************
693 5328 : SUBROUTINE setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
694 : TYPE(fist_neighbor_type), POINTER :: nonbonded
695 : TYPE(pair_potential_pp_type), POINTER :: potparm
696 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
697 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
698 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
699 : TYPE(cell_type), POINTER :: cell
700 :
701 : CHARACTER(LEN=*), PARAMETER :: routineN = 'setup_tersoff_arrays'
702 :
703 : INTEGER :: handle, i, iend, igrp, ikind, ilist, &
704 : ipair, istart, jkind, nkinds, npairs, &
705 : npairs_tot
706 5328 : INTEGER, DIMENSION(:), POINTER :: work_list, work_list2
707 5328 : INTEGER, DIMENSION(:, :), POINTER :: list
708 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi
709 5328 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rwork_list
710 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
711 : TYPE(pair_potential_single_type), POINTER :: pot
712 :
713 0 : CPASSERT(.NOT. ASSOCIATED(glob_loc_list))
714 5328 : CPASSERT(.NOT. ASSOCIATED(glob_loc_list_a))
715 5328 : CPASSERT(.NOT. ASSOCIATED(glob_cell_v))
716 5328 : CALL timeset(routineN, handle)
717 5328 : npairs_tot = 0
718 5328 : nkinds = SIZE(potparm%pot, 1)
719 202104 : DO ilist = 1, nonbonded%nlists
720 196776 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
721 196776 : npairs = neighbor_kind_pair%npairs
722 196776 : IF (npairs == 0) CYCLE
723 136074 : Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
724 66494 : istart = neighbor_kind_pair%grp_kind_start(igrp)
725 66494 : iend = neighbor_kind_pair%grp_kind_end(igrp)
726 66494 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
727 66494 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
728 66494 : pot => potparm%pot(ikind, jkind)%pot
729 66494 : npairs = iend - istart + 1
730 66494 : IF (pot%no_mb) CYCLE Kind_Group_Loop1
731 329624 : DO i = 1, SIZE(pot%type)
732 132920 : IF (pot%type(i) == tersoff_type) npairs_tot = npairs_tot + npairs
733 : END DO
734 : END DO Kind_Group_Loop1
735 : END DO
736 15984 : ALLOCATE (work_list(npairs_tot))
737 10656 : ALLOCATE (work_list2(npairs_tot))
738 15984 : ALLOCATE (glob_loc_list(2, npairs_tot))
739 15984 : ALLOCATE (glob_cell_v(3, npairs_tot))
740 : ! Fill arrays with data
741 5328 : npairs_tot = 0
742 202104 : DO ilist = 1, nonbonded%nlists
743 196776 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
744 196776 : npairs = neighbor_kind_pair%npairs
745 196776 : IF (npairs == 0) CYCLE
746 136074 : Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
747 66494 : istart = neighbor_kind_pair%grp_kind_start(igrp)
748 66494 : iend = neighbor_kind_pair%grp_kind_end(igrp)
749 66494 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
750 66494 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
751 66494 : list => neighbor_kind_pair%list
752 265976 : cvi = neighbor_kind_pair%cell_vector
753 66494 : pot => potparm%pot(ikind, jkind)%pot
754 66494 : npairs = iend - istart + 1
755 66494 : IF (pot%no_mb) CYCLE Kind_Group_Loop2
756 863486 : cell_v = MATMUL(cell%hmat, cvi)
757 329624 : DO i = 1, SIZE(pot%type)
758 : ! TERSOFF
759 132920 : IF (pot%type(i) == tersoff_type) THEN
760 9553152 : DO ipair = 1, npairs
761 56920440 : glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
762 38013372 : glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
763 : END DO
764 66412 : npairs_tot = npairs_tot + npairs
765 : END IF
766 : END DO
767 : END DO Kind_Group_Loop2
768 : END DO
769 : ! Order the arrays w.r.t. the first index of glob_loc_list
770 5328 : CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
771 9492068 : DO ipair = 1, npairs_tot
772 9492068 : work_list2(ipair) = glob_loc_list(2, work_list(ipair))
773 : END DO
774 18984136 : glob_loc_list(2, :) = work_list2
775 5328 : DEALLOCATE (work_list2)
776 15984 : ALLOCATE (rwork_list(3, npairs_tot))
777 9492068 : DO ipair = 1, npairs_tot
778 75899248 : rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
779 : END DO
780 75904576 : glob_cell_v = rwork_list
781 5328 : DEALLOCATE (rwork_list)
782 5328 : DEALLOCATE (work_list)
783 15984 : ALLOCATE (glob_loc_list_a(npairs_tot))
784 18984136 : glob_loc_list_a = glob_loc_list(1, :)
785 5328 : CALL timestop(handle)
786 10656 : END SUBROUTINE setup_tersoff_arrays
787 :
788 : ! **************************************************************************************************
789 : !> \brief ...
790 : !> \param glob_loc_list ...
791 : !> \param glob_cell_v ...
792 : !> \param glob_loc_list_a ...
793 : !> \par History
794 : !> Fast implementation of the tersoff potential - [tlaino] 2007
795 : !> \author Teodoro Laino - University of Zurich
796 : ! **************************************************************************************************
797 5328 : SUBROUTINE destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
798 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
799 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
800 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
801 :
802 5328 : IF (ASSOCIATED(glob_loc_list)) THEN
803 5328 : DEALLOCATE (glob_loc_list)
804 : END IF
805 5328 : IF (ASSOCIATED(glob_loc_list_a)) THEN
806 5328 : DEALLOCATE (glob_loc_list_a)
807 : END IF
808 5328 : IF (ASSOCIATED(glob_cell_v)) THEN
809 5328 : DEALLOCATE (glob_cell_v)
810 : END IF
811 :
812 5328 : END SUBROUTINE destroy_tersoff_arrays
813 :
814 : END MODULE manybody_tersoff
815 :
|
