Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> Efficient tersoff implementation
11 : !> \author CJM, I-Feng W. Kuo, Teodoro Laino
12 : ! **************************************************************************************************
13 : MODULE manybody_tersoff
14 :
15 : USE cell_types, ONLY: cell_type
16 : USE fist_neighbor_list_types, ONLY: fist_neighbor_type,&
17 : neighbor_kind_pairs_type
18 : USE fist_nonbond_env_types, ONLY: pos_type
19 : USE kinds, ONLY: dp
20 : USE mathconstants, ONLY: pi
21 : USE pair_potential_types, ONLY: pair_potential_pp_type,&
22 : pair_potential_single_type,&
23 : tersoff_pot_type,&
24 : tersoff_type
25 : USE util, ONLY: sort
26 : #include "./base/base_uses.f90"
27 :
28 : IMPLICIT NONE
29 :
30 : PRIVATE
31 : PUBLIC :: setup_tersoff_arrays, destroy_tersoff_arrays, &
32 : tersoff_forces, tersoff_energy
33 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_tersoff'
34 :
35 : CONTAINS
36 :
37 : ! **************************************************************************************************
38 : !> \brief ...
39 : !> \param pot_loc ...
40 : !> \param tersoff ...
41 : !> \param r_last_update_pbc ...
42 : !> \param atom_a ...
43 : !> \param atom_b ...
44 : !> \param nloc_size ...
45 : !> \param full_loc_list ...
46 : !> \param loc_cell_v ...
47 : !> \param cell_v ...
48 : !> \param drij ...
49 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
50 : ! **************************************************************************************************
51 250380 : SUBROUTINE tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
52 250380 : 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 250380 : full_loc_list, loc_cell_v, cell_v, tersoff%rcutsq)
67 250380 : f_C = ter_f_C(tersoff, drij)
68 250380 : f_A = ter_f_A(tersoff, drij)
69 250380 : f_R = ter_f_R(tersoff, drij)
70 250380 : pot_loc = f_C*(f_R + b_ij*f_A)
71 :
72 250380 : END SUBROUTINE tersoff_energy
73 :
74 : ! **************************************************************************************************
75 : !> \brief ...
76 : !> \param tersoff ...
77 : !> \param r ...
78 : !> \return ...
79 : !> \author I-Feng W. Kuo
80 : ! **************************************************************************************************
81 4601928 : 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 4601928 : bigR = tersoff%bigR
89 4601928 : bigD = tersoff%bigD
90 4601928 : RmD = tersoff%bigR - tersoff%bigD
91 4601928 : RpD = tersoff%bigR + tersoff%bigD
92 4601928 : ter_f_C = 0.0_dp
93 4601928 : IF (r < RmD) ter_f_C = 1.0_dp
94 4601928 : IF (r > RpD) ter_f_C = 0.0_dp
95 4601928 : IF ((r < RpD) .AND. (r > RmD)) THEN
96 1761012 : ter_f_C = 0.5_dp*(1.0_dp - SIN(0.5_dp*PI*(r - bigR)/(bigD)))
97 : END IF
98 4601928 : 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 1275672 : 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 1275672 : bigR = tersoff%bigR
115 1275672 : bigD = tersoff%bigD
116 1275672 : RmD = tersoff%bigR - tersoff%bigD
117 1275672 : 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 1275672 : IF ((r < RpD) .AND. (r > RmD)) THEN
122 465791 : ter_f_C_d = (0.25_dp*PI/bigD)*COS(0.5_dp*PI*(r - bigR)/(bigD))/r
123 : END IF
124 :
125 1275672 : 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 500760 : 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 500760 : A = tersoff%A
142 500760 : lambda1 = tersoff%lambda1
143 500760 : ter_f_R = 0.0_dp
144 500760 : ter_f_R = A*EXP(-lambda1*r)
145 :
146 500760 : 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 250380 : 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 250380 : A = tersoff%A
163 250380 : lambda1 = tersoff%lambda1
164 250380 : f_R = A*EXP(-lambda1*r)
165 250380 : ter_f_R_d = 0.0_dp
166 250380 : ter_f_R_d = lambda1*f_R/r
167 :
168 250380 : 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 500760 : 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 500760 : B = tersoff%B
185 500760 : lambda2 = tersoff%lambda2
186 500760 : ter_f_A = 0.0_dp
187 500760 : ter_f_A = -B*EXP(-lambda2*r)
188 :
189 500760 : 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 250380 : 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 250380 : B = tersoff%B
206 250380 : lambda2 = tersoff%lambda2
207 250380 : ter_f_A_d = 0.0_dp
208 250380 : ter_f_A_d = -B*lambda2*EXP(-lambda2*r)/r
209 :
210 250380 : 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 500760 : FUNCTION ter_b_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
248 500760 : 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 500760 : n = tersoff%n
261 500760 : beta = tersoff%beta
262 500760 : ter_b_ij = 0.0_dp
263 : zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, &
264 500760 : n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq)
265 500760 : ter_b_ij = (1.0_dp + (beta*zeta_ij)**n)**(-0.5_dp/n)
266 :
267 500760 : 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 250380 : FUNCTION ter_b_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
284 250380 : 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 250380 : n = tersoff%n
298 250380 : beta = tersoff%beta
299 250380 : beta_n = beta**n
300 : zeta_ij = ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
301 250380 : full_loc_list, loc_cell_v, cell_v, rcutsq)
302 250380 : zeta_ij_nm1 = 0.0_dp
303 250380 : IF (zeta_ij > 0.0_dp) zeta_ij_nm1 = zeta_ij**(n - 1.0_dp)
304 250380 : zeta_ij_n = zeta_ij**(n)
305 :
306 250380 : ter_b_ij_d = 0.0_dp
307 : ter_b_ij_d = -0.5_dp*beta_n*zeta_ij_nm1* &
308 250380 : ((1.0_dp + beta_n*zeta_ij_n)**((-0.5_dp/n) - 1.0_dp))
309 :
310 250380 : 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 751140 : FUNCTION ter_zeta_ij(tersoff, r_last_update_pbc, iparticle, jparticle, n_loc_size, &
329 751140 : 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 751140 : ter_zeta_ij = 0.0_dp
345 751140 : n = tersoff%n
346 751140 : lambda3 = tersoff%lambda3
347 751140 : rab2_max = rcutsq
348 3004560 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
349 3004560 : drij = SQRT(DOT_PRODUCT(rij, rij))
350 751140 : ter_zeta_ij = 0.0_dp
351 81945396 : DO ilist = 1, n_loc_size
352 81194256 : kparticle = full_loc_list(2, ilist)
353 81194256 : IF (kparticle == jparticle) CYCLE
354 318489000 : cell_v_2 = loc_cell_v(:, ilist)
355 318489000 : rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
356 318489000 : drik = DOT_PRODUCT(rik, rik)
357 79622250 : IF (drik > rab2_max) CYCLE
358 3075876 : drik = SQRT(drik)
359 12303504 : costheta = DOT_PRODUCT(rij, rik)/(drij*drik)
360 3075876 : IF (costheta < -1.0_dp) costheta = -1.0_dp
361 3075876 : IF (costheta > +1.0_dp) costheta = +1.0_dp
362 3075876 : f_C = ter_f_C(tersoff, drik)
363 3075876 : gterm = ter_g(tersoff, costheta)
364 3075876 : expterm = EXP((lambda3*(drij - drik))**3)
365 81945396 : ter_zeta_ij = ter_zeta_ij + f_C*gterm*expterm
366 : END DO
367 :
368 751140 : 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 250380 : SUBROUTINE ter_zeta_ij_d(tersoff, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, prefactor, &
390 250380 : 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 250380 : n = tersoff%n
412 250380 : lambda3 = tersoff%lambda3
413 250380 : rab2_max = rcutsq
414 :
415 1001520 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
416 1001520 : drij = SQRT(DOT_PRODUCT(rij, rij))
417 1001520 : rij_hat(:) = rij(:)/drij
418 :
419 27315132 : nparticle = SIZE(r_last_update_pbc)
420 27315132 : DO ilist = 1, n_loc_size
421 27064752 : kparticle = full_loc_list(2, ilist)
422 27064752 : IF (kparticle == jparticle) CYCLE
423 106163000 : cell_v_2 = loc_cell_v(:, ilist)
424 106163000 : rik(:) = r_last_update_pbc(kparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v_2
425 106163000 : drik = DOT_PRODUCT(rik, rik)
426 :
427 26540750 : IF (drik > rab2_max) CYCLE
428 1025292 : drik = SQRT(drik)
429 4101168 : rik_hat(:) = rik(:)/drik
430 4101168 : costheta = DOT_PRODUCT(rij, rik)/(drij*drik)
431 1025292 : IF (costheta < -1.0_dp) costheta = -1.0_dp
432 1025292 : IF (costheta > +1.0_dp) costheta = +1.0_dp
433 :
434 4101168 : dcosdrj(:) = (1.0_dp/(drij))*(rik_hat(:) - costheta*rij_hat(:))
435 4101168 : dcosdrk(:) = (1.0_dp/(drik))*(rij_hat(:) - costheta*rik_hat(:))
436 4101168 : dcosdri(:) = -(dcosdrj(:) + dcosdrk(:))
437 :
438 1025292 : f_C = ter_f_C(tersoff, drik)
439 1025292 : f_C_d = ter_f_C_d(tersoff, drik)
440 1025292 : gterm = ter_g(tersoff, costheta)
441 1025292 : gterm_d = ter_g_d(tersoff, costheta) !still need d(costheta)/dR term
442 1025292 : expterm = EXP((lambda3*(drij - drik))**3)
443 1025292 : 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 4101168 : + 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 4101168 : + f_C*gterm*expterm_d*(rij_hat)
452 :
453 : drk = f_C_d*gterm*expterm*(-rik) &
454 : + f_C*gterm_d*expterm*(dcosdrk) &
455 4101168 : + f_C*gterm*expterm_d*(-rik_hat)
456 :
457 1025292 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + prefactor*dri(1)
458 1025292 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + prefactor*dri(2)
459 1025292 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + prefactor*dri(3)
460 :
461 1025292 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) + prefactor*drj(1)
462 1025292 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) + prefactor*drj(2)
463 1025292 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) + prefactor*drj(3)
464 :
465 1025292 : f_nonbond(1, kparticle) = f_nonbond(1, kparticle) + prefactor*drk(1)
466 1025292 : f_nonbond(2, kparticle) = f_nonbond(2, kparticle) + prefactor*drk(2)
467 1025292 : f_nonbond(3, kparticle) = f_nonbond(3, kparticle) + prefactor*drk(3)
468 :
469 1275672 : IF (use_virial) THEN
470 284746 : pv_nonbond(1, 1) = pv_nonbond(1, 1) + prefactor*(rij(1)*drj(1) + rik(1)*drk(1))
471 284746 : pv_nonbond(1, 2) = pv_nonbond(1, 2) + prefactor*(rij(1)*drj(2) + rik(1)*drk(2))
472 284746 : pv_nonbond(1, 3) = pv_nonbond(1, 3) + prefactor*(rij(1)*drj(3) + rik(1)*drk(3))
473 :
474 284746 : pv_nonbond(2, 1) = pv_nonbond(2, 1) + prefactor*(rij(2)*drj(1) + rik(2)*drk(1))
475 284746 : pv_nonbond(2, 2) = pv_nonbond(2, 2) + prefactor*(rij(2)*drj(2) + rik(2)*drk(2))
476 284746 : pv_nonbond(2, 3) = pv_nonbond(2, 3) + prefactor*(rij(2)*drj(3) + rik(2)*drk(3))
477 :
478 284746 : pv_nonbond(3, 1) = pv_nonbond(3, 1) + prefactor*(rij(3)*drj(1) + rik(3)*drk(1))
479 284746 : pv_nonbond(3, 2) = pv_nonbond(3, 2) + prefactor*(rij(3)*drj(2) + rik(3)*drk(2))
480 284746 : pv_nonbond(3, 3) = pv_nonbond(3, 3) + prefactor*(rij(3)*drj(3) + rik(3)*drk(3))
481 : END IF
482 : END DO
483 250380 : 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 4101168 : 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 4101168 : c = tersoff%c
500 4101168 : d = tersoff%d
501 4101168 : h = tersoff%h
502 4101168 : c2 = c*c
503 4101168 : d2 = d*d
504 4101168 : ter_g = 0.0_dp
505 4101168 : ter_g = 1.0_dp + (c2/d2) - (c2)/(d2 + (h - costheta)**2)
506 :
507 4101168 : END FUNCTION ter_g
508 :
509 : ! **************************************************************************************************
510 : !> \brief ...
511 : !> \param tersoff ...
512 : !> \param costheta ...
513 : !> \return ...
514 : !> \author I-Feng W. Kuo
515 : ! **************************************************************************************************
516 1025292 : 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 1025292 : c = tersoff%c
524 1025292 : d = tersoff%d
525 1025292 : h = tersoff%h
526 1025292 : c2 = c*c
527 1025292 : d2 = d*d
528 1025292 : hc = h - costheta
529 :
530 : sintheta = SQRT(1.0 - costheta**2)
531 :
532 1025292 : ter_g_d = 0.0_dp
533 : ! Still need d(costheta)/dR
534 1025292 : ter_g_d = (-2.0_dp*c2*hc)/(d2 + hc**2)**2
535 1025292 : 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 500760 : SUBROUTINE tersoff_forces(tersoff, r_last_update_pbc, cell_v, n_loc_size, &
556 250380 : 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 250380 : CALL timeset(routineN, handle)
578 1001520 : rij(:) = r_last_update_pbc(jparticle)%r(:) - r_last_update_pbc(iparticle)%r(:) + cell_v
579 1001520 : drij = SQRT(DOT_PRODUCT(rij, rij))
580 1001520 : rij_hat(:) = rij(:)/drij
581 :
582 250380 : fac = -0.5_dp
583 250380 : 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 250380 : 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 250380 : f_A = ter_f_A(tersoff, drij)
586 250380 : f_A_d = ter_f_A_d(tersoff, drij)
587 250380 : f_C = ter_f_C(tersoff, drij)
588 250380 : f_C_d = ter_f_C_d(tersoff, drij)
589 250380 : f_R = ter_f_R(tersoff, drij)
590 250380 : 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 250380 : f_R1 = f_C_d*f_R*fac
595 250380 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_R1*rij(1)
596 250380 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_R1*rij(2)
597 250380 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_R1*rij(3)
598 250380 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_R1*rij(1)
599 250380 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_R1*rij(2)
600 250380 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_R1*rij(3)
601 :
602 250380 : IF (use_virial) THEN
603 74582 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_R1*rij(1)*rij(1)
604 74582 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_R1*rij(1)*rij(2)
605 74582 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_R1*rij(1)*rij(3)
606 74582 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_R1*rij(2)*rij(1)
607 74582 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_R1*rij(2)*rij(2)
608 74582 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_R1*rij(2)*rij(3)
609 74582 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_R1*rij(3)*rij(1)
610 74582 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_R1*rij(3)*rij(2)
611 74582 : pv_nonbond(3, 3) = pv_nonbond(3, 3) - f_R1*rij(3)*rij(3)
612 : END IF
613 :
614 250380 : f_R2 = f_C*f_R_d*fac
615 250380 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_R2*rij(1)
616 250380 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_R2*rij(2)
617 250380 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_R2*rij(3)
618 250380 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_R2*rij(1)
619 250380 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_R2*rij(2)
620 250380 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_R2*rij(3)
621 :
622 250380 : IF (use_virial) THEN
623 74582 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_R2*rij(1)*rij(1)
624 74582 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_R2*rij(1)*rij(2)
625 74582 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_R2*rij(1)*rij(3)
626 74582 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_R2*rij(2)*rij(1)
627 74582 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_R2*rij(2)*rij(2)
628 74582 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_R2*rij(2)*rij(3)
629 74582 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_R2*rij(3)*rij(1)
630 74582 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_R2*rij(3)*rij(2)
631 74582 : 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 250380 : f_A1 = f_C_d*b_ij*f_A*fac
636 250380 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A1*rij(1)
637 250380 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A1*rij(2)
638 250380 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A1*rij(3)
639 250380 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A1*rij(1)
640 250380 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A1*rij(2)
641 250380 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A1*rij(3)
642 :
643 250380 : IF (use_virial) THEN
644 74582 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_A1*rij(1)*rij(1)
645 74582 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_A1*rij(1)*rij(2)
646 74582 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_A1*rij(1)*rij(3)
647 74582 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_A1*rij(2)*rij(1)
648 74582 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_A1*rij(2)*rij(2)
649 74582 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_A1*rij(2)*rij(3)
650 74582 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_A1*rij(3)*rij(1)
651 74582 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_A1*rij(3)*rij(2)
652 74582 : 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 250380 : f_A2 = f_C*b_ij*f_A_d*fac
657 250380 : f_nonbond(1, iparticle) = f_nonbond(1, iparticle) + f_A2*rij(1)
658 250380 : f_nonbond(2, iparticle) = f_nonbond(2, iparticle) + f_A2*rij(2)
659 250380 : f_nonbond(3, iparticle) = f_nonbond(3, iparticle) + f_A2*rij(3)
660 250380 : f_nonbond(1, jparticle) = f_nonbond(1, jparticle) - f_A2*rij(1)
661 250380 : f_nonbond(2, jparticle) = f_nonbond(2, jparticle) - f_A2*rij(2)
662 250380 : f_nonbond(3, jparticle) = f_nonbond(3, jparticle) - f_A2*rij(3)
663 :
664 250380 : IF (use_virial) THEN
665 74582 : pv_nonbond(1, 1) = pv_nonbond(1, 1) - f_A2*rij(1)*rij(1)
666 74582 : pv_nonbond(1, 2) = pv_nonbond(1, 2) - f_A2*rij(1)*rij(2)
667 74582 : pv_nonbond(1, 3) = pv_nonbond(1, 3) - f_A2*rij(1)*rij(3)
668 74582 : pv_nonbond(2, 1) = pv_nonbond(2, 1) - f_A2*rij(2)*rij(1)
669 74582 : pv_nonbond(2, 2) = pv_nonbond(2, 2) - f_A2*rij(2)*rij(2)
670 74582 : pv_nonbond(2, 3) = pv_nonbond(2, 3) - f_A2*rij(2)*rij(3)
671 74582 : pv_nonbond(3, 1) = pv_nonbond(3, 1) - f_A2*rij(3)*rij(1)
672 74582 : pv_nonbond(3, 2) = pv_nonbond(3, 2) - f_A2*rij(3)*rij(2)
673 74582 : 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 250380 : 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 250380 : n_loc_size, full_loc_list, loc_cell_v, cell_v, rcutsq, use_virial)
680 250380 : CALL timestop(handle)
681 250380 : 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 5328 : 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 5328 : INTEGER, DIMENSION(:), POINTER :: work_list, work_list2
709 5328 : INTEGER, DIMENSION(:, :), POINTER :: list
710 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi
711 5328 : 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 5328 : CPASSERT(.NOT. ASSOCIATED(glob_loc_list_a))
717 5328 : CPASSERT(.NOT. ASSOCIATED(glob_cell_v))
718 5328 : CALL timeset(routineN, handle)
719 5328 : npairs_tot = 0
720 5328 : nkinds = SIZE(potparm%pot, 1)
721 202104 : DO ilist = 1, nonbonded%nlists
722 196776 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
723 196776 : npairs = neighbor_kind_pair%npairs
724 196776 : IF (npairs == 0) CYCLE
725 136598 : Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
726 66756 : istart = neighbor_kind_pair%grp_kind_start(igrp)
727 66756 : iend = neighbor_kind_pair%grp_kind_end(igrp)
728 66756 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
729 66756 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
730 66756 : pot => potparm%pot(ikind, jkind)%pot
731 66756 : npairs = iend - istart + 1
732 66756 : IF (pot%no_mb) CYCLE Kind_Group_Loop1
733 330148 : DO i = 1, SIZE(pot%type)
734 133444 : IF (pot%type(i) == tersoff_type) npairs_tot = npairs_tot + npairs
735 : END DO
736 : END DO Kind_Group_Loop1
737 : END DO
738 15984 : ALLOCATE (work_list(npairs_tot))
739 10656 : ALLOCATE (work_list2(npairs_tot))
740 15984 : ALLOCATE (glob_loc_list(2, npairs_tot))
741 15984 : ALLOCATE (glob_cell_v(3, npairs_tot))
742 : ! Fill arrays with data
743 5328 : npairs_tot = 0
744 202104 : DO ilist = 1, nonbonded%nlists
745 196776 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
746 196776 : npairs = neighbor_kind_pair%npairs
747 196776 : IF (npairs == 0) CYCLE
748 136598 : Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
749 66756 : istart = neighbor_kind_pair%grp_kind_start(igrp)
750 66756 : iend = neighbor_kind_pair%grp_kind_end(igrp)
751 66756 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
752 66756 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
753 66756 : list => neighbor_kind_pair%list
754 267024 : cvi = neighbor_kind_pair%cell_vector
755 66756 : pot => potparm%pot(ikind, jkind)%pot
756 66756 : npairs = iend - istart + 1
757 66756 : IF (pot%no_mb) CYCLE Kind_Group_Loop2
758 866892 : cell_v = MATMUL(cell%hmat, cvi)
759 330148 : DO i = 1, SIZE(pot%type)
760 : ! TERSOFF
761 133444 : IF (pot%type(i) == tersoff_type) THEN
762 9725342 : DO ipair = 1, npairs
763 57952008 : glob_loc_list(:, npairs_tot + ipair) = list(:, istart - 1 + ipair)
764 38701346 : glob_cell_v(1:3, npairs_tot + ipair) = cell_v(1:3)
765 : END DO
766 66674 : 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 5328 : CALL sort(glob_loc_list(1, :), npairs_tot, work_list)
773 9663996 : DO ipair = 1, npairs_tot
774 9663996 : work_list2(ipair) = glob_loc_list(2, work_list(ipair))
775 : END DO
776 19327992 : glob_loc_list(2, :) = work_list2
777 5328 : DEALLOCATE (work_list2)
778 15984 : ALLOCATE (rwork_list(3, npairs_tot))
779 9663996 : DO ipair = 1, npairs_tot
780 77274672 : rwork_list(:, ipair) = glob_cell_v(:, work_list(ipair))
781 : END DO
782 77280000 : glob_cell_v = rwork_list
783 5328 : DEALLOCATE (rwork_list)
784 5328 : DEALLOCATE (work_list)
785 15984 : ALLOCATE (glob_loc_list_a(npairs_tot))
786 19327992 : glob_loc_list_a = glob_loc_list(1, :)
787 5328 : CALL timestop(handle)
788 10656 : 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 5328 : 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 5328 : IF (ASSOCIATED(glob_loc_list)) THEN
805 5328 : DEALLOCATE (glob_loc_list)
806 : END IF
807 5328 : IF (ASSOCIATED(glob_loc_list_a)) THEN
808 5328 : DEALLOCATE (glob_loc_list_a)
809 : END IF
810 5328 : IF (ASSOCIATED(glob_cell_v)) THEN
811 5328 : DEALLOCATE (glob_cell_v)
812 : END IF
813 :
814 5328 : END SUBROUTINE destroy_tersoff_arrays
815 :
816 : END MODULE manybody_tersoff
817 :
|
