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