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