Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \par History
10 : !> JGH (15-Mar-2001) : New routine ewald_setup (former pme_setup)
11 : !> JGH (23-Mar-2001) : Get rid of global variable ewald_grp
12 : ! **************************************************************************************************
13 : MODULE ewalds
14 :
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE bibliography, ONLY: Ewald1921,&
18 : cite_reference
19 : USE cell_types, ONLY: cell_type
20 : USE dg_rho0_types, ONLY: dg_rho0_type
21 : USE dg_types, ONLY: dg_get,&
22 : dg_type
23 : USE distribution_1d_types, ONLY: distribution_1d_type
24 : USE ewald_environment_types, ONLY: ewald_env_get,&
25 : ewald_environment_type
26 : USE ewald_pw_types, ONLY: ewald_pw_get,&
27 : ewald_pw_type
28 : USE kinds, ONLY: dp
29 : USE mathconstants, ONLY: fourpi,&
30 : oorootpi,&
31 : pi
32 : USE message_passing, ONLY: mp_comm_type
33 : USE particle_types, ONLY: particle_type
34 : USE pw_grid_types, ONLY: pw_grid_type
35 : USE pw_poisson_types, ONLY: do_ewald_none
36 : USE pw_pool_types, ONLY: pw_pool_type
37 : USE shell_potential_types, ONLY: get_shell,&
38 : shell_kind_type
39 : USE structure_factor_types, ONLY: structure_factor_type
40 : USE structure_factors, ONLY: structure_factor_allocate,&
41 : structure_factor_deallocate,&
42 : structure_factor_evaluate
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 : LOGICAL, PRIVATE, PARAMETER :: debug_this_module = .TRUE.
47 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ewalds'
48 :
49 : PRIVATE
50 : PUBLIC :: ewald_evaluate, ewald_self, ewald_self_atom, ewald_print
51 :
52 : CONTAINS
53 :
54 : ! **************************************************************************************************
55 : !> \brief computes the potential and the force from the g-space part of
56 : !> the 1/r potential
57 : !> Ref.: J.-P. Hansen, Enrico Fermi School, 1985
58 : !> Note: Only the positive G-vectors are used in the sum.
59 : !> \param ewald_env ...
60 : !> \param ewald_pw ...
61 : !> \param cell ...
62 : !> \param atomic_kind_set ...
63 : !> \param particle_set ...
64 : !> \param local_particles ...
65 : !> \param fg_coulomb ...
66 : !> \param vg_coulomb ...
67 : !> \param pv_g ...
68 : !> \param use_virial ...
69 : !> \param charges ...
70 : !> \param e_coulomb ...
71 : !> \param pv_coulomb ...
72 : !> \par History
73 : !> JGH (21-Feb-2001) : changed name
74 : !> \author CJM
75 : ! **************************************************************************************************
76 28739 : SUBROUTINE ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
77 28739 : local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial, charges, e_coulomb, &
78 : pv_coulomb)
79 : TYPE(ewald_environment_type), POINTER :: ewald_env
80 : TYPE(ewald_pw_type), POINTER :: ewald_pw
81 : TYPE(cell_type), POINTER :: cell
82 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
83 : TYPE(particle_type), POINTER :: particle_set(:)
84 : TYPE(distribution_1d_type), POINTER :: local_particles
85 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fg_coulomb
86 : REAL(KIND=dp), INTENT(OUT) :: vg_coulomb
87 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: pv_g
88 : LOGICAL, INTENT(IN) :: use_virial
89 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges, e_coulomb
90 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL, &
91 : POINTER :: pv_coulomb
92 :
93 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_evaluate'
94 :
95 : COMPLEX(KIND=dp) :: snode
96 28739 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: summe
97 : INTEGER :: gpt, handle, iparticle, iparticle_kind, &
98 : iparticle_local, lp, mp, nnodes, node, &
99 : np, nparticle_kind, nparticle_local
100 28739 : INTEGER, DIMENSION(:, :), POINTER :: bds
101 : LOGICAL :: atenergy, atstress, use_charge_array
102 : REAL(KIND=dp) :: alpha, denom, e_igdotr, factor, &
103 : four_alpha_sq, gauss, pref, q
104 : REAL(KIND=dp), DIMENSION(3) :: vec
105 28739 : REAL(KIND=dp), DIMENSION(:), POINTER :: charge
106 28739 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0
107 : TYPE(atomic_kind_type), POINTER :: atomic_kind
108 : TYPE(dg_rho0_type), POINTER :: dg_rho0
109 : TYPE(dg_type), POINTER :: dg
110 : TYPE(mp_comm_type) :: group
111 : TYPE(pw_grid_type), POINTER :: pw_grid
112 : TYPE(pw_pool_type), POINTER :: pw_pool
113 : TYPE(structure_factor_type) :: exp_igr
114 :
115 28739 : CALL timeset(routineN, handle)
116 28739 : CALL cite_reference(Ewald1921)
117 28739 : use_charge_array = PRESENT(charges)
118 28739 : IF (use_charge_array) use_charge_array = ASSOCIATED(charges)
119 28739 : atenergy = PRESENT(e_coulomb)
120 28739 : IF (atenergy) atenergy = ASSOCIATED(e_coulomb)
121 29042 : IF (atenergy) e_coulomb = 0._dp
122 :
123 28739 : atstress = PRESENT(pv_coulomb)
124 28739 : IF (atstress) atstress = ASSOCIATED(pv_coulomb)
125 32678 : IF (atstress) pv_coulomb = 0._dp
126 :
127 : ! pointing
128 28739 : CALL ewald_env_get(ewald_env, alpha=alpha, group=group)
129 28739 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, dg=dg)
130 28739 : CALL dg_get(dg, dg_rho0=dg_rho0)
131 28739 : rho0 => dg_rho0%density%array
132 28739 : pw_grid => pw_pool%pw_grid
133 28739 : bds => pw_grid%bounds
134 :
135 : ! allocating
136 28739 : nparticle_kind = SIZE(atomic_kind_set)
137 28739 : nnodes = 0
138 128109 : DO iparticle_kind = 1, nparticle_kind
139 128109 : nnodes = nnodes + local_particles%n_el(iparticle_kind)
140 : END DO
141 :
142 28739 : CALL structure_factor_allocate(pw_grid%bounds, nnodes, exp_igr)
143 :
144 86217 : ALLOCATE (summe(1:pw_grid%ngpts_cut))
145 85684 : ALLOCATE (charge(1:nnodes))
146 :
147 : ! Initializing vg_coulomb and fg_coulomb
148 28739 : vg_coulomb = 0.0_dp
149 1718291 : fg_coulomb = 0.0_dp
150 31989 : IF (use_virial) pv_g = 0.0_dp
151 : ! defining four_alpha_sq
152 28739 : four_alpha_sq = 4.0_dp*alpha**2
153 : ! zero node count
154 28739 : node = 0
155 128109 : DO iparticle_kind = 1, nparticle_kind
156 99370 : nparticle_local = local_particles%n_el(iparticle_kind)
157 128109 : IF (use_charge_array) THEN
158 894 : DO iparticle_local = 1, nparticle_local
159 402 : node = node + 1
160 402 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
161 402 : charge(node) = charges(iparticle)
162 5226 : vec = MATMUL(cell%h_inv, particle_set(iparticle)%r)
163 : CALL structure_factor_evaluate(vec, exp_igr%lb, &
164 894 : exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
165 : END DO
166 : ELSE
167 98878 : atomic_kind => atomic_kind_set(iparticle_kind)
168 98878 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q)
169 520864 : DO iparticle_local = 1, nparticle_local
170 421986 : node = node + 1
171 421986 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
172 421986 : charge(node) = q
173 5485818 : vec = MATMUL(cell%h_inv, particle_set(iparticle)%r)
174 : CALL structure_factor_evaluate(vec, exp_igr%lb, &
175 520864 : exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
176 : END DO
177 : END IF
178 : END DO
179 :
180 63964050 : summe(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
181 : ! looping over the positive g-vectors
182 63964050 : DO gpt = 1, pw_grid%ngpts_cut_local
183 :
184 63935311 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
185 63935311 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
186 63935311 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
187 :
188 63935311 : lp = lp + bds(1, 1)
189 63935311 : mp = mp + bds(1, 2)
190 63935311 : np = np + bds(1, 3)
191 :
192 : ! initializing sum to be used in the energy and force
193 872097294 : DO node = 1, nnodes
194 : summe(gpt) = summe(gpt) + charge(node)* &
195 : (exp_igr%ex(lp, node) &
196 : *exp_igr%ey(mp, node) &
197 872068555 : *exp_igr%ez(np, node))
198 : END DO
199 : END DO
200 28739 : CALL group%sum(summe)
201 :
202 28739 : pref = fourpi/pw_grid%vol
203 :
204 : ! looping over the positive g-vectors
205 63964050 : DO gpt = 1, pw_grid%ngpts_cut_local
206 : ! computing the potential energy
207 63935311 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
208 63935311 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
209 63935311 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
210 :
211 63935311 : lp = lp + bds(1, 1)
212 63935311 : mp = mp + bds(1, 2)
213 63935311 : np = np + bds(1, 3)
214 :
215 63935311 : IF (pw_grid%gsq(gpt) <= 1.0E-10_dp) CYCLE
216 :
217 63906572 : gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt)
218 63906572 : factor = gauss*REAL(summe(gpt)*CONJG(summe(gpt)), KIND=dp)
219 63906572 : vg_coulomb = vg_coulomb + factor
220 :
221 : ! atomic energies
222 63906572 : IF (atenergy) THEN
223 2893145 : DO node = 1, nnodes
224 : snode = CONJG(exp_igr%ex(lp, node) &
225 : *exp_igr%ey(mp, node) &
226 1735887 : *exp_igr%ez(np, node))
227 2893145 : e_coulomb(node) = e_coulomb(node) + gauss*charge(node)*REAL(summe(gpt)*snode, KIND=dp)
228 : END DO
229 : END IF
230 :
231 : ! computing the force
232 : node = 0
233 871617428 : DO node = 1, nnodes
234 : e_igdotr = AIMAG(summe(gpt)*CONJG &
235 : (exp_igr%ex(lp, node) &
236 : *exp_igr%ey(mp, node) &
237 807710856 : *exp_igr%ez(np, node)))
238 : fg_coulomb(:, node) = fg_coulomb(:, node) &
239 3294749996 : + charge(node)*gauss*e_igdotr*pw_grid%g(:, gpt)
240 : END DO
241 :
242 : ! compute the virial P*V
243 63906572 : denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
244 63906572 : IF (use_virial) THEN
245 1353328 : pv_g(1, 1) = pv_g(1, 1) + factor*(1.0_dp - 2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(1, gpt)*denom)
246 1353328 : pv_g(1, 2) = pv_g(1, 2) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
247 1353328 : pv_g(1, 3) = pv_g(1, 3) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
248 1353328 : pv_g(2, 1) = pv_g(2, 1) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
249 1353328 : pv_g(2, 2) = pv_g(2, 2) + factor*(1.0_dp - 2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(2, gpt)*denom)
250 1353328 : pv_g(2, 3) = pv_g(2, 3) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
251 1353328 : pv_g(3, 1) = pv_g(3, 1) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
252 1353328 : pv_g(3, 2) = pv_g(3, 2) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
253 1353328 : pv_g(3, 3) = pv_g(3, 3) + factor*(1.0_dp - 2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(3, gpt)*denom)
254 : END IF
255 63935311 : IF (atstress) THEN
256 2893145 : DO node = 1, nnodes
257 : snode = CONJG(exp_igr%ex(lp, node) &
258 : *exp_igr%ey(mp, node) &
259 1735887 : *exp_igr%ez(np, node))
260 1735887 : factor = gauss*charge(node)*REAL(summe(gpt)*snode, KIND=dp)
261 1735887 : pv_coulomb(1, 1, node) = pv_coulomb(1, 1, node) + factor*(1.0_dp - 2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(1, gpt)*denom)
262 1735887 : pv_coulomb(1, 2, node) = pv_coulomb(1, 2, node) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
263 1735887 : pv_coulomb(1, 3, node) = pv_coulomb(1, 3, node) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
264 1735887 : pv_coulomb(2, 1, node) = pv_coulomb(2, 1, node) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
265 1735887 : pv_coulomb(2, 2, node) = pv_coulomb(2, 2, node) + factor*(1.0_dp - 2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(2, gpt)*denom)
266 1735887 : pv_coulomb(2, 3, node) = pv_coulomb(2, 3, node) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
267 1735887 : pv_coulomb(3, 1, node) = pv_coulomb(3, 1, node) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
268 1735887 : pv_coulomb(3, 2, node) = pv_coulomb(3, 2, node) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
269 2893145 : pv_coulomb(3, 3, node) = pv_coulomb(3, 3, node) + factor*(1.0_dp - 2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(3, gpt)*denom)
270 : END DO
271 : END IF
272 : END DO
273 :
274 28739 : vg_coulomb = vg_coulomb*pref
275 31989 : IF (use_virial) pv_g = pv_g*pref
276 29042 : IF (atenergy) e_coulomb = e_coulomb*pref
277 32678 : IF (atstress) pv_coulomb = pv_coulomb*pref
278 :
279 1718291 : fg_coulomb = fg_coulomb*(2.0_dp*pref)
280 :
281 28739 : CALL structure_factor_deallocate(exp_igr)
282 :
283 28739 : DEALLOCATE (charge, summe)
284 :
285 28739 : CALL timestop(handle)
286 :
287 114956 : END SUBROUTINE ewald_evaluate
288 :
289 : ! **************************************************************************************************
290 : !> \brief Computes the self interaction from g-space
291 : !> and the neutralizing background
292 : !> \param ewald_env ...
293 : !> \param cell ...
294 : !> \param atomic_kind_set ...
295 : !> \param local_particles ...
296 : !> \param e_self ...
297 : !> \param e_neut ...
298 : !> \param charges ...
299 : !> \par History
300 : !> none
301 : !> \author CJM
302 : ! **************************************************************************************************
303 60041 : SUBROUTINE ewald_self(ewald_env, cell, atomic_kind_set, local_particles, e_self, &
304 : e_neut, charges)
305 :
306 : TYPE(ewald_environment_type), POINTER :: ewald_env
307 : TYPE(cell_type), POINTER :: cell
308 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
309 : TYPE(distribution_1d_type), POINTER :: local_particles
310 : REAL(KIND=dp), INTENT(OUT) :: e_self, e_neut
311 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
312 :
313 : INTEGER :: ewald_type, iparticle_kind, &
314 : nparticle_kind, nparticle_local
315 : LOGICAL :: is_shell
316 : REAL(KIND=dp) :: alpha, mm_radius, q, q_neutg, q_self, &
317 : q_sum, qcore, qshell
318 : TYPE(atomic_kind_type), POINTER :: atomic_kind
319 : TYPE(mp_comm_type) :: group
320 : TYPE(shell_kind_type), POINTER :: shell
321 :
322 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, &
323 60041 : alpha=alpha, group=group)
324 60041 : q_neutg = 0.0_dp
325 60041 : q_self = 0.0_dp
326 60041 : q_sum = 0.0_dp
327 60041 : nparticle_kind = SIZE(atomic_kind_set)
328 60041 : IF (ASSOCIATED(charges)) THEN
329 2644 : q_self = DOT_PRODUCT(charges, charges)
330 2644 : q_sum = SUM(charges)
331 : ! check and abort..
332 1928 : DO iparticle_kind = 1, nparticle_kind
333 1300 : atomic_kind => atomic_kind_set(iparticle_kind)
334 1300 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius)
335 1928 : IF (mm_radius > 0.0_dp) THEN
336 0 : CPABORT("Array of charges not implemented for mm_radius>0.0 !!")
337 : END IF
338 : END DO
339 : ELSE
340 263059 : DO iparticle_kind = 1, nparticle_kind
341 203646 : atomic_kind => atomic_kind_set(iparticle_kind)
342 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius, &
343 203646 : qeff=q, shell_active=is_shell, shell=shell)
344 203646 : nparticle_local = local_particles%n_el(iparticle_kind)
345 263059 : IF (is_shell) THEN
346 15954 : CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell)
347 : ! MI: the core-shell ES interaction, when core and shell belong to the same ion, is excluded
348 : ! in the nonbond correction term. Therefore, here the self interaction is computed entirely
349 15954 : q_self = q_self + qcore*qcore*nparticle_local + qshell*qshell*nparticle_local
350 15954 : q_sum = q_sum + qcore*nparticle_local + qshell*nparticle_local
351 15954 : IF (mm_radius > 0) THEN
352 : ! the core is always a point charge
353 0 : q_neutg = q_neutg + 2.0_dp*qshell*mm_radius**2
354 : END IF
355 : ELSE
356 187692 : q_self = q_self + q*q*nparticle_local
357 187692 : q_sum = q_sum + q*nparticle_local
358 187692 : IF (mm_radius > 0) THEN
359 0 : q_neutg = q_neutg + 2.0_dp*q*mm_radius**2
360 : END IF
361 : END IF
362 : END DO
363 :
364 59413 : CALL group%sum(q_self)
365 59413 : CALL group%sum(q_sum)
366 : END IF
367 :
368 60041 : e_neut = 0.0_dp
369 60041 : e_self = 0.0_dp
370 60041 : IF (ewald_type /= do_ewald_none) THEN
371 60041 : e_self = -q_self*alpha*oorootpi
372 60041 : e_neut = -q_sum*pi/(2.0_dp*cell%deth)*(q_sum/alpha**2 - q_neutg)
373 : END IF
374 :
375 60041 : END SUBROUTINE ewald_self
376 :
377 : ! **************************************************************************************************
378 : !> \brief Computes the self interaction per atom
379 : !> \param ewald_env ...
380 : !> \param atomic_kind_set ...
381 : !> \param local_particles ...
382 : !> \param e_self ...
383 : !> \param charges ...
384 : !> \par History
385 : !> none
386 : !> \author JHU from ewald_self
387 : ! **************************************************************************************************
388 636 : SUBROUTINE ewald_self_atom(ewald_env, atomic_kind_set, local_particles, e_self, &
389 : charges)
390 :
391 : TYPE(ewald_environment_type), POINTER :: ewald_env
392 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set(:)
393 : TYPE(distribution_1d_type), POINTER :: local_particles
394 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: e_self(:)
395 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
396 :
397 : INTEGER :: ewald_type, ii, iparticle_kind, &
398 : iparticle_local, nparticle_kind, &
399 : nparticle_local
400 : LOGICAL :: is_shell
401 : REAL(KIND=dp) :: alpha, fself, q, qcore, qshell
402 : TYPE(atomic_kind_type), POINTER :: atomic_kind
403 : TYPE(shell_kind_type), POINTER :: shell
404 :
405 636 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha)
406 :
407 636 : fself = alpha*oorootpi
408 :
409 636 : IF (ewald_type /= do_ewald_none) THEN
410 636 : nparticle_kind = SIZE(atomic_kind_set)
411 636 : IF (ASSOCIATED(charges)) THEN
412 0 : CPABORT("Atomic energy not implemented for charges")
413 : ELSE
414 1908 : DO iparticle_kind = 1, nparticle_kind
415 1272 : atomic_kind => atomic_kind_set(iparticle_kind)
416 1272 : nparticle_local = local_particles%n_el(iparticle_kind)
417 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, &
418 1272 : shell_active=is_shell, shell=shell)
419 1908 : IF (is_shell) THEN
420 30 : CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell)
421 1110 : DO iparticle_local = 1, nparticle_local
422 1080 : ii = local_particles%list(iparticle_kind)%array(iparticle_local)
423 1110 : e_self(ii) = e_self(ii) - (qcore*qcore + qshell*qshell)*fself
424 : END DO
425 : ELSE
426 2871 : DO iparticle_local = 1, nparticle_local
427 1629 : ii = local_particles%list(iparticle_kind)%array(iparticle_local)
428 2871 : e_self(ii) = e_self(ii) - q*q*fself
429 : END DO
430 : END IF
431 : END DO
432 : END IF
433 : END IF
434 :
435 636 : END SUBROUTINE ewald_self_atom
436 :
437 : ! **************************************************************************************************
438 : !> \brief ...
439 : !> \param iw ...
440 : !> \param pot_nonbond ...
441 : !> \param e_gspace ...
442 : !> \param e_self ...
443 : !> \param e_neut ...
444 : !> \param e_bonded ...
445 : !> \par History
446 : !> none
447 : !> \author CJM
448 : ! **************************************************************************************************
449 60041 : SUBROUTINE ewald_print(iw, pot_nonbond, e_gspace, e_self, e_neut, e_bonded)
450 :
451 : INTEGER, INTENT(IN) :: iw
452 : REAL(KIND=dp), INTENT(IN) :: pot_nonbond, e_gspace, e_self, e_neut, &
453 : e_bonded
454 :
455 60041 : IF (iw > 0) THEN
456 208 : WRITE (iw, '( A, A )') ' *********************************', &
457 416 : '**********************************************'
458 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL GSPACE ENERGY', &
459 416 : '[hartree]', '= ', e_gspace
460 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL NONBONDED ENERGY', &
461 416 : '[hartree]', '= ', pot_nonbond
462 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' SELF ENERGY CORRECTION', &
463 416 : '[hartree]', '= ', e_self
464 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' NEUT. BACKGROUND', &
465 416 : '[hartree]', '= ', e_neut
466 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' BONDED CORRECTION', &
467 416 : '[hartree]', '= ', e_bonded
468 208 : WRITE (iw, '( A, A )') ' *********************************', &
469 416 : '**********************************************'
470 : END IF
471 60041 : END SUBROUTINE ewald_print
472 :
473 : END MODULE ewalds
|
