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