(git:0de0cc2)
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 ! **************************************************************************************************
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,&
25  USE util, ONLY: sort
26 #include "./base/base_uses.f90"
27 
28  IMPLICIT NONE
29 
30  PRIVATE
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  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 
816 END MODULE manybody_tersoff
817 
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
Handles all functions related to the CELL.
Definition: cell_types.F:15
Define the neighbor list data types and the corresponding functionality.
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition: list.F:24
subroutine, public setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
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)
...
subroutine, public destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
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 of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
integer, parameter, public tersoff_type
All kind of helpful little routines.
Definition: util.F:14