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