(git:ed6f26b)
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 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 END DO
544 END DO
545 END IF
546 END DO pairs
547!$OMP END DO
548 IF (use_virial) THEN
549 DO i = 1, 3
550 DO j = 1, 3
551!$OMP ATOMIC
552 pv_nonbond(j, i) = pv_nonbond(j, i) + pv_thread(j, i)
553 END DO
554 END DO
555 END IF
556!$OMP END PARALLEL
557 END DO kind_group_loop
558 END DO lists
559
560 !sample peak memory
561 CALL m_memory()
562
563 DEALLOCATE (mm_radius)
564 DEALLOCATE (qeff)
565 DEALLOCATE (qcore)
566 DEALLOCATE (qshell)
567 DEALLOCATE (is_shell_kind)
568
569 CALL timestop(handle)
570
571 END SUBROUTINE force_nonbond
572
573 ! **************************************************************************************************
574 !> \brief Adds a non-bonding contribution to the total force and optionally to
575 !> the virial.
576 ! **************************************************************************************************
577! **************************************************************************************************
578!> \brief ...
579!> \param f_nonbond_a ...
580!> \param f_nonbond_b ...
581!> \param pv ...
582!> \param fscalar ...
583!> \param rab ...
584!> \param use_virial ...
585! **************************************************************************************************
586 SUBROUTINE add_force_nonbond(f_nonbond_a, f_nonbond_b, pv, fscalar, rab, use_virial)
587
588 REAL(kind=dp), DIMENSION(3), INTENT(INOUT) :: f_nonbond_a, f_nonbond_b
589 REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv
590 REAL(kind=dp), INTENT(IN) :: fscalar
591 REAL(kind=dp), DIMENSION(3), INTENT(IN) :: rab
592 LOGICAL, INTENT(IN) :: use_virial
593
594 REAL(kind=dp), DIMENSION(3) :: fr
595
596 fr(1) = fscalar*rab(1)
597 fr(2) = fscalar*rab(2)
598 fr(3) = fscalar*rab(3)
599 f_nonbond_a(1) = f_nonbond_a(1) - fr(1)
600 f_nonbond_a(2) = f_nonbond_a(2) - fr(2)
601 f_nonbond_a(3) = f_nonbond_a(3) - fr(3)
602 f_nonbond_b(1) = f_nonbond_b(1) + fr(1)
603 f_nonbond_b(2) = f_nonbond_b(2) + fr(2)
604 f_nonbond_b(3) = f_nonbond_b(3) + fr(3)
605 IF (use_virial) THEN
606 pv(1, 1) = pv(1, 1) + rab(1)*fr(1)
607 pv(1, 2) = pv(1, 2) + rab(1)*fr(2)
608 pv(1, 3) = pv(1, 3) + rab(1)*fr(3)
609 pv(2, 1) = pv(2, 1) + rab(2)*fr(1)
610 pv(2, 2) = pv(2, 2) + rab(2)*fr(2)
611 pv(2, 3) = pv(2, 3) + rab(2)*fr(3)
612 pv(3, 1) = pv(3, 1) + rab(3)*fr(1)
613 pv(3, 2) = pv(3, 2) + rab(3)*fr(2)
614 pv(3, 3) = pv(3, 3) + rab(3)*fr(3)
615 END IF
616
617 END SUBROUTINE
618
619! **************************************************************************************************
620!> \brief corrects electrostatics for bonded terms
621!> \param fist_nonbond_env ...
622!> \param atomic_kind_set ...
623!> \param local_particles ...
624!> \param particle_set ...
625!> \param ewald_env ...
626!> \param v_bonded_corr ...
627!> \param pv_bc ...
628!> \param shell_particle_set ...
629!> \param core_particle_set ...
630!> \param atprop_env ...
631!> \param cell ...
632!> \param use_virial ...
633!> \par History
634!> Split routines to clean and to fix a bug with the tensor whose
635!> original definition was not correct for PBC.. [Teodoro Laino -06/2007]
636! **************************************************************************************************
637 SUBROUTINE bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
638 local_particles, particle_set, ewald_env, v_bonded_corr, pv_bc, &
639 shell_particle_set, core_particle_set, atprop_env, cell, use_virial)
640
641 TYPE(fist_nonbond_env_type), POINTER :: fist_nonbond_env
642 TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
643 TYPE(distribution_1d_type), POINTER :: local_particles
644 TYPE(particle_type), POINTER :: particle_set(:)
645 TYPE(ewald_environment_type), POINTER :: ewald_env
646 REAL(kind=dp), INTENT(OUT) :: v_bonded_corr
647 REAL(kind=dp), DIMENSION(:, :), INTENT(OUT) :: pv_bc
648 TYPE(particle_type), OPTIONAL, POINTER :: shell_particle_set(:), &
649 core_particle_set(:)
650 TYPE(atprop_type), POINTER :: atprop_env
651 TYPE(cell_type), POINTER :: cell
652 LOGICAL, INTENT(IN) :: use_virial
653
654 CHARACTER(LEN=*), PARAMETER :: routinen = 'bonded_correct_gaussian'
655
656 INTEGER :: atom_a, atom_b, handle, iatom, iend, igrp, ilist, ipair, istart, kind_a, kind_b, &
657 natoms_per_kind, nkind, npairs, shell_a, shell_b
658 INTEGER, DIMENSION(:, :), POINTER :: list
659 LOGICAL :: a_is_shell, b_is_shell, do_multipoles, &
660 full_nl, shell_adiabatic
661 REAL(kind=dp) :: alpha, const, fac_cor, fac_ei, qcore_a, &
662 qcore_b, qeff_a, qeff_b, qshell_a, &
663 qshell_b
664 REAL(kind=dp), DIMENSION(3) :: rca, rcb, rsa, rsb
665 REAL(kind=dp), DIMENSION(:, :), POINTER :: ij_kind_full_fac
666 TYPE(atomic_kind_type), POINTER :: atomic_kind
667 TYPE(fist_neighbor_type), POINTER :: nonbonded
668 TYPE(mp_comm_type) :: group
669 TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
670 TYPE(pair_potential_pp_type), POINTER :: potparm, potparm14
671 TYPE(pair_potential_single_type), POINTER :: pot
672 TYPE(shell_kind_type), POINTER :: shell_kind
673
674 CALL timeset(routinen, handle)
675
676 ! Initializing values
677 IF (use_virial) pv_bc = 0.0_dp
678 v_bonded_corr = 0.0_dp
679
680 CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, &
681 potparm14=potparm14, potparm=potparm, &
682 ij_kind_full_fac=ij_kind_full_fac)
683 CALL ewald_env_get(ewald_env, alpha=alpha, do_multipoles=do_multipoles, &
684 group=group)
685 ! Defining the constants
686 const = 2.0_dp*alpha*oorootpi
687
688 CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
689 shell_adiabatic=shell_adiabatic)
690
691 lists: DO ilist = 1, nonbonded%nlists
692 neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
693 npairs = neighbor_kind_pair%nscale
694 IF (npairs == 0) cycle
695 list => neighbor_kind_pair%list
696 kind_group_loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
697 istart = neighbor_kind_pair%grp_kind_start(igrp)
698 IF (istart > npairs) THEN
699 EXIT
700 END IF
701 iend = min(npairs, neighbor_kind_pair%grp_kind_end(igrp))
702
703 pairs: DO ipair = istart, iend
704 atom_a = list(1, ipair)
705 atom_b = list(2, ipair)
706 ! Get actual atomic kinds, since atom_a is not always of
707 ! kind_a and atom_b of kind_b, ie. they might be swapped.
708 kind_a = particle_set(atom_a)%atomic_kind%kind_number
709 kind_b = particle_set(atom_b)%atomic_kind%kind_number
710
711 ! take the proper potential, only for full_nl test
712 pot => potparm%pot(kind_a, kind_b)%pot
713 IF (ipair <= neighbor_kind_pair%nscale) THEN
714 IF (neighbor_kind_pair%is_onfo(ipair)) THEN
715 pot => potparm14%pot(kind_a, kind_b)%pot
716 END IF
717 END IF
718
719 ! Determine the scaling factors
720 fac_ei = ij_kind_full_fac(kind_a, kind_b)
721 full_nl = any(pot%type == tersoff_type) .OR. any(pot%type == siepmann_type) &
722 .OR. any(pot%type == gal_type) .OR. any(pot%type == gal21_type) &
723 .OR. any(pot%type == nequip_type) .OR. any(pot%type == allegro_type) &
724 .OR. any(pot%type == deepmd_type)
725 IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
726 fac_ei = fac_ei*0.5_dp
727 END IF
728 IF (ipair <= neighbor_kind_pair%nscale) THEN
729 fac_ei = fac_ei*neighbor_kind_pair%ei_scale(ipair)
730 END IF
731 ! The amount of correction is related to the
732 ! amount of scaling as follows:
733 fac_cor = 1.0_dp - fac_ei
734 IF (fac_cor <= 0.0_dp) cycle
735
736 ! Parameters for kind a
737 atomic_kind => atomic_kind_set(kind_a)
738 CALL get_atomic_kind(atomic_kind, qeff=qeff_a, shell=shell_kind)
739 IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_a = fist_nonbond_env%charges(atom_a)
740 a_is_shell = ASSOCIATED(shell_kind)
741 IF (a_is_shell) THEN
742 CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
743 charge_shell=qshell_a)
744 shell_a = particle_set(atom_a)%shell_index
745 rca = core_particle_set(shell_a)%r
746 rsa = shell_particle_set(shell_a)%r
747 ELSE
748 qcore_a = qeff_a
749 qshell_a = huge(0.0_dp)
750 shell_a = 0
751 rca = particle_set(atom_a)%r
752 rsa = 0.0_dp
753 END IF
754
755 ! Parameters for kind b
756 atomic_kind => atomic_kind_set(kind_b)
757 CALL get_atomic_kind(atomic_kind, qeff=qeff_b, shell=shell_kind)
758 IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_b = fist_nonbond_env%charges(atom_b)
759 b_is_shell = ASSOCIATED(shell_kind)
760 IF (b_is_shell) THEN
761 CALL get_shell(shell=shell_kind, charge_core=qcore_b, &
762 charge_shell=qshell_b)
763 shell_b = particle_set(atom_b)%shell_index
764 rcb = core_particle_set(shell_b)%r
765 rsb = shell_particle_set(shell_b)%r
766 ELSE
767 qcore_b = qeff_b
768 qshell_b = huge(0.0_dp)
769 shell_b = 0
770 rcb = particle_set(atom_b)%r
771 rsb = 0.0_dp
772 END IF
773
774 ! First part: take care of core/ion-core/ion correction
775 IF (a_is_shell .AND. b_is_shell) THEN
776 ! correct for core-core interaction
777 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
778 v_bonded_corr, core_particle_set, core_particle_set, &
779 shell_a, shell_b, .true., alpha, qcore_a, qcore_b, &
780 const, fac_cor, pv_bc, atprop_env, use_virial)
781 ELSE IF (a_is_shell) THEN
782 ! correct for core-ion interaction
783 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
784 v_bonded_corr, core_particle_set, particle_set, &
785 shell_a, atom_b, .true., alpha, qcore_a, qcore_b, &
786 const, fac_cor, pv_bc, atprop_env, use_virial)
787 ELSE IF (b_is_shell) THEN
788 ! correct for ion-core interaction
789 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
790 v_bonded_corr, particle_set, core_particle_set, &
791 atom_a, shell_b, .true., alpha, qcore_a, qcore_b, &
792 const, fac_cor, pv_bc, atprop_env, use_virial)
793 ELSE
794 ! correct for ion-ion interaction
795 CALL bonded_correct_gaussian_low(rca, rcb, cell, &
796 v_bonded_corr, particle_set, particle_set, &
797 atom_a, atom_b, .true., alpha, qcore_a, qcore_b, &
798 const, fac_cor, pv_bc, atprop_env, use_virial)
799 END IF
800
801 ! Second part: take care of shell-shell, shell-core/ion and
802 ! core/ion-shell corrections
803 IF (a_is_shell .AND. b_is_shell) THEN
804 ! correct for shell-shell interaction
805 CALL bonded_correct_gaussian_low(rsa, rsa, cell, &
806 v_bonded_corr, shell_particle_set, shell_particle_set, &
807 shell_a, shell_b, shell_adiabatic, alpha, qshell_a, &
808 qshell_b, const, fac_cor, pv_bc, atprop_env, use_virial)
809 END IF
810 IF (a_is_shell) THEN
811 IF (b_is_shell) THEN
812 ! correct for shell-core interaction
813 CALL bonded_correct_gaussian_low(rsa, rcb, cell, &
814 v_bonded_corr, shell_particle_set, core_particle_set, &
815 shell_a, shell_b, shell_adiabatic, alpha, qshell_a, qcore_b, &
816 const, fac_cor, pv_bc, atprop_env, use_virial)
817 ELSE
818 ! correct for shell-ion interaction
819 CALL bonded_correct_gaussian_low(rsa, rcb, cell, &
820 v_bonded_corr, shell_particle_set, particle_set, &
821 shell_a, atom_b, shell_adiabatic, alpha, qshell_a, qcore_b, &
822 const, fac_cor, pv_bc, atprop_env, use_virial)
823 END IF
824 END IF
825 IF (b_is_shell) THEN
826 IF (a_is_shell) THEN
827 ! correct for core-shell interaction
828 CALL bonded_correct_gaussian_low(rca, rsb, cell, &
829 v_bonded_corr, core_particle_set, shell_particle_set, &
830 shell_a, shell_b, shell_adiabatic, alpha, qcore_a, qshell_b, &
831 const, fac_cor, pv_bc, atprop_env, use_virial)
832 ELSE
833 ! correct for ion-shell interaction
834 CALL bonded_correct_gaussian_low(rca, rsb, cell, &
835 v_bonded_corr, particle_set, shell_particle_set, &
836 atom_a, shell_b, shell_adiabatic, alpha, qcore_a, qshell_b, &
837 const, fac_cor, pv_bc, atprop_env, use_virial)
838 END IF
839 END IF
840 END DO pairs
841 END DO kind_group_loop
842 END DO lists
843
844 ! Always correct core-shell interaction within one atom.
845 nkind = SIZE(atomic_kind_set)
846 DO kind_a = 1, nkind
847 ! parameters for kind a
848 atomic_kind => atomic_kind_set(kind_a)
849 CALL get_atomic_kind(atomic_kind, shell=shell_kind)
850 IF (ASSOCIATED(shell_kind)) THEN
851 CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
852 charge_shell=qshell_a)
853
854 natoms_per_kind = local_particles%n_el(kind_a)
855 DO iatom = 1, natoms_per_kind
856
857 ! Data for atom a
858 atom_a = local_particles%list(kind_a)%array(iatom)
859 shell_a = particle_set(atom_a)%shell_index
860 rca = core_particle_set(shell_a)%r
861 rsa = shell_particle_set(shell_a)%r
862
863 CALL bonded_correct_gaussian_low_sh(rca, rsa, cell, &
864 v_bonded_corr, core_particle_set, shell_particle_set, &
865 shell_a, shell_adiabatic, alpha, qcore_a, qshell_a, &
866 const, pv_bc, atprop_env, use_virial)
867
868 END DO
869 END IF
870 END DO
871
872 CALL group%sum(v_bonded_corr)
873
874 CALL timestop(handle)
875
876 END SUBROUTINE bonded_correct_gaussian
877
878! **************************************************************************************************
879!> \brief ...
880!> \param r1 ...
881!> \param r2 ...
882!> \param cell ...
883!> \param v_bonded_corr ...
884!> \param particle_set1 ...
885!> \param particle_set2 ...
886!> \param i ...
887!> \param j ...
888!> \param shell_adiabatic ...
889!> \param alpha ...
890!> \param q1 ...
891!> \param q2 ...
892!> \param const ...
893!> \param fac_cor ...
894!> \param pv_bc ...
895!> \param atprop_env ...
896!> \param use_virial ...
897!> \par History
898!> Split routines to clean and to fix a bug with the tensor whose
899!> original definition was not correct for PBC..
900!> \author Teodoro Laino
901! **************************************************************************************************
902 SUBROUTINE bonded_correct_gaussian_low(r1, r2, cell, v_bonded_corr, &
903 particle_set1, particle_set2, i, j, shell_adiabatic, alpha, q1, q2, &
904 const, fac_cor, pv_bc, atprop_env, use_virial)
905 REAL(kind=dp), DIMENSION(3) :: r1, r2
906 TYPE(cell_type), POINTER :: cell
907 REAL(kind=dp), INTENT(INOUT) :: v_bonded_corr
908 TYPE(particle_type), POINTER :: particle_set1(:), particle_set2(:)
909 INTEGER, INTENT(IN) :: i, j
910 LOGICAL, INTENT(IN) :: shell_adiabatic
911 REAL(kind=dp), INTENT(IN) :: alpha, q1, q2, const, fac_cor
912 REAL(kind=dp), INTENT(INOUT) :: pv_bc(3, 3)
913 TYPE(atprop_type), POINTER :: atprop_env
914 LOGICAL, INTENT(IN) :: use_virial
915
916 REAL(kind=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
917 ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, pc = 0.3275911_dp
918
919 INTEGER :: iatom, jatom
920 REAL(kind=dp) :: arg, dij, e_arg_arg, errf, fscalar, &
921 idij, rijsq, tc, vbc
922 REAL(kind=dp), DIMENSION(3) :: fij_com, rij
923 REAL(kind=dp), DIMENSION(3, 3) :: fbc
924
925 rij = r1 - r2
926 rij = pbc(rij, cell)
927 rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
928 idij = 1.0_dp/sqrt(rijsq)
929 dij = rijsq*idij
930 arg = alpha*dij
931 e_arg_arg = exp(-arg**2)
932 tc = 1.0_dp/(1.0_dp + pc*arg)
933
934 ! Defining errf=1-erfc
935 errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
936
937 ! Getting the potential
938 vbc = -q1*q2*idij*errf*fac_cor
939 v_bonded_corr = v_bonded_corr + vbc
940 IF (atprop_env%energy) THEN
941 iatom = particle_set1(i)%atom_index
942 atprop_env%atener(iatom) = atprop_env%atener(iatom) + 0.5_dp*vbc
943 jatom = particle_set2(j)%atom_index
944 atprop_env%atener(jatom) = atprop_env%atener(jatom) + 0.5_dp*vbc
945 END IF
946
947 ! Subtracting the force from the total force
948 fscalar = q1*q2*idij**2*(idij*errf - const*e_arg_arg)*fac_cor
949
950 particle_set1(i)%f(1) = particle_set1(i)%f(1) - fscalar*rij(1)
951 particle_set1(i)%f(2) = particle_set1(i)%f(2) - fscalar*rij(2)
952 particle_set1(i)%f(3) = particle_set1(i)%f(3) - fscalar*rij(3)
953
954 particle_set2(j)%f(1) = particle_set2(j)%f(1) + fscalar*rij(1)
955 particle_set2(j)%f(2) = particle_set2(j)%f(2) + fscalar*rij(2)
956 particle_set2(j)%f(3) = particle_set2(j)%f(3) + fscalar*rij(3)
957
958 IF (use_virial .AND. shell_adiabatic) THEN
959 fij_com = fscalar*rij
960 fbc(1, 1) = -fij_com(1)*rij(1)
961 fbc(1, 2) = -fij_com(1)*rij(2)
962 fbc(1, 3) = -fij_com(1)*rij(3)
963 fbc(2, 1) = -fij_com(2)*rij(1)
964 fbc(2, 2) = -fij_com(2)*rij(2)
965 fbc(2, 3) = -fij_com(2)*rij(3)
966 fbc(3, 1) = -fij_com(3)*rij(1)
967 fbc(3, 2) = -fij_com(3)*rij(2)
968 fbc(3, 3) = -fij_com(3)*rij(3)
969 pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
970 END IF
971
972 END SUBROUTINE bonded_correct_gaussian_low
973
974! **************************************************************************************************
975!> \brief specific for shell models cleans the interaction core-shell on the same
976!> atom
977!> \param r1 ...
978!> \param r2 ...
979!> \param cell ...
980!> \param v_bonded_corr ...
981!> \param core_particle_set ...
982!> \param shell_particle_set ...
983!> \param i ...
984!> \param shell_adiabatic ...
985!> \param alpha ...
986!> \param q1 ...
987!> \param q2 ...
988!> \param const ...
989!> \param pv_bc ...
990!> \param atprop_env ...
991!> \param use_virial ...
992!> \par History
993!> Split routines to clean and to fix a bug with the tensor whose
994!> original definition was not correct for PBC..
995!> \author Teodoro Laino
996! **************************************************************************************************
997 SUBROUTINE bonded_correct_gaussian_low_sh(r1, r2, cell, v_bonded_corr, &
998 core_particle_set, shell_particle_set, i, shell_adiabatic, alpha, q1, q2, &
999 const, pv_bc, atprop_env, use_virial)
1000 REAL(kind=dp), DIMENSION(3) :: r1, r2
1001 TYPE(cell_type), POINTER :: cell
1002 REAL(kind=dp), INTENT(INOUT) :: v_bonded_corr
1003 TYPE(particle_type), POINTER :: core_particle_set(:), &
1004 shell_particle_set(:)
1005 INTEGER, INTENT(IN) :: i
1006 LOGICAL, INTENT(IN) :: shell_adiabatic
1007 REAL(kind=dp), INTENT(IN) :: alpha, q1, q2, const
1008 REAL(kind=dp), INTENT(INOUT) :: pv_bc(3, 3)
1009 TYPE(atprop_type), POINTER :: atprop_env
1010 LOGICAL, INTENT(IN) :: use_virial
1011
1012 REAL(kind=dp), PARAMETER :: ac1 = 0.254829592_dp, ac2 = -0.284496736_dp, &
1013 ac3 = 1.421413741_dp, ac4 = -1.453152027_dp, ac5 = 1.061405429_dp, pc = 0.3275911_dp
1014
1015 INTEGER :: iatom
1016 REAL(kind=dp) :: arg, dij, e_arg_arg, efac, errf, ffac, &
1017 fscalar, idij, rijsq, tc, tc2, tc4, vbc
1018 REAL(kind=dp), DIMENSION(3) :: fr, rij
1019 REAL(kind=dp), DIMENSION(3, 3) :: fbc
1020
1021 rij = r1 - r2
1022 rij = pbc(rij, cell)
1023 rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
1024 dij = sqrt(rijsq)
1025 ! Two possible limiting cases according the value of dij
1026 arg = alpha*dij
1027 ! and this is a magic number.. it is related to the order expansion
1028 ! and to the value of the polynomial coefficients
1029 IF (arg > 0.355_dp) THEN
1030 idij = 1.0_dp/dij
1031 e_arg_arg = exp(-arg*arg)
1032 tc = 1.0_dp/(1.0_dp + pc*arg)
1033 ! defining errf = 1 - erfc
1034 errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
1035 efac = idij*errf
1036 ffac = idij**2*(efac - const*e_arg_arg)
1037 ELSE
1038 tc = arg*arg
1039 tc2 = tc*tc
1040 tc4 = tc2*tc2
1041 efac = const*(1.0_dp - tc/3.0_dp + tc2/10.0_dp - tc*tc2/42.0_dp + tc4/216.0_dp - &
1042 tc*tc4/1320.0_dp + tc2*tc4/9360.0_dp)
1043 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 + &
1044 tc4/132.0_dp - tc*tc4/780.0_dp)
1045 END IF
1046
1047 ! getting the potential
1048 vbc = -q1*q2*efac
1049 v_bonded_corr = v_bonded_corr + vbc
1050 IF (atprop_env%energy) THEN
1051 iatom = shell_particle_set(i)%atom_index
1052 atprop_env%atener(iatom) = atprop_env%atener(iatom) + vbc
1053 END IF
1054
1055 ! subtracting the force from the total force
1056 fscalar = q1*q2*ffac
1057 fr(:) = fscalar*rij(:)
1058
1059 core_particle_set(i)%f(1) = core_particle_set(i)%f(1) - fr(1)
1060 core_particle_set(i)%f(2) = core_particle_set(i)%f(2) - fr(2)
1061 core_particle_set(i)%f(3) = core_particle_set(i)%f(3) - fr(3)
1062
1063 shell_particle_set(i)%f(1) = shell_particle_set(i)%f(1) + fr(1)
1064 shell_particle_set(i)%f(2) = shell_particle_set(i)%f(2) + fr(2)
1065 shell_particle_set(i)%f(3) = shell_particle_set(i)%f(3) + fr(3)
1066
1067 IF (use_virial .AND. shell_adiabatic) THEN
1068 fbc(1, 1) = -fr(1)*rij(1)
1069 fbc(1, 2) = -fr(1)*rij(2)
1070 fbc(1, 3) = -fr(1)*rij(3)
1071 fbc(2, 1) = -fr(2)*rij(1)
1072 fbc(2, 2) = -fr(2)*rij(2)
1073 fbc(2, 3) = -fr(2)*rij(3)
1074 fbc(3, 1) = -fr(3)*rij(1)
1075 fbc(3, 2) = -fr(3)*rij(2)
1076 fbc(3, 3) = -fr(3)*rij(3)
1077 pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
1078 END IF
1079
1080 END SUBROUTINE bonded_correct_gaussian_low_sh
1081
1082END 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, 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:640
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: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: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: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:416
pair_potential_types::pair_potential_single_type
Definition pair_potential_types.F:383
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