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