Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculate ntrinsic atomic orbitals and analyze wavefunctions
10 : !> \par History
11 : !> 03.2023 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE iao_types
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind_set
17 : USE basis_set_types, ONLY: gto_basis_set_p_type
18 : USE cell_types, ONLY: cell_type,&
19 : pbc
20 : USE cp_blacs_env, ONLY: cp_blacs_env_type
21 : USE cp_dbcsr_api, ONLY: dbcsr_get_info,&
22 : dbcsr_p_type,&
23 : dbcsr_type
24 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm
25 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
26 : cp_fm_struct_release,&
27 : cp_fm_struct_type
28 : USE cp_fm_types, ONLY: cp_fm_create,&
29 : cp_fm_release,&
30 : cp_fm_to_fm_submat,&
31 : cp_fm_type,&
32 : cp_fm_write_formatted
33 : USE cp_log_handling, ONLY: cp_get_default_logger,&
34 : cp_logger_type
35 : USE cp_output_handling, ONLY: cp_p_file,&
36 : cp_print_key_finished_output,&
37 : cp_print_key_should_output,&
38 : cp_print_key_unit_nr
39 : USE cp_units, ONLY: cp_unit_from_cp2k
40 : USE input_constants, ONLY: do_iaoloc_enone,&
41 : do_iaoloc_pm2
42 : USE input_section_types, ONLY: section_vals_get,&
43 : section_vals_get_subs_vals,&
44 : section_vals_type,&
45 : section_vals_val_get
46 : USE kinds, ONLY: dp
47 : USE message_passing, ONLY: mp_para_env_type
48 : USE particle_types, ONLY: particle_type
49 : USE qs_environment_types, ONLY: get_qs_env,&
50 : qs_environment_type
51 : #include "./base/base_uses.f90"
52 :
53 : IMPLICIT NONE
54 : PRIVATE
55 :
56 : PUBLIC :: iao_env_type, iao_read_input, iao_set_default, &
57 : organise_printout_for_rose, print_fragment_association
58 :
59 : ! **************************************************************************************************
60 : TYPE iao_env_type
61 : LOGICAL :: do_iao = .FALSE.
62 : !
63 : REAL(KIND=dp) :: eps_svd = 0.0_dp
64 : REAL(KIND=dp) :: eps_occ = 0.0_dp
65 : ! chages
66 : LOGICAL :: do_charges = .FALSE.
67 : ! one-center expansion
68 : LOGICAL :: do_oce = .FALSE.
69 : INTEGER :: lmax_oce = 0
70 : INTEGER :: nbas_oce = 0
71 : LOGICAL :: append_oce = .FALSE.
72 : ! Bond orbitals
73 : LOGICAL :: do_bondorbitals = .FALSE.
74 : ! Wannier centers
75 : LOGICAL :: do_center = .FALSE.
76 : LOGICAL :: pos_periodic = .FALSE.
77 : INTEGER :: loc_operator = 0
78 : INTEGER :: eloc_function = 0
79 : REAL(KIND=dp) :: eloc_weight = 0.0_dp
80 : ! Molden
81 : LOGICAL :: molden_iao = .FALSE.
82 : LOGICAL :: molden_ibo = .FALSE.
83 : ! CUBE files
84 : LOGICAL :: cubes_iao = .FALSE.
85 : LOGICAL :: cubes_ibo = .FALSE.
86 : ! Input sections
87 : TYPE(section_vals_type), POINTER :: iao_cubes_section => NULL(), &
88 : iao_molden_section => NULL(), &
89 : ibo_cubes_section => NULL(), &
90 : ibo_molden_section => NULL(), &
91 : ibo_cc_section => NULL()
92 : ! Fragments
93 : LOGICAL :: do_fragments = .FALSE.
94 : INTEGER :: nfrags = 0
95 : END TYPE iao_env_type
96 :
97 : ! **************************************************************************************************
98 :
99 : CONTAINS
100 :
101 : ! **************************************************************************************************
102 : !> \brief ...
103 : !> \param iao_env ...
104 : ! **************************************************************************************************
105 64 : SUBROUTINE iao_set_default(iao_env)
106 : TYPE(iao_env_type), INTENT(INOUT) :: iao_env
107 :
108 : !iao
109 64 : iao_env%do_iao = .FALSE.
110 64 : iao_env%eps_svd = 0.0_dp
111 64 : iao_env%eps_occ = 0.0_dp
112 : ! charges
113 64 : iao_env%do_charges = .FALSE.
114 : ! one-center expansion
115 64 : iao_env%do_oce = .FALSE.
116 64 : iao_env%lmax_oce = 3
117 64 : iao_env%nbas_oce = 10
118 64 : iao_env%append_oce = .FALSE.
119 : ! Bond orbitals
120 64 : iao_env%do_bondorbitals = .FALSE.
121 : ! Wannier centers
122 64 : iao_env%do_center = .FALSE.
123 64 : iao_env%pos_periodic = .FALSE.
124 64 : iao_env%loc_operator = do_iaoloc_pm2
125 64 : iao_env%eloc_function = do_iaoloc_enone
126 64 : iao_env%eloc_weight = 0.0_dp
127 : ! i/o
128 64 : iao_env%molden_iao = .FALSE.
129 64 : iao_env%molden_ibo = .FALSE.
130 64 : iao_env%cubes_iao = .FALSE.
131 64 : iao_env%cubes_ibo = .FALSE.
132 : ! Input sections
133 64 : NULLIFY (iao_env%iao_cubes_section, iao_env%iao_molden_section)
134 64 : NULLIFY (iao_env%ibo_cubes_section, iao_env%ibo_molden_section)
135 64 : NULLIFY (iao_env%ibo_cc_section)
136 :
137 64 : END SUBROUTINE iao_set_default
138 :
139 : ! **************************************************************************************************
140 :
141 : ! **************************************************************************************************
142 : !> \brief ...
143 : !> \param iao_env ...
144 : !> \param iao_section ...
145 : !> \param cell ...
146 : ! **************************************************************************************************
147 68 : SUBROUTINE iao_read_input(iao_env, iao_section, cell)
148 : TYPE(iao_env_type), INTENT(INOUT) :: iao_env
149 : TYPE(section_vals_type), POINTER :: iao_section
150 : TYPE(cell_type), OPTIONAL :: cell
151 :
152 : LOGICAL :: explicit, iao_explicit
153 : TYPE(section_vals_type), POINTER :: subsection
154 :
155 34 : CALL iao_set_default(iao_env)
156 :
157 34 : CALL section_vals_get(iao_section, explicit=iao_explicit)
158 34 : IF (iao_explicit) THEN
159 6 : iao_env%do_iao = .TRUE.
160 : ! input options
161 6 : CALL section_vals_val_get(iao_section, "EPS_SVD", r_val=iao_env%eps_svd)
162 6 : CALL section_vals_val_get(iao_section, "EPS_OCC", r_val=iao_env%eps_occ)
163 6 : CALL section_vals_val_get(iao_section, "ATOMIC_CHARGES", l_val=iao_env%do_charges)
164 6 : iao_env%iao_molden_section => section_vals_get_subs_vals(iao_section, "IAO_MOLDEN")
165 6 : CALL section_vals_get(iao_env%iao_molden_section, explicit=iao_env%molden_iao)
166 6 : iao_env%iao_cubes_section => section_vals_get_subs_vals(iao_section, "IAO_CUBES")
167 6 : CALL section_vals_get(iao_env%iao_cubes_section, explicit=iao_env%cubes_iao)
168 6 : subsection => section_vals_get_subs_vals(iao_section, "ONE_CENTER_EXPANSION")
169 6 : CALL section_vals_get(subsection, explicit=iao_env%do_oce)
170 6 : IF (iao_env%do_oce) THEN
171 4 : subsection => section_vals_get_subs_vals(iao_section, "ONE_CENTER_EXPANSION")
172 4 : CALL section_vals_val_get(subsection, "LMAX", i_val=iao_env%lmax_oce)
173 4 : CALL section_vals_val_get(subsection, "NBAS", i_val=iao_env%nbas_oce)
174 4 : CALL section_vals_val_get(subsection, "APPEND", l_val=iao_env%append_oce)
175 : END IF
176 6 : subsection => section_vals_get_subs_vals(iao_section, "BOND_ORBITALS")
177 6 : CALL section_vals_get(subsection, explicit=iao_env%do_bondorbitals)
178 6 : IF (iao_env%do_bondorbitals) THEN
179 6 : subsection => section_vals_get_subs_vals(iao_section, "BOND_ORBITALS")
180 6 : CALL section_vals_val_get(subsection, "LOCALIZATION_OPERATOR", i_val=iao_env%loc_operator)
181 6 : CALL section_vals_val_get(subsection, "ENERGY_LOCALIZATION_FUNCTION", i_val=iao_env%eloc_function)
182 6 : CALL section_vals_val_get(subsection, "ENERGY_LOCALIZATION_WEIGHT", r_val=iao_env%eloc_weight)
183 6 : iao_env%ibo_molden_section => section_vals_get_subs_vals(subsection, "IBO_MOLDEN")
184 6 : CALL section_vals_get(iao_env%ibo_molden_section, explicit=iao_env%molden_ibo)
185 6 : iao_env%ibo_cubes_section => section_vals_get_subs_vals(subsection, "IBO_CUBES")
186 6 : CALL section_vals_get(iao_env%ibo_cubes_section, explicit=iao_env%cubes_ibo)
187 6 : iao_env%ibo_cc_section => section_vals_get_subs_vals(subsection, "CHARGE_CENTER")
188 6 : CALL section_vals_get(iao_env%ibo_cc_section, explicit=iao_env%do_center)
189 6 : IF (iao_env%do_center) THEN
190 : CALL section_vals_val_get(iao_env%ibo_cc_section, "POSITION_OPERATOR_BERRY", &
191 6 : l_val=iao_env%pos_periodic, explicit=explicit)
192 6 : IF (.NOT. explicit) THEN
193 : ! set default according to cell periodicity
194 2 : iao_env%pos_periodic = .TRUE.
195 2 : IF (PRESENT(cell)) THEN
196 8 : IF (ALL(cell%perd == 0)) iao_env%pos_periodic = .FALSE.
197 : END IF
198 : END IF
199 : END IF
200 : END IF
201 : END IF
202 :
203 34 : END SUBROUTINE iao_read_input
204 : ! **************************************************************************************************
205 : !> \brief Prints fragmented overlap matrices for interface with ROSE
206 : !> \param qs_env ...
207 : !> \param orb_basis_set_list ...
208 : ! **************************************************************************************************
209 2 : SUBROUTINE organise_printout_for_rose(qs_env, orb_basis_set_list)
210 : TYPE(qs_environment_type), POINTER :: qs_env
211 : TYPE(gto_basis_set_p_type), DIMENSION(:), &
212 : INTENT(IN), POINTER :: orb_basis_set_list
213 :
214 : CHARACTER(len=*), PARAMETER :: routineN = 'organise_printout_for_rose'
215 :
216 : INTEGER :: after, colskip, handle, iatom, ii, &
217 : ikind, iounit, jj, ncol, nfrags, nrow, &
218 : rowfrag, rowskip, totalcol, totalrow
219 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of
220 : LOGICAL :: lfirstcol, lfirstcoljj
221 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
222 : TYPE(cp_blacs_env_type), POINTER :: blacs_env
223 : TYPE(cp_fm_struct_type), POINTER :: fmstruct
224 : TYPE(cp_fm_type) :: fm_sorb, fm_sorb_frags12, fm_sorb_frags22
225 : TYPE(cp_logger_type), POINTER :: logger
226 2 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
227 : TYPE(dbcsr_type), POINTER :: smat
228 : TYPE(mp_para_env_type), POINTER :: para_env
229 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
230 :
231 2 : CALL timeset(routineN, handle)
232 :
233 : !get number of fragments
234 2 : NULLIFY (particle_set, atomic_kind_set)
235 : CALL get_qs_env(qs_env, particle_set=particle_set, atomic_kind_set=atomic_kind_set, &
236 2 : matrix_s_kp=matrix_s)
237 2 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
238 2 : smat => matrix_s(1, 1)%matrix
239 :
240 2 : nfrags = 0
241 16 : DO ii = 1, SIZE(particle_set)
242 14 : IF (nfrags < particle_set(ii)%fragment_index) &
243 2 : nfrags = particle_set(ii)%fragment_index
244 : END DO
245 :
246 2 : logger => cp_get_default_logger()
247 2 : IF (BTEST(cp_print_key_should_output(logger%iter_info, &
248 : qs_env%input, "DFT%PRINT%IAO_ANALYSIS%IAO_OVERLAP"), cp_p_file)) THEN
249 : iounit = cp_print_key_unit_nr(logger, qs_env%input, "DFT%PRINT%IAO_ANALYSIS%IAO_OVERLAP", &
250 : extension=".Log", file_form="FORMATTED", file_action="WRITE", &
251 2 : file_status="REPLACE")
252 2 : CALL section_vals_val_get(qs_env%input, "DFT%PRINT%IAO_ANALYSIS%IAO_OVERLAP%NDIGITS", i_val=after)
253 : after = MIN(MAX(after, 1), 16)
254 :
255 6 : DO ii = 1, nfrags
256 14 : DO jj = 1, nfrags
257 8 : NULLIFY (para_env, blacs_env)
258 8 : CALL get_qs_env(qs_env=qs_env, blacs_env=blacs_env, para_env=para_env)
259 8 : totalcol = 0
260 8 : totalrow = 0
261 8 : rowfrag = 0
262 8 : colskip = 0
263 8 : rowskip = 0
264 8 : lfirstcol = .TRUE.
265 8 : lfirstcoljj = .TRUE.
266 :
267 64 : DO iatom = 1, SIZE(particle_set)
268 56 : ikind = kind_of(iatom)
269 56 : IF (particle_set(iatom)%fragment_index == jj) THEN
270 28 : totalcol = totalcol + orb_basis_set_list(ikind)%gto_basis_set%nsgf
271 28 : lfirstcoljj = .FALSE.
272 28 : ELSEIF (lfirstcoljj) THEN
273 16 : colskip = colskip + orb_basis_set_list(ikind)%gto_basis_set%nsgf
274 : END IF
275 56 : IF (particle_set(iatom)%fragment_index == ii) THEN
276 28 : totalrow = totalrow + orb_basis_set_list(ikind)%gto_basis_set%nsgf
277 28 : lfirstcol = .FALSE.
278 28 : ELSEIF (lfirstcol) THEN
279 16 : rowskip = rowskip + orb_basis_set_list(ikind)%gto_basis_set%nsgf
280 : END IF
281 64 : rowfrag = rowfrag + orb_basis_set_list(ikind)%gto_basis_set%nsgf
282 : END DO
283 :
284 8 : CALL dbcsr_get_info(matrix_s(1, 1)%matrix, nfullrows_total=nrow, nfullcols_total=ncol)
285 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, context=blacs_env, &
286 8 : nrow_global=nrow, ncol_global=ncol)
287 8 : CALL cp_fm_create(matrix=fm_sorb, matrix_struct=fmstruct, name="Overlap matrix S_B1")
288 8 : CALL cp_fm_struct_release(fmstruct=fmstruct)
289 8 : CALL copy_dbcsr_to_fm(matrix_s(1, 1)%matrix, fm_sorb)
290 :
291 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, context=blacs_env, &
292 8 : nrow_global=totalrow, ncol_global=totalcol)
293 8 : CALL cp_fm_create(matrix=fm_sorb_frags22, matrix_struct=fmstruct, name="Overlap matrix S_B1^kk' (S22)")
294 8 : CALL cp_fm_struct_release(fmstruct=fmstruct)
295 : CALL cp_fm_to_fm_submat(fm_sorb, fm_sorb_frags22, &
296 8 : totalrow, totalcol, 1 + rowskip, 1 + colskip, 1, 1)
297 :
298 : CALL cp_fm_struct_create(fmstruct=fmstruct, para_env=para_env, context=blacs_env, &
299 8 : nrow_global=rowfrag, ncol_global=totalcol)
300 8 : CALL cp_fm_create(matrix=fm_sorb_frags12, matrix_struct=fmstruct, name="Overlap matrix S_B1^k' (S12)")
301 8 : CALL cp_fm_struct_release(fmstruct=fmstruct)
302 8 : CALL cp_fm_to_fm_submat(fm_sorb, fm_sorb_frags12, rowfrag, totalcol, 1, 1 + colskip, 1, 1)
303 :
304 8 : IF (iounit > 0) WRITE (iounit, *) "FRAGMENTS kk' NR. ", ii, jj
305 8 : CALL cp_fm_write_formatted(fm_sorb_frags22, iounit, "Overlap matrix S_B1^kk' (S22)")
306 8 : IF (ii /= jj) THEN
307 4 : IF (iounit > 0) WRITE (iounit, *) "FRAGMENTS kk' NR. ", ii, jj
308 4 : CALL cp_fm_write_formatted(fm_sorb_frags12, iounit, "Overlap matrix S_B1^k' (S12)")
309 : END IF
310 :
311 8 : CALL cp_fm_release(fm_sorb)
312 8 : CALL cp_fm_release(fm_sorb_frags22)
313 44 : CALL cp_fm_release(fm_sorb_frags12)
314 :
315 : END DO
316 : END DO
317 :
318 2 : CALL cp_print_key_finished_output(iounit, logger, qs_env%input, "DFT%PRINT%IAO_ANALYSIS%IAO_OVERLAP")
319 : END IF
320 :
321 2 : CALL timestop(handle)
322 :
323 4 : END SUBROUTINE organise_printout_for_rose
324 : ! **************************************************************************************************
325 : !> \brief Associates IBO centers to fragments
326 : !> \param qs_env ...
327 : !> \param moments ...
328 : !> \param unit_nr ...
329 : ! **************************************************************************************************
330 2 : SUBROUTINE print_fragment_association(qs_env, moments, unit_nr)
331 : TYPE(qs_environment_type), POINTER :: qs_env
332 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: moments
333 : INTEGER, INTENT(IN), OPTIONAL :: unit_nr
334 :
335 : CHARACTER(len=*), PARAMETER :: routineN = 'print_fragment_association'
336 :
337 : INTEGER :: handle, iatom, is, ispin, jcenter, &
338 : natom, ncenters, nspin
339 2 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: atomcenter
340 : REAL(KIND=dp) :: conv, dab, mindab
341 : REAL(KIND=dp), DIMENSION(3) :: rab
342 : TYPE(cell_type), POINTER :: cell
343 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
344 :
345 2 : CALL timeset(routineN, handle)
346 :
347 2 : NULLIFY (cell, particle_set)
348 2 : CALL get_qs_env(qs_env=qs_env, cell=cell, particle_set=particle_set, natom=natom)
349 2 : conv = cp_unit_from_cp2k(1.0_dp, TRIM("angstrom"))
350 2 : nspin = SIZE(moments, 3)
351 2 : ncenters = SIZE(moments, 2)
352 8 : ALLOCATE (atomcenter(ncenters, nspin))
353 2 : atomcenter = 0
354 4 : DO ispin = 1, nspin
355 22 : DO jcenter = 1, SIZE(moments, 2)
356 20 : mindab = 10.0_dp
357 162 : DO iatom = 1, natom
358 140 : rab(:) = pbc(particle_set(iatom)%r(:), moments(:, jcenter, ispin), cell)
359 140 : dab = SQRT(rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3))
360 160 : IF (dab <= mindab) THEN
361 40 : mindab = dab
362 40 : atomcenter(jcenter, ispin) = iatom
363 : END IF
364 : END DO
365 : END DO
366 :
367 4 : IF (unit_nr > 0) THEN
368 1 : WRITE (unit_nr, "(/,T2,A,i1)") "Intrinsic Bond Orbitals: Centers associated to fragments for Spin ", &
369 2 : ispin
370 1 : WRITE (unit_nr, "(T2,A6,T30,A6,T60,A4,T70,A8)") "Center", "Moment", "Atom", "Fragment"
371 11 : DO is = 1, ncenters
372 10 : WRITE (unit_nr, "(i7,3F15.8,4X,I7,4X,I7)") is, moments(1:3, is, ispin), atomcenter(is, ispin), &
373 21 : particle_set(atomcenter(is, ispin))%fragment_index
374 : END DO
375 1 : WRITE (unit_nr, "(/)")
376 : END IF
377 : END DO
378 :
379 2 : CALL timestop(handle)
380 :
381 4 : END SUBROUTINE print_fragment_association
382 :
383 0 : END MODULE iao_types
|
