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