Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 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 23817 : 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(:), OPTIONAL, 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 23817 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
78 23817 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
79 :
80 0 : CPASSERT(PRESENT(ifirst) .EQV. PRESENT(ilast))
81 23817 : NULLIFY (cell, particle_set, qs_kind_set, local_particles, para_env)
82 23817 : 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 2613 : local_particles=local_particles, para_env=para_env)
86 : END IF
87 23817 : IF (PRESENT(fist_env)) THEN
88 : CALL fist_env_get(fist_env, cell=cell, particle_set=particle_set, &
89 21204 : local_particles=local_particles, para_env=para_env)
90 : END IF
91 23817 : do_molecule = .FALSE.
92 23817 : IF (PRESENT(ifirst) .AND. PRESENT(ilast)) do_molecule = .TRUE.
93 23817 : IF (PRESENT(drpoint)) drpoint = 0.0_dp
94 23817 : IF (reference == use_mom_ref_user .AND. (.NOT. PRESENT(ref_point))) THEN
95 0 : CPABORT("User reference point not supplied!")
96 : END IF
97 0 : SELECT CASE (reference)
98 : CASE DEFAULT
99 0 : CPABORT("Type of reference point not implemented")
100 : CASE (use_mom_ref_com)
101 1484 : rpoint = 0._dp
102 1484 : mtot = 0._dp
103 1484 : center(:) = 0._dp
104 1484 : IF (do_molecule) THEN
105 6 : mass_low = -HUGE(mass_low)
106 : ! fold the molecule around the heaviest atom in the molecule
107 24 : DO iatom = ifirst, ilast
108 18 : atomic_kind => particle_set(iatom)%atomic_kind
109 18 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
110 24 : IF (mass > mass_low) THEN
111 6 : mass_low = mass
112 24 : center = particle_set(iatom)%r
113 : END IF
114 : END DO
115 24 : DO iatom = ifirst, ilast
116 72 : ria = particle_set(iatom)%r
117 126 : ria = pbc(ria - center, cell) + center
118 18 : atomic_kind => particle_set(iatom)%atomic_kind
119 18 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
120 72 : rpoint(:) = rpoint(:) + mass*ria(:)
121 18 : IF (PRESENT(drpoint)) drpoint = drpoint + mass*particle_set(iatom)%v
122 42 : mtot = mtot + mass
123 : END DO
124 : ELSE
125 5080 : DO ikind = 1, SIZE(local_particles%n_el)
126 7773 : DO ia = 1, local_particles%n_el(ikind)
127 2693 : iatom = local_particles%list(ikind)%array(ia)
128 10772 : ria = particle_set(iatom)%r
129 10772 : ria = pbc(ria, cell)
130 2693 : atomic_kind => particle_set(iatom)%atomic_kind
131 2693 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
132 10772 : rpoint(:) = rpoint(:) + mass*ria(:)
133 6653 : IF (PRESENT(drpoint)) drpoint = drpoint + mass*particle_set(iatom)%v
134 8988 : mtot = mtot + mass
135 : END DO
136 : END DO
137 1478 : CALL para_env%sum(rpoint)
138 1478 : CALL para_env%sum(mtot)
139 : END IF
140 1484 : IF (ABS(mtot) > 0._dp) THEN
141 5936 : rpoint(:) = rpoint(:)/mtot
142 3884 : IF (PRESENT(drpoint)) drpoint = drpoint/mtot
143 : END IF
144 : CASE (use_mom_ref_coac)
145 72 : rpoint = 0._dp
146 72 : ztot = 0._dp
147 72 : center(:) = 0._dp
148 72 : IF (do_molecule) THEN
149 0 : mass_low = -HUGE(mass_low)
150 : ! fold the molecule around the heaviest atom in the molecule
151 0 : DO iatom = ifirst, ilast
152 0 : atomic_kind => particle_set(iatom)%atomic_kind
153 0 : CALL get_atomic_kind(atomic_kind=atomic_kind, mass=mass)
154 0 : IF (mass > mass_low) THEN
155 0 : mass_low = mass
156 0 : center = particle_set(iatom)%r
157 : END IF
158 : END DO
159 0 : DO iatom = ifirst, ilast
160 0 : ria = particle_set(iatom)%r
161 0 : ria = pbc(ria - center, cell) + center
162 0 : atomic_kind => particle_set(iatom)%atomic_kind
163 0 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
164 0 : CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
165 0 : rpoint(:) = rpoint(:) + charge*ria(:)
166 0 : IF (PRESENT(drpoint)) drpoint = drpoint + charge*particle_set(iatom)%v
167 0 : ztot = ztot + charge
168 : END DO
169 : ELSE
170 202 : DO ikind = 1, SIZE(local_particles%n_el)
171 303 : DO ia = 1, local_particles%n_el(ikind)
172 101 : iatom = local_particles%list(ikind)%array(ia)
173 404 : ria = particle_set(iatom)%r
174 404 : ria = pbc(ria, cell)
175 101 : atomic_kind => particle_set(iatom)%atomic_kind
176 101 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
177 101 : CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
178 404 : rpoint(:) = rpoint(:) + charge*ria(:)
179 101 : IF (PRESENT(drpoint)) drpoint = drpoint + charge*particle_set(iatom)%v
180 332 : ztot = ztot + charge
181 : END DO
182 : END DO
183 72 : CALL para_env%sum(rpoint)
184 72 : CALL para_env%sum(ztot)
185 : END IF
186 72 : IF (ABS(ztot) > 0._dp) THEN
187 288 : rpoint(:) = rpoint(:)/ztot
188 72 : IF (PRESENT(drpoint)) drpoint = drpoint/ztot
189 : END IF
190 : CASE (use_mom_ref_user)
191 16 : rpoint = ref_point
192 : CASE (use_mom_ref_zero)
193 23817 : rpoint = 0._dp
194 : END SELECT
195 :
196 23817 : END SUBROUTINE get_reference_point
197 :
198 : END MODULE moments_utils
199 :
|
