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