76 SUBROUTINE ewald_evaluate(ewald_env, ewald_pw, cell, atomic_kind_set, particle_set, &
77 local_particles, fg_coulomb, vg_coulomb, pv_g, use_virial, charges, e_coulomb, &
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, &
93 CHARACTER(LEN=*),
PARAMETER :: routinen =
'ewald_evaluate'
95 COMPLEX(KIND=dp) :: snode
96 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 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 REAL(kind=
dp),
DIMENSION(:),
POINTER :: charge
106 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: rho0
115 CALL timeset(routinen, handle)
117 use_charge_array =
PRESENT(charges)
118 IF (use_charge_array) use_charge_array =
ASSOCIATED(charges)
119 atenergy =
PRESENT(e_coulomb)
120 IF (atenergy) atenergy =
ASSOCIATED(e_coulomb)
121 IF (atenergy) e_coulomb = 0._dp
123 atstress =
PRESENT(pv_coulomb)
124 IF (atstress) atstress =
ASSOCIATED(pv_coulomb)
125 IF (atstress) pv_coulomb = 0._dp
130 CALL dg_get(dg, dg_rho0=dg_rho0)
131 rho0 => dg_rho0%density%array
132 pw_grid => pw_pool%pw_grid
133 bds => pw_grid%bounds
136 nparticle_kind =
SIZE(atomic_kind_set)
138 DO iparticle_kind = 1, nparticle_kind
139 nnodes = nnodes + local_particles%n_el(iparticle_kind)
144 ALLOCATE (summe(1:pw_grid%ngpts_cut))
145 ALLOCATE (charge(1:nnodes))
150 IF (use_virial) pv_g = 0.0_dp
152 four_alpha_sq = 4.0_dp*alpha**2
155 DO iparticle_kind = 1, nparticle_kind
156 nparticle_local = local_particles%n_el(iparticle_kind)
157 IF (use_charge_array)
THEN
158 DO iparticle_local = 1, nparticle_local
160 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
161 charge(node) = charges(iparticle)
162 vec = matmul(cell%h_inv, particle_set(iparticle)%r)
164 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
167 atomic_kind => atomic_kind_set(iparticle_kind)
169 DO iparticle_local = 1, nparticle_local
171 iparticle = local_particles%list(iparticle_kind)%array(iparticle_local)
173 vec = matmul(cell%h_inv, particle_set(iparticle)%r)
175 exp_igr%ex(:, node), exp_igr%ey(:, node), exp_igr%ez(:, node))
180 summe(:) = cmplx(0.0_dp, 0.0_dp, kind=
dp)
182 DO gpt = 1, pw_grid%ngpts_cut_local
184 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
185 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
186 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
194 summe(gpt) = summe(gpt) + charge(node)* &
195 (exp_igr%ex(lp, node) &
196 *exp_igr%ey(mp, node) &
197 *exp_igr%ez(np, node))
200 CALL group%sum(summe)
205 DO gpt = 1, pw_grid%ngpts_cut_local
207 lp = pw_grid%mapl%pos(pw_grid%g_hat(1, gpt))
208 mp = pw_grid%mapm%pos(pw_grid%g_hat(2, gpt))
209 np = pw_grid%mapn%pos(pw_grid%g_hat(3, gpt))
215 IF (pw_grid%gsq(gpt) <= 1.0e-10_dp) cycle
217 gauss = (rho0(lp, mp, np)*pw_grid%vol)**2/pw_grid%gsq(gpt)
218 factor = gauss*real(summe(gpt)*conjg(summe(gpt)), kind=
dp)
219 vg_coulomb = vg_coulomb + factor
224 snode = conjg(exp_igr%ex(lp, node) &
225 *exp_igr%ey(mp, node) &
226 *exp_igr%ez(np, node))
227 e_coulomb(node) = e_coulomb(node) + gauss*charge(node)*real(summe(gpt)*snode, kind=
dp)
234 e_igdotr = aimag(summe(gpt)*conjg &
235 (exp_igr%ex(lp, node) &
236 *exp_igr%ey(mp, node) &
237 *exp_igr%ez(np, node)))
238 fg_coulomb(:, node) = fg_coulomb(:, node) &
239 + charge(node)*gauss*e_igdotr*pw_grid%g(:, gpt)
243 denom = 1.0_dp/four_alpha_sq + 1.0_dp/pw_grid%gsq(gpt)
245 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 pv_g(1, 2) = pv_g(1, 2) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(2, gpt)*denom)
247 pv_g(1, 3) = pv_g(1, 3) - factor*(2.0_dp*pw_grid%g(1, gpt)*pw_grid%g(3, gpt)*denom)
248 pv_g(2, 1) = pv_g(2, 1) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(1, gpt)*denom)
249 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 pv_g(2, 3) = pv_g(2, 3) - factor*(2.0_dp*pw_grid%g(2, gpt)*pw_grid%g(3, gpt)*denom)
251 pv_g(3, 1) = pv_g(3, 1) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(1, gpt)*denom)
252 pv_g(3, 2) = pv_g(3, 2) - factor*(2.0_dp*pw_grid%g(3, gpt)*pw_grid%g(2, gpt)*denom)
253 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)
257 snode = conjg(exp_igr%ex(lp, node) &
258 *exp_igr%ey(mp, node) &
259 *exp_igr%ez(np, node))
260 factor = gauss*charge(node)*real(summe(gpt)*snode, kind=
dp)
261 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 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 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 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 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 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 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 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 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)
274 vg_coulomb = vg_coulomb*pref
275 IF (use_virial) pv_g = pv_g*pref
276 IF (atenergy) e_coulomb = e_coulomb*pref
277 IF (atstress) pv_coulomb = pv_coulomb*pref
279 fg_coulomb = fg_coulomb*(2.0_dp*pref)
283 DEALLOCATE (charge, summe)
285 CALL timestop(handle)
303 SUBROUTINE ewald_self(ewald_env, cell, atomic_kind_set, local_particles, e_self, &
310 REAL(kind=
dp),
INTENT(OUT) :: e_self, e_neut
311 REAL(kind=
dp),
DIMENSION(:),
POINTER :: charges
313 INTEGER :: ewald_type, iparticle_kind, &
314 nparticle_kind, nparticle_local
316 REAL(kind=
dp) :: alpha, mm_radius, q, q_neutg, q_self, &
323 alpha=alpha, group=group)
327 nparticle_kind =
SIZE(atomic_kind_set)
328 IF (
ASSOCIATED(charges))
THEN
329 q_self = dot_product(charges, charges)
332 DO iparticle_kind = 1, nparticle_kind
333 atomic_kind => atomic_kind_set(iparticle_kind)
335 IF (mm_radius > 0.0_dp)
THEN
336 cpabort(
"Array of charges not implemented for mm_radius>0.0 !!")
340 DO iparticle_kind = 1, nparticle_kind
341 atomic_kind => atomic_kind_set(iparticle_kind)
343 qeff=q, shell_active=is_shell, shell=shell)
344 nparticle_local = local_particles%n_el(iparticle_kind)
346 CALL get_shell(shell=shell, charge_core=qcore, charge_shell=qshell)
349 q_self = q_self + qcore*qcore*nparticle_local + qshell*qshell*nparticle_local
350 q_sum = q_sum + qcore*nparticle_local + qshell*nparticle_local
351 IF (mm_radius > 0)
THEN
353 q_neutg = q_neutg + 2.0_dp*qshell*mm_radius**2
356 q_self = q_self + q*q*nparticle_local
357 q_sum = q_sum + q*nparticle_local
358 IF (mm_radius > 0)
THEN
359 q_neutg = q_neutg + 2.0_dp*q*mm_radius**2
364 CALL group%sum(q_self)
365 CALL group%sum(q_sum)
372 e_neut = -q_sum*
pi/(2.0_dp*cell%deth)*(q_sum/alpha**2 - q_neutg)