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 Set of routines handling the localization for molecular properties
10 : ! **************************************************************************************************
11 : MODULE molecular_moments
12 : USE atomic_kind_types, ONLY: atomic_kind_type,&
13 : get_atomic_kind
14 : USE cell_types, ONLY: cell_type,&
15 : pbc
16 : USE cp_control_types, ONLY: dft_control_type
17 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
18 : dbcsr_deallocate_matrix,&
19 : dbcsr_p_type,&
20 : dbcsr_set
21 : USE cp_dbcsr_operations, ONLY: cp_dbcsr_sm_fm_multiply
22 : USE cp_fm_basic_linalg, ONLY: cp_fm_schur_product
23 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
24 : cp_fm_struct_release,&
25 : cp_fm_struct_type
26 : USE cp_fm_types, ONLY: cp_fm_create,&
27 : cp_fm_get_info,&
28 : cp_fm_release,&
29 : cp_fm_to_fm,&
30 : cp_fm_type
31 : USE cp_log_handling, ONLY: cp_get_default_logger,&
32 : cp_logger_type
33 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
34 : cp_print_key_unit_nr
35 : USE distribution_1d_types, ONLY: distribution_1d_type
36 : USE input_constants, ONLY: use_mom_ref_com
37 : USE input_section_types, ONLY: section_vals_type,&
38 : section_vals_val_get
39 : USE kinds, ONLY: dp
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE molecule_kind_types, ONLY: get_molecule_kind,&
42 : molecule_kind_type
43 : USE molecule_types, ONLY: molecule_type
44 : USE moments_utils, ONLY: get_reference_point
45 : USE orbital_pointers, ONLY: indco,&
46 : ncoset
47 : USE particle_types, ONLY: particle_type
48 : USE qs_environment_types, ONLY: get_qs_env,&
49 : qs_environment_type
50 : USE qs_kind_types, ONLY: get_qs_kind,&
51 : qs_kind_type
52 : USE qs_loc_types, ONLY: qs_loc_env_type
53 : USE qs_moments, ONLY: build_local_moment_matrix
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : ! *** Public ***
61 : PUBLIC :: calculate_molecular_moments
62 :
63 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'molecular_moments'
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief Calculates electrical molecular moments using local operators r-r_ref
69 : !> r_ref: center of mass of the molecule
70 : !> Output is in atomic units
71 : !> \param qs_env the qs_env in which the qs_env lives
72 : !> \param qs_loc_env ...
73 : !> \param mo_local ...
74 : !> \param loc_print_key ...
75 : !> \param molecule_set ...
76 : ! **************************************************************************************************
77 2 : SUBROUTINE calculate_molecular_moments(qs_env, qs_loc_env, mo_local, loc_print_key, molecule_set)
78 : TYPE(qs_environment_type), POINTER :: qs_env
79 : TYPE(qs_loc_env_type), INTENT(IN) :: qs_loc_env
80 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: mo_local
81 : TYPE(section_vals_type), POINTER :: loc_print_key
82 : TYPE(molecule_type), POINTER :: molecule_set(:)
83 :
84 : INTEGER :: akind, first_atom, i, iatom, ikind, imol, imol_now, iounit, iproc, ispin, j, lx, &
85 : ly, lz, molkind, n, n1, n2, natom, ncol_global, nm, nmol, nmols, norder, nrow_global, ns, &
86 : nspins
87 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: states
88 : INTEGER, DIMENSION(2) :: nstates
89 : LOGICAL :: floating, ghost
90 : REAL(KIND=dp) :: zeff, zmom, zwfc
91 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: charge_set
92 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: moment_set
93 : REAL(KIND=dp), DIMENSION(3) :: rcc, ria
94 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: ref_point
95 : TYPE(atomic_kind_type), POINTER :: atomic_kind
96 : TYPE(cell_type), POINTER :: cell
97 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
98 : TYPE(cp_fm_type) :: momv, mvector, omvector
99 : TYPE(cp_logger_type), POINTER :: logger
100 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
101 : TYPE(dft_control_type), POINTER :: dft_control
102 : TYPE(distribution_1d_type), POINTER :: local_molecules
103 : TYPE(molecule_kind_type), POINTER :: molecule_kind
104 : TYPE(mp_para_env_type), POINTER :: para_env
105 2 : TYPE(particle_type), POINTER :: particle_set(:)
106 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
107 :
108 4 : logger => cp_get_default_logger()
109 :
110 2 : CALL get_qs_env(qs_env, dft_control=dft_control)
111 2 : nspins = dft_control%nspins
112 2 : zwfc = 3.0_dp - REAL(nspins, KIND=dp)
113 :
114 2 : CALL section_vals_val_get(loc_print_key, "MOLECULAR_MOMENTS%ORDER", i_val=norder)
115 2 : CPASSERT(norder >= 0)
116 2 : nm = ncoset(norder) - 1
117 :
118 2 : CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell)
119 2 : particle_set => qs_loc_env%particle_set
120 2 : para_env => qs_loc_env%para_env
121 2 : local_molecules => qs_loc_env%local_molecules
122 2 : molkind = SIZE(local_molecules%n_el)
123 2 : nmols = SIZE(molecule_set)
124 12 : ALLOCATE (charge_set(nmols), moment_set(nm, nmols))
125 6 : charge_set = 0.0_dp
126 42 : moment_set = 0.0_dp
127 :
128 2 : IF (norder > 0) THEN
129 2 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
130 6 : DO imol = 1, SIZE(molecule_set)
131 4 : molecule_kind => molecule_set(imol)%molecule_kind
132 4 : first_atom = molecule_set(imol)%first_atom
133 4 : CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom)
134 : ! Get reference point for this molecule
135 : CALL get_reference_point(rcc, qs_env=qs_env, reference=use_mom_ref_com, &
136 : ref_point=ref_point, ifirst=first_atom, &
137 4 : ilast=first_atom + natom - 1)
138 48 : ALLOCATE (moments(nm))
139 40 : DO i = 1, nm
140 36 : ALLOCATE (moments(i)%matrix)
141 36 : CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, 'MOM MAT')
142 40 : CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
143 : END DO
144 : !
145 4 : CALL build_local_moment_matrix(qs_env, moments, norder, rcc)
146 : !
147 12 : DO ispin = 1, nspins
148 8 : IF (ASSOCIATED(molecule_set(imol)%lmi)) THEN
149 4 : nstates(1) = molecule_set(imol)%lmi(ispin)%nstates
150 : ELSE
151 4 : nstates(1) = 0
152 : END IF
153 8 : nstates(2) = para_env%mepos
154 8 : CALL para_env%maxloc(nstates)
155 8 : IF (nstates(1) == 0) CYCLE
156 8 : ns = nstates(1)
157 8 : iproc = nstates(2)
158 24 : ALLOCATE (states(ns))
159 8 : IF (iproc == para_env%mepos) THEN
160 20 : states(:) = molecule_set(imol)%lmi(ispin)%states(:)
161 : ELSE
162 20 : states(:) = 0
163 : END IF
164 8 : CALL para_env%bcast(states, iproc)
165 : ! assemble local states for this molecule
166 : ASSOCIATE (mo_localized => mo_local(ispin))
167 8 : CALL cp_fm_get_info(mo_localized, ncol_global=ncol_global, nrow_global=nrow_global)
168 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_global, ncol_global=ns, &
169 : para_env=mo_localized%matrix_struct%para_env, &
170 8 : context=mo_localized%matrix_struct%context)
171 8 : CALL cp_fm_create(mvector, fm_struct, name="mvector")
172 8 : CALL cp_fm_create(omvector, fm_struct, name="omvector")
173 8 : CALL cp_fm_create(momv, fm_struct, name="omvector")
174 8 : CALL cp_fm_struct_release(fm_struct)
175 : !
176 40 : DO i = 1, ns
177 40 : CALL cp_fm_to_fm(mo_localized, mvector, 1, states(i), i)
178 : END DO
179 : END ASSOCIATE
180 80 : DO i = 1, nm
181 72 : CALL cp_dbcsr_sm_fm_multiply(moments(i)%matrix, mvector, omvector, ns)
182 72 : CALL cp_fm_schur_product(mvector, omvector, momv)
183 2096 : moment_set(i, imol) = moment_set(i, imol) - zwfc*SUM(momv%local_data)
184 : END DO
185 : !
186 8 : CALL cp_fm_release(mvector)
187 8 : CALL cp_fm_release(omvector)
188 8 : CALL cp_fm_release(momv)
189 28 : DEALLOCATE (states)
190 : END DO
191 40 : DO i = 1, nm
192 40 : CALL dbcsr_deallocate_matrix(moments(i)%matrix)
193 : END DO
194 10 : DEALLOCATE (moments)
195 : END DO
196 : END IF
197 : !
198 6 : DO ikind = 1, molkind ! loop over different molecules
199 4 : nmol = SIZE(local_molecules%list(ikind)%array)
200 8 : DO imol = 1, nmol ! all the molecules of the kind
201 2 : imol_now = local_molecules%list(ikind)%array(imol) ! index in the global array
202 2 : molecule_kind => molecule_set(imol_now)%molecule_kind
203 2 : first_atom = molecule_set(imol_now)%first_atom
204 2 : CALL get_molecule_kind(molecule_kind=molecule_kind, natom=natom)
205 : ! Get reference point for this molecule
206 : CALL get_reference_point(rcc, qs_env=qs_env, reference=use_mom_ref_com, &
207 : ref_point=ref_point, ifirst=first_atom, &
208 2 : ilast=first_atom + natom - 1)
209 : ! charge
210 8 : DO iatom = 1, natom
211 6 : i = first_atom + iatom - 1
212 6 : atomic_kind => particle_set(i)%atomic_kind
213 6 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
214 6 : CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
215 14 : IF (.NOT. ghost .AND. .NOT. floating) THEN
216 6 : CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
217 6 : charge_set(imol_now) = charge_set(imol_now) + zeff
218 : END IF
219 : END DO
220 6 : DO ispin = 1, nspins
221 6 : IF (ASSOCIATED(molecule_set(imol_now)%lmi(ispin)%states)) THEN
222 4 : ns = SIZE(molecule_set(imol_now)%lmi(ispin)%states)
223 4 : charge_set(imol_now) = charge_set(imol_now) - zwfc*ns
224 : END IF
225 : END DO
226 : !
227 8 : IF (norder > 0) THEN
228 : ! nuclear contribution
229 20 : DO i = 1, nm
230 18 : lx = indco(1, i + 1)
231 18 : ly = indco(2, i + 1)
232 18 : lz = indco(3, i + 1)
233 74 : DO iatom = 1, natom
234 54 : j = first_atom + iatom - 1
235 54 : atomic_kind => particle_set(j)%atomic_kind
236 54 : CALL get_atomic_kind(atomic_kind, kind_number=akind)
237 54 : CALL get_qs_kind(qs_kind_set(akind), ghost=ghost, floating=floating)
238 126 : IF (.NOT. ghost .AND. .NOT. floating) THEN
239 54 : CALL get_qs_kind(qs_kind_set(akind), core_charge=zeff)
240 216 : ria = particle_set(j)%r - rcc
241 216 : ria = pbc(ria, cell)
242 54 : zmom = zeff
243 54 : IF (lx /= 0) zmom = zmom*ria(1)**lx
244 54 : IF (ly /= 0) zmom = zmom*ria(2)**ly
245 54 : IF (lz /= 0) zmom = zmom*ria(3)**lz
246 54 : moment_set(i, imol_now) = moment_set(i, imol_now) + zmom
247 : END IF
248 : END DO
249 : END DO
250 : END IF
251 : END DO
252 : END DO
253 2 : CALL para_env%sum(moment_set)
254 2 : CALL para_env%sum(charge_set)
255 :
256 : iounit = cp_print_key_unit_nr(logger, loc_print_key, "MOLECULAR_MOMENTS", &
257 2 : extension=".MolMom", middle_name="MOLECULAR_MOMENTS")
258 2 : IF (iounit > 0) THEN
259 3 : DO i = 1, SIZE(charge_set)
260 2 : WRITE (UNIT=iounit, FMT='(A,I6,A,F12.6)') " # molecule nr:", i, " Charge:", charge_set(I)
261 7 : DO n = 1, norder
262 4 : n1 = ncoset(n - 1)
263 4 : n2 = ncoset(n) - 1
264 6 : WRITE (UNIT=iounit, FMT='(T4,A,I2,10(T16,6F12.6))') "Order:", n, moment_set(n1:n2, i)
265 : END DO
266 : END DO
267 : END IF
268 : CALL cp_print_key_finished_output(iounit, logger, loc_print_key, &
269 2 : "MOLECULAR_MOMENTS")
270 :
271 2 : DEALLOCATE (charge_set, moment_set)
272 :
273 6 : END SUBROUTINE calculate_molecular_moments
274 : !------------------------------------------------------------------------------
275 :
276 : END MODULE molecular_moments
277 :
|
