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