Line data Source code
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 82936 : SUBROUTINE force_nonbond(fist_nonbond_env, ewald_env, particle_set, cell, &
80 82936 : 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 82936 : INTEGER, DIMENSION(:, :), POINTER :: list
100 : LOGICAL :: all_terms, do_multipoles, full_nl, &
101 : shell_present
102 82936 : 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 82936 : 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 82936 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ij_kind_full_fac
114 82936 : 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 82936 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc, &
122 82936 : rcore_last_update_pbc, &
123 82936 : rshell_last_update_pbc
124 : TYPE(shell_kind_type), POINTER :: shell_kind
125 82936 : TYPE(spline_data_p_type), DIMENSION(:), POINTER :: spline_data
126 : TYPE(spline_factor_type), POINTER :: spl_f
127 :
128 82936 : CALL timeset(routineN, handle)
129 82936 : NULLIFY (logger)
130 82936 : logger => cp_get_default_logger()
131 82936 : 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 82936 : 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 82936 : interaction_cutoffs=ei_interaction_cutoffs)
141 :
142 : ! Initializing the potential energy, pressure tensor and force
143 82936 : pot_nonbond = 0.0_dp
144 33086652 : f_nonbond(:, :) = 0.0_dp
145 :
146 82936 : IF (use_virial) THEN
147 221884 : pv_nonbond(:, :) = 0.0_dp
148 : END IF
149 82936 : shell_present = .FALSE.
150 82936 : IF (PRESENT(fshell_nonbond)) THEN
151 9364 : CPASSERT(PRESENT(fcore_nonbond))
152 3057596 : fshell_nonbond = 0.0_dp
153 3057596 : fcore_nonbond = 0.0_dp
154 : shell_present = .TRUE.
155 : END IF
156 : ! Load atomic kind information
157 248808 : ALLOCATE (mm_radius(nkind))
158 248808 : ALLOCATE (qeff(nkind))
159 248808 : ALLOCATE (qcore(nkind))
160 248808 : ALLOCATE (qshell(nkind))
161 248808 : ALLOCATE (is_shell_kind(nkind))
162 316021 : DO ikind = 1, nkind
163 233085 : atomic_kind => atomic_kind_set(ikind)
164 : CALL get_atomic_kind(atomic_kind, &
165 : qeff=qeff(ikind), &
166 : mm_radius=mm_radius(ikind), &
167 233085 : shell=shell_kind)
168 233085 : is_shell_kind(ikind) = ASSOCIATED(shell_kind)
169 316021 : IF (ASSOCIATED(shell_kind)) THEN
170 : CALL get_shell(shell=shell_kind, &
171 : charge_core=qcore(ikind), &
172 15798 : charge_shell=qshell(ikind))
173 : ELSE
174 217287 : qcore(ikind) = 0.0_dp
175 217287 : qshell(ikind) = 0.0_dp
176 : END IF
177 : END DO
178 : ! Starting the force loop
179 9872066 : Lists: DO ilist = 1, nonbonded%nlists
180 9789130 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
181 9789130 : npairs = neighbor_kind_pair%npairs
182 9789130 : IF (npairs == 0) CYCLE
183 2630036 : list => neighbor_kind_pair%list
184 10520144 : cvi = neighbor_kind_pair%cell_vector
185 34190468 : cell_v = MATMUL(cell%hmat, cvi)
186 11259283 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
187 8546311 : istart = neighbor_kind_pair%grp_kind_start(igrp)
188 8546311 : 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 18335441 : !$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 82936 : CALL m_memory()
571 :
572 82936 : DEALLOCATE (mm_radius)
573 82936 : DEALLOCATE (qeff)
574 82936 : DEALLOCATE (qcore)
575 82936 : DEALLOCATE (qshell)
576 82936 : DEALLOCATE (is_shell_kind)
577 :
578 82936 : CALL timestop(handle)
579 :
580 248808 : 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 1006726708 : 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 1006726708 : fr(1) = fscalar*rab(1)
606 1006726708 : fr(2) = fscalar*rab(2)
607 1006726708 : fr(3) = fscalar*rab(3)
608 1006726708 : f_nonbond_a(1) = f_nonbond_a(1) - fr(1)
609 1006726708 : f_nonbond_a(2) = f_nonbond_a(2) - fr(2)
610 1006726708 : f_nonbond_a(3) = f_nonbond_a(3) - fr(3)
611 1006726708 : f_nonbond_b(1) = f_nonbond_b(1) + fr(1)
612 1006726708 : f_nonbond_b(2) = f_nonbond_b(2) + fr(2)
613 1006726708 : f_nonbond_b(3) = f_nonbond_b(3) + fr(3)
614 1006726708 : IF (use_virial) THEN
615 361579480 : pv(1, 1) = pv(1, 1) + rab(1)*fr(1)
616 361579480 : pv(1, 2) = pv(1, 2) + rab(1)*fr(2)
617 361579480 : pv(1, 3) = pv(1, 3) + rab(1)*fr(3)
618 361579480 : pv(2, 1) = pv(2, 1) + rab(2)*fr(1)
619 361579480 : pv(2, 2) = pv(2, 2) + rab(2)*fr(2)
620 361579480 : pv(2, 3) = pv(2, 3) + rab(2)*fr(3)
621 361579480 : pv(3, 1) = pv(3, 1) + rab(3)*fr(1)
622 361579480 : pv(3, 2) = pv(3, 2) + rab(3)*fr(2)
623 361579480 : pv(3, 3) = pv(3, 3) + rab(3)*fr(3)
624 : END IF
625 :
626 1006726708 : 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 238984 : SUBROUTINE bonded_correct_gaussian(fist_nonbond_env, atomic_kind_set, &
647 59746 : 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 59746 : 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 59746 : 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 59746 : CALL timeset(routineN, handle)
684 :
685 : ! Initializing values
686 240290 : IF (use_virial) pv_bc = 0.0_dp
687 59746 : v_bonded_corr = 0.0_dp
688 :
689 : CALL fist_nonbond_env_get(fist_nonbond_env, nonbonded=nonbonded, &
690 : potparm14=potparm14, potparm=potparm, &
691 59746 : ij_kind_full_fac=ij_kind_full_fac)
692 : CALL ewald_env_get(ewald_env, alpha=alpha, do_multipoles=do_multipoles, &
693 59746 : group=group)
694 : ! Defining the constants
695 59746 : const = 2.0_dp*alpha*oorootpi
696 :
697 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
698 59746 : shell_adiabatic=shell_adiabatic)
699 :
700 4996200 : Lists: DO ilist = 1, nonbonded%nlists
701 4936454 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
702 4936454 : npairs = neighbor_kind_pair%nscale
703 4936454 : IF (npairs == 0) CYCLE
704 66091 : list => neighbor_kind_pair%list
705 2264856 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
706 2190303 : istart = neighbor_kind_pair%grp_kind_start(igrp)
707 2190303 : IF (istart > npairs) THEN
708 : EXIT
709 : END IF
710 2139019 : iend = MIN(npairs, neighbor_kind_pair%grp_kind_end(igrp))
711 :
712 12106077 : Pairs: DO ipair = istart, iend
713 5030604 : atom_a = list(1, ipair)
714 5030604 : 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 5030604 : kind_a = particle_set(atom_a)%atomic_kind%kind_number
718 5030604 : kind_b = particle_set(atom_b)%atomic_kind%kind_number
719 :
720 : ! take the proper potential, only for full_nl test
721 5030604 : pot => potparm%pot(kind_a, kind_b)%pot
722 5030604 : IF (ipair <= neighbor_kind_pair%nscale) THEN
723 5030604 : IF (neighbor_kind_pair%is_onfo(ipair)) THEN
724 873150 : pot => potparm14%pot(kind_a, kind_b)%pot
725 : END IF
726 : END IF
727 :
728 : ! Determine the scaling factors
729 5030604 : 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 70428456 : .OR. ANY(pot%type == deepmd_type)
734 5030604 : IF ((.NOT. full_nl) .AND. (atom_a == atom_b)) THEN
735 0 : fac_ei = fac_ei*0.5_dp
736 : END IF
737 5030604 : IF (ipair <= neighbor_kind_pair%nscale) THEN
738 5030604 : 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 5030604 : fac_cor = 1.0_dp - fac_ei
743 5030604 : IF (fac_cor <= 0.0_dp) CYCLE
744 :
745 : ! Parameters for kind a
746 5028223 : atomic_kind => atomic_kind_set(kind_a)
747 5028223 : CALL get_atomic_kind(atomic_kind, qeff=qeff_a, shell=shell_kind)
748 5028223 : IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_a = fist_nonbond_env%charges(atom_a)
749 5028223 : a_is_shell = ASSOCIATED(shell_kind)
750 5028223 : IF (a_is_shell) THEN
751 : CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
752 8 : charge_shell=qshell_a)
753 8 : shell_a = particle_set(atom_a)%shell_index
754 32 : rca = core_particle_set(shell_a)%r
755 32 : rsa = shell_particle_set(shell_a)%r
756 : ELSE
757 5028215 : qcore_a = qeff_a
758 5028215 : qshell_a = HUGE(0.0_dp)
759 5028215 : shell_a = 0
760 20112860 : rca = particle_set(atom_a)%r
761 5028215 : rsa = 0.0_dp
762 : END IF
763 :
764 : ! Parameters for kind b
765 5028223 : atomic_kind => atomic_kind_set(kind_b)
766 5028223 : CALL get_atomic_kind(atomic_kind, qeff=qeff_b, shell=shell_kind)
767 5028223 : IF (ASSOCIATED(fist_nonbond_env%charges)) qeff_b = fist_nonbond_env%charges(atom_b)
768 5028223 : b_is_shell = ASSOCIATED(shell_kind)
769 5028223 : IF (b_is_shell) THEN
770 : CALL get_shell(shell=shell_kind, charge_core=qcore_b, &
771 492 : charge_shell=qshell_b)
772 492 : shell_b = particle_set(atom_b)%shell_index
773 1968 : rcb = core_particle_set(shell_b)%r
774 1968 : rsb = shell_particle_set(shell_b)%r
775 : ELSE
776 5027731 : qcore_b = qeff_b
777 5027731 : qshell_b = HUGE(0.0_dp)
778 5027731 : shell_b = 0
779 20110924 : rcb = particle_set(atom_b)%r
780 5027731 : rsb = 0.0_dp
781 : END IF
782 :
783 : ! First part: take care of core/ion-core/ion correction
784 5028223 : 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 0 : const, fac_cor, pv_bc, atprop_env, use_virial)
790 5028223 : 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 8 : const, fac_cor, pv_bc, atprop_env, use_virial)
796 5028215 : 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 492 : 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 5027723 : 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 5028223 : 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 0 : qshell_b, const, fac_cor, pv_bc, atprop_env, use_virial)
818 : END IF
819 5028223 : IF (a_is_shell) THEN
820 8 : 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 0 : 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 8 : const, fac_cor, pv_bc, atprop_env, use_virial)
832 : END IF
833 : END IF
834 17223688 : IF (b_is_shell) THEN
835 492 : 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 0 : 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 492 : 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 59746 : nkind = SIZE(atomic_kind_set)
855 264102 : DO kind_a = 1, nkind
856 : ! parameters for kind a
857 204356 : atomic_kind => atomic_kind_set(kind_a)
858 204356 : CALL get_atomic_kind(atomic_kind, shell=shell_kind)
859 264102 : IF (ASSOCIATED(shell_kind)) THEN
860 : CALL get_shell(shell=shell_kind, charge_core=qcore_a, &
861 15778 : charge_shell=qshell_a)
862 :
863 15778 : natoms_per_kind = local_particles%n_el(kind_a)
864 431567 : DO iatom = 1, natoms_per_kind
865 :
866 : ! Data for atom a
867 415789 : atom_a = local_particles%list(kind_a)%array(iatom)
868 415789 : shell_a = particle_set(atom_a)%shell_index
869 1663156 : rca = core_particle_set(shell_a)%r
870 1663156 : 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 431567 : const, pv_bc, atprop_env, use_virial)
876 :
877 : END DO
878 : END IF
879 : END DO
880 :
881 59746 : CALL group%sum(v_bonded_corr)
882 :
883 59746 : CALL timestop(handle)
884 :
885 59746 : 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 5028723 : 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 20114892 : rij = r1 - r2
935 20114892 : rij = pbc(rij, cell)
936 5028723 : rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
937 5028723 : idij = 1.0_dp/SQRT(rijsq)
938 5028723 : dij = rijsq*idij
939 5028723 : arg = alpha*dij
940 5028723 : e_arg_arg = EXP(-arg**2)
941 5028723 : tc = 1.0_dp/(1.0_dp + pc*arg)
942 :
943 : ! Defining errf=1-erfc
944 5028723 : errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
945 :
946 : ! Getting the potential
947 5028723 : vbc = -q1*q2*idij*errf*fac_cor
948 5028723 : v_bonded_corr = v_bonded_corr + vbc
949 5028723 : IF (atprop_env%energy) THEN
950 909 : iatom = particle_set1(i)%atom_index
951 909 : atprop_env%atener(iatom) = atprop_env%atener(iatom) + 0.5_dp*vbc
952 909 : jatom = particle_set2(j)%atom_index
953 909 : 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 5028723 : fscalar = q1*q2*idij**2*(idij*errf - const*e_arg_arg)*fac_cor
958 :
959 5028723 : particle_set1(i)%f(1) = particle_set1(i)%f(1) - fscalar*rij(1)
960 5028723 : particle_set1(i)%f(2) = particle_set1(i)%f(2) - fscalar*rij(2)
961 5028723 : particle_set1(i)%f(3) = particle_set1(i)%f(3) - fscalar*rij(3)
962 :
963 5028723 : particle_set2(j)%f(1) = particle_set2(j)%f(1) + fscalar*rij(1)
964 5028723 : particle_set2(j)%f(2) = particle_set2(j)%f(2) + fscalar*rij(2)
965 5028723 : particle_set2(j)%f(3) = particle_set2(j)%f(3) + fscalar*rij(3)
966 :
967 5028723 : IF (use_virial .AND. shell_adiabatic) THEN
968 2130648 : fij_com = fscalar*rij
969 532662 : fbc(1, 1) = -fij_com(1)*rij(1)
970 532662 : fbc(1, 2) = -fij_com(1)*rij(2)
971 532662 : fbc(1, 3) = -fij_com(1)*rij(3)
972 532662 : fbc(2, 1) = -fij_com(2)*rij(1)
973 532662 : fbc(2, 2) = -fij_com(2)*rij(2)
974 532662 : fbc(2, 3) = -fij_com(2)*rij(3)
975 532662 : fbc(3, 1) = -fij_com(3)*rij(1)
976 532662 : fbc(3, 2) = -fij_com(3)*rij(2)
977 532662 : fbc(3, 3) = -fij_com(3)*rij(3)
978 6924606 : pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
979 532662 : IF (atprop_env%stress) THEN
980 : ! Atomic stress tensors
981 909 : iatom = particle_set1(i)%atom_index
982 11817 : atprop_env%atstress(:, :, iatom) = atprop_env%atstress(:, :, iatom) + 0.5_dp*fbc(:, :)
983 909 : jatom = particle_set2(j)%atom_index
984 11817 : atprop_env%atstress(:, :, jatom) = atprop_env%atstress(:, :, jatom) + 0.5_dp*fbc(:, :)
985 : END IF
986 : END IF
987 :
988 5028723 : 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 415789 : 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 1663156 : rij = r1 - r2
1038 1663156 : rij = pbc(rij, cell)
1039 415789 : rijsq = rij(1)*rij(1) + rij(2)*rij(2) + rij(3)*rij(3)
1040 415789 : dij = SQRT(rijsq)
1041 : ! Two possible limiting cases according the value of dij
1042 415789 : 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 415789 : IF (arg > 0.355_dp) THEN
1046 0 : idij = 1.0_dp/dij
1047 0 : e_arg_arg = EXP(-arg*arg)
1048 0 : tc = 1.0_dp/(1.0_dp + pc*arg)
1049 : ! defining errf = 1 - erfc
1050 0 : errf = 1.0_dp - ((((ac5*tc + ac4)*tc + ac3)*tc + ac2)*tc + ac1)*tc*e_arg_arg
1051 0 : efac = idij*errf
1052 0 : ffac = idij**2*(efac - const*e_arg_arg)
1053 : ELSE
1054 415789 : tc = arg*arg
1055 415789 : tc2 = tc*tc
1056 415789 : 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 415789 : 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 415789 : tc4/132.0_dp - tc*tc4/780.0_dp)
1061 : END IF
1062 :
1063 : ! getting the potential
1064 415789 : vbc = -q1*q2*efac
1065 415789 : v_bonded_corr = v_bonded_corr + vbc
1066 415789 : IF (atprop_env%energy) THEN
1067 1080 : iatom = shell_particle_set(i)%atom_index
1068 1080 : atprop_env%atener(iatom) = atprop_env%atener(iatom) + vbc
1069 : END IF
1070 :
1071 : ! subtracting the force from the total force
1072 415789 : fscalar = q1*q2*ffac
1073 1663156 : fr(:) = fscalar*rij(:)
1074 :
1075 415789 : core_particle_set(i)%f(1) = core_particle_set(i)%f(1) - fr(1)
1076 415789 : core_particle_set(i)%f(2) = core_particle_set(i)%f(2) - fr(2)
1077 415789 : core_particle_set(i)%f(3) = core_particle_set(i)%f(3) - fr(3)
1078 :
1079 415789 : shell_particle_set(i)%f(1) = shell_particle_set(i)%f(1) + fr(1)
1080 415789 : shell_particle_set(i)%f(2) = shell_particle_set(i)%f(2) + fr(2)
1081 415789 : shell_particle_set(i)%f(3) = shell_particle_set(i)%f(3) + fr(3)
1082 :
1083 415789 : IF (use_virial .AND. shell_adiabatic) THEN
1084 341710 : fbc(1, 1) = -fr(1)*rij(1)
1085 341710 : fbc(1, 2) = -fr(1)*rij(2)
1086 341710 : fbc(1, 3) = -fr(1)*rij(3)
1087 341710 : fbc(2, 1) = -fr(2)*rij(1)
1088 341710 : fbc(2, 2) = -fr(2)*rij(2)
1089 341710 : fbc(2, 3) = -fr(2)*rij(3)
1090 341710 : fbc(3, 1) = -fr(3)*rij(1)
1091 341710 : fbc(3, 2) = -fr(3)*rij(2)
1092 341710 : fbc(3, 3) = -fr(3)*rij(3)
1093 4442230 : pv_bc(:, :) = pv_bc(:, :) + fbc(:, :)
1094 341710 : IF (atprop_env%stress) THEN
1095 : ! Atomic stress tensors
1096 0 : iatom = shell_particle_set(i)%atom_index
1097 0 : atprop_env%atstress(:, :, iatom) = atprop_env%atstress(:, :, iatom) + fbc(:, :)
1098 : END IF
1099 : END IF
1100 :
1101 415789 : END SUBROUTINE bonded_correct_gaussian_low_sh
1102 :
1103 : END MODULE fist_nonbond_force
|
