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 : !> \brief Treats the electrostatic for multipoles (up to quadrupoles)
10 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
11 : !> inclusion of optional electric field damping for the polarizable atoms
12 : !> Rodolphe Vuilleumier and Mathieu Salanne - 12.2009
13 : ! **************************************************************************************************
14 : MODULE ewalds_multipole
15 : USE atomic_kind_types, ONLY: atomic_kind_type
16 : USE bibliography, ONLY: Aguado2003, &
17 : Laino2008, &
18 : cite_reference
19 : USE cell_types, ONLY: cell_type, &
20 : pbc
21 : USE cp_log_handling, ONLY: cp_get_default_logger, &
22 : cp_logger_type
23 : USE damping_dipole_types, ONLY: damping_type, &
24 : no_damping, &
25 : tang_toennies
26 : USE dg_rho0_types, ONLY: dg_rho0_type
27 : USE dg_types, ONLY: dg_get, &
28 : dg_type
29 : USE distribution_1d_types, ONLY: distribution_1d_type
30 : USE ewald_environment_types, ONLY: ewald_env_get, &
31 : ewald_environment_type
32 : USE ewald_pw_types, ONLY: ewald_pw_get, &
33 : ewald_pw_type
34 : USE fist_neighbor_list_control, ONLY: list_control
35 : USE fist_neighbor_list_types, ONLY: fist_neighbor_type, &
36 : neighbor_kind_pairs_type
37 : USE fist_nonbond_env_types, ONLY: fist_nonbond_env_get, &
38 : fist_nonbond_env_type, &
39 : pos_type
40 : USE input_section_types, ONLY: section_vals_type
41 : USE kinds, ONLY: dp
42 : USE mathconstants, ONLY: fourpi, &
43 : oorootpi, &
44 : pi, &
45 : sqrthalf, &
46 : z_zero
47 : USE message_passing, ONLY: mp_comm_type
48 : USE parallel_rng_types, ONLY: UNIFORM, &
49 : rng_stream_type
50 : USE particle_types, ONLY: particle_type
51 : USE pw_grid_types, ONLY: pw_grid_type
52 : USE pw_pool_types, ONLY: pw_pool_type
53 : USE structure_factor_types, ONLY: structure_factor_type
54 : USE structure_factors, ONLY: structure_factor_allocate, &
55 : structure_factor_deallocate, &
56 : structure_factor_evaluate
57 : #include "./base/base_uses.f90"
58 :
59 : #:include "ewalds_multipole_sr.fypp"
60 :
61 : IMPLICIT NONE
62 : PRIVATE
63 :
64 : TYPE charge_mono_type
65 : REAL(KIND=dp), DIMENSION(:), &
66 : POINTER :: charge => NULL()
67 : REAL(KIND=dp), DIMENSION(:, :), &
68 : POINTER :: pos => NULL()
69 : END TYPE charge_mono_type
70 : TYPE multi_charge_type
71 : TYPE(charge_mono_type), DIMENSION(:), &
72 : POINTER :: charge_typ => NULL()
73 : END TYPE multi_charge_type
74 :
75 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .FALSE.
76 : LOGICAL, PRIVATE, PARAMETER :: debug_r_space = .FALSE.
77 : LOGICAL, PRIVATE, PARAMETER :: debug_g_space = .FALSE.
78 : LOGICAL, PRIVATE, PARAMETER :: debug_e_field = .FALSE.
79 : LOGICAL, PRIVATE, PARAMETER :: debug_e_field_en = .FALSE.
80 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewalds_multipole'
81 :
82 : PUBLIC :: ewald_multipole_evaluate
83 :
84 : CONTAINS
85 :
86 : ! **************************************************************************************************
87 : !> \brief Computes the potential and the force for a lattice sum of multipoles (up to quadrupole)
88 : !> \param ewald_env ...
89 : !> \param ewald_pw ...
90 : !> \param nonbond_env ...
91 : !> \param cell ...
92 : !> \param particle_set ...
93 : !> \param local_particles ...
94 : !> \param energy_local ...
95 : !> \param energy_glob ...
96 : !> \param e_neut ...
97 : !> \param e_self ...
98 : !> \param task ...
99 : !> \param do_correction_bonded ...
100 : !> \param do_forces ...
101 : !> \param do_stress ...
102 : !> \param do_efield ...
103 : !> \param radii ...
104 : !> \param charges ...
105 : !> \param dipoles ...
106 : !> \param quadrupoles ...
107 : !> \param forces_local ...
108 : !> \param forces_glob ...
109 : !> \param pv_local ...
110 : !> \param pv_glob ...
111 : !> \param efield0 ...
112 : !> \param efield1 ...
113 : !> \param efield2 ...
114 : !> \param iw ...
115 : !> \param do_debug ...
116 : !> \param atomic_kind_set ...
117 : !> \param mm_section ...
118 : !> \par Note
119 : !> atomic_kind_set and mm_section are between the arguments only
120 : !> for debug purpose (therefore optional) and can be avoided when this
121 : !> function is called in other part of the program
122 : !> \par Note
123 : !> When a gaussian multipole is used instead of point multipole, i.e.
124 : !> when radii(i)>0, the electrostatic fields (efield0, efield1, efield2)
125 : !> become derivatives of the electrostatic potential energy towards
126 : !> these gaussian multipoles.
127 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
128 : ! **************************************************************************************************
129 11208 : RECURSIVE SUBROUTINE ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, &
130 : cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, &
131 : task, do_correction_bonded, do_forces, do_stress, &
132 : do_efield, radii, charges, dipoles, &
133 7472 : quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1, &
134 3736 : efield2, iw, do_debug, atomic_kind_set, mm_section)
135 : TYPE(ewald_environment_type), POINTER :: ewald_env
136 : TYPE(ewald_pw_type), POINTER :: ewald_pw
137 : TYPE(fist_nonbond_env_type), POINTER :: nonbond_env
138 : TYPE(cell_type), POINTER :: cell
139 : TYPE(particle_type), POINTER :: particle_set(:)
140 : TYPE(distribution_1d_type), POINTER :: local_particles
141 : REAL(KIND=dp), INTENT(INOUT) :: energy_local, energy_glob
142 : REAL(KIND=dp), INTENT(OUT) :: e_neut, e_self
143 : LOGICAL, DIMENSION(3), INTENT(IN) :: task
144 : LOGICAL, INTENT(IN) :: do_correction_bonded, do_forces, &
145 : do_stress, do_efield
146 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges
147 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
148 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
149 : POINTER :: quadrupoles
150 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
151 : OPTIONAL :: forces_local, forces_glob, pv_local, &
152 : pv_glob
153 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT), OPTIONAL :: efield0
154 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
155 : OPTIONAL :: efield1, efield2
156 : INTEGER, INTENT(IN) :: iw
157 : LOGICAL, INTENT(IN) :: do_debug
158 : TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, &
159 : POINTER :: atomic_kind_set
160 : TYPE(section_vals_type), OPTIONAL, POINTER :: mm_section
161 :
162 : CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_evaluate'
163 :
164 : INTEGER :: handle, i, j, size1, size2
165 : LOGICAL :: check_debug, check_efield, check_forces, &
166 : do_task(3)
167 : LOGICAL, DIMENSION(3, 3) :: my_task
168 : REAL(KIND=dp) :: e_bonded, e_bonded_t, e_rspace, &
169 : e_rspace_t, energy_glob_t
170 3736 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0_lr, efield0_sr
171 3736 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1_lr, efield1_sr, efield2_lr, &
172 3736 : efield2_sr
173 : TYPE(mp_comm_type) :: group
174 :
175 3736 : CALL cite_reference(Aguado2003)
176 3736 : CALL cite_reference(Laino2008)
177 3736 : CALL timeset(routineN, handle)
178 3736 : CPASSERT(ASSOCIATED(nonbond_env))
179 : check_debug = (debug_this_module .OR. debug_r_space .OR. debug_g_space .OR. debug_e_field .OR. debug_e_field_en) &
180 3736 : .EQV. debug_this_module
181 : CPASSERT(check_debug)
182 3736 : check_forces = do_forces .EQV. (PRESENT(forces_local) .AND. PRESENT(forces_glob))
183 3736 : CPASSERT(check_forces)
184 3736 : check_efield = do_efield .EQV. (PRESENT(efield0) .OR. PRESENT(efield1) .OR. PRESENT(efield2))
185 3736 : CPASSERT(check_efield)
186 : ! Debugging this module
187 : IF (debug_this_module .AND. do_debug) THEN
188 : ! Debug specifically real space part
189 : IF (debug_r_space) THEN
190 : CALL debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, &
191 : particle_set, local_particles, iw, debug_r_space)
192 : CPABORT("Debug Multipole Requested: Real Part!")
193 : END IF
194 : ! Debug electric fields and gradients as pure derivatives
195 : IF (debug_e_field) THEN
196 : CPASSERT(PRESENT(atomic_kind_set))
197 : CPASSERT(PRESENT(mm_section))
198 : CALL debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, &
199 : cell, particle_set, local_particles, radii, charges, dipoles, &
200 : quadrupoles, task, iw, atomic_kind_set, mm_section)
201 : CPABORT("Debug Multipole Requested: POT+EFIELDS+GRAD!")
202 : END IF
203 : ! Debug the potential, electric fields and electric fields gradient in oder
204 : ! to retrieve the correct energy
205 : IF (debug_e_field_en) THEN
206 : CALL debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, &
207 : cell, particle_set, local_particles, radii, charges, dipoles, &
208 : quadrupoles, task, iw)
209 : CPABORT("Debug Multipole Requested: POT+EFIELDS+GRAD to give the correct energy!")
210 : END IF
211 : END IF
212 :
213 : ! Setup the tasks (needed to skip useless parts in the real-space term)
214 3736 : do_task = task
215 14944 : DO i = 1, 3
216 14944 : IF (do_task(i)) THEN
217 3270 : SELECT CASE (i)
218 : CASE (1)
219 5364 : do_task(1) = ANY(charges /= 0.0_dp)
220 : CASE (2)
221 128446 : do_task(2) = ANY(dipoles /= 0.0_dp)
222 : CASE (3)
223 39072 : do_task(3) = ANY(quadrupoles /= 0.0_dp)
224 : END SELECT
225 : END IF
226 : END DO
227 14944 : DO i = 1, 3
228 37360 : DO j = i, 3
229 22416 : my_task(j, i) = do_task(i) .AND. do_task(j)
230 33624 : my_task(i, j) = my_task(j, i)
231 : END DO
232 : END DO
233 :
234 : ! Allocate arrays for the evaluation of the potential, fields and electrostatic field gradients
235 3736 : NULLIFY (efield0_sr, efield0_lr, efield1_sr, efield1_lr, efield2_sr, efield2_lr)
236 3736 : IF (do_efield) THEN
237 2578 : IF (PRESENT(efield0)) THEN
238 1840 : size1 = SIZE(efield0)
239 5520 : ALLOCATE (efield0_sr(size1))
240 3680 : ALLOCATE (efield0_lr(size1))
241 18304 : efield0_sr = 0.0_dp
242 18304 : efield0_lr = 0.0_dp
243 : END IF
244 2578 : IF (PRESENT(efield1)) THEN
245 2578 : size1 = SIZE(efield1, 1)
246 2578 : size2 = SIZE(efield1, 2)
247 10312 : ALLOCATE (efield1_sr(size1, size2))
248 7734 : ALLOCATE (efield1_lr(size1, size2))
249 657034 : efield1_sr = 0.0_dp
250 657034 : efield1_lr = 0.0_dp
251 : END IF
252 2578 : IF (PRESENT(efield2)) THEN
253 2134 : size1 = SIZE(efield2, 1)
254 2134 : size2 = SIZE(efield2, 2)
255 8536 : ALLOCATE (efield2_sr(size1, size2))
256 6402 : ALLOCATE (efield2_lr(size1, size2))
257 1101534 : efield2_sr = 0.0_dp
258 1101534 : efield2_lr = 0.0_dp
259 : END IF
260 : END IF
261 :
262 3736 : e_rspace = 0.0_dp
263 3736 : e_bonded = 0.0_dp
264 3736 : IF ((.NOT. debug_g_space) .AND. (nonbond_env%do_nonbonded)) THEN
265 : ! Compute the Real Space (Short-Range) part of the Ewald sum.
266 : ! This contribution is only added when the nonbonded flag in the input
267 : ! is set, because these contributions depend. the neighborlists.
268 : CALL ewald_multipole_SR(nonbond_env, ewald_env, atomic_kind_set, &
269 : particle_set, cell, e_rspace, my_task, &
270 : do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, &
271 8852 : forces_glob, pv_glob, efield0_sr, efield1_sr, efield2_sr)
272 3736 : energy_glob = energy_glob + e_rspace
273 :
274 3736 : IF (do_correction_bonded) THEN
275 : ! The corrections for bonded interactions are stored in the Real Space
276 : ! (Short-Range) part of the fields array.
277 : CALL ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, &
278 : cell, e_bonded, my_task, do_forces, do_efield, do_stress, &
279 : charges, dipoles, quadrupoles, forces_glob, pv_glob, &
280 3372 : efield0_sr, efield1_sr, efield2_sr)
281 1896 : energy_glob = energy_glob + e_bonded
282 : END IF
283 : END IF
284 :
285 3736 : e_neut = 0.0_dp
286 3736 : e_self = 0.0_dp
287 3736 : energy_local = 0.0_dp
288 : IF (.NOT. debug_r_space) THEN
289 : ! Compute the Reciprocal Space (Long-Range) part of the Ewald sum
290 : CALL ewald_multipole_LR(ewald_env, ewald_pw, cell, particle_set, &
291 : local_particles, energy_local, my_task, do_forces, do_efield, do_stress, &
292 : charges, dipoles, quadrupoles, forces_local, pv_local, efield0_lr, efield1_lr, &
293 8852 : efield2_lr)
294 :
295 : ! Self-Interactions corrections
296 : CALL ewald_multipole_self(ewald_env, cell, local_particles, e_self, &
297 : e_neut, my_task, do_efield, radii, charges, dipoles, quadrupoles, &
298 3736 : efield0_lr, efield1_lr, efield2_lr)
299 : END IF
300 :
301 : ! Sumup energy contributions for possible IO
302 3736 : CALL ewald_env_get(ewald_env, group=group)
303 3736 : energy_glob_t = energy_glob
304 3736 : e_rspace_t = e_rspace
305 3736 : e_bonded_t = e_bonded
306 3736 : CALL group%sum(energy_glob_t)
307 3736 : CALL group%sum(e_rspace_t)
308 3736 : CALL group%sum(e_bonded_t)
309 : ! Print some info about energetics
310 3736 : CALL ewald_multipole_print(iw, energy_local, e_rspace_t, e_bonded_t, e_self, e_neut)
311 :
312 : ! Gather the components of the potential, fields and electrostatic field gradients
313 3736 : IF (do_efield) THEN
314 2578 : IF (PRESENT(efield0)) THEN
315 18304 : efield0 = efield0_sr + efield0_lr
316 34768 : CALL group%sum(efield0)
317 1840 : DEALLOCATE (efield0_sr)
318 1840 : DEALLOCATE (efield0_lr)
319 : END IF
320 2578 : IF (PRESENT(efield1)) THEN
321 657034 : efield1 = efield1_sr + efield1_lr
322 1311490 : CALL group%sum(efield1)
323 2578 : DEALLOCATE (efield1_sr)
324 2578 : DEALLOCATE (efield1_lr)
325 : END IF
326 2578 : IF (PRESENT(efield2)) THEN
327 1101534 : efield2 = efield2_sr + efield2_lr
328 2200934 : CALL group%sum(efield2)
329 2134 : DEALLOCATE (efield2_sr)
330 2134 : DEALLOCATE (efield2_lr)
331 : END IF
332 : END IF
333 3736 : CALL timestop(handle)
334 3736 : END SUBROUTINE ewald_multipole_evaluate
335 :
336 : ! **************************************************************************************************
337 : !> \brief computes the potential and the force for a lattice sum of multipoles
338 : !> up to quadrupole - Short Range (Real Space) Term
339 : !> \param nonbond_env ...
340 : !> \param ewald_env ...
341 : !> \param atomic_kind_set ...
342 : !> \param particle_set ...
343 : !> \param cell ...
344 : !> \param energy ...
345 : !> \param task ...
346 : !> \param do_forces ...
347 : !> \param do_efield ...
348 : !> \param do_stress ...
349 : !> \param radii ...
350 : !> \param charges ...
351 : !> \param dipoles ...
352 : !> \param quadrupoles ...
353 : !> \param forces ...
354 : !> \param pv ...
355 : !> \param efield0 ...
356 : !> \param efield1 ...
357 : !> \param efield2 ...
358 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
359 : ! **************************************************************************************************
360 7472 : SUBROUTINE ewald_multipole_SR(nonbond_env, ewald_env, atomic_kind_set, &
361 : particle_set, cell, energy, task, &
362 : do_forces, do_efield, do_stress, radii, charges, dipoles, quadrupoles, &
363 3736 : forces, pv, efield0, efield1, efield2)
364 : TYPE(fist_nonbond_env_type), POINTER :: nonbond_env
365 : TYPE(ewald_environment_type), POINTER :: ewald_env
366 : TYPE(atomic_kind_type), DIMENSION(:), OPTIONAL, &
367 : POINTER :: atomic_kind_set
368 : TYPE(particle_type), POINTER :: particle_set(:)
369 : TYPE(cell_type), POINTER :: cell
370 : REAL(KIND=dp), INTENT(INOUT) :: energy
371 : LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task
372 : LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress
373 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges
374 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
375 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
376 : POINTER :: quadrupoles
377 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
378 : OPTIONAL :: forces, pv
379 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0
380 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2
381 :
382 : CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_SR'
383 :
384 : INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, i, iend, igrp, ikind, ilist, ipair, &
385 : istart, itype_ij, itype_ji, jkind, k, kind_a, kind_b, kk, nkdamp_ij, nkdamp_ji, nkinds, &
386 : npairs
387 3736 : INTEGER, DIMENSION(:, :), POINTER :: list
388 : LOGICAL :: do_efield0, do_efield1, do_efield2, &
389 : force_eval
390 : REAL(KIND=dp) :: alpha, beta, ch_i, ch_j, dampa_ij, dampa_ji, dampaexpi, dampaexpj, &
391 : dampfac_ij, dampfac_ji, dampfuncdiffi, dampfuncdiffj, dampfunci, dampfuncj, dampsumfi, &
392 : dampsumfj, ef0_i, ef0_j, eloc, fac, fac_ij, factorial, ir, irab2, ptens11, ptens12, &
393 : ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, rab2_max, radius, &
394 : rcut, tij, tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, &
395 : tmp33, tmp_ij, tmp_ji, xf
396 : REAL(KIND=dp), DIMENSION(0:5) :: f
397 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi, damptij_a, damptji_a, dp_i, &
398 : dp_j, ef1_i, ef1_j, fr, rab, tij_a
399 : REAL(KIND=dp), DIMENSION(3, 3) :: damptij_ab, damptji_ab, ef2_i, ef2_j, &
400 : qp_i, qp_j, tij_ab
401 : REAL(KIND=dp), DIMENSION(3, 3, 3) :: tij_abc
402 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd
403 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
404 : TYPE(damping_type) :: damping_ij, damping_ji
405 : TYPE(fist_neighbor_type), POINTER :: nonbonded
406 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
407 3736 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc
408 :
409 3736 : CALL timeset(routineN, handle)
410 3736 : NULLIFY (nonbonded, r_last_update, r_last_update_pbc)
411 3736 : do_efield0 = do_efield .AND. ASSOCIATED(efield0)
412 3736 : do_efield1 = do_efield .AND. ASSOCIATED(efield1)
413 3736 : do_efield2 = do_efield .AND. ASSOCIATED(efield2)
414 : IF (do_stress) THEN
415 : ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
416 : ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
417 : ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
418 : END IF
419 : ! Get nonbond_env info
420 : CALL fist_nonbond_env_get(nonbond_env, nonbonded=nonbonded, natom_types=nkinds, &
421 3736 : r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc)
422 3736 : CALL ewald_env_get(ewald_env, alpha=alpha, rcut=rcut)
423 3736 : rab2_max = rcut**2
424 : IF (debug_r_space) THEN
425 : rab2_max = HUGE(0.0_dp)
426 : END IF
427 : ! Starting the force loop
428 5278820 : Lists: DO ilist = 1, nonbonded%nlists
429 5275084 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
430 5275084 : npairs = neighbor_kind_pair%npairs
431 5275084 : IF (npairs == 0) CYCLE Lists
432 1842978 : list => neighbor_kind_pair%list
433 7371912 : cvi = neighbor_kind_pair%cell_vector
434 23958714 : cell_v = MATMUL(cell%hmat, cvi)
435 4961652 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
436 3114938 : istart = neighbor_kind_pair%grp_kind_start(igrp)
437 3114938 : iend = neighbor_kind_pair%grp_kind_end(igrp)
438 3114938 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
439 3114938 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
440 :
441 3114938 : itype_ij = no_damping
442 3114938 : nkdamp_ij = 1
443 3114938 : dampa_ij = 1.0_dp
444 3114938 : dampfac_ij = 0.0_dp
445 :
446 3114938 : itype_ji = no_damping
447 3114938 : nkdamp_ji = 1
448 3114938 : dampa_ji = 1.0_dp
449 3114938 : dampfac_ji = 0.0_dp
450 3114938 : IF (PRESENT(atomic_kind_set)) THEN
451 3062577 : IF (ASSOCIATED(atomic_kind_set(jkind)%damping)) THEN
452 22135 : damping_ij = atomic_kind_set(jkind)%damping%damp(ikind)
453 22135 : itype_ij = damping_ij%itype
454 22135 : nkdamp_ij = damping_ij%order
455 22135 : dampa_ij = damping_ij%bij
456 22135 : dampfac_ij = damping_ij%cij
457 : END IF
458 :
459 3062577 : IF (ASSOCIATED(atomic_kind_set(ikind)%damping)) THEN
460 13035 : damping_ji = atomic_kind_set(ikind)%damping%damp(jkind)
461 13035 : itype_ji = damping_ji%itype
462 13035 : nkdamp_ji = damping_ji%order
463 13035 : dampa_ji = damping_ji%bij
464 13035 : dampfac_ji = damping_ji%cij
465 : END IF
466 : END IF
467 :
468 568112844 : Pairs: DO ipair = istart, iend
469 559722822 : IF (ipair <= neighbor_kind_pair%nscale) THEN
470 : ! scale the electrostatic interaction if needed
471 : ! (most often scaled to zero)
472 97950 : fac_ij = neighbor_kind_pair%ei_scale(ipair)
473 97950 : IF (fac_ij <= 0) CYCLE Pairs
474 : ELSE
475 : fac_ij = 1.0_dp
476 : END IF
477 559624872 : atom_a = list(1, ipair)
478 559624872 : atom_b = list(2, ipair)
479 559624872 : kind_a = particle_set(atom_a)%atomic_kind%kind_number
480 559624872 : kind_b = particle_set(atom_b)%atomic_kind%kind_number
481 559624872 : IF (atom_a == atom_b) fac_ij = 0.5_dp
482 2238499488 : rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
483 2238499488 : rab = rab + cell_v
484 559624872 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
485 562739810 : IF (rab2 <= rab2_max) THEN
486 21715939 : IF (PRESENT(radii)) THEN
487 21380381 : radius = SQRT(radii(atom_a)*radii(atom_a) + radii(atom_b)*radii(atom_b))
488 : ELSE
489 : radius = 0.0_dp
490 : END IF
491 21715939 : IF (radius > 0.0_dp) THEN
492 11 : beta = sqrthalf/radius
493 6727 : $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_GAUSS", damping=False, store_energy=True, store_forces=True)
494 : ELSE
495 14817322789 : $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERFC", damping=True, store_energy=True, store_forces=True )
496 : END IF
497 : END IF
498 : END DO Pairs
499 : END DO Kind_Group_Loop
500 : END DO Lists
501 3736 : IF (do_stress) THEN
502 38 : pv(1, 1) = pv(1, 1) + ptens11
503 38 : pv(1, 2) = pv(1, 2) + (ptens12 + ptens21)*0.5_dp
504 38 : pv(1, 3) = pv(1, 3) + (ptens13 + ptens31)*0.5_dp
505 38 : pv(2, 1) = pv(1, 2)
506 38 : pv(2, 2) = pv(2, 2) + ptens22
507 38 : pv(2, 3) = pv(2, 3) + (ptens23 + ptens32)*0.5_dp
508 38 : pv(3, 1) = pv(1, 3)
509 38 : pv(3, 2) = pv(2, 3)
510 38 : pv(3, 3) = pv(3, 3) + ptens33
511 : END IF
512 :
513 3736 : CALL timestop(handle)
514 3736 : END SUBROUTINE ewald_multipole_SR
515 :
516 : ! **************************************************************************************************
517 : !> \brief computes the bonded correction for the potential and the force for a
518 : !> lattice sum of multipoles up to quadrupole
519 : !> \param nonbond_env ...
520 : !> \param particle_set ...
521 : !> \param ewald_env ...
522 : !> \param cell ...
523 : !> \param energy ...
524 : !> \param task ...
525 : !> \param do_forces ...
526 : !> \param do_efield ...
527 : !> \param do_stress ...
528 : !> \param charges ...
529 : !> \param dipoles ...
530 : !> \param quadrupoles ...
531 : !> \param forces ...
532 : !> \param pv ...
533 : !> \param efield0 ...
534 : !> \param efield1 ...
535 : !> \param efield2 ...
536 : !> \author Teodoro Laino [tlaino] - 05.2009
537 : ! **************************************************************************************************
538 3792 : SUBROUTINE ewald_multipole_bonded(nonbond_env, particle_set, ewald_env, &
539 : cell, energy, task, do_forces, do_efield, do_stress, charges, &
540 1896 : dipoles, quadrupoles, forces, pv, efield0, efield1, efield2)
541 :
542 : TYPE(fist_nonbond_env_type), POINTER :: nonbond_env
543 : TYPE(particle_type), POINTER :: particle_set(:)
544 : TYPE(ewald_environment_type), POINTER :: ewald_env
545 : TYPE(cell_type), POINTER :: cell
546 : REAL(KIND=dp), INTENT(INOUT) :: energy
547 : LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task
548 : LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress
549 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
550 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
551 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
552 : POINTER :: quadrupoles
553 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
554 : OPTIONAL :: forces, pv
555 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0
556 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2
557 :
558 : CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_bonded'
559 :
560 : INTEGER :: a, atom_a, atom_b, b, c, d, e, handle, &
561 : i, iend, igrp, ilist, ipair, istart, &
562 : k, nscale
563 1896 : INTEGER, DIMENSION(:, :), POINTER :: list
564 : LOGICAL :: do_efield0, do_efield1, do_efield2, &
565 : force_eval
566 : REAL(KIND=dp) :: alpha, ch_i, ch_j, ef0_i, ef0_j, eloc, fac, fac_ij, ir, irab2, ptens11, &
567 : ptens12, ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33, r, rab2, tij, &
568 : tmp, tmp1, tmp11, tmp12, tmp13, tmp2, tmp21, tmp22, tmp23, tmp31, tmp32, tmp33, tmp_ij, &
569 : tmp_ji
570 : REAL(KIND=dp), DIMENSION(0:5) :: f
571 : REAL(KIND=dp), DIMENSION(3) :: dp_i, dp_j, ef1_i, ef1_j, fr, rab, tij_a
572 : REAL(KIND=dp), DIMENSION(3, 3) :: ef2_i, ef2_j, qp_i, qp_j, tij_ab
573 : REAL(KIND=dp), DIMENSION(3, 3, 3) :: tij_abc
574 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3) :: tij_abcd
575 : REAL(KIND=dp), DIMENSION(3, 3, 3, 3, 3) :: tij_abcde
576 : TYPE(fist_neighbor_type), POINTER :: nonbonded
577 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
578 :
579 1896 : CALL timeset(routineN, handle)
580 1896 : do_efield0 = do_efield .AND. ASSOCIATED(efield0)
581 1896 : do_efield1 = do_efield .AND. ASSOCIATED(efield1)
582 1896 : do_efield2 = do_efield .AND. ASSOCIATED(efield2)
583 : IF (do_stress) THEN
584 : ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
585 : ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
586 : ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
587 : END IF
588 1896 : CALL ewald_env_get(ewald_env, alpha=alpha)
589 1896 : CALL fist_nonbond_env_get(nonbond_env, nonbonded=nonbonded)
590 :
591 : ! Starting the force loop
592 5152080 : Lists: DO ilist = 1, nonbonded%nlists
593 5150184 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
594 5150184 : nscale = neighbor_kind_pair%nscale
595 5150184 : IF (nscale == 0) CYCLE Lists
596 1157 : list => neighbor_kind_pair%list
597 59169 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
598 56116 : istart = neighbor_kind_pair%grp_kind_start(igrp)
599 56116 : IF (istart > nscale) CYCLE Kind_Group_Loop
600 50773 : iend = MIN(neighbor_kind_pair%grp_kind_end(igrp), nscale)
601 5298907 : Pairs: DO ipair = istart, iend
602 : ! only use pairs that are (partially) excluded for electrostatics
603 97950 : fac_ij = -1.0_dp + neighbor_kind_pair%ei_scale(ipair)
604 97950 : IF (fac_ij >= 0) CYCLE Pairs
605 :
606 97950 : atom_a = list(1, ipair)
607 97950 : atom_b = list(2, ipair)
608 :
609 391800 : rab = particle_set(atom_b)%r - particle_set(atom_a)%r
610 391800 : rab = pbc(rab, cell)
611 97950 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
612 59621194 : $: ewalds_multipole_sr_macro(mode="SCREENED_COULOMB_ERF", store_energy=True, store_forces=True)
613 : END DO Pairs
614 : END DO Kind_Group_Loop
615 : END DO Lists
616 1896 : IF (do_stress) THEN
617 36 : pv(1, 1) = pv(1, 1) + ptens11
618 36 : pv(1, 2) = pv(1, 2) + (ptens12 + ptens21)*0.5_dp
619 36 : pv(1, 3) = pv(1, 3) + (ptens13 + ptens31)*0.5_dp
620 36 : pv(2, 1) = pv(1, 2)
621 36 : pv(2, 2) = pv(2, 2) + ptens22
622 36 : pv(2, 3) = pv(2, 3) + (ptens23 + ptens32)*0.5_dp
623 36 : pv(3, 1) = pv(1, 3)
624 36 : pv(3, 2) = pv(2, 3)
625 36 : pv(3, 3) = pv(3, 3) + ptens33
626 : END IF
627 :
628 1896 : CALL timestop(handle)
629 1896 : END SUBROUTINE ewald_multipole_bonded
630 :
631 : ! **************************************************************************************************
632 : !> \brief computes the potential and the force for a lattice sum of multipoles
633 : !> up to quadrupole - Long Range (Reciprocal Space) Term
634 : !> \param ewald_env ...
635 : !> \param ewald_pw ...
636 : !> \param cell ...
637 : !> \param particle_set ...
638 : !> \param local_particles ...
639 : !> \param energy ...
640 : !> \param task ...
641 : !> \param do_forces ...
642 : !> \param do_efield ...
643 : !> \param do_stress ...
644 : !> \param charges ...
645 : !> \param dipoles ...
646 : !> \param quadrupoles ...
647 : !> \param forces ...
648 : !> \param pv ...
649 : !> \param efield0 ...
650 : !> \param efield1 ...
651 : !> \param efield2 ...
652 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
653 : ! **************************************************************************************************
654 3736 : SUBROUTINE ewald_multipole_LR(ewald_env, ewald_pw, cell, particle_set, &
655 : local_particles, energy, task, do_forces, do_efield, do_stress, &
656 3736 : charges, dipoles, quadrupoles, forces, pv, efield0, efield1, efield2)
657 : TYPE(ewald_environment_type), POINTER :: ewald_env
658 : TYPE(ewald_pw_type), POINTER :: ewald_pw
659 : TYPE(cell_type), POINTER :: cell
660 : TYPE(particle_type), POINTER :: particle_set(:)
661 : TYPE(distribution_1d_type), POINTER :: local_particles
662 : REAL(KIND=dp), INTENT(INOUT) :: energy
663 : LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task
664 : LOGICAL, INTENT(IN) :: do_forces, do_efield, do_stress
665 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
666 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
667 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
668 : POINTER :: quadrupoles
669 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
670 : OPTIONAL :: forces, pv
671 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0
672 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2
673 :
674 : CHARACTER(len=*), PARAMETER :: routineN = 'ewald_multipole_LR'
675 :
676 : COMPLEX(KIND=dp) :: atm_factor, atm_factor_st(3), cnjg_fac, &
677 : fac, summe_tmp
678 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: summe_ef
679 3736 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: summe_st
680 : INTEGER :: gpt, handle, iparticle, iparticle_kind, iparticle_local, lp, mp, nnodes, &
681 : node, np, nparticle_kind, nparticle_local
682 3736 : INTEGER, DIMENSION(:, :), POINTER :: bds
683 : LOGICAL :: do_efield0, do_efield1, do_efield2
684 : REAL(KIND=dp) :: alpha, denom, dipole_t(3), f0, factor, &
685 : four_alpha_sq, gauss, pref, q_t, tmp, &
686 : trq_t
687 : REAL(KIND=dp), DIMENSION(3) :: tmp_v, vec
688 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_tmp
689 3736 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0
690 : TYPE(dg_rho0_type), POINTER :: dg_rho0
691 : TYPE(dg_type), POINTER :: dg
692 : TYPE(pw_grid_type), POINTER :: pw_grid
693 : TYPE(pw_pool_type), POINTER :: pw_pool
694 : TYPE(structure_factor_type) :: exp_igr
695 : TYPE(mp_comm_type) :: group
696 :
697 3736 : CALL timeset(routineN, handle)
698 3736 : do_efield0 = do_efield .AND. ASSOCIATED(efield0)
699 3736 : do_efield1 = do_efield .AND. ASSOCIATED(efield1)
700 3736 : do_efield2 = do_efield .AND. ASSOCIATED(efield2)
701 :
702 : ! Gathering data from the ewald environment
703 3736 : CALL ewald_env_get(ewald_env, alpha=alpha, group=group)
704 3736 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, dg=dg)
705 3736 : CALL dg_get(dg, dg_rho0=dg_rho0)
706 3736 : rho0 => dg_rho0%density%array
707 3736 : pw_grid => pw_pool%pw_grid
708 3736 : bds => pw_grid%bounds
709 :
710 : ! Allocation of working arrays
711 3736 : nparticle_kind = SIZE(local_particles%n_el)
712 3736 : nnodes = 0
713 11444 : DO iparticle_kind = 1, nparticle_kind
714 11444 : nnodes = nnodes + local_particles%n_el(iparticle_kind)
715 : END DO
716 3736 : CALL structure_factor_allocate(pw_grid%bounds, nnodes, exp_igr)
717 :
718 11208 : ALLOCATE (summe_ef(1:pw_grid%ngpts_cut))
719 3736 : summe_ef = z_zero
720 : ! Stress Tensor
721 3736 : IF (do_stress) THEN
722 38 : pv_tmp = 0.0_dp
723 114 : ALLOCATE (summe_st(3, 1:pw_grid%ngpts_cut))
724 38 : summe_st = z_zero
725 : END IF
726 :
727 : ! Defining four_alpha_sq
728 3736 : four_alpha_sq = 4.0_dp*alpha**2
729 3736 : dipole_t = 0.0_dp
730 3736 : q_t = 0.0_dp
731 3736 : trq_t = 0.0_dp
732 : ! Zero node count
733 3736 : node = 0
734 11444 : DO iparticle_kind = 1, nparticle_kind
735 7708 : nparticle_local = local_particles%n_el(iparticle_kind)
736 102983 : DO iparticle_local = 1, nparticle_local
737 91539 : node = node + 1
738 91539 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
739 1190007 : vec = MATMUL(cell%h_inv, particle_set(iparticle)%r)
740 : CALL structure_factor_evaluate(vec, exp_igr%lb, &
741 91539 : exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
742 :
743 : ! Computing the total charge, dipole and quadrupole trace (if any)
744 100704 : IF (ANY(task(1, :))) THEN
745 88484 : q_t = q_t + charges(iparticle)
746 : END IF
747 123835 : IF (ANY(task(2, :))) THEN
748 323728 : dipole_t = dipole_t + dipoles(:, iparticle)
749 : END IF
750 354682 : IF (ANY(task(3, :))) THEN
751 : trq_t = trq_t + quadrupoles(1, 1, iparticle) + &
752 : quadrupoles(2, 2, iparticle) + &
753 6627 : quadrupoles(3, 3, iparticle)
754 : END IF
755 : END DO
756 : END DO
757 :
758 : ! Looping over the positive g-vectors
759 143536764 : DO gpt = 1, pw_grid%ngpts_cut_local
760 143533028 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
761 143533028 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
762 143533028 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
763 :
764 143533028 : lp = lp + bds(1, 1)
765 143533028 : mp = mp + bds(1, 2)
766 143533028 : np = np + bds(1, 3)
767 :
768 : ! Initializing sum to be used in the energy and force
769 143533028 : node = 0
770 428239136 : DO iparticle_kind = 1, nparticle_kind
771 284702372 : nparticle_local = local_particles%n_el(iparticle_kind)
772 909275529 : DO iparticle_local = 1, nparticle_local
773 481040129 : node = node + 1
774 481040129 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
775 : ! Density for energy and forces
776 : CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
777 481040129 : dipoles, quadrupoles)
778 481040129 : summe_tmp = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node)
779 481040129 : summe_ef(gpt) = summe_ef(gpt) + atm_factor*summe_tmp
780 :
781 : ! Precompute pseudo-density for stress tensor calculation
782 765742501 : IF (do_stress) THEN
783 : CALL get_atom_factor_stress(atm_factor_st, pw_grid, gpt, iparticle, task, &
784 8802187 : dipoles, quadrupoles)
785 35208748 : summe_st(1:3, gpt) = summe_st(1:3, gpt) + atm_factor_st(1:3)*summe_tmp
786 : END IF
787 : END DO
788 : END DO
789 : END DO
790 3736 : CALL group%sum(q_t)
791 3736 : CALL group%sum(dipole_t)
792 3736 : CALL group%sum(trq_t)
793 3736 : CALL group%sum(summe_ef)
794 3736 : IF (do_stress) CALL group%sum(summe_st)
795 :
796 : ! Looping over the positive g-vectors
797 143536764 : DO gpt = 1, pw_grid%ngpts_cut_local
798 : ! computing the potential energy
799 143533028 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
800 143533028 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
801 143533028 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
802 :
803 143533028 : lp = lp + bds(1, 1)
804 143533028 : mp = mp + bds(1, 2)
805 143533028 : np = np + bds(1, 3)
806 :
807 143533028 : IF (pw_grid%gsq(gpt) == 0.0_dp) THEN
808 : ! G=0 vector for dipole-dipole and charge-quadrupole
809 : energy = energy + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t) &
810 14944 : - (1.0_dp/9.0_dp)*q_t*trq_t
811 : ! Stress tensor
812 3736 : IF (do_stress) THEN
813 152 : pv_tmp(1, 1) = pv_tmp(1, 1) + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t)
814 152 : pv_tmp(2, 2) = pv_tmp(2, 2) + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t)
815 152 : pv_tmp(3, 3) = pv_tmp(3, 3) + (1.0_dp/6.0_dp)*DOT_PRODUCT(dipole_t, dipole_t)
816 : END IF
817 : ! Corrections for G=0 to potential, field and field gradient
818 3736 : IF (do_efield .AND. (debug_e_field_en .OR. (.NOT. debug_this_module))) THEN
819 : ! This term is important and may give problems if one is debugging
820 : ! VS finite differences since it comes from a residual integral in
821 : ! the complex plane (cannot be reproduced with finite differences)
822 : node = 0
823 7920 : DO iparticle_kind = 1, nparticle_kind
824 5342 : nparticle_local = local_particles%n_el(iparticle_kind)
825 89727 : DO iparticle_local = 1, nparticle_local
826 81807 : node = node + 1
827 81807 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
828 :
829 : ! Potential
830 : IF (do_efield0) THEN
831 : efield0(iparticle) = efield0(iparticle)
832 : END IF
833 : ! Electrostatic field
834 81807 : IF (do_efield1) THEN
835 327228 : efield1(1:3, iparticle) = efield1(1:3, iparticle) - (1.0_dp/6.0_dp)*dipole_t(1:3)
836 : END IF
837 : ! Electrostatic field gradients
838 87149 : IF (do_efield2) THEN
839 54970 : efield2(1, iparticle) = efield2(1, iparticle) - (1.0_dp/(18.0_dp))*q_t
840 54970 : efield2(5, iparticle) = efield2(5, iparticle) - (1.0_dp/(18.0_dp))*q_t
841 54970 : efield2(9, iparticle) = efield2(9, iparticle) - (1.0_dp/(18.0_dp))*q_t
842 : END IF
843 : END DO
844 : END DO
845 : END IF
846 : CYCLE
847 : END IF
848 143529292 : gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt)
849 143529292 : factor = gauss*REAL(summe_ef(gpt)*CONJG(summe_ef(gpt)), KIND=dp)
850 143529292 : energy = energy + factor
851 :
852 143529292 : IF (do_forces .OR. do_efield) THEN
853 : node = 0
854 428223956 : DO iparticle_kind = 1, nparticle_kind
855 284694664 : nparticle_local = local_particles%n_el(iparticle_kind)
856 909172546 : DO iparticle_local = 1, nparticle_local
857 480948590 : node = node + 1
858 480948590 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
859 480948590 : fac = exp_igr%ex(lp, node)*exp_igr%ey(mp, node)*exp_igr%ez(np, node)
860 480948590 : cnjg_fac = CONJG(fac)
861 :
862 : ! Forces
863 480948590 : IF (do_forces) THEN
864 : CALL get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
865 222572226 : dipoles, quadrupoles)
866 :
867 222572226 : tmp = gauss*AIMAG(summe_ef(gpt)*(cnjg_fac*CONJG(atm_factor)))
868 222572226 : forces(1, node) = forces(1, node) + tmp*pw_grid%g(1, gpt)
869 222572226 : forces(2, node) = forces(2, node) + tmp*pw_grid%g(2, gpt)
870 222572226 : forces(3, node) = forces(3, node) + tmp*pw_grid%g(3, gpt)
871 : END IF
872 :
873 : ! Electric field
874 765643254 : IF (do_efield) THEN
875 : ! Potential
876 258810650 : IF (do_efield0) THEN
877 27316155 : efield0(iparticle) = efield0(iparticle) + gauss*REAL(fac*CONJG(summe_ef(gpt)), KIND=dp)
878 : END IF
879 : ! Electric field
880 258810650 : IF (do_efield1) THEN
881 258810650 : tmp = AIMAG(fac*CONJG(summe_ef(gpt)))*gauss
882 258810650 : efield1(1, iparticle) = efield1(1, iparticle) - tmp*pw_grid%g(1, gpt)
883 258810650 : efield1(2, iparticle) = efield1(2, iparticle) - tmp*pw_grid%g(2, gpt)
884 258810650 : efield1(3, iparticle) = efield1(3, iparticle) - tmp*pw_grid%g(3, gpt)
885 : END IF
886 : ! Electric field gradient
887 258810650 : IF (do_efield2) THEN
888 185990301 : tmp_v(1) = REAL(fac*CONJG(summe_ef(gpt)), KIND=dp)*pw_grid%g(1, gpt)*gauss
889 185990301 : tmp_v(2) = REAL(fac*CONJG(summe_ef(gpt)), KIND=dp)*pw_grid%g(2, gpt)*gauss
890 185990301 : tmp_v(3) = REAL(fac*CONJG(summe_ef(gpt)), KIND=dp)*pw_grid%g(3, gpt)*gauss
891 :
892 185990301 : efield2(1, iparticle) = efield2(1, iparticle) + tmp_v(1)*pw_grid%g(1, gpt)
893 185990301 : efield2(2, iparticle) = efield2(2, iparticle) + tmp_v(1)*pw_grid%g(2, gpt)
894 185990301 : efield2(3, iparticle) = efield2(3, iparticle) + tmp_v(1)*pw_grid%g(3, gpt)
895 185990301 : efield2(4, iparticle) = efield2(4, iparticle) + tmp_v(2)*pw_grid%g(1, gpt)
896 185990301 : efield2(5, iparticle) = efield2(5, iparticle) + tmp_v(2)*pw_grid%g(2, gpt)
897 185990301 : efield2(6, iparticle) = efield2(6, iparticle) + tmp_v(2)*pw_grid%g(3, gpt)
898 185990301 : efield2(7, iparticle) = efield2(7, iparticle) + tmp_v(3)*pw_grid%g(1, gpt)
899 185990301 : efield2(8, iparticle) = efield2(8, iparticle) + tmp_v(3)*pw_grid%g(2, gpt)
900 185990301 : efield2(9, iparticle) = efield2(9, iparticle) + tmp_v(3)*pw_grid%g(3, gpt)
901 : END IF
902 : END IF
903 : END DO
904 : END DO
905 : END IF
906 :
907 : ! Compute the virial P*V
908 143533028 : IF (do_stress) THEN
909 : ! The Stress Tensor can be decomposed in two main components.
910 : ! The first one is just a normal ewald component for reciprocal space
911 1841078 : denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
912 1841078 : pv_tmp(1, 1) = pv_tmp(1, 1) + factor*(1.0_dp - 2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(1, gpt)*denom)
913 1841078 : pv_tmp(1, 2) = pv_tmp(1, 2) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
914 1841078 : pv_tmp(1, 3) = pv_tmp(1, 3) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
915 1841078 : pv_tmp(2, 1) = pv_tmp(2, 1) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
916 1841078 : pv_tmp(2, 2) = pv_tmp(2, 2) + factor*(1.0_dp - 2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(2, gpt)*denom)
917 1841078 : pv_tmp(2, 3) = pv_tmp(2, 3) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
918 1841078 : pv_tmp(3, 1) = pv_tmp(3, 1) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
919 1841078 : pv_tmp(3, 2) = pv_tmp(3, 2) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
920 1841078 : pv_tmp(3, 3) = pv_tmp(3, 3) + factor*(1.0_dp - 2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(3, gpt)*denom)
921 : ! The second one can be written in the following way
922 1841078 : f0 = 2.0_dp*gauss
923 1841078 : pv_tmp(1, 1) = pv_tmp(1, 1) + f0*pw_grid%g(1, gpt)*REAL(summe_st(1, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
924 1841078 : pv_tmp(1, 2) = pv_tmp(1, 2) + f0*pw_grid%g(1, gpt)*REAL(summe_st(2, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
925 1841078 : pv_tmp(1, 3) = pv_tmp(1, 3) + f0*pw_grid%g(1, gpt)*REAL(summe_st(3, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
926 1841078 : pv_tmp(2, 1) = pv_tmp(2, 1) + f0*pw_grid%g(2, gpt)*REAL(summe_st(1, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
927 1841078 : pv_tmp(2, 2) = pv_tmp(2, 2) + f0*pw_grid%g(2, gpt)*REAL(summe_st(2, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
928 1841078 : pv_tmp(2, 3) = pv_tmp(2, 3) + f0*pw_grid%g(2, gpt)*REAL(summe_st(3, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
929 1841078 : pv_tmp(3, 1) = pv_tmp(3, 1) + f0*pw_grid%g(3, gpt)*REAL(summe_st(1, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
930 1841078 : pv_tmp(3, 2) = pv_tmp(3, 2) + f0*pw_grid%g(3, gpt)*REAL(summe_st(2, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
931 1841078 : pv_tmp(3, 3) = pv_tmp(3, 3) + f0*pw_grid%g(3, gpt)*REAL(summe_st(3, gpt)*CONJG(summe_ef(gpt)), KIND=dp)
932 : END IF
933 : END DO
934 3736 : pref = fourpi/pw_grid%vol
935 3736 : energy = energy*pref
936 :
937 3736 : CALL structure_factor_deallocate(exp_igr)
938 3736 : DEALLOCATE (summe_ef)
939 3736 : IF (do_stress) THEN
940 494 : pv_tmp = pv_tmp*pref
941 : ! Symmetrize the tensor
942 38 : pv(1, 1) = pv(1, 1) + pv_tmp(1, 1)
943 38 : pv(1, 2) = pv(1, 2) + (pv_tmp(1, 2) + pv_tmp(2, 1))*0.5_dp
944 38 : pv(1, 3) = pv(1, 3) + (pv_tmp(1, 3) + pv_tmp(3, 1))*0.5_dp
945 38 : pv(2, 1) = pv(1, 2)
946 38 : pv(2, 2) = pv(2, 2) + pv_tmp(2, 2)
947 38 : pv(2, 3) = pv(2, 3) + (pv_tmp(2, 3) + pv_tmp(3, 2))*0.5_dp
948 38 : pv(3, 1) = pv(1, 3)
949 38 : pv(3, 2) = pv(2, 3)
950 38 : pv(3, 3) = pv(3, 3) + pv_tmp(3, 3)
951 38 : DEALLOCATE (summe_st)
952 : END IF
953 3736 : IF (do_forces) THEN
954 40754 : forces = 2.0_dp*forces*pref
955 : END IF
956 3736 : IF (do_efield0) THEN
957 18304 : efield0 = 2.0_dp*efield0*pref
958 : END IF
959 3736 : IF (do_efield1) THEN
960 657034 : efield1 = 2.0_dp*efield1*pref
961 : END IF
962 3736 : IF (do_efield2) THEN
963 1101534 : efield2 = 2.0_dp*efield2*pref
964 : END IF
965 3736 : CALL timestop(handle)
966 :
967 22416 : END SUBROUTINE ewald_multipole_LR
968 :
969 : ! **************************************************************************************************
970 : !> \brief Computes the atom factor including charge, dipole and quadrupole
971 : !> \param atm_factor ...
972 : !> \param pw_grid ...
973 : !> \param gpt ...
974 : !> \param iparticle ...
975 : !> \param task ...
976 : !> \param charges ...
977 : !> \param dipoles ...
978 : !> \param quadrupoles ...
979 : !> \par History
980 : !> none
981 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
982 : ! **************************************************************************************************
983 703612355 : SUBROUTINE get_atom_factor(atm_factor, pw_grid, gpt, iparticle, task, charges, &
984 : dipoles, quadrupoles)
985 : COMPLEX(KIND=dp), INTENT(OUT) :: atm_factor
986 : TYPE(pw_grid_type), POINTER :: pw_grid
987 : INTEGER, INTENT(IN) :: gpt
988 : INTEGER :: iparticle
989 : LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task
990 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
991 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
992 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
993 : POINTER :: quadrupoles
994 :
995 : COMPLEX(KIND=dp) :: tmp
996 : INTEGER :: i, j
997 :
998 703612355 : atm_factor = z_zero
999 703612355 : IF (task(1, 1)) THEN
1000 : ! Charge
1001 485719082 : atm_factor = atm_factor + charges(iparticle)
1002 : END IF
1003 703612355 : IF (task(2, 2)) THEN
1004 : ! Dipole
1005 : tmp = z_zero
1006 2073876004 : DO i = 1, 3
1007 2073876004 : tmp = tmp + dipoles(i, iparticle)*pw_grid%g(i, gpt)
1008 : END DO
1009 518469001 : atm_factor = atm_factor + tmp*CMPLX(0.0_dp, -1.0_dp, KIND=dp)
1010 : END IF
1011 703612355 : IF (task(3, 3)) THEN
1012 : ! Quadrupole
1013 : tmp = z_zero
1014 939275996 : DO i = 1, 3
1015 3052646987 : DO j = 1, 3
1016 2817827988 : tmp = tmp + quadrupoles(j, i, iparticle)*pw_grid%g(j, gpt)*pw_grid%g(i, gpt)
1017 : END DO
1018 : END DO
1019 234818999 : atm_factor = atm_factor - 1.0_dp/3.0_dp*tmp
1020 : END IF
1021 :
1022 703612355 : END SUBROUTINE get_atom_factor
1023 :
1024 : ! **************************************************************************************************
1025 : !> \brief Computes the atom factor including charge, dipole and quadrupole
1026 : !> \param atm_factor ...
1027 : !> \param pw_grid ...
1028 : !> \param gpt ...
1029 : !> \param iparticle ...
1030 : !> \param task ...
1031 : !> \param dipoles ...
1032 : !> \param quadrupoles ...
1033 : !> \par History
1034 : !> none
1035 : !> \author Teodoro Laino [tlaino] - 12.2007 - University of Zurich
1036 : ! **************************************************************************************************
1037 8802187 : SUBROUTINE get_atom_factor_stress(atm_factor, pw_grid, gpt, iparticle, task, &
1038 : dipoles, quadrupoles)
1039 : COMPLEX(KIND=dp), INTENT(OUT) :: atm_factor(3)
1040 : TYPE(pw_grid_type), POINTER :: pw_grid
1041 : INTEGER, INTENT(IN) :: gpt
1042 : INTEGER :: iparticle
1043 : LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task
1044 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
1045 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
1046 : POINTER :: quadrupoles
1047 :
1048 : INTEGER :: i
1049 :
1050 8802187 : atm_factor = z_zero
1051 12459514 : IF (ANY(task(2, :))) THEN
1052 : ! Dipole
1053 31453744 : atm_factor = dipoles(:, iparticle)*CMPLX(0.0_dp, -1.0_dp, KIND=dp)
1054 : END IF
1055 32084455 : IF (ANY(task(3, :))) THEN
1056 : ! Quadrupole
1057 6049968 : DO i = 1, 3
1058 : atm_factor(1) = atm_factor(1) - 1.0_dp/3.0_dp* &
1059 : (quadrupoles(1, i, iparticle)*pw_grid%g(i, gpt) + &
1060 4537476 : quadrupoles(i, 1, iparticle)*pw_grid%g(i, gpt))
1061 : atm_factor(2) = atm_factor(2) - 1.0_dp/3.0_dp* &
1062 : (quadrupoles(2, i, iparticle)*pw_grid%g(i, gpt) + &
1063 4537476 : quadrupoles(i, 2, iparticle)*pw_grid%g(i, gpt))
1064 : atm_factor(3) = atm_factor(3) - 1.0_dp/3.0_dp* &
1065 : (quadrupoles(3, i, iparticle)*pw_grid%g(i, gpt) + &
1066 6049968 : quadrupoles(i, 3, iparticle)*pw_grid%g(i, gpt))
1067 : END DO
1068 : END IF
1069 :
1070 8802187 : END SUBROUTINE get_atom_factor_stress
1071 :
1072 : ! **************************************************************************************************
1073 : !> \brief Computes the self interaction from g-space and the neutralizing background
1074 : !> when using multipoles
1075 : !> \param ewald_env ...
1076 : !> \param cell ...
1077 : !> \param local_particles ...
1078 : !> \param e_self ...
1079 : !> \param e_neut ...
1080 : !> \param task ...
1081 : !> \param do_efield ...
1082 : !> \param radii ...
1083 : !> \param charges ...
1084 : !> \param dipoles ...
1085 : !> \param quadrupoles ...
1086 : !> \param efield0 ...
1087 : !> \param efield1 ...
1088 : !> \param efield2 ...
1089 : !> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007
1090 : ! **************************************************************************************************
1091 3736 : SUBROUTINE ewald_multipole_self(ewald_env, cell, local_particles, e_self, &
1092 : e_neut, task, do_efield, radii, charges, dipoles, quadrupoles, efield0, &
1093 : efield1, efield2)
1094 : TYPE(ewald_environment_type), POINTER :: ewald_env
1095 : TYPE(cell_type), POINTER :: cell
1096 : TYPE(distribution_1d_type), POINTER :: local_particles
1097 : REAL(KIND=dp), INTENT(OUT) :: e_self, e_neut
1098 : LOGICAL, DIMENSION(3, 3), INTENT(IN) :: task
1099 : LOGICAL, INTENT(IN) :: do_efield
1100 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges
1101 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
1102 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
1103 : POINTER :: quadrupoles
1104 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0
1105 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2
1106 :
1107 : REAL(KIND=dp), PARAMETER :: f23 = 2.0_dp/3.0_dp, &
1108 : f415 = 4.0_dp/15.0_dp
1109 :
1110 : INTEGER :: ewald_type, i, iparticle, &
1111 : iparticle_kind, iparticle_local, j, &
1112 : nparticle_local
1113 : LOGICAL :: do_efield0, do_efield1, do_efield2, &
1114 : lradii
1115 : REAL(KIND=dp) :: alpha, ch_qu_self, ch_qu_self_tmp, &
1116 : dipole_self, fac1, fac2, fac3, fac4, &
1117 : q, q_neutg, q_self, q_sum, qu_qu_self, &
1118 : radius
1119 : TYPE(mp_comm_type) :: group
1120 :
1121 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha, &
1122 3736 : group=group)
1123 :
1124 3736 : do_efield0 = do_efield .AND. ASSOCIATED(efield0)
1125 3736 : do_efield1 = do_efield .AND. ASSOCIATED(efield1)
1126 3736 : do_efield2 = do_efield .AND. ASSOCIATED(efield2)
1127 3736 : q_self = 0.0_dp
1128 3736 : q_sum = 0.0_dp
1129 3736 : dipole_self = 0.0_dp
1130 3736 : ch_qu_self = 0.0_dp
1131 3736 : qu_qu_self = 0.0_dp
1132 3736 : fac1 = 2.0_dp*alpha*oorootpi
1133 3736 : fac2 = 6.0_dp*(f23**2)*(alpha**3)*oorootpi
1134 3736 : fac3 = (2.0_dp*oorootpi)*f23*alpha**3
1135 3736 : fac4 = (4.0_dp*oorootpi)*f415*alpha**5
1136 3736 : lradii = PRESENT(radii)
1137 3736 : radius = 0.0_dp
1138 3736 : q_neutg = 0.0_dp
1139 11444 : DO iparticle_kind = 1, SIZE(local_particles%n_el)
1140 7708 : nparticle_local = local_particles%n_el(iparticle_kind)
1141 102983 : DO iparticle_local = 1, nparticle_local
1142 91539 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1143 100704 : IF (ANY(task(1, :))) THEN
1144 : ! Charge - Charge
1145 88484 : q = charges(iparticle)
1146 88484 : IF (lradii) radius = radii(iparticle)
1147 88484 : IF (radius > 0) THEN
1148 11 : q_neutg = q_neutg + 2.0_dp*q*radius**2
1149 : END IF
1150 88484 : q_self = q_self + q*q
1151 88484 : q_sum = q_sum + q
1152 : ! Potential
1153 88484 : IF (do_efield0) THEN
1154 5917 : efield0(iparticle) = efield0(iparticle) - q*fac1
1155 : END IF
1156 :
1157 88484 : IF (task(1, 3)) THEN
1158 : ! Charge - Quadrupole
1159 : ch_qu_self_tmp = 0.0_dp
1160 24644 : DO i = 1, 3
1161 24644 : ch_qu_self_tmp = ch_qu_self_tmp + quadrupoles(i, i, iparticle)*q
1162 : END DO
1163 6161 : ch_qu_self = ch_qu_self + ch_qu_self_tmp
1164 : ! Electric Field Gradient
1165 6161 : IF (do_efield2) THEN
1166 5811 : efield2(1, iparticle) = efield2(1, iparticle) + fac2*q
1167 5811 : efield2(5, iparticle) = efield2(5, iparticle) + fac2*q
1168 5811 : efield2(9, iparticle) = efield2(9, iparticle) + fac2*q
1169 : END IF
1170 : END IF
1171 : END IF
1172 123835 : IF (ANY(task(2, :))) THEN
1173 : ! Dipole - Dipole
1174 323728 : DO i = 1, 3
1175 323728 : dipole_self = dipole_self + dipoles(i, iparticle)**2
1176 : END DO
1177 : ! Electric Field
1178 80932 : IF (do_efield1) THEN
1179 72038 : efield1(1, iparticle) = efield1(1, iparticle) + fac3*dipoles(1, iparticle)
1180 72038 : efield1(2, iparticle) = efield1(2, iparticle) + fac3*dipoles(2, iparticle)
1181 72038 : efield1(3, iparticle) = efield1(3, iparticle) + fac3*dipoles(3, iparticle)
1182 : END IF
1183 : END IF
1184 354682 : IF (ANY(task(3, :))) THEN
1185 : ! Quadrupole - Quadrupole
1186 26508 : DO i = 1, 3
1187 86151 : DO j = 1, 3
1188 79524 : qu_qu_self = qu_qu_self + quadrupoles(j, i, iparticle)**2
1189 : END DO
1190 : END DO
1191 : ! Electric Field Gradient
1192 6627 : IF (do_efield2) THEN
1193 5811 : efield2(1, iparticle) = efield2(1, iparticle) + fac4*quadrupoles(1, 1, iparticle)
1194 5811 : efield2(2, iparticle) = efield2(2, iparticle) + fac4*quadrupoles(2, 1, iparticle)
1195 5811 : efield2(3, iparticle) = efield2(3, iparticle) + fac4*quadrupoles(3, 1, iparticle)
1196 5811 : efield2(4, iparticle) = efield2(4, iparticle) + fac4*quadrupoles(1, 2, iparticle)
1197 5811 : efield2(5, iparticle) = efield2(5, iparticle) + fac4*quadrupoles(2, 2, iparticle)
1198 5811 : efield2(6, iparticle) = efield2(6, iparticle) + fac4*quadrupoles(3, 2, iparticle)
1199 5811 : efield2(7, iparticle) = efield2(7, iparticle) + fac4*quadrupoles(1, 3, iparticle)
1200 5811 : efield2(8, iparticle) = efield2(8, iparticle) + fac4*quadrupoles(2, 3, iparticle)
1201 5811 : efield2(9, iparticle) = efield2(9, iparticle) + fac4*quadrupoles(3, 3, iparticle)
1202 : END IF
1203 : END IF
1204 : END DO
1205 : END DO
1206 :
1207 3736 : CALL group%sum(q_neutg)
1208 3736 : CALL group%sum(q_self)
1209 3736 : CALL group%sum(q_sum)
1210 3736 : CALL group%sum(dipole_self)
1211 3736 : CALL group%sum(ch_qu_self)
1212 3736 : CALL group%sum(qu_qu_self)
1213 :
1214 3736 : e_self = -(q_self + f23*(dipole_self - f23*ch_qu_self + f415*qu_qu_self*alpha**2)*alpha**2)*alpha*oorootpi
1215 3736 : fac1 = pi/(2.0_dp*cell%deth)
1216 3736 : e_neut = -q_sum*fac1*(q_sum/alpha**2 - q_neutg)
1217 :
1218 : ! Correcting Potential for the neutralizing background charge
1219 11444 : DO iparticle_kind = 1, SIZE(local_particles%n_el)
1220 7708 : nparticle_local = local_particles%n_el(iparticle_kind)
1221 102983 : DO iparticle_local = 1, nparticle_local
1222 91539 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
1223 108412 : IF (ANY(task(1, :))) THEN
1224 : ! Potential energy
1225 88484 : IF (do_efield0) THEN
1226 5917 : efield0(iparticle) = efield0(iparticle) - q_sum*2.0_dp*fac1/alpha**2
1227 5917 : IF (lradii) radius = radii(iparticle)
1228 5917 : IF (radius > 0) THEN
1229 0 : q = charges(iparticle)
1230 0 : efield0(iparticle) = efield0(iparticle) + fac1*radius**2*(q_sum + q)
1231 : END IF
1232 : END IF
1233 : END IF
1234 : END DO
1235 : END DO
1236 :
1237 3736 : END SUBROUTINE ewald_multipole_self
1238 :
1239 : ! **************************************************************************************************
1240 : !> \brief ...
1241 : !> \param iw ...
1242 : !> \param e_gspace ...
1243 : !> \param e_rspace ...
1244 : !> \param e_bonded ...
1245 : !> \param e_self ...
1246 : !> \param e_neut ...
1247 : !> \author Teodoro Laino [tlaino] - University of Zurich - 12.2007
1248 : ! **************************************************************************************************
1249 3736 : SUBROUTINE ewald_multipole_print(iw, e_gspace, e_rspace, e_bonded, e_self, e_neut)
1250 :
1251 : INTEGER, INTENT(IN) :: iw
1252 : REAL(KIND=dp), INTENT(IN) :: e_gspace, e_rspace, e_bonded, e_self, &
1253 : e_neut
1254 :
1255 3736 : IF (iw > 0) THEN
1256 642 : WRITE (iw, '( A, A )') ' *********************************', &
1257 1284 : '**********************************************'
1258 642 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL GSPACE ENERGY', &
1259 1284 : '[hartree]', '= ', e_gspace
1260 642 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL RSPACE ENERGY', &
1261 1284 : '[hartree]', '= ', e_rspace
1262 642 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' BONDED CORRECTION', &
1263 1284 : '[hartree]', '= ', e_bonded
1264 642 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' SELF ENERGY CORRECTION', &
1265 1284 : '[hartree]', '= ', e_self
1266 642 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' NEUTRALIZ. BCKGR. ENERGY', &
1267 1284 : '[hartree]', '= ', e_neut
1268 642 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' TOTAL ELECTROSTATIC EN.', &
1269 1284 : '[hartree]', '= ', e_rspace + e_bonded + e_gspace + e_self + e_neut
1270 642 : WRITE (iw, '( A, A )') ' *********************************', &
1271 1284 : '**********************************************'
1272 : END IF
1273 3736 : END SUBROUTINE ewald_multipole_print
1274 :
1275 : ! **************************************************************************************************
1276 : !> \brief Debug routines for multipoles
1277 : !> \param ewald_env ...
1278 : !> \param ewald_pw ...
1279 : !> \param nonbond_env ...
1280 : !> \param cell ...
1281 : !> \param particle_set ...
1282 : !> \param local_particles ...
1283 : !> \param iw ...
1284 : !> \param debug_r_space ...
1285 : !> \date 05.2008
1286 : !> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008
1287 : ! **************************************************************************************************
1288 0 : SUBROUTINE debug_ewald_multipoles(ewald_env, ewald_pw, nonbond_env, cell, &
1289 : particle_set, local_particles, iw, debug_r_space)
1290 : TYPE(ewald_environment_type), POINTER :: ewald_env
1291 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1292 : TYPE(fist_nonbond_env_type), POINTER :: nonbond_env
1293 : TYPE(cell_type), POINTER :: cell
1294 : TYPE(particle_type), DIMENSION(:), &
1295 : POINTER :: particle_set
1296 : TYPE(distribution_1d_type), POINTER :: local_particles
1297 : INTEGER, INTENT(IN) :: iw
1298 : LOGICAL, INTENT(IN) :: debug_r_space
1299 :
1300 : INTEGER :: nparticles
1301 : LOGICAL, DIMENSION(3) :: task
1302 : REAL(KIND=dp) :: e_neut, e_self, g_energy, &
1303 : r_energy, debug_energy
1304 0 : REAL(KIND=dp), POINTER, DIMENSION(:) :: charges
1305 : REAL(KIND=dp), POINTER, &
1306 0 : DIMENSION(:, :) :: dipoles, g_forces, g_pv, &
1307 0 : r_forces, r_pv, e_field1, &
1308 0 : e_field2
1309 : REAL(KIND=dp), POINTER, &
1310 0 : DIMENSION(:, :, :) :: quadrupoles
1311 : TYPE(rng_stream_type) :: random_stream
1312 : TYPE(multi_charge_type), DIMENSION(:), &
1313 0 : POINTER :: multipoles
1314 :
1315 0 : NULLIFY (multipoles, charges, dipoles, g_forces, g_pv, &
1316 0 : r_forces, r_pv, e_field1, e_field2)
1317 : random_stream = rng_stream_type(name="DEBUG_EWALD_MULTIPOLE", &
1318 0 : distribution_type=UNIFORM)
1319 : ! check: charge - charge
1320 0 : task = .FALSE.
1321 0 : nparticles = SIZE(particle_set)
1322 :
1323 : ! Allocate charges, dipoles, quadrupoles
1324 0 : ALLOCATE (charges(nparticles))
1325 0 : ALLOCATE (dipoles(3, nparticles))
1326 0 : ALLOCATE (quadrupoles(3, 3, nparticles))
1327 :
1328 : ! Allocate arrays for forces
1329 0 : ALLOCATE (r_forces(3, nparticles))
1330 0 : ALLOCATE (g_forces(3, nparticles))
1331 0 : ALLOCATE (e_field1(3, nparticles))
1332 0 : ALLOCATE (e_field2(3, nparticles))
1333 0 : ALLOCATE (g_pv(3, 3))
1334 0 : ALLOCATE (r_pv(3, 3))
1335 :
1336 : ! Debug CHARGES-CHARGES
1337 0 : task(1) = .TRUE.
1338 0 : charges = 0.0_dp
1339 0 : dipoles = 0.0_dp
1340 0 : quadrupoles = 0.0_dp
1341 0 : r_forces = 0.0_dp
1342 0 : g_forces = 0.0_dp
1343 0 : e_field1 = 0.0_dp
1344 0 : e_field2 = 0.0_dp
1345 0 : g_pv = 0.0_dp
1346 0 : r_pv = 0.0_dp
1347 0 : g_energy = 0.0_dp
1348 0 : r_energy = 0.0_dp
1349 : e_neut = 0.0_dp
1350 : e_self = 0.0_dp
1351 :
1352 : CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, &
1353 0 : random_stream=random_stream, charges=charges)
1354 : CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles, "CHARGE", echarge=1.0_dp, &
1355 0 : random_stream=random_stream, charges=charges)
1356 : CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
1357 0 : debug_r_space)
1358 :
1359 0 : WRITE (iw, *) "DEBUG ENERGY (CHARGE-CHARGE): ", debug_energy
1360 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, &
1361 : particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
1362 : task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., &
1363 : charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
1364 0 : forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.)
1365 0 : CALL release_multi_type(multipoles)
1366 :
1367 : ! Debug CHARGES-DIPOLES
1368 0 : task(1) = .TRUE.
1369 0 : task(2) = .TRUE.
1370 0 : charges = 0.0_dp
1371 0 : dipoles = 0.0_dp
1372 0 : quadrupoles = 0.0_dp
1373 0 : r_forces = 0.0_dp
1374 0 : g_forces = 0.0_dp
1375 0 : e_field1 = 0.0_dp
1376 0 : e_field2 = 0.0_dp
1377 0 : g_pv = 0.0_dp
1378 0 : r_pv = 0.0_dp
1379 0 : g_energy = 0.0_dp
1380 0 : r_energy = 0.0_dp
1381 : e_neut = 0.0_dp
1382 : e_self = 0.0_dp
1383 :
1384 : CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, &
1385 0 : random_stream=random_stream, charges=charges)
1386 : CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles, "DIPOLE", echarge=0.5_dp, &
1387 0 : random_stream=random_stream, dipoles=dipoles)
1388 0 : WRITE (iw, '("CHARGES",F15.9)') charges
1389 0 : WRITE (iw, '("DIPOLES",3F15.9)') dipoles
1390 : CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
1391 0 : debug_r_space)
1392 :
1393 0 : WRITE (iw, *) "DEBUG ENERGY (CHARGE-DIPOLE): ", debug_energy
1394 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, &
1395 : particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
1396 : task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., &
1397 : charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
1398 0 : forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.)
1399 0 : CALL release_multi_type(multipoles)
1400 :
1401 : ! Debug DIPOLES-DIPOLES
1402 0 : task(2) = .TRUE.
1403 0 : charges = 0.0_dp
1404 0 : dipoles = 0.0_dp
1405 0 : quadrupoles = 0.0_dp
1406 0 : r_forces = 0.0_dp
1407 0 : g_forces = 0.0_dp
1408 0 : e_field1 = 0.0_dp
1409 0 : e_field2 = 0.0_dp
1410 0 : g_pv = 0.0_dp
1411 0 : r_pv = 0.0_dp
1412 0 : g_energy = 0.0_dp
1413 0 : r_energy = 0.0_dp
1414 : e_neut = 0.0_dp
1415 : e_self = 0.0_dp
1416 :
1417 : CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "DIPOLE", echarge=10000.0_dp, &
1418 0 : random_stream=random_stream, dipoles=dipoles)
1419 : CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles, "DIPOLE", echarge=20000._dp, &
1420 0 : random_stream=random_stream, dipoles=dipoles)
1421 0 : WRITE (iw, '("DIPOLES",3F15.9)') dipoles
1422 : CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
1423 0 : debug_r_space)
1424 :
1425 0 : WRITE (iw, *) "DEBUG ENERGY (DIPOLE-DIPOLE): ", debug_energy
1426 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, &
1427 : particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
1428 : task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., &
1429 : charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
1430 0 : forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.)
1431 0 : CALL release_multi_type(multipoles)
1432 :
1433 : ! Debug CHARGES-QUADRUPOLES
1434 0 : task(1) = .TRUE.
1435 0 : task(3) = .TRUE.
1436 0 : charges = 0.0_dp
1437 0 : dipoles = 0.0_dp
1438 0 : quadrupoles = 0.0_dp
1439 0 : r_forces = 0.0_dp
1440 0 : g_forces = 0.0_dp
1441 0 : e_field1 = 0.0_dp
1442 0 : e_field2 = 0.0_dp
1443 0 : g_pv = 0.0_dp
1444 0 : r_pv = 0.0_dp
1445 0 : g_energy = 0.0_dp
1446 0 : r_energy = 0.0_dp
1447 : e_neut = 0.0_dp
1448 : e_self = 0.0_dp
1449 :
1450 : CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "CHARGE", echarge=-1.0_dp, &
1451 0 : random_stream=random_stream, charges=charges)
1452 : CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles, "QUADRUPOLE", echarge=10.0_dp, &
1453 0 : random_stream=random_stream, quadrupoles=quadrupoles)
1454 0 : WRITE (iw, '("CHARGES",F15.9)') charges
1455 0 : WRITE (iw, '("QUADRUPOLES",9F15.9)') quadrupoles
1456 : CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
1457 0 : debug_r_space)
1458 :
1459 0 : WRITE (iw, *) "DEBUG ENERGY (CHARGE-QUADRUPOLE): ", debug_energy
1460 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, &
1461 : particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
1462 : task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., &
1463 : charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
1464 0 : forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.)
1465 0 : CALL release_multi_type(multipoles)
1466 :
1467 : ! Debug DIPOLES-QUADRUPOLES
1468 0 : task(2) = .TRUE.
1469 0 : task(3) = .TRUE.
1470 0 : charges = 0.0_dp
1471 0 : dipoles = 0.0_dp
1472 0 : quadrupoles = 0.0_dp
1473 0 : r_forces = 0.0_dp
1474 0 : g_forces = 0.0_dp
1475 0 : e_field1 = 0.0_dp
1476 0 : e_field2 = 0.0_dp
1477 0 : g_pv = 0.0_dp
1478 0 : r_pv = 0.0_dp
1479 0 : g_energy = 0.0_dp
1480 0 : r_energy = 0.0_dp
1481 : e_neut = 0.0_dp
1482 : e_self = 0.0_dp
1483 :
1484 : CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "DIPOLE", echarge=10000.0_dp, &
1485 0 : random_stream=random_stream, dipoles=dipoles)
1486 : CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles, "QUADRUPOLE", echarge=10000.0_dp, &
1487 0 : random_stream=random_stream, quadrupoles=quadrupoles)
1488 0 : WRITE (iw, '("DIPOLES",3F15.9)') dipoles
1489 0 : WRITE (iw, '("QUADRUPOLES",9F15.9)') quadrupoles
1490 : CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
1491 0 : debug_r_space)
1492 :
1493 0 : WRITE (iw, *) "DEBUG ENERGY (DIPOLE-QUADRUPOLE): ", debug_energy
1494 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, &
1495 : particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
1496 : task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., &
1497 : charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
1498 0 : forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.)
1499 0 : CALL release_multi_type(multipoles)
1500 :
1501 : ! Debug QUADRUPOLES-QUADRUPOLES
1502 0 : task(3) = .TRUE.
1503 0 : charges = 0.0_dp
1504 0 : dipoles = 0.0_dp
1505 0 : quadrupoles = 0.0_dp
1506 0 : r_forces = 0.0_dp
1507 0 : g_forces = 0.0_dp
1508 0 : e_field1 = 0.0_dp
1509 0 : e_field2 = 0.0_dp
1510 0 : g_pv = 0.0_dp
1511 0 : r_pv = 0.0_dp
1512 0 : g_energy = 0.0_dp
1513 0 : r_energy = 0.0_dp
1514 : e_neut = 0.0_dp
1515 : e_self = 0.0_dp
1516 :
1517 : CALL create_multi_type(multipoles, nparticles, 1, nparticles/2, "QUADRUPOLE", echarge=-20000.0_dp, &
1518 0 : random_stream=random_stream, quadrupoles=quadrupoles)
1519 : CALL create_multi_type(multipoles, nparticles, nparticles/2 + 1, nparticles, "QUADRUPOLE", echarge=10000.0_dp, &
1520 0 : random_stream=random_stream, quadrupoles=quadrupoles)
1521 0 : WRITE (iw, '("QUADRUPOLES",9F15.9)') quadrupoles
1522 : CALL debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, debug_energy, &
1523 0 : debug_r_space)
1524 :
1525 0 : WRITE (iw, *) "DEBUG ENERGY (QUADRUPOLE-QUADRUPOLE): ", debug_energy
1526 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, &
1527 : particle_set, local_particles, g_energy, r_energy, e_neut, e_self, &
1528 : task, do_correction_bonded=.FALSE., do_forces=.TRUE., do_stress=.TRUE., do_efield=.FALSE., &
1529 : charges=charges, dipoles=dipoles, quadrupoles=quadrupoles, forces_local=g_forces, &
1530 0 : forces_glob=r_forces, pv_local=g_pv, pv_glob=r_pv, iw=iw, do_debug=.FALSE.)
1531 0 : CALL release_multi_type(multipoles)
1532 :
1533 0 : DEALLOCATE (charges)
1534 0 : DEALLOCATE (dipoles)
1535 0 : DEALLOCATE (quadrupoles)
1536 0 : DEALLOCATE (r_forces)
1537 0 : DEALLOCATE (g_forces)
1538 0 : DEALLOCATE (e_field1)
1539 0 : DEALLOCATE (e_field2)
1540 0 : DEALLOCATE (g_pv)
1541 0 : DEALLOCATE (r_pv)
1542 :
1543 : CONTAINS
1544 : ! **************************************************************************************************
1545 : !> \brief Debug routines for multipoles - low level - charge interactions
1546 : !> \param particle_set ...
1547 : !> \param cell ...
1548 : !> \param nonbond_env ...
1549 : !> \param multipoles ...
1550 : !> \param energy ...
1551 : !> \param debug_r_space ...
1552 : !> \date 05.2008
1553 : !> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008
1554 : ! **************************************************************************************************
1555 0 : SUBROUTINE debug_ewald_multipole_low(particle_set, cell, nonbond_env, multipoles, &
1556 : energy, debug_r_space)
1557 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1558 : TYPE(cell_type), POINTER :: cell
1559 : TYPE(fist_nonbond_env_type), POINTER :: nonbond_env
1560 : TYPE(multi_charge_type), DIMENSION(:), POINTER :: multipoles
1561 : REAL(KIND=dp), INTENT(OUT) :: energy
1562 : LOGICAL, INTENT(IN) :: debug_r_space
1563 :
1564 : INTEGER :: atom_a, atom_b, icell, iend, igrp, &
1565 : ikind, ilist, ipair, istart, jcell, &
1566 : jkind, k, k1, kcell, l, l1, ncells, &
1567 : nkinds, npairs
1568 0 : INTEGER, DIMENSION(:, :), POINTER :: list
1569 : REAL(KIND=dp) :: fac_ij, q, r, rab2, rab2_max
1570 : REAL(KIND=dp), DIMENSION(3) :: cell_v, cvi, rab, rab0, rm
1571 : TYPE(fist_neighbor_type), POINTER :: nonbonded
1572 : TYPE(neighbor_kind_pairs_type), POINTER :: neighbor_kind_pair
1573 0 : TYPE(pos_type), DIMENSION(:), POINTER :: r_last_update, r_last_update_pbc
1574 :
1575 0 : energy = 0.0_dp
1576 : CALL fist_nonbond_env_get(nonbond_env, nonbonded=nonbonded, natom_types=nkinds, &
1577 0 : r_last_update=r_last_update, r_last_update_pbc=r_last_update_pbc)
1578 0 : rab2_max = HUGE(0.0_dp)
1579 0 : IF (debug_r_space) THEN
1580 : ! This debugs the real space part of the multipole Ewald summation scheme
1581 : ! Starting the force loop
1582 0 : Lists: DO ilist = 1, nonbonded%nlists
1583 0 : neighbor_kind_pair => nonbonded%neighbor_kind_pairs(ilist)
1584 0 : npairs = neighbor_kind_pair%npairs
1585 0 : IF (npairs == 0) CYCLE Lists
1586 0 : list => neighbor_kind_pair%list
1587 0 : cvi = neighbor_kind_pair%cell_vector
1588 0 : cell_v = MATMUL(cell%hmat, cvi)
1589 0 : Kind_Group_Loop: DO igrp = 1, neighbor_kind_pair%ngrp_kind
1590 0 : istart = neighbor_kind_pair%grp_kind_start(igrp)
1591 0 : iend = neighbor_kind_pair%grp_kind_end(igrp)
1592 0 : ikind = neighbor_kind_pair%ij_kind(1, igrp)
1593 0 : jkind = neighbor_kind_pair%ij_kind(2, igrp)
1594 0 : Pairs: DO ipair = istart, iend
1595 0 : fac_ij = 1.0_dp
1596 0 : atom_a = list(1, ipair)
1597 0 : atom_b = list(2, ipair)
1598 0 : IF (atom_a == atom_b) fac_ij = 0.5_dp
1599 0 : rab = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
1600 0 : rab = rab + cell_v
1601 0 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
1602 0 : IF (rab2 <= rab2_max) THEN
1603 :
1604 0 : DO k = 1, SIZE(multipoles(atom_a)%charge_typ)
1605 0 : DO k1 = 1, SIZE(multipoles(atom_a)%charge_typ(k)%charge)
1606 :
1607 0 : DO l = 1, SIZE(multipoles(atom_b)%charge_typ)
1608 0 : DO l1 = 1, SIZE(multipoles(atom_b)%charge_typ(l)%charge)
1609 :
1610 0 : rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:, l1) - multipoles(atom_a)%charge_typ(k)%pos(:, k1)
1611 0 : r = NORM2(rm)
1612 0 : q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1)
1613 0 : energy = energy + q/r*fac_ij
1614 : END DO
1615 : END DO
1616 :
1617 : END DO
1618 : END DO
1619 :
1620 : END IF
1621 : END DO Pairs
1622 : END DO Kind_Group_Loop
1623 : END DO Lists
1624 : ELSE
1625 0 : ncells = 6
1626 : !Debugs the sum of real + space terms.. (Charge-Charge and Charge-Dipole should be anyway wrong but
1627 : !all the other terms should be correct)
1628 0 : DO atom_a = 1, SIZE(particle_set)
1629 0 : DO atom_b = atom_a, SIZE(particle_set)
1630 0 : fac_ij = 1.0_dp
1631 0 : IF (atom_a == atom_b) fac_ij = 0.5_dp
1632 0 : rab0 = r_last_update_pbc(atom_b)%r - r_last_update_pbc(atom_a)%r
1633 : ! Loop over cells
1634 0 : DO icell = -ncells, ncells
1635 0 : DO jcell = -ncells, ncells
1636 0 : DO kcell = -ncells, ncells
1637 0 : cell_v = MATMUL(cell%hmat, REAL([icell, jcell, kcell], KIND=dp))
1638 0 : IF (ALL(cell_v == 0.0_dp) .AND. (atom_a == atom_b)) CYCLE
1639 0 : rab = rab0 + cell_v
1640 0 : rab2 = rab(1)**2 + rab(2)**2 + rab(3)**2
1641 0 : IF (rab2 <= rab2_max) THEN
1642 :
1643 0 : DO k = 1, SIZE(multipoles(atom_a)%charge_typ)
1644 0 : DO k1 = 1, SIZE(multipoles(atom_a)%charge_typ(k)%charge)
1645 :
1646 0 : DO l = 1, SIZE(multipoles(atom_b)%charge_typ)
1647 0 : DO l1 = 1, SIZE(multipoles(atom_b)%charge_typ(l)%charge)
1648 :
1649 0 : rm = rab + multipoles(atom_b)%charge_typ(l)%pos(:, l1) - multipoles(atom_a)%charge_typ(k)%pos(:, k1)
1650 0 : r = NORM2(rm)
1651 0 : q = multipoles(atom_b)%charge_typ(l)%charge(l1)*multipoles(atom_a)%charge_typ(k)%charge(k1)
1652 0 : energy = energy + q/r*fac_ij
1653 : END DO
1654 : END DO
1655 :
1656 : END DO
1657 : END DO
1658 :
1659 : END IF
1660 : END DO
1661 : END DO
1662 : END DO
1663 : END DO
1664 : END DO
1665 : END IF
1666 0 : END SUBROUTINE debug_ewald_multipole_low
1667 :
1668 : ! **************************************************************************************************
1669 : !> \brief create multi_type for multipoles
1670 : !> \param multipoles ...
1671 : !> \param idim ...
1672 : !> \param istart ...
1673 : !> \param iend ...
1674 : !> \param label ...
1675 : !> \param echarge ...
1676 : !> \param random_stream ...
1677 : !> \param charges ...
1678 : !> \param dipoles ...
1679 : !> \param quadrupoles ...
1680 : !> \date 05.2008
1681 : !> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008
1682 : ! **************************************************************************************************
1683 0 : SUBROUTINE create_multi_type(multipoles, idim, istart, iend, label, echarge, &
1684 : random_stream, charges, dipoles, quadrupoles)
1685 : TYPE(multi_charge_type), DIMENSION(:), POINTER :: multipoles
1686 : INTEGER, INTENT(IN) :: idim, istart, iend
1687 : CHARACTER(LEN=*), INTENT(IN) :: label
1688 : REAL(KIND=dp), INTENT(IN) :: echarge
1689 : TYPE(rng_stream_type), INTENT(INOUT) :: random_stream
1690 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges
1691 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
1692 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
1693 : POINTER :: quadrupoles
1694 :
1695 : INTEGER :: i, isize, k, l, m
1696 : REAL(KIND=dp) :: dx, r2, rvec(3), rvec1(3), rvec2(3)
1697 :
1698 0 : IF (ASSOCIATED(multipoles)) THEN
1699 0 : CPASSERT(SIZE(multipoles) == idim)
1700 : ELSE
1701 0 : ALLOCATE (multipoles(idim))
1702 0 : DO i = 1, idim
1703 0 : NULLIFY (multipoles(i)%charge_typ)
1704 : END DO
1705 : END IF
1706 0 : DO i = istart, iend
1707 0 : IF (ASSOCIATED(multipoles(i)%charge_typ)) THEN
1708 : ! make a copy of the array and enlarge the array type by 1
1709 0 : isize = SIZE(multipoles(i)%charge_typ) + 1
1710 : ELSE
1711 0 : isize = 1
1712 : END IF
1713 0 : CALL reallocate_charge_type(multipoles(i)%charge_typ, 1, isize)
1714 0 : SELECT CASE (label)
1715 : CASE ("CHARGE")
1716 0 : CPASSERT(PRESENT(charges))
1717 0 : CPASSERT(ASSOCIATED(charges))
1718 0 : ALLOCATE (multipoles(i)%charge_typ(isize)%charge(1))
1719 0 : ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 1))
1720 :
1721 0 : multipoles(i)%charge_typ(isize)%charge(1) = echarge
1722 0 : multipoles(i)%charge_typ(isize)%pos(1:3, 1) = 0.0_dp
1723 0 : charges(i) = charges(i) + echarge
1724 : CASE ("DIPOLE")
1725 0 : dx = 1.0E-4_dp
1726 0 : CPASSERT(PRESENT(dipoles))
1727 0 : CPASSERT(ASSOCIATED(dipoles))
1728 0 : ALLOCATE (multipoles(i)%charge_typ(isize)%charge(2))
1729 0 : ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 2))
1730 0 : CALL random_stream%fill(rvec)
1731 0 : rvec = rvec/(2.0_dp*NORM2(rvec))*dx
1732 0 : multipoles(i)%charge_typ(isize)%charge(1) = echarge
1733 0 : multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec
1734 0 : multipoles(i)%charge_typ(isize)%charge(2) = -echarge
1735 0 : multipoles(i)%charge_typ(isize)%pos(1:3, 2) = -rvec
1736 :
1737 0 : dipoles(:, i) = dipoles(:, i) + 2.0_dp*echarge*rvec
1738 : CASE ("QUADRUPOLE")
1739 0 : dx = 1.0E-2_dp
1740 0 : CPASSERT(PRESENT(quadrupoles))
1741 0 : CPASSERT(ASSOCIATED(quadrupoles))
1742 0 : ALLOCATE (multipoles(i)%charge_typ(isize)%charge(4))
1743 0 : ALLOCATE (multipoles(i)%charge_typ(isize)%pos(3, 4))
1744 0 : CALL random_stream%fill(rvec1)
1745 0 : CALL random_stream%fill(rvec2)
1746 0 : rvec1 = rvec1/NORM2(rvec1)
1747 0 : rvec2 = rvec2 - DOT_PRODUCT(rvec2, rvec1)*rvec1
1748 0 : rvec2 = rvec2/NORM2(rvec2)
1749 : !
1750 0 : rvec1 = rvec1/2.0_dp*dx
1751 0 : rvec2 = rvec2/2.0_dp*dx
1752 : ! + (4) ^ - (1)
1753 : ! |rvec2
1754 : ! |
1755 : ! 0------> rvec1
1756 : !
1757 : !
1758 : ! - (3) + (2)
1759 0 : multipoles(i)%charge_typ(isize)%charge(1) = -echarge
1760 0 : multipoles(i)%charge_typ(isize)%pos(1:3, 1) = rvec1 + rvec2
1761 0 : multipoles(i)%charge_typ(isize)%charge(2) = echarge
1762 0 : multipoles(i)%charge_typ(isize)%pos(1:3, 2) = rvec1 - rvec2
1763 0 : multipoles(i)%charge_typ(isize)%charge(3) = -echarge
1764 0 : multipoles(i)%charge_typ(isize)%pos(1:3, 3) = -rvec1 - rvec2
1765 0 : multipoles(i)%charge_typ(isize)%charge(4) = echarge
1766 0 : multipoles(i)%charge_typ(isize)%pos(1:3, 4) = -rvec1 + rvec2
1767 :
1768 0 : DO k = 1, 4
1769 0 : r2 = DOT_PRODUCT(multipoles(i)%charge_typ(isize)%pos(:, k), multipoles(i)%charge_typ(isize)%pos(:, k))
1770 0 : DO l = 1, 3
1771 0 : DO m = 1, 3
1772 : quadrupoles(m, l, i) = quadrupoles(m, l, i) + 3.0_dp*0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)* &
1773 : multipoles(i)%charge_typ(isize)%pos(l, k)* &
1774 0 : multipoles(i)%charge_typ(isize)%pos(m, k)
1775 0 : IF (m == l) quadrupoles(m, l, i) = quadrupoles(m, l, i) - 0.5_dp*multipoles(i)%charge_typ(isize)%charge(k)*r2
1776 : END DO
1777 : END DO
1778 : END DO
1779 :
1780 : END SELECT
1781 : END DO
1782 0 : END SUBROUTINE create_multi_type
1783 :
1784 : ! **************************************************************************************************
1785 : !> \brief release multi_type for multipoles
1786 : !> \param multipoles ...
1787 : !> \date 05.2008
1788 : !> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008
1789 : ! **************************************************************************************************
1790 0 : SUBROUTINE release_multi_type(multipoles)
1791 : TYPE(multi_charge_type), DIMENSION(:), POINTER :: multipoles
1792 :
1793 : INTEGER :: i, j
1794 :
1795 0 : IF (ASSOCIATED(multipoles)) THEN
1796 0 : DO i = 1, SIZE(multipoles)
1797 0 : DO j = 1, SIZE(multipoles(i)%charge_typ)
1798 0 : DEALLOCATE (multipoles(i)%charge_typ(j)%charge)
1799 0 : DEALLOCATE (multipoles(i)%charge_typ(j)%pos)
1800 : END DO
1801 0 : DEALLOCATE (multipoles(i)%charge_typ)
1802 : END DO
1803 : END IF
1804 0 : END SUBROUTINE release_multi_type
1805 :
1806 : ! **************************************************************************************************
1807 : !> \brief reallocates multi_type for multipoles
1808 : !> \param charge_typ ...
1809 : !> \param istart ...
1810 : !> \param iend ...
1811 : !> \date 05.2008
1812 : !> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008
1813 : ! **************************************************************************************************
1814 0 : SUBROUTINE reallocate_charge_type(charge_typ, istart, iend)
1815 : TYPE(charge_mono_type), DIMENSION(:), POINTER :: charge_typ
1816 : INTEGER, INTENT(IN) :: istart, iend
1817 :
1818 : INTEGER :: i, isize, j, jsize, jsize1, jsize2
1819 0 : TYPE(charge_mono_type), DIMENSION(:), POINTER :: charge_typ_bk
1820 :
1821 0 : IF (ASSOCIATED(charge_typ)) THEN
1822 0 : isize = SIZE(charge_typ)
1823 0 : ALLOCATE (charge_typ_bk(1:isize))
1824 0 : DO j = 1, isize
1825 0 : jsize = SIZE(charge_typ(j)%charge)
1826 0 : ALLOCATE (charge_typ_bk(j)%charge(jsize))
1827 0 : jsize1 = SIZE(charge_typ(j)%pos, 1)
1828 0 : jsize2 = SIZE(charge_typ(j)%pos, 2)
1829 0 : ALLOCATE (charge_typ_bk(j)%pos(jsize1, jsize2))
1830 0 : charge_typ_bk(j)%pos = charge_typ(j)%pos
1831 0 : charge_typ_bk(j)%charge = charge_typ(j)%charge
1832 : END DO
1833 0 : DO j = 1, SIZE(charge_typ)
1834 0 : DEALLOCATE (charge_typ(j)%charge)
1835 0 : DEALLOCATE (charge_typ(j)%pos)
1836 : END DO
1837 0 : DEALLOCATE (charge_typ)
1838 : ! Reallocate
1839 0 : ALLOCATE (charge_typ_bk(istart:iend))
1840 0 : DO i = istart, isize
1841 0 : jsize = SIZE(charge_typ_bk(j)%charge)
1842 0 : ALLOCATE (charge_typ(j)%charge(jsize))
1843 0 : jsize1 = SIZE(charge_typ_bk(j)%pos, 1)
1844 0 : jsize2 = SIZE(charge_typ_bk(j)%pos, 2)
1845 0 : ALLOCATE (charge_typ(j)%pos(jsize1, jsize2))
1846 0 : charge_typ(j)%pos = charge_typ_bk(j)%pos
1847 0 : charge_typ(j)%charge = charge_typ_bk(j)%charge
1848 : END DO
1849 0 : DO j = 1, SIZE(charge_typ_bk)
1850 0 : DEALLOCATE (charge_typ_bk(j)%charge)
1851 0 : DEALLOCATE (charge_typ_bk(j)%pos)
1852 : END DO
1853 0 : DEALLOCATE (charge_typ_bk)
1854 : ELSE
1855 0 : ALLOCATE (charge_typ(istart:iend))
1856 : END IF
1857 :
1858 0 : END SUBROUTINE reallocate_charge_type
1859 :
1860 : END SUBROUTINE debug_ewald_multipoles
1861 :
1862 : ! **************************************************************************************************
1863 : !> \brief Routine to debug potential, field and electric field gradients
1864 : !> \param ewald_env ...
1865 : !> \param ewald_pw ...
1866 : !> \param nonbond_env ...
1867 : !> \param cell ...
1868 : !> \param particle_set ...
1869 : !> \param local_particles ...
1870 : !> \param radii ...
1871 : !> \param charges ...
1872 : !> \param dipoles ...
1873 : !> \param quadrupoles ...
1874 : !> \param task ...
1875 : !> \param iw ...
1876 : !> \param atomic_kind_set ...
1877 : !> \param mm_section ...
1878 : !> \date 05.2008
1879 : !> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008
1880 : ! **************************************************************************************************
1881 0 : SUBROUTINE debug_ewald_multipoles_fields(ewald_env, ewald_pw, nonbond_env, cell, &
1882 : particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw, &
1883 : atomic_kind_set, mm_section)
1884 : TYPE(ewald_environment_type), POINTER :: ewald_env
1885 : TYPE(ewald_pw_type), POINTER :: ewald_pw
1886 : TYPE(fist_nonbond_env_type), POINTER :: nonbond_env
1887 : TYPE(cell_type), POINTER :: cell
1888 : TYPE(particle_type), POINTER :: particle_set(:)
1889 : TYPE(distribution_1d_type), POINTER :: local_particles
1890 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges
1891 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
1892 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
1893 : POINTER :: quadrupoles
1894 : LOGICAL, DIMENSION(3), INTENT(IN) :: task
1895 : INTEGER, INTENT(IN) :: iw
1896 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
1897 : TYPE(section_vals_type), POINTER :: mm_section
1898 :
1899 : INTEGER :: i, iparticle_kind, j, k, &
1900 : nparticle_local, nparticles
1901 : REAL(KIND=dp) :: coord(3), dq, e_neut, e_self, efield1n(3), efield2n(3, 3), ene(2), &
1902 : energy_glob, energy_local, enev(3, 2), o_tot_ene, pot, pv_glob(3, 3), pv_local(3, 3), &
1903 : tot_ene
1904 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efield1, efield2, forces_glob, &
1905 : forces_local
1906 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0, lcharges
1907 : TYPE(cp_logger_type), POINTER :: logger
1908 0 : TYPE(particle_type), DIMENSION(:), POINTER :: core_particle_set, shell_particle_set
1909 :
1910 0 : NULLIFY (lcharges, shell_particle_set, core_particle_set)
1911 0 : NULLIFY (logger)
1912 0 : logger => cp_get_default_logger()
1913 :
1914 0 : nparticles = SIZE(particle_set)
1915 0 : nparticle_local = 0
1916 0 : DO iparticle_kind = 1, SIZE(local_particles%n_el)
1917 0 : nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind)
1918 : END DO
1919 0 : ALLOCATE (lcharges(nparticles))
1920 0 : ALLOCATE (forces_glob(3, nparticles))
1921 0 : ALLOCATE (forces_local(3, nparticle_local))
1922 0 : ALLOCATE (efield0(nparticles))
1923 0 : ALLOCATE (efield1(3, nparticles))
1924 0 : ALLOCATE (efield2(9, nparticles))
1925 0 : forces_glob = 0.0_dp
1926 0 : forces_local = 0.0_dp
1927 0 : efield0 = 0.0_dp
1928 0 : efield1 = 0.0_dp
1929 0 : efield2 = 0.0_dp
1930 0 : pv_local = 0.0_dp
1931 0 : pv_glob = 0.0_dp
1932 0 : energy_glob = 0.0_dp
1933 0 : energy_local = 0.0_dp
1934 : e_neut = 0.0_dp
1935 : e_self = 0.0_dp
1936 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, &
1937 : local_particles, energy_local, energy_glob, e_neut, e_self, task, .FALSE., .TRUE., .TRUE., &
1938 : .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, &
1939 0 : efield0, efield1, efield2, iw, do_debug=.FALSE.)
1940 0 : o_tot_ene = energy_local + energy_glob + e_neut + e_self
1941 0 : WRITE (iw, *) "TOTAL ENERGY :: ========>", o_tot_ene
1942 : ! Debug Potential
1943 0 : dq = 0.001_dp
1944 0 : tot_ene = 0.0_dp
1945 0 : DO i = 1, nparticles
1946 0 : DO k = 1, 2
1947 0 : lcharges = charges
1948 0 : lcharges(i) = charges(i) + (-1.0_dp)**k*dq
1949 0 : forces_glob = 0.0_dp
1950 0 : forces_local = 0.0_dp
1951 0 : pv_local = 0.0_dp
1952 0 : pv_glob = 0.0_dp
1953 0 : energy_glob = 0.0_dp
1954 0 : energy_local = 0.0_dp
1955 : e_neut = 0.0_dp
1956 : e_self = 0.0_dp
1957 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, &
1958 : local_particles, energy_local, energy_glob, e_neut, e_self, &
1959 : task, .FALSE., .FALSE., .FALSE., .FALSE., radii, &
1960 0 : lcharges, dipoles, quadrupoles, iw=iw, do_debug=.FALSE.)
1961 0 : ene(k) = energy_local + energy_glob + e_neut + e_self
1962 : END DO
1963 0 : pot = (ene(2) - ene(1))/(2.0_dp*dq)
1964 0 : WRITE (iw, '(A,I8,3(A,F15.9))') "POTENTIAL FOR ATOM: ", i, " NUMERICAL: ", pot, " ANALYTICAL: ", efield0(i), &
1965 0 : " ERROR: ", pot - efield0(i)
1966 0 : tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i)
1967 : END DO
1968 0 : WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
1969 0 : WRITE (iw, '(/,/,/)')
1970 : ! Debug Field
1971 0 : dq = 0.001_dp
1972 0 : DO i = 1, nparticles
1973 0 : coord = particle_set(i)%r
1974 0 : DO j = 1, 3
1975 0 : DO k = 1, 2
1976 0 : particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k*dq
1977 :
1978 : ! Rebuild neighbor lists
1979 : CALL list_control(atomic_kind_set, particle_set, local_particles, &
1980 : cell, nonbond_env, logger%para_env, mm_section, &
1981 0 : shell_particle_set, core_particle_set)
1982 :
1983 0 : forces_glob = 0.0_dp
1984 0 : forces_local = 0.0_dp
1985 0 : pv_local = 0.0_dp
1986 0 : pv_glob = 0.0_dp
1987 0 : energy_glob = 0.0_dp
1988 0 : energy_local = 0.0_dp
1989 : e_neut = 0.0_dp
1990 : e_self = 0.0_dp
1991 0 : efield0 = 0.0_dp
1992 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, &
1993 : local_particles, energy_local, energy_glob, e_neut, e_self, &
1994 : task, .FALSE., .TRUE., .TRUE., .TRUE., radii, &
1995 : charges, dipoles, quadrupoles, forces_local, forces_glob, &
1996 0 : pv_local, pv_glob, efield0, iw=iw, do_debug=.FALSE.)
1997 0 : ene(k) = efield0(i)
1998 0 : particle_set(i)%r(j) = coord(j)
1999 : END DO
2000 0 : efield1n(j) = -(ene(2) - ene(1))/(2.0_dp*dq)
2001 : END DO
2002 0 : WRITE (iw, '(/,A,I8)') "FIELD FOR ATOM: ", i
2003 0 : WRITE (iw, '(A,3F15.9)') " NUMERICAL: ", efield1n, " ANALYTICAL: ", efield1(:, i), &
2004 0 : " ERROR: ", efield1n - efield1(:, i)
2005 0 : IF (task(2)) THEN
2006 0 : tot_ene = tot_ene - 0.5_dp*DOT_PRODUCT(efield1(:, i), dipoles(:, i))
2007 : END IF
2008 : END DO
2009 0 : WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
2010 :
2011 : ! Debug Field Gradient
2012 0 : dq = 0.0001_dp
2013 0 : DO i = 1, nparticles
2014 0 : coord = particle_set(i)%r
2015 0 : DO j = 1, 3
2016 0 : DO k = 1, 2
2017 0 : particle_set(i)%r(j) = coord(j) + (-1.0_dp)**k*dq
2018 :
2019 : ! Rebuild neighbor lists
2020 : CALL list_control(atomic_kind_set, particle_set, local_particles, &
2021 : cell, nonbond_env, logger%para_env, mm_section, &
2022 0 : shell_particle_set, core_particle_set)
2023 :
2024 0 : forces_glob = 0.0_dp
2025 0 : forces_local = 0.0_dp
2026 0 : pv_local = 0.0_dp
2027 0 : pv_glob = 0.0_dp
2028 0 : energy_glob = 0.0_dp
2029 0 : energy_local = 0.0_dp
2030 : e_neut = 0.0_dp
2031 : e_self = 0.0_dp
2032 0 : efield1 = 0.0_dp
2033 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, &
2034 : local_particles, energy_local, energy_glob, e_neut, e_self, &
2035 : task, .FALSE., .TRUE., .TRUE., .TRUE., radii, &
2036 : charges, dipoles, quadrupoles, forces_local, forces_glob, &
2037 0 : pv_local, pv_glob, efield1=efield1, iw=iw, do_debug=.FALSE.)
2038 0 : enev(:, k) = efield1(:, i)
2039 0 : particle_set(i)%r(j) = coord(j)
2040 : END DO
2041 0 : efield2n(:, j) = (enev(:, 2) - enev(:, 1))/(2.0_dp*dq)
2042 : END DO
2043 0 : WRITE (iw, '(/,A,I8)') "FIELD GRADIENT FOR ATOM: ", i
2044 0 : WRITE (iw, '(A,9F15.9)') " NUMERICAL: ", efield2n, &
2045 0 : " ANALYTICAL: ", efield2(:, i), &
2046 0 : " ERROR: ", RESHAPE(efield2n, [9]) - efield2(:, i)
2047 : END DO
2048 0 : END SUBROUTINE debug_ewald_multipoles_fields
2049 :
2050 : ! **************************************************************************************************
2051 : !> \brief Routine to debug potential, field and electric field gradients
2052 : !> \param ewald_env ...
2053 : !> \param ewald_pw ...
2054 : !> \param nonbond_env ...
2055 : !> \param cell ...
2056 : !> \param particle_set ...
2057 : !> \param local_particles ...
2058 : !> \param radii ...
2059 : !> \param charges ...
2060 : !> \param dipoles ...
2061 : !> \param quadrupoles ...
2062 : !> \param task ...
2063 : !> \param iw ...
2064 : !> \date 05.2008
2065 : !> \author Teodoro Laino [tlaino] - University of Zurich - 05.2008
2066 : ! **************************************************************************************************
2067 0 : SUBROUTINE debug_ewald_multipoles_fields2(ewald_env, ewald_pw, nonbond_env, cell, &
2068 : particle_set, local_particles, radii, charges, dipoles, quadrupoles, task, iw)
2069 : TYPE(ewald_environment_type), POINTER :: ewald_env
2070 : TYPE(ewald_pw_type), POINTER :: ewald_pw
2071 : TYPE(fist_nonbond_env_type), POINTER :: nonbond_env
2072 : TYPE(cell_type), POINTER :: cell
2073 : TYPE(particle_type), POINTER :: particle_set(:)
2074 : TYPE(distribution_1d_type), POINTER :: local_particles
2075 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: radii, charges
2076 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: dipoles
2077 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
2078 : POINTER :: quadrupoles
2079 : LOGICAL, DIMENSION(3), INTENT(IN) :: task
2080 : INTEGER, INTENT(IN) :: iw
2081 :
2082 : INTEGER :: i, ind, iparticle_kind, j, k, &
2083 : nparticle_local, nparticles
2084 : REAL(KIND=dp) :: e_neut, e_self, energy_glob, &
2085 : energy_local, o_tot_ene, prod, &
2086 : pv_glob(3, 3), pv_local(3, 3), tot_ene
2087 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efield1, efield2, forces_glob, &
2088 : forces_local
2089 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0
2090 : TYPE(cp_logger_type), POINTER :: logger
2091 :
2092 0 : NULLIFY (logger)
2093 0 : logger => cp_get_default_logger()
2094 :
2095 0 : nparticles = SIZE(particle_set)
2096 0 : nparticle_local = 0
2097 0 : DO iparticle_kind = 1, SIZE(local_particles%n_el)
2098 0 : nparticle_local = nparticle_local + local_particles%n_el(iparticle_kind)
2099 : END DO
2100 0 : ALLOCATE (forces_glob(3, nparticles))
2101 0 : ALLOCATE (forces_local(3, nparticle_local))
2102 0 : ALLOCATE (efield0(nparticles))
2103 0 : ALLOCATE (efield1(3, nparticles))
2104 0 : ALLOCATE (efield2(9, nparticles))
2105 0 : forces_glob = 0.0_dp
2106 0 : forces_local = 0.0_dp
2107 0 : efield0 = 0.0_dp
2108 0 : efield1 = 0.0_dp
2109 0 : efield2 = 0.0_dp
2110 0 : pv_local = 0.0_dp
2111 0 : pv_glob = 0.0_dp
2112 0 : energy_glob = 0.0_dp
2113 0 : energy_local = 0.0_dp
2114 : e_neut = 0.0_dp
2115 : e_self = 0.0_dp
2116 : CALL ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, &
2117 : local_particles, energy_local, energy_glob, e_neut, e_self, task, .FALSE., .TRUE., .TRUE., &
2118 : .TRUE., radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, &
2119 0 : efield0, efield1, efield2, iw, do_debug=.FALSE.)
2120 0 : o_tot_ene = energy_local + energy_glob + e_neut + e_self
2121 0 : WRITE (iw, *) "TOTAL ENERGY :: ========>", o_tot_ene
2122 :
2123 : ! Debug Potential
2124 0 : tot_ene = 0.0_dp
2125 0 : IF (task(1)) THEN
2126 0 : DO i = 1, nparticles
2127 0 : tot_ene = tot_ene + 0.5_dp*efield0(i)*charges(i)
2128 : END DO
2129 0 : WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
2130 0 : WRITE (iw, '(/,/,/)')
2131 : END IF
2132 :
2133 : ! Debug Field
2134 0 : IF (task(2)) THEN
2135 0 : DO i = 1, nparticles
2136 0 : tot_ene = tot_ene - 0.5_dp*DOT_PRODUCT(efield1(:, i), dipoles(:, i))
2137 : END DO
2138 0 : WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
2139 0 : WRITE (iw, '(/,/,/)')
2140 : END IF
2141 :
2142 : ! Debug Field Gradient
2143 0 : IF (task(3)) THEN
2144 0 : DO i = 1, nparticles
2145 : ind = 0
2146 : prod = 0.0_dp
2147 0 : DO j = 1, 3
2148 0 : DO k = 1, 3
2149 0 : ind = ind + 1
2150 0 : prod = prod + efield2(ind, i)*quadrupoles(j, k, i)
2151 : END DO
2152 : END DO
2153 0 : tot_ene = tot_ene - 0.5_dp*(1.0_dp/3.0_dp)*prod
2154 : END DO
2155 0 : WRITE (iw, *) "ENERGIES: ", o_tot_ene, tot_ene, o_tot_ene - tot_ene
2156 0 : WRITE (iw, '(/,/,/)')
2157 : END IF
2158 :
2159 0 : END SUBROUTINE debug_ewald_multipoles_fields2
2160 :
2161 0 : END MODULE ewalds_multipole
|
