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