(git:1f285aa)
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 ! **************************************************************************************************
17 MODULE 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 
61 CONTAINS
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 
1103 END MODULE fist_nonbond_force
pbc
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
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