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