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