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