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