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 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 : ! **************************************************************************************************
16 : MODULE 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_ace, ONLY: ace_add_force_virial,&
30 : ace_energy_store_force_virial
31 : USE manybody_deepmd, ONLY: deepmd_add_force_virial,&
32 : deepmd_energy_store_force_virial
33 : USE manybody_e3nn, ONLY: e3nn_add_force_virial,&
34 : e3nn_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_siepmann, ONLY: destroy_siepmann_arrays,&
45 : print_nr_ions_siepmann,&
46 : setup_siepmann_arrays,&
47 : siepmann_energy,&
48 : siepmann_forces_v2,&
49 : siepmann_forces_v3
50 : USE manybody_tersoff, ONLY: destroy_tersoff_arrays,&
51 : setup_tersoff_arrays,&
52 : tersoff_energy,&
53 : tersoff_forces
54 : USE message_passing, ONLY: mp_para_env_type
55 : USE pair_potential_types, ONLY: &
56 : ace_type, allegro_type, deepmd_type, ea_type, eam_pot_type, gal21_pot_type, gal21_type, &
57 : gal_pot_type, gal_type, mace_type, nequip_type, pair_potential_pp_type, &
58 : pair_potential_single_type, siepmann_pot_type, siepmann_type, tersoff_pot_type, &
59 : tersoff_type
60 : USE particle_types, ONLY: particle_type
61 : USE util, ONLY: sort
62 : #include "./base/base_uses.f90"
63 :
64 : IMPLICIT NONE
65 :
66 : PRIVATE
67 : PUBLIC :: energy_manybody
68 : PUBLIC :: force_nonbond_manybody
69 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'manybody_potential'
70 :
71 : CONTAINS
72 :
73 : ! **************************************************************************************************
74 : !> \brief computes the embedding contribution to the energy
75 : !> \param fist_nonbond_env ...
76 : !> \param atomic_kind_set ...
77 : !> \param local_particles ...
78 : !> \param particle_set ...
79 : !> \param cell ...
80 : !> \param pot_manybody ...
81 : !> \param para_env ...
82 : !> \param mm_section ...
83 : !> \param use_virial ...
84 : !> \par History
85 : !> tlaino [2007] - New algorithm for tersoff potential
86 : !> \author CJM, I-Feng W. Kuo, Teodoro Laino
87 : ! **************************************************************************************************
88 77140 : SUBROUTINE energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, &
89 : particle_set, cell, pot_manybody, para_env, mm_section, use_virial)
90 :
91 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
92 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
93 : TYPE(distribution_1d_type), POINTER :: local_particles
94 : TYPE(particle_type), POINTER :: particle_set(:)
95 : TYPE(cell_type), POINTER :: cell
96 : REAL(dp), INTENT(INOUT) :: pot_manybody
97 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
98 : TYPE(section_vals_type), POINTER :: mm_section
99 : LOGICAL, INTENT(IN) :: use_virial
100 :
101 : CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_manybody'
102 :
103 : INTEGER :: atom_a, atom_b, handle, i, iend, ifirst, igrp, ikind, ilast, ilist, indexa, &
104 : ipair, iparticle, iparticle_local, istart, iunique, jkind, junique, mpair, nkinds, &
105 : nloc_size, npairs, nparticle, nparticle_local, nr_h3O, nr_o, nr_oh, nunique
106 77140 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list
107 77140 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list
108 : LOGICAL :: any_ace, any_allegro, any_deepmd, &
109 : any_gal, any_gal21, any_mace, &
110 : any_nequip, any_siepmann, any_tersoff
111 : REAL(KIND=dp) :: drij, embed, pot_ace, pot_allegro, &
112 : pot_deepmd, pot_loc, pot_mace, &
113 : pot_nequip, qr, rab2_max, rij(3)
114 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi
115 77140 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
116 77140 : REAL(KIND=dp), POINTER :: fembed(:)
117 : TYPE(eam_pot_type), POINTER :: eam
118 77140 : TYPE(eam_type), DIMENSION(:), POINTER :: eam_data
119 : TYPE(fist_neighbor_type), POINTER :: nonbonded
120 : TYPE(gal21_pot_type), POINTER :: gal21
121 : TYPE(gal_pot_type), POINTER :: gal
122 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
123 : TYPE(pair_potential_pp_type), POINTER :: potparm
124 : TYPE(pair_potential_single_type), POINTER :: pot
125 77140 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
126 : TYPE(siepmann_pot_type), POINTER :: siepmann
127 : TYPE(tersoff_pot_type), POINTER :: tersoff
128 :
129 77140 : NULLIFY (eam, siepmann, tersoff, gal, gal21)
130 77140 : any_tersoff = .FALSE.
131 77140 : any_siepmann = .FALSE.
132 77140 : any_gal = .FALSE.
133 77140 : any_gal21 = .FALSE.
134 77140 : any_allegro = .FALSE.
135 77140 : any_nequip = .FALSE.
136 77140 : any_mace = .FALSE.
137 77140 : any_ace = .FALSE.
138 77140 : any_deepmd = .FALSE.
139 77140 : CALL timeset(routineN, handle)
140 : CALL fist_nonbond_env_get(fist_nonbond_env, r_last_update_pbc=r_last_update_pbc, &
141 77140 : potparm=potparm, eam_data=eam_data)
142 : ! EAM requires a single loop
143 304504 : DO ikind = 1, SIZE(atomic_kind_set)
144 227364 : pot => potparm%pot(ikind, ikind)%pot
145 531916 : DO i = 1, SIZE(pot%type)
146 227412 : IF (pot%type(i) /= ea_type) CYCLE
147 488 : eam => pot%set(i)%eam
148 488 : nparticle = SIZE(particle_set)
149 1464 : ALLOCATE (fembed(nparticle))
150 14258 : fembed(:) = 0._dp
151 488 : CPASSERT(ASSOCIATED(eam_data))
152 : ! computation of embedding function and energy
153 488 : nparticle_local = local_particles%n_el(ikind)
154 4136 : DO iparticle_local = 1, nparticle_local
155 3648 : iparticle = local_particles%list(ikind)%array(iparticle_local)
156 3648 : indexa = INT(eam_data(iparticle)%rho/eam%drhoar) + 1
157 3648 : IF (indexa > eam%npoints - 1) indexa = eam%npoints - 1
158 3648 : qr = eam_data(iparticle)%rho - eam%rhoval(indexa)
159 :
160 3648 : embed = eam%frho(indexa) + qr*eam%frhop(indexa)
161 3648 : fembed(iparticle) = eam%frhop(indexa) + qr*(eam%frhop(indexa + 1) - eam%frhop(indexa))/eam%drhoar
162 :
163 4136 : pot_manybody = pot_manybody + embed
164 : END DO
165 : ! communicate data
166 28028 : CALL para_env%sum(fembed)
167 14258 : DO iparticle = 1, nparticle
168 14258 : IF (particle_set(iparticle)%atomic_kind%kind_number == ikind) THEN
169 7296 : eam_data(iparticle)%f_embed = fembed(iparticle)
170 : END IF
171 : END DO
172 :
173 454776 : DEALLOCATE (fembed)
174 : END DO
175 : END DO
176 : ! Other manybody potential
177 304504 : DO ikind = 1, SIZE(atomic_kind_set)
178 1773650 : DO jkind = ikind, SIZE(atomic_kind_set)
179 1469146 : pot => potparm%pot(ikind, jkind)%pot
180 2935368 : any_tersoff = any_tersoff .OR. ANY(pot%type == tersoff_type)
181 2938334 : any_nequip = any_nequip .OR. ANY(pot%type == nequip_type)
182 2938338 : any_mace = any_mace .OR. ANY(pot%type == mace_type)
183 2937722 : any_ace = any_ace .OR. ANY(pot%type == ace_type)
184 2938334 : any_allegro = any_allegro .OR. ANY(pot%type == allegro_type)
185 2938334 : any_deepmd = any_deepmd .OR. ANY(pot%type == deepmd_type)
186 2938319 : any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type)
187 2938339 : any_gal = any_gal .OR. ANY(pot%type == gal_type)
188 3165703 : any_gal21 = any_gal21 .OR. ANY(pot%type == gal21_type)
189 : END DO
190 : END DO
191 77140 : CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, natom_types=nkinds)
192 : ! NEQUIP
193 77140 : IF (any_nequip) THEN
194 2 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
195 : CALL e3nn_energy_store_force_virial(nonbonded, particle_set, local_particles, cell, &
196 : atomic_kind_set, potparm, r_last_update_pbc, &
197 : pot_nequip, fist_nonbond_env, &
198 2 : para_env, use_virial, nequip_type)
199 2 : pot_manybody = pot_manybody + pot_nequip
200 : END IF
201 : ! ALLEGRO
202 77140 : IF (any_allegro) THEN
203 2 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
204 : CALL e3nn_energy_store_force_virial(nonbonded, particle_set, local_particles, cell, &
205 : atomic_kind_set, potparm, r_last_update_pbc, &
206 : pot_allegro, fist_nonbond_env, &
207 2 : para_env, use_virial, allegro_type)
208 2 : pot_manybody = pot_manybody + pot_allegro
209 : END IF
210 : ! MACE (reuses the NequIP message-passing evaluation path)
211 77140 : IF (any_mace) THEN
212 2 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
213 : CALL e3nn_energy_store_force_virial(nonbonded, particle_set, local_particles, cell, &
214 : atomic_kind_set, potparm, r_last_update_pbc, &
215 : pot_mace, fist_nonbond_env, &
216 2 : para_env, use_virial, mace_type)
217 2 : pot_manybody = pot_manybody + pot_mace
218 : END IF
219 : ! ACE
220 77140 : IF (any_ace) THEN
221 : CALL ace_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, &
222 206 : fist_nonbond_env, pot_ace)
223 206 : pot_manybody = pot_manybody + pot_ace
224 : END IF
225 : ! DEEPMD
226 77140 : IF (any_deepmd) THEN
227 : CALL deepmd_energy_store_force_virial(particle_set, cell, atomic_kind_set, potparm, &
228 2 : fist_nonbond_env, pot_deepmd, para_env)
229 2 : pot_manybody = pot_manybody + pot_deepmd
230 : END IF
231 :
232 : ! TERSOFF
233 77140 : IF (any_tersoff) THEN
234 2664 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
235 2664 : CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
236 101052 : DO ilist = 1, nonbonded%nlists
237 98388 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
238 98388 : npairs = neighbor_kind_pair%npairs
239 98388 : IF (npairs == 0) CYCLE
240 68185 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
241 33321 : istart = neighbor_kind_pair%grp_kind_start(igrp)
242 33321 : iend = neighbor_kind_pair%grp_kind_end(igrp)
243 33321 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
244 33321 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
245 33321 : list => neighbor_kind_pair%list
246 133284 : cvi = neighbor_kind_pair%cell_vector
247 33321 : pot => potparm%pot(ikind, jkind)%pot
248 165032 : DO i = 1, SIZE(pot%type)
249 33323 : IF (pot%type(i) /= tersoff_type) CYCLE
250 33280 : rab2_max = pot%set(i)%tersoff%rcutsq
251 432640 : cell_v = MATMUL(cell%hmat, cvi)
252 33280 : pot => potparm%pot(ikind, jkind)%pot
253 33280 : tersoff => pot%set(i)%tersoff
254 33280 : npairs = iend - istart + 1
255 66601 : IF (npairs /= 0) THEN
256 166400 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
257 29009284 : sort_list = list(:, istart:iend)
258 : ! Sort the list of neighbors, this increases the efficiency for single
259 : ! potential contributions
260 33280 : CALL sort(sort_list(1, :), npairs, work_list)
261 4862614 : DO ipair = 1, npairs
262 4862614 : work_list(ipair) = sort_list(2, work_list(ipair))
263 : END DO
264 9691948 : sort_list(2, :) = work_list
265 : ! find number of unique elements of array index 1
266 33280 : nunique = 1
267 4829334 : DO ipair = 1, npairs - 1
268 4829334 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
269 : END DO
270 33280 : ipair = 1
271 33280 : junique = sort_list(1, ipair)
272 33280 : ifirst = 1
273 344350 : DO iunique = 1, nunique
274 311070 : atom_a = junique
275 311070 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
276 61587820 : DO mpair = ifirst, SIZE(glob_loc_list_a)
277 61587820 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
278 : END DO
279 70828815 : ifirst = mpair
280 70828815 : DO mpair = ifirst, SIZE(glob_loc_list_a)
281 70828815 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
282 : END DO
283 311070 : ilast = mpair - 1
284 311070 : nloc_size = 0
285 311070 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
286 5140404 : DO WHILE (ipair <= npairs)
287 5107124 : IF (sort_list(1, ipair) /= junique) EXIT
288 4829334 : atom_b = sort_list(2, ipair)
289 : ! Energy terms
290 4829334 : pot_loc = 0.0_dp
291 19317336 : rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
292 19317336 : drij = DOT_PRODUCT(rij, rij)
293 4829334 : ipair = ipair + 1
294 4829334 : IF (drij > rab2_max) CYCLE
295 250380 : drij = SQRT(drij)
296 : CALL tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
297 250380 : glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), cell_v, drij)
298 5107124 : pot_manybody = pot_manybody + 0.5_dp*pot_loc
299 : END DO
300 311070 : ifirst = ilast + 1
301 344350 : IF (ipair <= npairs) junique = sort_list(1, ipair)
302 : END DO
303 33280 : DEALLOCATE (sort_list, work_list)
304 : END IF
305 : END DO
306 : END DO Kind_Group_Loop
307 : END DO
308 2664 : CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
309 : END IF
310 :
311 : !SIEPMANN POTENTIAL
312 77140 : IF (any_siepmann) THEN
313 21 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
314 21 : nr_oh = 0
315 21 : nr_h3O = 0
316 21 : nr_o = 0
317 21 : CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
318 588 : DO ilist = 1, nonbonded%nlists
319 567 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
320 567 : npairs = neighbor_kind_pair%npairs
321 567 : IF (npairs == 0) CYCLE
322 918 : Kind_Group_Loop_2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
323 708 : istart = neighbor_kind_pair%grp_kind_start(igrp)
324 708 : iend = neighbor_kind_pair%grp_kind_end(igrp)
325 708 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
326 708 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
327 708 : list => neighbor_kind_pair%list
328 2832 : cvi = neighbor_kind_pair%cell_vector
329 708 : pot => potparm%pot(ikind, jkind)%pot
330 1983 : DO i = 1, SIZE(pot%type)
331 708 : IF (pot%type(i) /= siepmann_type) CYCLE
332 165 : rab2_max = pot%set(i)%siepmann%rcutsq
333 2145 : cell_v = MATMUL(cell%hmat, cvi)
334 165 : pot => potparm%pot(ikind, jkind)%pot
335 165 : siepmann => pot%set(i)%siepmann
336 165 : npairs = iend - istart + 1
337 873 : IF (npairs /= 0) THEN
338 825 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
339 109533 : sort_list = list(:, istart:iend)
340 : ! Sort the list of neighbors, this increases the efficiency for single
341 : ! potential contributions
342 165 : CALL sort(sort_list(1, :), npairs, work_list)
343 18393 : DO ipair = 1, npairs
344 18393 : work_list(ipair) = sort_list(2, work_list(ipair))
345 : END DO
346 36621 : sort_list(2, :) = work_list
347 : ! find number of unique elements of array index 1
348 165 : nunique = 1
349 18228 : DO ipair = 1, npairs - 1
350 18228 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
351 : END DO
352 165 : ipair = 1
353 165 : junique = sort_list(1, ipair)
354 165 : ifirst = 1
355 5340 : DO iunique = 1, nunique
356 5175 : atom_a = junique
357 5175 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
358 91602 : DO mpair = ifirst, SIZE(glob_loc_list_a)
359 91602 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
360 : END DO
361 62187 : ifirst = mpair
362 62187 : DO mpair = ifirst, SIZE(glob_loc_list_a)
363 62187 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
364 : END DO
365 5175 : ilast = mpair - 1
366 5175 : nloc_size = 0
367 5175 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
368 23403 : DO WHILE (ipair <= npairs)
369 23238 : IF (sort_list(1, ipair) /= junique) EXIT
370 18228 : atom_b = sort_list(2, ipair)
371 : ! Energy terms
372 18228 : pot_loc = 0.0_dp
373 72912 : rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
374 72912 : drij = DOT_PRODUCT(rij, rij)
375 18228 : ipair = ipair + 1
376 18228 : IF (drij > rab2_max) CYCLE
377 318 : drij = SQRT(drij)
378 : CALL siepmann_energy(pot_loc, siepmann, r_last_update_pbc, atom_a, atom_b, nloc_size, &
379 : glob_loc_list(:, ifirst:ilast), cell_v, cell, drij, &
380 318 : particle_set, nr_oh, nr_h3O, nr_o)
381 23238 : pot_manybody = pot_manybody + pot_loc
382 : END DO
383 5175 : ifirst = ilast + 1
384 5340 : IF (ipair <= npairs) junique = sort_list(1, ipair)
385 : END DO
386 165 : DEALLOCATE (sort_list, work_list)
387 : END IF
388 : END DO
389 : END DO Kind_Group_Loop_2
390 : END DO
391 21 : CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
392 : CALL print_nr_ions_siepmann(nr_oh, mm_section, para_env, print_oh=.TRUE., &
393 21 : print_h3o=.FALSE., print_o=.FALSE.)
394 : CALL print_nr_ions_siepmann(nr_h3o, mm_section, para_env, print_oh=.FALSE., &
395 21 : print_h3o=.TRUE., print_o=.FALSE.)
396 : CALL print_nr_ions_siepmann(nr_o, mm_section, para_env, print_oh=.FALSE., &
397 21 : print_h3o=.FALSE., print_o=.TRUE.)
398 : END IF
399 :
400 : !GAL19 POTENTIAL
401 77140 : IF (any_gal) THEN
402 1 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
403 1 : CALL setup_gal_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
404 28 : DO ilist = 1, nonbonded%nlists
405 27 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
406 27 : npairs = neighbor_kind_pair%npairs
407 27 : IF (npairs == 0) CYCLE
408 168 : Kind_Group_Loop_3: DO igrp = 1, neighbor_kind_pair%ngrp_kind
409 158 : istart = neighbor_kind_pair%grp_kind_start(igrp)
410 158 : iend = neighbor_kind_pair%grp_kind_end(igrp)
411 158 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
412 158 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
413 158 : list => neighbor_kind_pair%list
414 632 : cvi = neighbor_kind_pair%cell_vector
415 158 : pot => potparm%pot(ikind, jkind)%pot
416 343 : DO i = 1, SIZE(pot%type)
417 158 : IF (pot%type(i) /= gal_type) CYCLE
418 9 : rab2_max = pot%set(i)%gal%rcutsq
419 117 : cell_v = MATMUL(cell%hmat, cvi)
420 9 : pot => potparm%pot(ikind, jkind)%pot
421 9 : gal => pot%set(i)%gal
422 9 : npairs = iend - istart + 1
423 167 : IF (npairs /= 0) THEN
424 45 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
425 45609 : sort_list = list(:, istart:iend)
426 : ! Sort the list of neighbors, this increases the efficiency for single
427 : ! potential contributions
428 9 : CALL sort(sort_list(1, :), npairs, work_list)
429 7609 : DO ipair = 1, npairs
430 7609 : work_list(ipair) = sort_list(2, work_list(ipair))
431 : END DO
432 15209 : sort_list(2, :) = work_list
433 : ! find number of unique elements of array index 1
434 9 : nunique = 1
435 7600 : DO ipair = 1, npairs - 1
436 7600 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
437 : END DO
438 9 : ipair = 1
439 9 : junique = sort_list(1, ipair)
440 9 : ifirst = 1
441 659 : DO iunique = 1, nunique
442 650 : atom_a = junique
443 650 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
444 36198 : DO mpair = ifirst, SIZE(glob_loc_list_a)
445 36198 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
446 : END DO
447 24581 : ifirst = mpair
448 24581 : DO mpair = ifirst, SIZE(glob_loc_list_a)
449 24581 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
450 : END DO
451 650 : ilast = mpair - 1
452 650 : nloc_size = 0
453 650 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
454 8250 : DO WHILE (ipair <= npairs)
455 8241 : IF (sort_list(1, ipair) /= junique) EXIT
456 7600 : atom_b = sort_list(2, ipair)
457 : ! Energy terms
458 7600 : pot_loc = 0.0_dp
459 30400 : rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
460 30400 : drij = DOT_PRODUCT(rij, rij)
461 7600 : ipair = ipair + 1
462 7600 : IF (drij > rab2_max) CYCLE
463 2004 : drij = SQRT(drij)
464 : CALL gal_energy(pot_loc, gal, r_last_update_pbc, atom_a, atom_b, &
465 2004 : cell, particle_set, mm_section)
466 :
467 8241 : pot_manybody = pot_manybody + pot_loc
468 : END DO
469 650 : ifirst = ilast + 1
470 659 : IF (ipair <= npairs) junique = sort_list(1, ipair)
471 : END DO
472 9 : DEALLOCATE (sort_list, work_list)
473 : END IF
474 : END DO
475 : END DO Kind_Group_Loop_3
476 : END DO
477 1 : CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
478 : END IF
479 :
480 : !GAL21 POTENTIAL
481 77140 : IF (any_gal21) THEN
482 1 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
483 1 : CALL setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
484 28 : DO ilist = 1, nonbonded%nlists
485 27 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
486 27 : npairs = neighbor_kind_pair%npairs
487 27 : IF (npairs == 0) CYCLE
488 168 : Kind_Group_Loop_5: DO igrp = 1, neighbor_kind_pair%ngrp_kind
489 158 : istart = neighbor_kind_pair%grp_kind_start(igrp)
490 158 : iend = neighbor_kind_pair%grp_kind_end(igrp)
491 158 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
492 158 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
493 158 : list => neighbor_kind_pair%list
494 632 : cvi = neighbor_kind_pair%cell_vector
495 158 : pot => potparm%pot(ikind, jkind)%pot
496 343 : DO i = 1, SIZE(pot%type)
497 158 : IF (pot%type(i) /= gal21_type) CYCLE
498 9 : rab2_max = pot%set(i)%gal21%rcutsq
499 117 : cell_v = MATMUL(cell%hmat, cvi)
500 9 : pot => potparm%pot(ikind, jkind)%pot
501 9 : gal21 => pot%set(i)%gal21
502 9 : npairs = iend - istart + 1
503 167 : IF (npairs /= 0) THEN
504 45 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
505 52809 : sort_list = list(:, istart:iend)
506 : ! Sort the list of neighbors, this increases the efficiency for single
507 : ! potential contributions
508 9 : CALL sort(sort_list(1, :), npairs, work_list)
509 8809 : DO ipair = 1, npairs
510 8809 : work_list(ipair) = sort_list(2, work_list(ipair))
511 : END DO
512 17609 : sort_list(2, :) = work_list
513 : ! find number of unique elements of array index 1
514 9 : nunique = 1
515 8800 : DO ipair = 1, npairs - 1
516 8800 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
517 : END DO
518 9 : ipair = 1
519 9 : junique = sort_list(1, ipair)
520 9 : ifirst = 1
521 710 : DO iunique = 1, nunique
522 701 : atom_a = junique
523 701 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
524 42242 : DO mpair = ifirst, SIZE(glob_loc_list_a)
525 42242 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
526 : END DO
527 30069 : ifirst = mpair
528 30069 : DO mpair = ifirst, SIZE(glob_loc_list_a)
529 30069 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
530 : END DO
531 701 : ilast = mpair - 1
532 701 : nloc_size = 0
533 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
534 9501 : DO WHILE (ipair <= npairs)
535 9492 : IF (sort_list(1, ipair) /= junique) EXIT
536 8800 : atom_b = sort_list(2, ipair)
537 : ! Energy terms
538 8800 : pot_loc = 0.0_dp
539 35200 : rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
540 35200 : drij = DOT_PRODUCT(rij, rij)
541 8800 : ipair = ipair + 1
542 8800 : IF (drij > rab2_max) CYCLE
543 : drij = SQRT(drij)
544 : CALL gal21_energy(pot_loc, gal21, r_last_update_pbc, atom_a, atom_b, &
545 5732 : cell, particle_set, mm_section)
546 :
547 9492 : pot_manybody = pot_manybody + pot_loc
548 : END DO
549 701 : ifirst = ilast + 1
550 710 : IF (ipair <= npairs) junique = sort_list(1, ipair)
551 : END DO
552 9 : DEALLOCATE (sort_list, work_list)
553 : END IF
554 : END DO
555 : END DO Kind_Group_Loop_5
556 : END DO
557 1 : CALL destroy_gal21_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
558 : END IF
559 :
560 77140 : CALL timestop(handle)
561 77140 : END SUBROUTINE energy_manybody
562 :
563 : ! **************************************************************************************************
564 : !> \brief ...
565 : !> \param fist_nonbond_env ...
566 : !> \param particle_set ...
567 : !> \param cell ...
568 : !> \param f_nonbond ...
569 : !> \param pv_nonbond ...
570 : !> \param use_virial ...
571 : !> \par History
572 : !> Fast implementation of the tersoff potential - [tlaino] 2007
573 : !> \author I-Feng W. Kuo, Teodoro Laino
574 : ! **************************************************************************************************
575 66832 : SUBROUTINE force_nonbond_manybody(fist_nonbond_env, particle_set, cell, &
576 66832 : f_nonbond, pv_nonbond, use_virial)
577 :
578 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
579 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
580 : TYPE(cell_type), POINTER :: cell
581 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
582 : LOGICAL, INTENT(IN) :: use_virial
583 :
584 : CHARACTER(LEN=*), PARAMETER :: routineN = 'force_nonbond_manybody'
585 :
586 : INTEGER :: atom_a, atom_b, handle, i, i_a, i_b, iend, ifirst, igrp, ikind, ilast, ilist, &
587 : ipair, istart, iunique, jkind, junique, kind_a, kind_b, mpair, nkinds, nloc_size, npairs, &
588 : nunique
589 66832 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index
590 66832 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list
591 66832 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list
592 : LOGICAL :: any_ace, any_allegro, any_deepmd, &
593 : any_gal, any_gal21, any_mace, &
594 : any_nequip, any_siepmann, any_tersoff
595 : REAL(KIND=dp) :: f_eam, fac, fr(3), ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, &
596 : ptens31, ptens32, ptens33, rab(3), rab2, rab2_max, rtmp(3)
597 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi
598 66832 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
599 : TYPE(eam_pot_type), POINTER :: eam_a, eam_b
600 66832 : TYPE(eam_type), DIMENSION(:), POINTER :: eam_data
601 : TYPE(fist_neighbor_type), POINTER :: nonbonded
602 : TYPE(gal21_pot_type), POINTER :: gal21
603 : TYPE(gal_pot_type), POINTER :: gal
604 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
605 : TYPE(pair_potential_pp_type), POINTER :: potparm
606 : TYPE(pair_potential_single_type), POINTER :: pot
607 66832 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
608 : TYPE(siepmann_pot_type), POINTER :: siepmann
609 : TYPE(tersoff_pot_type), POINTER :: tersoff
610 :
611 66832 : any_tersoff = .FALSE.
612 66832 : any_allegro = .FALSE.
613 66832 : any_nequip = .FALSE.
614 66832 : any_mace = .FALSE.
615 66832 : any_siepmann = .FALSE.
616 66832 : any_ace = .FALSE.
617 66832 : any_deepmd = .FALSE.
618 66832 : any_gal = .FALSE.
619 66832 : any_gal21 = .FALSE.
620 66832 : CALL timeset(routineN, handle)
621 66832 : NULLIFY (eam_a, eam_b, tersoff, siepmann, gal, gal21)
622 :
623 : CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, potparm=potparm, &
624 66832 : natom_types=nkinds, eam_data=eam_data, r_last_update_pbc=r_last_update_pbc)
625 :
626 : ! Initializing the potential energy, pressure tensor and force
627 : IF (use_virial) THEN
628 : ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
629 : ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
630 : ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
631 : END IF
632 :
633 66832 : nkinds = SIZE(potparm%pot, 1)
634 267328 : ALLOCATE (eam_kinds_index(nkinds, nkinds))
635 2943336 : eam_kinds_index = -1
636 273598 : DO ikind = 1, nkinds
637 1711850 : DO jkind = ikind, nkinds
638 3083318 : DO i = 1, SIZE(potparm%pot(ikind, jkind)%pot%type)
639 2876552 : IF (potparm%pot(ikind, jkind)%pot%type(i) == ea_type) THEN
640 : ! At the moment we allow only 1 EAM per each kinds pair..
641 692 : CPASSERT(eam_kinds_index(ikind, jkind) == -1)
642 692 : CPASSERT(eam_kinds_index(jkind, ikind) == -1)
643 692 : eam_kinds_index(ikind, jkind) = i
644 692 : eam_kinds_index(jkind, ikind) = i
645 : END IF
646 : END DO
647 : END DO
648 : END DO
649 273598 : DO ikind = 1, nkinds
650 1711850 : DO jkind = ikind, nkinds
651 3082700 : any_ace = any_ace .OR. ANY(potparm%pot(ikind, jkind)%pot%type == ace_type)
652 : END DO
653 : END DO
654 : ! ACE
655 66832 : IF (any_ace) THEN
656 206 : CALL ace_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
657 : END IF
658 :
659 273598 : DO ikind = 1, nkinds
660 1711850 : DO jkind = ikind, nkinds
661 3083312 : any_deepmd = any_deepmd .OR. ANY(potparm%pot(ikind, jkind)%pot%type == deepmd_type)
662 : END DO
663 : END DO
664 : ! DEEPMD
665 66832 : IF (any_deepmd) THEN
666 2 : CALL deepmd_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
667 : END IF
668 :
669 : ! NEQUIP
670 273598 : DO ikind = 1, nkinds
671 1711850 : DO jkind = ikind, nkinds
672 3083312 : any_nequip = any_nequip .OR. ANY(potparm%pot(ikind, jkind)%pot%type == nequip_type)
673 : END DO
674 : END DO
675 66832 : IF (any_nequip) THEN
676 2 : CALL e3nn_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
677 : END IF
678 :
679 : ! ALLEGRO
680 273598 : DO ikind = 1, nkinds
681 1711850 : DO jkind = ikind, nkinds
682 3083312 : any_allegro = any_allegro .OR. ANY(potparm%pot(ikind, jkind)%pot%type == allegro_type)
683 : END DO
684 : END DO
685 66832 : IF (any_allegro) THEN
686 2 : CALL e3nn_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
687 : END IF
688 :
689 : ! MACE (reuses the NequIP force/virial accumulation)
690 273598 : DO ikind = 1, nkinds
691 1711850 : DO jkind = ikind, nkinds
692 3083316 : any_mace = any_mace .OR. ANY(potparm%pot(ikind, jkind)%pot%type == mace_type)
693 : END DO
694 : END DO
695 66832 : IF (any_mace) THEN
696 2 : CALL e3nn_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
697 : END IF
698 :
699 : ! starting the force loop
700 7332466 : DO ilist = 1, nonbonded%nlists
701 7265634 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
702 7265634 : npairs = neighbor_kind_pair%npairs
703 7265634 : IF (npairs == 0) CYCLE
704 9469821 : Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
705 7284498 : istart = neighbor_kind_pair%grp_kind_start(igrp)
706 7284498 : iend = neighbor_kind_pair%grp_kind_end(igrp)
707 7284498 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
708 7284498 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
709 7284498 : list => neighbor_kind_pair%list
710 29137992 : cvi = neighbor_kind_pair%cell_vector
711 7284498 : pot => potparm%pot(ikind, jkind)%pot
712 7284498 : IF (pot%no_mb) CYCLE Kind_Group_Loop1
713 63830 : rab2_max = pot%rcutsq
714 829790 : cell_v = MATMUL(cell%hmat, cvi)
715 94382 : any_tersoff = any_tersoff .OR. ANY(pot%type == tersoff_type)
716 127497 : any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type)
717 111148 : any_ace = any_ace .OR. ANY(pot%type == ace_type)
718 127660 : any_deepmd = any_deepmd .OR. ANY(pot%type == deepmd_type)
719 127653 : any_gal = any_gal .OR. ANY(pot%type == gal_type)
720 127653 : any_gal21 = any_gal21 .OR. ANY(pot%type == gal21_type)
721 127530 : any_nequip = any_nequip .OR. ANY(pot%type == nequip_type)
722 127530 : any_allegro = any_allegro .OR. ANY(pot%type == allegro_type)
723 63830 : i = eam_kinds_index(ikind, jkind)
724 63830 : IF (i == -1) CYCLE Kind_Group_Loop1
725 : ! EAM
726 13535 : CPASSERT(ASSOCIATED(eam_data))
727 7437070 : DO ipair = istart, iend
728 157901 : atom_a = list(1, ipair)
729 157901 : atom_b = list(2, ipair)
730 157901 : fac = 1.0_dp
731 157901 : IF (atom_a == atom_b) fac = 0.5_dp
732 157901 : kind_a = particle_set(atom_a)%atomic_kind%kind_number
733 157901 : kind_b = particle_set(atom_b)%atomic_kind%kind_number
734 157901 : i_a = eam_kinds_index(kind_a, kind_a)
735 157901 : i_b = eam_kinds_index(kind_b, kind_b)
736 157901 : eam_a => potparm%pot(kind_a, kind_a)%pot%set(i_a)%eam
737 157901 : eam_b => potparm%pot(kind_b, kind_b)%pot%set(i_b)%eam
738 :
739 : !set this outside the potential type in case need multiple potentials
740 : !Do everything necessary for EAM here
741 631604 : rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
742 631604 : rab = rab + cell_v
743 157901 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
744 7442399 : IF (rab2 <= rab2_max) THEN
745 97493 : CALL get_force_eam(rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam)
746 97493 : f_eam = f_eam*fac
747 :
748 97493 : fr(1) = -f_eam*rab(1)
749 97493 : fr(2) = -f_eam*rab(2)
750 97493 : fr(3) = -f_eam*rab(3)
751 97493 : f_nonbond(1, atom_a) = f_nonbond(1, atom_a) - fr(1)
752 97493 : f_nonbond(2, atom_a) = f_nonbond(2, atom_a) - fr(2)
753 97493 : f_nonbond(3, atom_a) = f_nonbond(3, atom_a) - fr(3)
754 :
755 97493 : f_nonbond(1, atom_b) = f_nonbond(1, atom_b) + fr(1)
756 97493 : f_nonbond(2, atom_b) = f_nonbond(2, atom_b) + fr(2)
757 97493 : f_nonbond(3, atom_b) = f_nonbond(3, atom_b) + fr(3)
758 97493 : IF (use_virial) THEN
759 4112 : ptens11 = ptens11 + rab(1)*fr(1)
760 4112 : ptens21 = ptens21 + rab(2)*fr(1)
761 4112 : ptens31 = ptens31 + rab(3)*fr(1)
762 4112 : ptens12 = ptens12 + rab(1)*fr(2)
763 4112 : ptens22 = ptens22 + rab(2)*fr(2)
764 4112 : ptens32 = ptens32 + rab(3)*fr(2)
765 4112 : ptens13 = ptens13 + rab(1)*fr(3)
766 4112 : ptens23 = ptens23 + rab(2)*fr(3)
767 4112 : ptens33 = ptens33 + rab(3)*fr(3)
768 : END IF
769 : END IF
770 : END DO
771 : END DO Kind_Group_Loop1
772 : END DO
773 66832 : DEALLOCATE (eam_kinds_index)
774 :
775 : ! Special way of handling the tersoff potential..
776 66832 : IF (any_tersoff) THEN
777 2664 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
778 2664 : CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
779 101052 : DO ilist = 1, nonbonded%nlists
780 98388 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
781 98388 : npairs = neighbor_kind_pair%npairs
782 98388 : IF (npairs == 0) CYCLE
783 68185 : Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
784 33321 : istart = neighbor_kind_pair%grp_kind_start(igrp)
785 33321 : iend = neighbor_kind_pair%grp_kind_end(igrp)
786 33321 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
787 33321 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
788 33321 : list => neighbor_kind_pair%list
789 133284 : cvi = neighbor_kind_pair%cell_vector
790 33321 : pot => potparm%pot(ikind, jkind)%pot
791 :
792 33321 : IF (pot%no_mb) CYCLE Kind_Group_Loop2
793 33285 : rab2_max = pot%rcutsq
794 432705 : cell_v = MATMUL(cell%hmat, cvi)
795 164960 : DO i = 1, SIZE(pot%type)
796 : ! TERSOFF
797 66608 : IF (pot%type(i) == tersoff_type) THEN
798 33280 : npairs = iend - istart + 1
799 33280 : tersoff => pot%set(i)%tersoff
800 166400 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
801 29042564 : sort_list = list(:, istart:iend)
802 : ! Sort the list of neighbors, this increases the efficiency for single
803 : ! potential contributions
804 33280 : CALL sort(sort_list(1, :), npairs, work_list)
805 4862614 : DO ipair = 1, npairs
806 4862614 : work_list(ipair) = sort_list(2, work_list(ipair))
807 : END DO
808 9725228 : sort_list(2, :) = work_list
809 : ! find number of unique elements of array index 1
810 33280 : nunique = 1
811 4829334 : DO ipair = 1, npairs - 1
812 4829334 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
813 : END DO
814 33280 : ipair = 1
815 33280 : junique = sort_list(1, ipair)
816 33280 : ifirst = 1
817 344350 : DO iunique = 1, nunique
818 311070 : atom_a = junique
819 311070 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
820 61587820 : DO mpair = ifirst, SIZE(glob_loc_list_a)
821 61587820 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
822 : END DO
823 70828815 : ifirst = mpair
824 70828815 : DO mpair = ifirst, SIZE(glob_loc_list_a)
825 70828815 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
826 : END DO
827 311070 : ilast = mpair - 1
828 311070 : nloc_size = 0
829 311070 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
830 5140404 : DO WHILE (ipair <= npairs)
831 5107124 : IF (sort_list(1, ipair) /= junique) EXIT
832 4829334 : atom_b = sort_list(2, ipair)
833 : ! Derivative terms
834 19317336 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
835 4829334 : ipair = ipair + 1
836 19628406 : IF (DOT_PRODUCT(rtmp, rtmp) <= tersoff%rcutsq) THEN
837 : CALL tersoff_forces(tersoff, r_last_update_pbc, cell_v, &
838 : nloc_size, glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), &
839 250380 : atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, tersoff%rcutsq)
840 : END IF
841 : END DO
842 311070 : ifirst = ilast + 1
843 344350 : IF (ipair <= npairs) junique = sort_list(1, ipair)
844 : END DO
845 33280 : DEALLOCATE (sort_list, work_list)
846 : END IF
847 : END DO
848 : END DO Kind_Group_Loop2
849 : END DO
850 2664 : CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
851 : END IF
852 : ! Special way of handling the siepmann potential..
853 66832 : IF (any_siepmann) THEN
854 21 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
855 21 : CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
856 588 : DO ilist = 1, nonbonded%nlists
857 567 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
858 567 : npairs = neighbor_kind_pair%npairs
859 567 : IF (npairs == 0) CYCLE
860 918 : Kind_Group_Loop3: DO igrp = 1, neighbor_kind_pair%ngrp_kind
861 708 : istart = neighbor_kind_pair%grp_kind_start(igrp)
862 708 : iend = neighbor_kind_pair%grp_kind_end(igrp)
863 708 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
864 708 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
865 708 : list => neighbor_kind_pair%list
866 2832 : cvi = neighbor_kind_pair%cell_vector
867 708 : pot => potparm%pot(ikind, jkind)%pot
868 :
869 708 : IF (pot%no_mb) CYCLE Kind_Group_Loop3
870 165 : rab2_max = pot%rcutsq
871 2145 : cell_v = MATMUL(cell%hmat, cvi)
872 897 : DO i = 1, SIZE(pot%type)
873 : ! SIEPMANN
874 873 : IF (pot%type(i) == siepmann_type) THEN
875 165 : npairs = iend - istart + 1
876 165 : siepmann => pot%set(i)%siepmann
877 825 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
878 109698 : sort_list = list(:, istart:iend)
879 : ! Sort the list of neighbors, this increases the efficiency for single
880 : ! potential contributions
881 165 : CALL sort(sort_list(1, :), npairs, work_list)
882 18393 : DO ipair = 1, npairs
883 18393 : work_list(ipair) = sort_list(2, work_list(ipair))
884 : END DO
885 36786 : sort_list(2, :) = work_list
886 : ! find number of unique elements of array index 1
887 165 : nunique = 1
888 18228 : DO ipair = 1, npairs - 1
889 18228 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
890 : END DO
891 165 : ipair = 1
892 165 : junique = sort_list(1, ipair)
893 165 : ifirst = 1
894 5340 : DO iunique = 1, nunique
895 5175 : atom_a = junique
896 5175 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
897 91602 : DO mpair = ifirst, SIZE(glob_loc_list_a)
898 91602 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
899 : END DO
900 62187 : ifirst = mpair
901 62187 : DO mpair = ifirst, SIZE(glob_loc_list_a)
902 62187 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
903 : END DO
904 5175 : ilast = mpair - 1
905 5175 : nloc_size = 0
906 5175 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
907 23403 : DO WHILE (ipair <= npairs)
908 23238 : IF (sort_list(1, ipair) /= junique) EXIT
909 18228 : atom_b = sort_list(2, ipair)
910 : ! Derivative terms
911 72912 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
912 18228 : ipair = ipair + 1
913 78087 : IF (DOT_PRODUCT(rtmp, rtmp) <= siepmann%rcutsq) THEN
914 : CALL siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, &
915 : atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, &
916 318 : particle_set)
917 : CALL siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, &
918 : nloc_size, glob_loc_list(:, ifirst:ilast), &
919 : atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, &
920 318 : cell, particle_set)
921 : END IF
922 : END DO
923 5175 : ifirst = ilast + 1
924 5340 : IF (ipair <= npairs) junique = sort_list(1, ipair)
925 : END DO
926 165 : DEALLOCATE (sort_list, work_list)
927 : END IF
928 : END DO
929 : END DO Kind_Group_Loop3
930 : END DO
931 21 : CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
932 : END IF
933 :
934 : ! GAL19 potential..
935 66832 : IF (any_gal) THEN
936 1 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
937 1 : CALL setup_gal_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
938 28 : DO ilist = 1, nonbonded%nlists
939 27 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
940 27 : npairs = neighbor_kind_pair%npairs
941 27 : IF (npairs == 0) CYCLE
942 168 : Kind_Group_Loop4: DO igrp = 1, neighbor_kind_pair%ngrp_kind
943 158 : istart = neighbor_kind_pair%grp_kind_start(igrp)
944 158 : iend = neighbor_kind_pair%grp_kind_end(igrp)
945 158 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
946 158 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
947 158 : list => neighbor_kind_pair%list
948 632 : cvi = neighbor_kind_pair%cell_vector
949 158 : pot => potparm%pot(ikind, jkind)%pot
950 :
951 158 : IF (pot%no_mb) CYCLE Kind_Group_Loop4
952 9 : rab2_max = pot%rcutsq
953 117 : cell_v = MATMUL(cell%hmat, cvi)
954 45 : DO i = 1, SIZE(pot%type)
955 : ! GAL19
956 167 : IF (pot%type(i) == gal_type) THEN
957 9 : npairs = iend - istart + 1
958 9 : gal => pot%set(i)%gal
959 45 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
960 45618 : sort_list = list(:, istart:iend)
961 : ! Sort the list of neighbors, this increases the efficiency for single
962 : ! potential contributions
963 9 : CALL sort(sort_list(1, :), npairs, work_list)
964 7609 : DO ipair = 1, npairs
965 7609 : work_list(ipair) = sort_list(2, work_list(ipair))
966 : END DO
967 15218 : sort_list(2, :) = work_list
968 : ! find number of unique elements of array index 1
969 9 : nunique = 1
970 7600 : DO ipair = 1, npairs - 1
971 7600 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
972 : END DO
973 9 : ipair = 1
974 9 : junique = sort_list(1, ipair)
975 9 : ifirst = 1
976 659 : DO iunique = 1, nunique
977 650 : atom_a = junique
978 650 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
979 36198 : DO mpair = ifirst, SIZE(glob_loc_list_a)
980 36198 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
981 : END DO
982 24581 : ifirst = mpair
983 24581 : DO mpair = ifirst, SIZE(glob_loc_list_a)
984 24581 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
985 : END DO
986 650 : ilast = mpair - 1
987 650 : nloc_size = 0
988 650 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
989 8250 : DO WHILE (ipair <= npairs)
990 8241 : IF (sort_list(1, ipair) /= junique) EXIT
991 7600 : atom_b = sort_list(2, ipair)
992 : ! Derivative terms
993 30400 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
994 7600 : ipair = ipair + 1
995 31050 : IF (DOT_PRODUCT(rtmp, rtmp) <= gal%rcutsq) THEN
996 : CALL gal_forces(gal, r_last_update_pbc, &
997 : atom_a, atom_b, f_nonbond, use_virial, &
998 2004 : cell, particle_set)
999 : END IF
1000 : END DO
1001 650 : ifirst = ilast + 1
1002 659 : IF (ipair <= npairs) junique = sort_list(1, ipair)
1003 : END DO
1004 9 : DEALLOCATE (sort_list, work_list)
1005 : END IF
1006 : END DO
1007 : END DO Kind_Group_Loop4
1008 : END DO
1009 1 : CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
1010 : END IF
1011 :
1012 : ! GAL21 potential..
1013 66832 : IF (any_gal21) THEN
1014 1 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
1015 1 : CALL setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
1016 28 : DO ilist = 1, nonbonded%nlists
1017 27 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
1018 27 : npairs = neighbor_kind_pair%npairs
1019 27 : IF (npairs == 0) CYCLE
1020 168 : Kind_Group_Loop6: DO igrp = 1, neighbor_kind_pair%ngrp_kind
1021 158 : istart = neighbor_kind_pair%grp_kind_start(igrp)
1022 158 : iend = neighbor_kind_pair%grp_kind_end(igrp)
1023 158 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
1024 158 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
1025 158 : list => neighbor_kind_pair%list
1026 632 : cvi = neighbor_kind_pair%cell_vector
1027 158 : pot => potparm%pot(ikind, jkind)%pot
1028 :
1029 158 : IF (pot%no_mb) CYCLE Kind_Group_Loop6
1030 9 : rab2_max = pot%rcutsq
1031 117 : cell_v = MATMUL(cell%hmat, cvi)
1032 45 : DO i = 1, SIZE(pot%type)
1033 : ! GAL21
1034 167 : IF (pot%type(i) == gal21_type) THEN
1035 9 : npairs = iend - istart + 1
1036 9 : gal21 => pot%set(i)%gal21
1037 45 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
1038 52818 : sort_list = list(:, istart:iend)
1039 : ! Sort the list of neighbors, this increases the efficiency for single
1040 : ! potential contributions
1041 9 : CALL sort(sort_list(1, :), npairs, work_list)
1042 8809 : DO ipair = 1, npairs
1043 8809 : work_list(ipair) = sort_list(2, work_list(ipair))
1044 : END DO
1045 17618 : sort_list(2, :) = work_list
1046 : ! find number of unique elements of array index 1
1047 9 : nunique = 1
1048 8800 : DO ipair = 1, npairs - 1
1049 8800 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
1050 : END DO
1051 9 : ipair = 1
1052 9 : junique = sort_list(1, ipair)
1053 9 : ifirst = 1
1054 710 : DO iunique = 1, nunique
1055 701 : atom_a = junique
1056 701 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
1057 42242 : DO mpair = ifirst, SIZE(glob_loc_list_a)
1058 42242 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
1059 : END DO
1060 30069 : ifirst = mpair
1061 30069 : DO mpair = ifirst, SIZE(glob_loc_list_a)
1062 30069 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
1063 : END DO
1064 701 : ilast = mpair - 1
1065 701 : nloc_size = 0
1066 701 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
1067 9501 : DO WHILE (ipair <= npairs)
1068 9492 : IF (sort_list(1, ipair) /= junique) EXIT
1069 8800 : atom_b = sort_list(2, ipair)
1070 : ! Derivative terms
1071 35200 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
1072 8800 : ipair = ipair + 1
1073 35901 : IF (DOT_PRODUCT(rtmp, rtmp) <= gal21%rcutsq) THEN
1074 : CALL gal21_forces(gal21, r_last_update_pbc, &
1075 : atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, &
1076 5732 : cell, particle_set)
1077 : END IF
1078 : END DO
1079 701 : ifirst = ilast + 1
1080 710 : IF (ipair <= npairs) junique = sort_list(1, ipair)
1081 : END DO
1082 9 : DEALLOCATE (sort_list, work_list)
1083 : END IF
1084 : END DO
1085 : END DO Kind_Group_Loop6
1086 : END DO
1087 1 : CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
1088 : END IF
1089 :
1090 66832 : IF (use_virial) THEN
1091 5242 : pv_nonbond(1, 1) = pv_nonbond(1, 1) + ptens11
1092 5242 : pv_nonbond(1, 2) = pv_nonbond(1, 2) + ptens12
1093 5242 : pv_nonbond(1, 3) = pv_nonbond(1, 3) + ptens13
1094 5242 : pv_nonbond(2, 1) = pv_nonbond(2, 1) + ptens21
1095 5242 : pv_nonbond(2, 2) = pv_nonbond(2, 2) + ptens22
1096 5242 : pv_nonbond(2, 3) = pv_nonbond(2, 3) + ptens23
1097 5242 : pv_nonbond(3, 1) = pv_nonbond(3, 1) + ptens31
1098 5242 : pv_nonbond(3, 2) = pv_nonbond(3, 2) + ptens32
1099 5242 : pv_nonbond(3, 3) = pv_nonbond(3, 3) + ptens33
1100 : END IF
1101 66832 : CALL timestop(handle)
1102 66832 : END SUBROUTINE force_nonbond_manybody
1103 :
1104 : END MODULE manybody_potential
1105 :
|
