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 Calculate intrinsic atomic orbitals and analyze wavefunctions
10 : !> \par History
11 : !> 03.2023 created [JGH]
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE iao_analysis
15 : USE ai_contraction, ONLY: block_add,&
16 : contraction
17 : USE ai_overlap, ONLY: overlap_ab
18 : USE atomic_charges, ONLY: print_atomic_charges
19 : USE atomic_kind_types, ONLY: atomic_kind_type,&
20 : get_atomic_kind
21 : USE auto_basis, ONLY: create_oce_basis
22 : USE basis_set_types, ONLY: deallocate_gto_basis_set,&
23 : get_gto_basis_set,&
24 : gto_basis_set_p_type,&
25 : gto_basis_set_type,&
26 : init_orb_basis_set
27 : USE bibliography, ONLY: Knizia2013,&
28 : cite_reference
29 : USE cell_types, ONLY: cell_type,&
30 : pbc
31 : USE cp_array_utils, ONLY: cp_2d_r_p_type
32 : USE cp_control_types, ONLY: dft_control_type
33 : USE cp_dbcsr_api, ONLY: &
34 : dbcsr_copy, dbcsr_create, dbcsr_distribution_type, dbcsr_get_block_p, dbcsr_get_info, &
35 : dbcsr_iterator_blocks_left, dbcsr_iterator_next_block, dbcsr_iterator_start, &
36 : dbcsr_iterator_stop, dbcsr_iterator_type, dbcsr_multiply, dbcsr_p_type, dbcsr_release, &
37 : dbcsr_replicate_all, dbcsr_set, dbcsr_type, dbcsr_type_no_symmetry
38 : USE cp_dbcsr_contrib, ONLY: dbcsr_reserve_diag_blocks,&
39 : dbcsr_trace
40 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
41 : copy_fm_to_dbcsr,&
42 : cp_dbcsr_plus_fm_fm_t,&
43 : cp_dbcsr_sm_fm_multiply,&
44 : dbcsr_deallocate_matrix_set
45 : USE cp_fm_basic_linalg, ONLY: cp_fm_column_scale,&
46 : cp_fm_geadd,&
47 : cp_fm_invert,&
48 : cp_fm_rot_cols,&
49 : cp_fm_rot_rows
50 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
51 : cp_fm_struct_release,&
52 : cp_fm_struct_type
53 : USE cp_fm_types, ONLY: cp_fm_create,&
54 : cp_fm_get_diag,&
55 : cp_fm_get_element,&
56 : cp_fm_get_info,&
57 : cp_fm_release,&
58 : cp_fm_set_all,&
59 : cp_fm_to_fm,&
60 : cp_fm_type
61 : USE cp_log_handling, ONLY: cp_get_default_logger,&
62 : cp_logger_type
63 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
64 : cp_print_key_unit_nr
65 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
66 : USE iao_types, ONLY: iao_env_type
67 : USE input_constants, ONLY: do_iaoloc_enone,&
68 : do_iaoloc_pm2
69 : USE input_section_types, ONLY: section_get_ivals,&
70 : section_vals_get,&
71 : section_vals_type,&
72 : section_vals_val_get
73 : USE kinds, ONLY: default_path_length,&
74 : default_string_length,&
75 : dp
76 : USE machine, ONLY: m_walltime
77 : USE mathconstants, ONLY: twopi
78 : USE mathlib, ONLY: invmat_symm,&
79 : jacobi
80 : USE message_passing, ONLY: mp_comm_type
81 : USE min_basis_set, ONLY: create_minbas_set
82 : USE molden_utils, ONLY: write_mos_molden
83 : USE orbital_pointers, ONLY: ncoset
84 : USE parallel_gemm_api, ONLY: parallel_gemm
85 : USE particle_list_types, ONLY: particle_list_type
86 : USE particle_methods, ONLY: get_particle_set
87 : USE particle_types, ONLY: particle_type
88 : USE pw_env_types, ONLY: pw_env_get,&
89 : pw_env_type
90 : USE pw_pool_types, ONLY: pw_pool_type
91 : USE pw_types, ONLY: pw_c1d_gs_type,&
92 : pw_r3d_rs_type
93 : USE qs_collocate_density, ONLY: calculate_wavefunction
94 : USE qs_environment_types, ONLY: get_qs_env,&
95 : qs_environment_type
96 : USE qs_interactions, ONLY: init_interaction_radii_orb_basis
97 : USE qs_kind_types, ONLY: get_qs_kind,&
98 : qs_kind_type
99 : USE qs_ks_types, ONLY: get_ks_env,&
100 : qs_ks_env_type
101 : USE qs_loc_utils, ONLY: compute_berry_operator
102 : USE qs_localization_methods, ONLY: initialize_weights,&
103 : rotate_orbitals,&
104 : scdm_qrfact
105 : USE qs_mo_methods, ONLY: make_basis_lowdin
106 : USE qs_mo_types, ONLY: allocate_mo_set,&
107 : deallocate_mo_set,&
108 : get_mo_set,&
109 : init_mo_set,&
110 : mo_set_type,&
111 : set_mo_set
112 : USE qs_moments, ONLY: build_local_moment_matrix
113 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type,&
114 : release_neighbor_list_sets
115 : USE qs_neighbor_lists, ONLY: setup_neighbor_list
116 : USE qs_overlap, ONLY: build_overlap_matrix_simple
117 : USE qs_subsys_types, ONLY: qs_subsys_get,&
118 : qs_subsys_type
119 : #include "./base/base_uses.f90"
120 :
121 : IMPLICIT NONE
122 : PRIVATE
123 :
124 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'iao_analysis'
125 :
126 : PUBLIC :: iao_wfn_analysis, iao_charges, iao_calculate_dmat, print_center_spread
127 :
128 : INTERFACE iao_calculate_dmat
129 : MODULE PROCEDURE iao_calculate_dmat_diag, & ! diagonal occupation matrix
130 : iao_calculate_dmat_full ! full occupation matrix
131 : END INTERFACE
132 :
133 : ! **************************************************************************************************
134 :
135 : CONTAINS
136 :
137 : ! **************************************************************************************************
138 : !> \brief ...
139 : !> \param qs_env ...
140 : !> \param iao_env ...
141 : !> \param unit_nr ...
142 : !> \param c_iao_coef ...
143 : !> \param mos ...
144 : !> \param bond_centers ...
145 : ! **************************************************************************************************
146 2 : SUBROUTINE iao_wfn_analysis(qs_env, iao_env, unit_nr, c_iao_coef, mos, bond_centers)
147 : TYPE(qs_environment_type), POINTER :: qs_env
148 : TYPE(iao_env_type), INTENT(IN) :: iao_env
149 : INTEGER, INTENT(IN) :: unit_nr
150 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:), &
151 : OPTIONAL :: c_iao_coef
152 : TYPE(mo_set_type), DIMENSION(:), OPTIONAL, POINTER :: mos
153 : REAL(KIND=dp), DIMENSION(:, :, :), OPTIONAL :: bond_centers
154 :
155 : CHARACTER(len=*), PARAMETER :: routineN = 'iao_wfn_analysis'
156 :
157 : CHARACTER(LEN=2) :: element_symbol
158 : CHARACTER(LEN=default_string_length) :: bname
159 : INTEGER :: handle, i, iatom, ikind, isgf, ispin, &
160 : nao, natom, nimages, nkind, no, norb, &
161 : nref, ns, nsgf, nspin, nvec, nx, order
162 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf, nsgfat
163 : INTEGER, DIMENSION(2) :: nocc
164 : LOGICAL :: cubes_iao, cubes_ibo, molden_iao, &
165 : molden_ibo, uniform_occupation
166 : REAL(KIND=dp) :: fin, fout, t1, t2, trace, zval
167 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fdiag
168 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: mcharge
169 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: moments
170 2 : REAL(KIND=dp), DIMENSION(:), POINTER :: occupation_numbers
171 2 : TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: smat_kind
172 2 : TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:, :) :: oce_atom
173 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
174 : TYPE(cp_fm_type) :: p_orb_ref, p_ref_orb, smo
175 2 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: c_loc_orb, ciao, cloc, cvec, iao_coef
176 2 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij_atom
177 : TYPE(cp_fm_type), POINTER :: mo_coeff
178 2 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: sref, sro, sxo
179 2 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
180 : TYPE(dbcsr_type) :: dmat
181 : TYPE(dbcsr_type), POINTER :: smat
182 : TYPE(dft_control_type), POINTER :: dft_control
183 2 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: oce_basis_set_list, orb_basis_set_list, &
184 2 : ref_basis_set_list
185 : TYPE(gto_basis_set_type), POINTER :: ocebasis, orbbasis, refbasis
186 2 : TYPE(mo_set_type), ALLOCATABLE, DIMENSION(:) :: mo_iao, mo_loc
187 2 : TYPE(mo_set_type), DIMENSION(:), POINTER :: my_mos
188 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
189 2 : POINTER :: sro_list, srr_list, sxo_list
190 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
191 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
192 : TYPE(qs_kind_type), POINTER :: qs_kind
193 : TYPE(qs_ks_env_type), POINTER :: ks_env
194 : TYPE(section_vals_type), POINTER :: iao_cubes_section, iao_molden_section, &
195 : ibo_cc_section, ibo_cubes_section, &
196 : ibo_molden_section
197 :
198 : ! only do IAO analysis if explicitly requested
199 0 : CPASSERT(iao_env%do_iao)
200 :
201 : ! k-points?
202 2 : CALL get_qs_env(qs_env, dft_control=dft_control)
203 2 : nspin = dft_control%nspins
204 2 : nimages = dft_control%nimages
205 2 : IF (nimages > 1) THEN
206 0 : IF (unit_nr > 0) THEN
207 : WRITE (UNIT=unit_nr, FMT="(T2,A)") &
208 0 : "K-Points: Intrinsic Atomic Orbitals Analysis not available."
209 : END IF
210 : END IF
211 0 : IF (nimages > 1) RETURN
212 :
213 2 : CALL timeset(routineN, handle)
214 :
215 2 : IF (unit_nr > 0) THEN
216 1 : WRITE (unit_nr, '(/,T2,A)') '!-----------------------------------------------------------------------------!'
217 1 : WRITE (UNIT=unit_nr, FMT="(T24,A)") "INTRINSIC ATOMIC ORBITALS ANALYSIS"
218 1 : WRITE (UNIT=unit_nr, FMT="(T13,A)") "G. Knizia, J. Chem. Theory Comput. 9, 4834-4843 (2013)"
219 1 : WRITE (unit_nr, '(T2,A)') '!-----------------------------------------------------------------------------!'
220 : END IF
221 2 : CALL cite_reference(Knizia2013)
222 :
223 : ! input sections
224 2 : iao_molden_section => iao_env%iao_molden_section
225 2 : iao_cubes_section => iao_env%iao_cubes_section
226 2 : ibo_molden_section => iao_env%ibo_molden_section
227 2 : ibo_cubes_section => iao_env%ibo_cubes_section
228 2 : ibo_cc_section => iao_env%ibo_cc_section
229 :
230 : !
231 2 : molden_iao = .FALSE.
232 2 : IF (ASSOCIATED(iao_molden_section)) THEN
233 0 : CALL section_vals_get(iao_molden_section, explicit=molden_iao)
234 : END IF
235 2 : cubes_iao = .FALSE.
236 2 : IF (ASSOCIATED(iao_cubes_section)) THEN
237 0 : CALL section_vals_get(iao_cubes_section, explicit=cubes_iao)
238 : END IF
239 2 : molden_ibo = .FALSE.
240 2 : IF (ASSOCIATED(ibo_molden_section)) THEN
241 0 : CALL section_vals_get(ibo_molden_section, explicit=molden_ibo)
242 : END IF
243 2 : cubes_ibo = .FALSE.
244 2 : IF (ASSOCIATED(ibo_cubes_section)) THEN
245 0 : CALL section_vals_get(ibo_cubes_section, explicit=cubes_ibo)
246 : END IF
247 :
248 : ! check for or generate reference basis
249 2 : CALL create_minbas_set(qs_env, unit_nr)
250 :
251 : ! overlap matrices
252 2 : CALL get_qs_env(qs_env, matrix_s_kp=matrix_s)
253 2 : smat => matrix_s(1, 1)%matrix
254 : !
255 2 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
256 2 : nkind = SIZE(qs_kind_set)
257 16 : ALLOCATE (ref_basis_set_list(nkind), orb_basis_set_list(nkind))
258 6 : DO ikind = 1, nkind
259 4 : qs_kind => qs_kind_set(ikind)
260 4 : NULLIFY (ref_basis_set_list(ikind)%gto_basis_set)
261 4 : NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
262 4 : NULLIFY (refbasis, orbbasis)
263 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type="ORB")
264 4 : IF (ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
265 4 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=refbasis, basis_type="MIN")
266 6 : IF (ASSOCIATED(refbasis)) ref_basis_set_list(ikind)%gto_basis_set => refbasis
267 : END DO
268 :
269 : ! neighbor lists
270 2 : NULLIFY (srr_list, sro_list)
271 2 : CALL setup_neighbor_list(srr_list, ref_basis_set_list, qs_env=qs_env)
272 2 : CALL setup_neighbor_list(sro_list, ref_basis_set_list, orb_basis_set_list, qs_env=qs_env)
273 :
274 : ! Srr and Sro overlap matrices
275 2 : NULLIFY (sref, sro)
276 2 : CALL get_qs_env(qs_env, ks_env=ks_env)
277 : CALL build_overlap_matrix_simple(ks_env, sref, &
278 2 : ref_basis_set_list, ref_basis_set_list, srr_list)
279 : CALL build_overlap_matrix_simple(ks_env, sro, &
280 2 : ref_basis_set_list, orb_basis_set_list, sro_list)
281 : !
282 2 : IF (PRESENT(mos)) THEN
283 2 : my_mos => mos
284 : ELSE
285 0 : CALL get_qs_env(qs_env, mos=my_mos)
286 : END IF
287 2 : CALL get_mo_set(my_mos(1), mo_coeff=mo_coeff)
288 2 : CALL dbcsr_get_info(sro(1)%matrix, nfullrows_total=nref, nfullcols_total=nao)
289 :
290 : ! Projectors
291 2 : NULLIFY (fm_struct)
292 : CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nref, &
293 2 : ncol_global=nao, para_env=mo_coeff%matrix_struct%para_env)
294 2 : CALL cp_fm_create(p_ref_orb, fm_struct)
295 2 : CALL cp_fm_struct_release(fm_struct)
296 2 : NULLIFY (fm_struct)
297 : CALL cp_fm_struct_create(fm_struct, context=mo_coeff%matrix_struct%context, nrow_global=nao, &
298 2 : ncol_global=nref, para_env=mo_coeff%matrix_struct%para_env)
299 2 : CALL cp_fm_create(p_orb_ref, fm_struct)
300 2 : CALL cp_fm_struct_release(fm_struct)
301 : !
302 2 : CALL iao_projectors(smat, sref(1)%matrix, sro(1)%matrix, p_orb_ref, p_ref_orb, iao_env%eps_svd)
303 :
304 : ! occupied orbitals, orbitals considered for IAO generation
305 2 : nocc(1:2) = 0
306 4 : DO ispin = 1, nspin
307 : CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nmo=norb, uniform_occupation=uniform_occupation, &
308 2 : occupation_numbers=occupation_numbers)
309 2 : IF (uniform_occupation) THEN
310 2 : nvec = norb
311 : ELSE
312 0 : nvec = 0
313 0 : DO i = 1, norb
314 0 : IF (occupation_numbers(i) > iao_env%eps_occ) THEN
315 0 : nvec = i
316 : ELSE
317 : EXIT
318 : END IF
319 : END DO
320 : END IF
321 6 : nocc(ispin) = nvec
322 : END DO
323 : ! generate IAOs
324 12 : ALLOCATE (iao_coef(nspin), cvec(nspin))
325 4 : DO ispin = 1, nspin
326 2 : CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff)
327 2 : nvec = nocc(ispin)
328 4 : IF (nvec > 0) THEN
329 2 : NULLIFY (fm_struct)
330 2 : CALL cp_fm_struct_create(fm_struct, ncol_global=nvec, template_fmstruct=mo_coeff%matrix_struct)
331 2 : CALL cp_fm_create(cvec(ispin), fm_struct)
332 2 : CALL cp_fm_struct_release(fm_struct)
333 2 : CALL cp_fm_to_fm(mo_coeff, cvec(ispin), nvec)
334 : !
335 2 : NULLIFY (fm_struct)
336 2 : CALL cp_fm_struct_create(fm_struct, ncol_global=nref, template_fmstruct=mo_coeff%matrix_struct)
337 2 : CALL cp_fm_create(iao_coef(ispin), fm_struct)
338 2 : CALL cp_fm_struct_release(fm_struct)
339 : !
340 2 : CALL intrinsic_ao_calc(smat, p_orb_ref, p_ref_orb, cvec(ispin), iao_coef(ispin))
341 : END IF
342 : END DO
343 : !
344 2 : IF (iao_env%do_charges) THEN
345 : ! MOs in IAO basis
346 6 : ALLOCATE (ciao(nspin))
347 4 : DO ispin = 1, nspin
348 2 : CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
349 2 : CALL cp_fm_create(smo, mo_coeff%matrix_struct)
350 2 : NULLIFY (fm_struct)
351 2 : CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
352 2 : CALL cp_fm_create(ciao(ispin), fm_struct)
353 2 : CALL cp_fm_struct_release(fm_struct)
354 : ! CIAO = A(T)*S*C
355 2 : CALL cp_dbcsr_sm_fm_multiply(smat, mo_coeff, smo, ncol=norb)
356 2 : CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
357 6 : CALL cp_fm_release(smo)
358 : END DO
359 : !
360 : ! population analysis
361 2 : IF (unit_nr > 0) THEN
362 1 : WRITE (unit_nr, '(/,T2,A)') 'Intrinsic AO Population Analysis '
363 : END IF
364 2 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
365 2 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
366 8 : ALLOCATE (mcharge(natom, nspin))
367 2 : CALL dbcsr_get_info(sref(1)%matrix, nfullrows_total=nref)
368 4 : DO ispin = 1, nspin
369 : CALL get_mo_set(my_mos(ispin), &
370 : uniform_occupation=uniform_occupation, &
371 2 : occupation_numbers=occupation_numbers)
372 : ! diagonal block matrix
373 2 : CALL dbcsr_create(dmat, template=sref(1)%matrix)
374 2 : CALL dbcsr_reserve_diag_blocks(dmat)
375 : !
376 2 : CALL iao_calculate_dmat(ciao(ispin), dmat, occupation_numbers, uniform_occupation)
377 2 : CALL dbcsr_trace(dmat, trace)
378 2 : IF (unit_nr > 0) THEN
379 1 : WRITE (unit_nr, '(T2,A,I2,T66,F15.4)') 'Number of Electrons: Trace(Piao) Spin ', ispin, trace
380 : END IF
381 2 : CALL iao_charges(dmat, mcharge(:, ispin))
382 8 : CALL dbcsr_release(dmat)
383 : END DO
384 : CALL print_atomic_charges(particle_set, qs_kind_set, unit_nr, "Intrinsic Atomic Orbital Charges", &
385 2 : electronic_charges=mcharge)
386 2 : DEALLOCATE (mcharge)
387 4 : CALL cp_fm_release(ciao)
388 : END IF
389 : !
390 : ! Deallocate the neighbor list structure
391 2 : CALL release_neighbor_list_sets(srr_list)
392 2 : CALL release_neighbor_list_sets(sro_list)
393 2 : IF (ASSOCIATED(sref)) CALL dbcsr_deallocate_matrix_set(sref)
394 2 : IF (ASSOCIATED(sro)) CALL dbcsr_deallocate_matrix_set(sro)
395 2 : CALL cp_fm_release(p_ref_orb)
396 2 : CALL cp_fm_release(p_orb_ref)
397 2 : CALL cp_fm_release(cvec)
398 2 : DEALLOCATE (orb_basis_set_list)
399 : !
400 2 : IF (iao_env%do_oce) THEN
401 : ! generate OCE basis
402 0 : IF (unit_nr > 0) THEN
403 0 : WRITE (unit_nr, '(T2,A)') "IAO One-Center Expansion: OCE Basis Set"
404 : END IF
405 0 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
406 0 : nkind = SIZE(qs_kind_set)
407 0 : ALLOCATE (oce_basis_set_list(nkind), orb_basis_set_list(nkind))
408 0 : DO ikind = 1, nkind
409 0 : qs_kind => qs_kind_set(ikind)
410 0 : NULLIFY (orb_basis_set_list(ikind)%gto_basis_set)
411 0 : NULLIFY (orbbasis)
412 0 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orbbasis, basis_type="ORB")
413 0 : IF (ASSOCIATED(orbbasis)) orb_basis_set_list(ikind)%gto_basis_set => orbbasis
414 : END DO
415 0 : DO ikind = 1, nkind
416 0 : orbbasis => orb_basis_set_list(ikind)%gto_basis_set
417 0 : CPASSERT(ASSOCIATED(orbbasis))
418 0 : NULLIFY (ocebasis)
419 0 : CALL create_oce_basis(ocebasis, orbbasis, iao_env%lmax_oce, iao_env%nbas_oce)
420 0 : CALL init_orb_basis_set(ocebasis)
421 0 : CALL init_interaction_radii_orb_basis(ocebasis, dft_control%qs_control%eps_pgf_orb)
422 0 : oce_basis_set_list(ikind)%gto_basis_set => ocebasis
423 0 : IF (unit_nr > 0) THEN
424 0 : qs_kind => qs_kind_set(ikind)
425 0 : CALL get_qs_kind(qs_kind, zeff=zval, element_symbol=element_symbol)
426 0 : CALL get_gto_basis_set(ocebasis, name=bname, nsgf=nsgf)
427 0 : WRITE (unit_nr, '(T2,A,A,T14,A,I4,T40,A,A30)') "Kind: ", element_symbol, "NBasFun: ", nsgf, &
428 0 : "OCE Basis: ", ADJUSTL(TRIM(bname))
429 : END IF
430 : END DO
431 0 : IF (unit_nr > 0) WRITE (unit_nr, *)
432 : ! OCE basis overlap
433 0 : ALLOCATE (smat_kind(nkind))
434 0 : DO ikind = 1, nkind
435 0 : ocebasis => oce_basis_set_list(ikind)%gto_basis_set
436 0 : CALL get_gto_basis_set(gto_basis_set=ocebasis, nsgf=nsgf)
437 0 : ALLOCATE (smat_kind(ikind)%array(nsgf, nsgf))
438 : END DO
439 0 : CALL oce_overlap_matrix(smat_kind, oce_basis_set_list)
440 : ! overlap between IAO OCE basis and orbital basis
441 0 : NULLIFY (sxo, sxo_list)
442 0 : CALL setup_neighbor_list(sxo_list, oce_basis_set_list, orb_basis_set_list, qs_env=qs_env)
443 0 : CALL get_qs_env(qs_env, ks_env=ks_env)
444 : CALL build_overlap_matrix_simple(ks_env, sxo, &
445 0 : oce_basis_set_list, orb_basis_set_list, sxo_list)
446 : !
447 : ! One Center Expansion of IAOs
448 0 : CALL get_qs_env(qs_env=qs_env, natom=natom)
449 0 : ALLOCATE (oce_atom(natom, nspin))
450 0 : CALL iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
451 : !
452 0 : DO ikind = 1, nkind
453 0 : ocebasis => oce_basis_set_list(ikind)%gto_basis_set
454 0 : CALL deallocate_gto_basis_set(ocebasis)
455 : END DO
456 0 : DEALLOCATE (oce_basis_set_list, orb_basis_set_list)
457 : !
458 0 : CALL release_neighbor_list_sets(sxo_list)
459 0 : IF (ASSOCIATED(sxo)) CALL dbcsr_deallocate_matrix_set(sxo)
460 0 : DO ikind = 1, nkind
461 0 : DEALLOCATE (smat_kind(ikind)%array)
462 : END DO
463 0 : DEALLOCATE (smat_kind)
464 0 : DO ispin = 1, nspin
465 0 : DO iatom = 1, natom
466 0 : DEALLOCATE (oce_atom(iatom, ispin)%array)
467 : END DO
468 : END DO
469 0 : DEALLOCATE (oce_atom)
470 : END IF
471 : !
472 2 : IF (molden_iao) THEN
473 : ! Print basis functions: molden file
474 0 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
475 0 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
476 0 : IF (unit_nr > 0) THEN
477 0 : WRITE (unit_nr, '(T2,A)') ' Write IAO in MOLDEN Format'
478 : END IF
479 0 : ALLOCATE (mo_iao(nspin))
480 0 : DO ispin = 1, nspin
481 0 : CALL get_mo_set(my_mos(ispin), nao=nao)
482 0 : CALL allocate_mo_set(mo_iao(ispin), nao, nref, nref, 0.0_dp, 1.0_dp, 0.0_dp)
483 0 : CALL init_mo_set(mo_iao(ispin), fm_ref=iao_coef(ispin), name="iao_set")
484 0 : CALL cp_fm_to_fm(iao_coef(ispin), mo_iao(ispin)%mo_coeff)
485 0 : CALL set_mo_set(mo_iao(ispin), homo=nref)
486 0 : mo_iao(ispin)%occupation_numbers = 1.0_dp
487 : END DO
488 : ! Print IAO basis functions: molden format
489 0 : CALL write_mos_molden(mo_iao, qs_kind_set, particle_set, iao_molden_section)
490 0 : DO ispin = 1, nspin
491 0 : CALL deallocate_mo_set(mo_iao(ispin))
492 : END DO
493 0 : DEALLOCATE (mo_iao)
494 : END IF
495 2 : IF (cubes_iao) THEN
496 0 : IF (unit_nr > 0) THEN
497 0 : WRITE (unit_nr, '(T2,A)') ' Write IAO as CUBE Files'
498 : END IF
499 : ! Print basis functions: cube file
500 0 : CALL print_iao_cubes(qs_env, iao_cubes_section, iao_coef, ref_basis_set_list)
501 : END IF
502 : !
503 : ! Intrinsic bond orbitals
504 2 : IF (iao_env%do_bondorbitals) THEN
505 0 : IF (unit_nr > 0) THEN
506 0 : WRITE (unit_nr, '(/,T2,A)') 'Intrinsic Bond Orbital Generation'
507 : END IF
508 : ! MOs in IAO basis -> ciao
509 0 : ALLOCATE (ciao(nspin))
510 0 : DO ispin = 1, nspin
511 0 : CALL get_mo_set(my_mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=norb)
512 0 : CALL cp_fm_create(smo, mo_coeff%matrix_struct)
513 0 : NULLIFY (fm_struct)
514 0 : CALL cp_fm_struct_create(fm_struct, nrow_global=nref, template_fmstruct=mo_coeff%matrix_struct)
515 0 : CALL cp_fm_create(ciao(ispin), fm_struct)
516 0 : CALL cp_fm_struct_release(fm_struct)
517 : ! CIAO = A(T)*S*C
518 0 : CALL cp_dbcsr_sm_fm_multiply(smat, mo_coeff, smo, ncol=norb)
519 0 : CALL parallel_gemm('T', 'N', nref, norb, nao, 1.0_dp, iao_coef(ispin), smo, 0.0_dp, ciao(ispin))
520 0 : CALL cp_fm_release(smo)
521 : END DO
522 : !
523 : ! localize occupied orbitals nocc(ispin), see IAO generation
524 : !
525 0 : ALLOCATE (cloc(nspin), c_loc_orb(nspin))
526 0 : DO ispin = 1, nspin
527 0 : NULLIFY (fm_struct)
528 : CALL cp_fm_struct_create(fm_struct, ncol_global=nocc(ispin), &
529 0 : template_fmstruct=ciao(ispin)%matrix_struct)
530 0 : CALL cp_fm_create(cloc(ispin), fm_struct)
531 0 : CALL cp_fm_struct_release(fm_struct)
532 0 : CALL cp_fm_to_fm(ciao(ispin), cloc(ispin), ncol=nocc(ispin))
533 : END DO
534 0 : CALL cp_fm_release(ciao)
535 0 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, natom=natom)
536 0 : ALLOCATE (first_sgf(natom), last_sgf(natom), nsgfat(natom))
537 0 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
538 : CALL get_particle_set(particle_set, qs_kind_set, first_sgf=first_sgf, last_sgf=last_sgf, &
539 0 : nsgf=nsgfat, basis=ref_basis_set_list)
540 :
541 0 : IF (iao_env%loc_operator /= do_iaoloc_pm2) THEN
542 0 : CPABORT("IAO localization operator NYA")
543 : END IF
544 0 : IF (iao_env%eloc_function /= do_iaoloc_enone) THEN
545 0 : CPABORT("IAO energy weight localization NYA")
546 : END IF
547 :
548 0 : DO ispin = 1, nspin
549 0 : nvec = nocc(ispin)
550 0 : IF (nvec > 0) THEN
551 0 : ALLOCATE (zij_atom(natom, 1), fdiag(nvec))
552 0 : NULLIFY (fm_struct)
553 : CALL cp_fm_struct_create(fm_struct, nrow_global=nvec, ncol_global=nvec, &
554 0 : template_fmstruct=cloc(ispin)%matrix_struct)
555 0 : DO iatom = 1, natom
556 0 : CALL cp_fm_create(zij_atom(iatom, 1), fm_struct)
557 0 : isgf = first_sgf(iatom)
558 : CALL parallel_gemm('T', 'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
559 : 0.0_dp, zij_atom(iatom, 1), &
560 0 : a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
561 : END DO
562 0 : CALL cp_fm_struct_release(fm_struct)
563 : !
564 0 : t1 = m_walltime()
565 0 : order = 4
566 0 : fin = 0.0_dp
567 0 : DO iatom = 1, natom
568 0 : CALL cp_fm_get_diag(zij_atom(iatom, 1), fdiag)
569 0 : fin = fin + SUM(fdiag**order)
570 : END DO
571 0 : fin = fin**(1._dp/order)
572 : ! QR localization
573 0 : CALL scdm_qrfact(cloc(ispin))
574 : !
575 0 : DO iatom = 1, natom
576 0 : isgf = first_sgf(iatom)
577 : CALL parallel_gemm('T', 'N', nvec, nvec, nsgfat(iatom), 1.0_dp, cloc(ispin), cloc(ispin), &
578 : 0.0_dp, zij_atom(iatom, 1), &
579 0 : a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
580 : END DO
581 0 : fout = 0.0_dp
582 0 : DO iatom = 1, natom
583 0 : CALL cp_fm_get_diag(zij_atom(iatom, 1), fdiag)
584 0 : fout = fout + SUM(fdiag**order)
585 : END DO
586 0 : fout = fout**(1._dp/order)
587 0 : DEALLOCATE (fdiag)
588 0 : t2 = m_walltime()
589 0 : IF (unit_nr > 0) THEN
590 : WRITE (unit_nr, '(T2,A,F14.8,A,F14.8,A,F8.2)') &
591 0 : 'SCDM pre-localization: fin=', fin, " fout=", fout, " Time=", t2 - t1
592 : END IF
593 : !
594 0 : CALL pm_localization(zij_atom, cloc(ispin), order, 1.E-12_dp, unit_nr)
595 : !
596 0 : CALL cp_fm_release(zij_atom)
597 0 : CALL get_mo_set(my_mos(ispin), nao=nao)
598 0 : NULLIFY (fm_struct)
599 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, &
600 0 : template_fmstruct=cloc(ispin)%matrix_struct)
601 0 : CALL cp_fm_create(c_loc_orb(ispin), fm_struct)
602 0 : CALL cp_fm_struct_release(fm_struct)
603 : CALL parallel_gemm('N', 'N', nao, nvec, nref, 1.0_dp, iao_coef(ispin), cloc(ispin), &
604 0 : 0.0_dp, c_loc_orb(ispin))
605 : END IF
606 : END DO
607 0 : DEALLOCATE (first_sgf, last_sgf, nsgfat)
608 0 : CALL cp_fm_release(cloc)
609 : !
610 0 : IF (iao_env%do_center) THEN
611 0 : IF (unit_nr > 0) THEN
612 0 : WRITE (unit_nr, '(T2,A)') ' Calculate Localized Orbital Centers and Spread'
613 0 : IF (iao_env%pos_periodic) THEN
614 0 : WRITE (unit_nr, '(T2,A)') ' Use Berry Phase Position Operator'
615 : ELSE
616 0 : WRITE (unit_nr, '(T2,A)') ' Use Local Position Operator'
617 : END IF
618 : END IF
619 0 : nvec = MAXVAL(nocc)
620 : ! x y z m2(1) m2(2)
621 0 : ALLOCATE (moments(5, nvec, nspin))
622 0 : moments = 0.0_dp
623 0 : IF (iao_env%pos_periodic) THEN
624 0 : CALL center_spread_berry(qs_env, c_loc_orb, moments)
625 : ELSE
626 0 : CALL center_spread_loc(qs_env, c_loc_orb, moments)
627 : END IF
628 : !
629 0 : IF (ASSOCIATED(ibo_cc_section)) THEN
630 0 : CALL print_center_spread(moments, nocc, ibo_cc_section)
631 : END IF
632 : !
633 0 : IF (PRESENT(bond_centers)) THEN
634 0 : nx = SIZE(bond_centers, 1)
635 0 : no = SIZE(bond_centers, 2)
636 0 : ns = SIZE(bond_centers, 3)
637 0 : CPASSERT(no >= nvec)
638 0 : CPASSERT(ns == nspin)
639 0 : CPASSERT(nx >= 3)
640 0 : nx = MIN(nx, 5)
641 0 : bond_centers(1:nx, 1:nvec, 1:ns) = moments(1:nx, 1:nvec, 1:ns)
642 : END IF
643 0 : DEALLOCATE (moments)
644 : END IF
645 : !
646 0 : IF (molden_ibo) THEN
647 0 : CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set)
648 0 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set)
649 0 : IF (unit_nr > 0) THEN
650 0 : WRITE (unit_nr, '(T2,A)') ' Write IBO in MOLDEN Format'
651 : END IF
652 0 : ALLOCATE (mo_loc(nspin))
653 0 : DO ispin = 1, nspin
654 0 : CALL get_mo_set(my_mos(ispin), nao=nao)
655 0 : nvec = nocc(ispin)
656 0 : CALL allocate_mo_set(mo_loc(ispin), nao, nvec, nvec, 0.0_dp, 1.0_dp, 0.0_dp)
657 0 : CALL init_mo_set(mo_loc(ispin), fm_ref=c_loc_orb(ispin), name="ibo_orb")
658 0 : CALL cp_fm_to_fm(c_loc_orb(ispin), mo_loc(ispin)%mo_coeff)
659 0 : CALL set_mo_set(mo_loc(ispin), homo=nvec)
660 0 : mo_loc(ispin)%occupation_numbers = 1.0_dp
661 : END DO
662 : ! Print IBO functions: molden format
663 0 : CALL write_mos_molden(mo_loc, qs_kind_set, particle_set, ibo_molden_section)
664 0 : DO ispin = 1, nspin
665 0 : CALL deallocate_mo_set(mo_loc(ispin))
666 : END DO
667 0 : DEALLOCATE (mo_loc)
668 : END IF
669 : !
670 0 : IF (cubes_ibo) THEN
671 0 : IF (unit_nr > 0) THEN
672 0 : WRITE (unit_nr, '(T2,A)') ' Write IBO on CUBE files'
673 : END IF
674 : ! Print localized orbital functions: cube file
675 0 : CALL print_ibo_cubes(qs_env, ibo_cubes_section, c_loc_orb)
676 : END IF
677 : !
678 0 : IF (PRESENT(mos) .AND. ALLOCATED(c_loc_orb)) THEN
679 : ! c_loc_orb -> mos
680 0 : DO ispin = 1, nspin
681 0 : CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
682 0 : nvec = nocc(ispin)
683 0 : CALL cp_fm_to_fm(c_loc_orb(ispin), mo_coeff, ncol=nvec)
684 : END DO
685 : END IF
686 0 : CALL cp_fm_release(c_loc_orb)
687 : END IF
688 2 : DEALLOCATE (ref_basis_set_list)
689 :
690 2 : IF (PRESENT(c_iao_coef)) THEN
691 2 : CALL cp_fm_release(c_iao_coef)
692 6 : ALLOCATE (c_iao_coef(nspin))
693 4 : DO ispin = 1, nspin
694 2 : CALL cp_fm_create(c_iao_coef(ispin), iao_coef(ispin)%matrix_struct)
695 4 : CALL cp_fm_to_fm(iao_coef(ispin), c_iao_coef(ispin))
696 : END DO
697 : END IF
698 2 : CALL cp_fm_release(iao_coef)
699 :
700 2 : IF (unit_nr > 0) THEN
701 : WRITE (unit_nr, '(T2,A)') &
702 1 : '!----------------------------END OF IAO ANALYSIS------------------------------!'
703 : END IF
704 :
705 2 : CALL timestop(handle)
706 :
707 12 : END SUBROUTINE iao_wfn_analysis
708 :
709 : ! **************************************************************************************************
710 : !> \brief Computes projector matrices for ref to orb basis and reverse
711 : !> \param smat ...
712 : !> \param sref ...
713 : !> \param s_r_o ...
714 : !> \param p_o_r ...
715 : !> \param p_r_o ...
716 : !> \param eps_svd ...
717 : ! **************************************************************************************************
718 14 : SUBROUTINE iao_projectors(smat, sref, s_r_o, p_o_r, p_r_o, eps_svd)
719 : TYPE(dbcsr_type), INTENT(INOUT) :: smat, sref, s_r_o
720 : TYPE(cp_fm_type), INTENT(INOUT) :: p_o_r, p_r_o
721 : REAL(KIND=dp), INTENT(IN) :: eps_svd
722 :
723 : CHARACTER(len=*), PARAMETER :: routineN = 'iao_projectors'
724 :
725 : INTEGER :: handle, norb, nref
726 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
727 : TYPE(cp_fm_type) :: fm_inv, fm_s, fm_sro
728 :
729 2 : CALL timeset(routineN, handle)
730 :
731 2 : CALL cp_fm_get_info(p_r_o, nrow_global=nref, ncol_global=norb)
732 2 : CALL cp_fm_create(fm_sro, p_r_o%matrix_struct)
733 2 : CALL copy_dbcsr_to_fm(s_r_o, fm_sro)
734 :
735 2 : NULLIFY (fm_struct)
736 : CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
737 2 : template_fmstruct=p_r_o%matrix_struct)
738 2 : CALL cp_fm_create(fm_s, fm_struct, name="smat")
739 2 : CALL cp_fm_create(fm_inv, fm_struct, name="sinv")
740 2 : CALL cp_fm_struct_release(fm_struct)
741 2 : CALL copy_dbcsr_to_fm(smat, fm_s)
742 2 : CALL cp_fm_invert(fm_s, fm_inv, eps_svd=eps_svd)
743 2 : CALL parallel_gemm('N', 'T', norb, nref, norb, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_o_r)
744 2 : CALL cp_fm_release(fm_s)
745 2 : CALL cp_fm_release(fm_inv)
746 :
747 2 : NULLIFY (fm_struct)
748 : CALL cp_fm_struct_create(fm_struct, nrow_global=nref, ncol_global=nref, &
749 2 : template_fmstruct=p_r_o%matrix_struct)
750 2 : CALL cp_fm_create(fm_s, fm_struct, name="sref")
751 2 : CALL cp_fm_create(fm_inv, fm_struct, name="sinv")
752 2 : CALL cp_fm_struct_release(fm_struct)
753 2 : CALL copy_dbcsr_to_fm(sref, fm_s)
754 2 : CALL cp_fm_invert(fm_s, fm_inv, eps_svd=eps_svd)
755 2 : CALL parallel_gemm('N', 'N', nref, norb, nref, 1.0_dp, fm_inv, fm_sro, 0.0_dp, p_r_o)
756 2 : CALL cp_fm_release(fm_s)
757 2 : CALL cp_fm_release(fm_inv)
758 :
759 2 : CALL cp_fm_release(fm_sro)
760 :
761 2 : CALL timestop(handle)
762 :
763 2 : END SUBROUTINE iao_projectors
764 :
765 : ! **************************************************************************************************
766 : !> \brief Computes intrinsic orbitals for a given MO vector set
767 : !> \param smat ...
768 : !> \param p_o_r ...
769 : !> \param p_r_o ...
770 : !> \param cvec ...
771 : !> \param avec ...
772 : ! **************************************************************************************************
773 18 : SUBROUTINE intrinsic_ao_calc(smat, p_o_r, p_r_o, cvec, avec)
774 : TYPE(dbcsr_type), INTENT(INOUT) :: smat
775 : TYPE(cp_fm_type), INTENT(INOUT) :: p_o_r, p_r_o, cvec, avec
776 :
777 : CHARACTER(len=*), PARAMETER :: routineN = 'intrinsic_ao_calc'
778 :
779 : INTEGER :: handle, nao, norb, nref
780 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
781 : TYPE(cp_fm_type) :: ctvec, pc, sc, sct, vec1
782 :
783 2 : CALL timeset(routineN, handle)
784 :
785 : ! number of orbitals
786 2 : CALL cp_fm_get_info(cvec, ncol_global=norb)
787 : ! basis set sizes
788 2 : CALL cp_fm_get_info(p_r_o, nrow_global=nref, ncol_global=nao)
789 : ! temp array
790 2 : NULLIFY (fm_struct)
791 : CALL cp_fm_struct_create(fm_struct, nrow_global=nref, ncol_global=norb, &
792 2 : template_fmstruct=cvec%matrix_struct)
793 2 : CALL cp_fm_create(pc, fm_struct)
794 2 : CALL cp_fm_struct_release(fm_struct)
795 : ! CT = orth(Por.Pro.C)
796 2 : CALL parallel_gemm('N', 'N', nref, norb, nao, 1.0_dp, p_r_o, cvec, 0.0_dp, pc)
797 2 : CALL cp_fm_create(ctvec, cvec%matrix_struct)
798 2 : CALL parallel_gemm('N', 'N', nao, norb, nref, 1.0_dp, p_o_r, pc, 0.0_dp, ctvec)
799 2 : CALL cp_fm_release(pc)
800 2 : CALL make_basis_lowdin(ctvec, norb, smat)
801 : ! S*C and S*CT
802 2 : CALL cp_fm_create(sc, cvec%matrix_struct)
803 2 : CALL cp_dbcsr_sm_fm_multiply(smat, cvec, sc, ncol=norb)
804 2 : CALL cp_fm_create(sct, cvec%matrix_struct)
805 2 : CALL cp_dbcsr_sm_fm_multiply(smat, ctvec, sct, ncol=norb)
806 : CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=nref, &
807 2 : template_fmstruct=cvec%matrix_struct)
808 2 : CALL cp_fm_create(pc, fm_struct)
809 2 : CALL cp_fm_struct_release(fm_struct)
810 : ! V1 = (CT*SCT(T))Por
811 2 : CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sct, p_o_r, 0.0_dp, pc)
812 2 : NULLIFY (fm_struct)
813 : CALL cp_fm_struct_create(fm_struct, nrow_global=nao, ncol_global=nref, &
814 2 : template_fmstruct=cvec%matrix_struct)
815 2 : CALL cp_fm_create(vec1, fm_struct)
816 2 : CALL cp_fm_struct_release(fm_struct)
817 2 : CALL parallel_gemm('N', 'N', nao, nref, norb, 1.0_dp, ctvec, pc, 0.0_dp, vec1)
818 : ! A = C*SC(T)*V1
819 2 : CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
820 2 : CALL parallel_gemm('N', 'N', nao, nref, norb, 1.0_dp, cvec, pc, 0.0_dp, avec)
821 : ! V1 = Por - V1
822 2 : CALL cp_fm_geadd(1.0_dp, 'N', p_o_r, -1.0_dp, vec1)
823 : ! A = A + V1
824 2 : CALL cp_fm_geadd(1.0_dp, 'N', vec1, 1.0_dp, avec)
825 : ! A = A - C*SC(T)*V1
826 2 : CALL parallel_gemm('T', 'N', norb, nref, nao, 1.0_dp, sc, vec1, 0.0_dp, pc)
827 2 : CALL parallel_gemm('N', 'N', nao, nref, norb, -1.0_dp, cvec, pc, 1.0_dp, avec)
828 : ! A = orth(A)
829 2 : CALL make_basis_lowdin(avec, nref, smat)
830 :
831 : ! clean up
832 2 : CALL cp_fm_release(pc)
833 2 : CALL cp_fm_release(vec1)
834 2 : CALL cp_fm_release(sc)
835 2 : CALL cp_fm_release(sct)
836 2 : CALL cp_fm_release(ctvec)
837 :
838 2 : CALL timestop(handle)
839 :
840 2 : END SUBROUTINE intrinsic_ao_calc
841 :
842 : ! **************************************************************************************************
843 : !> \brief Calculate the density matrix from fm vectors using occupation numbers
844 : !> \param cvec ...
845 : !> \param density_matrix ...
846 : !> \param occupation ...
847 : !> \param uniform_occupation ...
848 : ! **************************************************************************************************
849 8 : SUBROUTINE iao_calculate_dmat_diag(cvec, density_matrix, occupation, uniform_occupation)
850 :
851 : TYPE(cp_fm_type), INTENT(IN) :: cvec
852 : TYPE(dbcsr_type) :: density_matrix
853 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: occupation
854 : LOGICAL, INTENT(IN) :: uniform_occupation
855 :
856 : CHARACTER(len=*), PARAMETER :: routineN = 'iao_calculate_dmat_diag'
857 :
858 : INTEGER :: handle, ncol
859 : REAL(KIND=dp) :: alpha
860 : TYPE(cp_fm_type) :: fm_tmp
861 :
862 8 : CALL timeset(routineN, handle)
863 :
864 8 : CALL dbcsr_set(density_matrix, 0.0_dp)
865 :
866 8 : CALL cp_fm_get_info(cvec, ncol_global=ncol)
867 8 : IF (.NOT. uniform_occupation) THEN
868 6 : CALL cp_fm_create(fm_tmp, cvec%matrix_struct)
869 6 : CALL cp_fm_to_fm(cvec, fm_tmp)
870 6 : CALL cp_fm_column_scale(fm_tmp, occupation(1:ncol))
871 6 : alpha = 1.0_dp
872 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
873 : matrix_v=cvec, matrix_g=fm_tmp, &
874 6 : ncol=ncol, alpha=alpha)
875 6 : CALL cp_fm_release(fm_tmp)
876 : ELSE
877 2 : alpha = occupation(1)
878 : CALL cp_dbcsr_plus_fm_fm_t(sparse_matrix=density_matrix, &
879 2 : matrix_v=cvec, ncol=ncol, alpha=alpha)
880 : END IF
881 :
882 8 : CALL timestop(handle)
883 :
884 8 : END SUBROUTINE iao_calculate_dmat_diag
885 :
886 : ! **************************************************************************************************
887 : !> \brief Calculate the density matrix from fm vectors using an occupation matrix
888 : !> \param cvec ...
889 : !> \param density_matrix ...
890 : !> \param weight ...
891 : !> \param occmat ...
892 : ! **************************************************************************************************
893 0 : SUBROUTINE iao_calculate_dmat_full(cvec, density_matrix, weight, occmat)
894 :
895 : TYPE(cp_fm_type), INTENT(IN) :: cvec
896 : TYPE(dbcsr_type) :: density_matrix
897 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: weight
898 : TYPE(cp_fm_type), INTENT(IN) :: occmat
899 :
900 : CHARACTER(len=*), PARAMETER :: routineN = 'iao_calculate_dmat_full'
901 :
902 : INTEGER :: handle, ic, jc, ncol
903 : REAL(KIND=dp) :: alpha
904 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
905 : TYPE(cp_fm_type) :: fm1, fm2
906 :
907 0 : CALL timeset(routineN, handle)
908 :
909 0 : CALL dbcsr_set(density_matrix, 0.0_dp)
910 :
911 : CALL cp_fm_struct_create(fm_struct, ncol_global=1, &
912 0 : template_fmstruct=cvec%matrix_struct)
913 0 : CALL cp_fm_create(fm1, fm_struct)
914 0 : CALL cp_fm_create(fm2, fm_struct)
915 0 : CALL cp_fm_struct_release(fm_struct)
916 :
917 0 : CALL cp_fm_get_info(cvec, ncol_global=ncol)
918 0 : DO ic = 1, ncol
919 0 : CALL cp_fm_to_fm(cvec, fm1, ncol=1, source_start=ic, target_start=1)
920 0 : DO jc = 1, ncol
921 0 : CALL cp_fm_to_fm(cvec, fm2, ncol=1, source_start=jc, target_start=1)
922 0 : CALL cp_fm_get_element(occmat, ic, jc, alpha)
923 0 : alpha = weight(ic)*alpha
924 : CALL cp_dbcsr_plus_fm_fm_t(density_matrix, fm1, fm2, ncol=1, &
925 0 : alpha=alpha, symmetry_mode=1)
926 : END DO
927 : END DO
928 0 : CALL cp_fm_release(fm1)
929 0 : CALL cp_fm_release(fm2)
930 :
931 0 : CALL timestop(handle)
932 :
933 0 : END SUBROUTINE iao_calculate_dmat_full
934 :
935 : ! **************************************************************************************************
936 : !> \brief compute the atomic charges (orthogonal basis)
937 : !> \param p_matrix ...
938 : !> \param charges ...
939 : ! **************************************************************************************************
940 10 : SUBROUTINE iao_charges(p_matrix, charges)
941 : TYPE(dbcsr_type) :: p_matrix
942 : REAL(KIND=dp), DIMENSION(:) :: charges
943 :
944 : INTEGER :: i, iblock_col, iblock_row
945 : REAL(kind=dp) :: trace
946 10 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block
947 : TYPE(dbcsr_iterator_type) :: iter
948 : TYPE(mp_comm_type) :: group
949 :
950 40 : charges = 0.0_dp
951 :
952 10 : CALL dbcsr_iterator_start(iter, p_matrix)
953 25 : DO WHILE (dbcsr_iterator_blocks_left(iter))
954 15 : NULLIFY (p_block)
955 15 : CALL dbcsr_iterator_next_block(iter, iblock_row, iblock_col, p_block)
956 15 : CPASSERT(iblock_row == iblock_col)
957 15 : trace = 0.0_dp
958 45 : DO i = 1, SIZE(p_block, 1)
959 45 : trace = trace + p_block(i, i)
960 : END DO
961 15 : charges(iblock_row) = trace
962 : END DO
963 10 : CALL dbcsr_iterator_stop(iter)
964 :
965 10 : CALL dbcsr_get_info(p_matrix, group=group)
966 70 : CALL group%sum(charges)
967 :
968 10 : END SUBROUTINE iao_charges
969 :
970 : ! **************************************************************************************************
971 : !> \brief Computes and prints the Cube Files for the intrinsic atomic orbitals
972 : !> \param qs_env ...
973 : !> \param print_section ...
974 : !> \param iao_coef ...
975 : !> \param basis_set_list ...
976 : ! **************************************************************************************************
977 0 : SUBROUTINE print_iao_cubes(qs_env, print_section, iao_coef, basis_set_list)
978 : TYPE(qs_environment_type), POINTER :: qs_env
979 : TYPE(section_vals_type), POINTER :: print_section
980 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: iao_coef
981 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
982 :
983 : CHARACTER(LEN=default_path_length) :: filename, title
984 : INTEGER :: i, i_rep, ispin, ivec, iw, j, n_rep, &
985 : nat, natom, norb, nspins, nstart
986 0 : INTEGER, DIMENSION(:), POINTER :: atom_list, blk_sizes, first_bas, stride
987 : LOGICAL :: mpi_io
988 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
989 : TYPE(cell_type), POINTER :: cell
990 : TYPE(cp_logger_type), POINTER :: logger
991 : TYPE(dft_control_type), POINTER :: dft_control
992 : TYPE(particle_list_type), POINTER :: particles
993 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
994 : TYPE(pw_c1d_gs_type) :: wf_g
995 : TYPE(pw_env_type), POINTER :: pw_env
996 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
997 : TYPE(pw_r3d_rs_type) :: wf_r
998 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
999 : TYPE(qs_subsys_type), POINTER :: subsys
1000 :
1001 0 : logger => cp_get_default_logger()
1002 0 : stride => section_get_ivals(print_section, "STRIDE")
1003 0 : CALL section_vals_val_get(print_section, "ATOM_LIST", n_rep_val=n_rep)
1004 :
1005 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1006 0 : subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
1007 0 : CALL qs_subsys_get(subsys, particles=particles)
1008 :
1009 0 : CALL get_qs_env(qs_env=qs_env, natom=natom)
1010 0 : ALLOCATE (blk_sizes(natom), first_bas(0:natom))
1011 0 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=blk_sizes, basis=basis_set_list)
1012 0 : first_bas(0) = 0
1013 0 : DO i = 1, natom
1014 0 : first_bas(i) = first_bas(i - 1) + blk_sizes(i)
1015 : END DO
1016 :
1017 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1018 0 : CALL auxbas_pw_pool%create_pw(wf_r)
1019 0 : CALL auxbas_pw_pool%create_pw(wf_g)
1020 :
1021 0 : nspins = SIZE(iao_coef)
1022 0 : nstart = MIN(1, n_rep)
1023 :
1024 0 : DO ispin = 1, nspins
1025 0 : DO i_rep = nstart, n_rep
1026 0 : CALL cp_fm_get_info(iao_coef(ispin), ncol_global=norb)
1027 0 : IF (i_rep == 0) THEN
1028 0 : nat = natom
1029 : ELSE
1030 0 : CALL section_vals_val_get(print_section, "ATOM_LIST", i_rep_val=i_rep, i_vals=atom_list)
1031 0 : nat = SIZE(atom_list)
1032 : END IF
1033 0 : DO i = 1, nat
1034 0 : IF (i_rep == 0) THEN
1035 0 : j = i
1036 : ELSE
1037 0 : j = atom_list(i)
1038 : END IF
1039 0 : CPASSERT(j >= 1 .AND. j <= natom)
1040 0 : DO ivec = first_bas(j - 1) + 1, first_bas(j)
1041 0 : WRITE (filename, '(a4,I5.5,a1,I1.1)') "IAO_", ivec, "_", ispin
1042 0 : WRITE (title, *) "Intrinsic Atomic Orbitals ", ivec, " atom ", j, " spin ", ispin
1043 0 : mpi_io = .TRUE.
1044 : iw = cp_print_key_unit_nr(logger, print_section, "", extension=".cube", &
1045 : middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE., &
1046 0 : mpi_io=mpi_io)
1047 : CALL calculate_wavefunction(iao_coef(ispin), ivec, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1048 0 : cell, dft_control, particle_set, pw_env)
1049 0 : CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
1050 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "", mpi_io=mpi_io)
1051 : END DO
1052 : END DO
1053 : END DO
1054 : END DO
1055 0 : DEALLOCATE (blk_sizes, first_bas)
1056 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1057 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1058 :
1059 0 : END SUBROUTINE print_iao_cubes
1060 :
1061 : ! **************************************************************************************************
1062 : !> \brief Computes and prints the Cube Files for the intrinsic bond orbitals
1063 : !> \param qs_env ...
1064 : !> \param print_section ...
1065 : !> \param ibo_coef ...
1066 : ! **************************************************************************************************
1067 0 : SUBROUTINE print_ibo_cubes(qs_env, print_section, ibo_coef)
1068 : TYPE(qs_environment_type), POINTER :: qs_env
1069 : TYPE(section_vals_type), POINTER :: print_section
1070 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: ibo_coef
1071 :
1072 : CHARACTER(LEN=default_path_length) :: filename, title
1073 : INTEGER :: i, i_rep, ispin, iw, j, n_rep, norb, &
1074 : nspins, nstart, nstate
1075 0 : INTEGER, DIMENSION(:), POINTER :: state_list, stride
1076 : LOGICAL :: mpi_io
1077 0 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1078 : TYPE(cell_type), POINTER :: cell
1079 : TYPE(cp_logger_type), POINTER :: logger
1080 : TYPE(dft_control_type), POINTER :: dft_control
1081 : TYPE(particle_list_type), POINTER :: particles
1082 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1083 : TYPE(pw_c1d_gs_type) :: wf_g
1084 : TYPE(pw_env_type), POINTER :: pw_env
1085 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
1086 : TYPE(pw_r3d_rs_type) :: wf_r
1087 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1088 : TYPE(qs_subsys_type), POINTER :: subsys
1089 :
1090 0 : CPASSERT(ASSOCIATED(print_section))
1091 :
1092 0 : logger => cp_get_default_logger()
1093 0 : stride => section_get_ivals(print_section, "STRIDE")
1094 0 : CALL section_vals_val_get(print_section, "STATE_LIST", n_rep_val=n_rep)
1095 :
1096 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, &
1097 0 : subsys=subsys, cell=cell, particle_set=particle_set, pw_env=pw_env, dft_control=dft_control)
1098 0 : CALL qs_subsys_get(subsys, particles=particles)
1099 :
1100 0 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
1101 0 : CALL auxbas_pw_pool%create_pw(wf_r)
1102 0 : CALL auxbas_pw_pool%create_pw(wf_g)
1103 :
1104 0 : nspins = SIZE(ibo_coef)
1105 0 : nstart = MIN(1, n_rep)
1106 :
1107 0 : DO ispin = 1, nspins
1108 0 : DO i_rep = nstart, n_rep
1109 0 : CALL cp_fm_get_info(ibo_coef(ispin), ncol_global=norb)
1110 0 : IF (i_rep == 0) THEN
1111 0 : nstate = norb
1112 : ELSE
1113 0 : CALL section_vals_val_get(print_section, "STATE_LIST", i_rep_val=i_rep, i_vals=state_list)
1114 0 : nstate = SIZE(state_list)
1115 : END IF
1116 0 : DO i = 1, nstate
1117 0 : IF (i_rep == 0) THEN
1118 0 : j = i
1119 : ELSE
1120 0 : j = state_list(i)
1121 : END IF
1122 0 : CPASSERT(j >= 1 .AND. j <= norb)
1123 0 : WRITE (filename, '(a4,I5.5,a1,I1.1)') "IBO_", j, "_", ispin
1124 0 : WRITE (title, *) "Intrinsic Atomic Orbitals ", j, " spin ", ispin
1125 0 : mpi_io = .TRUE.
1126 : iw = cp_print_key_unit_nr(logger, print_section, "", extension=".cube", &
1127 : middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE., &
1128 0 : mpi_io=mpi_io)
1129 : CALL calculate_wavefunction(ibo_coef(ispin), j, wf_r, wf_g, atomic_kind_set, qs_kind_set, &
1130 0 : cell, dft_control, particle_set, pw_env)
1131 0 : CALL cp_pw_to_cube(wf_r, iw, title, particles=particles, stride=stride, mpi_io=mpi_io)
1132 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "", mpi_io=mpi_io)
1133 : END DO
1134 : END DO
1135 : END DO
1136 :
1137 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
1138 0 : CALL auxbas_pw_pool%give_back_pw(wf_g)
1139 :
1140 0 : END SUBROUTINE print_ibo_cubes
1141 :
1142 : ! **************************************************************************************************
1143 : !> \brief prints charge center and spreads for all orbitals
1144 : !> \param moments ...
1145 : !> \param nocc ...
1146 : !> \param print_section ...
1147 : !> \param iounit ...
1148 : ! **************************************************************************************************
1149 0 : SUBROUTINE print_center_spread(moments, nocc, print_section, iounit)
1150 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: moments
1151 : INTEGER, DIMENSION(:), INTENT(IN) :: nocc
1152 : TYPE(section_vals_type), OPTIONAL, POINTER :: print_section
1153 : INTEGER, INTENT(IN), OPTIONAL :: iounit
1154 :
1155 : CHARACTER(LEN=default_path_length) :: filename
1156 : INTEGER :: is, ispin, iw, nspin
1157 : TYPE(cp_logger_type), POINTER :: logger
1158 :
1159 0 : logger => cp_get_default_logger()
1160 0 : nspin = SIZE(moments, 3)
1161 :
1162 0 : IF (PRESENT(print_section)) THEN
1163 0 : WRITE (filename, '(A18,I1.1)') "IBO_CENTERS_SPREAD"
1164 : iw = cp_print_key_unit_nr(logger, print_section, "", extension=".csp", &
1165 0 : middle_name=TRIM(filename), file_position="REWIND", log_filename=.FALSE.)
1166 0 : ELSEIF (PRESENT(iounit)) THEN
1167 0 : iw = iounit
1168 : ELSE
1169 0 : iw = -1
1170 : END IF
1171 0 : IF (iw > 0) THEN
1172 0 : DO ispin = 1, nspin
1173 0 : WRITE (iw, "(T2,A,i1)") "Intrinsic Bond Orbitals: Centers and Spread for Spin ", ispin
1174 0 : WRITE (iw, "(A7,T30,A6,T68,A7)") "State", "Center", "Spreads"
1175 0 : DO is = 1, nocc(ispin)
1176 0 : WRITE (iw, "(i7,3F15.8,8X,2F10.5)") is, moments(:, is, ispin)
1177 : END DO
1178 : END DO
1179 : END IF
1180 0 : IF (PRESENT(print_section)) THEN
1181 0 : CALL cp_print_key_finished_output(iw, logger, print_section, "")
1182 : END IF
1183 :
1184 0 : END SUBROUTINE print_center_spread
1185 :
1186 : ! **************************************************************************************************
1187 : !> \brief ...
1188 : !> \param qs_env ...
1189 : !> \param c_loc_orb ...
1190 : !> \param moments ...
1191 : ! **************************************************************************************************
1192 0 : SUBROUTINE center_spread_loc(qs_env, c_loc_orb, moments)
1193 : TYPE(qs_environment_type), POINTER :: qs_env
1194 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c_loc_orb
1195 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: moments
1196 :
1197 : CHARACTER(len=*), PARAMETER :: routineN = 'center_spread_loc'
1198 : INTEGER, DIMENSION(6), PARAMETER :: list = (/1, 2, 3, 4, 7, 9/)
1199 :
1200 : INTEGER :: handle, i, iop, iorb, ispin, nao, norb, &
1201 : nspin
1202 : REAL(KIND=dp) :: xmii
1203 : REAL(KIND=dp), DIMENSION(3) :: rpoint
1204 : TYPE(cell_type), POINTER :: cell
1205 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1206 : TYPE(cp_fm_type) :: ccmat, ocvec
1207 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: dipmat
1208 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_s
1209 : TYPE(dbcsr_type), POINTER :: omat, smat
1210 :
1211 0 : CALL timeset(routineN, handle)
1212 :
1213 0 : CALL get_qs_env(qs_env=qs_env, matrix_s_kp=matrix_s)
1214 0 : smat => matrix_s(1, 1)%matrix
1215 0 : rpoint = 0.0_dp
1216 0 : nspin = SIZE(c_loc_orb)
1217 0 : moments = 0.0_dp
1218 :
1219 0 : ALLOCATE (dipmat(9))
1220 0 : DO i = 1, 9
1221 0 : ALLOCATE (dipmat(i)%matrix)
1222 0 : CALL dbcsr_copy(dipmat(i)%matrix, smat)
1223 0 : CALL dbcsr_set(dipmat(i)%matrix, 0.0_dp)
1224 : END DO
1225 :
1226 0 : CALL build_local_moment_matrix(qs_env, dipmat, 2, rpoint)
1227 :
1228 0 : DO ispin = 1, nspin
1229 0 : CALL cp_fm_get_info(c_loc_orb(ispin), nrow_global=nao, ncol_global=norb)
1230 0 : CALL cp_fm_create(ocvec, c_loc_orb(ispin)%matrix_struct)
1231 0 : NULLIFY (fm_struct)
1232 : CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
1233 0 : template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1234 0 : CALL cp_fm_create(ccmat, fm_struct)
1235 0 : CALL cp_fm_struct_release(fm_struct)
1236 0 : DO i = 1, 6
1237 0 : iop = list(i)
1238 0 : omat => dipmat(iop)%matrix
1239 0 : CALL cp_dbcsr_sm_fm_multiply(omat, c_loc_orb(ispin), ocvec, ncol=norb)
1240 : CALL parallel_gemm('T', 'N', norb, norb, nao, 1.0_dp, c_loc_orb(ispin), ocvec, &
1241 0 : 0.0_dp, ccmat)
1242 0 : DO iorb = 1, norb
1243 0 : CALL cp_fm_get_element(ccmat, iorb, iorb, xmii)
1244 0 : IF (iop <= 3) THEN
1245 0 : moments(iop, iorb, ispin) = moments(iop, iorb, ispin) + xmii
1246 0 : moments(4, iorb, ispin) = moments(4, iorb, ispin) - xmii**2
1247 : ELSE
1248 0 : moments(4, iorb, ispin) = moments(4, iorb, ispin) + xmii
1249 : END IF
1250 : END DO
1251 : END DO
1252 0 : CALL cp_fm_release(ocvec)
1253 0 : CALL cp_fm_release(ccmat)
1254 : END DO
1255 :
1256 0 : DO i = 1, 9
1257 0 : CALL dbcsr_release(dipmat(i)%matrix)
1258 0 : DEALLOCATE (dipmat(i)%matrix)
1259 : END DO
1260 0 : DEALLOCATE (dipmat)
1261 :
1262 0 : CALL get_qs_env(qs_env=qs_env, cell=cell)
1263 0 : DO ispin = 1, nspin
1264 0 : DO iorb = 1, norb
1265 0 : moments(1:3, iorb, ispin) = pbc(moments(1:3, iorb, ispin), cell)
1266 : END DO
1267 : END DO
1268 :
1269 0 : CALL timestop(handle)
1270 :
1271 0 : END SUBROUTINE center_spread_loc
1272 :
1273 : ! **************************************************************************************************
1274 : !> \brief ...
1275 : !> \param qs_env ...
1276 : !> \param c_loc_orb ...
1277 : !> \param moments ...
1278 : ! **************************************************************************************************
1279 0 : SUBROUTINE center_spread_berry(qs_env, c_loc_orb, moments)
1280 : TYPE(qs_environment_type), POINTER :: qs_env
1281 : TYPE(cp_fm_type), DIMENSION(:), INTENT(IN) :: c_loc_orb
1282 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(INOUT) :: moments
1283 :
1284 : CHARACTER(len=*), PARAMETER :: routineN = 'center_spread_berry'
1285 :
1286 : COMPLEX(KIND=dp) :: z
1287 : INTEGER :: dim_op, handle, i, idir, ispin, istate, &
1288 : j, jdir, nao, norb, nspin
1289 : REAL(dp), DIMENSION(3) :: c, cpbc
1290 : REAL(dp), DIMENSION(6) :: weights
1291 : REAL(KIND=dp) :: imagpart, realpart, spread_i, spread_ii
1292 : TYPE(cell_type), POINTER :: cell
1293 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1294 : TYPE(cp_fm_type) :: opvec
1295 0 : TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:, :) :: zij
1296 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
1297 0 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
1298 :
1299 0 : CALL timeset(routineN, handle)
1300 :
1301 0 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, cell=cell)
1302 :
1303 0 : IF (cell%orthorhombic) THEN
1304 0 : dim_op = 3
1305 : ELSE
1306 0 : dim_op = 6
1307 : END IF
1308 0 : ALLOCATE (op_sm_set(2, dim_op))
1309 0 : DO i = 1, dim_op
1310 0 : DO j = 1, 2
1311 0 : NULLIFY (op_sm_set(j, i)%matrix)
1312 0 : ALLOCATE (op_sm_set(j, i)%matrix)
1313 0 : CALL dbcsr_copy(op_sm_set(j, i)%matrix, matrix_s(1)%matrix)
1314 0 : CALL dbcsr_set(op_sm_set(j, i)%matrix, 0.0_dp)
1315 : END DO
1316 : END DO
1317 0 : CALL initialize_weights(cell, weights)
1318 :
1319 0 : CALL compute_berry_operator(qs_env, cell, op_sm_set, dim_op)
1320 :
1321 0 : nspin = SIZE(c_loc_orb, 1)
1322 0 : DO ispin = 1, nspin
1323 0 : CALL cp_fm_get_info(c_loc_orb(ispin), nrow_global=nao, ncol_global=norb)
1324 0 : CALL cp_fm_create(opvec, c_loc_orb(ispin)%matrix_struct)
1325 0 : NULLIFY (fm_struct)
1326 : CALL cp_fm_struct_create(fm_struct, nrow_global=norb, ncol_global=norb, &
1327 0 : template_fmstruct=c_loc_orb(ispin)%matrix_struct)
1328 0 : ALLOCATE (zij(2, dim_op))
1329 :
1330 : ! Compute zij here
1331 0 : DO i = 1, dim_op
1332 0 : DO j = 1, 2
1333 0 : CALL cp_fm_create(zij(j, i), fm_struct)
1334 0 : CALL cp_fm_set_all(zij(j, i), 0.0_dp)
1335 0 : CALL cp_dbcsr_sm_fm_multiply(op_sm_set(j, i)%matrix, c_loc_orb(ispin), opvec, ncol=norb)
1336 0 : CALL parallel_gemm("T", "N", norb, norb, nao, 1.0_dp, c_loc_orb(ispin), opvec, 0.0_dp, zij(j, i))
1337 : END DO
1338 : END DO
1339 :
1340 0 : CALL cp_fm_struct_release(fm_struct)
1341 0 : CALL cp_fm_release(opvec)
1342 :
1343 0 : DO istate = 1, norb
1344 0 : c = 0.0_dp
1345 0 : spread_i = 0.0_dp
1346 0 : spread_ii = 0.0_dp
1347 0 : DO jdir = 1, dim_op
1348 0 : CALL cp_fm_get_element(zij(1, jdir), istate, istate, realpart)
1349 0 : CALL cp_fm_get_element(zij(2, jdir), istate, istate, imagpart)
1350 0 : z = CMPLX(realpart, imagpart, dp)
1351 : spread_i = spread_i - weights(jdir)* &
1352 0 : LOG(realpart*realpart + imagpart*imagpart)/twopi/twopi
1353 : spread_ii = spread_ii + weights(jdir)* &
1354 0 : (1.0_dp - (realpart*realpart + imagpart*imagpart))/twopi/twopi
1355 0 : IF (jdir < 4) THEN
1356 0 : DO idir = 1, 3
1357 : c(idir) = c(idir) + &
1358 0 : (cell%hmat(idir, jdir)/twopi)*AIMAG(LOG(z))
1359 : END DO
1360 : END IF
1361 : END DO
1362 0 : cpbc = pbc(c, cell)
1363 0 : moments(1:3, istate, ispin) = cpbc(1:3)
1364 0 : moments(4, istate, ispin) = spread_i
1365 0 : moments(5, istate, ispin) = spread_ii
1366 : END DO
1367 :
1368 0 : CALL cp_fm_release(zij)
1369 :
1370 : END DO
1371 :
1372 0 : DO i = 1, dim_op
1373 0 : DO j = 1, 2
1374 0 : CALL dbcsr_release(op_sm_set(j, i)%matrix)
1375 0 : DEALLOCATE (op_sm_set(j, i)%matrix)
1376 : END DO
1377 : END DO
1378 0 : DEALLOCATE (op_sm_set)
1379 :
1380 0 : CALL timestop(handle)
1381 :
1382 0 : END SUBROUTINE center_spread_berry
1383 :
1384 : ! **************************************************************************************************
1385 : !> \brief ...
1386 : !> \param qs_env ...
1387 : !> \param oce_basis_set_list ...
1388 : !> \param smat_kind ...
1389 : !> \param sxo ...
1390 : !> \param iao_coef ...
1391 : !> \param oce_atom ...
1392 : !> \param unit_nr ...
1393 : ! **************************************************************************************************
1394 0 : SUBROUTINE iao_oce_expansion(qs_env, oce_basis_set_list, smat_kind, sxo, iao_coef, oce_atom, unit_nr)
1395 : TYPE(qs_environment_type), POINTER :: qs_env
1396 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: oce_basis_set_list
1397 : TYPE(cp_2d_r_p_type), DIMENSION(:) :: smat_kind
1398 : TYPE(dbcsr_p_type), DIMENSION(:) :: sxo
1399 : TYPE(cp_fm_type), DIMENSION(:) :: iao_coef
1400 : TYPE(cp_2d_r_p_type), DIMENSION(:, :) :: oce_atom
1401 : INTEGER, INTENT(IN) :: unit_nr
1402 :
1403 : CHARACTER(len=*), PARAMETER :: routineN = 'iao_oce_expansion'
1404 :
1405 : INTEGER :: handle, i, i1, i2, iao, iatom, ikind, &
1406 : iset, ishell, ispin, l, m, maxl, n, &
1407 : natom, nkind, noce, ns, nset, nsgf, &
1408 : nspin
1409 0 : INTEGER, DIMENSION(:), POINTER :: iao_blk_sizes, nshell, oce_blk_sizes, &
1410 0 : orb_blk_sizes
1411 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgf, last_sgf, lval
1412 : LOGICAL :: found
1413 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: ev, vector
1414 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: amat, bmat, filter, oce_comp, prol, vec
1415 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ablock
1416 0 : TYPE(cp_2d_r_p_type), ALLOCATABLE, DIMENSION(:) :: sinv_kind
1417 : TYPE(dbcsr_distribution_type) :: dbcsr_dist
1418 : TYPE(dbcsr_type) :: iao_vec, sx_vec
1419 : TYPE(gto_basis_set_type), POINTER :: oce_basis
1420 : TYPE(mp_comm_type) :: group
1421 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1422 0 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1423 : TYPE(qs_ks_env_type), POINTER :: ks_env
1424 :
1425 0 : CALL timeset(routineN, handle)
1426 :
1427 : ! basis sets: block sizes
1428 : CALL dbcsr_get_info(sxo(1)%matrix, row_blk_size=oce_blk_sizes, col_blk_size=orb_blk_sizes, &
1429 0 : distribution=dbcsr_dist)
1430 0 : CALL get_qs_env(qs_env=qs_env, ks_env=ks_env, natom=natom)
1431 0 : CALL get_ks_env(ks_env=ks_env, particle_set=particle_set, qs_kind_set=qs_kind_set)
1432 0 : ALLOCATE (iao_blk_sizes(natom))
1433 0 : DO iatom = 1, natom
1434 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1435 0 : CALL get_qs_kind(qs_kind_set(ikind), nsgf=ns, basis_type="MIN")
1436 0 : iao_blk_sizes(iatom) = ns
1437 : END DO
1438 :
1439 : CALL dbcsr_create(iao_vec, name="IAO_VEC", dist=dbcsr_dist, &
1440 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=orb_blk_sizes, &
1441 0 : col_blk_size=iao_blk_sizes)
1442 : CALL dbcsr_create(sx_vec, name="SX_VEC", dist=dbcsr_dist, &
1443 : matrix_type=dbcsr_type_no_symmetry, row_blk_size=oce_blk_sizes, &
1444 0 : col_blk_size=iao_blk_sizes)
1445 0 : CALL dbcsr_reserve_diag_blocks(matrix=sx_vec)
1446 :
1447 0 : nkind = SIZE(smat_kind)
1448 0 : ALLOCATE (sinv_kind(nkind))
1449 0 : DO ikind = 1, nkind
1450 0 : noce = SIZE(smat_kind(ikind)%array, 1)
1451 0 : ALLOCATE (sinv_kind(ikind)%array(noce, noce))
1452 0 : sinv_kind(ikind)%array = smat_kind(ikind)%array
1453 0 : CALL invmat_symm(sinv_kind(ikind)%array)
1454 : END DO
1455 0 : CALL dbcsr_get_info(iao_vec, group=group)
1456 :
1457 0 : nspin = SIZE(iao_coef, 1)
1458 0 : ALLOCATE (oce_comp(natom, nspin))
1459 0 : oce_comp = 0.0_dp
1460 0 : DO ispin = 1, nspin
1461 0 : CALL copy_fm_to_dbcsr(iao_coef(ispin), iao_vec)
1462 : CALL dbcsr_multiply("N", "N", 1.0_dp, sxo(1)%matrix, iao_vec, 0.0_dp, sx_vec, &
1463 0 : retain_sparsity=.TRUE.)
1464 0 : DO iatom = 1, natom
1465 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1466 : CALL dbcsr_get_block_p(matrix=sx_vec, &
1467 0 : row=iatom, col=iatom, BLOCK=ablock, found=found)
1468 0 : IF (found) THEN
1469 0 : n = SIZE(ablock, 1)
1470 0 : m = SIZE(ablock, 2)
1471 0 : ablock = MATMUL(sinv_kind(ikind)%array, ablock)
1472 0 : ALLOCATE (amat(n, m), bmat(m, m), ev(m), vec(m, m))
1473 0 : amat(1:n, 1:m) = MATMUL(smat_kind(ikind)%array(1:n, 1:n), ablock(1:n, 1:m))
1474 0 : bmat(1:m, 1:m) = MATMUL(TRANSPOSE(ablock(1:n, 1:m)), amat(1:n, 1:m))
1475 0 : CALL jacobi(bmat, ev, vec)
1476 0 : oce_comp(iatom, ispin) = SUM(ev(1:m))/REAL(m, KIND=dp)
1477 0 : DO i = 1, m
1478 0 : ev(i) = 1._dp/SQRT(ev(i))
1479 0 : bmat(1:m, i) = vec(1:m, i)*ev(i)
1480 : END DO
1481 0 : bmat(1:m, 1:m) = MATMUL(bmat(1:m, 1:m), TRANSPOSE(vec(1:m, 1:m)))
1482 0 : ablock(1:n, 1:m) = MATMUL(ablock(1:n, 1:m), bmat(1:m, 1:m))
1483 0 : DEALLOCATE (amat, bmat, ev, vec)
1484 : END IF
1485 : END DO
1486 0 : CALL dbcsr_replicate_all(sx_vec)
1487 0 : DO iatom = 1, natom
1488 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1489 : CALL dbcsr_get_block_p(matrix=sx_vec, &
1490 0 : row=iatom, col=iatom, BLOCK=ablock, found=found)
1491 0 : CPASSERT(found)
1492 0 : n = SIZE(ablock, 1)
1493 0 : m = SIZE(ablock, 2)
1494 0 : ALLOCATE (oce_atom(iatom, ispin)%array(n, m))
1495 0 : oce_atom(iatom, ispin)%array = ablock
1496 : END DO
1497 : END DO
1498 :
1499 0 : CALL group%sum(oce_comp)
1500 0 : IF (unit_nr > 0) THEN
1501 0 : DO iatom = 1, natom
1502 0 : WRITE (unit_nr, "(T4,A,I6,T30,A,2F12.4)") "Atom:", iatom, " Completeness: ", oce_comp(iatom, 1:nspin)
1503 0 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, kind_number=ikind)
1504 0 : oce_basis => oce_basis_set_list(ikind)%gto_basis_set
1505 : CALL get_gto_basis_set(oce_basis, nset=nset, nshell=nshell, nsgf=nsgf, maxl=maxl, &
1506 0 : l=lval, first_sgf=first_sgf, last_sgf=last_sgf)
1507 0 : ALLOCATE (filter(nsgf, 0:maxl), vector(nsgf), prol(0:maxl, nspin))
1508 0 : filter = 0.0_dp
1509 0 : DO iset = 1, nset
1510 0 : DO ishell = 1, nshell(iset)
1511 0 : l = lval(ishell, iset)
1512 0 : i1 = first_sgf(ishell, iset)
1513 0 : i2 = last_sgf(ishell, iset)
1514 0 : filter(i1:i2, l) = 1.0_dp
1515 : END DO
1516 : END DO
1517 : !
1518 0 : n = SIZE(oce_atom(iatom, 1)%array, 1)
1519 0 : m = SIZE(oce_atom(iatom, 1)%array, 2)
1520 0 : CPASSERT(n == nsgf)
1521 0 : DO iao = 1, m
1522 0 : prol = 0.0_dp
1523 0 : DO ispin = 1, nspin
1524 0 : DO l = 0, maxl
1525 0 : vector(1:n) = oce_atom(iatom, ispin)%array(1:n, iao)*filter(1:n, l)
1526 0 : prol(l, ispin) = SUM(MATMUL(smat_kind(ikind)%array(1:n, 1:n), vector(1:n))*vector(1:n))
1527 : END DO
1528 : END DO
1529 0 : WRITE (unit_nr, "(T4,I3,T15,A,T39,(6F7.4))") iao, " l-contributions:", (SUM(prol(l, :)), l=0, maxl)
1530 : END DO
1531 0 : DEALLOCATE (filter, vector, prol)
1532 : END DO
1533 0 : WRITE (unit_nr, *)
1534 : END IF
1535 :
1536 0 : DEALLOCATE (oce_comp)
1537 0 : DO ikind = 1, nkind
1538 0 : DEALLOCATE (sinv_kind(ikind)%array)
1539 : END DO
1540 0 : DEALLOCATE (sinv_kind)
1541 0 : DEALLOCATE (iao_blk_sizes)
1542 0 : CALL dbcsr_release(iao_vec)
1543 0 : CALL dbcsr_release(sx_vec)
1544 :
1545 0 : CALL timestop(handle)
1546 :
1547 0 : END SUBROUTINE iao_oce_expansion
1548 :
1549 : ! **************************************************************************************************
1550 : !> \brief ...
1551 : !> \param smat_kind ...
1552 : !> \param basis_set_list ...
1553 : ! **************************************************************************************************
1554 0 : SUBROUTINE oce_overlap_matrix(smat_kind, basis_set_list)
1555 : TYPE(cp_2d_r_p_type), DIMENSION(:) :: smat_kind
1556 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
1557 :
1558 : CHARACTER(len=*), PARAMETER :: routineN = 'oce_overlap_matrix'
1559 :
1560 : INTEGER :: handle, ikind, iset, jset, ldsab, m1, &
1561 : m2, n1, n2, ncoa, ncob, nkind, nseta, &
1562 : sgfa, sgfb
1563 0 : INTEGER, DIMENSION(:), POINTER :: la_max, la_min, npgfa, nsgfa
1564 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
1565 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: oint, owork
1566 : REAL(KIND=dp), DIMENSION(3) :: rab
1567 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, scon_a, smat, zeta
1568 : TYPE(gto_basis_set_type), POINTER :: basis_set
1569 :
1570 0 : CALL timeset(routineN, handle)
1571 :
1572 0 : rab(1:3) = 0.0_dp
1573 :
1574 0 : nkind = SIZE(smat_kind)
1575 0 : ldsab = 0
1576 0 : DO ikind = 1, nkind
1577 0 : basis_set => basis_set_list(ikind)%gto_basis_set
1578 0 : CPASSERT(ASSOCIATED(basis_set))
1579 0 : CALL get_gto_basis_set(gto_basis_set=basis_set, maxco=m1, nsgf=m2)
1580 0 : ldsab = MAX(m1, m2, ldsab)
1581 : END DO
1582 :
1583 0 : ALLOCATE (oint(ldsab, ldsab), owork(ldsab, ldsab))
1584 :
1585 0 : DO ikind = 1, nkind
1586 0 : basis_set => basis_set_list(ikind)%gto_basis_set
1587 0 : CPASSERT(ASSOCIATED(basis_set))
1588 0 : smat => smat_kind(ikind)%array
1589 0 : smat = 0.0_dp
1590 : ! basis ikind
1591 0 : first_sgfa => basis_set%first_sgf
1592 0 : la_max => basis_set%lmax
1593 0 : la_min => basis_set%lmin
1594 0 : npgfa => basis_set%npgf
1595 0 : nseta = basis_set%nset
1596 0 : nsgfa => basis_set%nsgf_set
1597 0 : rpgfa => basis_set%pgf_radius
1598 0 : scon_a => basis_set%scon
1599 0 : zeta => basis_set%zet
1600 :
1601 0 : DO iset = 1, nseta
1602 :
1603 0 : ncoa = npgfa(iset)*ncoset(la_max(iset))
1604 0 : n1 = npgfa(iset)*(ncoset(la_max(iset)) - ncoset(la_min(iset) - 1))
1605 0 : sgfa = first_sgfa(1, iset)
1606 :
1607 0 : DO jset = 1, nseta
1608 :
1609 0 : ncob = npgfa(jset)*ncoset(la_max(jset))
1610 0 : n2 = npgfa(jset)*(ncoset(la_max(jset)) - ncoset(la_min(jset) - 1))
1611 0 : sgfb = first_sgfa(1, jset)
1612 :
1613 : ! calculate integrals and derivatives
1614 : CALL overlap_ab(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
1615 : la_max(jset), la_min(jset), npgfa(jset), rpgfa(:, jset), zeta(:, jset), &
1616 0 : rab, sab=oint(:, :))
1617 : ! Contraction
1618 : CALL contraction(oint(:, :), owork, ca=scon_a(:, sgfa:), na=n1, ma=nsgfa(iset), &
1619 0 : cb=scon_a(:, sgfb:), nb=n2, mb=nsgfa(jset), fscale=1.0_dp, trans=.FALSE.)
1620 : CALL block_add("IN", owork, nsgfa(iset), nsgfa(jset), smat, &
1621 0 : sgfa, sgfb, trans=.FALSE.)
1622 :
1623 : END DO
1624 : END DO
1625 :
1626 : END DO
1627 0 : DEALLOCATE (oint, owork)
1628 :
1629 0 : CALL timestop(handle)
1630 :
1631 0 : END SUBROUTINE oce_overlap_matrix
1632 :
1633 : ! **************************************************************************************************
1634 : !> \brief 2 by 2 rotation optimization for the Pipek-Mezey operator
1635 : !> \param zij_fm_set ...
1636 : !> \param vectors ...
1637 : !> \param order ...
1638 : !> \param accuracy ...
1639 : !> \param unit_nr ...
1640 : !> \par History
1641 : !> 05-2005 created
1642 : !> 10-2023 adapted from jacobi_rotation_pipek [JGH]
1643 : !> \author MI
1644 : ! **************************************************************************************************
1645 0 : SUBROUTINE pm_localization(zij_fm_set, vectors, order, accuracy, unit_nr)
1646 :
1647 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
1648 : TYPE(cp_fm_type), INTENT(IN) :: vectors
1649 : INTEGER, INTENT(IN) :: order
1650 : REAL(KIND=dp), INTENT(IN) :: accuracy
1651 : INTEGER, INTENT(IN) :: unit_nr
1652 :
1653 : INTEGER, PARAMETER :: max_sweeps = 250
1654 :
1655 : INTEGER :: iatom, istate, jstate, natom, nstate, &
1656 : sweeps
1657 : REAL(KIND=dp) :: aij, bij, ct, ftarget, mii, mij, mjj, &
1658 : st, t1, t2, theta, tolerance, tt
1659 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fdiag
1660 : TYPE(cp_fm_type) :: rmat
1661 :
1662 0 : CALL cp_fm_create(rmat, zij_fm_set(1, 1)%matrix_struct)
1663 0 : CALL cp_fm_set_all(rmat, 0.0_dp, 1.0_dp)
1664 :
1665 0 : CALL cp_fm_get_info(rmat, nrow_global=nstate)
1666 0 : ALLOCATE (fdiag(nstate))
1667 0 : tolerance = 1.0e10_dp
1668 0 : sweeps = 0
1669 0 : natom = SIZE(zij_fm_set, 1)
1670 0 : IF (unit_nr > 0) THEN
1671 : WRITE (unit_nr, '(A,T30,A,T45,A,T63,A,T77,A)') &
1672 0 : " Jacobi Localization ", "Sweep", "Functional", "Gradient", "Time"
1673 : END IF
1674 : ! do jacobi sweeps until converged
1675 0 : DO sweeps = 1, max_sweeps
1676 0 : t1 = m_walltime()
1677 0 : tolerance = 0.0_dp
1678 0 : DO istate = 1, nstate
1679 0 : DO jstate = istate + 1, nstate
1680 : aij = 0.0_dp
1681 : bij = 0.0_dp
1682 0 : DO iatom = 1, natom
1683 0 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, mii)
1684 0 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, jstate, mij)
1685 0 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), jstate, jstate, mjj)
1686 0 : IF (order == 2) THEN
1687 0 : aij = aij + 4._dp*mij**2 - (mii - mjj)**2
1688 0 : bij = bij + 4._dp*mij*(mii - mjj)
1689 0 : ELSEIF (order == 4) THEN
1690 : aij = aij - mii**4 - mjj**4 + 6._dp*(mii**2 + mjj**2)*mij**2 + &
1691 0 : mii**3*mjj + mii*mjj**3
1692 0 : bij = bij + 4._dp*mij*(mii**3 - mjj**3)
1693 : ELSE
1694 0 : CPABORT("PM order")
1695 : END IF
1696 : END DO
1697 0 : IF ((aij**2 + bij**2) < 1.E-10_dp) CYCLE
1698 0 : tolerance = tolerance + bij**2
1699 0 : theta = 0.25_dp*ATAN2(bij, -aij)
1700 0 : ct = COS(theta)
1701 0 : st = SIN(theta)
1702 :
1703 0 : CALL cp_fm_rot_cols(rmat, istate, jstate, ct, st)
1704 :
1705 0 : DO iatom = 1, natom
1706 0 : CALL cp_fm_rot_cols(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1707 0 : CALL cp_fm_rot_rows(zij_fm_set(iatom, 1), istate, jstate, ct, st)
1708 : END DO
1709 : END DO
1710 : END DO
1711 0 : ftarget = 0.0_dp
1712 0 : DO iatom = 1, natom
1713 0 : CALL cp_fm_get_diag(zij_fm_set(iatom, 1), fdiag)
1714 0 : ftarget = ftarget + SUM(fdiag**order)
1715 : END DO
1716 0 : ftarget = ftarget**(1._dp/order)
1717 0 : tolerance = SQRT(tolerance)
1718 0 : t2 = m_walltime()
1719 0 : tt = t2 - t1
1720 0 : IF (unit_nr > 0) THEN
1721 0 : WRITE (unit_nr, '(T31,I4,T39,F16.8,T55,G16.8,T71,F10.2)') sweeps, ftarget, tolerance, tt
1722 : END IF
1723 0 : IF (tolerance < accuracy) EXIT
1724 : END DO
1725 0 : DEALLOCATE (fdiag)
1726 :
1727 0 : CALL rotate_orbitals(rmat, vectors)
1728 0 : CALL cp_fm_release(rmat)
1729 :
1730 0 : END SUBROUTINE pm_localization
1731 : ! **************************************************************************************************
1732 :
1733 0 : END MODULE iao_analysis
|
