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)
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
89 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ewald_evaluate'
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
111 CALL timeset(routinen, handle)
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
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
128 nparticle_kind =
SIZE(atomic_kind_set)
130 DO iparticle_kind = 1, nparticle_kind
131 nnodes = nnodes + local_particles%n_el(iparticle_kind)
136 ALLOCATE (summe(1:pw_grid%ngpts_cut))
137 ALLOCATE (charge(1:nnodes))
142 IF (use_virial) pv_g = 0.0_dp
144 four_alpha_sq = 4.0_dp*alpha**2
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
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)
156 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
159 atomic_kind => atomic_kind_set(iparticle_kind)
161 DO iparticle_local = 1, nparticle_local
163 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
165 vec = matmul(cell%h_inv, particle_set(iparticle)%r)
167 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
172 summe(:) = cmplx(0.0_dp, 0.0_dp, kind=
dp)
174 DO gpt = 1, pw_grid%ngpts_cut_local
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))
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))
192 CALL group%sum(summe)
197 DO gpt = 1, pw_grid%ngpts_cut_local
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))
207 IF (pw_grid%gsq(gpt) <= 1.0e-10_dp) cycle
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
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)
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)
235 denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
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)
249 vg_coulomb = vg_coulomb*pref
250 IF (use_virial) pv_g = pv_g*pref
251 IF (atenergy) e_coulomb = e_coulomb*pref
253 fg_coulomb = fg_coulomb*(2.0_dp*pref)
257 DEALLOCATE (charge, summe)
259 CALL timestop(handle)
277 SUBROUTINE ewald_self(ewald_env, cell, atomic_kind_set, local_particles, e_self, &
284 REAL(kind=
dp),
INTENT(OUT) :: e_self, e_neut
285 REAL(kind=
dp),
DIMENSION(:),
POINTER :: charges
287 INTEGER :: ewald_type, iparticle_kind, &
288 nparticle_kind, nparticle_local
290 REAL(kind=
dp) :: alpha, mm_radius, q, q_neutg, q_self, &
297 alpha=alpha, group=group)
301 nparticle_kind =
SIZE(atomic_kind_set)
302 IF (
ASSOCIATED(charges))
THEN
303 q_self = dot_product(charges, charges)
306 DO iparticle_kind = 1, nparticle_kind
307 atomic_kind => atomic_kind_set(iparticle_kind)
309 IF (mm_radius > 0.0_dp)
THEN
310 cpabort(
"Array of charges not implemented for mm_radius > 0.0")
314 DO iparticle_kind = 1, nparticle_kind
315 atomic_kind => atomic_kind_set(iparticle_kind)
317 qeff=q, shell_active=is_shell, shell=shell)
318 nparticle_local = local_particles%n_el(iparticle_kind)
320 CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell)
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
327 q_neutg = q_neutg + 2.0_dp*qshell*mm_radius**2
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
338 CALL group%sum(q_self)
339 CALL group%sum(q_sum)
346 e_neut = -q_sum*
pi/(2.0_dp*cell%deth)*(q_sum/alpha**2 - q_neutg)
368 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: e_self(:)
369 REAL(kind=
dp),
DIMENSION(:),
POINTER :: charges
371 INTEGER :: ewald_type, ii, iparticle_kind, &
372 iparticle_local, nparticle_kind, &
375 REAL(kind=
dp) :: alpha, fself, q, qcore, qshell
379 CALL ewald_env_get(ewald_env, ewald_type=ewald_type, alpha=alpha)
384 nparticle_kind =
SIZE(atomic_kind_set)
385 IF (
ASSOCIATED(charges))
THEN
386 cpabort(
"Atomic energy not implemented for charges")
388 DO iparticle_kind = 1, nparticle_kind
389 atomic_kind => atomic_kind_set(iparticle_kind)
390 nparticle_local = local_particles%n_el(iparticle_kind)
392 shell_active=is_shell, shell=shell)
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
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