LCOV - code coverage report
Current view: top level - src - qs_loc_methods.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 36.7 % 528 194
Test Date: 2025-07-25 12:55:17 Functions: 45.5 % 11 5

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

Generated by: LCOV version 2.0-1