(git:e5b1968)
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-2025 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, nequip_type, nosh_nosh, &
42 nosh_sh, pair_potential_pp_type, pair_potential_single_type, sh_sh, siepmann_type, &
43 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
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 == ace_type) .OR. any(pot%type == deepmd_type)
235 IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
236 fac_ei = 0.5_dp*fac_ei
237 fac_vdw = 0.5_dp*fac_vdw
238 END IF
239 ! decide which interactions to compute\b
240 IF (do_multipoles .OR. (.NOT. fist_nonbond_env%do_electrostatics)) THEN
241 fac_ei = 0.0_dp
242 END IF
243 IF (ipair <= neighbor_kind_pair%nscale) THEN
244 fac_ei = fac_ei*neighbor_kind_pair%ei_scale(ipair)
245 fac_vdw = fac_vdw*neighbor_kind_pair%vdw_scale(ipair)
246 END IF
247
248 IF (fac_ei > 0.0_dp) THEN
249 ! Get the electrostatic parameters for the atoms a and b
250 mm_radius_a = mm_radius(kind_a)
251 mm_radius_b = mm_radius(kind_b)
252 IF (ASSOCIATED(fist_nonbond_env%charges)) THEN
253 qeff_a = fist_nonbond_env%charges(atom_a)
254 qeff_b = fist_nonbond_env%charges(atom_b)
255 ELSE
256 qeff_a = qeff(kind_a)
257 qeff_b = qeff(kind_b)
258 END IF
259 IF (is_shell_kind(kind_a)) THEN
260 qcore_a = qcore(kind_a)
261 qshell_a = qshell(kind_a)
262 IF ((qcore_a == 0.0_dp) .AND. (qshell_a == 0.0_dp)) fac_ei = 0.0_dp
263 ELSE
264 qcore_a = qeff_a
265 qshell_a = huge(0.0_dp)
266 IF (qeff_a == 0.0_dp) fac_ei = 0.0_dp
267 END IF
268 IF (is_shell_kind(kind_b)) THEN
269 qcore_b = qcore(kind_b)
270 qshell_b = qshell(kind_b)
271 IF ((qcore_b == 0.0_dp) .AND. (qshell_b == 0.0_dp)) fac_ei = 0.0_dp
272 ELSE
273 qcore_b = qeff_b
274 qshell_b = huge(0.0_dp)
275 IF (qeff_b == 0.0_dp) fac_ei = 0.0_dp
276 END IF
277 ! Derive beta parameters
278 beta = 0.0_dp
279 beta_a = 0.0_dp
280 beta_b = 0.0_dp
281 IF (mm_radius_a > 0) THEN
282 beta_a = sqrthalf/mm_radius_a
283 END IF
284 IF (mm_radius_b > 0) THEN
285 beta_b = sqrthalf/mm_radius_b
286 END IF
287 IF ((mm_radius_a > 0) .OR. (mm_radius_b > 0)) THEN
288 beta = sqrthalf/sqrt(mm_radius_a*mm_radius_a + mm_radius_b*mm_radius_b)
289 END IF
290 END IF
291
292 ! In case we have only manybody potentials and no charges, this
293 ! pair of atom types can be ignored here.
294 IF (pot%no_pp .AND. (fac_ei == 0.0)) cycle
295
296 ! Setup spline_data set
297 spl_f => pot%spl_f
298 spline_data => pot%pair_spline_data
299 shell_type = pot%shell_type
300 IF (shell_type /= nosh_nosh) THEN
301 cpassert(.NOT. do_multipoles)
302 cpassert(shell_present)
303 END IF
304 rab2_max = pot%rcutsq
305
306 ! compute the relative vector(s) for this pair
307 IF (shell_type /= nosh_nosh) THEN
308 ! do shell
309 all_terms = .true.
310 IF (shell_type == sh_sh) THEN
311 shell_a = particle_set(atom_a)%shell_index
312 shell_b = particle_set(atom_b)%shell_index
313 rab_cc = rcore_last_update_pbc(shell_b)%r - rcore_last_update_pbc(shell_a)%r
314 rab_cs = rshell_last_update_pbc(shell_b)%r - rcore_last_update_pbc(shell_a)%r
315 rab_sc = rcore_last_update_pbc(shell_b)%r - rshell_last_update_pbc(shell_a)%r
316 rab_ss = rshell_last_update_pbc(shell_b)%r - rshell_last_update_pbc(shell_a)%r
317 rab_list(1:3, 1) = rab_cc(1:3) + cell_v(1:3)
318 rab_list(1:3, 2) = rab_cs(1:3) + cell_v(1:3)
319 rab_list(1:3, 3) = rab_sc(1:3) + cell_v(1:3)
320 rab_list(1:3, 4) = rab_ss(1:3) + cell_v(1:3)
321 ELSE IF ((shell_type == nosh_sh) .AND. (particle_set(atom_a)%shell_index /= 0)) THEN
322 shell_a = particle_set(atom_a)%shell_index
323 shell_b = 0
324 rab_cc = r_last_update_pbc(atom_b)%r - rcore_last_update_pbc(shell_a)%r
325 rab_sc = 0.0_dp
326 rab_cs = 0.0_dp
327 rab_ss = r_last_update_pbc(atom_b)%r - rshell_last_update_pbc(shell_a)%r
328 rab_list(1:3, 1) = rab_cc(1:3) + cell_v(1:3)
329 rab_list(1:3, 2) = 0.0_dp
330 rab_list(1:3, 3) = 0.0_dp
331 rab_list(1:3, 4) = rab_ss(1:3) + cell_v(1:3)
332 ELSE IF ((shell_type == nosh_sh) .AND. (particle_set(atom_b)%shell_index /= 0)) THEN
333 shell_b = particle_set(atom_b)%shell_index
334 shell_a = 0
335 rab_cc = rcore_last_update_pbc(shell_b)%r - r_last_update_pbc(atom_a)%r
336 rab_sc = 0.0_dp
337 rab_cs = 0.0_dp
338 rab_ss = rshell_last_update_pbc(shell_b)%r - r_last_update_pbc(atom_a)%r
339 rab_list(1:3, 1) = rab_cc(1:3) + cell_v(1:3)
340 rab_list(1:3, 2) = 0.0_dp
341 rab_list(1:3, 3) = 0.0_dp
342 rab_list(1:3, 4) = rab_ss(1:3) + cell_v(1:3)
343 ELSE
344 rab_list(:, :) = 0.0_dp
345 END IF
346 ! Compute the term only if all the pairs (cc,cs,sc,ss) are within the cut-off
347 check_terms: DO i = 1, 4
348 rab2_list(i) = rab_list(1, i)**2 + rab_list(2, i)**2 + rab_list(3, i)**2
349 IF (rab2_list(i) >= rab2_max) THEN
350 all_terms = .false.
351 EXIT check_terms
352 END IF
353 END DO check_terms
354 rab_com = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
355 ELSE
356 ! not do shell
357 rab_cc = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
358 rab_com = rab_cc
359 shell_a = 0
360 shell_b = 0
361 rab_list(:, :) = 0.0_dp
362 END IF
363 rab_com = rab_com + cell_v
364 rab2_com = rab_com(1)**2 + rab_com(2)**2 + rab_com(3)**2
365
366 ! compute the interactions for the current pair
367 etot = 0.0_dp
368 fatom_a(:) = 0.0_dp
369 fatom_b(:) = 0.0_dp
370 fcore_a(:) = 0.0_dp
371 fcore_b(:) = 0.0_dp
372 fshell_a(:) = 0.0_dp
373 fshell_b(:) = 0.0_dp
374 IF (use_virial) pv(:, :) = 0.0_dp
375 IF (shell_type /= nosh_nosh) THEN
376 ! do shell
377 IF ((rab2_com <= rab2_max) .AND. all_terms) THEN
378 IF (fac_ei > 0) THEN
379 ! core-core or core-ion/ion-core: Coulomb only
380 rab = rab_list(:, 1)
381 rab2 = rab2_list(1)
382 fscalar = 0.0_dp
383 IF (shell_a == 0) THEN
384 ! atom a is a plain ion and can have beta_a > 0
385 energy = potential_coulomb(rab2, fscalar, fac_ei*qeff_a*qcore_b, &
386 ewald_type, alpha, beta_a, &
387 ei_interaction_cutoffs(2, kind_a, kind_b))
388 CALL add_force_nonbond(fatom_a, fcore_b, pv, fscalar, rab, use_virial)
389 ELSE IF (shell_b == 0) THEN
390 ! atom b is a plain ion and can have beta_b > 0
391 energy = potential_coulomb(rab2, fscalar, fac_ei*qcore_a*qeff_b, &
392 ewald_type, alpha, beta_b, &
393 ei_interaction_cutoffs(2, kind_b, kind_a))
394 CALL add_force_nonbond(fcore_a, fatom_b, pv, fscalar, rab, use_virial)
395 ELSE
396 ! core-core interaction is always pure point charge
397 energy = potential_coulomb(rab2, fscalar, fac_ei*qcore_a*qcore_b, &
398 ewald_type, alpha, 0.0_dp, &
399 ei_interaction_cutoffs(1, kind_a, kind_b))
400 CALL add_force_nonbond(fcore_a, fcore_b, pv, fscalar, rab, use_virial)
401 END IF
402 etot = etot + energy
403 END IF
404
405 IF (shell_type == sh_sh) THEN
406 ! shell-shell: VDW + Coulomb
407 rab = rab_list(:, 4)
408 rab2 = rab2_list(4)
409 fscalar = 0.0_dp
410 IF (fac_vdw > 0) THEN
411 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
412 etot = etot + energy*fac_vdw
413 fscalar = fscalar*fac_vdw
414 END IF
415 IF (fac_ei > 0) THEN
416 ! note that potential_coulomb increments fscalar
417 energy = potential_coulomb(rab2, fscalar, fac_ei*qshell_a*qshell_b, &
418 ewald_type, alpha, beta, &
419 ei_interaction_cutoffs(3, kind_a, kind_b))
420 etot = etot + energy
421 END IF
422 CALL add_force_nonbond(fshell_a, fshell_b, pv, fscalar, rab, use_virial)
423
424 IF (fac_ei > 0) THEN
425 ! core-shell: Coulomb only
426 rab = rab_list(:, 2)
427 rab2 = rab2_list(2)
428 fscalar = 0.0_dp
429 ! swap kind_a and kind_b to get the right cutoff
430 energy = potential_coulomb(rab2, fscalar, fac_ei*qcore_a*qshell_b, &
431 ewald_type, alpha, beta_b, &
432 ei_interaction_cutoffs(2, kind_b, kind_a))
433 etot = etot + energy
434 CALL add_force_nonbond(fcore_a, fshell_b, pv, fscalar, rab, use_virial)
435
436 ! shell-core: Coulomb only
437 rab = rab_list(:, 3)
438 rab2 = rab2_list(3)
439 fscalar = 0.0_dp
440 energy = potential_coulomb(rab2, fscalar, fac_ei*qshell_a*qcore_b, &
441 ewald_type, alpha, beta_a, &
442 ei_interaction_cutoffs(2, kind_a, kind_b))
443 etot = etot + energy
444 CALL add_force_nonbond(fshell_a, fcore_b, pv, fscalar, rab, use_virial)
445 END IF
446 ELSE IF ((shell_type == nosh_sh) .AND. (shell_a == 0)) THEN
447 ! ion-shell: VDW + Coulomb
448 rab = rab_list(:, 4)
449 rab2 = rab2_list(4)
450 fscalar = 0.0_dp
451 IF (fac_vdw > 0) THEN
452 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
453 etot = etot + energy*fac_vdw
454 fscalar = fscalar*fac_vdw
455 END IF
456 IF (fac_ei > 0) THEN
457 ! note that potential_coulomb increments fscalar
458 energy = potential_coulomb(rab2, fscalar, fac_ei*qeff_a*qshell_b, &
459 ewald_type, alpha, beta, &
460 ei_interaction_cutoffs(3, kind_a, kind_b))
461 etot = etot + energy
462 END IF
463 CALL add_force_nonbond(fatom_a, fshell_b, pv, fscalar, rab, use_virial)
464 ELSE IF ((shell_type == nosh_sh) .AND. (shell_b == 0)) THEN
465 ! shell-ion : VDW + Coulomb
466 rab = rab_list(:, 4)
467 rab2 = rab2_list(4)
468 fscalar = 0.0_dp
469 IF (fac_vdw > 0) THEN
470 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
471 etot = etot + energy*fac_vdw
472 fscalar = fscalar*fac_vdw
473 END IF
474 IF (fac_ei > 0) THEN
475 ! note that potential_coulomb increments fscalar
476 energy = potential_coulomb(rab2, fscalar, fac_ei*qshell_a*qeff_b, &
477 ewald_type, alpha, beta, &
478 ei_interaction_cutoffs(3, kind_a, kind_b))
479 etot = etot + energy
480 END IF
481 CALL add_force_nonbond(fshell_a, fatom_b, pv, fscalar, rab, use_virial)
482 END IF
483 END IF
484 ELSE
485 IF (rab2_com <= rab2_max) THEN
486 ! NO SHELL MODEL...
487 ! Ion-Ion: no shell model, VDW + coulomb
488 rab = rab_com
489 rab2 = rab2_com
490 fscalar = 0.0_dp
491 IF (fac_vdw > 0) THEN
492 energy = potential_s(spline_data, rab2, fscalar, spl_f, logger)
493 etot = etot + energy*fac_vdw
494 fscalar = fscalar*fac_vdw
495 END IF
496 IF (fac_ei > 0) THEN
497 ! note that potential_coulomb increments fscalar
498 energy = potential_coulomb(rab2, fscalar, fac_ei*qeff_a*qeff_b, &
499 ewald_type, alpha, beta, &
500 ei_interaction_cutoffs(3, kind_a, kind_b))
501 etot = etot + energy
502 END IF
503 CALL add_force_nonbond(fatom_a, fatom_b, pv, fscalar, rab, use_virial)
504 END IF
505 END IF
506 ! Nonbonded energy
507!$OMP ATOMIC
508 pot_nonbond = pot_nonbond + etot
509 IF (atprop_env%energy) THEN
510 ! Update atomic energies
511!$OMP ATOMIC
512 atprop_env%atener(atom_a) = atprop_env%atener(atom_a) + 0.5_dp*etot
513!$OMP ATOMIC
514 atprop_env%atener(atom_b) = atprop_env%atener(atom_b) + 0.5_dp*etot
515 END IF
516 ! Nonbonded forces
517 DO i = 1, 3
518!$OMP ATOMIC
519 f_nonbond(i, atom_a) = f_nonbond(i, atom_a) + fatom_a(i)
520!$OMP ATOMIC
521 f_nonbond(i, atom_b) = f_nonbond(i, atom_b) + fatom_b(i)
522 END DO
523 IF (shell_a > 0) THEN
524 DO i = 1, 3
525!$OMP ATOMIC
526 fcore_nonbond(i, shell_a) = fcore_nonbond(i, shell_a) + fcore_a(i)
527!$OMP ATOMIC
528 fshell_nonbond(i, shell_a) = fshell_nonbond(i, shell_a) + fshell_a(i)
529 END DO
530 END IF
531 IF (shell_b > 0) THEN
532 DO i = 1, 3
533!$OMP ATOMIC
534 fcore_nonbond(i, shell_b) = fcore_nonbond(i, shell_b) + fcore_b(i)
535!$OMP ATOMIC
536 fshell_nonbond(i, shell_b) = fshell_nonbond(i, shell_b) + fshell_b(i)
537 END DO
538 END IF
539 ! Add the contribution of the current pair to the total pressure tensor
540 IF (use_virial) THEN
541 DO i = 1, 3
542 DO j = 1, 3
543 pv_thread(j, i) = pv_thread(j, i) + pv(j, i)
544 END DO
545 END DO
546 END IF
547 END DO pairs
548!$OMP END DO
549 IF (use_virial) THEN
550 DO i = 1, 3
551 DO j = 1, 3
552!$OMP ATOMIC
553 pv_nonbond(j, i) = pv_nonbond(j, i) + pv_thread(j, i)
554 END DO
555 END DO
556 END IF
557!$OMP END PARALLEL
558 END DO kind_group_loop
559 END DO lists
560
561 !sample peak memory
562 CALL m_memory()
563
564 DEALLOCATE (mm_radius)
565 DEALLOCATE (qeff)
566 DEALLOCATE (qcore)
567 DEALLOCATE (qshell)
568 DEALLOCATE (is_shell_kind)
569
570 CALL timestop(handle)
571
572 END SUBROUTINE force_nonbond
573
574 ! **************************************************************************************************
575 !> \brief Adds a non-bonding contribution to the total force and optionally to
576 !> the virial.
577 ! **************************************************************************************************
578! **************************************************************************************************
579!> \brief ...
580!> \param f_nonbond_a ...
581!> \param f_nonbond_b ...
582!> \param pv ...
583!> \param fscalar ...
584!> \param rab ...
585!> \param use_virial ...
586! **************************************************************************************************
587 SUBROUTINE add_force_nonbond(f_nonbond_a, f_nonbond_b, pv, fscalar, rab, use_virial)
588
589 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: f_nonbond_a, f_nonbond_b
590 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv
591 REAL(kind=dp), INTENT(IN) :: fscalar
592 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
593 LOGICAL, INTENT(IN) :: use_virial
594
595 REAL(kind=dp), DIMENSION(3) :: fr
596
597 fr(1) = fscalar*rab(1)
598 fr(2) = fscalar*rab(2)
599 fr(3) = fscalar*rab(3)
600 f_nonbond_a(1) = f_nonbond_a(1) - fr(1)
601 f_nonbond_a(2) = f_nonbond_a(2) - fr(2)
602 f_nonbond_a(3) = f_nonbond_a(3) - fr(3)
603 f_nonbond_b(1) = f_nonbond_b(1) + fr(1)
604 f_nonbond_b(2) = f_nonbond_b(2) + fr(2)
605 f_nonbond_b(3) = f_nonbond_b(3) + fr(3)
606 IF (use_virial) THEN
607 pv(1, 1) = pv(1, 1) + rab(1)*fr(1)
608 pv(1, 2) = pv(1, 2) + rab(1)*fr(2)
609 pv(1, 3) = pv(1, 3) + rab(1)*fr(3)
610 pv(2, 1) = pv(2, 1) + rab(2)*fr(1)
611 pv(2, 2) = pv(2, 2) + rab(2)*fr(2)
612 pv(2, 3) = pv(2, 3) + rab(2)*fr(3)
613 pv(3, 1) = pv(3, 1) + rab(3)*fr(1)
614 pv(3, 2) = pv(3, 2) + rab(3)*fr(2)
615 pv(3, 3) = pv(3, 3) + rab(3)*fr(3)
616 END IF
617
618 END SUBROUTINE
619
620! **************************************************************************************************
621!> \brief corrects electrostatics for bonded terms
622!> \param fist_nonbond_env ...
623!> \param atomic_kind_set ...
624!> \param local_particles ...
625!> \param particle_set ...
626!> \param ewald_env ...
627!> \param v_bonded_corr ...
628!> \param pv_bc ...
629!> \param shell_particle_set ...
630!> \param core_particle_set ...
631!> \param atprop_env ...
632!> \param cell ...
633!> \param use_virial ...
634!> \par History
635!> Split routines to clean and to fix a bug with the tensor whose
636!> original definition was not correct for PBC.. [Teodoro Laino -06/2007]
637! **************************************************************************************************
638 SUBROUTINE bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
639 local_particles, particle_set, ewald_env, v_bonded_corr, pv_bc, &
640 shell_particle_set, core_particle_set, atprop_env, cell, use_virial)
641
642 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
643 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
644 TYPE(distribution_1d_type), POINTER :: local_particles
645 TYPE(particle_type), POINTER :: particle_set(:)
646 TYPE(ewald_environment_type), POINTER :: ewald_env
647 REAL(kind=dp), INTENT(OUT) :: v_bonded_corr
648 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: pv_bc
649 TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
650 core_particle_set(:)
651 TYPE(atprop_type), POINTER :: atprop_env
652 TYPE(cell_type), POINTER :: cell
653 LOGICAL, INTENT(IN) :: use_virial
654
655 CHARACTER(LEN=*), PARAMETER :: routinen = 'bonded_correct_gaussian'
656
657 INTEGER :: atom_a, atom_b, handle, iatom, iend, igrp, ilist, ipair, istart, kind_a, kind_b, &
658 natoms_per_kind, nkind, npairs, shell_a, shell_b
659 INTEGER, DIMENSION(:, :), POINTER :: list
660 LOGICAL :: a_is_shell, b_is_shell, do_multipoles, &
661 full_nl, shell_adiabatic
662 REAL(kind=dp) :: alpha, const, fac_cor, fac_ei, qcore_a, &
663 qcore_b, qeff_a, qeff_b, qshell_a, &
664 qshell_b
665 REAL(kind=dp), DIMENSION(3) :: rca, rcb, rsa, rsb
666 REAL(kind=dp), DIMENSION(:, :), POINTER :: ij_kind_full_fac
667 TYPE(atomic_kind_type), POINTER :: atomic_kind
668 TYPE(fist_neighbor_type), POINTER :: nonbonded
669 TYPE(mp_comm_type) :: group
670 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
671 TYPE(pair_potential_pp_type), POINTER :: potparm, potparm14
672 TYPE(pair_potential_single_type), POINTER :: pot
673 TYPE(shell_kind_type), POINTER :: shell_kind
674
675 CALL timeset(routinen, handle)
676
677 ! Initializing values
678 IF (use_virial) pv_bc = 0.0_dp
679 v_bonded_corr = 0.0_dp
680
681 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, &
682 potparm14=potparm14, potparm=potparm, &
683 ij_kind_full_fac=ij_kind_full_fac)
684 CALL ewald_env_get(ewald_env, alpha=alpha, do_multipoles=do_multipoles, &
685 group=group)
686 ! Defining the constants
687 const = 2.0_dp*alpha*oorootpi
688
689 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
690 shell_adiabatic=shell_adiabatic)
691
692 lists: DO ilist = 1, nonbonded%nlists
693 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
694 npairs = neighbor_kind_pair%nscale
695 IF (npairs == 0) cycle
696 list => neighbor_kind_pair%list
697 kind_group_loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
698 istart = neighbor_kind_pair%grp_kind_start(igrp)
699 IF (istart > npairs) THEN
700 EXIT
701 END IF
702 iend = min(npairs, neighbor_kind_pair%grp_kind_end(igrp))
703
704 pairs: DO ipair = istart, iend
705 atom_a = list(1, ipair)
706 atom_b = list(2, ipair)
707 ! Get actual atomic kinds, since atom_a is not always of
708 ! kind_a and atom_b of kind_b, ie. they might be swapped.
709 kind_a = particle_set(atom_a)%atomic_kind%kind_number
710 kind_b = particle_set(atom_b)%atomic_kind%kind_number
711
712 ! take the proper potential, only for full_nl test
713 pot => potparm%pot(kind_a, kind_b)%pot
714 IF (ipair <= neighbor_kind_pair%nscale) THEN
715 IF (neighbor_kind_pair%is_onfo(ipair)) THEN
716 pot => potparm14%pot(kind_a, kind_b)%pot
717 END IF
718 END IF
719
720 ! Determine the scaling factors
721 fac_ei = ij_kind_full_fac(kind_a, kind_b)
722 full_nl = any(pot%type == tersoff_type) .OR. any(pot%type == siepmann_type) &
723 .OR. any(pot%type == gal_type) .OR. any(pot%type == gal21_type) &
724 .OR. any(pot%type == nequip_type) .OR. any(pot%type == allegro_type) &
725 .OR. any(pot%type == ace_type) .OR. any(pot%type == deepmd_type)
726 IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
727 fac_ei = fac_ei*0.5_dp
728 END IF
729 IF (ipair <= neighbor_kind_pair%nscale) THEN
730 fac_ei = fac_ei*neighbor_kind_pair%ei_scale(ipair)
731 END IF
732 ! The amount of correction is related to the
733 ! amount of scaling as follows:
734 fac_cor = 1.0_dp - fac_ei
735 IF (fac_cor <= 0.0_dp) cycle
736
737 ! Parameters for kind a
738 atomic_kind => atomic_kind_set(kind_a)
739 CALL get_atomic_kind(atomic_kind, qeff=qeff_a, shell=shell_kind)
740 IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_a = fist_nonbond_env%charges(atom_a)
741 a_is_shell = ASSOCIATED(shell_kind)
742 IF (a_is_shell) THEN
743 CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
744 charge_shell=qshell_a)
745 shell_a = particle_set(atom_a)%shell_index
746 rca = core_particle_set(shell_a)%r
747 rsa = shell_particle_set(shell_a)%r
748 ELSE
749 qcore_a = qeff_a
750 qshell_a = huge(0.0_dp)
751 shell_a = 0
752 rca = particle_set(atom_a)%r
753 rsa = 0.0_dp
754 END IF
755
756 ! Parameters for kind b
757 atomic_kind => atomic_kind_set(kind_b)
758 CALL get_atomic_kind(atomic_kind, qeff=qeff_b, shell=shell_kind)
759 IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_b = fist_nonbond_env%charges(atom_b)
760 b_is_shell = ASSOCIATED(shell_kind)
761 IF (b_is_shell) THEN
762 CALL get_shell(shell=shell_kind, charge_core=qcore_b, &
763 charge_shell=qshell_b)
764 shell_b = particle_set(atom_b)%shell_index
765 rcb = core_particle_set(shell_b)%r
766 rsb = shell_particle_set(shell_b)%r
767 ELSE
768 qcore_b = qeff_b
769 qshell_b = huge(0.0_dp)
770 shell_b = 0
771 rcb = particle_set(atom_b)%r
772 rsb = 0.0_dp
773 END IF
774
775 ! First part: take care of core/ion-core/ion correction
776 IF (a_is_shell .AND. b_is_shell) THEN
777 ! correct for core-core interaction
778 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
779 v_bonded_corr, core_particle_set, core_particle_set, &
780 shell_a, shell_b, .true., alpha, qcore_a, qcore_b, &
781 const, fac_cor, pv_bc, atprop_env, use_virial)
782 ELSE IF (a_is_shell) THEN
783 ! correct for core-ion interaction
784 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
785 v_bonded_corr, core_particle_set, particle_set, &
786 shell_a, atom_b, .true., alpha, qcore_a, qcore_b, &
787 const, fac_cor, pv_bc, atprop_env, use_virial)
788 ELSE IF (b_is_shell) THEN
789 ! correct for ion-core interaction
790 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
791 v_bonded_corr, particle_set, core_particle_set, &
792 atom_a, shell_b, .true., alpha, qcore_a, qcore_b, &
793 const, fac_cor, pv_bc, atprop_env, use_virial)
794 ELSE
795 ! correct for ion-ion interaction
796 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
797 v_bonded_corr, particle_set, particle_set, &
798 atom_a, atom_b, .true., alpha, qcore_a, qcore_b, &
799 const, fac_cor, pv_bc, atprop_env, use_virial)
800 END IF
801
802 ! Second part: take care of shell-shell, shell-core/ion and
803 ! core/ion-shell corrections
804 IF (a_is_shell .AND. b_is_shell) THEN
805 ! correct for shell-shell interaction
806 CALL bonded_correct_gaussian_low(rsa, rsa, cell, &
807 v_bonded_corr, shell_particle_set, shell_particle_set, &
808 shell_a, shell_b, shell_adiabatic, alpha, qshell_a, &
809 qshell_b, const, fac_cor, pv_bc, atprop_env, use_virial)
810 END IF
811 IF (a_is_shell) THEN
812 IF (b_is_shell) THEN
813 ! correct for shell-core interaction
814 CALL bonded_correct_gaussian_low(rsa, rcb, cell, &
815 v_bonded_corr, shell_particle_set, core_particle_set, &
816 shell_a, shell_b, shell_adiabatic, alpha, qshell_a, qcore_b, &
817 const, fac_cor, pv_bc, atprop_env, use_virial)
818 ELSE
819 ! correct for shell-ion interaction
820 CALL bonded_correct_gaussian_low(rsa, rcb, cell, &
821 v_bonded_corr, shell_particle_set, particle_set, &
822 shell_a, atom_b, shell_adiabatic, alpha, qshell_a, qcore_b, &
823 const, fac_cor, pv_bc, atprop_env, use_virial)
824 END IF
825 END IF
826 IF (b_is_shell) THEN
827 IF (a_is_shell) THEN
828 ! correct for core-shell interaction
829 CALL bonded_correct_gaussian_low(rca, rsb, cell, &
830 v_bonded_corr, core_particle_set, shell_particle_set, &
831 shell_a, shell_b, shell_adiabatic, alpha, qcore_a, qshell_b, &
832 const, fac_cor, pv_bc, atprop_env, use_virial)
833 ELSE
834 ! correct for ion-shell interaction
835 CALL bonded_correct_gaussian_low(rca, rsb, cell, &
836 v_bonded_corr, particle_set, shell_particle_set, &
837 atom_a, shell_b, shell_adiabatic, alpha, qcore_a, qshell_b, &
838 const, fac_cor, pv_bc, atprop_env, use_virial)
839 END IF
840 END IF
841 END DO pairs
842 END DO kind_group_loop
843 END DO lists
844
845 ! Always correct core-shell interaction within one atom.
846 nkind = SIZE(atomic_kind_set)
847 DO kind_a = 1, nkind
848 ! parameters for kind a
849 atomic_kind => atomic_kind_set(kind_a)
850 CALL get_atomic_kind(atomic_kind, shell=shell_kind)
851 IF (ASSOCIATED(shell_kind)) THEN
852 CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
853 charge_shell=qshell_a)
854
855 natoms_per_kind = local_particles%n_el(kind_a)
856 DO iatom = 1, natoms_per_kind
857
858 ! Data for atom a
859 atom_a = local_particles%list(kind_a)%array(iatom)
860 shell_a = particle_set(atom_a)%shell_index
861 rca = core_particle_set(shell_a)%r
862 rsa = shell_particle_set(shell_a)%r
863
864 CALL bonded_correct_gaussian_low_sh(rca, rsa, cell, &
865 v_bonded_corr, core_particle_set, shell_particle_set, &
866 shell_a, shell_adiabatic, alpha, qcore_a, qshell_a, &
867 const, pv_bc, atprop_env, use_virial)
868
869 END DO
870 END IF
871 END DO
872
873 CALL group%sum(v_bonded_corr)
874
875 CALL timestop(handle)
876
877 END SUBROUTINE bonded_correct_gaussian
878
879! **************************************************************************************************
880!> \brief ...
881!> \param r1 ...
882!> \param r2 ...
883!> \param cell ...
884!> \param v_bonded_corr ...
885!> \param particle_set1 ...
886!> \param particle_set2 ...
887!> \param i ...
888!> \param j ...
889!> \param shell_adiabatic ...
890!> \param alpha ...
891!> \param q1 ...
892!> \param q2 ...
893!> \param const ...
894!> \param fac_cor ...
895!> \param pv_bc ...
896!> \param atprop_env ...
897!> \param use_virial ...
898!> \par History
899!> Split routines to clean and to fix a bug with the tensor whose
900!> original definition was not correct for PBC..
901!> \author Teodoro Laino
902! **************************************************************************************************
903 SUBROUTINE bonded_correct_gaussian_low(r1, r2, cell, v_bonded_corr, &
904 particle_set1, particle_set2, i, j, shell_adiabatic, alpha, q1, q2, &
905 const, fac_cor, pv_bc, atprop_env, use_virial)
906 REAL(kind=dp), DIMENSION(3) :: r1, r2
907 TYPE(cell_type), POINTER :: cell
908 REAL(kind=dp), INTENT(INOUT) :: v_bonded_corr
909 TYPE(particle_type), POINTER :: particle_set1(:), particle_set2(:)
910 INTEGER, INTENT(IN) :: i, j
911 LOGICAL, INTENT(IN) :: shell_adiabatic
912 REAL(kind=dp), INTENT(IN) :: alpha, q1, q2, const, fac_cor
913 REAL(kind=dp), INTENT(INOUT) :: pv_bc(3, 3)
914 TYPE(atprop_type), POINTER :: atprop_env
915 LOGICAL, INTENT(IN) :: use_virial
916
917 REAL(kind=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
918 ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, pc = 0.3275911_dp
919
920 INTEGER :: iatom, jatom
921 REAL(kind=dp) :: arg, dij, e_arg_arg, errf, fscalar, &
922 idij, rijsq, tc, vbc
923 REAL(kind=dp), DIMENSION(3) :: fij_com, rij
924 REAL(kind=dp), DIMENSION(3, 3) :: fbc
925
926 rij = r1 - r2
927 rij = pbc(rij, cell)
928 rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
929 idij = 1.0_dp/sqrt(rijsq)
930 dij = rijsq*idij
931 arg = alpha*dij
932 e_arg_arg = exp(-arg**2)
933 tc = 1.0_dp/(1.0_dp + pc*arg)
934
935 ! Defining errf=1-erfc
936 errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
937
938 ! Getting the potential
939 vbc = -q1*q2*idij*errf*fac_cor
940 v_bonded_corr = v_bonded_corr + vbc
941 IF (atprop_env%energy) THEN
942 iatom = particle_set1(i)%atom_index
943 atprop_env%atener(iatom) = atprop_env%atener(iatom) + 0.5_dp*vbc
944 jatom = particle_set2(j)%atom_index
945 atprop_env%atener(jatom) = atprop_env%atener(jatom) + 0.5_dp*vbc
946 END IF
947
948 ! Subtracting the force from the total force
949 fscalar = q1*q2*idij**2*(idij*errf - const*e_arg_arg)*fac_cor
950
951 particle_set1(i)%f(1) = particle_set1(i)%f(1) - fscalar*rij(1)
952 particle_set1(i)%f(2) = particle_set1(i)%f(2) - fscalar*rij(2)
953 particle_set1(i)%f(3) = particle_set1(i)%f(3) - fscalar*rij(3)
954
955 particle_set2(j)%f(1) = particle_set2(j)%f(1) + fscalar*rij(1)
956 particle_set2(j)%f(2) = particle_set2(j)%f(2) + fscalar*rij(2)
957 particle_set2(j)%f(3) = particle_set2(j)%f(3) + fscalar*rij(3)
958
959 IF (use_virial .AND. shell_adiabatic) THEN
960 fij_com = fscalar*rij
961 fbc(1, 1) = -fij_com(1)*rij(1)
962 fbc(1, 2) = -fij_com(1)*rij(2)
963 fbc(1, 3) = -fij_com(1)*rij(3)
964 fbc(2, 1) = -fij_com(2)*rij(1)
965 fbc(2, 2) = -fij_com(2)*rij(2)
966 fbc(2, 3) = -fij_com(2)*rij(3)
967 fbc(3, 1) = -fij_com(3)*rij(1)
968 fbc(3, 2) = -fij_com(3)*rij(2)
969 fbc(3, 3) = -fij_com(3)*rij(3)
970 pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
971 END IF
972
973 END SUBROUTINE bonded_correct_gaussian_low
974
975! **************************************************************************************************
976!> \brief specific for shell models cleans the interaction core-shell on the same
977!> atom
978!> \param r1 ...
979!> \param r2 ...
980!> \param cell ...
981!> \param v_bonded_corr ...
982!> \param core_particle_set ...
983!> \param shell_particle_set ...
984!> \param i ...
985!> \param shell_adiabatic ...
986!> \param alpha ...
987!> \param q1 ...
988!> \param q2 ...
989!> \param const ...
990!> \param pv_bc ...
991!> \param atprop_env ...
992!> \param use_virial ...
993!> \par History
994!> Split routines to clean and to fix a bug with the tensor whose
995!> original definition was not correct for PBC..
996!> \author Teodoro Laino
997! **************************************************************************************************
998 SUBROUTINE bonded_correct_gaussian_low_sh(r1, r2, cell, v_bonded_corr, &
999 core_particle_set, shell_particle_set, i, shell_adiabatic, alpha, q1, q2, &
1000 const, pv_bc, atprop_env, use_virial)
1001 REAL(kind=dp), DIMENSION(3) :: r1, r2
1002 TYPE(cell_type), POINTER :: cell
1003 REAL(kind=dp), INTENT(INOUT) :: v_bonded_corr
1004 TYPE(particle_type), POINTER :: core_particle_set(:), &
1005 shell_particle_set(:)
1006 INTEGER, INTENT(IN) :: i
1007 LOGICAL, INTENT(IN) :: shell_adiabatic
1008 REAL(kind=dp), INTENT(IN) :: alpha, q1, q2, const
1009 REAL(kind=dp), INTENT(INOUT) :: pv_bc(3, 3)
1010 TYPE(atprop_type), POINTER :: atprop_env
1011 LOGICAL, INTENT(IN) :: use_virial
1012
1013 REAL(kind=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
1014 ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, pc = 0.3275911_dp
1015
1016 INTEGER :: iatom
1017 REAL(kind=dp) :: arg, dij, e_arg_arg, efac, errf, ffac, &
1018 fscalar, idij, rijsq, tc, tc2, tc4, vbc
1019 REAL(kind=dp), DIMENSION(3) :: fr, rij
1020 REAL(kind=dp), DIMENSION(3, 3) :: fbc
1021
1022 rij = r1 - r2
1023 rij = pbc(rij, cell)
1024 rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
1025 dij = sqrt(rijsq)
1026 ! Two possible limiting cases according the value of dij
1027 arg = alpha*dij
1028 ! and this is a magic number.. it is related to the order expansion
1029 ! and to the value of the polynomial coefficients
1030 IF (arg > 0.355_dp) THEN
1031 idij = 1.0_dp/dij
1032 e_arg_arg = exp(-arg*arg)
1033 tc = 1.0_dp/(1.0_dp + pc*arg)
1034 ! defining errf = 1 - erfc
1035 errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
1036 efac = idij*errf
1037 ffac = idij**2*(efac - const*e_arg_arg)
1038 ELSE
1039 tc = arg*arg
1040 tc2 = tc*tc
1041 tc4 = tc2*tc2
1042 efac = const*(1.0_dp - tc/3.0_dp + tc2/10.0_dp - tc*tc2/42.0_dp + tc4/216.0_dp - &
1043 tc*tc4/1320.0_dp + tc2*tc4/9360.0_dp)
1044 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 + &
1045 tc4/132.0_dp - tc*tc4/780.0_dp)
1046 END IF
1047
1048 ! getting the potential
1049 vbc = -q1*q2*efac
1050 v_bonded_corr = v_bonded_corr + vbc
1051 IF (atprop_env%energy) THEN
1052 iatom = shell_particle_set(i)%atom_index
1053 atprop_env%atener(iatom) = atprop_env%atener(iatom) + vbc
1054 END IF
1055
1056 ! subtracting the force from the total force
1057 fscalar = q1*q2*ffac
1058 fr(:) = fscalar*rij(:)
1059
1060 core_particle_set(i)%f(1) = core_particle_set(i)%f(1) - fr(1)
1061 core_particle_set(i)%f(2) = core_particle_set(i)%f(2) - fr(2)
1062 core_particle_set(i)%f(3) = core_particle_set(i)%f(3) - fr(3)
1063
1064 shell_particle_set(i)%f(1) = shell_particle_set(i)%f(1) + fr(1)
1065 shell_particle_set(i)%f(2) = shell_particle_set(i)%f(2) + fr(2)
1066 shell_particle_set(i)%f(3) = shell_particle_set(i)%f(3) + fr(3)
1067
1068 IF (use_virial .AND. shell_adiabatic) THEN
1069 fbc(1, 1) = -fr(1)*rij(1)
1070 fbc(1, 2) = -fr(1)*rij(2)
1071 fbc(1, 3) = -fr(1)*rij(3)
1072 fbc(2, 1) = -fr(2)*rij(1)
1073 fbc(2, 2) = -fr(2)*rij(2)
1074 fbc(2, 3) = -fr(2)*rij(3)
1075 fbc(3, 1) = -fr(3)*rij(1)
1076 fbc(3, 2) = -fr(3)*rij(2)
1077 fbc(3, 3) = -fr(3)*rij(3)
1078 pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
1079 END IF
1080
1081 END SUBROUTINE bonded_correct_gaussian_low_sh
1082
1083END MODULE fist_nonbond_force
cell_types::pbc
Definition cell_types.F:92
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:124
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, quip_data, nequip_data, allegro_data, deepmd_data, ace_data, charges)
sets a fist_nonbond_env
Definition fist_nonbond_env_types.F:182
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:641
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::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:55
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:56
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:104
fist_nonbond_env_types::pos_type
Definition fist_nonbond_env_types.F:48
message_passing::mp_comm_type
Definition message_passing.F:189
pair_potential_types::pair_potential_pp_type
Definition pair_potential_types.F:428
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:395
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