(git:a0df33d)
Loading...
Searching...
No Matches
fist_nonbond_force.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9!> \par History
10!> JGH (11 May 2001) : cleaning up of support structures
11!> CJM & HAF (27 July 2001): fixed bug with handling of cutoff larger than
12!> half the boxsize.
13!> 07.02.2005: getting rid of scaled_to_real calls in force loop (MK)
14!> 22.06.2013: OpenMP parallelisation of pair interaction loop (MK)
15!> \author CJM
16! **************************************************************************************************
17MODULE fist_nonbond_force
18 USE atomic_kind_types, ONLY: atomic_kind_type,&
19 get_atomic_kind,&
20 get_atomic_kind_set
21 USE atprop_types, ONLY: atprop_type
22 USE cell_types, ONLY: cell_type,&
23 pbc
24 USE cp_log_handling, ONLY: cp_get_default_logger,&
25 cp_logger_type
26 USE distribution_1d_types, ONLY: distribution_1d_type
27 USE ewald_environment_types, ONLY: ewald_env_get,&
28 ewald_environment_type
29 USE fist_neighbor_list_types, ONLY: fist_neighbor_type,&
30 neighbor_kind_pairs_type
31 USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get,&
32 fist_nonbond_env_type,&
33 pos_type
34 USE kinds, ONLY: dp
35 USE machine, ONLY: m_memory
36 USE mathconstants, ONLY: oorootpi,&
37 sqrthalf
38 USE message_passing, ONLY: mp_comm_type
39 USE pair_potential_coulomb, ONLY: potential_coulomb
40 USE pair_potential_types, ONLY: &
41 ace_type, allegro_type, deepmd_type, gal21_type, gal_type, mace_type, nequip_type, &
42 nosh_nosh, nosh_sh, pair_potential_pp_type, pair_potential_single_type, sh_sh, &
43 siepmann_type, tersoff_type
44 USE particle_types, ONLY: particle_type
45 USE shell_potential_types, ONLY: get_shell,&
46 shell_kind_type
47 USE splines_methods, ONLY: potential_s
48 USE splines_types, ONLY: spline_data_p_type,&
49 spline_factor_type
50#include "./base/base_uses.f90"
51
52 IMPLICIT NONE
53
54 PRIVATE
55
56 CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'fist_nonbond_force'
57 LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
58
59 PUBLIC :: force_nonbond, &
60 bonded_correct_gaussian
61
62CONTAINS
63
64! **************************************************************************************************
65!> \brief Calculates the force and the potential of the minimum image, and
66!> the pressure tensor
67!> \param fist_nonbond_env ...
68!> \param ewald_env ...
69!> \param particle_set ...
70!> \param cell ...
71!> \param pot_nonbond ...
72!> \param f_nonbond ...
73!> \param pv_nonbond ...
74!> \param fshell_nonbond ...
75!> \param fcore_nonbond ...
76!> \param atprop_env ...
77!> \param atomic_kind_set ...
78!> \param use_virial ...
79! **************************************************************************************************
80 SUBROUTINE force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
81 pot_nonbond, f_nonbond, pv_nonbond, fshell_nonbond, fcore_nonbond, &
82 atprop_env, atomic_kind_set, use_virial)
83
84 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
85 TYPE(ewald_environment_type), POINTER :: ewald_env
86 TYPE(particle_type), DIMENSION(:), INTENT(IN) :: particle_set
87 TYPE(cell_type), POINTER :: cell
88 REAL(kind=dp), INTENT(OUT) :: pot_nonbond
89 REAL(kind=dp), DIMENSION(:, :), INTENT(INOUT) :: f_nonbond, pv_nonbond
90 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT), &
91 OPTIONAL :: fshell_nonbond, fcore_nonbond
92 TYPE(atprop_type), POINTER :: atprop_env
93 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
94 LOGICAL, INTENT(IN) :: use_virial
95
96 CHARACTER(LEN=*), PARAMETER :: routinen = 'force_nonbond'
97
98 INTEGER :: atom_a, atom_b, ewald_type, handle, i, iend, igrp, ikind, ilist, ipair, istart, &
99 j, kind_a, kind_b, nkind, npairs, shell_a, shell_b, shell_type
100 INTEGER, DIMENSION(:, :), POINTER :: list
101 LOGICAL :: all_terms, do_multipoles, full_nl, &
102 shell_present
103 LOGICAL, ALLOCATABLE, DIMENSION(:) :: is_shell_kind
104 REAL(kind=dp) :: alpha, beta, beta_a, beta_b, energy, etot, fac_ei, fac_kind, fac_vdw, &
105 fscalar, mm_radius_a, mm_radius_b, qcore_a, qcore_b, qeff_a, qeff_b, qshell_a, qshell_b, &
106 rab2, rab2_com, rab2_max
107 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: mm_radius, qcore, qeff, qshell
108 REAL(kind=dp), DIMENSION(3) :: cell_v, cvi, fatom_a, fatom_b, fcore_a, &
109 fcore_b, fshell_a, fshell_b, rab, &
110 rab_cc, rab_com, rab_cs, rab_sc, rab_ss
111 REAL(kind=dp), DIMENSION(3, 3) :: pv, pv_thread
112 REAL(kind=dp), DIMENSION(3, 4) :: rab_list
113 REAL(kind=dp), DIMENSION(4) :: rab2_list
114 REAL(kind=dp), DIMENSION(:, :), POINTER :: ij_kind_full_fac
115 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: ei_interaction_cutoffs
116 TYPE(atomic_kind_type), POINTER :: atomic_kind
117 TYPE(cp_logger_type), POINTER :: logger
118 TYPE(fist_neighbor_type), POINTER :: nonbonded
119 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
120 TYPE(pair_potential_pp_type), POINTER :: potparm, potparm14
121 TYPE(pair_potential_single_type), POINTER :: pot
122 TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc, &
123 rcore_last_update_pbc, &
124 rshell_last_update_pbc
125 TYPE(shell_kind_type), POINTER :: shell_kind
126 TYPE(spline_data_p_type), DIMENSION(:), POINTER :: spline_data
127 TYPE(spline_factor_type), POINTER :: spl_f
128
129 CALL timeset(routinen, handle)
130 NULLIFY (logger)
131 logger => cp_get_default_logger()
132 NULLIFY (pot, rshell_last_update_pbc, spl_f, ij_kind_full_fac)
133 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, &
134 potparm14=potparm14, potparm=potparm, r_last_update=r_last_update, &
135 r_last_update_pbc=r_last_update_pbc, natom_types=nkind, &
136 rshell_last_update_pbc=rshell_last_update_pbc, &
137 rcore_last_update_pbc=rcore_last_update_pbc, &
138 ij_kind_full_fac=ij_kind_full_fac)
139 CALL ewald_env_get(ewald_env, alpha=alpha, ewald_type=ewald_type, &
140 do_multipoles=do_multipoles, &
141 interaction_cutoffs=ei_interaction_cutoffs)
142
143 ! Initializing the potential energy, pressure tensor and force
144 pot_nonbond = 0.0_dp
145 f_nonbond(:, :) = 0.0_dp
146
147 IF (use_virial) THEN
148 pv_nonbond(:, :) = 0.0_dp
149 END IF
150 shell_present = .false.
151 IF (PRESENT(fshell_nonbond)) THEN
152 cpassert(PRESENT(fcore_nonbond))
153 fshell_nonbond = 0.0_dp
154 fcore_nonbond = 0.0_dp
155 shell_present = .true.
156 END IF
157 ! Load atomic kind information
158 ALLOCATE (mm_radius(nkind))
159 ALLOCATE (qeff(nkind))
160 ALLOCATE (qcore(nkind))
161 ALLOCATE (qshell(nkind))
162 ALLOCATE (is_shell_kind(nkind))
163 DO ikind = 1, nkind
164 atomic_kind => atomic_kind_set(ikind)
165 CALL get_atomic_kind(atomic_kind, &
166 qeff=qeff(ikind), &
167 mm_radius=mm_radius(ikind), &
168 shell=shell_kind)
169 is_shell_kind(ikind) = ASSOCIATED(shell_kind)
170 IF (ASSOCIATED(shell_kind)) THEN
171 CALL get_shell(shell=shell_kind, &
172 charge_core=qcore(ikind), &
173 charge_shell=qshell(ikind))
174 ELSE
175 qcore(ikind) = 0.0_dp
176 qshell(ikind) = 0.0_dp
177 END IF
178 END DO
179 ! Starting the force loop
180 lists: DO ilist = 1, nonbonded%nlists
181 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
182 npairs = neighbor_kind_pair%npairs
183 IF (npairs == 0) cycle lists
184 list => neighbor_kind_pair%list
185 cvi = neighbor_kind_pair%cell_vector
186 cell_v = matmul(cell%hmat, cvi)
187 kind_group_loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
188 istart = neighbor_kind_pair%grp_kind_start(igrp)
189 iend = neighbor_kind_pair%grp_kind_end(igrp)
190!$OMP PARALLEL DEFAULT(NONE) &
191!$OMP PRIVATE(ipair,atom_a,atom_b,kind_a,kind_b,fac_kind,pot) &
192!$OMP PRIVATE(fac_ei,fac_vdw,atomic_kind,full_nl,qcore_a,qshell_a) &
193!$OMP PRIVATE(qeff_a,qcore_b,qshell_b,qeff_b,mm_radius_a,mm_radius_b) &
194!$OMP PRIVATE(shell_kind,beta,beta_a,beta_b,spl_f,spline_data) &
195!$OMP PRIVATE(shell_type,all_terms,rab_cc,rab_cs,rab_sc,rab_ss) &
196!$OMP PRIVATE(rab_list,rab2_list,rab_com,rab2_com,pv,pv_thread) &
197!$OMP PRIVATE(rab,rab2,rab2_max,fscalar,energy) &
198!$OMP PRIVATE(shell_a,shell_b,etot,fatom_a,fatom_b) &
199!$OMP PRIVATE(fcore_a,fcore_b,fshell_a,fshell_b,i,j) &
200!$OMP SHARED(shell_present) &
201!$OMP SHARED(istart,iend,list,particle_set,ij_kind_full_fac) &
202!$OMP SHARED(neighbor_kind_pair,atomic_kind_set,fist_nonbond_env) &
203!$OMP SHARED(potparm,potparm14,do_multipoles,r_last_update_pbc) &
204!$OMP SHARED(use_virial,ei_interaction_cutoffs,alpha,cell_v) &
205!$OMP SHARED(rcore_last_update_pbc,rshell_last_update_pbc) &
206!$OMP SHARED(f_nonbond,fcore_nonbond,fshell_nonbond,logger) &
207!$OMP SHARED(ewald_type,pot_nonbond,pv_nonbond,atprop_env) &
208!$OMP SHARED(is_shell_kind,mm_radius,qcore,qeff,qshell)
209 IF (use_virial) pv_thread(:, :) = 0.0_dp
210!$OMP DO
211 pairs: DO ipair = istart, iend
212 atom_a = list(1, ipair)
213 atom_b = list(2, ipair)
214 ! Get actual atomic kinds, since atom_a is not always of
215 ! kind_a and atom_b of kind_b, ie. they might be swapped.
216 kind_a = particle_set(atom_a)%atomic_kind%kind_number
217 kind_b = particle_set(atom_b)%atomic_kind%kind_number
218
219 fac_kind = ij_kind_full_fac(kind_a, kind_b)
220 ! take the proper potential
221 pot => potparm%pot(kind_a, kind_b)%pot
222 IF (ipair <= neighbor_kind_pair%nscale) THEN
223 IF (neighbor_kind_pair%is_onfo(ipair)) THEN
224 pot => potparm14%pot(kind_a, kind_b)%pot
225 END IF
226 END IF
227
228 ! Determine the scaling factors
229 fac_ei = fac_kind
230 fac_vdw = fac_kind
231 full_nl = any(pot%type == tersoff_type) .OR. any(pot%type == siepmann_type) &
232 .OR. any(pot%type == gal_type) .OR. any(pot%type == gal21_type) &
233 .OR. any(pot%type == nequip_type) .OR. any(pot%type == allegro_type) &
234 .OR. any(pot%type == mace_type) &
235 .OR. any(pot%type == ace_type) .OR. any(pot%type == deepmd_type)
236 IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
237 fac_ei = 0.5_dp*fac_ei
238 fac_vdw = 0.5_dp*fac_vdw
239 END IF
240 ! decide which interactions to compute\b
241 IF (do_multipoles .OR. (.NOT. fist_nonbond_env%do_electrostatics)) THEN
242 fac_ei = 0.0_dp
243 END IF
244 IF (ipair <= neighbor_kind_pair%nscale) THEN
245 fac_ei = fac_ei*neighbor_kind_pair%ei_scale(ipair)
246 fac_vdw = fac_vdw*neighbor_kind_pair%vdw_scale(ipair)
247 END IF
248
249 IF (fac_ei > 0.0_dp) THEN
250 ! Get the electrostatic parameters for the atoms a and b
251 mm_radius_a = mm_radius(kind_a)
252 mm_radius_b = mm_radius(kind_b)
253 IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
254 qeff_a = fist_nonbond_env%charges(atom_a)
255 qeff_b = fist_nonbond_env%charges(atom_b)
256 ELSE
257 qeff_a = qeff(kind_a)
258 qeff_b = qeff(kind_b)
259 END IF
260 IF (is_shell_kind(kind_a)) THEN
261 qcore_a = qcore(kind_a)
262 qshell_a = qshell(kind_a)
263 IF ((qcore_a == 0.0_dp) .AND. (qshell_a == 0.0_dp)) fac_ei = 0.0_dp
264 ELSE
265 qcore_a = qeff_a
266 qshell_a = huge(0.0_dp)
267 IF (qeff_a == 0.0_dp) fac_ei = 0.0_dp
268 END IF
269 IF (is_shell_kind(kind_b)) THEN
270 qcore_b = qcore(kind_b)
271 qshell_b = qshell(kind_b)
272 IF ((qcore_b == 0.0_dp) .AND. (qshell_b == 0.0_dp)) fac_ei = 0.0_dp
273 ELSE
274 qcore_b = qeff_b
275 qshell_b = huge(0.0_dp)
276 IF (qeff_b == 0.0_dp) fac_ei = 0.0_dp
277 END IF
278 ! Derive beta parameters
279 beta = 0.0_dp
280 beta_a = 0.0_dp
281 beta_b = 0.0_dp
282 IF (mm_radius_a > 0) THEN
283 beta_a = sqrthalf/mm_radius_a
284 END IF
285 IF (mm_radius_b > 0) THEN
286 beta_b = sqrthalf/mm_radius_b
287 END IF
288 IF ((mm_radius_a > 0) .OR. (mm_radius_b > 0)) THEN
289 beta = sqrthalf/sqrt(mm_radius_a*mm_radius_a + mm_radius_b*mm_radius_b)
290 END IF
291 END IF
292
293 ! In case we have only manybody potentials and no charges, this
294 ! pair of atom types can be ignored here.
295 IF (pot%no_pp .AND. (fac_ei == 0.0)) cycle pairs
296
297 ! Setup spline_data set
298 spl_f => pot%spl_f
299 spline_data => pot%pair_spline_data
300 shell_type = pot%shell_type
301 IF (shell_type /= nosh_nosh) THEN
302 cpassert(.NOT. do_multipoles)
303 cpassert(shell_present)
304 END IF
305 rab2_max = pot%rcutsq
306
307 ! compute the relative vector(s) for this pair
308 IF (shell_type /= nosh_nosh) THEN
309 ! do shell
310 all_terms = .true.
311 IF (shell_type == sh_sh) THEN
312 shell_a = particle_set(atom_a)%shell_index
313 shell_b = particle_set(atom_b)%shell_index
314 rab_cc = rcore_last_update_pbc(shell_b)%r - rcore_last_update_pbc(shell_a)%r
315 rab_cs = rshell_last_update_pbc(shell_b)%r - rcore_last_update_pbc(shell_a)%r
316 rab_sc = rcore_last_update_pbc(shell_b)%r - rshell_last_update_pbc(shell_a)%r
317 rab_ss = rshell_last_update_pbc(shell_b)%r - rshell_last_update_pbc(shell_a)%r
318 rab_list(1:3, 1) = rab_cc(1:3) + cell_v(1:3)
319 rab_list(1:3, 2) = rab_cs(1:3) + cell_v(1:3)
320 rab_list(1:3, 3) = rab_sc(1:3) + cell_v(1:3)
321 rab_list(1:3, 4) = rab_ss(1:3) + cell_v(1:3)
322 ELSE IF ((shell_type == nosh_sh) .AND. (particle_set(atom_a)%shell_index /= 0)) THEN
323 shell_a = particle_set(atom_a)%shell_index
324 shell_b = 0
325 rab_cc = r_last_update_pbc(atom_b)%r - rcore_last_update_pbc(shell_a)%r
326 rab_sc = 0.0_dp
327 rab_cs = 0.0_dp
328 rab_ss = r_last_update_pbc(atom_b)%r - rshell_last_update_pbc(shell_a)%r
329 rab_list(1:3, 1) = rab_cc(1:3) + cell_v(1:3)
330 rab_list(1:3, 2) = 0.0_dp
331 rab_list(1:3, 3) = 0.0_dp
332 rab_list(1:3, 4) = rab_ss(1:3) + cell_v(1:3)
333 ELSE IF ((shell_type == nosh_sh) .AND. (particle_set(atom_b)%shell_index /= 0)) THEN
334 shell_b = particle_set(atom_b)%shell_index
335 shell_a = 0
336 rab_cc = rcore_last_update_pbc(shell_b)%r - r_last_update_pbc(atom_a)%r
337 rab_sc = 0.0_dp
338 rab_cs = 0.0_dp
339 rab_ss = rshell_last_update_pbc(shell_b)%r - r_last_update_pbc(atom_a)%r
340 rab_list(1:3, 1) = rab_cc(1:3) + cell_v(1:3)
341 rab_list(1:3, 2) = 0.0_dp
342 rab_list(1:3, 3) = 0.0_dp
343 rab_list(1:3, 4) = rab_ss(1:3) + cell_v(1:3)
344 ELSE
345 rab_list(:, :) = 0.0_dp
346 END IF
347 ! Compute the term only if all the pairs (cc,cs,sc,ss) are within the cut-off
348 check_terms: DO i = 1, 4
349 rab2_list(i) = rab_list(1, i)**2 + rab_list(2, i)**2 + rab_list(3, i)**2
350 IF (rab2_list(i) >= rab2_max) THEN
351 all_terms = .false.
352 EXIT check_terms
353 END IF
354 END DO check_terms
355 rab_com = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
356 ELSE
357 ! not do shell
358 rab_cc = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
359 rab_com = rab_cc
360 shell_a = 0
361 shell_b = 0
362 rab_list(:, :) = 0.0_dp
363 END IF
364 rab_com = rab_com + cell_v
365 rab2_com = rab_com(1)**2 + rab_com(2)**2 + rab_com(3)**2
366
367 ! compute the interactions for the current pair
368 etot = 0.0_dp
369 fatom_a(:) = 0.0_dp
370 fatom_b(:) = 0.0_dp
371 fcore_a(:) = 0.0_dp
372 fcore_b(:) = 0.0_dp
373 fshell_a(:) = 0.0_dp
374 fshell_b(:) = 0.0_dp
375 IF (use_virial) pv(:, :) = 0.0_dp
376 IF (shell_type /= nosh_nosh) THEN
377 ! do shell
378 IF ((rab2_com <= rab2_max) .AND. all_terms) THEN
379 IF (fac_ei > 0) THEN
380 ! core-core or core-ion/ion-core: Coulomb only
381 rab = rab_list(:, 1)
382 rab2 = rab2_list(1)
383 fscalar = 0.0_dp
384 IF (shell_a == 0) THEN
385 ! atom a is a plain ion and can have beta_a > 0
386 energy = potential_coulomb(rab2, fscalar, fac_ei*qeff_a*qcore_b, &
387 ewald_type, alpha, beta_a, &
388 ei_interaction_cutoffs(2, kind_a, kind_b))
389 CALL add_force_nonbond(fatom_a, fcore_b, pv, fscalar, rab, use_virial)
390 ELSE IF (shell_b == 0) THEN
391 ! atom b is a plain ion and can have beta_b > 0
392 energy = potential_coulomb(rab2, fscalar, fac_ei*qcore_a*qeff_b, &
393 ewald_type, alpha, beta_b, &
394 ei_interaction_cutoffs(2, kind_b, kind_a))
395 CALL add_force_nonbond(fcore_a, fatom_b, pv, fscalar, rab, use_virial)
396 ELSE
397 ! core-core interaction is always pure point charge
398 energy = potential_coulomb(rab2, fscalar, fac_ei*qcore_a*qcore_b, &
399 ewald_type, alpha, 0.0_dp, &
400 ei_interaction_cutoffs(1, kind_a, kind_b))
401 CALL add_force_nonbond(fcore_a, fcore_b, pv, fscalar, rab, use_virial)
402 END IF
403 etot = etot + energy
404 END IF
405
406 IF (shell_type == sh_sh) THEN
407 ! shell-shell: VDW + Coulomb
408 rab = rab_list(:, 4)
409 rab2 = rab2_list(4)
410 fscalar = 0.0_dp
411 IF (fac_vdw > 0) THEN
412 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
413 etot = etot + energy*fac_vdw
414 fscalar = fscalar*fac_vdw
415 END IF
416 IF (fac_ei > 0) THEN
417 ! note that potential_coulomb increments fscalar
418 energy = potential_coulomb(rab2, fscalar, fac_ei*qshell_a*qshell_b, &
419 ewald_type, alpha, beta, &
420 ei_interaction_cutoffs(3, kind_a, kind_b))
421 etot = etot + energy
422 END IF
423 CALL add_force_nonbond(fshell_a, fshell_b, pv, fscalar, rab, use_virial)
424
425 IF (fac_ei > 0) THEN
426 ! core-shell: Coulomb only
427 rab = rab_list(:, 2)
428 rab2 = rab2_list(2)
429 fscalar = 0.0_dp
430 ! swap kind_a and kind_b to get the right cutoff
431 energy = potential_coulomb(rab2, fscalar, fac_ei*qcore_a*qshell_b, &
432 ewald_type, alpha, beta_b, &
433 ei_interaction_cutoffs(2, kind_b, kind_a))
434 etot = etot + energy
435 CALL add_force_nonbond(fcore_a, fshell_b, pv, fscalar, rab, use_virial)
436
437 ! shell-core: Coulomb only
438 rab = rab_list(:, 3)
439 rab2 = rab2_list(3)
440 fscalar = 0.0_dp
441 energy = potential_coulomb(rab2, fscalar, fac_ei*qshell_a*qcore_b, &
442 ewald_type, alpha, beta_a, &
443 ei_interaction_cutoffs(2, kind_a, kind_b))
444 etot = etot + energy
445 CALL add_force_nonbond(fshell_a, fcore_b, pv, fscalar, rab, use_virial)
446 END IF
447 ELSE IF ((shell_type == nosh_sh) .AND. (shell_a == 0)) THEN
448 ! ion-shell: VDW + Coulomb
449 rab = rab_list(:, 4)
450 rab2 = rab2_list(4)
451 fscalar = 0.0_dp
452 IF (fac_vdw > 0) THEN
453 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
454 etot = etot + energy*fac_vdw
455 fscalar = fscalar*fac_vdw
456 END IF
457 IF (fac_ei > 0) THEN
458 ! note that potential_coulomb increments fscalar
459 energy = potential_coulomb(rab2, fscalar, fac_ei*qeff_a*qshell_b, &
460 ewald_type, alpha, beta, &
461 ei_interaction_cutoffs(3, kind_a, kind_b))
462 etot = etot + energy
463 END IF
464 CALL add_force_nonbond(fatom_a, fshell_b, pv, fscalar, rab, use_virial)
465 ELSE IF ((shell_type == nosh_sh) .AND. (shell_b == 0)) THEN
466 ! shell-ion : VDW + Coulomb
467 rab = rab_list(:, 4)
468 rab2 = rab2_list(4)
469 fscalar = 0.0_dp
470 IF (fac_vdw > 0) THEN
471 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
472 etot = etot + energy*fac_vdw
473 fscalar = fscalar*fac_vdw
474 END IF
475 IF (fac_ei > 0) THEN
476 ! note that potential_coulomb increments fscalar
477 energy = potential_coulomb(rab2, fscalar, fac_ei*qshell_a*qeff_b, &
478 ewald_type, alpha, beta, &
479 ei_interaction_cutoffs(3, kind_a, kind_b))
480 etot = etot + energy
481 END IF
482 CALL add_force_nonbond(fshell_a, fatom_b, pv, fscalar, rab, use_virial)
483 END IF
484 END IF
485 ELSE
486 IF (rab2_com <= rab2_max) THEN
487 ! NO SHELL MODEL...
488 ! Ion-Ion: no shell model, VDW + coulomb
489 rab = rab_com
490 rab2 = rab2_com
491 fscalar = 0.0_dp
492 IF (fac_vdw > 0) THEN
493 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
494 etot = etot + energy*fac_vdw
495 fscalar = fscalar*fac_vdw
496 END IF
497 IF (fac_ei > 0) THEN
498 ! note that potential_coulomb increments fscalar
499 energy = potential_coulomb(rab2, fscalar, fac_ei*qeff_a*qeff_b, &
500 ewald_type, alpha, beta, &
501 ei_interaction_cutoffs(3, kind_a, kind_b))
502 etot = etot + energy
503 END IF
504 CALL add_force_nonbond(fatom_a, fatom_b, pv, fscalar, rab, use_virial)
505 END IF
506 END IF
507 ! Nonbonded energy
508!$OMP ATOMIC
509 pot_nonbond = pot_nonbond + etot
510 IF (atprop_env%energy) THEN
511 ! Update atomic energies
512!$OMP ATOMIC
513 atprop_env%atener(atom_a) = atprop_env%atener(atom_a) + 0.5_dp*etot
514!$OMP ATOMIC
515 atprop_env%atener(atom_b) = atprop_env%atener(atom_b) + 0.5_dp*etot
516 END IF
517 ! Nonbonded forces
518 DO i = 1, 3
519!$OMP ATOMIC
520 f_nonbond(i, atom_a) = f_nonbond(i, atom_a) + fatom_a(i)
521!$OMP ATOMIC
522 f_nonbond(i, atom_b) = f_nonbond(i, atom_b) + fatom_b(i)
523 END DO
524 IF (shell_a > 0) THEN
525 DO i = 1, 3
526!$OMP ATOMIC
527 fcore_nonbond(i, shell_a) = fcore_nonbond(i, shell_a) + fcore_a(i)
528!$OMP ATOMIC
529 fshell_nonbond(i, shell_a) = fshell_nonbond(i, shell_a) + fshell_a(i)
530 END DO
531 END IF
532 IF (shell_b > 0) THEN
533 DO i = 1, 3
534!$OMP ATOMIC
535 fcore_nonbond(i, shell_b) = fcore_nonbond(i, shell_b) + fcore_b(i)
536!$OMP ATOMIC
537 fshell_nonbond(i, shell_b) = fshell_nonbond(i, shell_b) + fshell_b(i)
538 END DO
539 END IF
540 ! Add the contribution of the current pair to the total pressure tensor
541 IF (use_virial) THEN
542 DO i = 1, 3
543 DO j = 1, 3
544 pv_thread(j, i) = pv_thread(j, i) + pv(j, i)
545 END DO
546 END DO
547 END IF
548 END DO pairs
549!$OMP END DO
550 IF (use_virial) THEN
551 DO i = 1, 3
552 DO j = 1, 3
553!$OMP ATOMIC
554 pv_nonbond(j, i) = pv_nonbond(j, i) + pv_thread(j, i)
555 END DO
556 END DO
557 END IF
558!$OMP END PARALLEL
559 END DO kind_group_loop
560 END DO lists
561
562 !sample peak memory
563 CALL m_memory()
564
565 DEALLOCATE (mm_radius)
566 DEALLOCATE (qeff)
567 DEALLOCATE (qcore)
568 DEALLOCATE (qshell)
569 DEALLOCATE (is_shell_kind)
570
571 CALL timestop(handle)
572
573 END SUBROUTINE force_nonbond
574
575 ! **************************************************************************************************
576 !> \brief Adds a non-bonding contribution to the total force and optionally to
577 !> the virial.
578 ! **************************************************************************************************
579! **************************************************************************************************
580!> \brief ...
581!> \param f_nonbond_a ...
582!> \param f_nonbond_b ...
583!> \param pv ...
584!> \param fscalar ...
585!> \param rab ...
586!> \param use_virial ...
587! **************************************************************************************************
588 SUBROUTINE add_force_nonbond(f_nonbond_a, f_nonbond_b, pv, fscalar, rab, use_virial)
589
590 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: f_nonbond_a, f_nonbond_b
591 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv
592 REAL(kind=dp), INTENT(IN) :: fscalar
593 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
594 LOGICAL, INTENT(IN) :: use_virial
595
596 REAL(kind=dp), DIMENSION(3) :: fr
597
598 fr(1) = fscalar*rab(1)
599 fr(2) = fscalar*rab(2)
600 fr(3) = fscalar*rab(3)
601 f_nonbond_a(1) = f_nonbond_a(1) - fr(1)
602 f_nonbond_a(2) = f_nonbond_a(2) - fr(2)
603 f_nonbond_a(3) = f_nonbond_a(3) - fr(3)
604 f_nonbond_b(1) = f_nonbond_b(1) + fr(1)
605 f_nonbond_b(2) = f_nonbond_b(2) + fr(2)
606 f_nonbond_b(3) = f_nonbond_b(3) + fr(3)
607 IF (use_virial) THEN
608 pv(1, 1) = pv(1, 1) + rab(1)*fr(1)
609 pv(1, 2) = pv(1, 2) + rab(1)*fr(2)
610 pv(1, 3) = pv(1, 3) + rab(1)*fr(3)
611 pv(2, 1) = pv(2, 1) + rab(2)*fr(1)
612 pv(2, 2) = pv(2, 2) + rab(2)*fr(2)
613 pv(2, 3) = pv(2, 3) + rab(2)*fr(3)
614 pv(3, 1) = pv(3, 1) + rab(3)*fr(1)
615 pv(3, 2) = pv(3, 2) + rab(3)*fr(2)
616 pv(3, 3) = pv(3, 3) + rab(3)*fr(3)
617 END IF
618
619 END SUBROUTINE add_force_nonbond
620
621! **************************************************************************************************
622!> \brief corrects electrostatics for bonded terms
623!> \param fist_nonbond_env ...
624!> \param atomic_kind_set ...
625!> \param local_particles ...
626!> \param particle_set ...
627!> \param ewald_env ...
628!> \param v_bonded_corr ...
629!> \param pv_bc ...
630!> \param shell_particle_set ...
631!> \param core_particle_set ...
632!> \param atprop_env ...
633!> \param cell ...
634!> \param use_virial ...
635!> \par History
636!> Split routines to clean and to fix a bug with the tensor whose
637!> original definition was not correct for PBC.. [Teodoro Laino -06/2007]
638! **************************************************************************************************
639 SUBROUTINE bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
640 local_particles, particle_set, ewald_env, v_bonded_corr, pv_bc, &
641 shell_particle_set, core_particle_set, atprop_env, cell, use_virial)
642
643 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
644 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
645 TYPE(distribution_1d_type), POINTER :: local_particles
646 TYPE(particle_type), POINTER :: particle_set(:)
647 TYPE(ewald_environment_type), POINTER :: ewald_env
648 REAL(kind=dp), INTENT(OUT) :: v_bonded_corr
649 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: pv_bc
650 TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
651 core_particle_set(:)
652 TYPE(atprop_type), POINTER :: atprop_env
653 TYPE(cell_type), POINTER :: cell
654 LOGICAL, INTENT(IN) :: use_virial
655
656 CHARACTER(LEN=*), PARAMETER :: routinen = 'bonded_correct_gaussian'
657
658 INTEGER :: atom_a, atom_b, handle, iatom, iend, igrp, ilist, ipair, istart, kind_a, kind_b, &
659 natoms_per_kind, nkind, npairs, shell_a, shell_b
660 INTEGER, DIMENSION(:, :), POINTER :: list
661 LOGICAL :: a_is_shell, b_is_shell, do_multipoles, &
662 full_nl, shell_adiabatic
663 REAL(kind=dp) :: alpha, const, fac_cor, fac_ei, qcore_a, &
664 qcore_b, qeff_a, qeff_b, qshell_a, &
665 qshell_b
666 REAL(kind=dp), DIMENSION(3) :: rca, rcb, rsa, rsb
667 REAL(kind=dp), DIMENSION(:, :), POINTER :: ij_kind_full_fac
668 TYPE(atomic_kind_type), POINTER :: atomic_kind
669 TYPE(fist_neighbor_type), POINTER :: nonbonded
670 TYPE(mp_comm_type) :: group
671 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
672 TYPE(pair_potential_pp_type), POINTER :: potparm, potparm14
673 TYPE(pair_potential_single_type), POINTER :: pot
674 TYPE(shell_kind_type), POINTER :: shell_kind
675
676 CALL timeset(routinen, handle)
677
678 ! Initializing values
679 IF (use_virial) pv_bc = 0.0_dp
680 v_bonded_corr = 0.0_dp
681
682 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, &
683 potparm14=potparm14, potparm=potparm, &
684 ij_kind_full_fac=ij_kind_full_fac)
685 CALL ewald_env_get(ewald_env, alpha=alpha, do_multipoles=do_multipoles, &
686 group=group)
687 ! Defining the constants
688 const = 2.0_dp*alpha*oorootpi
689
690 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
691 shell_adiabatic=shell_adiabatic)
692
693 lists: DO ilist = 1, nonbonded%nlists
694 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
695 npairs = neighbor_kind_pair%nscale
696 IF (npairs == 0) cycle lists
697 list => neighbor_kind_pair%list
698 kind_group_loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
699 istart = neighbor_kind_pair%grp_kind_start(igrp)
700 IF (istart > npairs) THEN
701 EXIT kind_group_loop
702 END IF
703 iend = min(npairs, neighbor_kind_pair%grp_kind_end(igrp))
704
705 pairs: DO ipair = istart, iend
706 atom_a = list(1, ipair)
707 atom_b = list(2, ipair)
708 ! Get actual atomic kinds, since atom_a is not always of
709 ! kind_a and atom_b of kind_b, ie. they might be swapped.
710 kind_a = particle_set(atom_a)%atomic_kind%kind_number
711 kind_b = particle_set(atom_b)%atomic_kind%kind_number
712
713 ! take the proper potential, only for full_nl test
714 pot => potparm%pot(kind_a, kind_b)%pot
715 IF (ipair <= neighbor_kind_pair%nscale) THEN
716 IF (neighbor_kind_pair%is_onfo(ipair)) THEN
717 pot => potparm14%pot(kind_a, kind_b)%pot
718 END IF
719 END IF
720
721 ! Determine the scaling factors
722 fac_ei = ij_kind_full_fac(kind_a, kind_b)
723 full_nl = any(pot%type == tersoff_type) .OR. any(pot%type == siepmann_type) &
724 .OR. any(pot%type == gal_type) .OR. any(pot%type == gal21_type) &
725 .OR. any(pot%type == nequip_type) .OR. any(pot%type == allegro_type) &
726 .OR. any(pot%type == mace_type) &
727 .OR. any(pot%type == ace_type) .OR. any(pot%type == deepmd_type)
728 IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
729 fac_ei = fac_ei*0.5_dp
730 END IF
731 IF (ipair <= neighbor_kind_pair%nscale) THEN
732 fac_ei = fac_ei*neighbor_kind_pair%ei_scale(ipair)
733 END IF
734 ! The amount of correction is related to the
735 ! amount of scaling as follows:
736 fac_cor = 1.0_dp - fac_ei
737 IF (fac_cor <= 0.0_dp) cycle pairs
738
739 ! Parameters for kind a
740 atomic_kind => atomic_kind_set(kind_a)
741 CALL get_atomic_kind(atomic_kind, qeff=qeff_a, shell=shell_kind)
742 IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_a = fist_nonbond_env%charges(atom_a)
743 a_is_shell = ASSOCIATED(shell_kind)
744 IF (a_is_shell) THEN
745 CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
746 charge_shell=qshell_a)
747 shell_a = particle_set(atom_a)%shell_index
748 rca = core_particle_set(shell_a)%r
749 rsa = shell_particle_set(shell_a)%r
750 ELSE
751 qcore_a = qeff_a
752 qshell_a = huge(0.0_dp)
753 shell_a = 0
754 rca = particle_set(atom_a)%r
755 rsa = 0.0_dp
756 END IF
757
758 ! Parameters for kind b
759 atomic_kind => atomic_kind_set(kind_b)
760 CALL get_atomic_kind(atomic_kind, qeff=qeff_b, shell=shell_kind)
761 IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_b = fist_nonbond_env%charges(atom_b)
762 b_is_shell = ASSOCIATED(shell_kind)
763 IF (b_is_shell) THEN
764 CALL get_shell(shell=shell_kind, charge_core=qcore_b, &
765 charge_shell=qshell_b)
766 shell_b = particle_set(atom_b)%shell_index
767 rcb = core_particle_set(shell_b)%r
768 rsb = shell_particle_set(shell_b)%r
769 ELSE
770 qcore_b = qeff_b
771 qshell_b = huge(0.0_dp)
772 shell_b = 0
773 rcb = particle_set(atom_b)%r
774 rsb = 0.0_dp
775 END IF
776
777 ! First part: take care of core/ion-core/ion correction
778 IF (a_is_shell .AND. b_is_shell) THEN
779 ! correct for core-core interaction
780 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
781 v_bonded_corr, core_particle_set, core_particle_set, &
782 shell_a, shell_b, .true., alpha, qcore_a, qcore_b, &
783 const, fac_cor, pv_bc, atprop_env, use_virial)
784 ELSE IF (a_is_shell) THEN
785 ! correct for core-ion interaction
786 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
787 v_bonded_corr, core_particle_set, particle_set, &
788 shell_a, atom_b, .true., alpha, qcore_a, qcore_b, &
789 const, fac_cor, pv_bc, atprop_env, use_virial)
790 ELSE IF (b_is_shell) THEN
791 ! correct for ion-core interaction
792 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
793 v_bonded_corr, particle_set, core_particle_set, &
794 atom_a, shell_b, .true., alpha, qcore_a, qcore_b, &
795 const, fac_cor, pv_bc, atprop_env, use_virial)
796 ELSE
797 ! correct for ion-ion interaction
798 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
799 v_bonded_corr, particle_set, particle_set, &
800 atom_a, atom_b, .true., alpha, qcore_a, qcore_b, &
801 const, fac_cor, pv_bc, atprop_env, use_virial)
802 END IF
803
804 ! Second part: take care of shell-shell, shell-core/ion and
805 ! core/ion-shell corrections
806 IF (a_is_shell .AND. b_is_shell) THEN
807 ! correct for shell-shell interaction
808 CALL bonded_correct_gaussian_low(rsa, rsa, cell, &
809 v_bonded_corr, shell_particle_set, shell_particle_set, &
810 shell_a, shell_b, shell_adiabatic, alpha, qshell_a, &
811 qshell_b, const, fac_cor, pv_bc, atprop_env, use_virial)
812 END IF
813 IF (a_is_shell) THEN
814 IF (b_is_shell) THEN
815 ! correct for shell-core interaction
816 CALL bonded_correct_gaussian_low(rsa, rcb, cell, &
817 v_bonded_corr, shell_particle_set, core_particle_set, &
818 shell_a, shell_b, shell_adiabatic, alpha, qshell_a, qcore_b, &
819 const, fac_cor, pv_bc, atprop_env, use_virial)
820 ELSE
821 ! correct for shell-ion interaction
822 CALL bonded_correct_gaussian_low(rsa, rcb, cell, &
823 v_bonded_corr, shell_particle_set, particle_set, &
824 shell_a, atom_b, shell_adiabatic, alpha, qshell_a, qcore_b, &
825 const, fac_cor, pv_bc, atprop_env, use_virial)
826 END IF
827 END IF
828 IF (b_is_shell) THEN
829 IF (a_is_shell) THEN
830 ! correct for core-shell interaction
831 CALL bonded_correct_gaussian_low(rca, rsb, cell, &
832 v_bonded_corr, core_particle_set, shell_particle_set, &
833 shell_a, shell_b, shell_adiabatic, alpha, qcore_a, qshell_b, &
834 const, fac_cor, pv_bc, atprop_env, use_virial)
835 ELSE
836 ! correct for ion-shell interaction
837 CALL bonded_correct_gaussian_low(rca, rsb, cell, &
838 v_bonded_corr, particle_set, shell_particle_set, &
839 atom_a, shell_b, shell_adiabatic, alpha, qcore_a, qshell_b, &
840 const, fac_cor, pv_bc, atprop_env, use_virial)
841 END IF
842 END IF
843 END DO pairs
844 END DO kind_group_loop
845 END DO lists
846
847 ! Always correct core-shell interaction within one atom.
848 nkind = SIZE(atomic_kind_set)
849 DO kind_a = 1, nkind
850 ! parameters for kind a
851 atomic_kind => atomic_kind_set(kind_a)
852 CALL get_atomic_kind(atomic_kind, shell=shell_kind)
853 IF (ASSOCIATED(shell_kind)) THEN
854 CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
855 charge_shell=qshell_a)
856
857 natoms_per_kind = local_particles%n_el(kind_a)
858 DO iatom = 1, natoms_per_kind
859
860 ! Data for atom a
861 atom_a = local_particles%list(kind_a)%array(iatom)
862 shell_a = particle_set(atom_a)%shell_index
863 rca = core_particle_set(shell_a)%r
864 rsa = shell_particle_set(shell_a)%r
865
866 CALL bonded_correct_gaussian_low_sh(rca, rsa, cell, &
867 v_bonded_corr, core_particle_set, shell_particle_set, &
868 shell_a, shell_adiabatic, alpha, qcore_a, qshell_a, &
869 const, pv_bc, atprop_env, use_virial)
870
871 END DO
872 END IF
873 END DO
874
875 CALL group%sum(v_bonded_corr)
876
877 CALL timestop(handle)
878
879 END SUBROUTINE bonded_correct_gaussian
880
881! **************************************************************************************************
882!> \brief ...
883!> \param r1 ...
884!> \param r2 ...
885!> \param cell ...
886!> \param v_bonded_corr ...
887!> \param particle_set1 ...
888!> \param particle_set2 ...
889!> \param i ...
890!> \param j ...
891!> \param shell_adiabatic ...
892!> \param alpha ...
893!> \param q1 ...
894!> \param q2 ...
895!> \param const ...
896!> \param fac_cor ...
897!> \param pv_bc ...
898!> \param atprop_env ...
899!> \param use_virial ...
900!> \par History
901!> Split routines to clean and to fix a bug with the tensor whose
902!> original definition was not correct for PBC..
903!> \author Teodoro Laino
904! **************************************************************************************************
905 SUBROUTINE bonded_correct_gaussian_low(r1, r2, cell, v_bonded_corr, &
906 particle_set1, particle_set2, i, j, shell_adiabatic, alpha, q1, q2, &
907 const, fac_cor, pv_bc, atprop_env, use_virial)
908 REAL(kind=dp), DIMENSION(3) :: r1, r2
909 TYPE(cell_type), POINTER :: cell
910 REAL(kind=dp), INTENT(INOUT) :: v_bonded_corr
911 TYPE(particle_type), POINTER :: particle_set1(:), particle_set2(:)
912 INTEGER, INTENT(IN) :: i, j
913 LOGICAL, INTENT(IN) :: shell_adiabatic
914 REAL(kind=dp), INTENT(IN) :: alpha, q1, q2, const, fac_cor
915 REAL(kind=dp), INTENT(INOUT) :: pv_bc(3, 3)
916 TYPE(atprop_type), POINTER :: atprop_env
917 LOGICAL, INTENT(IN) :: use_virial
918
919 REAL(kind=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
920 ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, pc = 0.3275911_dp
921
922 INTEGER :: iatom, jatom
923 REAL(kind=dp) :: arg, dij, e_arg_arg, errf, fscalar, &
924 idij, rijsq, tc, vbc
925 REAL(kind=dp), DIMENSION(3) :: fij_com, rij
926 REAL(kind=dp), DIMENSION(3, 3) :: fbc
927
928 rij = r1 - r2
929 rij = pbc(rij, cell)
930 rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
931 idij = 1.0_dp/sqrt(rijsq)
932 dij = rijsq*idij
933 arg = alpha*dij
934 e_arg_arg = exp(-arg**2)
935 tc = 1.0_dp/(1.0_dp + pc*arg)
936
937 ! Defining errf=1-erfc
938 errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
939
940 ! Getting the potential
941 vbc = -q1*q2*idij*errf*fac_cor
942 v_bonded_corr = v_bonded_corr + vbc
943 IF (atprop_env%energy) THEN
944 iatom = particle_set1(i)%atom_index
945 atprop_env%atener(iatom) = atprop_env%atener(iatom) + 0.5_dp*vbc
946 jatom = particle_set2(j)%atom_index
947 atprop_env%atener(jatom) = atprop_env%atener(jatom) + 0.5_dp*vbc
948 END IF
949
950 ! Subtracting the force from the total force
951 fscalar = q1*q2*idij**2*(idij*errf - const*e_arg_arg)*fac_cor
952
953 particle_set1(i)%f(1) = particle_set1(i)%f(1) - fscalar*rij(1)
954 particle_set1(i)%f(2) = particle_set1(i)%f(2) - fscalar*rij(2)
955 particle_set1(i)%f(3) = particle_set1(i)%f(3) - fscalar*rij(3)
956
957 particle_set2(j)%f(1) = particle_set2(j)%f(1) + fscalar*rij(1)
958 particle_set2(j)%f(2) = particle_set2(j)%f(2) + fscalar*rij(2)
959 particle_set2(j)%f(3) = particle_set2(j)%f(3) + fscalar*rij(3)
960
961 IF (use_virial .AND. shell_adiabatic) THEN
962 fij_com = fscalar*rij
963 fbc(1, 1) = -fij_com(1)*rij(1)
964 fbc(1, 2) = -fij_com(1)*rij(2)
965 fbc(1, 3) = -fij_com(1)*rij(3)
966 fbc(2, 1) = -fij_com(2)*rij(1)
967 fbc(2, 2) = -fij_com(2)*rij(2)
968 fbc(2, 3) = -fij_com(2)*rij(3)
969 fbc(3, 1) = -fij_com(3)*rij(1)
970 fbc(3, 2) = -fij_com(3)*rij(2)
971 fbc(3, 3) = -fij_com(3)*rij(3)
972 pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
973 END IF
974
975 END SUBROUTINE bonded_correct_gaussian_low
976
977! **************************************************************************************************
978!> \brief specific for shell models cleans the interaction core-shell on the same
979!> atom
980!> \param r1 ...
981!> \param r2 ...
982!> \param cell ...
983!> \param v_bonded_corr ...
984!> \param core_particle_set ...
985!> \param shell_particle_set ...
986!> \param i ...
987!> \param shell_adiabatic ...
988!> \param alpha ...
989!> \param q1 ...
990!> \param q2 ...
991!> \param const ...
992!> \param pv_bc ...
993!> \param atprop_env ...
994!> \param use_virial ...
995!> \par History
996!> Split routines to clean and to fix a bug with the tensor whose
997!> original definition was not correct for PBC..
998!> \author Teodoro Laino
999! **************************************************************************************************
1000 SUBROUTINE bonded_correct_gaussian_low_sh(r1, r2, cell, v_bonded_corr, &
1001 core_particle_set, shell_particle_set, i, shell_adiabatic, alpha, q1, q2, &
1002 const, pv_bc, atprop_env, use_virial)
1003 REAL(kind=dp), DIMENSION(3) :: r1, r2
1004 TYPE(cell_type), POINTER :: cell
1005 REAL(kind=dp), INTENT(INOUT) :: v_bonded_corr
1006 TYPE(particle_type), POINTER :: core_particle_set(:), &
1007 shell_particle_set(:)
1008 INTEGER, INTENT(IN) :: i
1009 LOGICAL, INTENT(IN) :: shell_adiabatic
1010 REAL(kind=dp), INTENT(IN) :: alpha, q1, q2, const
1011 REAL(kind=dp), INTENT(INOUT) :: pv_bc(3, 3)
1012 TYPE(atprop_type), POINTER :: atprop_env
1013 LOGICAL, INTENT(IN) :: use_virial
1014
1015 REAL(kind=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
1016 ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, pc = 0.3275911_dp
1017
1018 INTEGER :: iatom
1019 REAL(kind=dp) :: arg, dij, e_arg_arg, efac, errf, ffac, &
1020 fscalar, idij, rijsq, tc, tc2, tc4, vbc
1021 REAL(kind=dp), DIMENSION(3) :: fr, rij
1022 REAL(kind=dp), DIMENSION(3, 3) :: fbc
1023
1024 rij = r1 - r2
1025 rij = pbc(rij, cell)
1026 rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
1027 dij = sqrt(rijsq)
1028 ! Two possible limiting cases according the value of dij
1029 arg = alpha*dij
1030 ! and this is a magic number.. it is related to the order expansion
1031 ! and to the value of the polynomial coefficients
1032 IF (arg > 0.355_dp) THEN
1033 idij = 1.0_dp/dij
1034 e_arg_arg = exp(-arg*arg)
1035 tc = 1.0_dp/(1.0_dp + pc*arg)
1036 ! defining errf = 1 - erfc
1037 errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
1038 efac = idij*errf
1039 ffac = idij**2*(efac - const*e_arg_arg)
1040 ELSE
1041 tc = arg*arg
1042 tc2 = tc*tc
1043 tc4 = tc2*tc2
1044 efac = const*(1.0_dp - tc/3.0_dp + tc2/10.0_dp - tc*tc2/42.0_dp + tc4/216.0_dp - &
1045 tc*tc4/1320.0_dp + tc2*tc4/9360.0_dp)
1046 ffac = const*alpha**2*(2.0_dp/3.0_dp - 2.0_dp*tc/5.0_dp + tc2/7.0_dp - tc*tc2/27.0_dp + &
1047 tc4/132.0_dp - tc*tc4/780.0_dp)
1048 END IF
1049
1050 ! getting the potential
1051 vbc = -q1*q2*efac
1052 v_bonded_corr = v_bonded_corr + vbc
1053 IF (atprop_env%energy) THEN
1054 iatom = shell_particle_set(i)%atom_index
1055 atprop_env%atener(iatom) = atprop_env%atener(iatom) + vbc
1056 END IF
1057
1058 ! subtracting the force from the total force
1059 fscalar = q1*q2*ffac
1060 fr(:) = fscalar*rij(:)
1061
1062 core_particle_set(i)%f(1) = core_particle_set(i)%f(1) - fr(1)
1063 core_particle_set(i)%f(2) = core_particle_set(i)%f(2) - fr(2)
1064 core_particle_set(i)%f(3) = core_particle_set(i)%f(3) - fr(3)
1065
1066 shell_particle_set(i)%f(1) = shell_particle_set(i)%f(1) + fr(1)
1067 shell_particle_set(i)%f(2) = shell_particle_set(i)%f(2) + fr(2)
1068 shell_particle_set(i)%f(3) = shell_particle_set(i)%f(3) + fr(3)
1069
1070 IF (use_virial .AND. shell_adiabatic) THEN
1071 fbc(1, 1) = -fr(1)*rij(1)
1072 fbc(1, 2) = -fr(1)*rij(2)
1073 fbc(1, 3) = -fr(1)*rij(3)
1074 fbc(2, 1) = -fr(2)*rij(1)
1075 fbc(2, 2) = -fr(2)*rij(2)
1076 fbc(2, 3) = -fr(2)*rij(3)
1077 fbc(3, 1) = -fr(3)*rij(1)
1078 fbc(3, 2) = -fr(3)*rij(2)
1079 fbc(3, 3) = -fr(3)*rij(3)
1080 pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
1081 END IF
1082
1083 END SUBROUTINE bonded_correct_gaussian_low_sh
1084
1085END MODULE fist_nonbond_force
cell_types::pbc
Definition cell_types.F:103
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition atomic_kind_types.F:223
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
atprop_types
Holds information on atomic properties.
Definition atprop_types.F:14
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition distribution_1d_types.F:24
ewald_environment_types
Definition ewald_environment_types.F:14
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition ewald_environment_types.F:132
fist_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition fist_neighbor_list_types.F:11
fist_nonbond_env_types
Definition fist_nonbond_env_types.F:13
fist_nonbond_env_types::fist_nonbond_env_get
subroutine, public fist_nonbond_env_get(fist_nonbond_env, potparm14, potparm, nonbonded, rlist_cut, rlist_lowsq, aup, lup, ei_scale14, vdw_scale14, shift_cutoff, do_electrostatics, r_last_update, r_last_update_pbc, rshell_last_update_pbc, rcore_last_update_pbc, cell_last_update, num_update, last_update, counter, natom_types, long_range_correction, ij_kind_full_fac, eam_data, nequip_data, deepmd_data, ace_data, charges)
sets a fist_nonbond_env
Definition fist_nonbond_env_types.F:160
fist_nonbond_force
Definition fist_nonbond_force.F:17
fist_nonbond_force::force_nonbond
subroutine, public force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, pot_nonbond, f_nonbond, pv_nonbond, fshell_nonbond, fcore_nonbond, atprop_env, atomic_kind_set, use_virial)
Calculates the force and the potential of the minimum image, and the pressure tensor.
Definition fist_nonbond_force.F:83
fist_nonbond_force::bonded_correct_gaussian
subroutine, public bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, local_particles, particle_set, ewald_env, v_bonded_corr, pv_bc, shell_particle_set, core_particle_set, atprop_env, cell, use_virial)
corrects electrostatics for bonded terms
Definition fist_nonbond_force.F:642
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
list
An array-based list which grows on demand. When the internal array is full, a new array of twice the ...
Definition list.F:24
machine
Machine interface based on Fortran 2003 and POSIX.
Definition machine.F:17
machine::m_memory
subroutine, public m_memory(mem)
Returns the total amount of memory [bytes] in use, if known, zero otherwise.
Definition machine.F:440
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::oorootpi
real(kind=dp), parameter, public oorootpi
Definition mathconstants.F:115
mathconstants::sqrthalf
real(kind=dp), parameter, public sqrthalf
Definition mathconstants.F:127
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
pair_potential_coulomb
Definition pair_potential_coulomb.F:8
pair_potential_coulomb::potential_coulomb
real(kind=dp) function, public potential_coulomb(r2, fscalar, qfac, ewald_type, alpha, beta, interaction_cutoff)
Evaluates the electrostatic energy and force.
Definition pair_potential_coulomb.F:38
pair_potential_types
Definition pair_potential_types.F:14
pair_potential_types::sh_sh
integer, parameter, public sh_sh
Definition pair_potential_types.F:84
pair_potential_types::nosh_nosh
integer, parameter, public nosh_nosh
Definition pair_potential_types.F:84
pair_potential_types::allegro_type
integer, parameter, public allegro_type
Definition pair_potential_types.F:36
pair_potential_types::gal_type
integer, parameter, public gal_type
Definition pair_potential_types.F:36
pair_potential_types::nequip_type
integer, parameter, public nequip_type
Definition pair_potential_types.F:36
pair_potential_types::deepmd_type
integer, parameter, public deepmd_type
Definition pair_potential_types.F:36
pair_potential_types::siepmann_type
integer, parameter, public siepmann_type
Definition pair_potential_types.F:36
pair_potential_types::nosh_sh
integer, parameter, public nosh_sh
Definition pair_potential_types.F:84
pair_potential_types::ace_type
integer, parameter, public ace_type
Definition pair_potential_types.F:36
pair_potential_types::mace_type
integer, parameter, public mace_type
Definition pair_potential_types.F:36
pair_potential_types::gal21_type
integer, parameter, public gal21_type
Definition pair_potential_types.F:36
pair_potential_types::tersoff_type
integer, parameter, public tersoff_type
Definition pair_potential_types.F:36
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
shell_potential_types
Definition shell_potential_types.F:11
shell_potential_types::get_shell
elemental subroutine, public get_shell(shell, charge, charge_core, charge_shell, mass_core, mass_shell, k2_spring, k4_spring, max_dist, shell_cutoff)
...
Definition shell_potential_types.F:71
splines_methods
routines for handling splines
Definition splines_methods.F:14
splines_methods::potential_s
real(kind=dp) function, public potential_s(spl_p, xxi, y1, spl_f, logger)
calculates the potential interpolated with splines value at a given point and the first derivative....
Definition splines_methods.F:162
splines_types
routines for handling splines_types
Definition splines_types.F:14
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:64
fist_neighbor_list_types::fist_neighbor_type
Definition fist_neighbor_list_types.F:48
fist_neighbor_list_types::neighbor_kind_pairs_type
Definition fist_neighbor_list_types.F:26
fist_nonbond_env_types::fist_nonbond_env_type
Definition fist_nonbond_env_types.F:86
fist_nonbond_env_types::pos_type
Definition fist_nonbond_env_types.F:43
message_passing::mp_comm_type
Definition message_passing.F:158
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:404
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:371
particle_types::particle_type
Definition particle_types.F:35
shell_potential_types::shell_kind_type
Define the shell type.
Definition shell_potential_types.F:28
splines_types::spline_data_p_type
Definition splines_types.F:55
splines_types::spline_factor_type
Definition splines_types.F:71