Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 23165 : 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 23165 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
78 23165 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
79 :
80 0 : CPASSERT(PRESENT(ifirst) .EQV. PRESENT(ilast))
81 23165 : NULLIFY (cell, particle_set, qs_kind_set, local_particles, para_env)
82 23165 : 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 2063 : local_particles=local_particles, para_env=para_env)
86 : END IF
87 23165 : IF (PRESENT(fist_env)) THEN
88 : CALL fist_env_get(fist_env, cell=cell, particle_set=particle_set, &
89 21102 : local_particles=local_particles, para_env=para_env)
90 : END IF
91 23165 : do_molecule = .FALSE.
92 23165 : IF (PRESENT(ifirst) .AND. PRESENT(ilast)) do_molecule = .TRUE.
93 23165 : IF (PRESENT(drpoint)) drpoint = 0.0_dp
94 23165 : SELECT CASE (reference)
95 : CASE DEFAULT
96 0 : CPABORT("Type of reference point not implemented")
97 : CASE (use_mom_ref_com)
98 1172 : rpoint = 0._dp
99 1172 : mtot = 0._dp
100 1172 : center(:) = 0._dp
101 1172 : 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 3832 : DO ikind = 1, SIZE(local_particles%n_el)
123 5901 : DO ia = 1, local_particles%n_el(ikind)
124 2069 : iatom = local_particles%list(ikind)%array(ia)
125 8276 : ria = particle_set(iatom)%r
126 8276 : ria = pbc(ria, cell)
127 2069 : atomic_kind => particle_set(iatom)%atomic_kind
128 2069 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
129 8276 : rpoint(:) = rpoint(:) + mass*ria(:)
130 4157 : IF (PRESENT(drpoint)) drpoint = drpoint + mass*particle_set(iatom)%v
131 6804 : mtot = mtot + mass
132 : END DO
133 : END DO
134 1166 : CALL para_env%sum(rpoint)
135 1166 : CALL para_env%sum(mtot)
136 : END IF
137 1172 : IF (ABS(mtot) > 0._dp) THEN
138 4688 : rpoint(:) = rpoint(:)/mtot
139 2324 : IF (PRESENT(drpoint)) drpoint = drpoint/mtot
140 : END IF
141 : CASE (use_mom_ref_coac)
142 6 : rpoint = 0._dp
143 6 : ztot = 0._dp
144 6 : center(:) = 0._dp
145 6 : 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 18 : DO ikind = 1, SIZE(local_particles%n_el)
168 27 : DO ia = 1, local_particles%n_el(ikind)
169 9 : iatom = local_particles%list(ikind)%array(ia)
170 36 : ria = particle_set(iatom)%r
171 36 : ria = pbc(ria, cell)
172 9 : atomic_kind => particle_set(iatom)%atomic_kind
173 9 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
174 9 : CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
175 36 : rpoint(:) = rpoint(:) + charge*ria(:)
176 9 : IF (PRESENT(drpoint)) drpoint = drpoint + charge*particle_set(iatom)%v
177 30 : ztot = ztot + charge
178 : END DO
179 : END DO
180 6 : CALL para_env%sum(rpoint)
181 6 : CALL para_env%sum(ztot)
182 : END IF
183 6 : IF (ABS(ztot) > 0._dp) THEN
184 24 : rpoint(:) = rpoint(:)/ztot
185 6 : 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 23165 : rpoint = 0._dp
191 : END SELECT
192 :
193 23165 : END SUBROUTINE get_reference_point
194 :
195 : END MODULE moments_utils
196 :
|
