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 Driver for the localization that should be general
10 : !> for all the methods available and all the definition of the
11 : !> spread functional
12 : !> Write centers, spread and cubes only if required and for the
13 : !> selected states
14 : !> The localized functions are copied in the standard mos array
15 : !> for the next use
16 : !> \par History
17 : !> 01.2008 Teodoro Laino [tlaino] - University of Zurich
18 : !> - Merging the two localization codes and updating to new structures
19 : !> \author MI (04.2005)
20 : ! **************************************************************************************************
21 : MODULE qs_loc_methods
22 : USE atomic_kind_types, ONLY: atomic_kind_type,&
23 : deallocate_atomic_kind_set,&
24 : set_atomic_kind
25 : USE cell_types, ONLY: cell_type,&
26 : pbc
27 : USE cp_cfm_types, ONLY: cp_cfm_create,&
28 : cp_cfm_get_element,&
29 : cp_cfm_get_info,&
30 : cp_cfm_p_type,&
31 : cp_cfm_release,&
32 : cp_cfm_to_fm,&
33 : cp_cfm_type,&
34 : cp_fm_to_cfm
35 : USE cp_control_types, ONLY: dft_control_type
36 : USE cp_dbcsr_api, ONLY: dbcsr_copy,&
37 : dbcsr_deallocate_matrix,&
38 : dbcsr_p_type,&
39 : dbcsr_set
40 : USE cp_dbcsr_operations, ONLY: copy_dbcsr_to_fm,&
41 : cp_dbcsr_sm_fm_multiply
42 : USE cp_fm_basic_linalg, ONLY: cp_fm_rot_cols,&
43 : cp_fm_rot_rows,&
44 : cp_fm_schur_product
45 : USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
46 : fm_pool_create_fm
47 : USE cp_fm_struct, ONLY: cp_fm_struct_create,&
48 : cp_fm_struct_release,&
49 : cp_fm_struct_type
50 : USE cp_fm_types, ONLY: cp_fm_create,&
51 : cp_fm_get_element,&
52 : cp_fm_get_info,&
53 : cp_fm_maxabsval,&
54 : cp_fm_release,&
55 : cp_fm_set_all,&
56 : cp_fm_to_fm,&
57 : cp_fm_type
58 : USE cp_log_handling, ONLY: cp_get_default_logger,&
59 : cp_logger_type,&
60 : cp_to_string
61 : USE cp_output_handling, ONLY: cp_iter_string,&
62 : cp_p_file,&
63 : cp_print_key_finished_output,&
64 : cp_print_key_should_output,&
65 : cp_print_key_unit_nr
66 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
67 : USE cp_units, ONLY: cp_unit_from_cp2k
68 : USE input_constants, ONLY: &
69 : do_loc_crazy, do_loc_cs, do_loc_direct, do_loc_gapo, do_loc_jacobi, do_loc_l1_norm_sd, &
70 : do_loc_none, do_loc_scdm, dump_dcd, dump_dcd_aligned_cell, dump_xmol
71 : USE input_section_types, ONLY: section_get_ival,&
72 : section_get_ivals,&
73 : section_vals_get_subs_vals,&
74 : section_vals_type,&
75 : section_vals_val_get
76 : USE kinds, ONLY: default_path_length,&
77 : default_string_length,&
78 : dp,&
79 : sp
80 : USE machine, ONLY: m_flush
81 : USE mathconstants, ONLY: pi,&
82 : twopi
83 : USE message_passing, ONLY: mp_para_env_type
84 : USE molecule_types, ONLY: molecule_type
85 : USE motion_utils, ONLY: get_output_format
86 : USE orbital_pointers, ONLY: ncoset
87 : USE parallel_gemm_api, ONLY: parallel_gemm
88 : USE particle_list_types, ONLY: particle_list_type
89 : USE particle_methods, ONLY: get_particle_set,&
90 : write_particle_coordinates
91 : USE particle_types, ONLY: allocate_particle_set,&
92 : deallocate_particle_set,&
93 : particle_type
94 : USE physcon, ONLY: angstrom
95 : USE pw_env_types, ONLY: pw_env_get,&
96 : pw_env_type
97 : USE pw_pool_types, ONLY: pw_pool_type
98 : USE pw_types, ONLY: pw_c1d_gs_type,&
99 : pw_r3d_rs_type
100 : USE qs_collocate_density, ONLY: calculate_wavefunction
101 : USE qs_environment_types, ONLY: get_qs_env,&
102 : qs_environment_type
103 : USE qs_kind_types, ONLY: qs_kind_type
104 : USE qs_loc_types, ONLY: get_qs_loc_env,&
105 : qs_loc_env_type
106 : USE qs_localization_methods, ONLY: &
107 : approx_l1_norm_sd, cardoso_souloumiac, cardoso_souloumiac_pipek, crazy_rotations, &
108 : direct_mini, jacobi_cg_edf_ls, jacobi_rotations, rotate_orbitals, scdm_qrfact, zij_matrix
109 : USE qs_matrix_pools, ONLY: mpools_get
110 : USE qs_moments, ONLY: build_local_moment_matrix
111 : USE qs_subsys_types, ONLY: qs_subsys_get,&
112 : qs_subsys_type
113 : USE string_utilities, ONLY: xstring
114 : #include "./base/base_uses.f90"
115 :
116 : IMPLICIT NONE
117 :
118 : PRIVATE
119 :
120 : ! *** Global parameters ***
121 :
122 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_methods'
123 :
124 : ! *** Public ***
125 : PUBLIC :: qs_print_cubes, centers_spreads_berry, centers_second_moments_loc, &
126 : centers_second_moments_berry, jacobi_rotation_pipek, optimize_loc_berry, &
127 : optimize_loc_pipek
128 :
129 : CONTAINS
130 :
131 : ! **************************************************************************************************
132 : !> \brief Calculate and optimize the spread functional as calculated from
133 : !> the selected mos and by the definition using the berry phase
134 : !> as given by silvestrelli et al
135 : !> If required the centers and the spreads for each mos selected
136 : !> are calculated from z_ij and printed to file.
137 : !> The centers files is appended if already exists
138 : !> \param method indicates localization algorithm
139 : !> \param qs_loc_env new environment for the localization calculations
140 : !> \param vectors selected mos to be localized
141 : !> \param op_sm_set sparse matrices containing the integrals of the kind <mi e{iGr} nu>
142 : !> \param zij_fm_set set of full matrix of size nmoloc x nmoloc, will contain the z_ij numbers
143 : !> as defined by Silvestrelli et al
144 : !> \param para_env ...
145 : !> \param cell ...
146 : !> \param weights ...
147 : !> \param ispin ...
148 : !> \param print_loc_section ...
149 : !> \param restricted ...
150 : !> \param nextra ...
151 : !> \param nmo ...
152 : !> \param vectors_2 ...
153 : !> \param guess_mos ...
154 : !> \par History
155 : !> 04.2005 created [MI]
156 : !> \author MI
157 : !> \note
158 : !> This definition need the use of complex numbers, therefore the
159 : !> optimization routines are specific for this case
160 : !> The file for the centers and the spreads have a xyz format
161 : ! **************************************************************************************************
162 500 : SUBROUTINE optimize_loc_berry(method, qs_loc_env, vectors, op_sm_set, &
163 500 : zij_fm_set, para_env, cell, weights, ispin, print_loc_section, &
164 : restricted, nextra, nmo, vectors_2, guess_mos)
165 :
166 : INTEGER, INTENT(IN) :: method
167 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
168 : TYPE(cp_fm_type), INTENT(IN) :: vectors
169 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: op_sm_set
170 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
171 : TYPE(mp_para_env_type), POINTER :: para_env
172 : TYPE(cell_type), POINTER :: cell
173 : REAL(dp), DIMENSION(:) :: weights
174 : INTEGER, INTENT(IN) :: ispin
175 : TYPE(section_vals_type), POINTER :: print_loc_section
176 : INTEGER :: restricted
177 : INTEGER, INTENT(IN), OPTIONAL :: nextra, nmo
178 : TYPE(cp_fm_type), INTENT(IN), OPTIONAL :: vectors_2, guess_mos
179 :
180 : CHARACTER(len=*), PARAMETER :: routineN = 'optimize_loc_berry'
181 :
182 : INTEGER :: handle, idim1, idim2, max_iter, nao, &
183 : nmoloc, out_each, output_unit, sweeps
184 : LOGICAL :: converged, crazy_use_diag, &
185 : do_jacobi_refinement, my_do_mixed
186 : REAL(dp) :: crazy_scale, eps_localization, &
187 : max_crazy_angle, start_time, &
188 : target_time
189 500 : TYPE(cp_cfm_type), ALLOCATABLE, DIMENSION(:, :) :: my_zij_cfm_set
190 : TYPE(cp_cfm_type), POINTER :: complex_vectors
191 : TYPE(cp_logger_type), POINTER :: logger
192 :
193 500 : NULLIFY (complex_vectors)
194 :
195 500 : CALL timeset(routineN, handle)
196 500 : logger => cp_get_default_logger()
197 : output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
198 500 : extension=".locInfo")
199 :
200 : ! get rows and cols of the input
201 500 : CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc)
202 :
203 500 : CALL zij_matrix(vectors, op_sm_set, zij_fm_set)
204 :
205 500 : max_iter = qs_loc_env%localized_wfn_control%max_iter
206 500 : max_crazy_angle = qs_loc_env%localized_wfn_control%max_crazy_angle
207 500 : crazy_use_diag = qs_loc_env%localized_wfn_control%crazy_use_diag
208 500 : crazy_scale = qs_loc_env%localized_wfn_control%crazy_scale
209 500 : eps_localization = qs_loc_env%localized_wfn_control%eps_localization
210 500 : out_each = qs_loc_env%localized_wfn_control%out_each
211 500 : target_time = qs_loc_env%target_time
212 500 : start_time = qs_loc_env%start_time
213 500 : do_jacobi_refinement = qs_loc_env%localized_wfn_control%jacobi_refinement
214 500 : my_do_mixed = qs_loc_env%localized_wfn_control%do_mixed
215 :
216 : CALL centers_spreads_berry(qs_loc_env, nmoloc, cell, weights, &
217 500 : ispin, print_loc_section, zij=zij_fm_set, only_initial_out=.TRUE.)
218 :
219 500 : sweeps = 0
220 :
221 822 : SELECT CASE (method)
222 : CASE (do_loc_jacobi)
223 : CALL jacobi_rotations(weights, zij_fm_set, vectors, para_env, max_iter=max_iter, &
224 : eps_localization=eps_localization, sweeps=sweeps, &
225 : out_each=out_each, target_time=target_time, start_time=start_time, &
226 322 : restricted=restricted)
227 : CASE (do_loc_gapo)
228 2 : IF (my_do_mixed) THEN
229 2 : IF (nextra > 0) THEN
230 2 : IF (PRESENT(guess_mos)) THEN
231 : CALL jacobi_cg_edf_ls(para_env, weights, zij_fm_set, vectors, max_iter, &
232 : eps_localization, sweeps, out_each, nextra, &
233 : qs_loc_env%localized_wfn_control%do_cg_po, &
234 2 : nmo=nmo, vectors_2=vectors_2, mos_guess=guess_mos)
235 : ELSE
236 : CALL jacobi_cg_edf_ls(para_env, weights, zij_fm_set, vectors, max_iter, &
237 : eps_localization, sweeps, out_each, nextra, &
238 : qs_loc_env%localized_wfn_control%do_cg_po, &
239 0 : nmo=nmo, vectors_2=vectors_2)
240 : END IF
241 : ELSE
242 : CALL jacobi_cg_edf_ls(para_env, weights, zij_fm_set, vectors, max_iter, &
243 : eps_localization, sweeps, out_each, 0, &
244 0 : qs_loc_env%localized_wfn_control%do_cg_po)
245 : END IF
246 : ELSE
247 0 : CPABORT("GAPO works only with STATES MIXED")
248 : END IF
249 : CASE (do_loc_scdm)
250 : ! Decomposition
251 2 : CALL scdm_qrfact(vectors)
252 : ! Calculation of Zij
253 2 : CALL zij_matrix(vectors, op_sm_set, zij_fm_set)
254 2 : IF (do_jacobi_refinement) THEN
255 : ! Intermediate spread and centers
256 : CALL centers_spreads_berry(qs_loc_env, nmoloc, cell, weights, &
257 0 : ispin, print_loc_section, zij=zij_fm_set, only_initial_out=.TRUE.)
258 : CALL jacobi_rotations(weights, zij_fm_set, vectors, para_env, max_iter=max_iter, &
259 : eps_localization=eps_localization, sweeps=sweeps, &
260 : out_each=out_each, target_time=target_time, start_time=start_time, &
261 0 : restricted=restricted)
262 : END IF
263 : CASE (do_loc_crazy)
264 : CALL crazy_rotations(weights, zij_fm_set, vectors, max_iter=max_iter, max_crazy_angle=max_crazy_angle, &
265 : crazy_scale=crazy_scale, crazy_use_diag=crazy_use_diag, &
266 136 : eps_localization=eps_localization, iterations=sweeps, converged=converged)
267 : ! Possibly fallback to jacobi if the crazy rotation fails
268 136 : IF (.NOT. converged) THEN
269 68 : IF (qs_loc_env%localized_wfn_control%jacobi_fallback) THEN
270 68 : IF (output_unit > 0) WRITE (output_unit, "(T4,A,I6,/,T4,A)") &
271 34 : " Crazy Wannier localization not converged after ", sweeps, &
272 68 : " iterations, switching to jacobi rotations"
273 : CALL jacobi_rotations(weights, zij_fm_set, vectors, para_env, max_iter=max_iter, &
274 : eps_localization=eps_localization, sweeps=sweeps, &
275 : out_each=out_each, target_time=target_time, start_time=start_time, &
276 68 : restricted=restricted)
277 : ELSE
278 0 : IF (output_unit > 0) WRITE (output_unit, "(T4,A,I6,/,T4,A)") &
279 0 : " Crazy Wannier localization not converged after ", sweeps, &
280 0 : " iterations, and jacobi_fallback switched off "
281 : END IF
282 : END IF
283 : CASE (do_loc_direct)
284 : CALL direct_mini(weights, zij_fm_set, vectors, max_iter=max_iter, &
285 2 : eps_localization=eps_localization, iterations=sweeps)
286 : CASE (do_loc_l1_norm_sd)
287 30 : IF (.NOT. cell%orthorhombic) THEN
288 0 : CPABORT("Non-orthorhombic cell with the selected method NYI")
289 : ELSE
290 30 : CALL approx_l1_norm_sd(vectors, max_iter, eps_localization, converged, sweeps)
291 : ! here we need to set zij for the computation of the centers and spreads
292 30 : CALL zij_matrix(vectors, op_sm_set, zij_fm_set)
293 : END IF
294 : CASE (do_loc_cs)
295 78 : ALLOCATE (my_zij_cfm_set(SIZE(zij_fm_set, 1), SIZE(zij_fm_set, 2)))
296 18 : DO idim1 = 1, SIZE(zij_fm_set, 1)
297 54 : DO idim2 = 1, SIZE(zij_fm_set, 2)
298 36 : CALL cp_cfm_create(my_zij_cfm_set(idim1, idim2), zij_fm_set(idim1, idim2)%matrix_struct)
299 48 : CALL cp_fm_to_cfm(msourcer=zij_fm_set(idim1, idim2), mtarget=my_zij_cfm_set(idim1, idim2))
300 : END DO
301 : END DO
302 6 : ALLOCATE (complex_vectors)
303 6 : CALL cp_cfm_create(complex_vectors, vectors%matrix_struct)
304 6 : CALL cp_fm_to_cfm(msourcer=vectors, mtarget=complex_vectors)
305 :
306 : CALL cardoso_souloumiac(weights, my_zij_cfm_set, max_iter, eps_localization, sweeps, &
307 6 : out_each, complex_vectors)
308 :
309 6 : CALL cp_cfm_to_fm(complex_vectors, mtargetr=vectors)
310 6 : CALL cp_cfm_release(complex_vectors)
311 6 : DEALLOCATE (complex_vectors)
312 18 : DO idim1 = 1, SIZE(my_zij_cfm_set, 1)
313 54 : DO idim2 = 1, SIZE(my_zij_cfm_set, 2)
314 48 : CALL cp_cfm_release(my_zij_cfm_set(idim1, idim2))
315 : END DO
316 : END DO
317 12 : DEALLOCATE (my_zij_cfm_set)
318 :
319 : CASE (do_loc_none)
320 0 : IF (output_unit > 0) THEN
321 0 : WRITE (output_unit, '(T4,A,I6,A)') " No MOS localization applied "
322 : END IF
323 : CASE DEFAULT
324 500 : CPABORT("Unknown localization method")
325 : END SELECT
326 500 : IF (output_unit > 0) THEN
327 438 : IF (sweeps <= max_iter) WRITE (output_unit, '(T4,A,I3,A,I6,A)') " Localization for spin ", ispin, &
328 376 : " converged in ", sweeps, " iterations"
329 : END IF
330 :
331 : CALL centers_spreads_berry(qs_loc_env, nmoloc, cell, weights, &
332 500 : ispin, print_loc_section, zij=zij_fm_set)
333 :
334 500 : CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
335 :
336 500 : CALL timestop(handle)
337 :
338 500 : END SUBROUTINE optimize_loc_berry
339 :
340 : ! **************************************************************************************************
341 : !> \brief ...
342 : !> \param qs_env ...
343 : !> \param method ...
344 : !> \param qs_loc_env ...
345 : !> \param vectors ...
346 : !> \param zij_fm_set ...
347 : !> \param ispin ...
348 : !> \param print_loc_section ...
349 : !> \par History
350 : !> 04.2005 created [MI]
351 : !> \author MI
352 : ! **************************************************************************************************
353 6 : SUBROUTINE optimize_loc_pipek(qs_env, method, qs_loc_env, vectors, zij_fm_set, &
354 : ispin, print_loc_section)
355 : TYPE(qs_environment_type), POINTER :: qs_env
356 : INTEGER, INTENT(IN) :: method
357 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
358 : TYPE(cp_fm_type), INTENT(IN) :: vectors
359 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
360 : INTEGER, INTENT(IN) :: ispin
361 : TYPE(section_vals_type), POINTER :: print_loc_section
362 :
363 : CHARACTER(len=*), PARAMETER :: routineN = 'optimize_loc_pipek'
364 :
365 : INTEGER :: handle, iatom, idim1, idim2, isgf, ldz, &
366 : max_iter, nao, natom, ncol, nmoloc, &
367 : out_each, output_unit, sweeps
368 : INTEGER, ALLOCATABLE, DIMENSION(:) :: first_sgf, last_sgf, nsgf
369 : REAL(dp) :: eps
370 6 : TYPE(cp_cfm_type), POINTER :: complex_vectors, my_zij_cfm_set(:, :)
371 6 : TYPE(cp_fm_pool_p_type), DIMENSION(:), POINTER :: ao_ao_fm_pools
372 : TYPE(cp_fm_type) :: opvec
373 : TYPE(cp_fm_type), POINTER :: ov_fm
374 : TYPE(cp_logger_type), POINTER :: logger
375 6 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
376 6 : TYPE(molecule_type), DIMENSION(:), POINTER :: mols
377 6 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
378 6 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
379 :
380 6 : NULLIFY (mols, ov_fm)
381 :
382 6 : CALL timeset(routineN, handle)
383 6 : logger => cp_get_default_logger()
384 : output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
385 6 : extension=".locInfo")
386 :
387 6 : NULLIFY (particle_set)
388 : ! get rows and cols of the input
389 6 : CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc)
390 :
391 : ! replicate the input kind of matrix
392 6 : CALL cp_fm_create(opvec, vectors%matrix_struct)
393 6 : CALL cp_fm_set_all(opvec, 0.0_dp)
394 :
395 : CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, &
396 6 : particle_set=particle_set, qs_kind_set=qs_kind_set, molecule_set=mols)
397 :
398 : ! localization options
399 6 : max_iter = qs_loc_env%localized_wfn_control%max_iter
400 6 : eps = qs_loc_env%localized_wfn_control%eps_localization
401 6 : out_each = qs_loc_env%localized_wfn_control%out_each
402 :
403 6 : natom = SIZE(particle_set, 1)
404 18 : ALLOCATE (first_sgf(natom))
405 12 : ALLOCATE (last_sgf(natom))
406 12 : ALLOCATE (nsgf(natom))
407 :
408 : ! construction of
409 : CALL get_particle_set(particle_set, qs_kind_set, &
410 6 : first_sgf=first_sgf, last_sgf=last_sgf, nsgf=nsgf)
411 :
412 : ! Copy the overlap sparse matrix in a full matrix
413 6 : CALL mpools_get(qs_env%mpools, ao_ao_fm_pools=ao_ao_fm_pools)
414 6 : ALLOCATE (ov_fm)
415 6 : CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, ov_fm, name=" ")
416 6 : CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, ov_fm)
417 :
418 : ! Compute zij here
419 78 : DO iatom = 1, natom
420 72 : CALL cp_fm_set_all(zij_fm_set(iatom, 1), 0.0_dp)
421 72 : CALL cp_fm_get_info(zij_fm_set(iatom, 1), ncol_global=ldz)
422 72 : isgf = first_sgf(iatom)
423 72 : ncol = nsgf(iatom)
424 :
425 : ! multiply fmxfm, using only part of the ao : Ct x S
426 : CALL parallel_gemm('N', 'N', nao, nmoloc, nao, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, &
427 72 : a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1)
428 :
429 : CALL parallel_gemm('T', 'N', nmoloc, nmoloc, ncol, 0.5_dp, vectors, opvec, &
430 : 0.0_dp, zij_fm_set(iatom, 1), &
431 72 : a_first_col=1, a_first_row=isgf, b_first_col=1, b_first_row=isgf)
432 :
433 : CALL parallel_gemm('N', 'N', nao, nmoloc, ncol, 1.0_dp, ov_fm, vectors, 0.0_dp, opvec, &
434 72 : a_first_col=isgf, a_first_row=1, b_first_col=1, b_first_row=isgf)
435 :
436 : CALL parallel_gemm('T', 'N', nmoloc, nmoloc, nao, 0.5_dp, vectors, opvec, &
437 : 1.0_dp, zij_fm_set(iatom, 1), &
438 150 : a_first_col=1, a_first_row=1, b_first_col=1, b_first_row=1)
439 :
440 : END DO ! iatom
441 :
442 : ! And now perform the optimization and rotate the orbitals
443 6 : SELECT CASE (method)
444 : CASE (do_loc_jacobi)
445 0 : CALL jacobi_rotation_pipek(zij_fm_set, vectors, sweeps)
446 : CASE (do_loc_gapo)
447 0 : CPABORT("GAPO and Pipek not implemented.")
448 : CASE (do_loc_crazy)
449 0 : CPABORT("Crazy and Pipek not implemented.")
450 : CASE (do_loc_l1_norm_sd)
451 0 : CPABORT("L1 norm and Pipek not implemented.")
452 : CASE (do_loc_direct)
453 0 : CPABORT("Direct and Pipek not implemented.")
454 : CASE (do_loc_cs)
455 102 : ALLOCATE (my_zij_cfm_set(SIZE(zij_fm_set, 1), SIZE(zij_fm_set, 2)))
456 78 : DO idim1 = 1, SIZE(zij_fm_set, 1)
457 150 : DO idim2 = 1, SIZE(zij_fm_set, 2)
458 72 : CALL cp_cfm_create(my_zij_cfm_set(idim1, idim2), zij_fm_set(idim1, idim2)%matrix_struct)
459 144 : CALL cp_fm_to_cfm(msourcer=zij_fm_set(idim1, idim2), mtarget=my_zij_cfm_set(idim1, idim2))
460 : END DO
461 : END DO
462 6 : ALLOCATE (complex_vectors)
463 6 : CALL cp_cfm_create(complex_vectors, vectors%matrix_struct)
464 6 : CALL cp_fm_to_cfm(msourcer=vectors, mtarget=complex_vectors)
465 :
466 6 : CALL cardoso_souloumiac_pipek(my_zij_cfm_set, complex_vectors, sweeps, max_iter, eps, out_each)
467 :
468 6 : CALL cp_cfm_to_fm(complex_vectors, mtargetr=vectors)
469 6 : CALL cp_cfm_release(complex_vectors)
470 6 : DEALLOCATE (complex_vectors)
471 78 : DO idim1 = 1, SIZE(my_zij_cfm_set, 1)
472 150 : DO idim2 = 1, SIZE(my_zij_cfm_set, 2)
473 144 : CALL cp_cfm_release(my_zij_cfm_set(idim1, idim2))
474 : END DO
475 : END DO
476 12 : DEALLOCATE (my_zij_cfm_set)
477 :
478 : CASE (do_loc_none)
479 0 : IF (output_unit > 0) WRITE (output_unit, '(A,I6,A)') " No MOS localization applied "
480 : CASE DEFAULT
481 6 : CPABORT("Unknown localization method")
482 : END SELECT
483 :
484 9 : IF (output_unit > 0) WRITE (output_unit, '(A,I3,A,I6,A)') " Localization for spin ", ispin, &
485 6 : " converged in ", sweeps, " iterations"
486 :
487 : CALL centers_spreads_pipek(qs_loc_env, zij_fm_set, particle_set, ispin, &
488 6 : print_loc_section)
489 :
490 6 : DEALLOCATE (first_sgf, last_sgf, nsgf)
491 :
492 6 : CALL cp_fm_release(opvec)
493 6 : CALL cp_fm_release(ov_fm)
494 6 : DEALLOCATE (ov_fm)
495 6 : CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
496 :
497 6 : CALL timestop(handle)
498 :
499 18 : END SUBROUTINE optimize_loc_pipek
500 :
501 : ! **************************************************************************************************
502 : !> \brief 2by2 rotation for the pipek operator
503 : !> in this case the z_ij numbers are reals
504 : !> \param zij_fm_set ...
505 : !> \param vectors ...
506 : !> \param sweeps ...
507 : !> \par History
508 : !> 05-2005 created
509 : !> \author MI
510 : ! **************************************************************************************************
511 0 : SUBROUTINE jacobi_rotation_pipek(zij_fm_set, vectors, sweeps)
512 :
513 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
514 : TYPE(cp_fm_type), INTENT(IN) :: vectors
515 : INTEGER :: sweeps
516 :
517 : INTEGER :: iatom, istate, jstate, natom, nstate
518 : REAL(KIND=dp) :: aij, bij, ct, mii, mij, mjj, ratio, st, &
519 : theta, tolerance
520 : TYPE(cp_fm_type) :: rmat
521 :
522 0 : CALL cp_fm_create(rmat, zij_fm_set(1, 1)%matrix_struct)
523 0 : CALL cp_fm_set_all(rmat, 0.0_dp, 1.0_dp)
524 :
525 0 : CALL cp_fm_get_info(rmat, nrow_global=nstate)
526 0 : tolerance = 1.0e10_dp
527 0 : sweeps = 0
528 0 : natom = SIZE(zij_fm_set, 1)
529 : ! do jacobi sweeps until converged
530 0 : DO WHILE (tolerance >= 1.0e-4_dp)
531 0 : sweeps = sweeps + 1
532 0 : DO istate = 1, nstate
533 0 : DO jstate = istate + 1, nstate
534 : aij = 0.0_dp
535 : bij = 0.0_dp
536 0 : DO iatom = 1, natom
537 0 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, mii)
538 0 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, jstate, mij)
539 0 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), jstate, jstate, mjj)
540 0 : aij = aij + mij*(mii - mjj)
541 0 : bij = bij + mij*mij - 0.25_dp*(mii - mjj)*(mii - mjj)
542 : END DO
543 0 : IF (ABS(bij) > 1.E-10_dp) THEN
544 0 : ratio = -aij/bij
545 0 : theta = 0.25_dp*ATAN(ratio)
546 : ELSE
547 0 : bij = 0.0_dp
548 : theta = 0.0_dp
549 : END IF
550 : ! Check max or min
551 : ! To minimize the spread
552 0 : IF (theta > pi*0.5_dp) THEN
553 0 : theta = theta - pi*0.25_dp
554 0 : ELSE IF (theta < -pi*0.5_dp) THEN
555 0 : theta = theta + pi*0.25_dp
556 : END IF
557 :
558 0 : ct = COS(theta)
559 0 : st = SIN(theta)
560 :
561 0 : CALL cp_fm_rot_cols(rmat, istate, jstate, ct, st)
562 :
563 0 : DO iatom = 1, natom
564 0 : CALL cp_fm_rot_cols(zij_fm_set(iatom, 1), istate, jstate, ct, st)
565 0 : CALL cp_fm_rot_rows(zij_fm_set(iatom, 1), istate, jstate, ct, st)
566 : END DO
567 : END DO
568 : END DO
569 0 : CALL check_tolerance_real(zij_fm_set, tolerance)
570 : END DO
571 :
572 0 : CALL rotate_orbitals(rmat, vectors)
573 0 : CALL cp_fm_release(rmat)
574 :
575 0 : END SUBROUTINE jacobi_rotation_pipek
576 :
577 : ! **************************************************************************************************
578 : !> \brief ...
579 : !> \param zij_fm_set ...
580 : !> \param tolerance ...
581 : !> \par History
582 : !> 04.2005 created [MI]
583 : !> \author MI
584 : ! **************************************************************************************************
585 0 : SUBROUTINE check_tolerance_real(zij_fm_set, tolerance)
586 :
587 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
588 : REAL(dp), INTENT(OUT) :: tolerance
589 :
590 : INTEGER :: iatom, istate, jstate, natom, &
591 : ncol_local, nrow_global, nrow_local
592 0 : INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
593 : REAL(dp) :: grad_ij, zii, zij, zjj
594 0 : REAL(dp), DIMENSION(:, :), POINTER :: diag
595 : TYPE(cp_fm_type) :: force
596 :
597 0 : CALL cp_fm_create(force, zij_fm_set(1, 1)%matrix_struct)
598 0 : CALL cp_fm_set_all(force, 0._dp)
599 :
600 0 : NULLIFY (diag, col_indices, row_indices)
601 0 : natom = SIZE(zij_fm_set, 1)
602 : CALL cp_fm_get_info(zij_fm_set(1, 1), nrow_local=nrow_local, &
603 : ncol_local=ncol_local, nrow_global=nrow_global, &
604 0 : row_indices=row_indices, col_indices=col_indices)
605 0 : ALLOCATE (diag(nrow_global, natom))
606 :
607 0 : DO iatom = 1, natom
608 0 : DO istate = 1, nrow_global
609 0 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, diag(istate, iatom))
610 : END DO
611 : END DO
612 :
613 0 : DO istate = 1, nrow_local
614 0 : DO jstate = 1, ncol_local
615 : grad_ij = 0.0_dp
616 0 : DO iatom = 1, natom
617 0 : zii = diag(row_indices(istate), iatom)
618 0 : zjj = diag(col_indices(jstate), iatom)
619 0 : zij = zij_fm_set(iatom, 1)%local_data(istate, jstate)
620 0 : grad_ij = grad_ij + 4.0_dp*zij*(zjj - zii)
621 : END DO
622 0 : force%local_data(istate, jstate) = grad_ij
623 : END DO
624 : END DO
625 :
626 0 : DEALLOCATE (diag)
627 :
628 0 : CALL cp_fm_maxabsval(force, tolerance)
629 0 : CALL cp_fm_release(force)
630 :
631 0 : END SUBROUTINE check_tolerance_real
632 : ! **************************************************************************************************
633 : !> \brief ...
634 : !> \param qs_loc_env ...
635 : !> \param nmoloc ...
636 : !> \param cell ...
637 : !> \param weights ...
638 : !> \param ispin ...
639 : !> \param print_loc_section ...
640 : !> \param zij ...
641 : !> \param c_zij ...
642 : !> \param only_initial_out ...
643 : !> \par History
644 : !> 04.2005 created [MI]
645 : !> \author MI
646 : ! **************************************************************************************************
647 1112 : SUBROUTINE centers_spreads_berry(qs_loc_env, nmoloc, cell, weights, ispin, &
648 1112 : print_loc_section, zij, c_zij, only_initial_out)
649 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
650 : INTEGER, INTENT(IN) :: nmoloc
651 : TYPE(cell_type), POINTER :: cell
652 : REAL(dp), DIMENSION(:) :: weights
653 : INTEGER, INTENT(IN) :: ispin
654 : TYPE(section_vals_type), POINTER :: print_loc_section
655 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN), &
656 : OPTIONAL :: zij
657 : TYPE(cp_cfm_p_type), INTENT(INOUT), OPTIONAL :: c_zij(:, :)
658 : LOGICAL, INTENT(IN), OPTIONAL :: only_initial_out
659 :
660 : CHARACTER(len=default_path_length) :: file_tmp, iter
661 : COMPLEX(KIND=dp) :: z
662 : INTEGER :: idir, istate, jdir, nstates, &
663 : output_unit, unit_out_s
664 : LOGICAL :: my_only_init, zij_is_complex
665 : REAL(dp) :: avg_spread_ii, spread_i, spread_ii, &
666 : sum_spread_i, sum_spread_ii
667 : REAL(dp), DIMENSION(3) :: c, c2, cpbc
668 1112 : REAL(dp), DIMENSION(:, :), POINTER :: centers
669 : REAL(KIND=dp) :: imagpart, realpart
670 : TYPE(cp_logger_type), POINTER :: logger
671 : TYPE(section_vals_type), POINTER :: print_key
672 :
673 1112 : NULLIFY (centers, logger, print_key)
674 2224 : logger => cp_get_default_logger()
675 1112 : my_only_init = .FALSE.
676 1112 : IF (PRESENT(only_initial_out)) my_only_init = only_initial_out
677 1112 : IF (PRESENT(zij)) THEN
678 1112 : zij_is_complex = .FALSE.
679 1112 : CPASSERT(.NOT. PRESENT(c_zij))
680 1112 : CALL cp_fm_get_info(zij(1, 1), nrow_global=nstates)
681 : ELSE
682 0 : zij_is_complex = .TRUE.
683 0 : CPASSERT(PRESENT(c_zij))
684 0 : CALL cp_cfm_get_info(c_zij(1, 1)%matrix, nrow_global=nstates)
685 : END IF
686 :
687 1112 : file_tmp = TRIM(qs_loc_env%tag_mo)//"_spreads_s"//TRIM(ADJUSTL(cp_to_string(ispin)))
688 : output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
689 1112 : extension=".locInfo")
690 : unit_out_s = cp_print_key_unit_nr(logger, print_loc_section, "WANNIER_SPREADS", &
691 1112 : middle_name=file_tmp, extension=".data")
692 1112 : iter = cp_iter_string(logger%iter_info)
693 1112 : IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(i6,/,A)') nmoloc, TRIM(iter)
694 :
695 1112 : CPASSERT(nstates >= nmoloc)
696 :
697 1112 : centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
698 1112 : CPASSERT(ASSOCIATED(centers))
699 1112 : CPASSERT(SIZE(centers, 2) == nmoloc)
700 1112 : sum_spread_i = 0.0_dp
701 1112 : sum_spread_ii = 0.0_dp
702 1112 : avg_spread_ii = 0.0_dp
703 9078 : DO istate = 1, nmoloc
704 7966 : c = 0.0_dp
705 : c2 = 0.0_dp
706 7966 : spread_i = 0.0_dp
707 7966 : spread_ii = 0.0_dp
708 7966 : IF (.NOT. zij_is_complex) THEN
709 32416 : DO jdir = 1, SIZE(zij, 2)
710 24450 : CALL cp_fm_get_element(zij(1, jdir), istate, istate, realpart)
711 24450 : CALL cp_fm_get_element(zij(2, jdir), istate, istate, imagpart)
712 24450 : z = CMPLX(realpart, imagpart, dp)
713 : spread_i = spread_i - weights(jdir)* &
714 24450 : LOG(realpart*realpart + imagpart*imagpart)/twopi/twopi
715 : spread_ii = spread_ii + weights(jdir)* &
716 24450 : (1.0_dp - (realpart*realpart + imagpart*imagpart))/twopi/twopi
717 32416 : IF (jdir < 4) THEN
718 95592 : DO idir = 1, 3
719 : c(idir) = c(idir) + &
720 95592 : (cell%hmat(idir, jdir)/twopi)*AIMAG(LOG(z))
721 : END DO
722 : END IF
723 : END DO
724 : ELSE
725 0 : DO jdir = 1, SIZE(c_zij, 2)
726 0 : CALL cp_cfm_get_element(c_zij(1, jdir)%matrix, istate, istate, z)
727 0 : realpart = REAL(z, dp)
728 0 : imagpart = AIMAG(z)
729 : spread_i = spread_i - weights(jdir)* &
730 0 : LOG(realpart*realpart + imagpart*imagpart)/twopi/twopi
731 : spread_ii = spread_ii + weights(jdir)* &
732 0 : (1.0_dp - (realpart*realpart + imagpart*imagpart))/twopi/twopi
733 0 : IF (jdir < 4) THEN
734 0 : DO idir = 1, 3
735 : c(idir) = c(idir) + &
736 0 : (cell%hmat(idir, jdir)/twopi)*AIMAG(LOG(z))
737 : END DO
738 : END IF
739 : END DO
740 : END IF
741 7966 : cpbc = pbc(c, cell)
742 31864 : centers(1:3, istate) = cpbc(1:3)
743 7966 : centers(4, istate) = spread_i
744 7966 : centers(5, istate) = spread_ii
745 7966 : sum_spread_i = sum_spread_i + centers(4, istate)
746 7966 : sum_spread_ii = sum_spread_ii + centers(5, istate)
747 9134 : IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(I6,2F16.8)') istate, centers(4:5, istate)
748 : END DO
749 1112 : avg_spread_ii = sum_spread_ii/REAL(nmoloc, KIND=dp)
750 :
751 : ! Print of wannier centers
752 1112 : print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CENTERS")
753 1112 : IF (.NOT. my_only_init) CALL print_wannier_centers(qs_loc_env, print_key, centers, logger, ispin)
754 :
755 1112 : IF (output_unit > 0) THEN
756 500 : WRITE (output_unit, '(T4, A, 2x, 2A26,/,T23, A28)') " Spread Functional ", "sum_in -w_i ln(|z_in|^2)", &
757 1000 : "sum_in w_i(1-|z_in|^2)", "sum_in w_i(1-|z_in|^2)/n"
758 500 : IF (my_only_init) THEN
759 250 : WRITE (output_unit, '(T4,A,T38,2F20.10,/,T38,F20.10)') " Initial Spread (Berry) : ", &
760 500 : sum_spread_i, sum_spread_ii, avg_spread_ii
761 : ELSE
762 250 : WRITE (output_unit, '(T4,A,T38,2F20.10,/,T38,F20.10)') " Total Spread (Berry) : ", &
763 500 : sum_spread_i, sum_spread_ii, avg_spread_ii
764 : END IF
765 : END IF
766 :
767 1112 : IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(A,2F16.10)') " Total ", sum_spread_i, sum_spread_ii
768 :
769 1112 : CALL cp_print_key_finished_output(unit_out_s, logger, print_loc_section, "WANNIER_SPREADS")
770 1112 : CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
771 :
772 1112 : END SUBROUTINE centers_spreads_berry
773 :
774 : ! **************************************************************************************************
775 : !> \brief define and print the centers and spread
776 : !> when the pipek operator is used
777 : !> \param qs_loc_env ...
778 : !> \param zij_fm_set matrix elements that define the populations on atoms
779 : !> \param particle_set ...
780 : !> \param ispin spin 1 or 2
781 : !> \param print_loc_section ...
782 : !> \par History
783 : !> 05.2005 created [MI]
784 : ! **************************************************************************************************
785 6 : SUBROUTINE centers_spreads_pipek(qs_loc_env, zij_fm_set, particle_set, ispin, &
786 : print_loc_section)
787 :
788 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
789 : TYPE(cp_fm_type), DIMENSION(:, :), INTENT(IN) :: zij_fm_set
790 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
791 : INTEGER, INTENT(IN) :: ispin
792 : TYPE(section_vals_type), POINTER :: print_loc_section
793 :
794 : CHARACTER(len=default_path_length) :: file_tmp, iter
795 : INTEGER :: iatom, istate, natom, nstate, unit_out_s
796 6 : INTEGER, DIMENSION(:), POINTER :: atom_of_state
797 : REAL(dp) :: r(3)
798 6 : REAL(dp), ALLOCATABLE, DIMENSION(:) :: Qii, ziimax
799 6 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: diag
800 6 : REAL(dp), DIMENSION(:, :), POINTER :: centers
801 : TYPE(cp_logger_type), POINTER :: logger
802 : TYPE(section_vals_type), POINTER :: print_key
803 :
804 6 : NULLIFY (centers, logger, print_key)
805 12 : logger => cp_get_default_logger()
806 :
807 6 : CALL cp_fm_get_info(zij_fm_set(1, 1), nrow_global=nstate)
808 6 : natom = SIZE(zij_fm_set, 1)
809 :
810 6 : centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
811 6 : CPASSERT(ASSOCIATED(centers))
812 6 : CPASSERT(SIZE(centers, 2) == nstate)
813 :
814 6 : file_tmp = TRIM(qs_loc_env%tag_mo)//"_spreads_s"//TRIM(ADJUSTL(cp_to_string(ispin)))
815 : unit_out_s = cp_print_key_unit_nr(logger, print_loc_section, "WANNIER_SPREADS", &
816 6 : middle_name=file_tmp, extension=".data", log_filename=.FALSE.)
817 6 : iter = cp_iter_string(logger%iter_info)
818 6 : IF (unit_out_s > 0) WRITE (unit_out_s, '(i6,/,A)') TRIM(iter)
819 :
820 18 : ALLOCATE (atom_of_state(nstate))
821 198 : atom_of_state = 0
822 :
823 24 : ALLOCATE (diag(nstate, natom))
824 6 : diag = 0.0_dp
825 :
826 78 : DO iatom = 1, natom
827 2382 : DO istate = 1, nstate
828 2376 : CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, diag(istate, iatom))
829 : END DO
830 : END DO
831 :
832 24 : ALLOCATE (Qii(nstate), ziimax(nstate))
833 6 : ziimax = 0.0_dp
834 6 : Qii = 0.0_dp
835 :
836 78 : DO iatom = 1, natom
837 2382 : DO istate = 1, nstate
838 2304 : Qii(istate) = Qii(istate) + diag(istate, iatom)*diag(istate, iatom)
839 2376 : IF (ABS(diag(istate, iatom)) > ziimax(istate)) THEN
840 706 : ziimax(istate) = ABS(diag(istate, iatom))
841 706 : atom_of_state(istate) = iatom
842 : END IF
843 : END DO
844 : END DO
845 :
846 198 : DO istate = 1, nstate
847 192 : iatom = atom_of_state(istate)
848 768 : r(1:3) = particle_set(iatom)%r(1:3)
849 768 : centers(1:3, istate) = r(1:3)
850 192 : centers(4, istate) = 1.0_dp/Qii(istate)
851 198 : IF (unit_out_s > 0) WRITE (unit_out_s, '(I6,F16.8)') istate, angstrom*centers(4, istate)
852 : END DO
853 :
854 : ! Print the wannier centers
855 6 : print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CENTERS")
856 6 : CALL print_wannier_centers(qs_loc_env, print_key, centers, logger, ispin)
857 :
858 6 : CALL cp_print_key_finished_output(unit_out_s, logger, print_loc_section, "WANNIER_SPREADS")
859 :
860 6 : DEALLOCATE (Qii, ziimax, atom_of_state, diag)
861 :
862 6 : END SUBROUTINE centers_spreads_pipek
863 :
864 : ! **************************************************************************************************
865 : !> \brief write the cube files for a set of selected states
866 : !> \param qs_env ...
867 : !> \param mo_coeff set mos from which the states to be printed are extracted
868 : !> \param nstates number of states to be printed
869 : !> \param state_list list of the indexes of the states to be printed
870 : !> \param centers centers and spread, all=0 if they hva not been calculated
871 : !> \param print_key ...
872 : !> \param root initial part of the cube file names
873 : !> \param ispin ...
874 : !> \param idir ...
875 : !> \param state0 ...
876 : !> \param file_position ...
877 : !> \par History
878 : !> 08.2005 created [MI]
879 : !> \author MI
880 : !> \note
881 : !> This routine should not be in this module
882 : ! **************************************************************************************************
883 12 : SUBROUTINE qs_print_cubes(qs_env, mo_coeff, nstates, state_list, centers, &
884 : print_key, root, ispin, idir, state0, file_position)
885 :
886 : TYPE(qs_environment_type), POINTER :: qs_env
887 : TYPE(cp_fm_type), INTENT(IN) :: mo_coeff
888 : INTEGER, INTENT(IN) :: nstates
889 : INTEGER, DIMENSION(:), POINTER :: state_list
890 : REAL(dp), DIMENSION(:, :), POINTER :: centers
891 : TYPE(section_vals_type), POINTER :: print_key
892 : CHARACTER(LEN=*) :: root
893 : INTEGER, INTENT(IN), OPTIONAL :: ispin, idir
894 : INTEGER, OPTIONAL :: state0
895 : CHARACTER(LEN=default_string_length), OPTIONAL :: file_position
896 :
897 : CHARACTER(len=*), PARAMETER :: routineN = 'qs_print_cubes'
898 : CHARACTER, DIMENSION(3), PARAMETER :: labels = ['x', 'y', 'z']
899 :
900 : CHARACTER :: label
901 : CHARACTER(LEN=default_path_length) :: file_tmp, filename, my_pos
902 : CHARACTER(LEN=default_string_length) :: title
903 : INTEGER :: handle, ia, ie, istate, ivector, &
904 : my_ispin, my_state0, unit_out_c
905 : LOGICAL :: add_idir, add_spin, mpi_io
906 12 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
907 : TYPE(cell_type), POINTER :: cell
908 : TYPE(cp_logger_type), POINTER :: logger
909 : TYPE(dft_control_type), POINTER :: dft_control
910 : TYPE(particle_list_type), POINTER :: particles
911 12 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
912 : TYPE(pw_c1d_gs_type) :: wf_g
913 : TYPE(pw_env_type), POINTER :: pw_env
914 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
915 : TYPE(pw_r3d_rs_type) :: wf_r
916 12 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
917 : TYPE(qs_subsys_type), POINTER :: subsys
918 :
919 12 : CALL timeset(routineN, handle)
920 12 : logger => cp_get_default_logger()
921 :
922 12 : CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
923 12 : CALL qs_subsys_get(subsys, particles=particles)
924 :
925 12 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
926 :
927 12 : CALL auxbas_pw_pool%create_pw(wf_r)
928 12 : CALL auxbas_pw_pool%create_pw(wf_g)
929 :
930 12 : my_state0 = 0
931 12 : IF (PRESENT(state0)) my_state0 = state0
932 :
933 12 : add_spin = .FALSE.
934 12 : my_ispin = 1
935 12 : IF (PRESENT(ispin)) THEN
936 8 : add_spin = .TRUE.
937 8 : my_ispin = ispin
938 : END IF
939 12 : add_idir = .FALSE.
940 12 : IF (PRESENT(idir)) THEN
941 0 : add_idir = .TRUE.
942 0 : label = labels(idir)
943 : END IF
944 :
945 12 : my_pos = "REWIND"
946 12 : IF (PRESENT(file_position)) THEN
947 12 : my_pos = file_position
948 : END IF
949 :
950 62 : DO istate = 1, nstates
951 50 : ivector = state_list(istate) - my_state0
952 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, cell=cell, &
953 50 : dft_control=dft_control, particle_set=particle_set, pw_env=pw_env)
954 :
955 : CALL calculate_wavefunction(mo_coeff, ivector, wf_r, wf_g, atomic_kind_set, &
956 50 : qs_kind_set, cell, dft_control, particle_set, pw_env)
957 :
958 : ! Formatting the middle part of the name
959 50 : ivector = state_list(istate)
960 50 : CALL xstring(root, ia, ie)
961 50 : IF (add_idir) THEN
962 0 : filename = root(ia:ie)//"_"//label//"_w"//TRIM(ADJUSTL(cp_to_string(ivector)))
963 : ELSE
964 50 : filename = root(ia:ie)//"_w"//TRIM(ADJUSTL(cp_to_string(ivector)))
965 : END IF
966 50 : IF (add_spin) THEN
967 38 : file_tmp = filename
968 38 : CALL xstring(file_tmp, ia, ie)
969 38 : filename = file_tmp(ia:ie)//"_s"//TRIM(ADJUSTL(cp_to_string(ispin)))
970 : END IF
971 :
972 : ! Using the print_key tools to open the file with the proper name
973 50 : mpi_io = .TRUE.
974 : unit_out_c = cp_print_key_unit_nr(logger, print_key, "", middle_name=filename, &
975 : extension=".cube", file_position=my_pos, log_filename=.FALSE., &
976 50 : mpi_io=mpi_io)
977 50 : IF (SIZE(centers, 1) == 6) THEN
978 50 : WRITE (title, '(A7,I5.5,A2,I1.1,A1,6F10.4)') "WFN ", ivector, "_s", my_ispin, " ", &
979 400 : centers(1:3, istate)*angstrom, centers(4:6, istate)*angstrom
980 : ELSE
981 0 : WRITE (title, '(A7,I5.5,A2,I1.1,A1,3F10.4)') "WFN ", ivector, "_s", my_ispin, " ", &
982 0 : centers(1:3, istate)*angstrom
983 : END IF
984 : CALL cp_pw_to_cube(wf_r, unit_out_c, title, &
985 : particles=particles, &
986 50 : stride=section_get_ivals(print_key, "STRIDE"), mpi_io=mpi_io)
987 62 : CALL cp_print_key_finished_output(unit_out_c, logger, print_key, "", mpi_io=mpi_io)
988 : END DO
989 :
990 12 : CALL auxbas_pw_pool%give_back_pw(wf_r)
991 12 : CALL auxbas_pw_pool%give_back_pw(wf_g)
992 12 : CALL timestop(handle)
993 12 : END SUBROUTINE qs_print_cubes
994 :
995 : ! **************************************************************************************************
996 : !> \brief Prints the electronic density
997 : !> \param qs_env ...
998 : !> \param nstates ...
999 : !> \param state_list ...
1000 : !> \param centers ...
1001 : !> \param print_key ...
1002 : !> \param root ...
1003 : !> \param mo_coeff ...
1004 : !> \param complex_mo_coeff ...
1005 : !> \param ispin ...
1006 : !> \param file_position ...
1007 : ! **************************************************************************************************
1008 :
1009 0 : SUBROUTINE qs_print_density_cubes(qs_env, nstates, state_list, centers, print_key, root, &
1010 : mo_coeff, complex_mo_coeff, ispin, file_position)
1011 :
1012 : TYPE(qs_environment_type), POINTER :: qs_env
1013 : INTEGER, INTENT(IN) :: nstates
1014 : INTEGER, DIMENSION(:), POINTER :: state_list
1015 : REAL(dp), DIMENSION(:, :), POINTER :: centers
1016 : TYPE(section_vals_type), POINTER :: print_key
1017 : CHARACTER(LEN=*) :: root
1018 : TYPE(cp_fm_type), OPTIONAL, POINTER :: mo_coeff
1019 : TYPE(cp_cfm_type), OPTIONAL, POINTER :: complex_mo_coeff
1020 : INTEGER, INTENT(IN), OPTIONAL :: ispin
1021 : CHARACTER(LEN=default_string_length), OPTIONAL :: file_position
1022 :
1023 : TYPE(cp_fm_type), POINTER :: squaredmat
1024 :
1025 0 : IF (PRESENT(mo_coeff)) THEN
1026 0 : CPASSERT(.NOT. PRESENT(complex_mo_coeff))
1027 0 : CALL cp_fm_create(squaredmat, mo_coeff%matrix_struct)
1028 0 : squaredmat%local_data = mo_coeff%local_data**2.0
1029 : ELSE
1030 0 : CPASSERT(PRESENT(complex_mo_coeff))
1031 0 : CALL cp_fm_create(squaredmat, complex_mo_coeff%matrix_struct)
1032 : squaredmat%local_data = REAL(CONJG(complex_mo_coeff%local_data) &
1033 0 : *complex_mo_coeff%local_data, dp)
1034 : END IF
1035 : CALL qs_print_cubes(qs_env, squaredmat, nstates, state_list, centers, print_key, root, &
1036 0 : ispin=ispin, file_position=file_position)
1037 0 : CALL cp_fm_release(squaredmat)
1038 :
1039 0 : END SUBROUTINE qs_print_density_cubes
1040 :
1041 : ! **************************************************************************************************
1042 : !> \brief Prints wannier centers
1043 : !> \param qs_loc_env ...
1044 : !> \param print_key ...
1045 : !> \param center ...
1046 : !> \param logger ...
1047 : !> \param ispin ...
1048 : ! **************************************************************************************************
1049 506 : SUBROUTINE print_wannier_centers(qs_loc_env, print_key, center, logger, ispin)
1050 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1051 : TYPE(section_vals_type), POINTER :: print_key
1052 : REAL(KIND=dp), INTENT(IN) :: center(:, :)
1053 : TYPE(cp_logger_type), POINTER :: logger
1054 : INTEGER, INTENT(IN) :: ispin
1055 :
1056 : CHARACTER(default_string_length) :: iter, unit_str
1057 : CHARACTER(LEN=default_string_length) :: my_ext, my_form
1058 : INTEGER :: iunit, l, nstates
1059 : LOGICAL :: first_time, init_traj
1060 : REAL(KIND=dp) :: unit_conv
1061 :
1062 506 : nstates = SIZE(center, 2)
1063 506 : my_form = "formatted"
1064 506 : my_ext = ".data"
1065 506 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, first_time=first_time) &
1066 : , cp_p_file)) THEN
1067 : ! Find out if we want to print IONS+CENTERS or ONLY CENTERS..
1068 140 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, "/IONS+CENTERS") &
1069 : , cp_p_file)) THEN
1070 38 : CALL get_output_format(print_key, my_form=my_form, my_ext=my_ext)
1071 : END IF
1072 140 : IF (first_time .OR. (.NOT. qs_loc_env%first_time)) THEN
1073 : iunit = cp_print_key_unit_nr(logger, print_key, "", extension=my_ext, file_form=my_form, &
1074 : middle_name=TRIM(qs_loc_env%tag_mo)//"_centers_s"//TRIM(ADJUSTL(cp_to_string(ispin))), &
1075 140 : log_filename=.FALSE., on_file=.TRUE., is_new_file=init_traj)
1076 140 : IF (iunit > 0) THEN
1077 : ! Gather units of measure for output (if available)
1078 70 : CALL section_vals_val_get(print_key, "UNIT", c_val=unit_str)
1079 70 : unit_conv = cp_unit_from_cp2k(1.0_dp, TRIM(unit_str))
1080 :
1081 70 : IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, "/IONS+CENTERS"), cp_p_file)) THEN
1082 : ! Different possible formats
1083 19 : CALL print_wannier_traj(qs_loc_env, print_key, center, iunit, init_traj, unit_conv)
1084 : ELSE
1085 : ! Default print format
1086 51 : iter = cp_iter_string(logger%iter_info)
1087 51 : WRITE (iunit, '(i6,/,a)') nstates, TRIM(iter)
1088 423 : DO l = 1, nstates
1089 1911 : WRITE (iunit, '(A,4F16.8)') "X", unit_conv*center(1:4, l)
1090 : END DO
1091 : END IF
1092 : END IF
1093 140 : CALL cp_print_key_finished_output(iunit, logger, print_key, on_file=.TRUE.)
1094 : END IF
1095 : END IF
1096 506 : END SUBROUTINE print_wannier_centers
1097 :
1098 : ! **************************************************************************************************
1099 : !> \brief computes spread and centers
1100 : !> \param qs_loc_env ...
1101 : !> \param print_key ...
1102 : !> \param center ...
1103 : !> \param iunit ...
1104 : !> \param init_traj ...
1105 : !> \param unit_conv ...
1106 : ! **************************************************************************************************
1107 19 : SUBROUTINE print_wannier_traj(qs_loc_env, print_key, center, iunit, init_traj, unit_conv)
1108 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1109 : TYPE(section_vals_type), POINTER :: print_key
1110 : REAL(KIND=dp), INTENT(IN) :: center(:, :)
1111 : INTEGER, INTENT(IN) :: iunit
1112 : LOGICAL, INTENT(IN) :: init_traj
1113 : REAL(KIND=dp), INTENT(IN) :: unit_conv
1114 :
1115 : CHARACTER(len=default_string_length) :: iter, remark1, remark2, title
1116 : INTEGER :: i, iskip, natom, ntot, outformat
1117 19 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1118 : TYPE(atomic_kind_type), POINTER :: atomic_kind
1119 : TYPE(cp_logger_type), POINTER :: logger
1120 19 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1121 :
1122 19 : NULLIFY (particle_set, atomic_kind_set, atomic_kind, logger)
1123 38 : logger => cp_get_default_logger()
1124 19 : natom = SIZE(qs_loc_env%particle_set)
1125 19 : ntot = natom + SIZE(center, 2)
1126 19 : CALL allocate_particle_set(particle_set, ntot)
1127 38 : ALLOCATE (atomic_kind_set(1))
1128 19 : atomic_kind => atomic_kind_set(1)
1129 : CALL set_atomic_kind(atomic_kind=atomic_kind, kind_number=0, &
1130 19 : name="X", element_symbol="X", mass=0.0_dp)
1131 : ! Particles
1132 283 : DO i = 1, natom
1133 264 : particle_set(i)%atomic_kind => qs_loc_env%particle_set(i)%atomic_kind
1134 1075 : particle_set(i)%r = pbc(qs_loc_env%particle_set(i)%r, qs_loc_env%cell)
1135 : END DO
1136 : ! Wannier Centers
1137 452 : DO i = natom + 1, ntot
1138 433 : particle_set(i)%atomic_kind => atomic_kind
1139 1751 : particle_set(i)%r = pbc(center(1:3, i - natom), qs_loc_env%cell)
1140 : END DO
1141 : ! Dump the structure
1142 19 : CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
1143 :
1144 : ! Header file
1145 0 : SELECT CASE (outformat)
1146 : CASE (dump_dcd, dump_dcd_aligned_cell)
1147 0 : IF (init_traj) THEN
1148 : !Lets write the header for the coordinate dcd
1149 : ! Note (TL) : even the new DCD format is unfortunately too poor
1150 : ! for our capabilities.. for example here the printing
1151 : ! of the geometry could be nested inside several iteration
1152 : ! levels.. this cannot be exactly reproduced with DCD.
1153 : ! Just as a compromise let's pick-up the value of the MD iteration
1154 : ! level. In any case this is not any sensible information for the standard.
1155 0 : iskip = section_get_ival(print_key, "EACH%MD")
1156 0 : WRITE (iunit) "CORD", 0, -1, iskip, &
1157 0 : 0, 0, 0, 0, 0, 0, REAL(0, KIND=sp), 1, 0, 0, 0, 0, 0, 0, 0, 0, 24
1158 0 : remark1 = "REMARK FILETYPE CORD DCD GENERATED BY CP2K"
1159 0 : remark2 = "REMARK Support new DCD format with cell information"
1160 0 : WRITE (iunit) 2, remark1, remark2
1161 0 : WRITE (iunit) ntot
1162 0 : CALL m_flush(iunit)
1163 : END IF
1164 : CASE (dump_xmol)
1165 19 : iter = cp_iter_string(logger%iter_info)
1166 19 : WRITE (UNIT=title, FMT="(A)") " Particles+Wannier centers. Iteration:"//TRIM(iter)
1167 : CASE DEFAULT
1168 19 : title = ""
1169 : END SELECT
1170 : CALL write_particle_coordinates(particle_set, iunit, outformat, "POS", title, qs_loc_env%cell, &
1171 19 : unit_conv=unit_conv)
1172 19 : CALL m_flush(iunit)
1173 19 : CALL deallocate_particle_set(particle_set)
1174 19 : CALL deallocate_atomic_kind_set(atomic_kind_set)
1175 38 : END SUBROUTINE print_wannier_traj
1176 :
1177 : ! **************************************************************************************************
1178 : !> \brief Compute the second moments of the centers using the local (non-periodic) pos operators
1179 : !> \param qs_env ...
1180 : !> \param qs_loc_env ...
1181 : !> \param print_loc_section ...
1182 : !> \param ispin ...
1183 : !> \par History
1184 : !> 07.2020 created [H. Elgabarty]
1185 : !> \author Hossam Elgabarty
1186 : ! **************************************************************************************************
1187 0 : SUBROUTINE centers_second_moments_loc(qs_env, qs_loc_env, print_loc_section, ispin)
1188 :
1189 : ! I might not actually need the qs_env
1190 : TYPE(qs_environment_type), POINTER :: qs_env
1191 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1192 : TYPE(section_vals_type), POINTER :: print_loc_section
1193 : INTEGER, INTENT(IN) :: ispin
1194 :
1195 : INTEGER, PARAMETER :: norder = 2
1196 : LOGICAL, PARAMETER :: debug_this_subroutine = .FALSE.
1197 :
1198 : CHARACTER(len=default_path_length) :: file_tmp
1199 : INTEGER :: imoment, istate, ncol_global, nm, &
1200 : nmoloc, nrow_global, output_unit, &
1201 : output_unit_sm
1202 0 : REAL(dp), DIMENSION(:, :), POINTER :: centers
1203 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: moment_set
1204 : REAL(KIND=dp), DIMENSION(3) :: rcc
1205 : REAL(KIND=dp), DIMENSION(6) :: cov
1206 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1207 : TYPE(cp_fm_type) :: momv, mvector, omvector
1208 0 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
1209 : TYPE(cp_fm_type), POINTER :: mo_localized
1210 : TYPE(cp_logger_type), POINTER :: logger
1211 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
1212 : TYPE(mp_para_env_type), POINTER :: para_env
1213 :
1214 0 : logger => cp_get_default_logger()
1215 :
1216 : output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
1217 0 : extension=".locInfo")
1218 :
1219 0 : IF (output_unit > 0) THEN
1220 : WRITE (output_unit, '(/,T2,A)') &
1221 0 : '!-----------------------------------------------------------------------------!'
1222 0 : WRITE (output_unit, '(T12,A)') "Computing second moments of Wannier functions..."
1223 : WRITE (output_unit, '(T2,A)') &
1224 0 : '!-----------------------------------------------------------------------------!'
1225 : END IF
1226 :
1227 0 : file_tmp = "_r_loc_cov_matrix_"//TRIM(ADJUSTL(cp_to_string(ispin)))
1228 : output_unit_sm = cp_print_key_unit_nr(logger, print_loc_section, "WANNIER_SPREADS", &
1229 0 : middle_name=file_tmp, extension=".data")
1230 :
1231 0 : CALL get_qs_loc_env(qs_loc_env=qs_loc_env, moloc_coeff=moloc_coeff)
1232 0 : centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
1233 :
1234 0 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1235 :
1236 0 : nm = ncoset(norder) - 1
1237 :
1238 : !CALL get_qs_env(qs_env, qs_kind_set=qs_kind_set, cell=cell)
1239 : !particle_set => qs_loc_env%particle_set
1240 0 : para_env => qs_loc_env%para_env
1241 :
1242 0 : CALL cp_fm_get_info(moloc_coeff(ispin), ncol_global=nmoloc)
1243 0 : ALLOCATE (moment_set(nm, nmoloc))
1244 0 : moment_set = 0.0_dp
1245 :
1246 0 : mo_localized => moloc_coeff(ispin)
1247 :
1248 0 : DO istate = 1, nmoloc
1249 0 : rcc = centers(1:3, istate)
1250 : !rcc = 0.0_dp
1251 :
1252 0 : ALLOCATE (moments(nm))
1253 0 : DO imoment = 1, nm
1254 0 : ALLOCATE (moments(imoment)%matrix)
1255 0 : CALL dbcsr_copy(moments(imoment)%matrix, matrix_s(ispin)%matrix, 'MOM MAT')
1256 0 : CALL dbcsr_set(moments(imoment)%matrix, 0.0_dp)
1257 : END DO
1258 : !
1259 :
1260 0 : CALL build_local_moment_matrix(qs_env, moments, norder, rcc)
1261 :
1262 0 : CALL cp_fm_get_info(mo_localized, ncol_global=ncol_global, nrow_global=nrow_global)
1263 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_global, ncol_global=1, &
1264 : para_env=mo_localized%matrix_struct%para_env, &
1265 0 : context=mo_localized%matrix_struct%context)
1266 0 : CALL cp_fm_create(mvector, fm_struct, name="mvector")
1267 0 : CALL cp_fm_create(omvector, fm_struct, name="omvector")
1268 0 : CALL cp_fm_create(momv, fm_struct, name="omvector")
1269 0 : CALL cp_fm_struct_release(fm_struct)
1270 :
1271 : !cp_fm_to_fm_columns(msource, mtarget, ncol, source_start,target_start)
1272 0 : CALL cp_fm_to_fm(mo_localized, mvector, 1, istate, 1)
1273 :
1274 0 : DO imoment = 1, nm
1275 0 : CALL cp_dbcsr_sm_fm_multiply(moments(imoment)%matrix, mvector, omvector, 1)
1276 0 : CALL cp_fm_schur_product(mvector, omvector, momv)
1277 : !moment_set(imoment, istate) = moment_set(imoment, istate) + SUM(momv%local_data)
1278 0 : moment_set(imoment, istate) = SUM(momv%local_data)
1279 : END DO
1280 : !
1281 0 : CALL cp_fm_release(mvector)
1282 0 : CALL cp_fm_release(omvector)
1283 0 : CALL cp_fm_release(momv)
1284 :
1285 0 : DO imoment = 1, nm
1286 0 : CALL dbcsr_deallocate_matrix(moments(imoment)%matrix)
1287 : END DO
1288 0 : DEALLOCATE (moments)
1289 : END DO
1290 :
1291 0 : CALL para_env%sum(moment_set)
1292 :
1293 0 : IF (output_unit_sm > 0) THEN
1294 0 : WRITE (UNIT=output_unit_sm, FMT='(A,T13,A,I1)') "#", &
1295 0 : "SECOND MOMENTS OF WANNIER CENTERS FOR SPIN ", ispin
1296 0 : WRITE (UNIT=output_unit_sm, FMT='(A,T29,A2,5(14x,A2))') "#", "XX", "XY", "XZ", "YY", "YZ", "ZZ"
1297 0 : DO istate = 1, nmoloc
1298 : cov = 0.0_dp
1299 0 : cov(1) = moment_set(4, istate) - moment_set(1, istate)*moment_set(1, istate)
1300 0 : cov(2) = moment_set(5, istate) - moment_set(1, istate)*moment_set(2, istate)
1301 0 : cov(3) = moment_set(6, istate) - moment_set(1, istate)*moment_set(3, istate)
1302 0 : cov(4) = moment_set(7, istate) - moment_set(2, istate)*moment_set(2, istate)
1303 0 : cov(5) = moment_set(8, istate) - moment_set(2, istate)*moment_set(3, istate)
1304 0 : cov(6) = moment_set(9, istate) - moment_set(3, istate)*moment_set(3, istate)
1305 0 : WRITE (UNIT=output_unit_sm, FMT='(T4,A,I6,6(T20,6F16.8))') "Center:", istate, cov(1:6)
1306 0 : IF (debug_this_subroutine) THEN
1307 : WRITE (UNIT=output_unit_sm, FMT='(T4,A,I5,9(T20,9F12.6))') "Center:", istate, moment_set(1:, istate)
1308 : END IF
1309 : END DO
1310 : END IF
1311 0 : CALL cp_print_key_finished_output(output_unit_sm, logger, print_loc_section, "WANNIER_SPREADS")
1312 0 : CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
1313 :
1314 0 : DEALLOCATE (moment_set)
1315 :
1316 0 : END SUBROUTINE centers_second_moments_loc
1317 :
1318 : ! **************************************************************************************************
1319 : !> \brief Compute the second moments of the centers using a periodic quadrupole operator
1320 : !> \brief c.f. Wheeler et al. PRB 100 245135 2019, and Kang et al. PRB 100 245134 2019
1321 : !> \brief The calculation is based on a Taylor expansion of the exp(i k_i r_i r_j k_j) operator to
1322 : !> \brief to first order in the cosine and the sine parts
1323 : !> \param qs_env ...
1324 : !> \param qs_loc_env ...
1325 : !> \param print_loc_section ...
1326 : !> \param ispin ...
1327 : !> \par History
1328 : !> 08.2020 created [H. Elgabarty]
1329 : !> \author Hossam Elgabarty
1330 : ! **************************************************************************************************
1331 0 : SUBROUTINE centers_second_moments_berry(qs_env, qs_loc_env, print_loc_section, ispin)
1332 :
1333 : TYPE(qs_environment_type), POINTER :: qs_env
1334 : TYPE(qs_loc_env_type), POINTER :: qs_loc_env
1335 : TYPE(section_vals_type), POINTER :: print_loc_section
1336 : INTEGER, INTENT(IN) :: ispin
1337 :
1338 : INTEGER, PARAMETER :: taylor_order = 1
1339 : LOGICAL, PARAMETER :: debug_this_subroutine = .FALSE.
1340 :
1341 : CHARACTER(len=default_path_length) :: file_tmp
1342 : COMPLEX(KIND=dp) :: z
1343 : INTEGER :: icov, imoment, istate, ncol_global, nm, &
1344 : nmoloc, norder, nrow_global, &
1345 : output_unit, output_unit_sm
1346 0 : REAL(dp), DIMENSION(:, :), POINTER :: centers
1347 : REAL(KIND=dp) :: imagpart, Lx, Ly, Lz, realpart
1348 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: moment_set
1349 : REAL(KIND=dp), DIMENSION(3) :: rcc
1350 : REAL(KIND=dp), DIMENSION(6) :: cov
1351 : TYPE(cell_type), POINTER :: cell
1352 : TYPE(cp_fm_struct_type), POINTER :: fm_struct
1353 : TYPE(cp_fm_type) :: momv, mvector, omvector
1354 0 : TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
1355 : TYPE(cp_fm_type), POINTER :: mo_localized
1356 : TYPE(cp_logger_type), POINTER :: logger
1357 0 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
1358 : TYPE(mp_para_env_type), POINTER :: para_env
1359 :
1360 : ! two pointers, one to each spin channel's coeff. matrix (nao x nmoloc)
1361 :
1362 0 : logger => cp_get_default_logger()
1363 :
1364 : output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
1365 0 : extension=".locInfo")
1366 :
1367 0 : IF (output_unit > 0) THEN
1368 : WRITE (output_unit, '(/,T2,A)') &
1369 0 : '!-----------------------------------------------------------------------------!'
1370 0 : WRITE (output_unit, '(T12,A)') "Computing second moments of Wannier functions..."
1371 : WRITE (output_unit, '(T2,A)') &
1372 0 : '!-----------------------------------------------------------------------------!'
1373 : END IF
1374 :
1375 0 : file_tmp = "_r_berry_cov_matrix_"//TRIM(ADJUSTL(cp_to_string(ispin)))
1376 : output_unit_sm = cp_print_key_unit_nr(logger, print_loc_section, "WANNIER_SPREADS", &
1377 0 : middle_name=file_tmp, extension=".data")
1378 :
1379 0 : norder = 2*(2*taylor_order + 1)
1380 0 : nm = (6 + 11*norder + 6*norder**2 + norder**3)/6 - 1
1381 :
1382 0 : NULLIFY (cell)
1383 0 : CALL get_qs_loc_env(qs_loc_env=qs_loc_env, moloc_coeff=moloc_coeff, cell=cell)
1384 0 : Lx = cell%hmat(1, 1)
1385 0 : Ly = cell%hmat(2, 2)
1386 0 : Lz = cell%hmat(3, 3)
1387 :
1388 0 : centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
1389 :
1390 0 : para_env => qs_loc_env%para_env
1391 :
1392 0 : CALL get_qs_env(qs_env, matrix_s=matrix_s)
1393 :
1394 0 : CALL cp_fm_get_info(moloc_coeff(ispin), ncol_global=nmoloc)
1395 :
1396 0 : ALLOCATE (moment_set(nm, nmoloc))
1397 0 : moment_set = 0.0_dp
1398 :
1399 0 : mo_localized => moloc_coeff(ispin)
1400 :
1401 0 : DO istate = 1, nmoloc
1402 0 : rcc = centers(1:3, istate)
1403 : !rcc = 0.0_dp
1404 :
1405 0 : ALLOCATE (moments(nm))
1406 0 : DO imoment = 1, nm
1407 0 : ALLOCATE (moments(imoment)%matrix)
1408 0 : CALL dbcsr_copy(moments(imoment)%matrix, matrix_s(ispin)%matrix, 'MOM MAT')
1409 0 : CALL dbcsr_set(moments(imoment)%matrix, 0.0_dp)
1410 : END DO
1411 : !
1412 :
1413 0 : CALL build_local_moment_matrix(qs_env, moments, norder, rcc)
1414 :
1415 0 : CALL cp_fm_get_info(mo_localized, ncol_global=ncol_global, nrow_global=nrow_global)
1416 : CALL cp_fm_struct_create(fm_struct, nrow_global=nrow_global, ncol_global=1, &
1417 : para_env=mo_localized%matrix_struct%para_env, &
1418 0 : context=mo_localized%matrix_struct%context)
1419 0 : CALL cp_fm_create(mvector, fm_struct, name="mvector")
1420 0 : CALL cp_fm_create(omvector, fm_struct, name="omvector")
1421 0 : CALL cp_fm_create(momv, fm_struct, name="omvector")
1422 0 : CALL cp_fm_struct_release(fm_struct)
1423 :
1424 : ! cp_fm_to_fm_columns(msource, mtarget, ncol, source_start,target_start)
1425 0 : CALL cp_fm_to_fm(mo_localized, mvector, 1, istate, 1)
1426 :
1427 0 : DO imoment = 1, nm
1428 0 : CALL cp_dbcsr_sm_fm_multiply(moments(imoment)%matrix, mvector, omvector, 1)
1429 0 : CALL cp_fm_schur_product(mvector, omvector, momv)
1430 0 : moment_set(imoment, istate) = SUM(momv%local_data)
1431 : END DO
1432 : !
1433 0 : CALL cp_fm_release(mvector)
1434 0 : CALL cp_fm_release(omvector)
1435 0 : CALL cp_fm_release(momv)
1436 :
1437 0 : DO imoment = 1, nm
1438 0 : CALL dbcsr_deallocate_matrix(moments(imoment)%matrix)
1439 : END DO
1440 0 : DEALLOCATE (moments)
1441 : END DO
1442 :
1443 0 : CALL para_env%sum(moment_set)
1444 :
1445 0 : IF (output_unit_sm > 0) THEN
1446 0 : WRITE (UNIT=output_unit_sm, FMT='(A,T13,A,I1)') "#", &
1447 0 : "SECOND MOMENTS OF WANNIER CENTERS FOR SPIN ", ispin
1448 0 : WRITE (UNIT=output_unit_sm, FMT='(A,T29,A2,5(14x,A2))') "#", "XX", "XY", "XZ", "YY", "YZ", "ZZ"
1449 0 : DO istate = 1, nmoloc
1450 0 : cov = 0.0_dp
1451 0 : DO icov = 1, 6
1452 0 : realpart = 0.0_dp
1453 0 : imagpart = 0.0_dp
1454 0 : z = CMPLX(realpart, imagpart, dp)
1455 0 : SELECT CASE (icov)
1456 :
1457 : !! XX
1458 : CASE (1)
1459 : realpart = 1.0 - 0.5*twopi*twopi/Lx/Lx/Lx/Lx &
1460 0 : *moment_set(20, istate)
1461 : imagpart = twopi/Lx/Lx*moment_set(4, istate) &
1462 0 : - twopi*twopi*twopi/Lx/Lx/Lx/Lx/Lx/Lx/6.0*moment_set(56, istate)
1463 :
1464 : ! third order
1465 : !realpart = realpart + twopi*twopi*twopi*twopi/Lx/Lx/Lx/Lx/Lx/Lx/Lx/Lx/24.0 * moment_set(120, istate)
1466 : !imagpart = imagpart + twopi*twopi*twopi*twopi*twopi/Lx/Lx/Lx/Lx/Lx/Lx/Lx/Lx/Lx/Lx/120 &
1467 : ! * moment_set(220, istate)
1468 :
1469 0 : z = CMPLX(realpart, imagpart, dp)
1470 0 : cov(1) = Lx*Lx/twopi*AIMAG(LOG(z)) - Lx*Lx/twopi/twopi*moment_set(1, istate)*moment_set(1, istate)
1471 :
1472 : !! XY
1473 : CASE (2)
1474 : realpart = 1.0 - 0.5*twopi*twopi/Lx/Ly/Lx/Ly &
1475 0 : *moment_set(23, istate)
1476 : imagpart = twopi/Lx/Ly*moment_set(5, istate) &
1477 0 : - twopi*twopi*twopi/Lx/Lx/Lx/Ly/Ly/Ly/6.0*moment_set(62, istate)
1478 :
1479 : ! third order
1480 : !realpart = realpart + twopi*twopi*twopi*twopi/Lx/Lx/Lx/Lx/Ly/Ly/Ly/Ly/24.0 * moment_set(130, istate)
1481 : !imagpart = imagpart + twopi*twopi*twopi*twopi*twopi/Lx/Lx/Lx/Lx/Lx/Ly/Ly/Ly/Ly/Ly/120 &
1482 : ! * moment_set(235, istate)
1483 :
1484 0 : z = CMPLX(realpart, imagpart, dp)
1485 0 : cov(2) = Lx*Ly/twopi*AIMAG(LOG(z)) - Lx*Ly/twopi/twopi*moment_set(1, istate)*moment_set(2, istate)
1486 :
1487 : !! XZ
1488 : CASE (3)
1489 : realpart = 1.0 - 0.5*twopi*twopi/Lx/Lz/Lx/Lz &
1490 0 : *moment_set(25, istate)
1491 : imagpart = twopi/Lx/Lz*moment_set(6, istate) &
1492 0 : - twopi*twopi*twopi/Lx/Lx/Lx/Lz/Lz/Lz/6.0*moment_set(65, istate)
1493 :
1494 : ! third order
1495 : !realpart = realpart + twopi*twopi*twopi*twopi/Lx/Lx/Lx/Lx/Lz/Lz/Lz/Lz/24.0 * moment_set(134, istate)
1496 : !imagpart = imagpart + twopi*twopi*twopi*twopi*twopi/Lx/Lx/Lx/Lx/Lx/Ly/Ly/Ly/Ly/Ly/120 &
1497 : ! * moment_set(240, istate)
1498 :
1499 0 : z = CMPLX(realpart, imagpart, dp)
1500 0 : cov(3) = Lx*Lz/twopi*AIMAG(LOG(z)) - Lx*Lz/twopi/twopi*moment_set(1, istate)*moment_set(3, istate)
1501 :
1502 : !! YY
1503 : CASE (4)
1504 : realpart = 1.0 - 0.5*twopi*twopi/Ly/Ly/Ly/Ly &
1505 0 : *moment_set(30, istate)
1506 : imagpart = twopi/Ly/Ly*moment_set(7, istate) &
1507 0 : - twopi*twopi*twopi/Ly/Ly/Ly/Ly/Ly/Ly/6.0*moment_set(77, istate)
1508 :
1509 : ! third order
1510 : !realpart = realpart + twopi*twopi*twopi*twopi/Ly/Ly/Ly/Ly/Ly/Ly/Ly/Ly/24.0 * moment_set(156, istate)
1511 : !imagpart = imagpart + twopi*twopi*twopi*twopi*twopi/Ly/Ly/Ly/Ly/Ly/Ly/Ly/Ly/Ly/Ly/120 &
1512 : ! * moment_set(275, istate)
1513 :
1514 0 : z = CMPLX(realpart, imagpart, dp)
1515 0 : cov(4) = Ly*Ly/twopi*AIMAG(LOG(z)) - Ly*Ly/twopi/twopi*moment_set(2, istate)*moment_set(2, istate)
1516 :
1517 : !! YZ
1518 : CASE (5)
1519 : realpart = 1.0 - 0.5*twopi*twopi/Ly/Lz/Ly/Lz &
1520 0 : *moment_set(32, istate)
1521 : imagpart = twopi/Ly/Lz*moment_set(8, istate) &
1522 0 : - twopi*twopi*twopi/Ly/Ly/Ly/Lz/Lz/Lz/6.0*moment_set(80, istate)
1523 :
1524 : ! third order
1525 : !realpart = realpart + twopi*twopi*twopi*twopi/Lz/Lz/Lz/Lz/Ly/Ly/Ly/Ly/24.0 * moment_set(160, istate)
1526 : !imagpart = imagpart + twopi*twopi*twopi*twopi*twopi/Lz/Lz/Lz/Lz/Lz/Ly/Ly/Ly/Ly/Ly/120 &
1527 : ! * moment_set(280, istate)
1528 :
1529 0 : z = CMPLX(realpart, imagpart, dp)
1530 0 : cov(5) = Ly*Lz/twopi*AIMAG(LOG(z)) - Ly*Lz/twopi/twopi*moment_set(2, istate)*moment_set(3, istate)
1531 :
1532 : !! ZZ
1533 : CASE (6)
1534 : realpart = 1.0 - 0.5*twopi*twopi/Lz/Lz/Lz/Lz &
1535 0 : *moment_set(34, istate)
1536 : imagpart = twopi/Lz/Lz*moment_set(9, istate) &
1537 0 : - twopi*twopi*twopi/Lz/Lz/Lz/Lz/Lz/Lz/6.0*moment_set(83, istate)
1538 :
1539 : ! third order
1540 : !realpart = realpart + twopi*twopi*twopi*twopi/Lz/Lz/Lz/Lz/Lz/Lz/Lz/Lz/24.0 * moment_set(164, istate)
1541 : !imagpart = imagpart + twopi*twopi*twopi*twopi*twopi/Lz/Lz/Lz/Lz/Lz/Lz/Lz/Lz/Lz/Lz/120 &
1542 : ! * moment_set(285, istate)
1543 :
1544 0 : z = CMPLX(realpart, imagpart, dp)
1545 0 : cov(6) = Lz*Lz/twopi*AIMAG(LOG(z)) - Lz*Lz/twopi/twopi*moment_set(3, istate)*moment_set(3, istate)
1546 :
1547 : END SELECT
1548 : END DO
1549 0 : WRITE (UNIT=output_unit_sm, FMT='(T4,A,I6,6(T20,6F16.8))') "Center:", istate, cov(1:6)
1550 0 : IF (debug_this_subroutine) THEN
1551 : WRITE (UNIT=output_unit_sm, FMT='(T4,A,I5,9(T20,9F12.6))') "Center:", istate, moment_set(1:, istate)
1552 : END IF
1553 : END DO
1554 : END IF
1555 0 : CALL cp_print_key_finished_output(output_unit_sm, logger, print_loc_section, "WANNIER_SPREADS")
1556 0 : CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
1557 :
1558 0 : DEALLOCATE (moment_set)
1559 :
1560 0 : END SUBROUTINE centers_second_moments_berry
1561 :
1562 : END MODULE qs_loc_methods
|
