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