LCOV - code coverage report
Current view: top level - src - moments_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 79.3 % 87 69
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 1 1

            Line data    Source code
       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              : !> \brief Calculates the moment integrals <a|r^m|b>
      10              : !> \par History
      11              : !>      12.2007 [tlaino] - Splitting common routines to QS and FIST
      12              : !>      06.2009 [tlaino] - Extending to molecular dipoles (interval of atoms)
      13              : !> \author JGH (20.07.2006)
      14              : ! **************************************************************************************************
      15              : MODULE moments_utils
      16              :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      17              :                                               get_atomic_kind
      18              :    USE cell_types,                      ONLY: cell_type,&
      19              :                                               pbc
      20              :    USE distribution_1d_types,           ONLY: distribution_1d_type
      21              :    USE fist_environment_types,          ONLY: fist_env_get,&
      22              :                                               fist_environment_type
      23              :    USE input_constants,                 ONLY: use_mom_ref_coac,&
      24              :                                               use_mom_ref_com,&
      25              :                                               use_mom_ref_user,&
      26              :                                               use_mom_ref_zero
      27              :    USE kinds,                           ONLY: dp
      28              :    USE message_passing,                 ONLY: mp_para_env_type
      29              :    USE particle_types,                  ONLY: particle_type
      30              :    USE qs_environment_types,            ONLY: get_qs_env,&
      31              :                                               qs_environment_type
      32              :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      33              :                                               qs_kind_type
      34              : #include "./base/base_uses.f90"
      35              : 
      36              :    IMPLICIT NONE
      37              : 
      38              :    PRIVATE
      39              : 
      40              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'moments_utils'
      41              : 
      42              : ! *** Public subroutines ***
      43              : 
      44              :    PUBLIC :: get_reference_point
      45              : 
      46              : CONTAINS
      47              : 
      48              : ! **************************************************************************************************
      49              : !> \brief ...
      50              : !> \param rpoint ...
      51              : !> \param drpoint ...
      52              : !> \param qs_env ...
      53              : !> \param fist_env ...
      54              : !> \param reference ...
      55              : !> \param ref_point ...
      56              : !> \param ifirst ...
      57              : !> \param ilast ...
      58              : ! **************************************************************************************************
      59        23819 :    SUBROUTINE get_reference_point(rpoint, drpoint, qs_env, fist_env, reference, ref_point, &
      60              :                                   ifirst, ilast)
      61              :       REAL(dp), DIMENSION(3), INTENT(OUT)                :: rpoint
      62              :       REAL(dp), DIMENSION(3), INTENT(OUT), OPTIONAL      :: drpoint
      63              :       TYPE(qs_environment_type), OPTIONAL, POINTER       :: qs_env
      64              :       TYPE(fist_environment_type), OPTIONAL, POINTER     :: fist_env
      65              :       INTEGER, INTENT(IN)                                :: reference
      66              :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
      67              :       INTEGER, INTENT(IN), OPTIONAL                      :: ifirst, ilast
      68              : 
      69              :       INTEGER                                            :: akind, ia, iatom, ikind
      70              :       LOGICAL                                            :: do_molecule
      71              :       REAL(dp)                                           :: charge, mass, mass_low, mtot, ztot
      72              :       REAL(dp), DIMENSION(3)                             :: center, ria
      73              :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
      74              :       TYPE(cell_type), POINTER                           :: cell
      75              :       TYPE(distribution_1d_type), POINTER                :: local_particles
      76              :       TYPE(mp_para_env_type), POINTER                    :: para_env
      77        23819 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      78        23819 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      79              : 
      80            0 :       CPASSERT(PRESENT(ifirst) .EQV. PRESENT(ilast))
      81        23819 :       NULLIFY (cell, particle_set, qs_kind_set, local_particles, para_env)
      82        23819 :       IF (PRESENT(qs_env)) THEN
      83              :          CALL get_qs_env(qs_env, cell=cell, particle_set=particle_set, &
      84              :                          qs_kind_set=qs_kind_set, &
      85         2615 :                          local_particles=local_particles, para_env=para_env)
      86              :       END IF
      87        23819 :       IF (PRESENT(fist_env)) THEN
      88              :          CALL fist_env_get(fist_env, cell=cell, particle_set=particle_set, &
      89        21204 :                            local_particles=local_particles, para_env=para_env)
      90              :       END IF
      91        23819 :       do_molecule = .FALSE.
      92        23819 :       IF (PRESENT(ifirst) .AND. PRESENT(ilast)) do_molecule = .TRUE.
      93        23819 :       IF (PRESENT(drpoint)) drpoint = 0.0_dp
      94        23819 :       SELECT CASE (reference)
      95              :       CASE DEFAULT
      96            0 :          CPABORT("Type of reference point not implemented")
      97              :       CASE (use_mom_ref_com)
      98         1496 :          rpoint = 0._dp
      99         1496 :          mtot = 0._dp
     100         1496 :          center(:) = 0._dp
     101         1496 :          IF (do_molecule) THEN
     102            6 :             mass_low = -HUGE(mass_low)
     103              :             ! fold the molecule around the heaviest atom in the molecule
     104           24 :             DO iatom = ifirst, ilast
     105           18 :                atomic_kind => particle_set(iatom)%atomic_kind
     106           18 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     107           24 :                IF (mass > mass_low) THEN
     108            6 :                   mass_low = mass
     109           24 :                   center = particle_set(iatom)%r
     110              :                END IF
     111              :             END DO
     112           24 :             DO iatom = ifirst, ilast
     113           72 :                ria = particle_set(iatom)%r
     114          126 :                ria = pbc(ria - center, cell) + center
     115           18 :                atomic_kind => particle_set(iatom)%atomic_kind
     116           18 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     117           72 :                rpoint(:) = rpoint(:) + mass*ria(:)
     118           18 :                IF (PRESENT(drpoint)) drpoint = drpoint + mass*particle_set(iatom)%v
     119           42 :                mtot = mtot + mass
     120              :             END DO
     121              :          ELSE
     122         5128 :             DO ikind = 1, SIZE(local_particles%n_el)
     123         7842 :                DO ia = 1, local_particles%n_el(ikind)
     124         2714 :                   iatom = local_particles%list(ikind)%array(ia)
     125        10856 :                   ria = particle_set(iatom)%r
     126        10856 :                   ria = pbc(ria, cell)
     127         2714 :                   atomic_kind => particle_set(iatom)%atomic_kind
     128         2714 :                   CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     129        10856 :                   rpoint(:) = rpoint(:) + mass*ria(:)
     130         6737 :                   IF (PRESENT(drpoint)) drpoint = drpoint + mass*particle_set(iatom)%v
     131         9066 :                   mtot = mtot + mass
     132              :                END DO
     133              :             END DO
     134         1490 :             CALL para_env%sum(rpoint)
     135         1490 :             CALL para_env%sum(mtot)
     136              :          END IF
     137         1496 :          IF (ABS(mtot) > 0._dp) THEN
     138         5984 :             rpoint(:) = rpoint(:)/mtot
     139         3944 :             IF (PRESENT(drpoint)) drpoint = drpoint/mtot
     140              :          END IF
     141              :       CASE (use_mom_ref_coac)
     142           54 :          rpoint = 0._dp
     143           54 :          ztot = 0._dp
     144           54 :          center(:) = 0._dp
     145           54 :          IF (do_molecule) THEN
     146            0 :             mass_low = -HUGE(mass_low)
     147              :             ! fold the molecule around the heaviest atom in the molecule
     148            0 :             DO iatom = ifirst, ilast
     149            0 :                atomic_kind => particle_set(iatom)%atomic_kind
     150            0 :                CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
     151            0 :                IF (mass > mass_low) THEN
     152            0 :                   mass_low = mass
     153            0 :                   center = particle_set(iatom)%r
     154              :                END IF
     155              :             END DO
     156            0 :             DO iatom = ifirst, ilast
     157            0 :                ria = particle_set(iatom)%r
     158            0 :                ria = pbc(ria - center, cell) + center
     159            0 :                atomic_kind => particle_set(iatom)%atomic_kind
     160            0 :                CALL get_atomic_kind(atomic_kind, kind_number=akind)
     161            0 :                CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
     162            0 :                rpoint(:) = rpoint(:) + charge*ria(:)
     163            0 :                IF (PRESENT(drpoint)) drpoint = drpoint + charge*particle_set(iatom)%v
     164            0 :                ztot = ztot + charge
     165              :             END DO
     166              :          ELSE
     167          148 :             DO ikind = 1, SIZE(local_particles%n_el)
     168          222 :                DO ia = 1, local_particles%n_el(ikind)
     169           74 :                   iatom = local_particles%list(ikind)%array(ia)
     170          296 :                   ria = particle_set(iatom)%r
     171          296 :                   ria = pbc(ria, cell)
     172           74 :                   atomic_kind => particle_set(iatom)%atomic_kind
     173           74 :                   CALL get_atomic_kind(atomic_kind, kind_number=akind)
     174           74 :                   CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
     175          296 :                   rpoint(:) = rpoint(:) + charge*ria(:)
     176           74 :                   IF (PRESENT(drpoint)) drpoint = drpoint + charge*particle_set(iatom)%v
     177          242 :                   ztot = ztot + charge
     178              :                END DO
     179              :             END DO
     180           54 :             CALL para_env%sum(rpoint)
     181           54 :             CALL para_env%sum(ztot)
     182              :          END IF
     183           54 :          IF (ABS(ztot) > 0._dp) THEN
     184          216 :             rpoint(:) = rpoint(:)/ztot
     185           54 :             IF (PRESENT(drpoint)) drpoint = drpoint/ztot
     186              :          END IF
     187              :       CASE (use_mom_ref_user)
     188           16 :          rpoint = ref_point
     189              :       CASE (use_mom_ref_zero)
     190        23819 :          rpoint = 0._dp
     191              :       END SELECT
     192              : 
     193        23819 :    END SUBROUTINE get_reference_point
     194              : 
     195              : END MODULE moments_utils
     196              : 
        

Generated by: LCOV version 2.0-1