Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 330372 : SUBROUTINE tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
52 330372 : 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 330372 : full_loc_list, loc_cell_v, cell_v, tersoff%rcutsq)
67 330372 : f_C = ter_f_C(tersoff, drij)
68 330372 : f_A = ter_f_A(tersoff, drij)
69 330372 : f_R = ter_f_R(tersoff, drij)
70 330372 : pot_loc = f_C*(f_R + b_ij*f_A)
71 :
72 330372 : END SUBROUTINE tersoff_energy
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param tersoff ...
77 : !> \param r ...
78 : !> \return ...
79 : !> \author I-Feng W. Kuo
80 : ! **************************************************************************************************
81 6383136 : 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 6383136 : bigR = tersoff%bigR
89 6383136 : bigD = tersoff%bigD
90 6383136 : RmD = tersoff%bigR - tersoff%bigD
91 6383136 : RpD = tersoff%bigR + tersoff%bigD
92 6383136 : ter_f_C = 0.0_dp
93 6383136 : IF (r < RmD) ter_f_C = 1.0_dp
94 6383136 : IF (r > RpD) ter_f_C = 0.0_dp
95 6383136 : IF ((r < RpD) .AND. (r > RmD)) THEN
96 1876824 : ter_f_C = 0.5_dp*(1.0_dp - SIN(0.5_dp*PI*(r - bigR)/(bigD)))
97 : END IF
98 6383136 : 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 1760970 : 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 1760970 : bigR = tersoff%bigR
115 1760970 : bigD = tersoff%bigD
116 1760970 : RmD = tersoff%bigR - tersoff%bigD
117 1760970 : 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 1760970 : IF ((r < RpD) .AND. (r > RmD)) THEN
122 496966 : ter_f_C_d = (0.25_dp*PI/bigD)*COS(0.5_dp*PI*(r - bigR)/(bigD))/r
123 : END IF
124 :
125 1760970 : 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 660744 : 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 660744 : A = tersoff%A
142 660744 : lambda1 = tersoff%lambda1
143 660744 : ter_f_R = 0.0_dp
144 660744 : ter_f_R = A*EXP(-lambda1*r)
145 :
146 660744 : 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 330372 : 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 330372 : A = tersoff%A
163 330372 : lambda1 = tersoff%lambda1
164 330372 : f_R = A*EXP(-lambda1*r)
165 330372 : ter_f_R_d = 0.0_dp
166 330372 : ter_f_R_d = lambda1*f_R/r
167 :
168 330372 : 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 660744 : 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 660744 : B = tersoff%B
185 660744 : lambda2 = tersoff%lambda2
186 660744 : ter_f_A = 0.0_dp
187 660744 : ter_f_A = -B*EXP(-lambda2*r)
188 :
189 660744 : 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 330372 : 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 330372 : B = tersoff%B
206 330372 : lambda2 = tersoff%lambda2
207 330372 : ter_f_A_d = 0.0_dp
208 330372 : ter_f_A_d = -B*lambda2*EXP(-lambda2*r)/r
209 :
210 330372 : 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 660744 : FUNCTION ter_b_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
248 660744 : 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 660744 : n = tersoff%n
261 660744 : beta = tersoff%beta
262 660744 : ter_b_ij = 0.0_dp
263 : zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, &
264 660744 : n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq)
265 660744 : ter_b_ij = (1.0_dp + (beta*zeta_ij)**n)**(-0.5_dp/n)
266 :
267 660744 : 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 330372 : FUNCTION ter_b_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
284 330372 : 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 330372 : n = tersoff%n
298 330372 : beta = tersoff%beta
299 330372 : beta_n = beta**n
300 : zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
301 330372 : full_loc_list, loc_cell_v, cell_v, rcutsq)
302 330372 : zeta_ij_nm1 = 0.0_dp
303 330372 : IF (zeta_ij > 0.0_dp) zeta_ij_nm1 = zeta_ij**(n - 1.0_dp)
304 330372 : zeta_ij_n = zeta_ij**(n)
305 :
306 330372 : ter_b_ij_d = 0.0_dp
307 : ter_b_ij_d = -0.5_dp*beta_n*zeta_ij_nm1* &
308 330372 : ((1.0_dp + beta_n*zeta_ij_n)**((-0.5_dp/n) - 1.0_dp))
309 :
310 330372 : 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 991116 : FUNCTION ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
329 991116 : 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 991116 : ter_zeta_ij = 0.0_dp
345 991116 : n = tersoff%n
346 991116 : lambda3 = tersoff%lambda3
347 991116 : rab2_max = rcutsq
348 3964464 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
349 3964464 : drij = SQRT(DOT_PRODUCT(rij, rij))
350 991116 : ter_zeta_ij = 0.0_dp
351 132213636 : DO ilist = 1, n_loc_size
352 131222520 : kparticle = full_loc_list(2, ilist)
353 131222520 : IF (kparticle == jparticle) CYCLE
354 517521288 : cell_v_2 = loc_cell_v(:, ilist)
355 517521288 : rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
356 517521288 : drik = DOT_PRODUCT(rik, rik)
357 129380322 : IF (drik > rab2_max) CYCLE
358 4291794 : drik = SQRT(drik)
359 17167176 : costheta = DOT_PRODUCT(rij, rik)/(drij*drik)
360 4291794 : IF (costheta < -1.0_dp) costheta = -1.0_dp
361 4291794 : IF (costheta > +1.0_dp) costheta = +1.0_dp
362 4291794 : f_C = ter_f_C(tersoff, drik)
363 4291794 : gterm = ter_g(tersoff, costheta)
364 4291794 : expterm = EXP((lambda3*(drij - drik))**3)
365 132213636 : ter_zeta_ij = ter_zeta_ij + f_C*gterm*expterm
366 : END DO
367 :
368 991116 : 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 330372 : SUBROUTINE ter_zeta_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, prefactor, &
390 330372 : 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 330372 : n = tersoff%n
412 330372 : lambda3 = tersoff%lambda3
413 330372 : rab2_max = rcutsq
414 :
415 1321488 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
416 1321488 : drij = SQRT(DOT_PRODUCT(rij, rij))
417 1321488 : rij_hat(:) = rij(:)/drij
418 :
419 44071212 : nparticle = SIZE(r_last_update_pbc)
420 44071212 : DO ilist = 1, n_loc_size
421 43740840 : kparticle = full_loc_list(2, ilist)
422 43740840 : IF (kparticle == jparticle) CYCLE
423 172507096 : cell_v_2 = loc_cell_v(:, ilist)
424 172507096 : rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
425 172507096 : drik = DOT_PRODUCT(rik, rik)
426 :
427 43126774 : IF (drik > rab2_max) CYCLE
428 1430598 : drik = SQRT(drik)
429 5722392 : rik_hat(:) = rik(:)/drik
430 5722392 : costheta = DOT_PRODUCT(rij, rik)/(drij*drik)
431 1430598 : IF (costheta < -1.0_dp) costheta = -1.0_dp
432 1430598 : IF (costheta > +1.0_dp) costheta = +1.0_dp
433 :
434 5722392 : dcosdrj(:) = (1.0_dp/(drij))*(rik_hat(:) - costheta*rij_hat(:))
435 5722392 : dcosdrk(:) = (1.0_dp/(drik))*(rij_hat(:) - costheta*rik_hat(:))
436 5722392 : dcosdri(:) = -(dcosdrj(:) + dcosdrk(:))
437 :
438 1430598 : f_C = ter_f_C(tersoff, drik)
439 1430598 : f_C_d = ter_f_C_d(tersoff, drik)
440 1430598 : gterm = ter_g(tersoff, costheta)
441 1430598 : gterm_d = ter_g_d(tersoff, costheta) !still need d(costheta)/dR term
442 1430598 : expterm = EXP((lambda3*(drij - drik))**3)
443 1430598 : 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 5722392 : + 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 5722392 : + f_C*gterm*expterm_d*(rij_hat)
452 :
453 : drk = f_C_d*gterm*expterm*(-rik) &
454 : + f_C*gterm_d*expterm*(dcosdrk) &
455 5722392 : + f_C*gterm*expterm_d*(-rik_hat)
456 :
457 1430598 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
458 1430598 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
459 1430598 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
460 :
461 1430598 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
462 1430598 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
463 1430598 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
464 :
465 1430598 : f_nonbond(1, kparticle) = f_nonbond(1, kparticle) + prefactor*drk(1)
466 1430598 : f_nonbond(2, kparticle) = f_nonbond(2, kparticle) + prefactor*drk(2)
467 1430598 : f_nonbond(3, kparticle) = f_nonbond(3, kparticle) + prefactor*drk(3)
468 :
469 1760970 : IF (use_virial) THEN
470 690052 : pv_nonbond(1, 1) = pv_nonbond(1, 1) + prefactor*(rij(1)*drj(1) + rik(1)*drk(1))
471 690052 : pv_nonbond(1, 2) = pv_nonbond(1, 2) + prefactor*(rij(1)*drj(2) + rik(1)*drk(2))
472 690052 : pv_nonbond(1, 3) = pv_nonbond(1, 3) + prefactor*(rij(1)*drj(3) + rik(1)*drk(3))
473 :
474 690052 : pv_nonbond(2, 1) = pv_nonbond(2, 1) + prefactor*(rij(2)*drj(1) + rik(2)*drk(1))
475 690052 : pv_nonbond(2, 2) = pv_nonbond(2, 2) + prefactor*(rij(2)*drj(2) + rik(2)*drk(2))
476 690052 : pv_nonbond(2, 3) = pv_nonbond(2, 3) + prefactor*(rij(2)*drj(3) + rik(2)*drk(3))
477 :
478 690052 : pv_nonbond(3, 1) = pv_nonbond(3, 1) + prefactor*(rij(3)*drj(1) + rik(3)*drk(1))
479 690052 : pv_nonbond(3, 2) = pv_nonbond(3, 2) + prefactor*(rij(3)*drj(2) + rik(3)*drk(2))
480 690052 : pv_nonbond(3, 3) = pv_nonbond(3, 3) + prefactor*(rij(3)*drj(3) + rik(3)*drk(3))
481 : END IF
482 : END DO
483 330372 : 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 5722392 : 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 5722392 : c = tersoff%c
500 5722392 : d = tersoff%d
501 5722392 : h = tersoff%h
502 5722392 : c2 = c*c
503 5722392 : d2 = d*d
504 5722392 : ter_g = 0.0_dp
505 5722392 : ter_g = 1.0_dp + (c2/d2) - (c2)/(d2 + (h - costheta)**2)
506 :
507 5722392 : END FUNCTION ter_g
508 :
509 : ! **************************************************************************************************
510 : !> \brief ...
511 : !> \param tersoff ...
512 : !> \param costheta ...
513 : !> \return ...
514 : !> \author I-Feng W. Kuo
515 : ! **************************************************************************************************
516 1430598 : 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 1430598 : c = tersoff%c
524 1430598 : d = tersoff%d
525 1430598 : h = tersoff%h
526 1430598 : c2 = c*c
527 1430598 : d2 = d*d
528 1430598 : hc = h - costheta
529 :
530 : sintheta = SQRT(1.0 - costheta**2)
531 :
532 1430598 : ter_g_d = 0.0_dp
533 : ! Still need d(costheta)/dR
534 1430598 : ter_g_d = (-2.0_dp*c2*hc)/(d2 + hc**2)**2
535 1430598 : 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 660744 : SUBROUTINE tersoff_forces(tersoff, r_last_update_pbc, cell_v, n_loc_size, &
556 330372 : 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 330372 : CALL timeset(routineN, handle)
578 1321488 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
579 1321488 : drij = SQRT(DOT_PRODUCT(rij, rij))
580 1321488 : rij_hat(:) = rij(:)/drij
581 :
582 330372 : fac = -0.5_dp
583 330372 : 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 330372 : 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 330372 : f_A = ter_f_A(tersoff, drij)
586 330372 : f_A_d = ter_f_A_d(tersoff, drij)
587 330372 : f_C = ter_f_C(tersoff, drij)
588 330372 : f_C_d = ter_f_C_d(tersoff, drij)
589 330372 : f_R = ter_f_R(tersoff, drij)
590 330372 : 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 330372 : f_R1 = f_C_d*f_R*fac
595 330372 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_R1*rij(1)
596 330372 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_R1*rij(2)
597 330372 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_R1*rij(3)
598 330372 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_R1*rij(1)
599 330372 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_R1*rij(2)
600 330372 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_R1*rij(3)
601 :
602 330372 : IF (use_virial) THEN
603 154574 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_R1*rij(1)*rij(1)
604 154574 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_R1*rij(1)*rij(2)
605 154574 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_R1*rij(1)*rij(3)
606 154574 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_R1*rij(2)*rij(1)
607 154574 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_R1*rij(2)*rij(2)
608 154574 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_R1*rij(2)*rij(3)
609 154574 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_R1*rij(3)*rij(1)
610 154574 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_R1*rij(3)*rij(2)
611 154574 : pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_R1*rij(3)*rij(3)
612 : END IF
613 :
614 330372 : f_R2 = f_C*f_R_d*fac
615 330372 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_R2*rij(1)
616 330372 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_R2*rij(2)
617 330372 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_R2*rij(3)
618 330372 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_R2*rij(1)
619 330372 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_R2*rij(2)
620 330372 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_R2*rij(3)
621 :
622 330372 : IF (use_virial) THEN
623 154574 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_R2*rij(1)*rij(1)
624 154574 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_R2*rij(1)*rij(2)
625 154574 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_R2*rij(1)*rij(3)
626 154574 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_R2*rij(2)*rij(1)
627 154574 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_R2*rij(2)*rij(2)
628 154574 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_R2*rij(2)*rij(3)
629 154574 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_R2*rij(3)*rij(1)
630 154574 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_R2*rij(3)*rij(2)
631 154574 : 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 330372 : f_A1 = f_C_d*b_ij*f_A*fac
636 330372 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A1*rij(1)
637 330372 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A1*rij(2)
638 330372 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A1*rij(3)
639 330372 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A1*rij(1)
640 330372 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A1*rij(2)
641 330372 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A1*rij(3)
642 :
643 330372 : IF (use_virial) THEN
644 154574 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_A1*rij(1)*rij(1)
645 154574 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_A1*rij(1)*rij(2)
646 154574 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_A1*rij(1)*rij(3)
647 154574 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_A1*rij(2)*rij(1)
648 154574 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_A1*rij(2)*rij(2)
649 154574 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_A1*rij(2)*rij(3)
650 154574 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_A1*rij(3)*rij(1)
651 154574 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_A1*rij(3)*rij(2)
652 154574 : 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 330372 : f_A2 = f_C*b_ij*f_A_d*fac
657 330372 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A2*rij(1)
658 330372 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A2*rij(2)
659 330372 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A2*rij(3)
660 330372 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A2*rij(1)
661 330372 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A2*rij(2)
662 330372 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A2*rij(3)
663 :
664 330372 : IF (use_virial) THEN
665 154574 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_A2*rij(1)*rij(1)
666 154574 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_A2*rij(1)*rij(2)
667 154574 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_A2*rij(1)*rij(3)
668 154574 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_A2*rij(2)*rij(1)
669 154574 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_A2*rij(2)*rij(2)
670 154574 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_A2*rij(2)*rij(3)
671 154574 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_A2*rij(3)*rij(1)
672 154574 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_A2*rij(3)*rij(2)
673 154574 : 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 330372 : 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 330372 : n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq, use_virial)
680 330372 : CALL timestop(handle)
681 330372 : 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 6248 : 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 6248 : INTEGER, DIMENSION(:), POINTER :: work_list, work_list2
709 6248 : INTEGER, DIMENSION(:, :), POINTER :: list
710 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi
711 6248 : 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 0 : CPASSERT(.NOT. ASSOCIATED(glob_loc_list))
716 6248 : CPASSERT(.NOT. ASSOCIATED(glob_loc_list_a))
717 6248 : CPASSERT(.NOT. ASSOCIATED(glob_cell_v))
718 6248 : CALL timeset(routineN, handle)
719 6248 : npairs_tot = 0
720 6248 : nkinds = SIZE(potparm%pot, 1)
721 231000 : DO ilist = 1, nonbonded%nlists
722 224752 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
723 224752 : npairs = neighbor_kind_pair%npairs
724 224752 : IF (npairs == 0) CYCLE
725 175462 : Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
726 85728 : istart = neighbor_kind_pair%grp_kind_start(igrp)
727 85728 : iend = neighbor_kind_pair%grp_kind_end(igrp)
728 85728 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
729 85728 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
730 85728 : pot => potparm%pot(ikind, jkind)%pot
731 85728 : npairs = iend - istart + 1
732 85728 : IF (pot%no_mb) CYCLE
733 396068 : DO i = 1, SIZE(pot%type)
734 171388 : IF (pot%type(i) == tersoff_type) npairs_tot = npairs_tot + npairs
735 : END DO
736 : END DO Kind_Group_Loop1
737 : END DO
738 18744 : ALLOCATE (work_list(npairs_tot))
739 12496 : ALLOCATE (work_list2(npairs_tot))
740 18744 : ALLOCATE (glob_loc_list(2, npairs_tot))
741 18744 : ALLOCATE (glob_cell_v(3, npairs_tot))
742 : ! Fill arrays with data
743 6248 : npairs_tot = 0
744 231000 : DO ilist = 1, nonbonded%nlists
745 224752 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
746 224752 : npairs = neighbor_kind_pair%npairs
747 224752 : IF (npairs == 0) CYCLE
748 175462 : Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
749 85728 : istart = neighbor_kind_pair%grp_kind_start(igrp)
750 85728 : iend = neighbor_kind_pair%grp_kind_end(igrp)
751 85728 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
752 85728 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
753 85728 : list => neighbor_kind_pair%list
754 342912 : cvi = neighbor_kind_pair%cell_vector
755 85728 : pot => potparm%pot(ikind, jkind)%pot
756 85728 : npairs = iend - istart + 1
757 85728 : IF (pot%no_mb) CYCLE
758 1113528 : cell_v = MATMUL(cell%hmat, cvi)
759 396068 : DO i = 1, SIZE(pot%type)
760 : ! TERSOFF
761 171388 : IF (pot%type(i) == tersoff_type) THEN
762 15378806 : DO ipair = 1, npairs
763 91758960 : glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
764 61258286 : glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
765 : END DO
766 85646 : 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 6248 : CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
773 15299408 : DO ipair = 1, npairs_tot
774 15299408 : work_list2(ipair) = glob_loc_list(2, work_list(ipair))
775 : END DO
776 30598816 : glob_loc_list(2, :) = work_list2
777 6248 : DEALLOCATE (work_list2)
778 18744 : ALLOCATE (rwork_list(3, npairs_tot))
779 15299408 : DO ipair = 1, npairs_tot
780 122351528 : rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
781 : END DO
782 122357776 : glob_cell_v = rwork_list
783 6248 : DEALLOCATE (rwork_list)
784 6248 : DEALLOCATE (work_list)
785 18744 : ALLOCATE (glob_loc_list_a(npairs_tot))
786 30598816 : glob_loc_list_a = glob_loc_list(1, :)
787 6248 : CALL timestop(handle)
788 12496 : 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 6248 : 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 6248 : IF (ASSOCIATED(glob_loc_list)) THEN
805 6248 : DEALLOCATE (glob_loc_list)
806 : END IF
807 6248 : IF (ASSOCIATED(glob_loc_list_a)) THEN
808 6248 : DEALLOCATE (glob_loc_list_a)
809 : END IF
810 6248 : IF (ASSOCIATED(glob_cell_v)) THEN
811 6248 : DEALLOCATE (glob_cell_v)
812 : END IF
813 :
814 6248 : END SUBROUTINE destroy_tersoff_arrays
815 :
816 : END MODULE manybody_tersoff
817 :
|
