(git:374b731)
Loading...
Searching...
No Matches
manybody_tersoff.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!> Efficient tersoff implementation
11!> \author CJM, I-Feng W. Kuo, Teodoro Laino
12! **************************************************************************************************
13MODULE 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
35CONTAINS
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 SUBROUTINE tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
52 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 full_loc_list, loc_cell_v, cell_v, tersoff%rcutsq)
67 f_c = ter_f_c(tersoff, drij)
68 f_a = ter_f_a(tersoff, drij)
69 f_r = ter_f_r(tersoff, drij)
70 pot_loc = f_c*(f_r + b_ij*f_a)
71
72 END SUBROUTINE tersoff_energy
73
74! **************************************************************************************************
75!> \brief ...
76!> \param tersoff ...
77!> \param r ...
78!> \return ...
79!> \author I-Feng W. Kuo
80! **************************************************************************************************
81 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 bigr = tersoff%bigR
89 bigd = tersoff%bigD
90 rmd = tersoff%bigR - tersoff%bigD
91 rpd = tersoff%bigR + tersoff%bigD
92 ter_f_c = 0.0_dp
93 IF (r < rmd) ter_f_c = 1.0_dp
94 IF (r > rpd) ter_f_c = 0.0_dp
95 IF ((r < rpd) .AND. (r > rmd)) THEN
96 ter_f_c = 0.5_dp*(1.0_dp - sin(0.5_dp*pi*(r - bigr)/(bigd)))
97 END IF
98 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 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 bigr = tersoff%bigR
115 bigd = tersoff%bigD
116 rmd = tersoff%bigR - tersoff%bigD
117 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 IF ((r < rpd) .AND. (r > rmd)) THEN
122 ter_f_c_d = (0.25_dp*pi/bigd)*cos(0.5_dp*pi*(r - bigr)/(bigd))/r
123 END IF
124
125 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 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 a = tersoff%A
142 lambda1 = tersoff%lambda1
143 ter_f_r = 0.0_dp
144 ter_f_r = a*exp(-lambda1*r)
145
146 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 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 a = tersoff%A
163 lambda1 = tersoff%lambda1
164 f_r = a*exp(-lambda1*r)
165 ter_f_r_d = 0.0_dp
166 ter_f_r_d = lambda1*f_r/r
167
168 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 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 b = tersoff%B
185 lambda2 = tersoff%lambda2
186 ter_f_a = 0.0_dp
187 ter_f_a = -b*exp(-lambda2*r)
188
189 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 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 b = tersoff%B
206 lambda2 = tersoff%lambda2
207 ter_f_a_d = 0.0_dp
208 ter_f_a_d = -b*lambda2*exp(-lambda2*r)/r
209
210 END FUNCTION ter_f_a_d
211
212! **************************************************************************************************
213!> \brief ...
214!> \param tersoff ...
215!> \return ...
216!> \author I-Feng W. Kuo
217! **************************************************************************************************
218 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 n = tersoff%n
225 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 ter_a_ij = 1.0_dp
230
231 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 FUNCTION ter_b_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
248 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 n = tersoff%n
261 beta = tersoff%beta
262 ter_b_ij = 0.0_dp
263 zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, &
264 n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq)
265 ter_b_ij = (1.0_dp + (beta*zeta_ij)**n)**(-0.5_dp/n)
266
267 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 FUNCTION ter_b_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
284 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 n = tersoff%n
298 beta = tersoff%beta
299 beta_n = beta**n
300 zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
301 full_loc_list, loc_cell_v, cell_v, rcutsq)
302 zeta_ij_nm1 = 0.0_dp
303 IF (zeta_ij > 0.0_dp) zeta_ij_nm1 = zeta_ij**(n - 1.0_dp)
304 zeta_ij_n = zeta_ij**(n)
305
306 ter_b_ij_d = 0.0_dp
307 ter_b_ij_d = -0.5_dp*beta_n*zeta_ij_nm1* &
308 ((1.0_dp + beta_n*zeta_ij_n)**((-0.5_dp/n) - 1.0_dp))
309
310 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 FUNCTION ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
329 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 ter_zeta_ij = 0.0_dp
345 n = tersoff%n
346 lambda3 = tersoff%lambda3
347 rab2_max = rcutsq
348 rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
349 drij = sqrt(dot_product(rij, rij))
350 ter_zeta_ij = 0.0_dp
351 DO ilist = 1, n_loc_size
352 kparticle = full_loc_list(2, ilist)
353 IF (kparticle == jparticle) cycle
354 cell_v_2 = loc_cell_v(:, ilist)
355 rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
356 drik = dot_product(rik, rik)
357 IF (drik > rab2_max) cycle
358 drik = sqrt(drik)
359 costheta = dot_product(rij, rik)/(drij*drik)
360 IF (costheta < -1.0_dp) costheta = -1.0_dp
361 IF (costheta > +1.0_dp) costheta = +1.0_dp
362 f_c = ter_f_c(tersoff, drik)
363 gterm = ter_g(tersoff, costheta)
364 expterm = exp((lambda3*(drij - drik))**3)
365 ter_zeta_ij = ter_zeta_ij + f_c*gterm*expterm
366 END DO
367
368 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 SUBROUTINE ter_zeta_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, prefactor, &
390 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 n = tersoff%n
412 lambda3 = tersoff%lambda3
413 rab2_max = rcutsq
414
415 rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
416 drij = sqrt(dot_product(rij, rij))
417 rij_hat(:) = rij(:)/drij
418
419 nparticle = SIZE(r_last_update_pbc)
420 DO ilist = 1, n_loc_size
421 kparticle = full_loc_list(2, ilist)
422 IF (kparticle == jparticle) cycle
423 cell_v_2 = loc_cell_v(:, ilist)
424 rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
425 drik = dot_product(rik, rik)
426
427 IF (drik > rab2_max) cycle
428 drik = sqrt(drik)
429 rik_hat(:) = rik(:)/drik
430 costheta = dot_product(rij, rik)/(drij*drik)
431 IF (costheta < -1.0_dp) costheta = -1.0_dp
432 IF (costheta > +1.0_dp) costheta = +1.0_dp
433
434 dcosdrj(:) = (1.0_dp/(drij))*(rik_hat(:) - costheta*rij_hat(:))
435 dcosdrk(:) = (1.0_dp/(drik))*(rij_hat(:) - costheta*rik_hat(:))
436 dcosdri(:) = -(dcosdrj(:) + dcosdrk(:))
437
438 f_c = ter_f_c(tersoff, drik)
439 f_c_d = ter_f_c_d(tersoff, drik)
440 gterm = ter_g(tersoff, costheta)
441 gterm_d = ter_g_d(tersoff, costheta) !still need d(costheta)/dR term
442 expterm = exp((lambda3*(drij - drik))**3)
443 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 + 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 + f_c*gterm*expterm_d*(rij_hat)
452
453 drk = f_c_d*gterm*expterm*(-rik) &
454 + f_c*gterm_d*expterm*(dcosdrk) &
455 + f_c*gterm*expterm_d*(-rik_hat)
456
457 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
458 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
459 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
460
461 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
462 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
463 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
464
465 f_nonbond(1, kparticle) = f_nonbond(1, kparticle) + prefactor*drk(1)
466 f_nonbond(2, kparticle) = f_nonbond(2, kparticle) + prefactor*drk(2)
467 f_nonbond(3, kparticle) = f_nonbond(3, kparticle) + prefactor*drk(3)
468
469 IF (use_virial) THEN
470 pv_nonbond(1, 1) = pv_nonbond(1, 1) + prefactor*(rij(1)*drj(1) + rik(1)*drk(1))
471 pv_nonbond(1, 2) = pv_nonbond(1, 2) + prefactor*(rij(1)*drj(2) + rik(1)*drk(2))
472 pv_nonbond(1, 3) = pv_nonbond(1, 3) + prefactor*(rij(1)*drj(3) + rik(1)*drk(3))
473
474 pv_nonbond(2, 1) = pv_nonbond(2, 1) + prefactor*(rij(2)*drj(1) + rik(2)*drk(1))
475 pv_nonbond(2, 2) = pv_nonbond(2, 2) + prefactor*(rij(2)*drj(2) + rik(2)*drk(2))
476 pv_nonbond(2, 3) = pv_nonbond(2, 3) + prefactor*(rij(2)*drj(3) + rik(2)*drk(3))
477
478 pv_nonbond(3, 1) = pv_nonbond(3, 1) + prefactor*(rij(3)*drj(1) + rik(3)*drk(1))
479 pv_nonbond(3, 2) = pv_nonbond(3, 2) + prefactor*(rij(3)*drj(2) + rik(3)*drk(2))
480 pv_nonbond(3, 3) = pv_nonbond(3, 3) + prefactor*(rij(3)*drj(3) + rik(3)*drk(3))
481 END IF
482 END DO
483 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 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 c = tersoff%c
500 d = tersoff%d
501 h = tersoff%h
502 c2 = c*c
503 d2 = d*d
504 ter_g = 0.0_dp
505 ter_g = 1.0_dp + (c2/d2) - (c2)/(d2 + (h - costheta)**2)
506
507 END FUNCTION ter_g
508
509! **************************************************************************************************
510!> \brief ...
511!> \param tersoff ...
512!> \param costheta ...
513!> \return ...
514!> \author I-Feng W. Kuo
515! **************************************************************************************************
516 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, sintheta
522
523 c = tersoff%c
524 d = tersoff%d
525 h = tersoff%h
526 c2 = c*c
527 d2 = d*d
528 hc = h - costheta
529
530 sintheta = sqrt(1.0 - costheta**2)
531
532 ter_g_d = 0.0_dp
533 ! Still need d(costheta)/dR
534 ter_g_d = (-2.0_dp*c2*hc)/(d2 + hc**2)**2
535 END FUNCTION ter_g_d
536
537! **************************************************************************************************
538!> \brief ...
539!> \param tersoff ...
540!> \param r_last_update_pbc ...
541!> \param cell_v ...
542!> \param n_loc_size ...
543!> \param full_loc_list ...
544!> \param loc_cell_v ...
545!> \param iparticle ...
546!> \param jparticle ...
547!> \param f_nonbond ...
548!> \param pv_nonbond ...
549!> \param use_virial ...
550!> \param rcutsq ...
551!> \par History
552!> Using a local list of neighbors - [tlaino] 2007
553!> \author I-Feng W. Kuo, Teodoro Laino
554! **************************************************************************************************
555 SUBROUTINE tersoff_forces(tersoff, r_last_update_pbc, cell_v, n_loc_size, &
556 full_loc_list, loc_cell_v, iparticle, jparticle, f_nonbond, pv_nonbond, &
557 use_virial, rcutsq)
558 TYPE(tersoff_pot_type), POINTER :: tersoff
559 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
560 REAL(kind=dp), DIMENSION(3) :: cell_v
561 INTEGER, INTENT(IN) :: n_loc_size
562 INTEGER, DIMENSION(2, 1:n_loc_size) :: full_loc_list
563 REAL(kind=dp), DIMENSION(3, 1:n_loc_size) :: loc_cell_v
564 INTEGER, INTENT(IN) :: iparticle, jparticle
565 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
566 LOGICAL, INTENT(IN) :: use_virial
567 REAL(kind=dp), INTENT(IN) :: rcutsq
568
569 CHARACTER(LEN=*), PARAMETER :: routinen = 'tersoff_forces'
570
571 INTEGER :: handle
572 REAL(kind=dp) :: b_ij, b_ij_d, drij, f_a, f_a1, f_a2, &
573 f_a_d, f_c, f_c_d, f_r, f_r1, f_r2, &
574 f_r_d, fac, prefactor, rij(3), &
575 rij_hat(3)
576
577 CALL timeset(routinen, handle)
578 rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
579 drij = sqrt(dot_product(rij, rij))
580 rij_hat(:) = rij(:)/drij
581
582 fac = -0.5_dp
583 b_ij = ter_b_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq)
584 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)
585 f_a = ter_f_a(tersoff, drij)
586 f_a_d = ter_f_a_d(tersoff, drij)
587 f_c = ter_f_c(tersoff, drij)
588 f_c_d = ter_f_c_d(tersoff, drij)
589 f_r = ter_f_r(tersoff, drij)
590 f_r_d = ter_f_r_d(tersoff, drij)
591
592 ! Lets do the easy one first, the repulsive term
593 ! Note a_ij = 1.0_dp so just going to ignore it...
594 f_r1 = f_c_d*f_r*fac
595 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_r1*rij(1)
596 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_r1*rij(2)
597 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_r1*rij(3)
598 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_r1*rij(1)
599 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_r1*rij(2)
600 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_r1*rij(3)
601
602 IF (use_virial) THEN
603 pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_r1*rij(1)*rij(1)
604 pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_r1*rij(1)*rij(2)
605 pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_r1*rij(1)*rij(3)
606 pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_r1*rij(2)*rij(1)
607 pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_r1*rij(2)*rij(2)
608 pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_r1*rij(2)*rij(3)
609 pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_r1*rij(3)*rij(1)
610 pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_r1*rij(3)*rij(2)
611 pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_r1*rij(3)*rij(3)
612 END IF
613
614 f_r2 = f_c*f_r_d*fac
615 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_r2*rij(1)
616 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_r2*rij(2)
617 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_r2*rij(3)
618 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_r2*rij(1)
619 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_r2*rij(2)
620 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_r2*rij(3)
621
622 IF (use_virial) THEN
623 pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_r2*rij(1)*rij(1)
624 pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_r2*rij(1)*rij(2)
625 pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_r2*rij(1)*rij(3)
626 pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_r2*rij(2)*rij(1)
627 pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_r2*rij(2)*rij(2)
628 pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_r2*rij(2)*rij(3)
629 pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_r2*rij(3)*rij(1)
630 pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_r2*rij(3)*rij(2)
631 pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_r2*rij(3)*rij(3)
632 END IF
633
634 ! Lets do the f_A1 piece derivative of F_C
635 f_a1 = f_c_d*b_ij*f_a*fac
636 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_a1*rij(1)
637 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_a1*rij(2)
638 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_a1*rij(3)
639 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_a1*rij(1)
640 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_a1*rij(2)
641 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_a1*rij(3)
642
643 IF (use_virial) THEN
644 pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_a1*rij(1)*rij(1)
645 pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_a1*rij(1)*rij(2)
646 pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_a1*rij(1)*rij(3)
647 pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_a1*rij(2)*rij(1)
648 pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_a1*rij(2)*rij(2)
649 pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_a1*rij(2)*rij(3)
650 pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_a1*rij(3)*rij(1)
651 pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_a1*rij(3)*rij(2)
652 pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_a1*rij(3)*rij(3)
653 END IF
654
655 ! Lets do the f_A2 piece derivative of F_A
656 f_a2 = f_c*b_ij*f_a_d*fac
657 f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_a2*rij(1)
658 f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_a2*rij(2)
659 f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_a2*rij(3)
660 f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_a2*rij(1)
661 f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_a2*rij(2)
662 f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_a2*rij(3)
663
664 IF (use_virial) THEN
665 pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_a2*rij(1)*rij(1)
666 pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_a2*rij(1)*rij(2)
667 pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_a2*rij(1)*rij(3)
668 pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_a2*rij(2)*rij(1)
669 pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_a2*rij(2)*rij(2)
670 pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_a2*rij(2)*rij(3)
671 pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_a2*rij(3)*rij(1)
672 pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_a2*rij(3)*rij(2)
673 pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_a2*rij(3)*rij(3)
674 END IF
675
676 ! Lets do the f_A3 piece derivative of b_ij
677 prefactor = f_c*b_ij_d*f_a*fac ! Note need to do d(Zeta_ij)/dR
678 CALL ter_zeta_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, prefactor, &
679 n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq, use_virial)
680 CALL timestop(handle)
681 END SUBROUTINE tersoff_forces
682
683! **************************************************************************************************
684!> \brief ...
685!> \param nonbonded ...
686!> \param potparm ...
687!> \param glob_loc_list ...
688!> \param glob_cell_v ...
689!> \param glob_loc_list_a ...
690!> \param cell ...
691!> \par History
692!> Fast implementation of the tersoff potential - [tlaino] 2007
693!> \author Teodoro Laino - University of Zurich
694! **************************************************************************************************
695 SUBROUTINE setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
696 TYPE(fist_neighbor_type), POINTER :: nonbonded
697 TYPE(pair_potential_pp_type), POINTER :: potparm
698 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
699 REAL(kind=dp), DIMENSION(:, :), POINTER :: glob_cell_v
700 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
701 TYPE(cell_type), POINTER :: cell
702
703 CHARACTER(LEN=*), PARAMETER :: routinen = 'setup_tersoff_arrays'
704
705 INTEGER :: handle, i, iend, igrp, ikind, ilist, &
706 ipair, istart, jkind, nkinds, npairs, &
707 npairs_tot
708 INTEGER, DIMENSION(:), POINTER :: work_list, work_list2
709 INTEGER, DIMENSION(:, :), POINTER :: list
710 REAL(kind=dp), DIMENSION(3) :: cell_v, cvi
711 REAL(kind=dp), DIMENSION(:, :), POINTER :: rwork_list
712 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
713 TYPE(pair_potential_single_type), POINTER :: pot
714
715 cpassert(.NOT. ASSOCIATED(glob_loc_list))
716 cpassert(.NOT. ASSOCIATED(glob_loc_list_a))
717 cpassert(.NOT. ASSOCIATED(glob_cell_v))
718 CALL timeset(routinen, handle)
719 npairs_tot = 0
720 nkinds = SIZE(potparm%pot, 1)
721 DO ilist = 1, nonbonded%nlists
722 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
723 npairs = neighbor_kind_pair%npairs
724 IF (npairs == 0) cycle
725 kind_group_loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
726 istart = neighbor_kind_pair%grp_kind_start(igrp)
727 iend = neighbor_kind_pair%grp_kind_end(igrp)
728 ikind = neighbor_kind_pair%ij_kind(1, igrp)
729 jkind = neighbor_kind_pair%ij_kind(2, igrp)
730 pot => potparm%pot(ikind, jkind)%pot
731 npairs = iend - istart + 1
732 IF (pot%no_mb) cycle
733 DO i = 1, SIZE(pot%type)
734 IF (pot%type(i) == tersoff_type) npairs_tot = npairs_tot + npairs
735 END DO
736 END DO kind_group_loop1
737 END DO
738 ALLOCATE (work_list(npairs_tot))
739 ALLOCATE (work_list2(npairs_tot))
740 ALLOCATE (glob_loc_list(2, npairs_tot))
741 ALLOCATE (glob_cell_v(3, npairs_tot))
742 ! Fill arrays with data
743 npairs_tot = 0
744 DO ilist = 1, nonbonded%nlists
745 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
746 npairs = neighbor_kind_pair%npairs
747 IF (npairs == 0) cycle
748 kind_group_loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
749 istart = neighbor_kind_pair%grp_kind_start(igrp)
750 iend = neighbor_kind_pair%grp_kind_end(igrp)
751 ikind = neighbor_kind_pair%ij_kind(1, igrp)
752 jkind = neighbor_kind_pair%ij_kind(2, igrp)
753 list => neighbor_kind_pair%list
754 cvi = neighbor_kind_pair%cell_vector
755 pot => potparm%pot(ikind, jkind)%pot
756 npairs = iend - istart + 1
757 IF (pot%no_mb) cycle
758 cell_v = matmul(cell%hmat, cvi)
759 DO i = 1, SIZE(pot%type)
760 ! TERSOFF
761 IF (pot%type(i) == tersoff_type) THEN
762 DO ipair = 1, npairs
763 glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
764 glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
765 END DO
766 npairs_tot = npairs_tot + npairs
767 END IF
768 END DO
769 END DO kind_group_loop2
770 END DO
771 ! Order the arrays w.r.t. the first index of glob_loc_list
772 CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
773 DO ipair = 1, npairs_tot
774 work_list2(ipair) = glob_loc_list(2, work_list(ipair))
775 END DO
776 glob_loc_list(2, :) = work_list2
777 DEALLOCATE (work_list2)
778 ALLOCATE (rwork_list(3, npairs_tot))
779 DO ipair = 1, npairs_tot
780 rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
781 END DO
782 glob_cell_v = rwork_list
783 DEALLOCATE (rwork_list)
784 DEALLOCATE (work_list)
785 ALLOCATE (glob_loc_list_a(npairs_tot))
786 glob_loc_list_a = glob_loc_list(1, :)
787 CALL timestop(handle)
788 END SUBROUTINE setup_tersoff_arrays
789
790! **************************************************************************************************
791!> \brief ...
792!> \param glob_loc_list ...
793!> \param glob_cell_v ...
794!> \param glob_loc_list_a ...
795!> \par History
796!> Fast implementation of the tersoff potential - [tlaino] 2007
797!> \author Teodoro Laino - University of Zurich
798! **************************************************************************************************
799 SUBROUTINE destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
800 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list
801 REAL(kind=dp), DIMENSION(:, :), POINTER :: glob_cell_v
802 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a
803
804 IF (ASSOCIATED(glob_loc_list)) THEN
805 DEALLOCATE (glob_loc_list)
806 END IF
807 IF (ASSOCIATED(glob_loc_list_a)) THEN
808 DEALLOCATE (glob_loc_list_a)
809 END IF
810 IF (ASSOCIATED(glob_cell_v)) THEN
811 DEALLOCATE (glob_cell_v)
812 END IF
813
814 END SUBROUTINE destroy_tersoff_arrays
815
816END MODULE manybody_tersoff
817
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
fist_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition fist_neighbor_list_types.F:11
fist_nonbond_env_types
Definition fist_nonbond_env_types.F:13
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
manybody_tersoff
Definition manybody_tersoff.F:13
manybody_tersoff::setup_tersoff_arrays
subroutine, public setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition manybody_tersoff.F:696
manybody_tersoff::tersoff_forces
subroutine, public tersoff_forces(tersoff, r_last_update_pbc, cell_v, n_loc_size, full_loc_list, loc_cell_v, iparticle, jparticle, f_nonbond, pv_nonbond, use_virial, rcutsq)
...
Definition manybody_tersoff.F:558
manybody_tersoff::destroy_tersoff_arrays
subroutine, public destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition manybody_tersoff.F:800
manybody_tersoff::tersoff_energy
subroutine, public tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, full_loc_list, loc_cell_v, cell_v, drij)
...
Definition manybody_tersoff.F:53
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
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::tersoff_type
integer, parameter, public tersoff_type
Definition pair_potential_types.F:35
util
All kind of helpful little routines.
Definition util.F:14
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
fist_neighbor_list_types::fist_neighbor_type
Definition fist_neighbor_list_types.F:48
fist_neighbor_list_types::neighbor_kind_pairs_type
Definition fist_neighbor_list_types.F:26
fist_nonbond_env_types::pos_type
Definition fist_nonbond_env_types.F:44
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:412
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:379
pair_potential_types::tersoff_pot_type
Definition pair_potential_types.F:250