LCOV - code coverage report
Current view: top level - src - pme.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 248 344 72.1 %
Date: 2024-03-28 07:31:50 Functions: 3 3 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \par History
      10             : !>      JGH (21-Mar-2001) : Complete rewrite
      11             : !> \author CJM and APSI
      12             : ! **************************************************************************************************
      13             : MODULE pme
      14             : 
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      16             :                                               get_atomic_kind
      17             :    USE atprop_types,                    ONLY: atprop_type
      18             :    USE bibliography,                    ONLY: cite_reference,&
      19             :                                               darden1993
      20             :    USE cell_types,                      ONLY: cell_type
      21             :    USE dg_rho0_types,                   ONLY: dg_rho0_type
      22             :    USE dg_types,                        ONLY: dg_get,&
      23             :                                               dg_type
      24             :    USE dgs,                             ONLY: dg_get_patch,&
      25             :                                               dg_get_strucfac,&
      26             :                                               dg_sum_patch,&
      27             :                                               dg_sum_patch_force_1d,&
      28             :                                               dg_sum_patch_force_3d
      29             :    USE ewald_environment_types,         ONLY: ewald_env_get,&
      30             :                                               ewald_environment_type
      31             :    USE ewald_pw_types,                  ONLY: ewald_pw_get,&
      32             :                                               ewald_pw_type
      33             :    USE kinds,                           ONLY: dp
      34             :    USE mathconstants,                   ONLY: fourpi
      35             :    USE message_passing,                 ONLY: mp_comm_type
      36             :    USE particle_types,                  ONLY: particle_type
      37             :    USE pme_tools,                       ONLY: get_center,&
      38             :                                               set_list
      39             :    USE pw_grid_types,                   ONLY: pw_grid_type
      40             :    USE pw_methods,                      ONLY: pw_copy,&
      41             :                                               pw_derive,&
      42             :                                               pw_integral_a2b,&
      43             :                                               pw_transfer
      44             :    USE pw_poisson_methods,              ONLY: pw_poisson_solve
      45             :    USE pw_poisson_types,                ONLY: pw_poisson_type
      46             :    USE pw_pool_types,                   ONLY: pw_pool_type
      47             :    USE pw_types,                        ONLY: pw_c1d_gs_type,&
      48             :                                               pw_r3d_rs_type
      49             :    USE realspace_grid_types,            ONLY: realspace_grid_desc_type,&
      50             :                                               realspace_grid_type,&
      51             :                                               rs_grid_create,&
      52             :                                               rs_grid_release,&
      53             :                                               rs_grid_set_box,&
      54             :                                               rs_grid_zero,&
      55             :                                               transfer_pw2rs,&
      56             :                                               transfer_rs2pw
      57             :    USE shell_potential_types,           ONLY: shell_kind_type
      58             :    USE structure_factor_types,          ONLY: structure_factor_type
      59             :    USE structure_factors,               ONLY: structure_factor_allocate,&
      60             :                                               structure_factor_deallocate,&
      61             :                                               structure_factor_init
      62             : #include "./base/base_uses.f90"
      63             : 
      64             :    IMPLICIT NONE
      65             : 
      66             :    PRIVATE
      67             :    PUBLIC :: pme_evaluate
      68             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'pme'
      69             : 
      70             : CONTAINS
      71             : 
      72             : ! **************************************************************************************************
      73             : !> \brief ...
      74             : !> \param ewald_env ...
      75             : !> \param ewald_pw ...
      76             : !> \param box ...
      77             : !> \param particle_set ...
      78             : !> \param vg_coulomb ...
      79             : !> \param fg_coulomb ...
      80             : !> \param pv_g ...
      81             : !> \param shell_particle_set ...
      82             : !> \param core_particle_set ...
      83             : !> \param fgshell_coulomb ...
      84             : !> \param fgcore_coulomb ...
      85             : !> \param use_virial ...
      86             : !> \param charges ...
      87             : !> \param atprop ...
      88             : !> \par History
      89             : !>      JGH (15-Mar-2001) : New electrostatic calculation and pressure tensor
      90             : !>      JGH (21-Mar-2001) : Complete rewrite
      91             : !>      JGH (21-Mar-2001) : Introduced real space density type for future
      92             : !>                          parallelisation
      93             : !> \author CJM and APSI
      94             : ! **************************************************************************************************
      95        1818 :    SUBROUTINE pme_evaluate(ewald_env, ewald_pw, box, particle_set, vg_coulomb, &
      96        1818 :                            fg_coulomb, pv_g, shell_particle_set, core_particle_set, &
      97        1818 :                            fgshell_coulomb, fgcore_coulomb, use_virial, charges, atprop)
      98             :       TYPE(ewald_environment_type), POINTER              :: ewald_env
      99             :       TYPE(ewald_pw_type), POINTER                       :: ewald_pw
     100             :       TYPE(cell_type), POINTER                           :: box
     101             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     102             :       REAL(KIND=dp), INTENT(OUT)                         :: vg_coulomb
     103             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
     104             :          OPTIONAL                                        :: fg_coulomb, pv_g
     105             :       TYPE(particle_type), DIMENSION(:), OPTIONAL, &
     106             :          POINTER                                         :: shell_particle_set, core_particle_set
     107             :       REAL(KIND=dp), DIMENSION(:, :), INTENT(OUT), &
     108             :          OPTIONAL                                        :: fgshell_coulomb, fgcore_coulomb
     109             :       LOGICAL, INTENT(IN)                                :: use_virial
     110             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: charges
     111             :       TYPE(atprop_type), POINTER                         :: atprop
     112             : 
     113             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'pme_evaluate'
     114             : 
     115             :       INTEGER                                            :: handle, i, ig, ipart, j, nd(3), npart, &
     116             :                                                             nshell, p1, p2
     117             :       LOGICAL                                            :: is1_core, is2_core
     118             :       REAL(KIND=dp)                                      :: alpha, dvols, fat1, ffa, ffb
     119             :       REAL(KIND=dp), DIMENSION(3)                        :: fat
     120             :       REAL(KIND=dp), DIMENSION(3, 3)                     :: f_stress, h_stress
     121             :       TYPE(dg_rho0_type), POINTER                        :: dg_rho0
     122             :       TYPE(dg_type), POINTER                             :: dg
     123             :       TYPE(mp_comm_type)                                 :: group
     124             :       TYPE(pw_c1d_gs_type)                               :: phi_g, rhob_g
     125        7272 :       TYPE(pw_c1d_gs_type), DIMENSION(3)                 :: dphi_g
     126             :       TYPE(pw_grid_type), POINTER                        :: grid_b, grid_s
     127             :       TYPE(pw_poisson_type), POINTER                     :: poisson_env
     128             :       TYPE(pw_pool_type), POINTER                        :: pw_big_pool, pw_small_pool
     129             :       TYPE(pw_r3d_rs_type)                               :: phi_r, rhob_r, rhos1, rhos2
     130             :       TYPE(realspace_grid_desc_type), POINTER            :: rs_desc
     131       65448 :       TYPE(realspace_grid_type), DIMENSION(3)            :: drpot
     132             :       TYPE(realspace_grid_type), POINTER                 :: rden, rpot
     133             :       TYPE(structure_factor_type)                        :: exp_igr
     134             : 
     135        1818 :       CALL timeset(routineN, handle)
     136        1818 :       NULLIFY (poisson_env, rden)
     137        1818 :       CALL cite_reference(Darden1993)
     138        1818 :       CALL ewald_env_get(ewald_env, alpha=alpha, group=group)
     139             :       CALL ewald_pw_get(ewald_pw, pw_big_pool=pw_big_pool, &
     140             :                         pw_small_pool=pw_small_pool, rs_desc=rs_desc, &
     141        1818 :                         poisson_env=poisson_env, dg=dg)
     142             : 
     143        1818 :       grid_b => pw_big_pool%pw_grid
     144        1818 :       grid_s => pw_small_pool%pw_grid
     145             : 
     146        1818 :       CALL dg_get(dg, dg_rho0=dg_rho0)
     147             : 
     148        1818 :       npart = SIZE(particle_set)
     149             : 
     150        1818 :       CALL structure_factor_init(exp_igr)
     151             : 
     152        1818 :       IF (PRESENT(shell_particle_set)) THEN
     153           0 :          CPASSERT(ASSOCIATED(shell_particle_set))
     154           0 :          CPASSERT(ASSOCIATED(core_particle_set))
     155           0 :          nshell = SIZE(shell_particle_set)
     156             :          CALL structure_factor_allocate(grid_s%bounds, npart, exp_igr, &
     157             :                                         allocate_centre=.TRUE., allocate_shell_e=.TRUE., &
     158           0 :                                         allocate_shell_centre=.TRUE., nshell=nshell)
     159             : 
     160             :       ELSE
     161             :          CALL structure_factor_allocate(grid_s%bounds, npart, exp_igr, &
     162        1818 :                                         allocate_centre=.TRUE.)
     163             :       END IF
     164             : 
     165        1818 :       CALL pw_small_pool%create_pw(rhos1)
     166        1818 :       CALL pw_small_pool%create_pw(rhos2)
     167             : 
     168       38178 :       ALLOCATE (rden)
     169        1818 :       CALL rs_grid_create(rden, rs_desc)
     170        1818 :       CALL rs_grid_set_box(grid_b, rs=rden)
     171        1818 :       CALL rs_grid_zero(rden)
     172             : 
     173        1818 :       CPASSERT(ASSOCIATED(box))
     174             : 
     175        1818 :       IF (rden%desc%parallel .AND. rden%desc%distributed) THEN
     176           0 :          CALL get_center(particle_set, box, exp_igr%centre, exp_igr%delta, grid_b%npts, 1)
     177             :       END IF
     178        1818 :       IF (PRESENT(shell_particle_set) .AND. rden%desc%parallel .AND. rden%desc%distributed) THEN
     179           0 :          CALL get_center(shell_particle_set, box, exp_igr%shell_centre, exp_igr%shell_delta, grid_b%npts, 1)
     180           0 :          CALL get_center(core_particle_set, box, exp_igr%core_centre, exp_igr%core_delta, grid_b%npts, 1)
     181             :       END IF
     182             : 
     183             :       !-------------- DENSITY CALCULATION ----------------
     184             : 
     185        1818 :       ipart = 0
     186       29400 :       DO
     187             : 
     188       14700 :          CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
     189       14700 :          CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
     190       14700 :          IF (p1 == 0 .AND. p2 == 0) EXIT
     191             : 
     192       12882 :          is1_core = (particle_set(p1)%shell_index /= 0)
     193       12882 :          IF (p2 /= 0) THEN
     194       11757 :             is2_core = (particle_set(p2)%shell_index /= 0)
     195             :          ELSE
     196        1125 :             is2_core = .FALSE.
     197             :          END IF
     198             : 
     199             :          ! calculate function on small boxes (we use double packing in FFT)
     200       14700 :          IF (is1_core .OR. is2_core) THEN
     201             :             CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
     202             :                            rhos1, rhos2, is1_core=is1_core, is2_core=is2_core, &
     203           0 :                            core_particle_set=core_particle_set, charges=charges)
     204             : 
     205             :             ! add boxes to real space grid (big box)
     206           0 :             IF (is1_core) THEN
     207           0 :                CALL dg_sum_patch(rden, rhos1, exp_igr%core_centre(:, particle_set(p1)%shell_index))
     208             :             ELSE
     209           0 :                CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1))
     210             :             END IF
     211           0 :             IF (p2 /= 0 .AND. is2_core) THEN
     212           0 :                CALL dg_sum_patch(rden, rhos2, exp_igr%core_centre(:, particle_set(p2)%shell_index))
     213           0 :             ELSE IF (p2 /= 0) THEN
     214           0 :                CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2))
     215             :             END IF
     216             :          ELSE
     217             :             CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
     218       12882 :                            rhos1, rhos2, charges=charges)
     219             :             ! add boxes to real space grid (big box)
     220       12882 :             CALL dg_sum_patch(rden, rhos1, exp_igr%centre(:, p1))
     221       12882 :             IF (p2 /= 0) CALL dg_sum_patch(rden, rhos2, exp_igr%centre(:, p2))
     222             :          END IF
     223             : 
     224             :       END DO
     225        1818 :       IF (PRESENT(shell_particle_set)) THEN
     226           0 :          ipart = 0
     227           0 :          DO
     228           0 :             CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rpot, ipart)
     229           0 :             CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rpot, ipart)
     230           0 :             IF (p1 == 0 .AND. p2 == 0) EXIT
     231             :             ! calculate function on small boxes (we use double packing in FFT)
     232             :             CALL get_patch(dg, shell_particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
     233           0 :                            rhos1, rhos2, is1_shell=.TRUE., is2_shell=.TRUE., charges=charges)
     234             :             ! add boxes to real space grid (big box)
     235           0 :             CALL dg_sum_patch(rpot, rhos1, exp_igr%shell_centre(:, p1))
     236           0 :             IF (p2 /= 0) CALL dg_sum_patch(rpot, rhos2, exp_igr%shell_centre(:, p2))
     237             :          END DO
     238             :       END IF
     239             : 
     240        1818 :       CALL pw_big_pool%create_pw(rhob_r)
     241        1818 :       CALL transfer_rs2pw(rden, rhob_r)
     242             : 
     243             :       !-------------- ELECTROSTATIC CALCULATION -----------
     244             : 
     245             :       ! allocate intermediate arrays
     246        7272 :       DO i = 1, 3
     247        7272 :          CALL pw_big_pool%create_pw(dphi_g(i))
     248             :       END DO
     249        1818 :       CALL pw_big_pool%create_pw(phi_r)
     250             : 
     251        1818 :       CALL pw_poisson_solve(poisson_env, rhob_r, vg_coulomb, phi_r, dphi_g, h_stress)
     252             : 
     253             :       ! atomic energies
     254        1818 :       IF (atprop%energy .OR. atprop%stress) THEN
     255         202 :          dvols = rhos1%pw_grid%dvol
     256        4242 :          ALLOCATE (rpot)
     257         202 :          CALL rs_grid_create(rpot, rs_desc)
     258         202 :          CALL transfer_pw2rs(rpot, phi_r)
     259         202 :          ipart = 0
     260         808 :          DO
     261         404 :             CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
     262         404 :             CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
     263         404 :             IF (p1 == 0 .AND. p2 == 0) EXIT
     264             :             ! integrate box and potential
     265             :             CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
     266         202 :                            rhos1, rhos2, charges=charges)
     267             :             ! add boxes to real space grid (big box)
     268         202 :             CALL dg_sum_patch_force_1d(rpot, rhos1, exp_igr%centre(:, p1), fat1)
     269         202 :             IF (atprop%energy) THEN
     270         202 :                atprop%atener(p1) = atprop%atener(p1) + 0.5_dp*fat1*dvols
     271             :             END IF
     272         202 :             IF (atprop%stress) THEN
     273         202 :                atprop%atstress(1, 1, p1) = atprop%atstress(1, 1, p1) + 0.5_dp*fat1*dvols
     274         202 :                atprop%atstress(2, 2, p1) = atprop%atstress(2, 2, p1) + 0.5_dp*fat1*dvols
     275         202 :                atprop%atstress(3, 3, p1) = atprop%atstress(3, 3, p1) + 0.5_dp*fat1*dvols
     276             :             END IF
     277         404 :             IF (p2 /= 0) THEN
     278         101 :                CALL dg_sum_patch_force_1d(rpot, rhos2, exp_igr%centre(:, p2), fat1)
     279         101 :                IF (atprop%energy) THEN
     280         101 :                   atprop%atener(p2) = atprop%atener(p2) + 0.5_dp*fat1*dvols
     281             :                END IF
     282         101 :                IF (atprop%stress) THEN
     283         101 :                   atprop%atstress(1, 1, p2) = atprop%atstress(1, 1, p2) + 0.5_dp*fat1*dvols
     284         101 :                   atprop%atstress(2, 2, p2) = atprop%atstress(2, 2, p2) + 0.5_dp*fat1*dvols
     285         101 :                   atprop%atstress(3, 3, p2) = atprop%atstress(3, 3, p2) + 0.5_dp*fat1*dvols
     286             :                END IF
     287             :             END IF
     288             :          END DO
     289         202 :          IF (atprop%stress) THEN
     290         202 :             CALL pw_big_pool%create_pw(phi_g)
     291         202 :             CALL pw_big_pool%create_pw(rhob_g)
     292         202 :             ffa = (0.5_dp/dg_rho0%zet(1))**2
     293         202 :             ffb = 1.0_dp/fourpi
     294         808 :             DO i = 1, 3
     295     4325022 :                DO ig = grid_b%first_gne0, grid_b%ngpts_cut_local
     296     4324416 :                   phi_g%array(ig) = ffb*dphi_g(i)%array(ig)*(ffa*grid_b%gsq(ig) + 1.0_dp)
     297     4325022 :                   phi_g%array(ig) = phi_g%array(ig)*poisson_env%green_fft%influence_fn%array(ig)
     298             :                END DO
     299         606 :                IF (grid_b%have_g0) phi_g%array(1) = 0.0_dp
     300        2020 :                DO j = 1, i
     301        1212 :                   CALL pw_copy(phi_g, rhob_g)
     302        1212 :                   nd = 0
     303        1212 :                   nd(j) = 1
     304        1212 :                   CALL pw_derive(rhob_g, nd)
     305        1212 :                   CALL pw_transfer(rhob_g, rhob_r)
     306        1212 :                   CALL transfer_pw2rs(rpot, rhob_r)
     307             : 
     308        1212 :                   ipart = 0
     309         606 :                   DO
     310        2424 :                      CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
     311        2424 :                      CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
     312        2424 :                      IF (p1 == 0 .AND. p2 == 0) EXIT
     313             :                      ! integrate box and potential
     314             :                      CALL get_patch(dg, particle_set, exp_igr, box, p1, p2, grid_b, grid_s, &
     315        1212 :                                     rhos1, rhos2, charges=charges)
     316             :                      ! add boxes to real space grid (big box)
     317        1212 :                      CALL dg_sum_patch_force_1d(rpot, rhos1, exp_igr%centre(:, p1), fat1)
     318        1212 :                      atprop%atstress(i, j, p1) = atprop%atstress(i, j, p1) + fat1*dvols
     319        1212 :                      IF (i /= j) atprop%atstress(j, i, p1) = atprop%atstress(j, i, p1) + fat1*dvols
     320        2424 :                      IF (p2 /= 0) THEN
     321         606 :                         CALL dg_sum_patch_force_1d(rpot, rhos2, exp_igr%centre(:, p2), fat1)
     322         606 :                         atprop%atstress(i, j, p2) = atprop%atstress(i, j, p2) + fat1*dvols
     323         606 :                         IF (i /= j) atprop%atstress(j, i, p2) = atprop%atstress(j, i, p2) + fat1*dvols
     324             :                      END IF
     325             :                   END DO
     326             :                END DO
     327             :             END DO
     328         202 :             CALL pw_big_pool%give_back_pw(phi_g)
     329         202 :             CALL pw_big_pool%give_back_pw(rhob_g)
     330             :          END IF
     331         202 :          CALL rs_grid_release(rpot)
     332         202 :          DEALLOCATE (rpot)
     333             :       END IF
     334             : 
     335        1818 :       CALL pw_big_pool%give_back_pw(phi_r)
     336             : 
     337             :       !---------- END OF ELECTROSTATIC CALCULATION --------
     338             : 
     339             :       !------------- STRESS TENSOR CALCULATION ------------
     340             : 
     341        1818 :       IF ((use_virial) .AND. (PRESENT(pv_g))) THEN
     342         808 :          DO i = 1, 3
     343        2020 :             DO j = i, 3
     344        1212 :                f_stress(i, j) = pw_integral_a2b(dphi_g(i), dphi_g(j))
     345        1818 :                f_stress(j, i) = f_stress(i, j)
     346             :             END DO
     347             :          END DO
     348         202 :          ffa = (1.0_dp/fourpi)*(0.5_dp/dg_rho0%zet(1))**2
     349        2626 :          f_stress = -ffa*f_stress
     350        2626 :          pv_g = h_stress + f_stress
     351             :       END IF
     352             : 
     353             :       !--------END OF STRESS TENSOR CALCULATION -----------
     354             : 
     355        7272 :       DO i = 1, 3
     356        5454 :          CALL rs_grid_create(drpot(i), rs_desc)
     357        5454 :          CALL rs_grid_set_box(grid_b, rs=drpot(i))
     358        5454 :          CALL pw_transfer(dphi_g(i), rhob_r)
     359        5454 :          CALL pw_big_pool%give_back_pw(dphi_g(i))
     360        7272 :          CALL transfer_pw2rs(drpot(i), rhob_r)
     361             :       END DO
     362             : 
     363        1818 :       CALL pw_big_pool%give_back_pw(rhob_r)
     364             : 
     365             :       !----------------- FORCE CALCULATION ----------------
     366             : 
     367             :       ! initialize the forces
     368        1818 :       IF (PRESENT(fg_coulomb)) THEN
     369      199026 :          fg_coulomb = 0.0_dp
     370        1818 :          dvols = rhos1%pw_grid%dvol
     371             : 
     372        1818 :          ipart = 0
     373       29400 :          DO
     374             : 
     375       14700 :             CALL set_list(particle_set, npart, exp_igr%centre, p1, rden, ipart, exp_igr%core_centre)
     376       14700 :             CALL set_list(particle_set, npart, exp_igr%centre, p2, rden, ipart, exp_igr%core_centre)
     377       14700 :             IF (p1 == 0 .AND. p2 == 0) EXIT
     378             : 
     379       12882 :             is1_core = (particle_set(p1)%shell_index /= 0)
     380       12882 :             IF (p2 /= 0) THEN
     381       11757 :                is2_core = (particle_set(p2)%shell_index /= 0)
     382             :             ELSE
     383        1125 :                is2_core = .FALSE.
     384             :             END IF
     385             : 
     386             :             ! calculate function on small boxes (we use double packing in FFT)
     387             : 
     388             :             CALL get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, &
     389       12882 :                                  is1_core=is1_core, is2_core=is2_core, charges=charges)
     390             : 
     391             :             ! sum boxes on real space grids (big box)
     392       12882 :             IF (is1_core) THEN
     393             :                CALL dg_sum_patch_force_3d(drpot, rhos1, &
     394           0 :                                           exp_igr%core_centre(:, particle_set(p1)%shell_index), fat)
     395             :                fgcore_coulomb(1, particle_set(p1)%shell_index) = &
     396           0 :                   fgcore_coulomb(1, particle_set(p1)%shell_index) - fat(1)*dvols
     397             :                fgcore_coulomb(2, particle_set(p1)%shell_index) = &
     398           0 :                   fgcore_coulomb(2, particle_set(p1)%shell_index) - fat(2)*dvols
     399             :                fgcore_coulomb(3, particle_set(p1)%shell_index) = &
     400           0 :                   fgcore_coulomb(3, particle_set(p1)%shell_index) - fat(3)*dvols
     401             :             ELSE
     402       12882 :                CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%centre(:, p1), fat)
     403       12882 :                fg_coulomb(1, p1) = fg_coulomb(1, p1) - fat(1)*dvols
     404       12882 :                fg_coulomb(2, p1) = fg_coulomb(2, p1) - fat(2)*dvols
     405       12882 :                fg_coulomb(3, p1) = fg_coulomb(3, p1) - fat(3)*dvols
     406             :             END IF
     407       14700 :             IF (p2 /= 0 .AND. is2_core) THEN
     408             :                CALL dg_sum_patch_force_3d(drpot, rhos1, &
     409           0 :                                           exp_igr%core_centre(:, particle_set(p2)%shell_index), fat)
     410             :                fgcore_coulomb(1, particle_set(p2)%shell_index) = &
     411           0 :                   fgcore_coulomb(1, particle_set(p2)%shell_index) - fat(1)*dvols
     412             :                fgcore_coulomb(2, particle_set(p2)%shell_index) = &
     413           0 :                   fgcore_coulomb(2, particle_set(p2)%shell_index) - fat(2)*dvols
     414             :                fgcore_coulomb(3, particle_set(p2)%shell_index) = &
     415           0 :                   fgcore_coulomb(3, particle_set(p2)%shell_index) - fat(3)*dvols
     416       12882 :             ELSEIF (p2 /= 0) THEN
     417       11757 :                CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%centre(:, p2), fat)
     418       11757 :                fg_coulomb(1, p2) = fg_coulomb(1, p2) - fat(1)*dvols
     419       11757 :                fg_coulomb(2, p2) = fg_coulomb(2, p2) - fat(2)*dvols
     420       11757 :                fg_coulomb(3, p2) = fg_coulomb(3, p2) - fat(3)*dvols
     421             :             END IF
     422             : 
     423             :          END DO
     424             :       END IF
     425        1818 :       IF (PRESENT(fgshell_coulomb)) THEN
     426           0 :          fgshell_coulomb = 0.0_dp
     427           0 :          dvols = rhos1%pw_grid%dvol
     428             : 
     429           0 :          ipart = 0
     430           0 :          DO
     431           0 :             CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p1, rden, ipart)
     432           0 :             CALL set_list(shell_particle_set, nshell, exp_igr%shell_centre, p2, rden, ipart)
     433           0 :             IF (p1 == 0 .AND. p2 == 0) EXIT
     434             : 
     435             :             ! calculate function on small boxes (we use double packing in FFT)
     436             :             CALL get_patch_again(dg, shell_particle_set, exp_igr, p1, p2, rhos1, rhos2, &
     437           0 :                                  is1_shell=.TRUE., is2_shell=.TRUE., charges=charges)
     438             : 
     439             :             ! sum boxes on real space grids (big box)
     440           0 :             CALL dg_sum_patch_force_3d(drpot, rhos1, exp_igr%shell_centre(:, p1), fat)
     441           0 :             fgshell_coulomb(1, p1) = fgshell_coulomb(1, p1) - fat(1)*dvols
     442           0 :             fgshell_coulomb(2, p1) = fgshell_coulomb(2, p1) - fat(2)*dvols
     443           0 :             fgshell_coulomb(3, p1) = fgshell_coulomb(3, p1) - fat(3)*dvols
     444           0 :             IF (p2 /= 0) THEN
     445           0 :                CALL dg_sum_patch_force_3d(drpot, rhos2, exp_igr%shell_centre(:, p2), fat)
     446           0 :                fgshell_coulomb(1, p2) = fgshell_coulomb(1, p2) - fat(1)*dvols
     447           0 :                fgshell_coulomb(2, p2) = fgshell_coulomb(2, p2) - fat(2)*dvols
     448           0 :                fgshell_coulomb(3, p2) = fgshell_coulomb(3, p2) - fat(3)*dvols
     449             :             END IF
     450             :          END DO
     451             : 
     452             :       END IF
     453             :       !--------------END OF FORCE CALCULATION -------------
     454             : 
     455             :       !------------------CLEANING UP ----------------------
     456             : 
     457        1818 :       CALL rs_grid_release(rden)
     458        1818 :       DEALLOCATE (rden)
     459        7272 :       DO i = 1, 3
     460        7272 :          CALL rs_grid_release(drpot(i))
     461             :       END DO
     462             : 
     463        1818 :       CALL pw_small_pool%give_back_pw(rhos1)
     464        1818 :       CALL pw_small_pool%give_back_pw(rhos2)
     465        1818 :       CALL structure_factor_deallocate(exp_igr)
     466             : 
     467        1818 :       CALL timestop(handle)
     468             : 
     469       29088 :    END SUBROUTINE pme_evaluate
     470             : 
     471             : ! **************************************************************************************************
     472             : !> \brief Calculates local density in a small box
     473             : !> \param dg ...
     474             : !> \param particle_set ...
     475             : !> \param exp_igr ...
     476             : !> \param box ...
     477             : !> \param p1 ...
     478             : !> \param p2 ...
     479             : !> \param grid_b ...
     480             : !> \param grid_s ...
     481             : !> \param rhos1 ...
     482             : !> \param rhos2 ...
     483             : !> \param is1_core ...
     484             : !> \param is2_core ...
     485             : !> \param is1_shell ...
     486             : !> \param is2_shell ...
     487             : !> \param core_particle_set ...
     488             : !> \param charges ...
     489             : !> \par History
     490             : !>      JGH (23-Mar-2001) : Switch to integer from particle list pointers
     491             : !> \author JGH (21-Mar-2001)
     492             : ! **************************************************************************************************
     493       14296 :    SUBROUTINE get_patch(dg, particle_set, exp_igr, box, p1, p2, &
     494             :                         grid_b, grid_s, rhos1, rhos2, is1_core, is2_core, is1_shell, &
     495             :                         is2_shell, core_particle_set, charges)
     496             : 
     497             :       TYPE(dg_type), POINTER                             :: dg
     498             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     499             :       TYPE(structure_factor_type)                        :: exp_igr
     500             :       TYPE(cell_type), POINTER                           :: box
     501             :       INTEGER, INTENT(IN)                                :: p1, p2
     502             :       TYPE(pw_grid_type), INTENT(IN)                     :: grid_b, grid_s
     503             :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: rhos1, rhos2
     504             :       LOGICAL, OPTIONAL                                  :: is1_core, is2_core, is1_shell, is2_shell
     505             :       TYPE(particle_type), DIMENSION(:), OPTIONAL, &
     506             :          POINTER                                         :: core_particle_set
     507             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: charges
     508             : 
     509       14296 :       COMPLEX(KIND=dp), DIMENSION(:), POINTER            :: ex1, ex2, ey1, ey2, ez1, ez2
     510       14296 :       INTEGER, DIMENSION(:), POINTER                     :: center1, center2
     511             :       LOGICAL                                            :: my_is1_core, my_is1_shell, my_is2_core, &
     512             :                                                             my_is2_shell, use_charge_array
     513             :       REAL(KIND=dp)                                      :: q1, q2
     514             :       REAL(KIND=dp), DIMENSION(3)                        :: r1, r2
     515             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     516             :       TYPE(dg_rho0_type), POINTER                        :: dg_rho0
     517             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho0
     518             :       TYPE(shell_kind_type), POINTER                     :: shell
     519             : 
     520       14296 :       NULLIFY (shell)
     521       14296 :       use_charge_array = PRESENT(charges)
     522       14296 :       IF (use_charge_array) use_charge_array = ASSOCIATED(charges)
     523       14296 :       my_is1_core = .FALSE.
     524       14296 :       my_is2_core = .FALSE.
     525       14296 :       IF (PRESENT(is1_core)) my_is1_core = is1_core
     526       14296 :       IF (PRESENT(is2_core)) my_is2_core = is2_core
     527       14296 :       IF (my_is1_core .OR. my_is2_core) THEN
     528           0 :          CPASSERT(PRESENT(core_particle_set))
     529             :       END IF
     530       14296 :       my_is1_shell = .FALSE.
     531       14296 :       my_is2_shell = .FALSE.
     532       14296 :       IF (PRESENT(is1_shell)) my_is1_shell = is1_shell
     533       14296 :       IF (PRESENT(is2_shell)) my_is2_shell = is2_shell
     534       14296 :       IF (my_is1_core .AND. my_is1_shell) THEN
     535           0 :          CPABORT("Shell-model: cannot be core and shell simultaneously")
     536             :       END IF
     537             : 
     538       14296 :       CALL dg_get(dg, dg_rho0=dg_rho0)
     539       14296 :       rho0 => dg_rho0%density
     540             : 
     541       14296 :       IF (my_is1_core) THEN
     542           0 :          r1 = core_particle_set(particle_set(p1)%shell_index)%r
     543             :       ELSE
     544       57184 :          r1 = particle_set(p1)%r
     545             :       END IF
     546       14296 :       atomic_kind => particle_set(p1)%atomic_kind
     547       14296 :       IF (my_is1_core) THEN
     548           0 :          CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     549           0 :          q1 = shell%charge_core
     550       14296 :       ELSE IF (my_is1_shell) THEN
     551           0 :          CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     552           0 :          q1 = shell%charge_shell
     553             :       ELSE
     554       14296 :          CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1)
     555             :       END IF
     556       14296 :       IF (use_charge_array) q1 = charges(p1)
     557             : 
     558       14296 :       IF (my_is1_shell) THEN
     559           0 :          center1 => exp_igr%shell_centre(:, p1)
     560           0 :          ex1 => exp_igr%shell_ex(:, p1)
     561           0 :          ey1 => exp_igr%shell_ey(:, p1)
     562           0 :          ez1 => exp_igr%shell_ez(:, p1)
     563       14296 :       ELSEIF (my_is1_core) THEN
     564           0 :          center1 => exp_igr%core_centre(:, particle_set(p1)%shell_index)
     565           0 :          ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index)
     566           0 :          ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index)
     567           0 :          ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index)
     568             :       ELSE
     569       14296 :          center1 => exp_igr%centre(:, p1)
     570       14296 :          ex1 => exp_igr%ex(:, p1)
     571       14296 :          ey1 => exp_igr%ey(:, p1)
     572       14296 :          ez1 => exp_igr%ez(:, p1)
     573             :       END IF
     574             : 
     575       14296 :       CPASSERT(ASSOCIATED(box))
     576             : 
     577             :       CALL dg_get_strucfac(box%hmat, r1, grid_s%npts, grid_b%npts, center1, &
     578       14296 :                            exp_igr%lb, ex1, ey1, ez1)
     579             : 
     580       14296 :       IF (p2 /= 0) THEN
     581       12464 :          IF (my_is2_core) THEN
     582           0 :             r2 = core_particle_set(particle_set(p2)%shell_index)%r
     583             :          ELSE
     584       49856 :             r2 = particle_set(p2)%r
     585             :          END IF
     586       12464 :          atomic_kind => particle_set(p2)%atomic_kind
     587       12464 :          IF (my_is2_core) THEN
     588           0 :             CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     589           0 :             q2 = shell%charge_core
     590       12464 :          ELSE IF (my_is2_shell) THEN
     591           0 :             CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     592           0 :             q2 = shell%charge_shell
     593             :          ELSE
     594       12464 :             CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2)
     595             :          END IF
     596       12464 :          IF (use_charge_array) q2 = charges(p2)
     597             : 
     598       12464 :          IF (my_is2_shell) THEN
     599           0 :             center2 => exp_igr%shell_centre(:, p2)
     600           0 :             ex2 => exp_igr%shell_ex(:, p2)
     601           0 :             ey2 => exp_igr%shell_ey(:, p2)
     602           0 :             ez2 => exp_igr%shell_ez(:, p2)
     603       12464 :          ELSEIF (my_is2_core) THEN
     604           0 :             center2 => exp_igr%core_centre(:, particle_set(p2)%shell_index)
     605           0 :             ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index)
     606           0 :             ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index)
     607           0 :             ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index)
     608             :          ELSE
     609       12464 :             center2 => exp_igr%centre(:, p2)
     610       12464 :             ex2 => exp_igr%ex(:, p2)
     611       12464 :             ey2 => exp_igr%ey(:, p2)
     612       12464 :             ez2 => exp_igr%ez(:, p2)
     613             :          END IF
     614             :          CALL dg_get_strucfac(box%hmat, r2, grid_s%npts, grid_b%npts, center2, &
     615       12464 :                               exp_igr%lb, ex2, ey2, ez2)
     616             :       END IF
     617             : 
     618       14296 :       IF (p2 == 0) THEN
     619        1832 :          CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1)
     620             :       ELSE
     621       12464 :          CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, ex1, ey1, ez1, ex2, ey2, ez2)
     622             :       END IF
     623             : 
     624       14296 :    END SUBROUTINE get_patch
     625             : 
     626             : ! **************************************************************************************************
     627             : !> \brief ...
     628             : !> \param dg ...
     629             : !> \param particle_set ...
     630             : !> \param exp_igr ...
     631             : !> \param p1 ...
     632             : !> \param p2 ...
     633             : !> \param rhos1 ...
     634             : !> \param rhos2 ...
     635             : !> \param is1_core ...
     636             : !> \param is2_core ...
     637             : !> \param is1_shell ...
     638             : !> \param is2_shell ...
     639             : !> \param charges ...
     640             : ! **************************************************************************************************
     641       12882 :    SUBROUTINE get_patch_again(dg, particle_set, exp_igr, p1, p2, rhos1, rhos2, is1_core, &
     642             :                               is2_core, is1_shell, is2_shell, charges)
     643             : 
     644             :       TYPE(dg_type), POINTER                             :: dg
     645             :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     646             :       TYPE(structure_factor_type)                        :: exp_igr
     647             :       INTEGER, INTENT(IN)                                :: p1, p2
     648             :       TYPE(pw_r3d_rs_type), INTENT(INOUT)                :: rhos1, rhos2
     649             :       LOGICAL, OPTIONAL                                  :: is1_core, is2_core, is1_shell, is2_shell
     650             :       REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER     :: charges
     651             : 
     652       12882 :       COMPLEX(KIND=dp), DIMENSION(:), POINTER            :: ex1, ex2, ey1, ey2, ez1, ez2
     653             :       LOGICAL                                            :: my_is1_core, my_is1_shell, my_is2_core, &
     654             :                                                             my_is2_shell, use_charge_array
     655             :       REAL(KIND=dp)                                      :: q1, q2
     656             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     657             :       TYPE(dg_rho0_type), POINTER                        :: dg_rho0
     658             :       TYPE(pw_r3d_rs_type), POINTER                      :: rho0
     659             :       TYPE(shell_kind_type), POINTER                     :: shell
     660             : 
     661       12882 :       NULLIFY (shell)
     662       12882 :       use_charge_array = PRESENT(charges)
     663       12882 :       IF (use_charge_array) use_charge_array = ASSOCIATED(charges)
     664       12882 :       my_is1_core = .FALSE.
     665       12882 :       my_is2_core = .FALSE.
     666       12882 :       IF (PRESENT(is1_core)) my_is1_core = is1_core
     667       12882 :       IF (PRESENT(is2_core)) my_is2_core = is2_core
     668       12882 :       my_is1_shell = .FALSE.
     669       12882 :       my_is2_shell = .FALSE.
     670       12882 :       IF (PRESENT(is1_shell)) my_is1_shell = is1_shell
     671       12882 :       IF (PRESENT(is2_shell)) my_is2_shell = is2_shell
     672             : 
     673       12882 :       CALL dg_get(dg, dg_rho0=dg_rho0)
     674       12882 :       rho0 => dg_rho0%density
     675             : 
     676       12882 :       atomic_kind => particle_set(p1)%atomic_kind
     677       12882 :       IF (my_is1_core) THEN
     678           0 :          CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     679           0 :          q1 = shell%charge_core
     680       12882 :       ELSE IF (my_is1_shell) THEN
     681           0 :          CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     682           0 :          q1 = shell%charge_shell
     683             :       ELSE
     684       12882 :          CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q1)
     685             :       END IF
     686       12882 :       IF (use_charge_array) q1 = charges(p1)
     687       12882 :       IF (my_is1_core) THEN
     688           0 :          ex1 => exp_igr%core_ex(:, particle_set(p1)%shell_index)
     689           0 :          ey1 => exp_igr%core_ey(:, particle_set(p1)%shell_index)
     690           0 :          ez1 => exp_igr%core_ez(:, particle_set(p1)%shell_index)
     691       12882 :       ELSEIF (my_is1_shell) THEN
     692           0 :          ex1 => exp_igr%shell_ex(:, p1)
     693           0 :          ey1 => exp_igr%shell_ey(:, p1)
     694           0 :          ez1 => exp_igr%shell_ez(:, p1)
     695             :       ELSE
     696       12882 :          ex1 => exp_igr%ex(:, p1)
     697       12882 :          ey1 => exp_igr%ey(:, p1)
     698       12882 :          ez1 => exp_igr%ez(:, p1)
     699             :       END IF
     700             : 
     701       12882 :       IF (p2 /= 0) THEN
     702       11757 :          atomic_kind => particle_set(p2)%atomic_kind
     703       11757 :          IF (my_is2_core) THEN
     704           0 :             CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     705           0 :             q2 = shell%charge_core
     706       11757 :          ELSE IF (my_is2_shell) THEN
     707           0 :             CALL get_atomic_kind(atomic_kind=atomic_kind, shell=shell)
     708           0 :             q2 = shell%charge_shell
     709             :          ELSE
     710       11757 :             CALL get_atomic_kind(atomic_kind=atomic_kind, qeff=q2)
     711             :          END IF
     712       11757 :          IF (use_charge_array) q2 = charges(p2)
     713       11757 :          IF (my_is2_core) THEN
     714           0 :             ex2 => exp_igr%core_ex(:, particle_set(p2)%shell_index)
     715           0 :             ey2 => exp_igr%core_ey(:, particle_set(p2)%shell_index)
     716           0 :             ez2 => exp_igr%core_ez(:, particle_set(p2)%shell_index)
     717       11757 :          ELSEIF (my_is2_shell) THEN
     718           0 :             ex2 => exp_igr%shell_ex(:, p2)
     719           0 :             ey2 => exp_igr%shell_ey(:, p2)
     720           0 :             ez2 => exp_igr%shell_ez(:, p2)
     721             :          ELSE
     722       11757 :             ex2 => exp_igr%ex(:, p2)
     723       11757 :             ey2 => exp_igr%ey(:, p2)
     724       11757 :             ez2 => exp_igr%ez(:, p2)
     725             :          END IF
     726             :       END IF
     727             : 
     728       12882 :       IF (p2 == 0) THEN
     729        1125 :          CALL dg_get_patch(rho0, rhos1, q1, ex1, ey1, ez1)
     730             :       ELSE
     731             :          CALL dg_get_patch(rho0, rhos1, rhos2, q1, q2, &
     732       11757 :                            ex1, ey1, ez1, ex2, ey2, ez2)
     733             :       END IF
     734             : 
     735       12882 :    END SUBROUTINE get_patch_again
     736             : 
     737             : END MODULE pme
     738             : 

Generated by: LCOV version 1.15