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