(git:0de0cc2)
qs_loc_methods.F
Go to the documentation of this file.
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 ! **************************************************************************************************
22  USE atomic_kind_types, ONLY: atomic_kind_type,&
25  USE cell_types, ONLY: cell_type,&
26  pbc
27  USE cp_control_types, ONLY: dft_control_type
33  USE cp_fm_pool_types, ONLY: cp_fm_pool_p_type,&
37  cp_fm_struct_type
38  USE cp_fm_types, ONLY: cp_fm_create,&
42  cp_fm_release,&
44  cp_fm_to_fm,&
45  cp_fm_type
47  cp_logger_type,&
48  cp_to_string
50  cp_p_file,&
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: &
66  section_vals_type,&
68  USE kinds, ONLY: default_path_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
77  USE orbital_pointers, ONLY: ncoset
78  USE parallel_gemm_api, ONLY: parallel_gemm
79  USE particle_list_types, ONLY: particle_list_type
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
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
99  direct_mini,&
103  scdm_qrfact,&
104  zij_matrix
105  USE qs_matrix_pools, ONLY: mpools_get
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 ***
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  SUBROUTINE optimize_loc_berry(method, qs_loc_env, vectors, op_sm_set, &
159  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  CALL timeset(routinen, handle)
188  logger => cp_get_default_logger()
189  output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
190  extension=".locInfo")
191 
192  ! get rows and cols of the input
193  CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc)
194 
195  CALL zij_matrix(vectors, op_sm_set, zij_fm_set)
196 
197  max_iter = qs_loc_env%localized_wfn_control%max_iter
198  max_crazy_angle = qs_loc_env%localized_wfn_control%max_crazy_angle
199  crazy_use_diag = qs_loc_env%localized_wfn_control%crazy_use_diag
200  crazy_scale = qs_loc_env%localized_wfn_control%crazy_scale
201  eps_localization = qs_loc_env%localized_wfn_control%eps_localization
202  out_each = qs_loc_env%localized_wfn_control%out_each
203  target_time = qs_loc_env%target_time
204  start_time = qs_loc_env%start_time
205  do_jacobi_refinement = qs_loc_env%localized_wfn_control%jacobi_refinement
206  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  ispin, print_loc_section, only_initial_out=.true.)
210 
211  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  restricted=restricted)
217  CASE (do_loc_gapo)
218  IF (my_do_mixed) THEN
219  IF (nextra > 0) THEN
220  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  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  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  qs_loc_env%localized_wfn_control%do_cg_po)
235  END IF
236  ELSE
237  cpabort("GAPO works only with STATES MIXED")
238  END IF
239  CASE (do_loc_scdm)
240  ! Decomposition
241  CALL scdm_qrfact(vectors)
242  ! Calculation of Zij
243  CALL zij_matrix(vectors, op_sm_set, zij_fm_set)
244  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  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  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  eps_localization=eps_localization, iterations=sweeps, converged=converged)
257  ! Possibly fallback to jacobi if the crazy rotation fails
258  IF (.NOT. converged) THEN
259  IF (qs_loc_env%localized_wfn_control%jacobi_fallback) THEN
260  IF (output_unit > 0) WRITE (output_unit, "(T4,A,I6,/,T4,A)") &
261  " Crazy Wannier localization not converged after ", sweeps, &
262  " 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  restricted=restricted)
267  ELSE
268  IF (output_unit > 0) WRITE (output_unit, "(T4,A,I6,/,T4,A)") &
269  " Crazy Wannier localization not converged after ", sweeps, &
270  " 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  eps_localization=eps_localization, iterations=sweeps)
276  CASE (do_loc_l1_norm_sd)
277  IF (.NOT. cell%orthorhombic) THEN
278  cpabort("Non-orthorhombic cell with the selected method NYI")
279  ELSE
280  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  CALL zij_matrix(vectors, op_sm_set, zij_fm_set)
283  END IF
284  CASE (do_loc_none)
285  IF (output_unit > 0) THEN
286  WRITE (output_unit, '(T4,A,I6,A)') " No MOS localization applied "
287  END IF
288  CASE DEFAULT
289  cpabort("Unknown localization method")
290  END SELECT
291  IF (output_unit > 0) THEN
292  IF (sweeps <= max_iter) WRITE (output_unit, '(T4,A,I3,A,I6,A)') " Localization for spin ", ispin, &
293  " converged in ", sweeps, " iterations"
294  END IF
295 
296  CALL centers_spreads_berry(qs_loc_env, zij_fm_set, nmoloc, cell, weights, &
297  ispin, print_loc_section)
298 
299  CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
300 
301  CALL timestop(handle)
302 
303  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  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  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  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s
338  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
339  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
340 
341  CALL timeset(routinen, handle)
342  logger => cp_get_default_logger()
343  output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
344  extension=".locInfo")
345 
346  NULLIFY (particle_set)
347  ! get rows and cols of the input
348  CALL cp_fm_get_info(vectors, nrow_global=nao, ncol_global=nmoloc)
349 
350  ! replicate the input kind of matrix
351  CALL cp_fm_create(opvec, vectors%matrix_struct)
352  CALL cp_fm_set_all(opvec, 0.0_dp)
353 
354  CALL get_qs_env(qs_env=qs_env, matrix_s=matrix_s, &
355  particle_set=particle_set, qs_kind_set=qs_kind_set)
356 
357  natom = SIZE(particle_set, 1)
358  ALLOCATE (first_sgf(natom))
359  ALLOCATE (last_sgf(natom))
360  ALLOCATE (nsgf(natom))
361 
362  ! construction of
363  CALL get_particle_set(particle_set, qs_kind_set, &
364  first_sgf=first_sgf, last_sgf=last_sgf, nsgf=nsgf)
365 
366  ! Copy the overlap sparse matrix in a full matrix
367  CALL mpools_get(qs_env%mpools, ao_ao_fm_pools=ao_ao_fm_pools)
368  ALLOCATE (ov_fm)
369  CALL fm_pool_create_fm(ao_ao_fm_pools(1)%pool, ov_fm, name=" ")
370  CALL copy_dbcsr_to_fm(matrix_s(1)%matrix, ov_fm)
371 
372  ! Compute zij here
373  DO iatom = 1, natom
374  CALL cp_fm_set_all(zij_fm_set(iatom, 1), 0.0_dp)
375  CALL cp_fm_get_info(zij_fm_set(iatom, 1), ncol_global=ldz)
376  isgf = first_sgf(iatom)
377  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  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  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  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  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  SELECT CASE (method)
398  CASE (do_loc_jacobi)
399  CALL jacobi_rotation_pipek(zij_fm_set, vectors, sweeps)
400  CASE (do_loc_gapo)
401  cpabort("GAPO and Pipek not implemented.")
402  CASE (do_loc_crazy)
403  cpabort("Crazy and Pipek not implemented.")
404  CASE (do_loc_l1_norm_sd)
405  cpabort("L1 norm and Pipek not implemented.")
406  CASE (do_loc_direct)
407  cpabort("Direct and Pipek not implemented.")
408  CASE (do_loc_none)
409  IF (output_unit > 0) WRITE (output_unit, '(A,I6,A)') " No MOS localization applied "
410  CASE DEFAULT
411  cpabort("Unknown localization method")
412  END SELECT
413 
414  IF (output_unit > 0) WRITE (output_unit, '(A,I3,A,I6,A)') " Localization for spin ", ispin, &
415  " converged in ", sweeps, " iterations"
416 
417  CALL centers_spreads_pipek(qs_loc_env, zij_fm_set, particle_set, ispin, &
418  print_loc_section)
419 
420  DEALLOCATE (first_sgf, last_sgf, nsgf)
421 
422  CALL cp_fm_release(opvec)
423  CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
424 
425  CALL timestop(handle)
426 
427  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  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  CALL cp_fm_create(rmat, zij_fm_set(1, 1)%matrix_struct)
451  CALL cp_fm_set_all(rmat, 0.0_dp, 1.0_dp)
452 
453  CALL cp_fm_get_info(rmat, nrow_global=nstate)
454  tolerance = 1.0e10_dp
455  sweeps = 0
456  natom = SIZE(zij_fm_set, 1)
457  ! do jacobi sweeps until converged
458  DO WHILE (tolerance >= 1.0e-4_dp)
459  sweeps = sweeps + 1
460  DO istate = 1, nstate
461  DO jstate = istate + 1, nstate
462  aij = 0.0_dp
463  bij = 0.0_dp
464  DO iatom = 1, natom
465  CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, mii)
466  CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, jstate, mij)
467  CALL cp_fm_get_element(zij_fm_set(iatom, 1), jstate, jstate, mjj)
468  aij = aij + mij*(mii - mjj)
469  bij = bij + mij*mij - 0.25_dp*(mii - mjj)*(mii - mjj)
470  END DO
471  IF (abs(bij) > 1.e-10_dp) THEN
472  ratio = -aij/bij
473  theta = 0.25_dp*atan(ratio)
474  ELSE
475  bij = 0.0_dp
476  theta = 0.0_dp
477  END IF
478  ! Check max or min
479  ! To minimize the spread
480  IF (theta > pi*0.5_dp) THEN
481  theta = theta - pi*0.25_dp
482  ELSE IF (theta < -pi*0.5_dp) THEN
483  theta = theta + pi*0.25_dp
484  END IF
485 
486  ct = cos(theta)
487  st = sin(theta)
488 
489  CALL cp_fm_rot_cols(rmat, istate, jstate, ct, st)
490 
491  DO iatom = 1, natom
492  CALL cp_fm_rot_cols(zij_fm_set(iatom, 1), istate, jstate, ct, st)
493  CALL cp_fm_rot_rows(zij_fm_set(iatom, 1), istate, jstate, ct, st)
494  END DO
495  END DO
496  END DO
497  CALL check_tolerance_real(zij_fm_set, tolerance)
498  END DO
499 
500  CALL rotate_orbitals(rmat, vectors)
501  CALL cp_fm_release(rmat)
502 
503  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  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  INTEGER, DIMENSION(:), POINTER :: col_indices, row_indices
521  REAL(dp) :: grad_ij, zii, zij, zjj
522  REAL(dp), DIMENSION(:, :), POINTER :: diag
523  TYPE(cp_fm_type) :: force
524 
525  CALL cp_fm_create(force, zij_fm_set(1, 1)%matrix_struct)
526  CALL cp_fm_set_all(force, 0._dp)
527 
528  NULLIFY (diag, col_indices, row_indices)
529  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  row_indices=row_indices, col_indices=col_indices)
533  ALLOCATE (diag(nrow_global, natom))
534 
535  DO iatom = 1, natom
536  DO istate = 1, nrow_global
537  CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, diag(istate, iatom))
538  END DO
539  END DO
540 
541  DO istate = 1, nrow_local
542  DO jstate = 1, ncol_local
543  grad_ij = 0.0_dp
544  DO iatom = 1, natom
545  zii = diag(row_indices(istate), iatom)
546  zjj = diag(col_indices(jstate), iatom)
547  zij = zij_fm_set(iatom, 1)%local_data(istate, jstate)
548  grad_ij = grad_ij + 4.0_dp*zij*(zjj - zii)
549  END DO
550  force%local_data(istate, jstate) = grad_ij
551  END DO
552  END DO
553 
554  DEALLOCATE (diag)
555 
556  CALL cp_fm_maxabsval(force, tolerance)
557  CALL cp_fm_release(force)
558 
559  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  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  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  NULLIFY (centers, logger, print_key)
599  logger => cp_get_default_logger()
600  my_only_init = .false.
601  IF (PRESENT(only_initial_out)) my_only_init = only_initial_out
602 
603  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  extension=".locInfo")
606  unit_out_s = cp_print_key_unit_nr(logger, print_loc_section, "WANNIER_SPREADS", &
607  middle_name=file_tmp, extension=".data")
608  iter = cp_iter_string(logger%iter_info)
609  IF (unit_out_s > 0 .AND. .NOT. my_only_init) WRITE (unit_out_s, '(i6,/,A)') nmoloc, trim(iter)
610 
611  CALL cp_fm_get_info(zij(1, 1), nrow_global=nstates)
612  cpassert(nstates >= nmoloc)
613 
614  centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
615  cpassert(ASSOCIATED(centers))
616  cpassert(SIZE(centers, 2) == nmoloc)
617  sum_spread_i = 0.0_dp
618  sum_spread_ii = 0.0_dp
619  avg_spread_ii = 0.0_dp
620  DO istate = 1, nmoloc
621  c = 0.0_dp
622  c2 = 0.0_dp
623  spread_i = 0.0_dp
624  spread_ii = 0.0_dp
625  DO jdir = 1, SIZE(zij, 2)
626  CALL cp_fm_get_element(zij(1, jdir), istate, istate, realpart)
627  CALL cp_fm_get_element(zij(2, jdir), istate, istate, imagpart)
628  z = cmplx(realpart, imagpart, dp)
629  spread_i = spread_i - weights(jdir)* &
630  log(realpart*realpart + imagpart*imagpart)/twopi/twopi
631  spread_ii = spread_ii + weights(jdir)* &
632  (1.0_dp - (realpart*realpart + imagpart*imagpart))/twopi/twopi
633  IF (jdir < 4) THEN
634  DO idir = 1, 3
635  c(idir) = c(idir) + &
636  (cell%hmat(idir, jdir)/twopi)*aimag(log(z))
637  END DO
638  END IF
639  END DO
640  cpbc = pbc(c, cell)
641  centers(1:3, istate) = cpbc(1:3)
642  centers(4, istate) = spread_i
643  centers(5, istate) = spread_ii
644  sum_spread_i = sum_spread_i + centers(4, istate)
645  sum_spread_ii = sum_spread_ii + centers(5, istate)
646  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  avg_spread_ii = sum_spread_ii/real(nmoloc, kind=dp)
649 
650  ! Print of wannier centers
651  print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CENTERS")
652  IF (.NOT. my_only_init) CALL print_wannier_centers(qs_loc_env, print_key, centers, logger, ispin)
653 
654  IF (output_unit > 0) THEN
655  WRITE (output_unit, '(T4, A, 2x, 2A26,/,T23, A28)') " Spread Functional ", "sum_in -w_i ln(|z_in|^2)", &
656  "sum_in w_i(1-|z_in|^2)", "sum_in w_i(1-|z_in|^2)/n"
657  IF (my_only_init) THEN
658  WRITE (output_unit, '(T4,A,T38,2F20.10,/,T38,F20.10)') " Initial Spread (Berry) : ", &
659  sum_spread_i, sum_spread_ii, avg_spread_ii
660  ELSE
661  WRITE (output_unit, '(T4,A,T38,2F20.10,/,T38,F20.10)') " Total Spread (Berry) : ", &
662  sum_spread_i, sum_spread_ii, avg_spread_ii
663  END IF
664  END IF
665 
666  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  CALL cp_print_key_finished_output(unit_out_s, logger, print_loc_section, "WANNIER_SPREADS")
669  CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
670 
671  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  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  INTEGER, DIMENSION(:), POINTER :: atom_of_state
696  REAL(dp) :: r(3)
697  REAL(dp), ALLOCATABLE, DIMENSION(:) :: qii, ziimax
698  REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: diag
699  REAL(dp), DIMENSION(:, :), POINTER :: centers
700  TYPE(cp_logger_type), POINTER :: logger
701  TYPE(section_vals_type), POINTER :: print_key
702 
703  NULLIFY (centers, logger, print_key)
704  logger => cp_get_default_logger()
705 
706  CALL cp_fm_get_info(zij_fm_set(1, 1), nrow_global=nstate)
707  natom = SIZE(zij_fm_set, 1)
708 
709  centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
710  cpassert(ASSOCIATED(centers))
711  cpassert(SIZE(centers, 2) == nstate)
712 
713  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  middle_name=file_tmp, extension=".data", log_filename=.false.)
716  iter = cp_iter_string(logger%iter_info)
717  IF (unit_out_s > 0) WRITE (unit_out_s, '(i6,/,A)') trim(iter)
718 
719  ALLOCATE (atom_of_state(nstate))
720  atom_of_state = 0
721 
722  ALLOCATE (diag(nstate, natom))
723  diag = 0.0_dp
724 
725  DO iatom = 1, natom
726  DO istate = 1, nstate
727  CALL cp_fm_get_element(zij_fm_set(iatom, 1), istate, istate, diag(istate, iatom))
728  END DO
729  END DO
730 
731  ALLOCATE (qii(nstate), ziimax(nstate))
732  ziimax = 0.0_dp
733  qii = 0.0_dp
734 
735  DO iatom = 1, natom
736  DO istate = 1, nstate
737  qii(istate) = qii(istate) + diag(istate, iatom)*diag(istate, iatom)
738  IF (abs(diag(istate, iatom)) > ziimax(istate)) THEN
739  ziimax(istate) = abs(diag(istate, iatom))
740  atom_of_state(istate) = iatom
741  END IF
742  END DO
743  END DO
744 
745  DO istate = 1, nstate
746  iatom = atom_of_state(istate)
747  r(1:3) = particle_set(iatom)%r(1:3)
748  centers(1:3, istate) = r(1:3)
749  centers(4, istate) = 1.0_dp/qii(istate)
750  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  print_key => section_vals_get_subs_vals(print_loc_section, "WANNIER_CENTERS")
755  CALL print_wannier_centers(qs_loc_env, print_key, centers, logger, ispin)
756 
757  CALL cp_print_key_finished_output(unit_out_s, logger, print_loc_section, "WANNIER_SPREADS")
758 
759  DEALLOCATE (qii, ziimax, atom_of_state, diag)
760 
761  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  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  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  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  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
816  TYPE(qs_subsys_type), POINTER :: subsys
817 
818  CALL timeset(routinen, handle)
819  NULLIFY (auxbas_pw_pool, pw_env, logger)
820  logger => cp_get_default_logger()
821 
822  CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, subsys=subsys)
823  CALL qs_subsys_get(subsys, particles=particles)
824 
825  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
826 
827  CALL auxbas_pw_pool%create_pw(wf_r)
828  CALL auxbas_pw_pool%create_pw(wf_g)
829 
830  my_state0 = 0
831  IF (PRESENT(state0)) my_state0 = state0
832 
833  add_spin = .false.
834  my_ispin = 1
835  IF (PRESENT(ispin)) THEN
836  add_spin = .true.
837  my_ispin = ispin
838  END IF
839  add_idir = .false.
840  IF (PRESENT(idir)) THEN
841  add_idir = .true.
842  label = labels(idir)
843  END IF
844 
845  my_pos = "REWIND"
846  IF (PRESENT(file_position)) THEN
847  my_pos = file_position
848  END IF
849 
850  DO istate = 1, nstates
851  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  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  qs_kind_set, cell, dft_control, particle_set, pw_env)
857 
858  ! Formatting the middle part of the name
859  ivector = state_list(istate)
860  CALL xstring(root, ia, ie)
861  IF (add_idir) THEN
862  filename = root(ia:ie)//"_"//label//"_w"//trim(adjustl(cp_to_string(ivector)))
863  ELSE
864  filename = root(ia:ie)//"_w"//trim(adjustl(cp_to_string(ivector)))
865  END IF
866  IF (add_spin) THEN
867  file_tmp = filename
868  CALL xstring(file_tmp, ia, ie)
869  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  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  mpi_io=mpi_io)
877  IF (SIZE(centers, 1) == 6) THEN
878  WRITE (title, '(A7,I5.5,A2,I1.1,A1,6F10.4)') "WFN ", ivector, "_s", my_ispin, " ", &
879  centers(1:3, istate)*angstrom, centers(4:6, istate)*angstrom
880  ELSE
881  WRITE (title, '(A7,I5.5,A2,I1.1,A1,3F10.4)') "WFN ", ivector, "_s", my_ispin, " ", &
882  centers(1:3, istate)*angstrom
883  END IF
884  CALL cp_pw_to_cube(wf_r, unit_out_c, title, &
885  particles=particles, &
886  stride=section_get_ivals(print_key, "STRIDE"), mpi_io=mpi_io)
887  CALL cp_print_key_finished_output(unit_out_c, logger, print_key, "", mpi_io=mpi_io)
888  END DO
889 
890  CALL auxbas_pw_pool%give_back_pw(wf_r)
891  CALL auxbas_pw_pool%give_back_pw(wf_g)
892  CALL timestop(handle)
893  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  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  nstates = SIZE(center, 2)
917  my_form = "formatted"
918  my_ext = ".data"
919  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  IF (btest(cp_print_key_should_output(logger%iter_info, print_key, "/IONS+CENTERS") &
923  , cp_p_file)) THEN
924  CALL get_output_format(print_key, my_form=my_form, my_ext=my_ext)
925  END IF
926  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  log_filename=.false., on_file=.true., is_new_file=init_traj)
930  IF (iunit > 0) THEN
931  ! Gather units of measure for output (if available)
932  CALL section_vals_val_get(print_key, "UNIT", c_val=unit_str)
933  unit_conv = cp_unit_from_cp2k(1.0_dp, trim(unit_str))
934 
935  IF (btest(cp_print_key_should_output(logger%iter_info, print_key, "/IONS+CENTERS"), cp_p_file)) THEN
936  ! Different possible formats
937  CALL print_wannier_traj(qs_loc_env, print_key, center, iunit, init_traj, unit_conv)
938  ELSE
939  ! Default print format
940  iter = cp_iter_string(logger%iter_info)
941  WRITE (iunit, '(i6,/,a)') nstates, trim(iter)
942  DO l = 1, nstates
943  WRITE (iunit, '(A,4F16.8)') "X", unit_conv*center(1:4, l)
944  END DO
945  END IF
946  END IF
947  CALL cp_print_key_finished_output(iunit, logger, print_key, on_file=.true.)
948  END IF
949  END IF
950  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  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  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  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
975 
976  NULLIFY (particle_set, atomic_kind_set, atomic_kind, logger)
977  logger => cp_get_default_logger()
978  natom = SIZE(qs_loc_env%particle_set)
979  ntot = natom + SIZE(center, 2)
980  CALL allocate_particle_set(particle_set, ntot)
981  ALLOCATE (atomic_kind_set(1))
982  atomic_kind => atomic_kind_set(1)
983  CALL set_atomic_kind(atomic_kind=atomic_kind, kind_number=0, &
984  name="X", element_symbol="X", mass=0.0_dp)
985  ! Particles
986  DO i = 1, natom
987  particle_set(i)%atomic_kind => qs_loc_env%particle_set(i)%atomic_kind
988  particle_set(i)%r = pbc(qs_loc_env%particle_set(i)%r, qs_loc_env%cell)
989  END DO
990  ! Wannier Centers
991  DO i = natom + 1, ntot
992  particle_set(i)%atomic_kind => atomic_kind
993  particle_set(i)%r = pbc(center(1:3, i - natom), qs_loc_env%cell)
994  END DO
995  ! Dump the structure
996  CALL section_vals_val_get(print_key, "FORMAT", i_val=outformat)
997 
998  ! Header file
999  SELECT CASE (outformat)
1001  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  iskip = section_get_ival(print_key, "EACH%MD")
1010  WRITE (iunit) "CORD", 0, -1, iskip, &
1011  0, 0, 0, 0, 0, 0, real(0, kind=sp), 1, 0, 0, 0, 0, 0, 0, 0, 0, 24
1012  remark1 = "REMARK FILETYPE CORD DCD GENERATED BY CP2K"
1013  remark2 = "REMARK Support new DCD format with cell information"
1014  WRITE (iunit) 2, remark1, remark2
1015  WRITE (iunit) ntot
1016  CALL m_flush(iunit)
1017  END IF
1018  CASE (dump_xmol)
1019  iter = cp_iter_string(logger%iter_info)
1020  WRITE (unit=title, fmt="(A)") " Particles+Wannier centers. Iteration:"//trim(iter)
1021  CASE DEFAULT
1022  title = ""
1023  END SELECT
1024  CALL write_particle_coordinates(particle_set, iunit, outformat, "POS", title, qs_loc_env%cell, &
1025  unit_conv=unit_conv)
1026  CALL m_flush(iunit)
1027  CALL deallocate_particle_set(particle_set)
1028  CALL deallocate_atomic_kind_set(atomic_kind_set)
1029  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  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  REAL(dp), DIMENSION(:, :), POINTER :: centers
1057  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  TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
1063  TYPE(cp_fm_type), POINTER :: mo_localized
1064  TYPE(cp_logger_type), POINTER :: logger
1065  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_s, moments
1066  TYPE(mp_para_env_type), POINTER :: para_env
1067 
1068  logger => cp_get_default_logger()
1069 
1070  output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
1071  extension=".locInfo")
1072 
1073  IF (output_unit > 0) THEN
1074  WRITE (output_unit, '(/,T2,A)') &
1075  '!-----------------------------------------------------------------------------!'
1076  WRITE (output_unit, '(T12,A)') "Computing second moments of Wannier functions..."
1077  WRITE (output_unit, '(T2,A)') &
1078  '!-----------------------------------------------------------------------------!'
1079  END IF
1080 
1081  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  middle_name=file_tmp, extension=".data")
1084 
1085  CALL get_qs_loc_env(qs_loc_env=qs_loc_env, moloc_coeff=moloc_coeff)
1086  centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
1087 
1088  CALL get_qs_env(qs_env, matrix_s=matrix_s)
1089 
1090  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  para_env => qs_loc_env%para_env
1095 
1096  CALL cp_fm_get_info(moloc_coeff(ispin), ncol_global=nmoloc)
1097  ALLOCATE (moment_set(nm, nmoloc))
1098  moment_set = 0.0_dp
1099 
1100  mo_localized => moloc_coeff(ispin)
1101 
1102  DO istate = 1, nmoloc
1103  rcc = centers(1:3, istate)
1104  !rcc = 0.0_dp
1105 
1106  ALLOCATE (moments(nm))
1107  DO imoment = 1, nm
1108  ALLOCATE (moments(imoment)%matrix)
1109  CALL dbcsr_copy(moments(imoment)%matrix, matrix_s(ispin)%matrix, 'MOM MAT')
1110  CALL dbcsr_set(moments(imoment)%matrix, 0.0_dp)
1111  END DO
1112  !
1113 
1114  CALL build_local_moment_matrix(qs_env, moments, norder, rcc)
1115 
1116  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  context=mo_localized%matrix_struct%context)
1120  CALL cp_fm_create(mvector, fm_struct, name="mvector")
1121  CALL cp_fm_create(omvector, fm_struct, name="omvector")
1122  CALL cp_fm_create(momv, fm_struct, name="omvector")
1123  CALL cp_fm_struct_release(fm_struct)
1124 
1125  !cp_fm_to_fm_columns(msource, mtarget, ncol, source_start,target_start)
1126  CALL cp_fm_to_fm(mo_localized, mvector, 1, istate, 1)
1127 
1128  DO imoment = 1, nm
1129  CALL cp_dbcsr_sm_fm_multiply(moments(imoment)%matrix, mvector, omvector, 1)
1130  CALL cp_fm_schur_product(mvector, omvector, momv)
1131  !moment_set(imoment, istate) = moment_set(imoment, istate) + SUM(momv%local_data)
1132  moment_set(imoment, istate) = sum(momv%local_data)
1133  END DO
1134  !
1135  CALL cp_fm_release(mvector)
1136  CALL cp_fm_release(omvector)
1137  CALL cp_fm_release(momv)
1138 
1139  DO imoment = 1, nm
1140  CALL dbcsr_deallocate_matrix(moments(imoment)%matrix)
1141  END DO
1142  DEALLOCATE (moments)
1143  END DO
1144 
1145  CALL para_env%sum(moment_set)
1146 
1147  IF (output_unit_sm > 0) THEN
1148  WRITE (unit=output_unit_sm, fmt='(A,T13,A,I1)') "#", &
1149  "SECOND MOMENTS OF WANNIER CENTERS FOR SPIN ", ispin
1150  WRITE (unit=output_unit_sm, fmt='(A,T29,A2,5(14x,A2))') "#", "XX", "XY", "XZ", "YY", "YZ", "ZZ"
1151  DO istate = 1, nmoloc
1152  cov = 0.0_dp
1153  cov(1) = moment_set(4, istate) - moment_set(1, istate)*moment_set(1, istate)
1154  cov(2) = moment_set(5, istate) - moment_set(1, istate)*moment_set(2, istate)
1155  cov(3) = moment_set(6, istate) - moment_set(1, istate)*moment_set(3, istate)
1156  cov(4) = moment_set(7, istate) - moment_set(2, istate)*moment_set(2, istate)
1157  cov(5) = moment_set(8, istate) - moment_set(2, istate)*moment_set(3, istate)
1158  cov(6) = moment_set(9, istate) - moment_set(3, istate)*moment_set(3, istate)
1159  WRITE (unit=output_unit_sm, fmt='(T4,A,I6,6(T20,6F16.8))') "Center:", istate, cov(1:6)
1160  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  CALL cp_print_key_finished_output(output_unit_sm, logger, print_loc_section, "WANNIER_SPREADS")
1166  CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
1167 
1168  DEALLOCATE (moment_set)
1169 
1170  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  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  REAL(dp), DIMENSION(:, :), POINTER :: centers
1201  REAL(kind=dp) :: imagpart, lx, ly, lz, realpart
1202  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  TYPE(cp_fm_type), DIMENSION(:), POINTER :: moloc_coeff
1209  TYPE(cp_fm_type), POINTER :: mo_localized
1210  TYPE(cp_logger_type), POINTER :: logger
1211  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  logger => cp_get_default_logger()
1217 
1218  output_unit = cp_print_key_unit_nr(logger, print_loc_section, "PROGRAM_RUN_INFO", &
1219  extension=".locInfo")
1220 
1221  IF (output_unit > 0) THEN
1222  WRITE (output_unit, '(/,T2,A)') &
1223  '!-----------------------------------------------------------------------------!'
1224  WRITE (output_unit, '(T12,A)') "Computing second moments of Wannier functions..."
1225  WRITE (output_unit, '(T2,A)') &
1226  '!-----------------------------------------------------------------------------!'
1227  END IF
1228 
1229  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  middle_name=file_tmp, extension=".data")
1232 
1233  norder = 2*(2*taylor_order + 1)
1234  nm = (6 + 11*norder + 6*norder**2 + norder**3)/6 - 1
1235 
1236  NULLIFY (cell)
1237  CALL get_qs_loc_env(qs_loc_env=qs_loc_env, moloc_coeff=moloc_coeff, cell=cell)
1238  lx = cell%hmat(1, 1)
1239  ly = cell%hmat(2, 2)
1240  lz = cell%hmat(3, 3)
1241 
1242  centers => qs_loc_env%localized_wfn_control%centers_set(ispin)%array
1243 
1244  para_env => qs_loc_env%para_env
1245 
1246  CALL get_qs_env(qs_env, matrix_s=matrix_s)
1247 
1248  CALL cp_fm_get_info(moloc_coeff(ispin), ncol_global=nmoloc)
1249 
1250  ALLOCATE (moment_set(nm, nmoloc))
1251  moment_set = 0.0_dp
1252 
1253  mo_localized => moloc_coeff(ispin)
1254 
1255  DO istate = 1, nmoloc
1256  rcc = centers(1:3, istate)
1257  !rcc = 0.0_dp
1258 
1259  ALLOCATE (moments(nm))
1260  DO imoment = 1, nm
1261  ALLOCATE (moments(imoment)%matrix)
1262  CALL dbcsr_copy(moments(imoment)%matrix, matrix_s(ispin)%matrix, 'MOM MAT')
1263  CALL dbcsr_set(moments(imoment)%matrix, 0.0_dp)
1264  END DO
1265  !
1266 
1267  CALL build_local_moment_matrix(qs_env, moments, norder, rcc)
1268 
1269  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  context=mo_localized%matrix_struct%context)
1273  CALL cp_fm_create(mvector, fm_struct, name="mvector")
1274  CALL cp_fm_create(omvector, fm_struct, name="omvector")
1275  CALL cp_fm_create(momv, fm_struct, name="omvector")
1276  CALL cp_fm_struct_release(fm_struct)
1277 
1278  ! cp_fm_to_fm_columns(msource, mtarget, ncol, source_start,target_start)
1279  CALL cp_fm_to_fm(mo_localized, mvector, 1, istate, 1)
1280 
1281  DO imoment = 1, nm
1282  CALL cp_dbcsr_sm_fm_multiply(moments(imoment)%matrix, mvector, omvector, 1)
1283  CALL cp_fm_schur_product(mvector, omvector, momv)
1284  moment_set(imoment, istate) = sum(momv%local_data)
1285  END DO
1286  !
1287  CALL cp_fm_release(mvector)
1288  CALL cp_fm_release(omvector)
1289  CALL cp_fm_release(momv)
1290 
1291  DO imoment = 1, nm
1292  CALL dbcsr_deallocate_matrix(moments(imoment)%matrix)
1293  END DO
1294  DEALLOCATE (moments)
1295  END DO
1296 
1297  CALL para_env%sum(moment_set)
1298 
1299  IF (output_unit_sm > 0) THEN
1300  WRITE (unit=output_unit_sm, fmt='(A,T13,A,I1)') "#", &
1301  "SECOND MOMENTS OF WANNIER CENTERS FOR SPIN ", ispin
1302  WRITE (unit=output_unit_sm, fmt='(A,T29,A2,5(14x,A2))') "#", "XX", "XY", "XZ", "YY", "YZ", "ZZ"
1303  DO istate = 1, nmoloc
1304  cov = 0.0_dp
1305  DO icov = 1, 6
1306  realpart = 0.0_dp
1307  imagpart = 0.0_dp
1308  z = cmplx(realpart, imagpart, dp)
1309  SELECT CASE (icov)
1310 
1311  !! XX
1312  CASE (1)
1313  realpart = 1.0 - 0.5*twopi*twopi/lx/lx/lx/lx &
1314  *moment_set(20, istate)
1315  imagpart = twopi/lx/lx*moment_set(4, istate) &
1316  - 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  z = cmplx(realpart, imagpart, dp)
1324  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  *moment_set(23, istate)
1330  imagpart = twopi/lx/ly*moment_set(5, istate) &
1331  - 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  z = cmplx(realpart, imagpart, dp)
1339  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  *moment_set(25, istate)
1345  imagpart = twopi/lx/lz*moment_set(6, istate) &
1346  - 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  z = cmplx(realpart, imagpart, dp)
1354  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  *moment_set(30, istate)
1360  imagpart = twopi/ly/ly*moment_set(7, istate) &
1361  - 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  z = cmplx(realpart, imagpart, dp)
1369  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  *moment_set(32, istate)
1375  imagpart = twopi/ly/lz*moment_set(8, istate) &
1376  - 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  z = cmplx(realpart, imagpart, dp)
1384  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  *moment_set(34, istate)
1390  imagpart = twopi/lz/lz*moment_set(9, istate) &
1391  - 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  z = cmplx(realpart, imagpart, dp)
1399  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  WRITE (unit=output_unit_sm, fmt='(T4,A,I6,6(T20,6F16.8))') "Center:", istate, cov(1:6)
1404  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  CALL cp_print_key_finished_output(output_unit_sm, logger, print_loc_section, "WANNIER_SPREADS")
1410  CALL cp_print_key_finished_output(output_unit, logger, print_loc_section, "PROGRAM_RUN_INFO")
1411 
1412  DEALLOCATE (moment_set)
1413 
1414  END SUBROUTINE centers_second_moments_berry
1415 
1416 END MODULE qs_loc_methods
subroutine pbc(r, r_pbc, s, s_pbc, a, b, c, alpha, beta, gamma, debug, info, pbc0, h, hinv)
...
Definition: dumpdcd.F:1203
Define the atomic kind types and their sub types.
subroutine, public set_atomic_kind(atomic_kind, element_symbol, name, mass, kind_number, natom, atom_list, fist_potential, shell, shell_active, damping)
Set the components of an atomic kind data set.
subroutine, public deallocate_atomic_kind_set(atomic_kind_set)
Destructor routine for a set of atomic kinds.
Handles all functions related to the CELL.
Definition: cell_types.F:15
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_sm_fm_multiply(matrix, fm_in, fm_out, ncol, alpha, beta)
multiply a dbcsr with a fm matrix
subroutine, public copy_dbcsr_to_fm(matrix, fm)
Copy a DBCSR matrix to a BLACS matrix.
basic linear algebra operations for full matrices
subroutine, public cp_fm_rot_rows(matrix, irow, jrow, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th rows.
subroutine, public cp_fm_rot_cols(matrix, icol, jcol, cs, sn)
Applies a planar rotation defined by cs and sn to the i'th and j'th columnns.
subroutine, public cp_fm_schur_product(matrix_a, matrix_b, matrix_c)
computes the schur product of two matrices c_ij = a_ij * b_ij
pool for for elements that are retained and released
subroutine, public fm_pool_create_fm(pool, element, name)
returns an element, allocating it if none is in the pool
represent the structure of a full matrix
Definition: cp_fm_struct.F:14
subroutine, public cp_fm_struct_create(fmstruct, para_env, context, nrow_global, ncol_global, nrow_block, ncol_block, descriptor, first_p_pos, local_leading_dimension, template_fmstruct, square_blocks, force_block)
allocates and initializes a full matrix structure
Definition: cp_fm_struct.F:125
subroutine, public cp_fm_struct_release(fmstruct)
releases a full matrix structure
Definition: cp_fm_struct.F:320
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_get_element(matrix, irow_global, icol_global, alpha, local)
returns an element of a fm this value is valid on every cpu using this call is expensive
Definition: cp_fm_types.F:643
subroutine, public cp_fm_maxabsval(matrix, a_max, ir_max, ic_max)
find the maximum absolute value of the matrix element maxval(abs(matrix))
Definition: cp_fm_types.F:1064
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
various routines to log and control the output. The idea is that decisions about where to log should ...
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
character(len=default_string_length) function, public cp_iter_string(iter_info, print_key, for_file)
returns the iteration string, a string that is useful to create unique filenames (once you trim it)
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
integer, parameter, public cp_p_file
integer function, public cp_print_key_should_output(iteration_info, basis_section, print_key_path, used_print_key, first_time)
returns what should be done with the given property if btest(res,cp_p_store) then the property should...
A wrapper around pw_to_cube() which accepts particle_list_type.
subroutine, public cp_pw_to_cube(pw, unit_nr, title, particles, stride, zero_tails, silent, mpi_io)
...
unit conversion facility
Definition: cp_units.F:30
real(kind=dp) function, public cp_unit_from_cp2k(value, unit_str, defaults, power)
converts from the internal cp2k units to the given unit
Definition: cp_units.F:1179
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public do_loc_jacobi
integer, parameter, public do_loc_l1_norm_sd
integer, parameter, public do_loc_none
integer, parameter, public dump_xmol
integer, parameter, public do_loc_scdm
integer, parameter, public do_loc_crazy
integer, parameter, public do_loc_gapo
integer, parameter, public dump_dcd_aligned_cell
integer, parameter, public dump_dcd
integer, parameter, public do_loc_direct
objects that represent the structure of input sections and the data contained in an input section
integer function, dimension(:), pointer, public section_get_ivals(section_vals, keyword_name)
...
integer function, public section_get_ival(section_vals, keyword_name)
...
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
integer, parameter, public default_string_length
Definition: kinds.F:57
integer, parameter, public default_path_length
Definition: kinds.F:58
integer, parameter, public sp
Definition: kinds.F:33
Machine interface based on Fortran 2003 and POSIX.
Definition: machine.F:17
subroutine, public m_flush(lunit)
flushes units if the &GLOBAL flag is set accordingly
Definition: machine.F:106
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public twopi
Interface to the message passing library MPI.
Output Utilities for MOTION_SECTION.
Definition: motion_utils.F:13
subroutine, public get_output_format(section, path, my_form, my_ext)
Info on the unit to be opened to dump MD informations.
Definition: motion_utils.F:485
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
basic linear algebra operations for full matrixes
represent a simple array based list of the given type
Define methods related to particle_type.
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
subroutine, public write_particle_coordinates(particle_set, iunit, output_format, content, title, cell, array, unit_conv, charge_occup, charge_beta, charge_extended, print_kind)
Should be able to write a few formats e.g. xmol, and some binary format (dcd) some format can be used...
Define the data structure for the particle information.
subroutine, public deallocate_particle_set(particle_set)
Deallocate a particle set.
subroutine, public allocate_particle_set(particle_set, nparticle)
Allocate a particle set.
Definition of physical constants:
Definition: physcon.F:68
real(kind=dp), parameter, public angstrom
Definition: physcon.F:144
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition: pw_pool_types.F:24
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_wavefunction(mo_vectors, ivector, rho, rho_gspace, atomic_kind_set, qs_kind_set, cell, dft_control, particle_set, pw_env, basis_type, external_vector)
maps a given wavefunction on the grid
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
Driver for the localization that should be general for all the methods available and all the definiti...
subroutine, public optimize_loc_pipek(qs_env, method, qs_loc_env, vectors, zij_fm_set, ispin, print_loc_section)
...
subroutine, public centers_second_moments_berry(qs_env, qs_loc_env, print_loc_section, ispin)
Compute the second moments of the centers using a periodic quadrupole operator.
subroutine, public centers_second_moments_loc(qs_env, qs_loc_env, print_loc_section, ispin)
Compute the second moments of the centers using the local (non-periodic) pos operators.
subroutine, public optimize_loc_berry(method, qs_loc_env, vectors, op_sm_set, zij_fm_set, para_env, cell, weights, ispin, print_loc_section, restricted, nextra, nmo, vectors_2, guess_mos)
Calculate and optimize the spread functional as calculated from the selected mos and by the definitio...
subroutine, public centers_spreads_berry(qs_loc_env, zij, nmoloc, cell, weights, ispin, print_loc_section, only_initial_out)
...
subroutine, public qs_print_cubes(qs_env, mo_coeff, nstates, state_list, centers, print_key, root, ispin, idir, state0, file_position)
write the cube files for a set of selected states
subroutine, public jacobi_rotation_pipek(zij_fm_set, vectors, sweeps)
2by2 rotation for the pipek operator in this case the z_ij numbers are reals
New version of the module for the localization of the molecular orbitals This should be able to use d...
Definition: qs_loc_types.F:25
subroutine, public get_qs_loc_env(qs_loc_env, cell, local_molecules, localized_wfn_control, moloc_coeff, op_sm_set, op_fm_set, para_env, particle_set, weights, dim_op)
...
Definition: qs_loc_types.F:317
Localization methods such as 2x2 Jacobi rotations Steepest Decents Conjugate Gradient.
subroutine, public rotate_orbitals(rmat, vectors)
...
subroutine, public jacobi_rotations(weights, zij, vectors, para_env, max_iter, eps_localization, sweeps, out_each, target_time, start_time, restricted)
wrapper for the jacobi routines, should be removed if jacobi_rot_para can deal with serial para_envs.
subroutine, public zij_matrix(vectors, op_sm_set, zij_fm_set)
...
subroutine, public direct_mini(weights, zij, vectors, max_iter, eps_localization, iterations)
use the exponential parametrization as described in to perform a direct mini Gerd Berghold et al....
subroutine, public crazy_rotations(weights, zij, vectors, max_iter, max_crazy_angle, crazy_scale, crazy_use_diag, eps_localization, iterations, converged)
yet another crazy try, computes the angles needed to rotate the orbitals first and rotates them all a...
subroutine, public approx_l1_norm_sd(C, iterations, eps, converged, sweeps)
...
subroutine, public scdm_qrfact(vectors)
...
subroutine, public jacobi_cg_edf_ls(para_env, weights, zij, vectors, max_iter, eps_localization, iter, out_each, nextra, do_cg, nmo, vectors_2, mos_guess)
combine jacobi rotations (serial) and conjugate gradient with golden section line search for partiall...
wrapper for the pools of matrixes
subroutine, public mpools_get(mpools, ao_mo_fm_pools, ao_ao_fm_pools, mo_mo_fm_pools, ao_mosub_fm_pools, mosub_mosub_fm_pools, maxao_maxmo_fm_pool, maxao_maxao_fm_pool, maxmo_maxmo_fm_pool)
returns various attributes of the mpools (notably the pools contained in it)
Calculates the moment integrals <a|r^m|b> and <a|r x d/dr|b>
Definition: qs_moments.F:14
subroutine, public build_local_moment_matrix(qs_env, moments, nmoments, ref_point, ref_points, basis_type)
...
Definition: qs_moments.F:558
types that represent a quickstep subsys
subroutine, public qs_subsys_get(subsys, atomic_kinds, atomic_kind_set, particles, particle_set, local_particles, molecules, molecule_set, molecule_kinds, molecule_kind_set, local_molecules, para_env, colvar_p, shell_particles, core_particles, gci, multipoles, natom, nparticle, ncore, nshell, nkind, atprop, virial, results, cell, cell_ref, use_ref_cell, energy, force, qs_kind_set, cp_subsys, nelectron_total, nelectron_spin)
...
Utilities for string manipulations.
elemental subroutine, public xstring(string, ia, ib)
...