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