File: moments_utils.F

package info (click to toggle)
cp2k 6.1-2
  • links: PTS, VCS
  • area: main
  • in suites: buster
  • size: 204,532 kB
  • sloc: fortran: 835,196; f90: 59,605; python: 9,861; sh: 7,882; cpp: 4,868; ansic: 2,807; xml: 2,185; lisp: 733; pascal: 612; perl: 547; makefile: 497; csh: 16
file content (196 lines) | stat: -rw-r--r-- 9,058 bytes parent folder | download
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