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_allegro, ONLY: allegro_add_force_virial,&
32 : allegro_energy_store_force_virial,&
33 : destroy_allegro_arrays,&
34 : setup_allegro_arrays
35 : USE manybody_deepmd, ONLY: deepmd_add_force_virial,&
36 : deepmd_energy_store_force_virial
37 : USE manybody_eam, ONLY: get_force_eam
38 : USE manybody_gal, ONLY: destroy_gal_arrays,&
39 : gal_energy,&
40 : gal_forces,&
41 : setup_gal_arrays
42 : USE manybody_gal21, ONLY: destroy_gal21_arrays,&
43 : gal21_energy,&
44 : gal21_forces,&
45 : setup_gal21_arrays
46 : USE manybody_nequip, ONLY: destroy_nequip_arrays,&
47 : nequip_add_force_virial,&
48 : nequip_energy_store_force_virial,&
49 : setup_nequip_arrays
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 : ace_type, allegro_pot_type, allegro_type, deepmd_type, ea_type, eam_pot_type, &
63 : gal21_pot_type, gal21_type, gal_pot_type, gal_type, nequip_pot_type, nequip_type, &
64 : pair_potential_pp_type, pair_potential_single_type, siepmann_pot_type, siepmann_type, &
65 : tersoff_pot_type, 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 :
77 : CONTAINS
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 : !> \param use_virial ...
90 : !> \par History
91 : !> tlaino [2007] - New algorithm for tersoff potential
92 : !> \author CJM, I-Feng W. Kuo, Teodoro Laino
93 : ! **************************************************************************************************
94 78998 : SUBROUTINE energy_manybody(fist_nonbond_env, atomic_kind_set, local_particles, &
95 : particle_set, cell, pot_manybody, para_env, mm_section, use_virial)
96 :
97 : TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
98 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
99 : TYPE(distribution_1d_type), POINTER :: local_particles
100 : TYPE(particle_type), POINTER :: particle_set(:)
101 : TYPE(cell_type), POINTER :: cell
102 : REAL(dp), INTENT(INOUT) :: pot_manybody
103 : TYPE(mp_para_env_type), OPTIONAL, POINTER :: para_env
104 : TYPE(section_vals_type), POINTER :: mm_section
105 : LOGICAL, INTENT(IN) :: use_virial
106 :
107 : CHARACTER(LEN=*), PARAMETER :: routineN = 'energy_manybody'
108 :
109 : INTEGER :: atom_a, atom_b, handle, i, iend, ifirst, igrp, ikind, ilast, ilist, indexa, &
110 : ipair, iparticle, iparticle_local, istart, iunique, jkind, junique, mpair, nkinds, &
111 : nloc_size, npairs, nparticle, nparticle_local, nr_h3O, nr_o, nr_oh, nunique
112 78998 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, unique_list_a, work_list
113 78998 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list
114 : LOGICAL :: any_ace, any_allegro, any_deepmd, &
115 : any_gal, any_gal21, any_nequip, &
116 : any_siepmann, any_tersoff
117 : REAL(KIND=dp) :: drij, embed, pot_ace, pot_allegro, &
118 : pot_deepmd, pot_loc, pot_nequip, qr, &
119 : rab2_max, rij(3)
120 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi
121 78998 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
122 78998 : REAL(KIND=dp), POINTER :: fembed(:)
123 : TYPE(allegro_pot_type), POINTER :: allegro
124 : TYPE(eam_pot_type), POINTER :: eam
125 78998 : TYPE(eam_type), DIMENSION(:), POINTER :: eam_data
126 : TYPE(fist_neighbor_type), POINTER :: nonbonded
127 : TYPE(gal21_pot_type), POINTER :: gal21
128 : TYPE(gal_pot_type), POINTER :: gal
129 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
130 : TYPE(nequip_pot_type), POINTER :: nequip
131 : TYPE(pair_potential_pp_type), POINTER :: potparm
132 : TYPE(pair_potential_single_type), POINTER :: pot
133 78998 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update_pbc
134 : TYPE(siepmann_pot_type), POINTER :: siepmann
135 : TYPE(tersoff_pot_type), POINTER :: tersoff
136 :
137 78998 : NULLIFY (eam, siepmann, tersoff, gal, gal21)
138 78998 : any_tersoff = .FALSE.
139 78998 : any_siepmann = .FALSE.
140 78998 : any_gal = .FALSE.
141 78998 : any_gal21 = .FALSE.
142 78998 : any_nequip = .FALSE.
143 78998 : any_allegro = .FALSE.
144 78998 : any_ace = .FALSE.
145 78998 : any_deepmd = .FALSE.
146 78998 : CALL timeset(routineN, handle)
147 : CALL fist_nonbond_env_get(fist_nonbond_env, r_last_update_pbc=r_last_update_pbc, &
148 78998 : potparm=potparm, eam_data=eam_data)
149 : ! EAM requires a single loop
150 308790 : DO ikind = 1, SIZE(atomic_kind_set)
151 229792 : pot => potparm%pot(ikind, ikind)%pot
152 538630 : DO i = 1, SIZE(pot%type)
153 229840 : IF (pot%type(i) /= ea_type) CYCLE
154 488 : eam => pot%set(i)%eam
155 488 : nparticle = SIZE(particle_set)
156 1464 : ALLOCATE (fembed(nparticle))
157 14258 : fembed(:) = 0._dp
158 488 : CPASSERT(ASSOCIATED(eam_data))
159 : ! computation of embedding function and energy
160 488 : nparticle_local = local_particles%n_el(ikind)
161 4136 : DO iparticle_local = 1, nparticle_local
162 3648 : iparticle = local_particles%list(ikind)%array(iparticle_local)
163 3648 : indexa = INT(eam_data(iparticle)%rho/eam%drhoar) + 1
164 3648 : IF (indexa > eam%npoints - 1) indexa = eam%npoints - 1
165 3648 : qr = eam_data(iparticle)%rho - eam%rhoval(indexa)
166 :
167 3648 : embed = eam%frho(indexa) + qr*eam%frhop(indexa)
168 3648 : fembed(iparticle) = eam%frhop(indexa) + qr*(eam%frhop(indexa + 1) - eam%frhop(indexa))/eam%drhoar
169 :
170 4136 : pot_manybody = pot_manybody + embed
171 : END DO
172 : ! communicate data
173 28028 : CALL para_env%sum(fembed)
174 14258 : DO iparticle = 1, nparticle
175 14258 : IF (particle_set(iparticle)%atomic_kind%kind_number == ikind) THEN
176 7296 : eam_data(iparticle)%f_embed = fembed(iparticle)
177 : END IF
178 : END DO
179 :
180 459632 : DEALLOCATE (fembed)
181 : END DO
182 : END DO
183 : ! Other manybody potential
184 308790 : DO ikind = 1, SIZE(atomic_kind_set)
185 1780482 : DO jkind = ikind, SIZE(atomic_kind_set)
186 1471692 : pot => potparm%pot(ikind, jkind)%pot
187 2940000 : any_tersoff = any_tersoff .OR. ANY(pot%type == tersoff_type)
188 2943420 : any_nequip = any_nequip .OR. ANY(pot%type == nequip_type)
189 2943424 : any_allegro = any_allegro .OR. ANY(pot%type == allegro_type)
190 2942814 : any_ace = any_ace .OR. ANY(pot%type == ace_type)
191 2943426 : any_deepmd = any_deepmd .OR. ANY(pot%type == deepmd_type)
192 2943411 : any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type)
193 2943431 : any_gal = any_gal .OR. ANY(pot%type == gal_type)
194 3173223 : any_gal21 = any_gal21 .OR. ANY(pot%type == gal21_type)
195 : END DO
196 : END DO
197 78998 : CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, natom_types=nkinds)
198 : ! NEQUIP
199 78998 : IF (any_nequip) THEN
200 4 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
201 4 : CALL setup_nequip_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
202 : CALL nequip_energy_store_force_virial(nonbonded, particle_set, cell, atomic_kind_set, potparm, &
203 : nequip, glob_loc_list_a, r_last_update_pbc, pot_nequip, &
204 4 : fist_nonbond_env, para_env, use_virial)
205 4 : pot_manybody = pot_manybody + pot_nequip
206 4 : CALL destroy_nequip_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
207 : END IF
208 : ! ALLEGRO
209 78998 : IF (any_allegro) THEN
210 4 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a, unique_list_a)
211 : CALL setup_allegro_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, &
212 4 : unique_list_a, cell)
213 : CALL allegro_energy_store_force_virial(nonbonded, particle_set, cell, atomic_kind_set, potparm, &
214 : allegro, glob_loc_list_a, r_last_update_pbc, pot_allegro, &
215 4 : fist_nonbond_env, unique_list_a, para_env, use_virial)
216 4 : pot_manybody = pot_manybody + pot_allegro
217 4 : CALL destroy_allegro_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a, unique_list_a)
218 : END IF
219 : ! ACE
220 78998 : 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 78998 : 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 78998 : IF (any_tersoff) THEN
234 3124 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
235 3124 : CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
236 115500 : DO ilist = 1, nonbonded%nlists
237 112376 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
238 112376 : npairs = neighbor_kind_pair%npairs
239 112376 : IF (npairs == 0) CYCLE
240 87731 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
241 42864 : istart = neighbor_kind_pair%grp_kind_start(igrp)
242 42864 : iend = neighbor_kind_pair%grp_kind_end(igrp)
243 42864 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
244 42864 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
245 42864 : list => neighbor_kind_pair%list
246 171456 : cvi = neighbor_kind_pair%cell_vector
247 42864 : pot => potparm%pot(ikind, jkind)%pot
248 198106 : DO i = 1, SIZE(pot%type)
249 42866 : IF (pot%type(i) /= tersoff_type) CYCLE
250 42823 : rab2_max = pot%set(i)%tersoff%rcutsq
251 556699 : cell_v = MATMUL(cell%hmat, cvi)
252 42823 : pot => potparm%pot(ikind, jkind)%pot
253 42823 : tersoff => pot%set(i)%tersoff
254 42823 : npairs = iend - istart + 1
255 85687 : IF (npairs /= 0) THEN
256 214115 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
257 45922303 : sort_list = list(:, istart:iend)
258 : ! Sort the list of neighbors, this increases the efficiency for single
259 : ! potential contributions
260 42823 : CALL sort(sort_list(1, :), npairs, work_list)
261 7689403 : DO ipair = 1, npairs
262 7689403 : work_list(ipair) = sort_list(2, work_list(ipair))
263 : END DO
264 15335983 : sort_list(2, :) = work_list
265 : ! find number of unique elements of array index 1
266 42823 : nunique = 1
267 7646580 : DO ipair = 1, npairs - 1
268 7646580 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
269 : END DO
270 42823 : ipair = 1
271 42823 : junique = sort_list(1, ipair)
272 42823 : ifirst = 1
273 494916 : DO iunique = 1, nunique
274 452093 : atom_a = junique
275 452093 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
276 93196653 : DO mpair = ifirst, SIZE(glob_loc_list_a)
277 93196653 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
278 : END DO
279 110993591 : ifirst = mpair
280 110993591 : DO mpair = ifirst, SIZE(glob_loc_list_a)
281 110993591 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
282 : END DO
283 452093 : ilast = mpair - 1
284 452093 : nloc_size = 0
285 452093 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
286 8098673 : DO WHILE (ipair <= npairs)
287 8055850 : IF (sort_list(1, ipair) /= junique) EXIT
288 7646580 : atom_b = sort_list(2, ipair)
289 : ! Energy terms
290 7646580 : pot_loc = 0.0_dp
291 30586320 : rij(:) = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
292 30586320 : drij = DOT_PRODUCT(rij, rij)
293 7646580 : ipair = ipair + 1
294 7646580 : IF (drij > rab2_max) CYCLE
295 330372 : drij = SQRT(drij)
296 : CALL tersoff_energy(pot_loc, tersoff, r_last_update_pbc, atom_a, atom_b, nloc_size, &
297 330372 : glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), cell_v, drij)
298 8055850 : pot_manybody = pot_manybody + 0.5_dp*pot_loc
299 : END DO
300 452093 : ifirst = ilast + 1
301 494916 : IF (ipair <= npairs) junique = sort_list(1, ipair)
302 : END DO
303 42823 : DEALLOCATE (sort_list, work_list)
304 : END IF
305 : END DO
306 : END DO Kind_Group_Loop
307 : END DO
308 3124 : CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
309 : END IF
310 :
311 : !SIEPMANN POTENTIAL
312 78998 : 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 78998 : 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 78998 : 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 701 : 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 5732 : 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 78998 : CALL timestop(handle)
561 78998 : 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 69618 : SUBROUTINE force_nonbond_manybody(fist_nonbond_env, particle_set, cell, &
576 69618 : 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 69618 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: eam_kinds_index
590 69618 : INTEGER, DIMENSION(:), POINTER :: glob_loc_list_a, work_list
591 69618 : INTEGER, DIMENSION(:, :), POINTER :: glob_loc_list, list, sort_list
592 : LOGICAL :: any_ace, any_allegro, any_deepmd, &
593 : any_gal, any_gal21, any_nequip, &
594 : 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 69618 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: glob_cell_v
599 : TYPE(eam_pot_type), POINTER :: eam_a, eam_b
600 69618 : 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 69618 : 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 69618 : any_tersoff = .FALSE.
612 69618 : any_nequip = .FALSE.
613 69618 : any_allegro = .FALSE.
614 69618 : any_siepmann = .FALSE.
615 69618 : any_ace = .FALSE.
616 69618 : any_deepmd = .FALSE.
617 69618 : any_gal = .FALSE.
618 69618 : any_gal21 = .FALSE.
619 69618 : CALL timeset(routineN, handle)
620 69618 : NULLIFY (eam_a, eam_b, tersoff, siepmann, gal, gal21)
621 :
622 : CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, potparm=potparm, &
623 69618 : natom_types=nkinds, eam_data=eam_data, r_last_update_pbc=r_last_update_pbc)
624 :
625 : ! Initializing the potential energy, pressure tensor and force
626 : IF (use_virial) THEN
627 : ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
628 : ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
629 : ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
630 : END IF
631 :
632 69618 : nkinds = SIZE(potparm%pot, 1)
633 278472 : ALLOCATE (eam_kinds_index(nkinds, nkinds))
634 2956782 : eam_kinds_index = -1
635 280668 : DO ikind = 1, nkinds
636 1724250 : DO jkind = ikind, nkinds
637 3098262 : DO i = 1, SIZE(potparm%pot(ikind, jkind)%pot%type)
638 2887212 : IF (potparm%pot(ikind, jkind)%pot%type(i) == ea_type) THEN
639 : ! At the moment we allow only 1 EAM per each kinds pair..
640 692 : CPASSERT(eam_kinds_index(ikind, jkind) == -1)
641 692 : CPASSERT(eam_kinds_index(jkind, ikind) == -1)
642 692 : eam_kinds_index(ikind, jkind) = i
643 692 : eam_kinds_index(jkind, ikind) = i
644 : END IF
645 : END DO
646 : END DO
647 : END DO
648 280668 : DO ikind = 1, nkinds
649 1724250 : DO jkind = ikind, nkinds
650 3097644 : any_ace = any_ace .OR. ANY(potparm%pot(ikind, jkind)%pot%type == ace_type)
651 : END DO
652 : END DO
653 : ! ACE
654 69618 : IF (any_ace) &
655 206 : CALL ace_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
656 :
657 280668 : DO ikind = 1, nkinds
658 1724250 : DO jkind = ikind, nkinds
659 3098256 : any_deepmd = any_deepmd .OR. ANY(potparm%pot(ikind, jkind)%pot%type == deepmd_type)
660 : END DO
661 : END DO
662 : ! DEEPMD
663 69618 : IF (any_deepmd) &
664 2 : CALL deepmd_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
665 :
666 : ! starting the force loop
667 8176358 : DO ilist = 1, nonbonded%nlists
668 8106740 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
669 8106740 : npairs = neighbor_kind_pair%npairs
670 8106740 : IF (npairs == 0) CYCLE
671 10164424 : Kind_Group_Loop1: DO igrp = 1, neighbor_kind_pair%ngrp_kind
672 7772534 : istart = neighbor_kind_pair%grp_kind_start(igrp)
673 7772534 : iend = neighbor_kind_pair%grp_kind_end(igrp)
674 7772534 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
675 7772534 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
676 7772534 : list => neighbor_kind_pair%list
677 31090136 : cvi = neighbor_kind_pair%cell_vector
678 7772534 : pot => potparm%pot(ikind, jkind)%pot
679 7772534 : IF (pot%no_mb) CYCLE Kind_Group_Loop1
680 73342 : rab2_max = pot%rcutsq
681 953446 : cell_v = MATMUL(cell%hmat, cvi)
682 103863 : any_tersoff = any_tersoff .OR. ANY(pot%type == tersoff_type)
683 146521 : any_siepmann = any_siepmann .OR. ANY(pot%type == siepmann_type)
684 130172 : any_ace = any_ace .OR. ANY(pot%type == ace_type)
685 146684 : any_deepmd = any_deepmd .OR. ANY(pot%type == deepmd_type)
686 146677 : any_gal = any_gal .OR. ANY(pot%type == gal_type)
687 146677 : any_gal21 = any_gal21 .OR. ANY(pot%type == gal21_type)
688 146570 : any_nequip = any_nequip .OR. ANY(pot%type == nequip_type)
689 146517 : any_allegro = any_allegro .OR. ANY(pot%type == allegro_type)
690 73342 : i = eam_kinds_index(ikind, jkind)
691 73342 : IF (i == -1) CYCLE Kind_Group_Loop1
692 : ! EAM
693 13535 : CPASSERT(ASSOCIATED(eam_data))
694 8278176 : DO ipair = istart, iend
695 157901 : atom_a = list(1, ipair)
696 157901 : atom_b = list(2, ipair)
697 157901 : fac = 1.0_dp
698 157901 : IF (atom_a == atom_b) fac = 0.5_dp
699 157901 : kind_a = particle_set(atom_a)%atomic_kind%kind_number
700 157901 : kind_b = particle_set(atom_b)%atomic_kind%kind_number
701 157901 : i_a = eam_kinds_index(kind_a, kind_a)
702 157901 : i_b = eam_kinds_index(kind_b, kind_b)
703 157901 : eam_a => potparm%pot(kind_a, kind_a)%pot%set(i_a)%eam
704 157901 : eam_b => potparm%pot(kind_b, kind_b)%pot%set(i_b)%eam
705 :
706 : !set this outside the potential type in case need multiple potentials
707 : !Do everything necessary for EAM here
708 631604 : rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
709 631604 : rab = rab + cell_v
710 157901 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
711 7930435 : IF (rab2 <= rab2_max) THEN
712 97493 : CALL get_force_eam(rab2, eam_a, eam_b, eam_data, atom_a, atom_b, f_eam)
713 97493 : f_eam = f_eam*fac
714 :
715 97493 : fr(1) = -f_eam*rab(1)
716 97493 : fr(2) = -f_eam*rab(2)
717 97493 : fr(3) = -f_eam*rab(3)
718 97493 : f_nonbond(1, atom_a) = f_nonbond(1, atom_a) - fr(1)
719 97493 : f_nonbond(2, atom_a) = f_nonbond(2, atom_a) - fr(2)
720 97493 : f_nonbond(3, atom_a) = f_nonbond(3, atom_a) - fr(3)
721 :
722 97493 : f_nonbond(1, atom_b) = f_nonbond(1, atom_b) + fr(1)
723 97493 : f_nonbond(2, atom_b) = f_nonbond(2, atom_b) + fr(2)
724 97493 : f_nonbond(3, atom_b) = f_nonbond(3, atom_b) + fr(3)
725 97493 : IF (use_virial) THEN
726 4112 : ptens11 = ptens11 + rab(1)*fr(1)
727 4112 : ptens21 = ptens21 + rab(2)*fr(1)
728 4112 : ptens31 = ptens31 + rab(3)*fr(1)
729 4112 : ptens12 = ptens12 + rab(1)*fr(2)
730 4112 : ptens22 = ptens22 + rab(2)*fr(2)
731 4112 : ptens32 = ptens32 + rab(3)*fr(2)
732 4112 : ptens13 = ptens13 + rab(1)*fr(3)
733 4112 : ptens23 = ptens23 + rab(2)*fr(3)
734 4112 : ptens33 = ptens33 + rab(3)*fr(3)
735 : END IF
736 : END IF
737 : END DO
738 : END DO Kind_Group_Loop1
739 : END DO
740 69618 : DEALLOCATE (eam_kinds_index)
741 :
742 : ! Special way of handling the tersoff potential..
743 69618 : IF (any_tersoff) THEN
744 3124 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
745 3124 : CALL setup_tersoff_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
746 115500 : DO ilist = 1, nonbonded%nlists
747 112376 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
748 112376 : npairs = neighbor_kind_pair%npairs
749 112376 : IF (npairs == 0) CYCLE
750 87731 : Kind_Group_Loop2: DO igrp = 1, neighbor_kind_pair%ngrp_kind
751 42864 : istart = neighbor_kind_pair%grp_kind_start(igrp)
752 42864 : iend = neighbor_kind_pair%grp_kind_end(igrp)
753 42864 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
754 42864 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
755 42864 : list => neighbor_kind_pair%list
756 171456 : cvi = neighbor_kind_pair%cell_vector
757 42864 : pot => potparm%pot(ikind, jkind)%pot
758 :
759 42864 : IF (pot%no_mb) CYCLE Kind_Group_Loop2
760 42828 : rab2_max = pot%rcutsq
761 556764 : cell_v = MATMUL(cell%hmat, cvi)
762 198034 : DO i = 1, SIZE(pot%type)
763 : ! TERSOFF
764 85694 : IF (pot%type(i) == tersoff_type) THEN
765 42823 : npairs = iend - istart + 1
766 42823 : tersoff => pot%set(i)%tersoff
767 214115 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
768 45965126 : sort_list = list(:, istart:iend)
769 : ! Sort the list of neighbors, this increases the efficiency for single
770 : ! potential contributions
771 42823 : CALL sort(sort_list(1, :), npairs, work_list)
772 7689403 : DO ipair = 1, npairs
773 7689403 : work_list(ipair) = sort_list(2, work_list(ipair))
774 : END DO
775 15378806 : sort_list(2, :) = work_list
776 : ! find number of unique elements of array index 1
777 42823 : nunique = 1
778 7646580 : DO ipair = 1, npairs - 1
779 7646580 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
780 : END DO
781 42823 : ipair = 1
782 42823 : junique = sort_list(1, ipair)
783 42823 : ifirst = 1
784 494916 : DO iunique = 1, nunique
785 452093 : atom_a = junique
786 452093 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
787 93196653 : DO mpair = ifirst, SIZE(glob_loc_list_a)
788 93196653 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
789 : END DO
790 110993591 : ifirst = mpair
791 110993591 : DO mpair = ifirst, SIZE(glob_loc_list_a)
792 110993591 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
793 : END DO
794 452093 : ilast = mpair - 1
795 452093 : nloc_size = 0
796 452093 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
797 8098673 : DO WHILE (ipair <= npairs)
798 8055850 : IF (sort_list(1, ipair) /= junique) EXIT
799 7646580 : atom_b = sort_list(2, ipair)
800 : ! Derivative terms
801 30586320 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
802 7646580 : ipair = ipair + 1
803 31038413 : IF (DOT_PRODUCT(rtmp, rtmp) <= tersoff%rcutsq) THEN
804 : CALL tersoff_forces(tersoff, r_last_update_pbc, cell_v, &
805 : nloc_size, glob_loc_list(:, ifirst:ilast), glob_cell_v(:, ifirst:ilast), &
806 330372 : atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, tersoff%rcutsq)
807 : END IF
808 : END DO
809 452093 : ifirst = ilast + 1
810 494916 : IF (ipair <= npairs) junique = sort_list(1, ipair)
811 : END DO
812 42823 : DEALLOCATE (sort_list, work_list)
813 : END IF
814 : END DO
815 : END DO Kind_Group_Loop2
816 : END DO
817 3124 : CALL destroy_tersoff_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
818 : END IF
819 : ! Special way of handling the siepmann potential..
820 69618 : IF (any_siepmann) THEN
821 21 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
822 21 : CALL setup_siepmann_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
823 588 : DO ilist = 1, nonbonded%nlists
824 567 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
825 567 : npairs = neighbor_kind_pair%npairs
826 567 : IF (npairs == 0) CYCLE
827 918 : Kind_Group_Loop3: DO igrp = 1, neighbor_kind_pair%ngrp_kind
828 708 : istart = neighbor_kind_pair%grp_kind_start(igrp)
829 708 : iend = neighbor_kind_pair%grp_kind_end(igrp)
830 708 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
831 708 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
832 708 : list => neighbor_kind_pair%list
833 2832 : cvi = neighbor_kind_pair%cell_vector
834 708 : pot => potparm%pot(ikind, jkind)%pot
835 :
836 708 : IF (pot%no_mb) CYCLE Kind_Group_Loop3
837 165 : rab2_max = pot%rcutsq
838 2145 : cell_v = MATMUL(cell%hmat, cvi)
839 897 : DO i = 1, SIZE(pot%type)
840 : ! SIEPMANN
841 873 : IF (pot%type(i) == siepmann_type) THEN
842 165 : npairs = iend - istart + 1
843 165 : siepmann => pot%set(i)%siepmann
844 825 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
845 109698 : sort_list = list(:, istart:iend)
846 : ! Sort the list of neighbors, this increases the efficiency for single
847 : ! potential contributions
848 165 : CALL sort(sort_list(1, :), npairs, work_list)
849 18393 : DO ipair = 1, npairs
850 18393 : work_list(ipair) = sort_list(2, work_list(ipair))
851 : END DO
852 36786 : sort_list(2, :) = work_list
853 : ! find number of unique elements of array index 1
854 165 : nunique = 1
855 18228 : DO ipair = 1, npairs - 1
856 18228 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
857 : END DO
858 165 : ipair = 1
859 165 : junique = sort_list(1, ipair)
860 165 : ifirst = 1
861 5340 : DO iunique = 1, nunique
862 5175 : atom_a = junique
863 5175 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
864 91602 : DO mpair = ifirst, SIZE(glob_loc_list_a)
865 91602 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
866 : END DO
867 62187 : ifirst = mpair
868 62187 : DO mpair = ifirst, SIZE(glob_loc_list_a)
869 62187 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
870 : END DO
871 5175 : ilast = mpair - 1
872 5175 : nloc_size = 0
873 5175 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
874 23403 : DO WHILE (ipair <= npairs)
875 23238 : IF (sort_list(1, ipair) /= junique) EXIT
876 18228 : atom_b = sort_list(2, ipair)
877 : ! Derivative terms
878 72912 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
879 18228 : ipair = ipair + 1
880 78087 : IF (DOT_PRODUCT(rtmp, rtmp) <= siepmann%rcutsq) THEN
881 : CALL siepmann_forces_v2(siepmann, r_last_update_pbc, cell_v, cell, &
882 : atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, &
883 318 : particle_set)
884 : CALL siepmann_forces_v3(siepmann, r_last_update_pbc, cell_v, &
885 : nloc_size, glob_loc_list(:, ifirst:ilast), &
886 : atom_a, atom_b, f_nonbond, use_virial, siepmann%rcutsq, &
887 318 : cell, particle_set)
888 : END IF
889 : END DO
890 5175 : ifirst = ilast + 1
891 5340 : IF (ipair <= npairs) junique = sort_list(1, ipair)
892 : END DO
893 165 : DEALLOCATE (sort_list, work_list)
894 : END IF
895 : END DO
896 : END DO Kind_Group_Loop3
897 : END DO
898 21 : CALL destroy_siepmann_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
899 : END IF
900 :
901 : ! GAL19 potential..
902 69618 : IF (any_gal) THEN
903 1 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
904 1 : CALL setup_gal_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
905 28 : DO ilist = 1, nonbonded%nlists
906 27 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
907 27 : npairs = neighbor_kind_pair%npairs
908 27 : IF (npairs == 0) CYCLE
909 168 : Kind_Group_Loop4: DO igrp = 1, neighbor_kind_pair%ngrp_kind
910 158 : istart = neighbor_kind_pair%grp_kind_start(igrp)
911 158 : iend = neighbor_kind_pair%grp_kind_end(igrp)
912 158 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
913 158 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
914 158 : list => neighbor_kind_pair%list
915 632 : cvi = neighbor_kind_pair%cell_vector
916 158 : pot => potparm%pot(ikind, jkind)%pot
917 :
918 158 : IF (pot%no_mb) CYCLE Kind_Group_Loop4
919 9 : rab2_max = pot%rcutsq
920 117 : cell_v = MATMUL(cell%hmat, cvi)
921 45 : DO i = 1, SIZE(pot%type)
922 : ! GAL19
923 167 : IF (pot%type(i) == gal_type) THEN
924 9 : npairs = iend - istart + 1
925 9 : gal => pot%set(i)%gal
926 45 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
927 45618 : sort_list = list(:, istart:iend)
928 : ! Sort the list of neighbors, this increases the efficiency for single
929 : ! potential contributions
930 9 : CALL sort(sort_list(1, :), npairs, work_list)
931 7609 : DO ipair = 1, npairs
932 7609 : work_list(ipair) = sort_list(2, work_list(ipair))
933 : END DO
934 15218 : sort_list(2, :) = work_list
935 : ! find number of unique elements of array index 1
936 9 : nunique = 1
937 7600 : DO ipair = 1, npairs - 1
938 7600 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
939 : END DO
940 9 : ipair = 1
941 9 : junique = sort_list(1, ipair)
942 9 : ifirst = 1
943 659 : DO iunique = 1, nunique
944 650 : atom_a = junique
945 650 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
946 36198 : DO mpair = ifirst, SIZE(glob_loc_list_a)
947 36198 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
948 : END DO
949 24581 : ifirst = mpair
950 24581 : DO mpair = ifirst, SIZE(glob_loc_list_a)
951 24581 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
952 : END DO
953 650 : ilast = mpair - 1
954 650 : nloc_size = 0
955 650 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
956 8250 : DO WHILE (ipair <= npairs)
957 8241 : IF (sort_list(1, ipair) /= junique) EXIT
958 7600 : atom_b = sort_list(2, ipair)
959 : ! Derivative terms
960 30400 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
961 7600 : ipair = ipair + 1
962 31050 : IF (DOT_PRODUCT(rtmp, rtmp) <= gal%rcutsq) THEN
963 : CALL gal_forces(gal, r_last_update_pbc, &
964 : atom_a, atom_b, f_nonbond, use_virial, &
965 2004 : cell, particle_set)
966 : END IF
967 : END DO
968 650 : ifirst = ilast + 1
969 659 : IF (ipair <= npairs) junique = sort_list(1, ipair)
970 : END DO
971 9 : DEALLOCATE (sort_list, work_list)
972 : END IF
973 : END DO
974 : END DO Kind_Group_Loop4
975 : END DO
976 1 : CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
977 : END IF
978 :
979 : ! GAL21 potential..
980 69618 : IF (any_gal21) THEN
981 1 : NULLIFY (glob_loc_list, glob_cell_v, glob_loc_list_a)
982 1 : CALL setup_gal21_arrays(nonbonded, potparm, glob_loc_list, glob_cell_v, glob_loc_list_a, cell)
983 28 : DO ilist = 1, nonbonded%nlists
984 27 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
985 27 : npairs = neighbor_kind_pair%npairs
986 27 : IF (npairs == 0) CYCLE
987 168 : Kind_Group_Loop6: DO igrp = 1, neighbor_kind_pair%ngrp_kind
988 158 : istart = neighbor_kind_pair%grp_kind_start(igrp)
989 158 : iend = neighbor_kind_pair%grp_kind_end(igrp)
990 158 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
991 158 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
992 158 : list => neighbor_kind_pair%list
993 632 : cvi = neighbor_kind_pair%cell_vector
994 158 : pot => potparm%pot(ikind, jkind)%pot
995 :
996 158 : IF (pot%no_mb) CYCLE Kind_Group_Loop6
997 9 : rab2_max = pot%rcutsq
998 117 : cell_v = MATMUL(cell%hmat, cvi)
999 45 : DO i = 1, SIZE(pot%type)
1000 : ! GAL21
1001 167 : IF (pot%type(i) == gal21_type) THEN
1002 9 : npairs = iend - istart + 1
1003 9 : gal21 => pot%set(i)%gal21
1004 45 : ALLOCATE (sort_list(2, npairs), work_list(npairs))
1005 52818 : sort_list = list(:, istart:iend)
1006 : ! Sort the list of neighbors, this increases the efficiency for single
1007 : ! potential contributions
1008 9 : CALL sort(sort_list(1, :), npairs, work_list)
1009 8809 : DO ipair = 1, npairs
1010 8809 : work_list(ipair) = sort_list(2, work_list(ipair))
1011 : END DO
1012 17618 : sort_list(2, :) = work_list
1013 : ! find number of unique elements of array index 1
1014 9 : nunique = 1
1015 8800 : DO ipair = 1, npairs - 1
1016 8800 : IF (sort_list(1, ipair + 1) /= sort_list(1, ipair)) nunique = nunique + 1
1017 : END DO
1018 9 : ipair = 1
1019 9 : junique = sort_list(1, ipair)
1020 9 : ifirst = 1
1021 710 : DO iunique = 1, nunique
1022 701 : atom_a = junique
1023 701 : IF (glob_loc_list_a(ifirst) > atom_a) CYCLE
1024 42242 : DO mpair = ifirst, SIZE(glob_loc_list_a)
1025 42242 : IF (glob_loc_list_a(mpair) == atom_a) EXIT
1026 : END DO
1027 30069 : ifirst = mpair
1028 30069 : DO mpair = ifirst, SIZE(glob_loc_list_a)
1029 30069 : IF (glob_loc_list_a(mpair) /= atom_a) EXIT
1030 : END DO
1031 701 : ilast = mpair - 1
1032 701 : nloc_size = 0
1033 701 : IF (ifirst /= 0) nloc_size = ilast - ifirst + 1
1034 9501 : DO WHILE (ipair <= npairs)
1035 9492 : IF (sort_list(1, ipair) /= junique) EXIT
1036 8800 : atom_b = sort_list(2, ipair)
1037 : ! Derivative terms
1038 35200 : rtmp = r_last_update_pbc(atom_b)%r(:) - r_last_update_pbc(atom_a)%r(:) + cell_v
1039 8800 : ipair = ipair + 1
1040 35901 : IF (DOT_PRODUCT(rtmp, rtmp) <= gal21%rcutsq) THEN
1041 : CALL gal21_forces(gal21, r_last_update_pbc, &
1042 : atom_a, atom_b, f_nonbond, pv_nonbond, use_virial, &
1043 5732 : cell, particle_set)
1044 : END IF
1045 : END DO
1046 701 : ifirst = ilast + 1
1047 710 : IF (ipair <= npairs) junique = sort_list(1, ipair)
1048 : END DO
1049 9 : DEALLOCATE (sort_list, work_list)
1050 : END IF
1051 : END DO
1052 : END DO Kind_Group_Loop6
1053 : END DO
1054 1 : CALL destroy_gal_arrays(glob_loc_list, glob_cell_v, glob_loc_list_a)
1055 : END IF
1056 :
1057 : ! NEQUIP
1058 69618 : IF (any_nequip) THEN
1059 4 : CALL nequip_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
1060 : END IF
1061 :
1062 : ! ALLEGRO
1063 69618 : IF (any_allegro) THEN
1064 4 : CALL allegro_add_force_virial(fist_nonbond_env, f_nonbond, pv_nonbond, use_virial)
1065 : END IF
1066 :
1067 69618 : IF (use_virial) THEN
1068 9344 : pv_nonbond(1, 1) = pv_nonbond(1, 1) + ptens11
1069 9344 : pv_nonbond(1, 2) = pv_nonbond(1, 2) + ptens12
1070 9344 : pv_nonbond(1, 3) = pv_nonbond(1, 3) + ptens13
1071 9344 : pv_nonbond(2, 1) = pv_nonbond(2, 1) + ptens21
1072 9344 : pv_nonbond(2, 2) = pv_nonbond(2, 2) + ptens22
1073 9344 : pv_nonbond(2, 3) = pv_nonbond(2, 3) + ptens23
1074 9344 : pv_nonbond(3, 1) = pv_nonbond(3, 1) + ptens31
1075 9344 : pv_nonbond(3, 2) = pv_nonbond(3, 2) + ptens32
1076 9344 : pv_nonbond(3, 3) = pv_nonbond(3, 3) + ptens33
1077 : END IF
1078 69618 : CALL timestop(handle)
1079 69618 : END SUBROUTINE force_nonbond_manybody
1080 :
1081 : END MODULE manybody_potential
1082 :
|
