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