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