93 fg_coulomb, vg_coulomb, pv_g, shell_particle_set, core_particle_set, &
94 fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop)
100 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: fg_coulomb
101 REAL(kind=
dp),
INTENT(OUT) :: vg_coulomb
102 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: pv_g
104 POINTER :: shell_particle_set, core_particle_set
105 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT), &
106 OPTIONAL :: fgshell_coulomb, fgcore_coulomb
107 LOGICAL,
INTENT(IN) :: use_virial
108 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL,
POINTER :: charges
111 CHARACTER(len=*),
PARAMETER :: routinen =
'spme_evaluate'
113 INTEGER :: handle, i, ig, ipart, j, n, ncore, &
114 npart, nshell, o_spline, p1, p1_shell
115 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: center, core_center, shell_center
116 INTEGER,
DIMENSION(3) :: nd, npts
118 REAL(kind=
dp) :: alpha, dvols, fat1, ffa, ffb
119 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: core_delta, delta, shell_delta
120 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: rhos
121 REAL(kind=
dp),
DIMENSION(3) :: fat
122 REAL(kind=
dp),
DIMENSION(3, 3) :: f_stress, h_stress
135 NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, rden, rpot)
136 CALL timeset(routinen, handle)
140 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, &
142 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, &
143 poisson_env=poisson_env)
145 green => poisson_env%green_fft
146 grid_spme => pw_pool%pw_grid
148 npart =
SIZE(particle_set)
153 ALLOCATE (rhos(n, n, n))
159 ALLOCATE (center(3, npart), delta(3, npart))
160 CALL get_center(particle_set, box, center, delta, npts, n)
162 IF (
PRESENT(shell_particle_set) .AND. (
PRESENT(core_particle_set)))
THEN
163 cpassert(
ASSOCIATED(shell_particle_set))
164 cpassert(
ASSOCIATED(core_particle_set))
165 nshell =
SIZE(shell_particle_set)
166 ncore =
SIZE(core_particle_set)
167 cpassert(nshell == ncore)
168 ALLOCATE (shell_center(3, nshell), shell_delta(3, nshell))
169 CALL get_center(shell_particle_set, box, shell_center, shell_delta, npts, n)
170 ALLOCATE (core_center(3, nshell), core_delta(3, nshell))
171 CALL get_center(core_particle_set, box, core_center, core_delta, npts, n)
178 CALL set_list(particle_set, npart, center, p1, rden, ipart)
181 do_shell = (particle_set(p1)%shell_index /= 0)
184 CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.false., &
185 is_shell=.false., unit_charge=.false., charges=charges)
191 IF (
PRESENT(shell_particle_set) .AND.
PRESENT(core_particle_set))
THEN
195 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
197 IF (p1_shell == 0)
EXIT
198 CALL get_patch(shell_particle_set, shell_delta, green, p1_shell, rhos, &
199 is_core=.false., is_shell=.true., unit_charge=.false.)
202 CALL dg_sum_patch(rden, rhos, shell_center(:, p1_shell))
207 CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
209 IF (p1_shell == 0)
EXIT
210 CALL get_patch(core_particle_set, core_delta, green, p1_shell, rhos, &
211 is_core=.true., is_shell=.false., unit_charge=.false.)
221 CALL pw_pool%create_pw(rhob_r)
226 CALL pw_pool%create_pw(rhob_g)
234 CALL pw_pool%create_pw(dphi_g(i))
238 CALL pw_pool%create_pw(phi_g)
244 IF (atprop%energy .OR. atprop%stress)
THEN
253 CALL set_list(particle_set, npart, center, p1, rden, ipart)
255 do_shell = (particle_set(p1)%shell_index /= 0)
258 CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.false., &
259 is_shell=.false., unit_charge=.false., charges=charges)
262 IF (atprop%energy)
THEN
263 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
265 IF (atprop%stress)
THEN
266 atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*fat1*dvols
267 atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*fat1*dvols
268 atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*fat1*dvols
272 IF (
PRESENT(shell_particle_set) .AND.
PRESENT(core_particle_set))
THEN
275 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, rden, ipart)
276 IF (p1_shell == 0)
EXIT
277 CALL get_patch(shell_particle_set, shell_delta, green, p1_shell, rhos, &
278 is_core=.false., is_shell=.true., unit_charge=.false.)
280 p1 = shell_particle_set(p1_shell)%atom_index
281 IF (atprop%energy)
THEN
282 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
287 CALL set_list(core_particle_set, nshell, core_center, p1_shell, rden, ipart)
288 IF (p1_shell == 0)
EXIT
289 CALL get_patch(core_particle_set, core_delta, green, p1_shell, rhos, &
290 is_core=.true., is_shell=.false., unit_charge=.false.)
292 p1 = core_particle_set(p1_shell)%atom_index
293 IF (atprop%energy)
THEN
294 atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
298 IF (atprop%stress)
THEN
299 ffa = (0.5_dp/alpha)**2
302 DO ig = grid_spme%first_gne0, grid_spme%ngpts_cut_local
303 phi_g%array(ig) = ffb*dphi_g(i)%array(ig)*(ffa*grid_spme%gsq(ig) + 1.0_dp)
304 phi_g%array(ig) = phi_g%array(ig)*poisson_env%green_fft%influence_fn%array(ig)
306 IF (grid_spme%have_g0) phi_g%array(1) = 0.0_dp
318 CALL set_list(particle_set, npart, center, p1, rden, ipart)
321 CALL get_patch(particle_set, delta, green, p1, rhos, &
322 is_core=.false., is_shell=.false., unit_charge=.false., charges=charges)
325 atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fat1*dvols
326 IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fat1*dvols
336 CALL pw_pool%give_back_pw(phi_g)
337 CALL pw_pool%give_back_pw(rhob_g)
338 DEALLOCATE (phi_g, rhob_g)
345 f_stress(j, i) = f_stress(i, j)
348 ffa = (1.0_dp/
fourpi)*(0.5_dp/alpha)**2
349 f_stress = -ffa*f_stress
350 pv_g = h_stress + f_stress
360 CALL pw_pool%give_back_pw(dphi_g(i))
364 CALL pw_pool%give_back_pw(rhob_r)
372 CALL set_list(particle_set, npart, center, p1, rden, ipart)
375 do_shell = (particle_set(p1)%shell_index /= 0)
378 CALL get_patch(particle_set, delta, green, p1, rhos, is_core=.false., &
379 is_shell=.false., unit_charge=.false., charges=charges)
383 fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols
384 fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols
385 fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols
388 IF (
PRESENT(shell_particle_set) .AND. (
PRESENT(core_particle_set)))
THEN
389 IF (
PRESENT(fgshell_coulomb))
THEN
391 fgshell_coulomb = 0.0_dp
394 CALL set_list(shell_particle_set, nshell, shell_center, p1_shell, &
396 IF (p1_shell == 0)
EXIT
398 CALL get_patch(shell_particle_set, shell_delta, green, &
399 p1_shell, rhos, is_core=.false., is_shell=.true., unit_charge=.false.)
403 fgshell_coulomb(1, p1_shell) = fgshell_coulomb(1, p1_shell) - fat(1)*dvols
404 fgshell_coulomb(2, p1_shell) = fgshell_coulomb(2, p1_shell) - fat(2)*dvols
405 fgshell_coulomb(3, p1_shell) = fgshell_coulomb(3, p1_shell) - fat(3)*dvols
409 IF (
PRESENT(fgcore_coulomb))
THEN
411 fgcore_coulomb = 0.0_dp
414 CALL set_list(core_particle_set, nshell, core_center, p1_shell, &
416 IF (p1_shell == 0)
EXIT
418 CALL get_patch(core_particle_set, core_delta, green, &
419 p1_shell, rhos, is_core=.true., is_shell=.false., unit_charge=.false.)
423 fgcore_coulomb(1, p1_shell) = fgcore_coulomb(1, p1_shell) - fat(1)*dvols
424 fgcore_coulomb(2, p1_shell) = fgcore_coulomb(2, p1_shell) - fat(2)*dvols
425 fgcore_coulomb(3, p1_shell) = fgcore_coulomb(3, p1_shell) - fat(3)*dvols
440 DEALLOCATE (center, delta)
441 IF (
ALLOCATED(shell_center))
THEN
442 DEALLOCATE (shell_center, shell_delta)
444 IF (
ALLOCATED(core_center))
THEN
445 DEALLOCATE (core_center, core_delta)
447 CALL timestop(handle)
466 particle_set_b, potential)
471 TYPE(
particle_type),
DIMENSION(:),
INTENT(IN) :: particle_set_a
472 REAL(kind=
dp),
DIMENSION(:),
POINTER :: charges_a
473 TYPE(
particle_type),
DIMENSION(:),
INTENT(IN) :: particle_set_b
474 REAL(kind=
dp),
DIMENSION(:),
INTENT(INOUT) :: potential
476 CHARACTER(len=*),
PARAMETER :: routinen =
'spme_potential'
478 INTEGER :: handle, ipart, n, npart_a, npart_b, &
480 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: center
481 INTEGER,
DIMENSION(3) :: npts
482 REAL(kind=
dp) :: alpha, dvols, fat1
483 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: delta
484 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: rhos
495 NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, pw_pool, &
497 CALL timeset(routinen, handle)
501 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
502 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
504 green => poisson_env%green_fft
505 grid_spme => pw_pool%pw_grid
507 npart_a =
SIZE(particle_set_a)
508 npart_b =
SIZE(particle_set_b)
513 ALLOCATE (rhos(n, n, n))
519 ALLOCATE (center(3, npart_a), delta(3, npart_a))
520 CALL get_center(particle_set_a, box, center, delta, npts, n)
526 CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
530 CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
535 DEALLOCATE (center, delta)
540 CALL pw_pool%create_pw(rhob_r)
545 CALL pw_pool%create_pw(rhob_g)
554 CALL pw_pool%create_pw(phi_g)
559 ALLOCATE (center(3, npart_b), delta(3, npart_b))
560 CALL get_center(particle_set_b, box, center, delta, npts, n)
567 CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
570 CALL get_patch(particle_set_b, delta, green, p1, rhos, &
571 is_core=.false., is_shell=.false., unit_charge=.true.)
574 potential(p1) = potential(p1) + fat1*dvols
578 CALL pw_pool%give_back_pw(phi_g)
579 CALL pw_pool%give_back_pw(rhob_g)
580 CALL pw_pool%give_back_pw(rhob_r)
581 DEALLOCATE (phi_g, rhob_g, rhob_r)
586 DEALLOCATE (center, delta)
587 CALL timestop(handle)
606 SUBROUTINE spme_forces(ewald_env, ewald_pw, box, particle_set_a, charges_a, &
607 particle_set_b, charges_b, forces_b)
612 TYPE(
particle_type),
DIMENSION(:),
INTENT(IN) :: particle_set_a
613 REAL(kind=
dp),
DIMENSION(:),
POINTER :: charges_a
614 TYPE(
particle_type),
DIMENSION(:),
INTENT(IN) :: particle_set_b
615 REAL(kind=
dp),
DIMENSION(:),
POINTER :: charges_b
616 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(INOUT) :: forces_b
618 CHARACTER(len=*),
PARAMETER :: routinen =
'spme_forces'
620 INTEGER :: handle, i, ipart, n, npart_a, npart_b, &
622 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: center
623 INTEGER,
DIMENSION(3) :: npts
624 REAL(kind=
dp) :: alpha, dvols
625 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: delta
626 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: rhos
627 REAL(kind=
dp),
DIMENSION(3) :: fat
640 NULLIFY (grid_spme, green, poisson_env, phi_g, rhob_g, rhob_r, &
642 CALL timeset(routinen, handle)
646 CALL ewald_env_get(ewald_env, alpha=alpha, o_spline=o_spline, group=group)
647 CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_pool, rs_desc=rs_desc, poisson_env=poisson_env)
649 green => poisson_env%green_fft
650 grid_spme => pw_pool%pw_grid
652 npart_a =
SIZE(particle_set_a)
653 npart_b =
SIZE(particle_set_b)
658 ALLOCATE (rhos(n, n, n))
664 ALLOCATE (center(3, npart_a), delta(3, npart_a))
665 CALL get_center(particle_set_a, box, center, delta, npts, n)
671 CALL set_list(particle_set_a, npart_a, center, p1, rden, ipart)
675 CALL get_patch(particle_set_a, delta, green, p1, rhos, charges_a)
680 DEALLOCATE (center, delta)
685 CALL pw_pool%create_pw(rhob_r)
690 CALL pw_pool%create_pw(rhob_g)
698 CALL pw_pool%create_pw(dphi_g(i))
702 CALL pw_pool%create_pw(phi_g)
713 CALL pw_pool%give_back_pw(dphi_g(i))
717 ALLOCATE (center(3, npart_b), delta(3, npart_b))
718 CALL get_center(particle_set_b, box, center, delta, npts, n)
721 CALL set_list(particle_set_b, npart_b, center, p1, rden, ipart)
724 CALL get_patch(particle_set_b, delta, green, p1, rhos, &
725 is_core=.false., is_shell=.false., unit_charge=.false., charges=charges_b)
728 forces_b(1, p1) = forces_b(1, p1) - fat(1)*dvols
729 forces_b(2, p1) = forces_b(2, p1) - fat(2)*dvols
730 forces_b(3, p1) = forces_b(3, p1) - fat(3)*dvols
736 CALL pw_pool%give_back_pw(phi_g)
737 CALL pw_pool%give_back_pw(rhob_g)
738 CALL pw_pool%give_back_pw(rhob_r)
739 DEALLOCATE (phi_g, rhob_g, rhob_r)
744 DEALLOCATE (center, delta)
745 CALL timestop(handle)