Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \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 : !> \par History
72 : !> JGH (21-Feb-2001) : changed name
73 : !> \author CJM
74 : ! **************************************************************************************************
75 28801 : SUBROUTINE ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
76 28801 : local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial, charges, e_coulomb)
77 : TYPE(ewald_environment_type), POINTER :: ewald_env
78 : TYPE(ewald_pw_type), POINTER :: ewald_pw
79 : TYPE(cell_type), POINTER :: cell
80 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
81 : TYPE(particle_type), POINTER :: particle_set(:)
82 : TYPE(distribution_1d_type), POINTER :: local_particles
83 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: fg_coulomb
84 : REAL(KIND=dp), INTENT(OUT) :: vg_coulomb
85 : REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT) :: pv_g
86 : LOGICAL, INTENT(IN) :: use_virial
87 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: charges, e_coulomb
88 :
89 : CHARACTER(LEN=*), PARAMETER :: routineN = 'ewald_evaluate'
90 :
91 : COMPLEX(KIND=dp) :: snode
92 : COMPLEX(KIND=dp), ALLOCATABLE, DIMENSION(:) :: summe
93 : INTEGER :: gpt, handle, iparticle, iparticle_kind, &
94 : iparticle_local, lp, mp, nnodes, node, &
95 : np, nparticle_kind, nparticle_local
96 28801 : INTEGER, DIMENSION(:, :), POINTER :: bds
97 : LOGICAL :: atenergy, use_charge_array
98 : REAL(KIND=dp) :: alpha, denom, e_igdotr, factor, &
99 : four_alpha_sq, gauss, pref, q
100 : REAL(KIND=dp), DIMENSION(3) :: vec
101 28801 : REAL(KIND=dp), DIMENSION(:), POINTER :: charge
102 28801 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: rho0
103 : TYPE(atomic_kind_type), POINTER :: atomic_kind
104 : TYPE(dg_rho0_type), POINTER :: dg_rho0
105 : TYPE(dg_type), POINTER :: dg
106 : TYPE(mp_comm_type) :: group
107 : TYPE(pw_grid_type), POINTER :: pw_grid
108 : TYPE(pw_pool_type), POINTER :: pw_pool
109 : TYPE(structure_factor_type) :: exp_igr
110 :
111 28801 : CALL timeset(routineN, handle)
112 28801 : CALL cite_reference(Ewald1921)
113 28801 : use_charge_array = .FALSE.
114 28801 : IF (PRESENT(charges)) use_charge_array = ASSOCIATED(charges)
115 28801 : atenergy = PRESENT(e_coulomb)
116 28801 : IF (atenergy) atenergy = ASSOCIATED(e_coulomb)
117 29104 : IF (atenergy) e_coulomb = 0._dp
118 :
119 : ! pointing
120 28801 : CALL ewald_env_get(ewald_env, alpha=alpha, group=group)
121 28801 : CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, dg=dg)
122 28801 : CALL dg_get(dg, dg_rho0=dg_rho0)
123 28801 : rho0 => dg_rho0%density%array
124 28801 : pw_grid => pw_pool%pw_grid
125 28801 : bds => pw_grid%bounds
126 :
127 : ! allocating
128 28801 : nparticle_kind = SIZE(atomic_kind_set)
129 28801 : nnodes = 0
130 128295 : DO iparticle_kind = 1, nparticle_kind
131 128295 : nnodes = nnodes + local_particles%n_el(iparticle_kind)
132 : END DO
133 :
134 28801 : CALL structure_factor_allocate(pw_grid%bounds, nnodes, exp_igr)
135 :
136 86403 : ALLOCATE (summe(1:pw_grid%ngpts_cut))
137 85870 : ALLOCATE (charge(1:nnodes))
138 :
139 : ! Initializing vg_coulomb and fg_coulomb
140 28801 : vg_coulomb = 0.0_dp
141 1723561 : fg_coulomb = 0.0_dp
142 32051 : IF (use_virial) pv_g = 0.0_dp
143 : ! defining four_alpha_sq
144 28801 : four_alpha_sq = 4.0_dp*alpha**2
145 : ! zero node count
146 28801 : node = 0
147 128295 : DO iparticle_kind = 1, nparticle_kind
148 99494 : nparticle_local = local_particles%n_el(iparticle_kind)
149 128295 : IF (use_charge_array) THEN
150 894 : DO iparticle_local = 1, nparticle_local
151 402 : node = node + 1
152 402 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
153 402 : charge(node) = charges(iparticle)
154 5226 : vec = MATMUL(cell%h_inv, particle_set(iparticle)%r)
155 : CALL structure_factor_evaluate(vec, exp_igr%lb, &
156 894 : exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
157 : END DO
158 : ELSE
159 99002 : atomic_kind => atomic_kind_set(iparticle_kind)
160 99002 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q)
161 522290 : DO iparticle_local = 1, nparticle_local
162 423288 : node = node + 1
163 423288 : iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
164 423288 : charge(node) = q
165 5502744 : vec = MATMUL(cell%h_inv, particle_set(iparticle)%r)
166 : CALL structure_factor_evaluate(vec, exp_igr%lb, &
167 522290 : exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
168 : END DO
169 : END IF
170 : END DO
171 :
172 63979178 : summe(:) = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
173 : ! looping over the positive g-vectors
174 63979178 : DO gpt = 1, pw_grid%ngpts_cut_local
175 :
176 63950377 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
177 63950377 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
178 63950377 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
179 :
180 63950377 : lp = lp + bds(1, 1)
181 63950377 : mp = mp + bds(1, 2)
182 63950377 : np = np + bds(1, 3)
183 :
184 : ! initializing sum to be used in the energy and force
185 872428808 : DO node = 1, nnodes
186 : summe(gpt) = summe(gpt) + charge(node)* &
187 : (exp_igr%ex(lp, node) &
188 : *exp_igr%ey(mp, node) &
189 872400007 : *exp_igr%ez(np, node))
190 : END DO
191 : END DO
192 28801 : CALL group%sum(summe)
193 :
194 28801 : pref = fourpi/pw_grid%vol
195 :
196 : ! looping over the positive g-vectors
197 63979178 : DO gpt = 1, pw_grid%ngpts_cut_local
198 : ! computing the potential energy
199 63950377 : lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
200 63950377 : mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
201 63950377 : np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
202 :
203 63950377 : lp = lp + bds(1, 1)
204 63950377 : mp = mp + bds(1, 2)
205 63950377 : np = np + bds(1, 3)
206 :
207 63950377 : IF (pw_grid%gsq(gpt) <= 1.0E-10_dp) CYCLE
208 :
209 63921576 : gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt)
210 63921576 : factor = gauss*REAL(summe(gpt)*CONJG(summe(gpt)), KIND=dp)
211 63921576 : vg_coulomb = vg_coulomb + factor
212 :
213 : ! atomic energies
214 63921576 : IF (atenergy) THEN
215 2893145 : DO node = 1, nnodes
216 : snode = CONJG(exp_igr%ex(lp, node) &
217 : *exp_igr%ey(mp, node) &
218 1735887 : *exp_igr%ez(np, node))
219 2893145 : e_coulomb(node) = e_coulomb(node) + gauss*charge(node)*REAL(summe(gpt)*snode, KIND=dp)
220 : END DO
221 : END IF
222 :
223 : ! computing the force
224 : node = 0
225 871947516 : DO node = 1, nnodes
226 : e_igdotr = AIMAG(summe(gpt)*CONJG &
227 : (exp_igr%ex(lp, node) &
228 : *exp_igr%ey(mp, node) &
229 808025940 : *exp_igr%ez(np, node)))
230 : fg_coulomb(:, node) = fg_coulomb(:, node) &
231 3296025336 : + charge(node)*gauss*e_igdotr*pw_grid%g(:, gpt)
232 : END DO
233 :
234 : ! compute the virial P*V
235 63921576 : denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
236 63950377 : IF (use_virial) THEN
237 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)
238 1353328 : pv_g(1, 2) = pv_g(1, 2) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
239 1353328 : pv_g(1, 3) = pv_g(1, 3) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
240 1353328 : pv_g(2, 1) = pv_g(2, 1) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
241 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)
242 1353328 : pv_g(2, 3) = pv_g(2, 3) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
243 1353328 : pv_g(3, 1) = pv_g(3, 1) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
244 1353328 : pv_g(3, 2) = pv_g(3, 2) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
245 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)
246 : END IF
247 : END DO
248 :
249 28801 : vg_coulomb = vg_coulomb*pref
250 32051 : IF (use_virial) pv_g = pv_g*pref
251 29104 : IF (atenergy) e_coulomb = e_coulomb*pref
252 :
253 1723561 : fg_coulomb = fg_coulomb*(2.0_dp*pref)
254 :
255 28801 : CALL structure_factor_deallocate(exp_igr)
256 :
257 28801 : DEALLOCATE (charge, summe)
258 :
259 28801 : CALL timestop(handle)
260 :
261 144005 : END SUBROUTINE ewald_evaluate
262 :
263 : ! **************************************************************************************************
264 : !> \brief Computes the self interaction from g-space
265 : !> and the neutralizing background
266 : !> \param ewald_env ...
267 : !> \param cell ...
268 : !> \param atomic_kind_set ...
269 : !> \param local_particles ...
270 : !> \param e_self ...
271 : !> \param e_neut ...
272 : !> \param charges ...
273 : !> \par History
274 : !> none
275 : !> \author CJM
276 : ! **************************************************************************************************
277 60275 : SUBROUTINE ewald_self(ewald_env, cell, atomic_kind_set, local_particles, e_self, &
278 : e_neut, charges)
279 :
280 : TYPE(ewald_environment_type), POINTER :: ewald_env
281 : TYPE(cell_type), POINTER :: cell
282 : TYPE(atomic_kind_type), POINTER :: atomic_kind_set(:)
283 : TYPE(distribution_1d_type), POINTER :: local_particles
284 : REAL(KIND=dp), INTENT(OUT) :: e_self, e_neut
285 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
286 :
287 : INTEGER :: ewald_type, iparticle_kind, &
288 : nparticle_kind, nparticle_local
289 : LOGICAL :: is_shell
290 : REAL(KIND=dp) :: alpha, mm_radius, q, q_neutg, q_self, &
291 : q_sum, qcore, qshell
292 : TYPE(atomic_kind_type), POINTER :: atomic_kind
293 : TYPE(mp_comm_type) :: group
294 : TYPE(shell_kind_type), POINTER :: shell
295 :
296 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, &
297 60275 : alpha=alpha, group=group)
298 60275 : q_neutg = 0.0_dp
299 60275 : q_self = 0.0_dp
300 60275 : q_sum = 0.0_dp
301 60275 : nparticle_kind = SIZE(atomic_kind_set)
302 60275 : IF (ASSOCIATED(charges)) THEN
303 2644 : q_self = DOT_PRODUCT(charges, charges)
304 2644 : q_sum = SUM(charges)
305 : ! check and abort..
306 1928 : DO iparticle_kind = 1, nparticle_kind
307 1300 : atomic_kind => atomic_kind_set(iparticle_kind)
308 1300 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius)
309 1928 : IF (mm_radius > 0.0_dp) THEN
310 0 : CPABORT("Array of charges not implemented for mm_radius > 0.0")
311 : END IF
312 : END DO
313 : ELSE
314 263761 : DO iparticle_kind = 1, nparticle_kind
315 204114 : atomic_kind => atomic_kind_set(iparticle_kind)
316 : CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius, &
317 204114 : qeff=q, shell_active=is_shell, shell=shell)
318 204114 : nparticle_local = local_particles%n_el(iparticle_kind)
319 263761 : IF (is_shell) THEN
320 15956 : CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell)
321 : ! MI: the core-shell ES interaction, when core and shell belong to the same ion, is excluded
322 : ! in the nonbond correction term. Therefore, here the self interaction is computed entirely
323 15956 : q_self = q_self + qcore*qcore*nparticle_local + qshell*qshell*nparticle_local
324 15956 : q_sum = q_sum + qcore*nparticle_local + qshell*nparticle_local
325 15956 : IF (mm_radius > 0) THEN
326 : ! the core is always a point charge
327 0 : q_neutg = q_neutg + 2.0_dp*qshell*mm_radius**2
328 : END IF
329 : ELSE
330 188158 : q_self = q_self + q*q*nparticle_local
331 188158 : q_sum = q_sum + q*nparticle_local
332 188158 : IF (mm_radius > 0) THEN
333 0 : q_neutg = q_neutg + 2.0_dp*q*mm_radius**2
334 : END IF
335 : END IF
336 : END DO
337 :
338 59647 : CALL group%sum(q_self)
339 59647 : CALL group%sum(q_sum)
340 : END IF
341 :
342 60275 : e_neut = 0.0_dp
343 60275 : e_self = 0.0_dp
344 60275 : IF (ewald_type /= do_ewald_none) THEN
345 60275 : e_self = -q_self*alpha*oorootpi
346 60275 : e_neut = -q_sum*pi/(2.0_dp*cell%deth)*(q_sum/alpha**2 - q_neutg)
347 : END IF
348 :
349 60275 : END SUBROUTINE ewald_self
350 :
351 : ! **************************************************************************************************
352 : !> \brief Computes the self interaction per atom
353 : !> \param ewald_env ...
354 : !> \param atomic_kind_set ...
355 : !> \param local_particles ...
356 : !> \param e_self ...
357 : !> \param charges ...
358 : !> \par History
359 : !> none
360 : !> \author JHU from ewald_self
361 : ! **************************************************************************************************
362 636 : SUBROUTINE ewald_self_atom(ewald_env, atomic_kind_set, local_particles, e_self, &
363 : charges)
364 :
365 : TYPE(ewald_environment_type), POINTER :: ewald_env
366 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set(:)
367 : TYPE(distribution_1d_type), POINTER :: local_particles
368 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: e_self(:)
369 : REAL(KIND=dp), DIMENSION(:), POINTER :: charges
370 :
371 : INTEGER :: ewald_type, ii, iparticle_kind, &
372 : iparticle_local, nparticle_kind, &
373 : nparticle_local
374 : LOGICAL :: is_shell
375 : REAL(KIND=dp) :: alpha, fself, q, qcore, qshell
376 : TYPE(atomic_kind_type), POINTER :: atomic_kind
377 : TYPE(shell_kind_type), POINTER :: shell
378 :
379 636 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha)
380 :
381 636 : fself = alpha*oorootpi
382 :
383 636 : IF (ewald_type /= do_ewald_none) THEN
384 636 : nparticle_kind = SIZE(atomic_kind_set)
385 636 : IF (ASSOCIATED(charges)) THEN
386 0 : CPABORT("Atomic energy not implemented for charges")
387 : ELSE
388 1908 : DO iparticle_kind = 1, nparticle_kind
389 1272 : atomic_kind => atomic_kind_set(iparticle_kind)
390 1272 : nparticle_local = local_particles%n_el(iparticle_kind)
391 : CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, &
392 1272 : shell_active=is_shell, shell=shell)
393 1908 : IF (is_shell) THEN
394 30 : CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell)
395 1110 : DO iparticle_local = 1, nparticle_local
396 1080 : ii = local_particles%list(iparticle_kind)%array(iparticle_local)
397 1110 : e_self(ii) = e_self(ii) - (qcore*qcore + qshell*qshell)*fself
398 : END DO
399 : ELSE
400 2871 : DO iparticle_local = 1, nparticle_local
401 1629 : ii = local_particles%list(iparticle_kind)%array(iparticle_local)
402 2871 : e_self(ii) = e_self(ii) - q*q*fself
403 : END DO
404 : END IF
405 : END DO
406 : END IF
407 : END IF
408 :
409 636 : END SUBROUTINE ewald_self_atom
410 :
411 : ! **************************************************************************************************
412 : !> \brief ...
413 : !> \param iw ...
414 : !> \param pot_nonbond ...
415 : !> \param e_gspace ...
416 : !> \param e_self ...
417 : !> \param e_neut ...
418 : !> \param e_bonded ...
419 : !> \par History
420 : !> none
421 : !> \author CJM
422 : ! **************************************************************************************************
423 60275 : SUBROUTINE ewald_print(iw, pot_nonbond, e_gspace, e_self, e_neut, e_bonded)
424 :
425 : INTEGER, INTENT(IN) :: iw
426 : REAL(KIND=dp), INTENT(IN) :: pot_nonbond, e_gspace, e_self, e_neut, &
427 : e_bonded
428 :
429 60275 : IF (iw > 0) THEN
430 208 : WRITE (iw, '( A, A )') ' *********************************', &
431 416 : '**********************************************'
432 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL GSPACE ENERGY', &
433 416 : '[hartree]', '= ', e_gspace
434 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL NONBONDED ENERGY', &
435 416 : '[hartree]', '= ', pot_nonbond
436 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' SELF ENERGY CORRECTION', &
437 416 : '[hartree]', '= ', e_self
438 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' NEUT. BACKGROUND', &
439 416 : '[hartree]', '= ', e_neut
440 208 : WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' BONDED CORRECTION', &
441 416 : '[hartree]', '= ', e_bonded
442 208 : WRITE (iw, '( A, A )') ' *********************************', &
443 416 : '**********************************************'
444 : END IF
445 60275 : END SUBROUTINE ewald_print
446 :
447 : END MODULE ewalds
|
