(git:374b731)
Loading...
Searching...
No Matches
manybody_potential.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 and general "lifting" of manybody_potential module
11!> 12.2007 [tlaino] - Splitting manybody module : In this module we should only
12!> keep the main routines for computing energy and forces of
13!> manybody potentials. Each potential should have his own module!
14!> \author CJM, I-Feng W. Kuo, Teodoro Laino
15! **************************************************************************************************
16MODULE manybody_potential
17
18 USE atomic_kind_types, ONLY: atomic_kind_type
19 USE cell_types, ONLY: cell_type
20 USE distribution_1d_types, ONLY: distribution_1d_type
21 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,&
22 neighbor_kind_pairs_type
23 USE fist_nonbond_env_types, ONLY: eam_type,&
24 fist_nonbond_env_get,&
25 fist_nonbond_env_type,&
26 pos_type
27 USE input_section_types, ONLY: section_vals_type
28 USE kinds, ONLY: dp
29 USE manybody_allegro, ONLY: allegro_add_force_virial,&
30 allegro_energy_store_force_virial,&
31 destroy_allegro_arrays,&
32 setup_allegro_arrays
33 USE manybody_deepmd, ONLY: deepmd_add_force_virial,&
34 deepmd_energy_store_force_virial
35 USE manybody_eam, ONLY: get_force_eam
36 USE manybody_gal, ONLY: destroy_gal_arrays,&
37 gal_energy,&
38 gal_forces,&
39 setup_gal_arrays
40 USE manybody_gal21, ONLY: destroy_gal21_arrays,&
41 gal21_energy,&
42 gal21_forces,&
43 setup_gal21_arrays
44 USE manybody_nequip, ONLY: destroy_nequip_arrays,&
45 nequip_add_force_virial,&
46 nequip_energy_store_force_virial,&
47 setup_nequip_arrays
48 USE manybody_quip, ONLY: quip_add_force_virial,&
49 quip_energy_store_force_virial
50 USE manybody_siepmann, ONLY: destroy_siepmann_arrays,&
51 print_nr_ions_siepmann,&
52 setup_siepmann_arrays,&
53 siepmann_energy,&
54 siepmann_forces_v2,&
55 siepmann_forces_v3
56 USE manybody_tersoff, ONLY: destroy_tersoff_arrays,&
57 setup_tersoff_arrays,&
58 tersoff_energy,&
59 tersoff_forces
60 USE message_passing, ONLY: mp_para_env_type
61 USE pair_potential_types, ONLY: &
62 allegro_pot_type, allegro_type, deepmd_type, ea_type, eam_pot_type, gal21_pot_type, &
63 gal21_type, gal_pot_type, gal_type, nequip_pot_type, nequip_type, pair_potential_pp_type, &
64 pair_potential_single_type, quip_type, siepmann_pot_type, siepmann_type, tersoff_pot_type, &
65 tersoff_type
66 USE particle_types, ONLY: particle_type
67 USE util, ONLY: sort
68#include "./base/base_uses.f90"
69
70 IMPLICIT NONE
71
72 PRIVATE
73 PUBLIC :: energy_manybody
74 PUBLIC :: force_nonbond_manybody
75 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_potential'
76
77CONTAINS
78
79! **************************************************************************************************
80!> \brief computes the embedding contribution to the energy
81!> \param fist_nonbond_env ...
82!> \param atomic_kind_set ...
83!> \param local_particles ...
84!> \param particle_set ...
85!> \param cell ...
86!> \param pot_manybody ...
87!> \param para_env ...
88!> \param mm_section ...
89!> \par History
90!> tlaino [2007] - New algorithm for tersoff potential
91!> \author CJM, I-Feng W. Kuo, Teodoro Laino
92! **************************************************************************************************
93 SUBROUTINE energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, &
94 particle_set, cell, pot_manybody, para_env, mm_section)
95
96 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
97 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
98 TYPE(distribution_1d_type), POINTER :: local_particles
99 TYPE(particle_type), POINTER :: particle_set(:)
100 TYPE(cell_type), POINTER :: cell
101 REAL(dp), INTENT(INOUT) :: pot_manybody
102 TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
103 TYPE(section_vals_type), POINTER :: mm_section
104
105 CHARACTER(LEN=*), PARAMETER :: routinen = 'energy_manybody'
106
107 INTEGER :: atom_a, atom_b, handle, i, iend, ifirst, igrp, ikind, ilast, ilist, indexa, &
108 ipair, iparticle, iparticle_local, istart, iunique, jkind, junique, mpair, nkinds, &
109 nloc_size, npairs, nparticle, nparticle_local, nr_h3o, nr_o, nr_oh, nunique
110 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, unique_list_a, work_list
111 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list
112 LOGICAL :: any_allegro, any_deepmd, any_gal, &
113 any_gal21, any_nequip, any_quip, &
114 any_siepmann, any_tersoff
115 REAL(kind=dp) :: drij, embed, pot_allegro, pot_deepmd, &
116 pot_loc, pot_nequip, pot_quip, qr, &
117 rab2_max, rij(3)
118 REAL(kind=dp), DIMENSION(3) :: cell_v, cvi
119 REAL(kind=dp), DIMENSION(:, :), POINTER :: glob_cell_v
120 REAL(kind=dp), POINTER :: fembed(:)
121 TYPE(allegro_pot_type), POINTER :: allegro
122 TYPE(eam_pot_type), POINTER :: eam
123 TYPE(eam_type), DIMENSION(:), POINTER :: eam_data
124 TYPE(fist_neighbor_type), POINTER :: nonbonded
125 TYPE(gal21_pot_type), POINTER :: gal21
126 TYPE(gal_pot_type), POINTER :: gal
127 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
128 TYPE(nequip_pot_type), POINTER :: nequip
129 TYPE(pair_potential_pp_type), POINTER :: potparm
130 TYPE(pair_potential_single_type), POINTER :: pot
131 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
132 TYPE(siepmann_pot_type), POINTER :: siepmann
133 TYPE(tersoff_pot_type), POINTER :: tersoff
134
135 NULLIFY (eam, siepmann, tersoff, gal, gal21)
136 any_tersoff = .false.
137 any_siepmann = .false.
138 any_gal = .false.
139 any_gal21 = .false.
140 any_quip = .false.
141 any_nequip = .false.
142 any_allegro = .false.
143 any_deepmd = .false.
144 CALL timeset(routinen, handle)
145 CALL fist_nonbond_env_get(fist_nonbond_env, r_last_update_pbc=r_last_update_pbc, &
146 potparm=potparm, eam_data=eam_data)
147 ! EAM requires a single loop
148 DO ikind = 1, SIZE(atomic_kind_set)
149 pot => potparm%pot(ikind, ikind)%pot
150 DO i = 1, SIZE(pot%type)
151 IF (pot%type(i) /= ea_type) cycle
152 eam => pot%set(i)%eam
153 nparticle = SIZE(particle_set)
154 ALLOCATE (fembed(nparticle))
155 fembed(:) = 0._dp
156 cpassert(ASSOCIATED(eam_data))
157 ! computation of embedding function and energy
158 nparticle_local = local_particles%n_el(ikind)
159 DO iparticle_local = 1, nparticle_local
160 iparticle = local_particles%list(ikind)%array(iparticle_local)
161 indexa = int(eam_data(iparticle)%rho/eam%drhoar) + 1
162 IF (indexa > eam%npoints - 1) indexa = eam%npoints - 1
163 qr = eam_data(iparticle)%rho - eam%rhoval(indexa)
164
165 embed = eam%frho(indexa) + qr*eam%frhop(indexa)
166 fembed(iparticle) = eam%frhop(indexa) + qr*(eam%frhop(indexa + 1) - eam%frhop(indexa))/eam%drhoar
167
168 pot_manybody = pot_manybody + embed
169 END DO
170 ! communicate data
171 CALL para_env%sum(fembed)
172 DO iparticle = 1, nparticle
173 IF (particle_set(iparticle)%atomic_kind%kind_number == ikind) THEN
174 eam_data(iparticle)%f_embed = fembed(iparticle)
175 END IF
176 END DO
177
178 DEALLOCATE (fembed)
179 END DO
180 END DO
181 ! Other manybody potential
182 DO ikind = 1, SIZE(atomic_kind_set)
183 DO jkind = ikind, SIZE(atomic_kind_set)
184 pot => potparm%pot(ikind, jkind)%pot
185 any_tersoff = any_tersoff .OR. any(pot%type == tersoff_type)
186 any_quip = any_quip .OR. any(pot%type == quip_type)
187 any_nequip = any_nequip .OR. any(pot%type == nequip_type)
188 any_allegro = any_allegro .OR. any(pot%type == allegro_type)
189 any_deepmd = any_deepmd .OR. any(pot%type == deepmd_type)
190 any_siepmann = any_siepmann .OR. any(pot%type == siepmann_type)
191 any_gal = any_gal .OR. any(pot%type == gal_type)
192 any_gal21 = any_gal21 .OR. any(pot%type == gal21_type)
193 END DO
194 END DO
195 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, natom_types=nkinds)
196 ! QUIP
197 IF (any_quip) THEN
198 CALL quip_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, &
199 fist_nonbond_env, pot_quip, para_env)
200 pot_manybody = pot_manybody + pot_quip
201 END IF
202 ! NEQUIP
203 IF (any_nequip) THEN
204 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
205 CALL setup_nequip_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
206 CALL nequip_energy_store_force_virial(nonbonded, particle_set, cell, atomic_kind_set, potparm, &
207 nequip, glob_loc_list_a, r_last_update_pbc, pot_nequip, &
208 fist_nonbond_env, para_env)
209 pot_manybody = pot_manybody + pot_nequip
210 CALL destroy_nequip_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
211 END IF
212 ! ALLEGRO
213 IF (any_allegro) THEN
214 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a, unique_list_a)
215 CALL setup_allegro_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, &
216 unique_list_a, cell)
217 CALL allegro_energy_store_force_virial(nonbonded, particle_set, cell, atomic_kind_set, potparm, &
218 allegro, glob_loc_list_a, r_last_update_pbc, pot_allegro, &
219 fist_nonbond_env, unique_list_a)
220 pot_manybody = pot_manybody + pot_allegro
221 CALL destroy_allegro_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a, unique_list_a)
222 END IF
223 ! DEEPMD
224 IF (any_deepmd) THEN
225 CALL deepmd_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, &
226 fist_nonbond_env, pot_deepmd, para_env)
227 pot_manybody = pot_manybody + pot_deepmd
228 END IF
229
230 ! TERSOFF
231 IF (any_tersoff) THEN
232 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
233 CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
234 DO ilist = 1, nonbonded%nlists
235 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
236 npairs = neighbor_kind_pair%npairs
237 IF (npairs == 0) cycle
238 kind_group_loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
239 istart = neighbor_kind_pair%grp_kind_start(igrp)
240 iend = neighbor_kind_pair%grp_kind_end(igrp)
241 ikind = neighbor_kind_pair%ij_kind(1, igrp)
242 jkind = neighbor_kind_pair%ij_kind(2, igrp)
243 list => neighbor_kind_pair%list
244 cvi = neighbor_kind_pair%cell_vector
245 pot => potparm%pot(ikind, jkind)%pot
246 DO i = 1, SIZE(pot%type)
247 IF (pot%type(i) /= tersoff_type) cycle
248 rab2_max = pot%set(i)%tersoff%rcutsq
249 cell_v = matmul(cell%hmat, cvi)
250 pot => potparm%pot(ikind, jkind)%pot
251 tersoff => pot%set(i)%tersoff
252 npairs = iend - istart + 1
253 IF (npairs /= 0) THEN
254 ALLOCATE (sort_list(2, npairs), work_list(npairs))
255 sort_list = list(:, istart:iend)
256 ! Sort the list of neighbors, this increases the efficiency for single
257 ! potential contributions
258 CALL sort(sort_list(1, :), npairs, work_list)
259 DO ipair = 1, npairs
260 work_list(ipair) = sort_list(2, work_list(ipair))
261 END DO
262 sort_list(2, :) = work_list
263 ! find number of unique elements of array index 1
264 nunique = 1
265 DO ipair = 1, npairs - 1
266 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
267 END DO
268 ipair = 1
269 junique = sort_list(1, ipair)
270 ifirst = 1
271 DO iunique = 1, nunique
272 atom_a = junique
273 IF (glob_loc_list_a(ifirst) > atom_a) cycle
274 DO mpair = ifirst, SIZE(glob_loc_list_a)
275 IF (glob_loc_list_a(mpair) == atom_a) EXIT
276 END DO
277 ifirst = mpair
278 DO mpair = ifirst, SIZE(glob_loc_list_a)
279 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
280 END DO
281 ilast = mpair - 1
282 nloc_size = 0
283 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
284 DO WHILE (ipair <= npairs)
285 IF (sort_list(1, ipair) /= junique) EXIT
286 atom_b = sort_list(2, ipair)
287 ! Energy terms
288 pot_loc = 0.0_dp
289 rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
290 drij = dot_product(rij, rij)
291 ipair = ipair + 1
292 IF (drij > rab2_max) cycle
293 drij = sqrt(drij)
294 CALL tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
295 glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), cell_v, drij)
296 pot_manybody = pot_manybody + 0.5_dp*pot_loc
297 END DO
298 ifirst = ilast + 1
299 IF (ipair <= npairs) junique = sort_list(1, ipair)
300 END DO
301 DEALLOCATE (sort_list, work_list)
302 END IF
303 END DO
304 END DO kind_group_loop
305 END DO
306 CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
307 END IF
308
309 !SIEPMANN POTENTIAL
310 IF (any_siepmann) THEN
311 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
312 nr_oh = 0
313 nr_h3o = 0
314 nr_o = 0
315 CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
316 DO ilist = 1, nonbonded%nlists
317 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
318 npairs = neighbor_kind_pair%npairs
319 IF (npairs == 0) cycle
320 kind_group_loop_2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
321 istart = neighbor_kind_pair%grp_kind_start(igrp)
322 iend = neighbor_kind_pair%grp_kind_end(igrp)
323 ikind = neighbor_kind_pair%ij_kind(1, igrp)
324 jkind = neighbor_kind_pair%ij_kind(2, igrp)
325 list => neighbor_kind_pair%list
326 cvi = neighbor_kind_pair%cell_vector
327 pot => potparm%pot(ikind, jkind)%pot
328 DO i = 1, SIZE(pot%type)
329 IF (pot%type(i) /= siepmann_type) cycle
330 rab2_max = pot%set(i)%siepmann%rcutsq
331 cell_v = matmul(cell%hmat, cvi)
332 pot => potparm%pot(ikind, jkind)%pot
333 siepmann => pot%set(i)%siepmann
334 npairs = iend - istart + 1
335 IF (npairs /= 0) THEN
336 ALLOCATE (sort_list(2, npairs), work_list(npairs))
337 sort_list = list(:, istart:iend)
338 ! Sort the list of neighbors, this increases the efficiency for single
339 ! potential contributions
340 CALL sort(sort_list(1, :), npairs, work_list)
341 DO ipair = 1, npairs
342 work_list(ipair) = sort_list(2, work_list(ipair))
343 END DO
344 sort_list(2, :) = work_list
345 ! find number of unique elements of array index 1
346 nunique = 1
347 DO ipair = 1, npairs - 1
348 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
349 END DO
350 ipair = 1
351 junique = sort_list(1, ipair)
352 ifirst = 1
353 DO iunique = 1, nunique
354 atom_a = junique
355 IF (glob_loc_list_a(ifirst) > atom_a) cycle
356 DO mpair = ifirst, SIZE(glob_loc_list_a)
357 IF (glob_loc_list_a(mpair) == atom_a) EXIT
358 END DO
359 ifirst = mpair
360 DO mpair = ifirst, SIZE(glob_loc_list_a)
361 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
362 END DO
363 ilast = mpair - 1
364 nloc_size = 0
365 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
366 DO WHILE (ipair <= npairs)
367 IF (sort_list(1, ipair) /= junique) EXIT
368 atom_b = sort_list(2, ipair)
369 ! Energy terms
370 pot_loc = 0.0_dp
371 rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
372 drij = dot_product(rij, rij)
373 ipair = ipair + 1
374 IF (drij > rab2_max) cycle
375 drij = sqrt(drij)
376 CALL siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size, &
377 glob_loc_list(:, ifirst:ilast), cell_v, cell, drij, &
378 particle_set, nr_oh, nr_h3o, nr_o)
379 pot_manybody = pot_manybody + pot_loc
380 END DO
381 ifirst = ilast + 1
382 IF (ipair <= npairs) junique = sort_list(1, ipair)
383 END DO
384 DEALLOCATE (sort_list, work_list)
385 END IF
386 END DO
387 END DO kind_group_loop_2
388 END DO
389 CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
390 CALL print_nr_ions_siepmann(nr_oh, mm_section, para_env, print_oh=.true., &
391 print_h3o=.false., print_o=.false.)
392 CALL print_nr_ions_siepmann(nr_h3o, mm_section, para_env, print_oh=.false., &
393 print_h3o=.true., print_o=.false.)
394 CALL print_nr_ions_siepmann(nr_o, mm_section, para_env, print_oh=.false., &
395 print_h3o=.false., print_o=.true.)
396 END IF
397
398 !GAL19 POTENTIAL
399 IF (any_gal) THEN
400 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
401 CALL setup_gal_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
402 DO ilist = 1, nonbonded%nlists
403 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
404 npairs = neighbor_kind_pair%npairs
405 IF (npairs == 0) cycle
406 kind_group_loop_3: DO igrp = 1, neighbor_kind_pair%ngrp_kind
407 istart = neighbor_kind_pair%grp_kind_start(igrp)
408 iend = neighbor_kind_pair%grp_kind_end(igrp)
409 ikind = neighbor_kind_pair%ij_kind(1, igrp)
410 jkind = neighbor_kind_pair%ij_kind(2, igrp)
411 list => neighbor_kind_pair%list
412 cvi = neighbor_kind_pair%cell_vector
413 pot => potparm%pot(ikind, jkind)%pot
414 DO i = 1, SIZE(pot%type)
415 IF (pot%type(i) /= gal_type) cycle
416 rab2_max = pot%set(i)%gal%rcutsq
417 cell_v = matmul(cell%hmat, cvi)
418 pot => potparm%pot(ikind, jkind)%pot
419 gal => pot%set(i)%gal
420 npairs = iend - istart + 1
421 IF (npairs /= 0) THEN
422 ALLOCATE (sort_list(2, npairs), work_list(npairs))
423 sort_list = list(:, istart:iend)
424 ! Sort the list of neighbors, this increases the efficiency for single
425 ! potential contributions
426 CALL sort(sort_list(1, :), npairs, work_list)
427 DO ipair = 1, npairs
428 work_list(ipair) = sort_list(2, work_list(ipair))
429 END DO
430 sort_list(2, :) = work_list
431 ! find number of unique elements of array index 1
432 nunique = 1
433 DO ipair = 1, npairs - 1
434 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
435 END DO
436 ipair = 1
437 junique = sort_list(1, ipair)
438 ifirst = 1
439 DO iunique = 1, nunique
440 atom_a = junique
441 IF (glob_loc_list_a(ifirst) > atom_a) cycle
442 DO mpair = ifirst, SIZE(glob_loc_list_a)
443 IF (glob_loc_list_a(mpair) == atom_a) EXIT
444 END DO
445 ifirst = mpair
446 DO mpair = ifirst, SIZE(glob_loc_list_a)
447 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
448 END DO
449 ilast = mpair - 1
450 nloc_size = 0
451 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
452 DO WHILE (ipair <= npairs)
453 IF (sort_list(1, ipair) /= junique) EXIT
454 atom_b = sort_list(2, ipair)
455 ! Energy terms
456 pot_loc = 0.0_dp
457 rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
458 drij = dot_product(rij, rij)
459 ipair = ipair + 1
460 IF (drij > rab2_max) cycle
461 drij = sqrt(drij)
462 CALL gal_energy(pot_loc, gal, r_last_update_pbc, atom_a, atom_b, &
463 cell, particle_set, mm_section)
464
465 pot_manybody = pot_manybody + pot_loc
466 END DO
467 ifirst = ilast + 1
468 IF (ipair <= npairs) junique = sort_list(1, ipair)
469 END DO
470 DEALLOCATE (sort_list, work_list)
471 END IF
472 END DO
473 END DO kind_group_loop_3
474 END DO
475 CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
476 END IF
477
478 !GAL21 POTENTIAL
479 IF (any_gal21) THEN
480 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
481 CALL setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
482 DO ilist = 1, nonbonded%nlists
483 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
484 npairs = neighbor_kind_pair%npairs
485 IF (npairs == 0) cycle
486 kind_group_loop_5: DO igrp = 1, neighbor_kind_pair%ngrp_kind
487 istart = neighbor_kind_pair%grp_kind_start(igrp)
488 iend = neighbor_kind_pair%grp_kind_end(igrp)
489 ikind = neighbor_kind_pair%ij_kind(1, igrp)
490 jkind = neighbor_kind_pair%ij_kind(2, igrp)
491 list => neighbor_kind_pair%list
492 cvi = neighbor_kind_pair%cell_vector
493 pot => potparm%pot(ikind, jkind)%pot
494 DO i = 1, SIZE(pot%type)
495 IF (pot%type(i) /= gal21_type) cycle
496 rab2_max = pot%set(i)%gal21%rcutsq
497 cell_v = matmul(cell%hmat, cvi)
498 pot => potparm%pot(ikind, jkind)%pot
499 gal21 => pot%set(i)%gal21
500 npairs = iend - istart + 1
501 IF (npairs /= 0) THEN
502 ALLOCATE (sort_list(2, npairs), work_list(npairs))
503 sort_list = list(:, istart:iend)
504 ! Sort the list of neighbors, this increases the efficiency for single
505 ! potential contributions
506 CALL sort(sort_list(1, :), npairs, work_list)
507 DO ipair = 1, npairs
508 work_list(ipair) = sort_list(2, work_list(ipair))
509 END DO
510 sort_list(2, :) = work_list
511 ! find number of unique elements of array index 1
512 nunique = 1
513 DO ipair = 1, npairs - 1
514 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
515 END DO
516 ipair = 1
517 junique = sort_list(1, ipair)
518 ifirst = 1
519 DO iunique = 1, nunique
520 atom_a = junique
521 IF (glob_loc_list_a(ifirst) > atom_a) cycle
522 DO mpair = ifirst, SIZE(glob_loc_list_a)
523 IF (glob_loc_list_a(mpair) == atom_a) EXIT
524 END DO
525 ifirst = mpair
526 DO mpair = ifirst, SIZE(glob_loc_list_a)
527 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
528 END DO
529 ilast = mpair - 1
530 nloc_size = 0
531 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
532 DO WHILE (ipair <= npairs)
533 IF (sort_list(1, ipair) /= junique) EXIT
534 atom_b = sort_list(2, ipair)
535 ! Energy terms
536 pot_loc = 0.0_dp
537 rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
538 drij = dot_product(rij, rij)
539 ipair = ipair + 1
540 IF (drij > rab2_max) cycle
541 drij = sqrt(drij)
542 CALL gal21_energy(pot_loc, gal21, r_last_update_pbc, atom_a, atom_b, &
543 cell, particle_set, mm_section)
544
545 pot_manybody = pot_manybody + pot_loc
546 END DO
547 ifirst = ilast + 1
548 IF (ipair <= npairs) junique = sort_list(1, ipair)
549 END DO
550 DEALLOCATE (sort_list, work_list)
551 END IF
552 END DO
553 END DO kind_group_loop_5
554 END DO
555 CALL destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
556 END IF
557
558 CALL timestop(handle)
559 END SUBROUTINE energy_manybody
560
561! **************************************************************************************************
562!> \brief ...
563!> \param fist_nonbond_env ...
564!> \param particle_set ...
565!> \param cell ...
566!> \param f_nonbond ...
567!> \param pv_nonbond ...
568!> \param use_virial ...
569!> \par History
570!> Fast implementation of the tersoff potential - [tlaino] 2007
571!> \author I-Feng W. Kuo, Teodoro Laino
572! **************************************************************************************************
573 SUBROUTINE force_nonbond_manybody(fist_nonbond_env, particle_set, cell, &
574 f_nonbond, pv_nonbond, use_virial)
575
576 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
577 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
578 TYPE(cell_type), POINTER :: cell
579 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
580 LOGICAL, INTENT(IN) :: use_virial
581
582 CHARACTER(LEN=*), PARAMETER :: routinen = 'force_nonbond_manybody'
583
584 INTEGER :: atom_a, atom_b, handle, i, i_a, i_b, iend, ifirst, igrp, ikind, ilast, ilist, &
585 ipair, istart, iunique, jkind, junique, kind_a, kind_b, mpair, nkinds, nloc_size, npairs, &
586 nunique
587 INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index
588 INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list
589 INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list
590 LOGICAL :: any_allegro, any_deepmd, any_gal, &
591 any_gal21, any_nequip, any_quip, &
592 any_siepmann, any_tersoff
593 REAL(kind=dp) :: f_eam, fac, fr(3), ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, &
594 ptens31, ptens32, ptens33, rab(3), rab2, rab2_max, rtmp(3)
595 REAL(kind=dp), DIMENSION(3) :: cell_v, cvi
596 REAL(kind=dp), DIMENSION(:, :), POINTER :: glob_cell_v
597 TYPE(eam_pot_type), POINTER :: eam_a, eam_b
598 TYPE(eam_type), DIMENSION(:), POINTER :: eam_data
599 TYPE(fist_neighbor_type), POINTER :: nonbonded
600 TYPE(gal21_pot_type), POINTER :: gal21
601 TYPE(gal_pot_type), POINTER :: gal
602 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
603 TYPE(pair_potential_pp_type), POINTER :: potparm
604 TYPE(pair_potential_single_type), POINTER :: pot
605 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
606 TYPE(siepmann_pot_type), POINTER :: siepmann
607 TYPE(tersoff_pot_type), POINTER :: tersoff
608
609 any_tersoff = .false.
610 any_quip = .false.
611 any_nequip = .false.
612 any_allegro = .false.
613 any_siepmann = .false.
614 any_deepmd = .false.
615 any_gal = .false.
616 any_gal21 = .false.
617 CALL timeset(routinen, handle)
618 NULLIFY (eam_a, eam_b, tersoff, siepmann, gal, gal21)
619
620 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, potparm=potparm, &
621 natom_types=nkinds, eam_data=eam_data, r_last_update_pbc=r_last_update_pbc)
622
623 ! Initializing the potential energy, pressure tensor and force
624 IF (use_virial) THEN
625 ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
626 ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
627 ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
628 END IF
629
630 nkinds = SIZE(potparm%pot, 1)
631 ALLOCATE (eam_kinds_index(nkinds, nkinds))
632 eam_kinds_index = -1
633 DO ikind = 1, nkinds
634 DO jkind = ikind, nkinds
635 DO i = 1, SIZE(potparm%pot(ikind, jkind)%pot%type)
636 IF (potparm%pot(ikind, jkind)%pot%type(i) == ea_type) THEN
637 ! At the moment we allow only 1 EAM per each kinds pair..
638 cpassert(eam_kinds_index(ikind, jkind) == -1)
639 cpassert(eam_kinds_index(jkind, ikind) == -1)
640 eam_kinds_index(ikind, jkind) = i
641 eam_kinds_index(jkind, ikind) = i
642 END IF
643 END DO
644 END DO
645 END DO
646 DO ikind = 1, nkinds
647 DO jkind = ikind, nkinds
648 any_deepmd = any_deepmd .OR. any(potparm%pot(ikind, jkind)%pot%type == deepmd_type)
649 END DO
650 END DO
651 ! DEEPMD
652 IF (any_deepmd) &
653 CALL deepmd_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
654
655 ! QUIP
656 IF (use_virial) &
657 CALL quip_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond)
658
659 ! starting the force loop
660 DO ilist = 1, nonbonded%nlists
661 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
662 npairs = neighbor_kind_pair%npairs
663 IF (npairs == 0) cycle
664 kind_group_loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
665 istart = neighbor_kind_pair%grp_kind_start(igrp)
666 iend = neighbor_kind_pair%grp_kind_end(igrp)
667 ikind = neighbor_kind_pair%ij_kind(1, igrp)
668 jkind = neighbor_kind_pair%ij_kind(2, igrp)
669 list => neighbor_kind_pair%list
670 cvi = neighbor_kind_pair%cell_vector
671 pot => potparm%pot(ikind, jkind)%pot
672 IF (pot%no_mb) cycle
673 rab2_max = pot%rcutsq
674 cell_v = matmul(cell%hmat, cvi)
675 any_tersoff = any_tersoff .OR. any(pot%type == tersoff_type)
676 any_siepmann = any_siepmann .OR. any(pot%type == siepmann_type)
677 any_deepmd = any_deepmd .OR. any(pot%type == deepmd_type)
678 any_gal = any_gal .OR. any(pot%type == gal_type)
679 any_gal21 = any_gal21 .OR. any(pot%type == gal21_type)
680 any_nequip = any_nequip .OR. any(pot%type == nequip_type)
681 any_allegro = any_allegro .OR. any(pot%type == allegro_type)
682 i = eam_kinds_index(ikind, jkind)
683 IF (i == -1) cycle
684 ! EAM
685 cpassert(ASSOCIATED(eam_data))
686 DO ipair = istart, iend
687 atom_a = list(1, ipair)
688 atom_b = list(2, ipair)
689 fac = 1.0_dp
690 IF (atom_a == atom_b) fac = 0.5_dp
691 kind_a = particle_set(atom_a)%atomic_kind%kind_number
692 kind_b = particle_set(atom_b)%atomic_kind%kind_number
693 i_a = eam_kinds_index(kind_a, kind_a)
694 i_b = eam_kinds_index(kind_b, kind_b)
695 eam_a => potparm%pot(kind_a, kind_a)%pot%set(i_a)%eam
696 eam_b => potparm%pot(kind_b, kind_b)%pot%set(i_b)%eam
697
698 !set this outside the potential type in case need multiple potentials
699 !Do everything necessary for EAM here
700 rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
701 rab = rab + cell_v
702 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
703 IF (rab2 <= rab2_max) THEN
704 CALL get_force_eam(rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam)
705 f_eam = f_eam*fac
706
707 fr(1) = -f_eam*rab(1)
708 fr(2) = -f_eam*rab(2)
709 fr(3) = -f_eam*rab(3)
710 f_nonbond(1, atom_a) = f_nonbond(1, atom_a) - fr(1)
711 f_nonbond(2, atom_a) = f_nonbond(2, atom_a) - fr(2)
712 f_nonbond(3, atom_a) = f_nonbond(3, atom_a) - fr(3)
713
714 f_nonbond(1, atom_b) = f_nonbond(1, atom_b) + fr(1)
715 f_nonbond(2, atom_b) = f_nonbond(2, atom_b) + fr(2)
716 f_nonbond(3, atom_b) = f_nonbond(3, atom_b) + fr(3)
717 IF (use_virial) THEN
718 ptens11 = ptens11 + rab(1)*fr(1)
719 ptens21 = ptens21 + rab(2)*fr(1)
720 ptens31 = ptens31 + rab(3)*fr(1)
721 ptens12 = ptens12 + rab(1)*fr(2)
722 ptens22 = ptens22 + rab(2)*fr(2)
723 ptens32 = ptens32 + rab(3)*fr(2)
724 ptens13 = ptens13 + rab(1)*fr(3)
725 ptens23 = ptens23 + rab(2)*fr(3)
726 ptens33 = ptens33 + rab(3)*fr(3)
727 END IF
728 END IF
729 END DO
730 END DO kind_group_loop1
731 END DO
732 DEALLOCATE (eam_kinds_index)
733
734 ! Special way of handling the tersoff potential..
735 IF (any_tersoff) THEN
736 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
737 CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
738 DO ilist = 1, nonbonded%nlists
739 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
740 npairs = neighbor_kind_pair%npairs
741 IF (npairs == 0) cycle
742 kind_group_loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
743 istart = neighbor_kind_pair%grp_kind_start(igrp)
744 iend = neighbor_kind_pair%grp_kind_end(igrp)
745 ikind = neighbor_kind_pair%ij_kind(1, igrp)
746 jkind = neighbor_kind_pair%ij_kind(2, igrp)
747 list => neighbor_kind_pair%list
748 cvi = neighbor_kind_pair%cell_vector
749 pot => potparm%pot(ikind, jkind)%pot
750
751 IF (pot%no_mb) cycle
752 rab2_max = pot%rcutsq
753 cell_v = matmul(cell%hmat, cvi)
754 DO i = 1, SIZE(pot%type)
755 ! TERSOFF
756 IF (pot%type(i) == tersoff_type) THEN
757 npairs = iend - istart + 1
758 tersoff => pot%set(i)%tersoff
759 ALLOCATE (sort_list(2, npairs), work_list(npairs))
760 sort_list = list(:, istart:iend)
761 ! Sort the list of neighbors, this increases the efficiency for single
762 ! potential contributions
763 CALL sort(sort_list(1, :), npairs, work_list)
764 DO ipair = 1, npairs
765 work_list(ipair) = sort_list(2, work_list(ipair))
766 END DO
767 sort_list(2, :) = work_list
768 ! find number of unique elements of array index 1
769 nunique = 1
770 DO ipair = 1, npairs - 1
771 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
772 END DO
773 ipair = 1
774 junique = sort_list(1, ipair)
775 ifirst = 1
776 DO iunique = 1, nunique
777 atom_a = junique
778 IF (glob_loc_list_a(ifirst) > atom_a) cycle
779 DO mpair = ifirst, SIZE(glob_loc_list_a)
780 IF (glob_loc_list_a(mpair) == atom_a) EXIT
781 END DO
782 ifirst = mpair
783 DO mpair = ifirst, SIZE(glob_loc_list_a)
784 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
785 END DO
786 ilast = mpair - 1
787 nloc_size = 0
788 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
789 DO WHILE (ipair <= npairs)
790 IF (sort_list(1, ipair) /= junique) EXIT
791 atom_b = sort_list(2, ipair)
792 ! Derivative terms
793 rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
794 ipair = ipair + 1
795 IF (dot_product(rtmp, rtmp) <= tersoff%rcutsq) THEN
796 CALL tersoff_forces(tersoff, r_last_update_pbc, cell_v, &
797 nloc_size, glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), &
798 atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, tersoff%rcutsq)
799 END IF
800 END DO
801 ifirst = ilast + 1
802 IF (ipair <= npairs) junique = sort_list(1, ipair)
803 END DO
804 DEALLOCATE (sort_list, work_list)
805 END IF
806 END DO
807 END DO kind_group_loop2
808 END DO
809 CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
810 END IF
811 ! Special way of handling the siepmann potential..
812 IF (any_siepmann) THEN
813 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
814 CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
815 DO ilist = 1, nonbonded%nlists
816 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
817 npairs = neighbor_kind_pair%npairs
818 IF (npairs == 0) cycle
819 kind_group_loop3: DO igrp = 1, neighbor_kind_pair%ngrp_kind
820 istart = neighbor_kind_pair%grp_kind_start(igrp)
821 iend = neighbor_kind_pair%grp_kind_end(igrp)
822 ikind = neighbor_kind_pair%ij_kind(1, igrp)
823 jkind = neighbor_kind_pair%ij_kind(2, igrp)
824 list => neighbor_kind_pair%list
825 cvi = neighbor_kind_pair%cell_vector
826 pot => potparm%pot(ikind, jkind)%pot
827
828 IF (pot%no_mb) cycle
829 rab2_max = pot%rcutsq
830 cell_v = matmul(cell%hmat, cvi)
831 DO i = 1, SIZE(pot%type)
832 ! SIEPMANN
833 IF (pot%type(i) == siepmann_type) THEN
834 npairs = iend - istart + 1
835 siepmann => pot%set(i)%siepmann
836 ALLOCATE (sort_list(2, npairs), work_list(npairs))
837 sort_list = list(:, istart:iend)
838 ! Sort the list of neighbors, this increases the efficiency for single
839 ! potential contributions
840 CALL sort(sort_list(1, :), npairs, work_list)
841 DO ipair = 1, npairs
842 work_list(ipair) = sort_list(2, work_list(ipair))
843 END DO
844 sort_list(2, :) = work_list
845 ! find number of unique elements of array index 1
846 nunique = 1
847 DO ipair = 1, npairs - 1
848 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
849 END DO
850 ipair = 1
851 junique = sort_list(1, ipair)
852 ifirst = 1
853 DO iunique = 1, nunique
854 atom_a = junique
855 IF (glob_loc_list_a(ifirst) > atom_a) cycle
856 DO mpair = ifirst, SIZE(glob_loc_list_a)
857 IF (glob_loc_list_a(mpair) == atom_a) EXIT
858 END DO
859 ifirst = mpair
860 DO mpair = ifirst, SIZE(glob_loc_list_a)
861 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
862 END DO
863 ilast = mpair - 1
864 nloc_size = 0
865 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
866 DO WHILE (ipair <= npairs)
867 IF (sort_list(1, ipair) /= junique) EXIT
868 atom_b = sort_list(2, ipair)
869 ! Derivative terms
870 rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
871 ipair = ipair + 1
872 IF (dot_product(rtmp, rtmp) <= siepmann%rcutsq) THEN
873 CALL siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, &
874 atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, &
875 particle_set)
876 CALL siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, &
877 nloc_size, glob_loc_list(:, ifirst:ilast), &
878 atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, &
879 cell, particle_set)
880 END IF
881 END DO
882 ifirst = ilast + 1
883 IF (ipair <= npairs) junique = sort_list(1, ipair)
884 END DO
885 DEALLOCATE (sort_list, work_list)
886 END IF
887 END DO
888 END DO kind_group_loop3
889 END DO
890 CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
891 END IF
892
893 ! GAL19 potential..
894 IF (any_gal) THEN
895 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
896 CALL setup_gal_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
897 DO ilist = 1, nonbonded%nlists
898 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
899 npairs = neighbor_kind_pair%npairs
900 IF (npairs == 0) cycle
901 kind_group_loop4: DO igrp = 1, neighbor_kind_pair%ngrp_kind
902 istart = neighbor_kind_pair%grp_kind_start(igrp)
903 iend = neighbor_kind_pair%grp_kind_end(igrp)
904 ikind = neighbor_kind_pair%ij_kind(1, igrp)
905 jkind = neighbor_kind_pair%ij_kind(2, igrp)
906 list => neighbor_kind_pair%list
907 cvi = neighbor_kind_pair%cell_vector
908 pot => potparm%pot(ikind, jkind)%pot
909
910 IF (pot%no_mb) cycle
911 rab2_max = pot%rcutsq
912 cell_v = matmul(cell%hmat, cvi)
913 DO i = 1, SIZE(pot%type)
914 ! GAL19
915 IF (pot%type(i) == gal_type) THEN
916 npairs = iend - istart + 1
917 gal => pot%set(i)%gal
918 ALLOCATE (sort_list(2, npairs), work_list(npairs))
919 sort_list = list(:, istart:iend)
920 ! Sort the list of neighbors, this increases the efficiency for single
921 ! potential contributions
922 CALL sort(sort_list(1, :), npairs, work_list)
923 DO ipair = 1, npairs
924 work_list(ipair) = sort_list(2, work_list(ipair))
925 END DO
926 sort_list(2, :) = work_list
927 ! find number of unique elements of array index 1
928 nunique = 1
929 DO ipair = 1, npairs - 1
930 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
931 END DO
932 ipair = 1
933 junique = sort_list(1, ipair)
934 ifirst = 1
935 DO iunique = 1, nunique
936 atom_a = junique
937 IF (glob_loc_list_a(ifirst) > atom_a) cycle
938 DO mpair = ifirst, SIZE(glob_loc_list_a)
939 IF (glob_loc_list_a(mpair) == atom_a) EXIT
940 END DO
941 ifirst = mpair
942 DO mpair = ifirst, SIZE(glob_loc_list_a)
943 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
944 END DO
945 ilast = mpair - 1
946 nloc_size = 0
947 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
948 DO WHILE (ipair <= npairs)
949 IF (sort_list(1, ipair) /= junique) EXIT
950 atom_b = sort_list(2, ipair)
951 ! Derivative terms
952 rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
953 ipair = ipair + 1
954 IF (dot_product(rtmp, rtmp) <= gal%rcutsq) THEN
955 CALL gal_forces(gal, r_last_update_pbc, &
956 atom_a, atom_b, f_nonbond, use_virial, &
957 cell, particle_set)
958 END IF
959 END DO
960 ifirst = ilast + 1
961 IF (ipair <= npairs) junique = sort_list(1, ipair)
962 END DO
963 DEALLOCATE (sort_list, work_list)
964 END IF
965 END DO
966 END DO kind_group_loop4
967 END DO
968 CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
969 END IF
970
971 ! GAL21 potential..
972 IF (any_gal21) THEN
973 NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
974 CALL setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
975 DO ilist = 1, nonbonded%nlists
976 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
977 npairs = neighbor_kind_pair%npairs
978 IF (npairs == 0) cycle
979 kind_group_loop6: DO igrp = 1, neighbor_kind_pair%ngrp_kind
980 istart = neighbor_kind_pair%grp_kind_start(igrp)
981 iend = neighbor_kind_pair%grp_kind_end(igrp)
982 ikind = neighbor_kind_pair%ij_kind(1, igrp)
983 jkind = neighbor_kind_pair%ij_kind(2, igrp)
984 list => neighbor_kind_pair%list
985 cvi = neighbor_kind_pair%cell_vector
986 pot => potparm%pot(ikind, jkind)%pot
987
988 IF (pot%no_mb) cycle
989 rab2_max = pot%rcutsq
990 cell_v = matmul(cell%hmat, cvi)
991 DO i = 1, SIZE(pot%type)
992 ! GAL21
993 IF (pot%type(i) == gal21_type) THEN
994 npairs = iend - istart + 1
995 gal21 => pot%set(i)%gal21
996 ALLOCATE (sort_list(2, npairs), work_list(npairs))
997 sort_list = list(:, istart:iend)
998 ! Sort the list of neighbors, this increases the efficiency for single
999 ! potential contributions
1000 CALL sort(sort_list(1, :), npairs, work_list)
1001 DO ipair = 1, npairs
1002 work_list(ipair) = sort_list(2, work_list(ipair))
1003 END DO
1004 sort_list(2, :) = work_list
1005 ! find number of unique elements of array index 1
1006 nunique = 1
1007 DO ipair = 1, npairs - 1
1008 IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
1009 END DO
1010 ipair = 1
1011 junique = sort_list(1, ipair)
1012 ifirst = 1
1013 DO iunique = 1, nunique
1014 atom_a = junique
1015 IF (glob_loc_list_a(ifirst) > atom_a) cycle
1016 DO mpair = ifirst, SIZE(glob_loc_list_a)
1017 IF (glob_loc_list_a(mpair) == atom_a) EXIT
1018 END DO
1019 ifirst = mpair
1020 DO mpair = ifirst, SIZE(glob_loc_list_a)
1021 IF (glob_loc_list_a(mpair) /= atom_a) EXIT
1022 END DO
1023 ilast = mpair - 1
1024 nloc_size = 0
1025 IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
1026 DO WHILE (ipair <= npairs)
1027 IF (sort_list(1, ipair) /= junique) EXIT
1028 atom_b = sort_list(2, ipair)
1029 ! Derivative terms
1030 rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
1031 ipair = ipair + 1
1032 IF (dot_product(rtmp, rtmp) <= gal21%rcutsq) THEN
1033 CALL gal21_forces(gal21, r_last_update_pbc, &
1034 atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, &
1035 cell, particle_set)
1036 END IF
1037 END DO
1038 ifirst = ilast + 1
1039 IF (ipair <= npairs) junique = sort_list(1, ipair)
1040 END DO
1041 DEALLOCATE (sort_list, work_list)
1042 END IF
1043 END DO
1044 END DO kind_group_loop6
1045 END DO
1046 CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
1047 END IF
1048
1049 ! NEQUIP
1050 IF (any_nequip) THEN
1051 CALL nequip_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
1052 END IF
1053
1054 ! ALLEGRO
1055 IF (any_allegro) THEN
1056 CALL allegro_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
1057 END IF
1058
1059 IF (use_virial) THEN
1060 pv_nonbond(1, 1) = pv_nonbond(1, 1) + ptens11
1061 pv_nonbond(1, 2) = pv_nonbond(1, 2) + ptens12
1062 pv_nonbond(1, 3) = pv_nonbond(1, 3) + ptens13
1063 pv_nonbond(2, 1) = pv_nonbond(2, 1) + ptens21
1064 pv_nonbond(2, 2) = pv_nonbond(2, 2) + ptens22
1065 pv_nonbond(2, 3) = pv_nonbond(2, 3) + ptens23
1066 pv_nonbond(3, 1) = pv_nonbond(3, 1) + ptens31
1067 pv_nonbond(3, 2) = pv_nonbond(3, 2) + ptens32
1068 pv_nonbond(3, 3) = pv_nonbond(3, 3) + ptens33
1069 END IF
1070 CALL timestop(handle)
1071 END SUBROUTINE force_nonbond_manybody
1072
1073END MODULE manybody_potential
1074
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition grid_common.h:48
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
fist_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition fist_neighbor_list_types.F:11
fist_nonbond_env_types
Definition fist_nonbond_env_types.F:13
fist_nonbond_env_types::fist_nonbond_env_get
subroutine, public fist_nonbond_env_get(fist_nonbond_env, potparm14, potparm, nonbonded, rlist_cut, rlist_lowsq, aup, lup, ei_scale14, vdw_scale14, shift_cutoff, do_electrostatics, r_last_update, r_last_update_pbc, rshell_last_update_pbc, rcore_last_update_pbc, cell_last_update, num_update, last_update, counter, natom_types, long_range_correction, ij_kind_full_fac, eam_data, quip_data, nequip_data, allegro_data, deepmd_data, charges)
sets a fist_nonbond_env
Definition fist_nonbond_env_types.F:157
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
manybody_allegro
Definition manybody_allegro.F:13
manybody_allegro::setup_allegro_arrays
subroutine, public setup_allegro_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, unique_list_a, cell)
...
Definition manybody_allegro.F:67
manybody_allegro::allegro_add_force_virial
subroutine, public allegro_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
...
Definition manybody_allegro.F:470
manybody_allegro::destroy_allegro_arrays
subroutine, public destroy_allegro_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a, unique_list_a)
...
Definition manybody_allegro.F:185
manybody_allegro::allegro_energy_store_force_virial
subroutine, public allegro_energy_store_force_virial(nonbonded, particle_set, cell, atomic_kind_set, potparm, allegro, glob_loc_list_a, r_last_update_pbc, pot_allegro, fist_nonbond_env, unique_list_a)
...
Definition manybody_allegro.F:224
manybody_deepmd
Definition manybody_deepmd.F:8
manybody_deepmd::deepmd_energy_store_force_virial
subroutine, public deepmd_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, fist_nonbond_env, pot_deepmd, para_env)
...
Definition manybody_deepmd.F:53
manybody_deepmd::deepmd_add_force_virial
subroutine, public deepmd_add_force_virial(fist_nonbond_env, force, pv_nonbond, use_virial)
...
Definition manybody_deepmd.F:194
manybody_eam
Definition manybody_eam.F:13
manybody_eam::get_force_eam
subroutine, public get_force_eam(rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam)
...
Definition manybody_eam.F:221
manybody_gal21
Implementation of the GAL21 potential.
Definition manybody_gal21.F:13
manybody_gal21::destroy_gal21_arrays
subroutine, public destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition manybody_gal21.F:844
manybody_gal21::gal21_forces
subroutine, public gal21_forces(gal21, r_last_update_pbc, iparticle, jparticle, f_nonbond, pv_nonbond, use_virial, cell, particle_set)
forces generated by the GAL2119 potential
Definition manybody_gal21.F:398
manybody_gal21::setup_gal21_arrays
subroutine, public setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition manybody_gal21.F:743
manybody_gal21::gal21_energy
subroutine, public gal21_energy(pot_loc, gal21, r_last_update_pbc, iparticle, jparticle, cell, particle_set, mm_section)
Main part of the energy evaluation of GAL2119.
Definition manybody_gal21.F:61
manybody_gal
Implementation of the GAL19 potential.
Definition manybody_gal.F:13
manybody_gal::destroy_gal_arrays
subroutine, public destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition manybody_gal.F:744
manybody_gal::setup_gal_arrays
subroutine, public setup_gal_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition manybody_gal.F:643
manybody_gal::gal_energy
subroutine, public gal_energy(pot_loc, gal, r_last_update_pbc, iparticle, jparticle, cell, particle_set, mm_section)
Main part of the energy evaluation of GAL19.
Definition manybody_gal.F:61
manybody_gal::gal_forces
subroutine, public gal_forces(gal, r_last_update_pbc, iparticle, jparticle, f_nonbond, use_virial, cell, particle_set)
forces generated by the GAL19 potential
Definition manybody_gal.F:373
manybody_nequip
Definition manybody_nequip.F:13
manybody_nequip::nequip_add_force_virial
subroutine, public nequip_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
...
Definition manybody_nequip.F:474
manybody_nequip::nequip_energy_store_force_virial
subroutine, public nequip_energy_store_force_virial(nonbonded, particle_set, cell, atomic_kind_set, potparm, nequip, glob_loc_list_a, r_last_update_pbc, pot_nequip, fist_nonbond_env, para_env)
...
Definition manybody_nequip.F:206
manybody_nequip::destroy_nequip_arrays
subroutine, public destroy_nequip_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition manybody_nequip.F:170
manybody_nequip::setup_nequip_arrays
subroutine, public setup_nequip_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition manybody_nequip.F:66
manybody_potential
Definition manybody_potential.F:16
manybody_potential::energy_manybody
subroutine, public energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, particle_set, cell, pot_manybody, para_env, mm_section)
computes the embedding contribution to the energy
Definition manybody_potential.F:95
manybody_potential::force_nonbond_manybody
subroutine, public force_nonbond_manybody(fist_nonbond_env, particle_set, cell, f_nonbond, pv_nonbond, use_virial)
...
Definition manybody_potential.F:575
manybody_quip
Definition manybody_quip.F:8
manybody_quip::quip_add_force_virial
subroutine, public quip_add_force_virial(fist_nonbond_env, force, virial)
...
Definition manybody_quip.F:195
manybody_quip::quip_energy_store_force_virial
subroutine, public quip_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, fist_nonbond_env, pot_quip, para_env)
...
Definition manybody_quip.F:56
manybody_siepmann
implementation of dipole and three-body part of Siepmann-Sprik potential dipole term: 3rd term in Eq....
Definition manybody_siepmann.F:17
manybody_siepmann::siepmann_energy
subroutine, public siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size, full_loc_list, cell_v, cell, drij, particle_set, nr_oh, nr_h3o, nr_o)
energy of two-body dipole term and three-body term
Definition manybody_siepmann.F:71
manybody_siepmann::setup_siepmann_arrays
subroutine, public setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition manybody_siepmann.F:707
manybody_siepmann::siepmann_forces_v2
subroutine, public siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, iparticle, jparticle, f_nonbond, use_virial, rcutsq, particle_set)
forces generated by the dipole term
Definition manybody_siepmann.F:628
manybody_siepmann::siepmann_forces_v3
subroutine, public siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, n_loc_size, full_loc_list, iparticle, jparticle, f_nonbond, use_virial, rcutsq, cell, particle_set)
forces generated by the three-body term
Definition manybody_siepmann.F:539
manybody_siepmann::print_nr_ions_siepmann
subroutine, public print_nr_ions_siepmann(nr_ions, mm_section, para_env, print_oh, print_h3o, print_o)
prints the number of OH- ions or H3O+ ions near surface
Definition manybody_siepmann.F:835
manybody_siepmann::destroy_siepmann_arrays
subroutine, public destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition manybody_siepmann.F:808
manybody_tersoff
Definition manybody_tersoff.F:13
manybody_tersoff::setup_tersoff_arrays
subroutine, public setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
...
Definition manybody_tersoff.F:696
manybody_tersoff::tersoff_forces
subroutine, public tersoff_forces(tersoff, r_last_update_pbc, cell_v, n_loc_size, full_loc_list, loc_cell_v, iparticle, jparticle, f_nonbond, pv_nonbond, use_virial, rcutsq)
...
Definition manybody_tersoff.F:558
manybody_tersoff::destroy_tersoff_arrays
subroutine, public destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
...
Definition manybody_tersoff.F:800
manybody_tersoff::tersoff_energy
subroutine, public tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, full_loc_list, loc_cell_v, cell_v, drij)
...
Definition manybody_tersoff.F:53
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
pair_potential_types
Definition pair_potential_types.F:14
pair_potential_types::allegro_type
integer, parameter, public allegro_type
Definition pair_potential_types.F:35
pair_potential_types::gal_type
integer, parameter, public gal_type
Definition pair_potential_types.F:35
pair_potential_types::nequip_type
integer, parameter, public nequip_type
Definition pair_potential_types.F:35
pair_potential_types::deepmd_type
integer, parameter, public deepmd_type
Definition pair_potential_types.F:35
pair_potential_types::quip_type
integer, parameter, public quip_type
Definition pair_potential_types.F:35
pair_potential_types::siepmann_type
integer, parameter, public siepmann_type
Definition pair_potential_types.F:35
pair_potential_types::gal21_type
integer, parameter, public gal21_type
Definition pair_potential_types.F:35
pair_potential_types::ea_type
integer, parameter, public ea_type
Definition pair_potential_types.F:35
pair_potential_types::tersoff_type
integer, parameter, public tersoff_type
Definition pair_potential_types.F:35
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
util
All kind of helpful little routines.
Definition util.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
fist_neighbor_list_types::fist_neighbor_type
Definition fist_neighbor_list_types.F:48
fist_neighbor_list_types::neighbor_kind_pairs_type
Definition fist_neighbor_list_types.F:26
fist_nonbond_env_types::eam_type
Definition fist_nonbond_env_types.F:48
fist_nonbond_env_types::fist_nonbond_env_type
Definition fist_nonbond_env_types.F:81
fist_nonbond_env_types::pos_type
Definition fist_nonbond_env_types.F:44
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
pair_potential_types::allegro_pot_type
Definition pair_potential_types.F:203
pair_potential_types::eam_pot_type
Definition pair_potential_types.F:170
pair_potential_types::gal21_pot_type
Definition pair_potential_types.F:305
pair_potential_types::gal_pot_type
Definition pair_potential_types.F:282
pair_potential_types::nequip_pot_type
Definition pair_potential_types.F:192
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:412
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:379
pair_potential_types::siepmann_pot_type
Definition pair_potential_types.F:269
pair_potential_types::tersoff_pot_type
Definition pair_potential_types.F:250
particle_types::particle_type
Definition particle_types.F:35