Loading [MathJax]/extensions/tex2jax.js
 (git:aabdcc8)
All Data Structures Namespaces Files Functions Variables Typedefs Enumerations Enumerator Macros
ewalds.F
Go to the documentation of this file.
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! **************************************************************************************************
13MODULE ewalds
14
17 USE bibliography, ONLY: ewald1921,&
18 cite_reference
19 USE cell_types, ONLY: cell_type
21 USE dg_types, ONLY: dg_get,&
26 USE ewald_pw_types, ONLY: ewald_pw_get,&
28 USE kinds, ONLY: dp
29 USE mathconstants, ONLY: fourpi,&
30 oorootpi,&
31 pi
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
51
52CONTAINS
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 SUBROUTINE ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
76 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 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 REAL(kind=dp), DIMENSION(:), POINTER :: charge
102 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 CALL timeset(routinen, handle)
112 CALL cite_reference(ewald1921)
113 use_charge_array = .false.
114 IF (PRESENT(charges)) use_charge_array = ASSOCIATED(charges)
115 atenergy = PRESENT(e_coulomb)
116 IF (atenergy) atenergy = ASSOCIATED(e_coulomb)
117 IF (atenergy) e_coulomb = 0._dp
118
119 ! pointing
120 CALL ewald_env_get(ewald_env, alpha=alpha, group=group)
121 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, dg=dg)
122 CALL dg_get(dg, dg_rho0=dg_rho0)
123 rho0 => dg_rho0%density%array
124 pw_grid => pw_pool%pw_grid
125 bds => pw_grid%bounds
126
127 ! allocating
128 nparticle_kind = SIZE(atomic_kind_set)
129 nnodes = 0
130 DO iparticle_kind = 1, nparticle_kind
131 nnodes = nnodes + local_particles%n_el(iparticle_kind)
132 END DO
133
134 CALL structure_factor_allocate(pw_grid%bounds, nnodes, exp_igr)
135
136 ALLOCATE (summe(1:pw_grid%ngpts_cut))
137 ALLOCATE (charge(1:nnodes))
138
139 ! Initializing vg_coulomb and fg_coulomb
140 vg_coulomb = 0.0_dp
141 fg_coulomb = 0.0_dp
142 IF (use_virial) pv_g = 0.0_dp
143 ! defining four_alpha_sq
144 four_alpha_sq = 4.0_dp*alpha**2
145 ! zero node count
146 node = 0
147 DO iparticle_kind = 1, nparticle_kind
148 nparticle_local = local_particles%n_el(iparticle_kind)
149 IF (use_charge_array) THEN
150 DO iparticle_local = 1, nparticle_local
151 node = node + 1
152 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
153 charge(node) = charges(iparticle)
154 vec = matmul(cell%h_inv, particle_set(iparticle)%r)
155 CALL structure_factor_evaluate(vec, exp_igr%lb, &
156 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
157 END DO
158 ELSE
159 atomic_kind => atomic_kind_set(iparticle_kind)
160 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q)
161 DO iparticle_local = 1, nparticle_local
162 node = node + 1
163 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
164 charge(node) = q
165 vec = matmul(cell%h_inv, particle_set(iparticle)%r)
166 CALL structure_factor_evaluate(vec, exp_igr%lb, &
167 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
168 END DO
169 END IF
170 END DO
171
172 summe(:) = cmplx(0.0_dp, 0.0_dp, kind=dp)
173 ! looping over the positive g-vectors
174 DO gpt = 1, pw_grid%ngpts_cut_local
175
176 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
177 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
178 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
179
180 lp = lp + bds(1, 1)
181 mp = mp + bds(1, 2)
182 np = np + bds(1, 3)
183
184 ! initializing sum to be used in the energy and force
185 DO node = 1, nnodes
186 summe(gpt) = summe(gpt) + charge(node)* &
187 (exp_igr%ex(lp, node) &
188 *exp_igr%ey(mp, node) &
189 *exp_igr%ez(np, node))
190 END DO
191 END DO
192 CALL group%sum(summe)
193
194 pref = fourpi/pw_grid%vol
195
196 ! looping over the positive g-vectors
197 DO gpt = 1, pw_grid%ngpts_cut_local
198 ! computing the potential energy
199 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
200 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
201 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
202
203 lp = lp + bds(1, 1)
204 mp = mp + bds(1, 2)
205 np = np + bds(1, 3)
206
207 IF (pw_grid%gsq(gpt) <= 1.0e-10_dp) cycle
208
209 gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt)
210 factor = gauss*real(summe(gpt)*conjg(summe(gpt)), kind=dp)
211 vg_coulomb = vg_coulomb + factor
212
213 ! atomic energies
214 IF (atenergy) THEN
215 DO node = 1, nnodes
216 snode = conjg(exp_igr%ex(lp, node) &
217 *exp_igr%ey(mp, node) &
218 *exp_igr%ez(np, node))
219 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 DO node = 1, nnodes
226 e_igdotr = aimag(summe(gpt)*conjg &
227 (exp_igr%ex(lp, node) &
228 *exp_igr%ey(mp, node) &
229 *exp_igr%ez(np, node)))
230 fg_coulomb(:, node) = fg_coulomb(:, node) &
231 + charge(node)*gauss*e_igdotr*pw_grid%g(:, gpt)
232 END DO
233
234 ! compute the virial P*V
235 denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
236 IF (use_virial) THEN
237 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 pv_g(1, 2) = pv_g(1, 2) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
239 pv_g(1, 3) = pv_g(1, 3) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
240 pv_g(2, 1) = pv_g(2, 1) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
241 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 pv_g(2, 3) = pv_g(2, 3) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
243 pv_g(3, 1) = pv_g(3, 1) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
244 pv_g(3, 2) = pv_g(3, 2) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
245 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 vg_coulomb = vg_coulomb*pref
250 IF (use_virial) pv_g = pv_g*pref
251 IF (atenergy) e_coulomb = e_coulomb*pref
252
253 fg_coulomb = fg_coulomb*(2.0_dp*pref)
254
255 CALL structure_factor_deallocate(exp_igr)
256
257 DEALLOCATE (charge, summe)
258
259 CALL timestop(handle)
260
261 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 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 alpha=alpha, group=group)
298 q_neutg = 0.0_dp
299 q_self = 0.0_dp
300 q_sum = 0.0_dp
301 nparticle_kind = SIZE(atomic_kind_set)
302 IF (ASSOCIATED(charges)) THEN
303 q_self = dot_product(charges, charges)
304 q_sum = sum(charges)
305 ! check and abort..
306 DO iparticle_kind = 1, nparticle_kind
307 atomic_kind => atomic_kind_set(iparticle_kind)
308 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius)
309 IF (mm_radius > 0.0_dp) THEN
310 cpabort("Array of charges not implemented for mm_radius > 0.0")
311 END IF
312 END DO
313 ELSE
314 DO iparticle_kind = 1, nparticle_kind
315 atomic_kind => atomic_kind_set(iparticle_kind)
316 CALL get_atomic_kind(atomic_kind=atomic_kind, mm_radius=mm_radius, &
317 qeff=q, shell_active=is_shell, shell=shell)
318 nparticle_local = local_particles%n_el(iparticle_kind)
319 IF (is_shell) THEN
320 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 q_self = q_self + qcore*qcore*nparticle_local + qshell*qshell*nparticle_local
324 q_sum = q_sum + qcore*nparticle_local + qshell*nparticle_local
325 IF (mm_radius > 0) THEN
326 ! the core is always a point charge
327 q_neutg = q_neutg + 2.0_dp*qshell*mm_radius**2
328 END IF
329 ELSE
330 q_self = q_self + q*q*nparticle_local
331 q_sum = q_sum + q*nparticle_local
332 IF (mm_radius > 0) THEN
333 q_neutg = q_neutg + 2.0_dp*q*mm_radius**2
334 END IF
335 END IF
336 END DO
337
338 CALL group%sum(q_self)
339 CALL group%sum(q_sum)
340 END IF
341
342 e_neut = 0.0_dp
343 e_self = 0.0_dp
344 IF (ewald_type /= do_ewald_none) THEN
345 e_self = -q_self*alpha*oorootpi
346 e_neut = -q_sum*pi/(2.0_dp*cell%deth)*(q_sum/alpha**2 - q_neutg)
347 END IF
348
349 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 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 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha)
380
381 fself = alpha*oorootpi
382
383 IF (ewald_type /= do_ewald_none) THEN
384 nparticle_kind = SIZE(atomic_kind_set)
385 IF (ASSOCIATED(charges)) THEN
386 cpabort("Atomic energy not implemented for charges")
387 ELSE
388 DO iparticle_kind = 1, nparticle_kind
389 atomic_kind => atomic_kind_set(iparticle_kind)
390 nparticle_local = local_particles%n_el(iparticle_kind)
391 CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q, &
392 shell_active=is_shell, shell=shell)
393 IF (is_shell) THEN
394 CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell)
395 DO iparticle_local = 1, nparticle_local
396 ii = local_particles%list(iparticle_kind)%array(iparticle_local)
397 e_self(ii) = e_self(ii) - (qcore*qcore + qshell*qshell)*fself
398 END DO
399 ELSE
400 DO iparticle_local = 1, nparticle_local
401 ii = local_particles%list(iparticle_kind)%array(iparticle_local)
402 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 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 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 IF (iw > 0) THEN
430 WRITE (iw, '( A, A )') ' *********************************', &
431 '**********************************************'
432 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL GSPACE ENERGY', &
433 '[hartree]', '= ', e_gspace
434 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' INITIAL NONBONDED ENERGY', &
435 '[hartree]', '= ', pot_nonbond
436 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' SELF ENERGY CORRECTION', &
437 '[hartree]', '= ', e_self
438 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' NEUT. BACKGROUND', &
439 '[hartree]', '= ', e_neut
440 WRITE (iw, '( A, A, T35, A, T56, E25.15 )') ' BONDED CORRECTION', &
441 '[hartree]', '= ', e_bonded
442 WRITE (iw, '( A, A )') ' *********************************', &
443 '**********************************************'
444 END IF
445 END SUBROUTINE ewald_print
446
447END MODULE ewalds
Define the atomic kind types and their sub types.
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
collects all references to literature in CP2K as new algorithms / method are included from literature...
integer, save, public ewald1921
Handles all functions related to the CELL.
Definition cell_types.F:15
subroutine, public dg_get(dg, dg_rho0)
Get the dg_type.
Definition dg_types.F:44
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
subroutine, public ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg)
get the ewald_pw environment to the correct program.
subroutine, public ewald_self_atom(ewald_env, atomic_kind_set, local_particles, e_self, charges)
Computes the self interaction per atom.
Definition ewalds.F:364
subroutine, public ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial, charges, e_coulomb)
computes the potential and the force from the g-space part of the 1/r potential Ref....
Definition ewalds.F:77
subroutine, public ewald_print(iw, pot_nonbond, e_gspace, e_self, e_neut, e_bonded)
...
Definition ewalds.F:424
subroutine, public ewald_self(ewald_env, cell, atomic_kind_set, local_particles, e_self, e_neut, charges)
Computes the self interaction from g-space and the neutralizing background.
Definition ewalds.F:279
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Definition of mathematical constants and functions.
real(kind=dp), parameter, public oorootpi
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public fourpi
Interface to the message passing library MPI.
Define the data structure for the particle information.
functions related to the poisson solver on regular grids
integer, parameter, public do_ewald_none
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
elemental subroutine, public get_shell(shell, charge, charge_core, charge_shell, mass_core, mass_shell, k2_spring, k4_spring, max_dist, shell_cutoff)
...
subroutine, public structure_factor_deallocate(exp_igr)
...
subroutine, public structure_factor_allocate(bds, nparts, exp_igr, allocate_centre, allocate_shell_e, allocate_shell_centre, nshell)
...
subroutine, public structure_factor_evaluate(delta, lb, ex, ey, ez)
...
Provides all information about an atomic kind.
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
Type for Gaussian Densities type = type of gaussian (PME) grid = grid number gcc = Gaussian contracti...
structure to store local (to a processor) ordered lists of integers.
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...