LCOV - code coverage report
Current view: top level - src - qs_dcdr_utils.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 90.3 % 475 429
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 8 8

            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 Calculate the derivatives of the MO coefficients wrt nuclear coordinates
      10              : !> \author Sandra Luber, Edward Ditler
      11              : ! **************************************************************************************************
      12              : 
      13              : MODULE qs_dcdr_utils
      14              : !#include "./common/cp_common_uses.f90"
      15              :    USE cell_types,                      ONLY: cell_type,&
      16              :                                               get_cell
      17              :    USE cp_control_types,                ONLY: dft_control_type
      18              :    USE cp_dbcsr_api,                    ONLY: dbcsr_copy,&
      19              :                                               dbcsr_create,&
      20              :                                               dbcsr_init_p,&
      21              :                                               dbcsr_p_type,&
      22              :                                               dbcsr_set,&
      23              :                                               dbcsr_type,&
      24              :                                               dbcsr_type_no_symmetry,&
      25              :                                               dbcsr_type_symmetric
      26              :    USE cp_dbcsr_cp2k_link,              ONLY: cp_dbcsr_alloc_block_from_nbl
      27              :    USE cp_dbcsr_operations,             ONLY: cp_dbcsr_sm_fm_multiply,&
      28              :                                               dbcsr_allocate_matrix_set,&
      29              :                                               dbcsr_deallocate_matrix_set
      30              :    USE cp_files,                        ONLY: close_file,&
      31              :                                               open_file
      32              :    USE cp_fm_struct,                    ONLY: cp_fm_struct_create,&
      33              :                                               cp_fm_struct_release
      34              :    USE cp_fm_types,                     ONLY: cp_fm_create,&
      35              :                                               cp_fm_get_info,&
      36              :                                               cp_fm_get_submatrix,&
      37              :                                               cp_fm_release,&
      38              :                                               cp_fm_set_all,&
      39              :                                               cp_fm_set_submatrix,&
      40              :                                               cp_fm_to_fm,&
      41              :                                               cp_fm_type
      42              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      43              :                                               cp_logger_get_default_io_unit,&
      44              :                                               cp_logger_type,&
      45              :                                               cp_to_string
      46              :    USE cp_output_handling,              ONLY: cp_p_file,&
      47              :                                               cp_print_key_finished_output,&
      48              :                                               cp_print_key_generate_filename,&
      49              :                                               cp_print_key_should_output,&
      50              :                                               cp_print_key_unit_nr
      51              :    USE cp_result_methods,               ONLY: get_results
      52              :    USE cp_result_types,                 ONLY: cp_result_type
      53              :    USE input_constants,                 ONLY: current_orb_center_wannier,&
      54              :                                               use_mom_ref_user
      55              :    USE input_section_types,             ONLY: section_vals_get_subs_vals,&
      56              :                                               section_vals_type,&
      57              :                                               section_vals_val_get
      58              :    USE kinds,                           ONLY: default_path_length,&
      59              :                                               default_string_length,&
      60              :                                               dp
      61              :    USE memory_utilities,                ONLY: reallocate
      62              :    USE message_passing,                 ONLY: mp_para_env_type
      63              :    USE molecule_types,                  ONLY: molecule_type
      64              :    USE moments_utils,                   ONLY: get_reference_point
      65              :    USE parallel_gemm_api,               ONLY: parallel_gemm
      66              :    USE particle_types,                  ONLY: particle_type
      67              :    USE qs_environment_types,            ONLY: get_qs_env,&
      68              :                                               qs_environment_type
      69              :    USE qs_kinetic,                      ONLY: build_kinetic_matrix
      70              :    USE qs_ks_types,                     ONLY: qs_ks_env_type
      71              :    USE qs_linres_types,                 ONLY: dcdr_env_type,&
      72              :                                               linres_control_type
      73              :    USE qs_loc_types,                    ONLY: get_qs_loc_env,&
      74              :                                               localized_wfn_control_type,&
      75              :                                               qs_loc_env_type
      76              :    USE qs_mo_types,                     ONLY: get_mo_set,&
      77              :                                               mo_set_type
      78              :    USE qs_moments,                      ONLY: build_local_moment_matrix
      79              :    USE qs_neighbor_list_types,          ONLY: neighbor_list_set_p_type
      80              :    USE qs_overlap,                      ONLY: build_overlap_matrix
      81              :    USE string_utilities,                ONLY: xstring
      82              : #include "./base/base_uses.f90"
      83              : 
      84              :    IMPLICIT NONE
      85              : 
      86              :    PRIVATE
      87              :    PUBLIC :: dcdr_env_cleanup, dcdr_env_init, dcdr_print, &
      88              :              get_loc_setting, shift_wannier_into_cell, &
      89              :              dcdr_write_restart, dcdr_read_restart, &
      90              :              multiply_localization
      91              : 
      92              :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_dcdr_utils'
      93              : 
      94              :    REAL(dp), DIMENSION(3, 3, 3), PARAMETER  :: Levi_Civita = RESHAPE((/ &
      95              :                                                           0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, &
      96              :                                                           0.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, -1.0_dp, 0.0_dp, 0.0_dp, &
      97              :                                                         0.0_dp, -1.0_dp, 0.0_dp, 1.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp, 0.0_dp/), &
      98              :                                                                      (/3, 3, 3/))
      99              : 
     100              : CONTAINS
     101              : 
     102              : ! **************************************************************************************************
     103              : !> \brief Multiply (ao_matrix @ mo_coeff) and store the column icenter in res
     104              : !> \param ao_matrix ...
     105              : !> \param mo_coeff ...
     106              : !> \param work Working space
     107              : !> \param nmo ...
     108              : !> \param icenter ...
     109              : !> \param res ...
     110              : !> \author Edward Ditler
     111              : ! **************************************************************************************************
     112          864 :    SUBROUTINE multiply_localization(ao_matrix, mo_coeff, work, nmo, icenter, res)
     113              :       TYPE(dbcsr_type), INTENT(IN), POINTER              :: ao_matrix
     114              :       TYPE(cp_fm_type), INTENT(IN)                       :: mo_coeff
     115              :       TYPE(cp_fm_type), INTENT(INOUT)                    :: work
     116              :       INTEGER, INTENT(IN)                                :: nmo, icenter
     117              :       TYPE(cp_fm_type), INTENT(IN)                       :: res
     118              : 
     119              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'multiply_localization'
     120              : 
     121              :       INTEGER                                            :: handle
     122              : 
     123          864 :       CALL timeset(routineN, handle)
     124              : 
     125              :       ! Multiply by the MO coefficients
     126          864 :       CALL cp_dbcsr_sm_fm_multiply(ao_matrix, mo_coeff, work, ncol=nmo)
     127              : 
     128              :       ! Only keep the icenter-th column
     129          864 :       CALL cp_fm_to_fm(work, res, 1, icenter, icenter)
     130              : 
     131              :       ! Reset the matrices
     132          864 :       CALL cp_fm_set_all(work, 0.0_dp)
     133              : 
     134          864 :       CALL timestop(handle)
     135          864 :    END SUBROUTINE multiply_localization
     136              : 
     137              : ! **************************************************************************************************
     138              : !> \brief Copied from linres_read_restart
     139              : !> \param qs_env ...
     140              : !> \param linres_section ...
     141              : !> \param vec ...
     142              : !> \param lambda ...
     143              : !> \param beta ...
     144              : !> \param tag ...
     145              : !> \note Adapted from linres_read_restart (ED)
     146              : !>       Would be nice not to crash but to start from zero if the present file doesn't match.
     147              : ! **************************************************************************************************
     148           18 :    SUBROUTINE dcdr_read_restart(qs_env, linres_section, vec, lambda, beta, tag)
     149              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     150              :       TYPE(section_vals_type), POINTER                   :: linres_section
     151              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: vec
     152              :       INTEGER, INTENT(IN)                                :: lambda, beta
     153              :       CHARACTER(LEN=*)                                   :: tag
     154              : 
     155              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'dcdr_read_restart'
     156              : 
     157              :       CHARACTER(LEN=default_path_length)                 :: filename
     158              :       CHARACTER(LEN=default_string_length)               :: my_middle
     159              :       INTEGER :: beta_tmp, handle, i, i_block, ia, ie, iostat, iounit, ispin, j, lambda_tmp, &
     160              :          max_block, n_rep_val, nao, nao_tmp, nmo, nmo_tmp, nspins, nspins_tmp, rst_unit
     161              :       LOGICAL                                            :: file_exists
     162           18 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
     163              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     164              :       TYPE(cp_logger_type), POINTER                      :: logger
     165           18 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     166              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     167              :       TYPE(section_vals_type), POINTER                   :: print_key
     168              : 
     169           18 :       file_exists = .FALSE.
     170              : 
     171           18 :       CALL timeset(routineN, handle)
     172              : 
     173           18 :       NULLIFY (mos, para_env, logger, print_key, vecbuffer)
     174           18 :       logger => cp_get_default_logger()
     175              : 
     176              :       iounit = cp_print_key_unit_nr(logger, linres_section, &
     177           18 :                                     "PRINT%PROGRAM_RUN_INFO", extension=".Log")
     178              : 
     179              :       CALL get_qs_env(qs_env=qs_env, &
     180              :                       para_env=para_env, &
     181           18 :                       mos=mos)
     182              : 
     183           18 :       nspins = SIZE(mos)
     184              : 
     185           18 :       rst_unit = -1
     186           18 :       IF (para_env%is_source()) THEN
     187              :          CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", &
     188            9 :                                    n_rep_val=n_rep_val)
     189              : 
     190            9 :          CALL XSTRING(tag, ia, ie)
     191              :          my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
     192            9 :                      //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
     193              : 
     194            9 :          IF (n_rep_val > 0) THEN
     195            0 :             CALL section_vals_val_get(linres_section, "WFN_RESTART_FILE_NAME", c_val=filename)
     196            0 :             CALL xstring(filename, ia, ie)
     197            0 :             filename = filename(ia:ie)//TRIM(my_middle)//".lr"
     198              :          ELSE
     199              :             ! try to read from the filename that is generated automatically from the printkey
     200            9 :             print_key => section_vals_get_subs_vals(linres_section, "PRINT%RESTART")
     201              :             filename = cp_print_key_generate_filename(logger, print_key, &
     202            9 :                                                       extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
     203              :          END IF
     204            9 :          INQUIRE (FILE=filename, exist=file_exists)
     205              :          !
     206              :          ! open file
     207            9 :          IF (file_exists) THEN
     208              :             CALL open_file(file_name=TRIM(filename), &
     209              :                            file_action="READ", &
     210              :                            file_form="UNFORMATTED", &
     211              :                            file_position="REWIND", &
     212              :                            file_status="OLD", &
     213            0 :                            unit_number=rst_unit)
     214              : 
     215            0 :             IF (iounit > 0) WRITE (iounit, "(T2,A)") &
     216            0 :                "LINRES| Reading response wavefunctions from the restart file <"//TRIM(ADJUSTL(filename))//">"
     217              :          ELSE
     218            9 :             IF (iounit > 0) WRITE (iounit, "(T2,A)") &
     219            9 :                "LINRES| Restart file  <"//TRIM(ADJUSTL(filename))//"> not found"
     220              :          END IF
     221              :       END IF
     222              : 
     223           18 :       CALL para_env%bcast(file_exists)
     224              : 
     225           18 :       IF (file_exists) THEN
     226              : 
     227            0 :          CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
     228            0 :          CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
     229              : 
     230            0 :          ALLOCATE (vecbuffer(nao, max_block))
     231              :          !
     232              :          ! read headers
     233            0 :          IF (rst_unit > 0) READ (rst_unit, IOSTAT=iostat) lambda_tmp, beta_tmp, nspins_tmp, nao_tmp
     234            0 :          CALL para_env%bcast(iostat)
     235            0 :          CALL para_env%bcast(beta_tmp)
     236            0 :          CALL para_env%bcast(lambda_tmp)
     237            0 :          CALL para_env%bcast(nspins_tmp)
     238            0 :          CALL para_env%bcast(nao_tmp)
     239              : 
     240              :          ! check that the number nao, nmo and nspins are
     241              :          ! the same as in the current mos
     242            0 :          IF (nspins_tmp .NE. nspins) THEN
     243            0 :             CPABORT("nspins not consistent")
     244              :          END IF
     245            0 :          IF (nao_tmp .NE. nao) CPABORT("nao not consistent")
     246              :          ! check that it's the right file
     247              :          ! the same as in the current mos
     248            0 :          IF (lambda_tmp .NE. lambda) CPABORT("lambda not consistent")
     249            0 :          IF (beta_tmp .NE. beta) CPABORT("beta not consistent")
     250              :          !
     251            0 :          DO ispin = 1, nspins
     252            0 :             CALL get_mo_set(mos(ispin), mo_coeff=mo_coeff)
     253            0 :             CALL cp_fm_get_info(mo_coeff, ncol_global=nmo)
     254              :             !
     255            0 :             IF (rst_unit > 0) READ (rst_unit) nmo_tmp
     256            0 :             CALL para_env%bcast(nmo_tmp)
     257            0 :             IF (nmo_tmp .NE. nmo) CPABORT("nmo not consistent")
     258              :             !
     259              :             ! read the response
     260            0 :             DO i = 1, nmo, MAX(max_block, 1)
     261            0 :                i_block = MIN(max_block, nmo - i + 1)
     262            0 :                DO j = 1, i_block
     263            0 :                   IF (rst_unit > 0) READ (rst_unit) vecbuffer(1:nao, j)
     264              :                END DO
     265            0 :                CALL para_env%bcast(vecbuffer)
     266            0 :                CALL cp_fm_set_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
     267              :             END DO
     268              :          END DO
     269              : 
     270            0 :          IF (iostat /= 0) THEN
     271            0 :             IF (iounit > 0) WRITE (iounit, "(T2,A)") &
     272            0 :                "LINRES| Restart file <"//TRIM(ADJUSTL(filename))//"> not found"
     273              :          END IF
     274              : 
     275            0 :          DEALLOCATE (vecbuffer)
     276              : 
     277              :       END IF
     278              : 
     279           18 :       IF (para_env%is_source()) THEN
     280            9 :          IF (file_exists) CALL close_file(unit_number=rst_unit)
     281              :       END IF
     282              : 
     283           18 :       CALL timestop(handle)
     284              : 
     285           18 :    END SUBROUTINE dcdr_read_restart
     286              : 
     287              : ! **************************************************************************************************
     288              : !> \brief Copied from linres_write_restart
     289              : !> \param qs_env ...
     290              : !> \param linres_section ...
     291              : !> \param vec ...
     292              : !> \param lambda ...
     293              : !> \param beta ...
     294              : !> \param tag ...
     295              : !> \note Adapted from linres_read_restart (ED)
     296              : !>       Would be nice not to crash but to start from zero if the present file doesn't match.
     297              : ! **************************************************************************************************
     298           18 :    SUBROUTINE dcdr_write_restart(qs_env, linres_section, vec, lambda, beta, tag)
     299              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     300              :       TYPE(section_vals_type), POINTER                   :: linres_section
     301              :       TYPE(cp_fm_type), DIMENSION(:), INTENT(IN)         :: vec
     302              :       INTEGER, INTENT(IN)                                :: lambda, beta
     303              :       CHARACTER(LEN=*)                                   :: tag
     304              : 
     305              :       CHARACTER(LEN=*), PARAMETER :: routineN = 'dcdr_write_restart'
     306              : 
     307              :       CHARACTER(LEN=default_path_length)                 :: filename
     308              :       CHARACTER(LEN=default_string_length)               :: my_middle, my_pos, my_status
     309              :       INTEGER                                            :: handle, i, i_block, ia, ie, iounit, &
     310              :                                                             ispin, j, max_block, nao, nmo, nspins, &
     311              :                                                             rst_unit
     312           18 :       REAL(KIND=dp), DIMENSION(:, :), POINTER            :: vecbuffer
     313              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     314              :       TYPE(cp_logger_type), POINTER                      :: logger
     315           18 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     316              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     317              :       TYPE(section_vals_type), POINTER                   :: print_key
     318              : 
     319           18 :       NULLIFY (logger, mo_coeff, mos, para_env, print_key, vecbuffer)
     320              : 
     321           18 :       CALL timeset(routineN, handle)
     322              : 
     323           18 :       logger => cp_get_default_logger()
     324              : 
     325           18 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, linres_section, "PRINT%RESTART", &
     326              :                                            used_print_key=print_key), &
     327              :                 cp_p_file)) THEN
     328              : 
     329              :          iounit = cp_print_key_unit_nr(logger, linres_section, &
     330           18 :                                        "PRINT%PROGRAM_RUN_INFO", extension=".Log")
     331              : 
     332              :          CALL get_qs_env(qs_env=qs_env, &
     333              :                          mos=mos, &
     334           18 :                          para_env=para_env)
     335              : 
     336           18 :          nspins = SIZE(mos)
     337              : 
     338           18 :          my_status = "REPLACE"
     339           18 :          my_pos = "REWIND"
     340           18 :          CALL XSTRING(tag, ia, ie)
     341              :          my_middle = "RESTART-"//tag(ia:ie)//TRIM("-")//TRIM(ADJUSTL(cp_to_string(beta))) &
     342           18 :                      //TRIM("-")//TRIM(ADJUSTL(cp_to_string(lambda)))
     343              :          rst_unit = cp_print_key_unit_nr(logger, linres_section, "PRINT%RESTART", &
     344              :                                          extension=".lr", middle_name=TRIM(my_middle), file_status=TRIM(my_status), &
     345           18 :                                          file_position=TRIM(my_pos), file_action="WRITE", file_form="UNFORMATTED")
     346              : 
     347              :          filename = cp_print_key_generate_filename(logger, print_key, &
     348           18 :                                                    extension=".lr", middle_name=TRIM(my_middle), my_local=.FALSE.)
     349              : 
     350           18 :          IF (iounit > 0) THEN
     351              :             WRITE (UNIT=iounit, FMT="(T2,A)") &
     352            9 :                "LINRES| Writing response functions to the restart file <"//TRIM(ADJUSTL(filename))//">"
     353              :          END IF
     354              : 
     355              :          !
     356              :          ! write data to file
     357              :          ! use the scalapack block size as a default for buffering columns
     358           18 :          CALL get_mo_set(mos(1), mo_coeff=mo_coeff)
     359           18 :          CALL cp_fm_get_info(mo_coeff, nrow_global=nao, ncol_block=max_block)
     360           72 :          ALLOCATE (vecbuffer(nao, max_block))
     361              : 
     362           18 :          IF (rst_unit > 0) WRITE (rst_unit) lambda, beta, nspins, nao
     363              : 
     364           36 :          DO ispin = 1, nspins
     365           18 :             CALL cp_fm_get_info(vec(ispin), ncol_global=nmo)
     366              : 
     367           18 :             IF (rst_unit > 0) WRITE (rst_unit) nmo
     368              : 
     369           72 :             DO i = 1, nmo, MAX(max_block, 1)
     370           18 :                i_block = MIN(max_block, nmo - i + 1)
     371           18 :                CALL cp_fm_get_submatrix(vec(ispin), vecbuffer, 1, i, nao, i_block)
     372              :                ! doing this in one write would increase efficiency, but breaks RESTART compatibility.
     373              :                ! to old ones, and in cases where max_block is different between runs, as might happen during
     374              :                ! restarts with a different number of CPUs
     375          108 :                DO j = 1, i_block
     376           90 :                   IF (rst_unit > 0) WRITE (rst_unit) vecbuffer(1:nao, j)
     377              :                END DO
     378              :             END DO
     379              :          END DO
     380              : 
     381           18 :          DEALLOCATE (vecbuffer)
     382              : 
     383              :          CALL cp_print_key_finished_output(rst_unit, logger, linres_section, &
     384           36 :                                            "PRINT%RESTART")
     385              :       END IF
     386              : 
     387           18 :       CALL timestop(handle)
     388              : 
     389           18 :    END SUBROUTINE dcdr_write_restart
     390              : 
     391              : ! **************************************************************************************************
     392              : !> \brief Print the APT and sum rules
     393              : !> \param dcdr_env ...
     394              : !> \param qs_env ...
     395              : !> \author Edward Ditler
     396              : ! **************************************************************************************************
     397           22 :    SUBROUTINE dcdr_print(dcdr_env, qs_env)
     398              :       TYPE(dcdr_env_type)                                :: dcdr_env
     399              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     400              : 
     401              :       CHARACTER(LEN=default_string_length)               :: description
     402              :       INTEGER                                            :: alpha, beta, delta, gamma, i, k, l, &
     403              :                                                             lambda, natom, nsubset, output_unit
     404           22 :       REAL(dp), DIMENSION(:, :, :), POINTER              :: apt_el_dcdr, apt_nuc_dcdr, apt_total_dcdr
     405           22 :       REAL(dp), DIMENSION(:, :, :, :), POINTER           :: apt_center_dcdr, apt_subset_dcdr
     406              :       REAL(kind=dp), DIMENSION(3, 3)                     :: sum_rule_0, sum_rule_1, sum_rule_2
     407              :       TYPE(cp_logger_type), POINTER                      :: logger
     408              :       TYPE(cp_result_type), POINTER                      :: results
     409           22 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     410           22 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     411              :       TYPE(section_vals_type), POINTER                   :: dcdr_section
     412              : 
     413           22 :       NULLIFY (logger)
     414              : 
     415           44 :       logger => cp_get_default_logger()
     416           22 :       output_unit = cp_logger_get_default_io_unit(logger)
     417              : 
     418           22 :       dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
     419              : 
     420           22 :       NULLIFY (particle_set)
     421           22 :       CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, molecule_set=molecule_set)
     422           22 :       natom = SIZE(particle_set)
     423           22 :       nsubset = SIZE(molecule_set)
     424              : 
     425           22 :       apt_el_dcdr => dcdr_env%apt_el_dcdr
     426           22 :       apt_nuc_dcdr => dcdr_env%apt_nuc_dcdr
     427           22 :       apt_total_dcdr => dcdr_env%apt_total_dcdr
     428           22 :       apt_subset_dcdr => dcdr_env%apt_el_dcdr_per_subset
     429           22 :       apt_center_dcdr => dcdr_env%apt_el_dcdr_per_center
     430              : 
     431           22 :       IF (dcdr_env%localized_psi0) THEN
     432            4 :          IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the final apt matrix per atom per subset'
     433           16 :          DO k = 1, natom
     434           52 :             DO l = 1, nsubset
     435           36 :                IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, *) 'APT | Subset', l
     436          156 :                DO i = 1, 3
     437          108 :                   IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,F15.6,F15.6,F15.6)") &
     438           90 :                      'APT | apt_subset ', i, apt_subset_dcdr(i, :, k, l)
     439              :                END DO
     440              :             END DO
     441              :          END DO
     442              :       END IF
     443              : 
     444           22 :       IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") &
     445           11 :          'APT | Write the final apt matrix per atom (Position perturbation)'
     446           88 :       DO l = 1, natom
     447           66 :          IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,I3,A,F15.6)") &
     448           33 :             'APT | Atom', l, ' - GAPT ', &
     449              :             (apt_total_dcdr(1, 1, l) &
     450              :              + apt_total_dcdr(2, 2, l) &
     451           66 :              + apt_total_dcdr(3, 3, l))/3._dp
     452          286 :          DO i = 1, 3
     453          264 :             IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,F15.6,F15.6,F15.6)") "APT | ", apt_total_dcdr(i, :, l)
     454              :          END DO
     455              :       END DO
     456              : 
     457           22 :       IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | Write the total apt matrix'
     458           88 :       DO i = 1, 3
     459           66 :          IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
     460          451 :                                               "(A,F15.6,F15.6,F15.6)") "APT | ", SUM(apt_total_dcdr(i, :, :), dim=2)
     461              :       END DO
     462           22 :       IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") 'APT | End Write the final apt matrix'
     463              : 
     464              :       ! Get the dipole
     465           22 :       CALL get_qs_env(qs_env, results=results)
     466           22 :       description = "[DIPOLE]"
     467           22 :       CALL get_results(results=results, description=description, values=dcdr_env%dipole_pos(1:3))
     468              : 
     469              :       ! Sum rules [for all alpha, beta]
     470           22 :       sum_rule_0 = 0._dp
     471           22 :       sum_rule_1 = 0._dp
     472           22 :       sum_rule_2 = 0._dp
     473              : 
     474           88 :       DO alpha = 1, 3
     475          286 :          DO beta = 1, 3
     476              :             ! 0: sum_lambda apt(alpha, beta, lambda)
     477          792 :             DO lambda = 1, natom
     478              :                sum_rule_0(alpha, beta) = sum_rule_0(alpha, beta) &
     479          792 :                                          + apt_total_dcdr(alpha, beta, lambda)
     480              :             END DO
     481              : 
     482              :             ! 1: sum_gamma epsilon_(alpha beta gamma) mu_gamma
     483          792 :             DO gamma = 1, 3
     484              :                sum_rule_1(alpha, beta) = sum_rule_1(alpha, beta) &
     485          792 :                                          + Levi_Civita(alpha, beta, gamma)*dcdr_env%dipole_pos(gamma)
     486              :             END DO
     487              : 
     488              :             ! 2: sum_(lambda gamma delta) R^lambda_gamma apt(delta, alpha, lambda)
     489          858 :             DO lambda = 1, natom
     490         2574 :                DO gamma = 1, 3
     491         7722 :                   DO delta = 1, 3
     492              :                      sum_rule_2(alpha, beta) = sum_rule_2(alpha, beta) &
     493              :                                                + Levi_Civita(beta, gamma, delta) &
     494              :                                                *particle_set(lambda)%r(gamma) &
     495         7128 :                                                *apt_total_dcdr(delta, alpha, lambda)
     496              :                   END DO
     497              :                END DO
     498              :             END DO
     499              : 
     500              :          END DO ! beta
     501              :       END DO   ! alpha
     502              : 
     503           22 :       IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A)") "APT | Position perturbation sum rules"
     504           22 :       IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, "(A,T18,A,T34,A,T49,A)") &
     505           11 :          "APT |", "Total APT", "Dipole", "R * APT"
     506           88 :       DO alpha = 1, 3
     507          286 :          DO beta = 1, 3
     508          198 :             IF (dcdr_env%output_unit > 0) WRITE (dcdr_env%output_unit, &
     509              :                                                  "(A,I3,I3,F15.6,F15.6,F15.6)") &
     510           99 :                "APT | ", &
     511           99 :                alpha, beta, &
     512           99 :                sum_rule_0(alpha, beta), &
     513           99 :                sum_rule_1(alpha, beta), &
     514          264 :                sum_rule_2(alpha, beta)
     515              :          END DO
     516              :       END DO
     517              : 
     518           22 :    END SUBROUTINE dcdr_print
     519              : 
     520              : ! **************************************************************************************************
     521              : !> \brief ...
     522              : !> \param r ...
     523              : !> \param cell ...
     524              : !> \param r_shifted ...
     525              : ! **************************************************************************************************
     526          144 :    SUBROUTINE shift_wannier_into_cell(r, cell, r_shifted)
     527              :       REAL(dp), DIMENSION(3), INTENT(in)                 :: r
     528              :       TYPE(cell_type), INTENT(in), POINTER               :: cell
     529              :       REAL(dp), DIMENSION(3), INTENT(out)                :: r_shifted
     530              : 
     531              :       INTEGER                                            :: i
     532              :       REAL(kind=dp), DIMENSION(3)                        :: abc
     533              : 
     534              :       ! Only orthorombic cell for now
     535          144 :       CALL get_cell(cell, abc=abc)
     536              : 
     537          576 :       DO i = 1, 3
     538          576 :          IF (r(i) < 0._dp) THEN
     539          234 :             r_shifted(i) = r(i) + abc(i)
     540          198 :          ELSE IF (r(i) > abc(i)) THEN
     541            0 :             r_shifted(i) = r(i) - abc(i)
     542              :          ELSE
     543          198 :             r_shifted(i) = r(i)
     544              :          END IF
     545              :       END DO
     546          144 :    END SUBROUTINE shift_wannier_into_cell
     547              : 
     548              : ! **************************************************************************************************
     549              : !> \brief ...
     550              : !> \param dcdr_env ...
     551              : !> \param qs_env ...
     552              : ! **************************************************************************************************
     553            4 :    SUBROUTINE get_loc_setting(dcdr_env, qs_env)
     554              :       TYPE(dcdr_env_type)                                :: dcdr_env
     555              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     556              : 
     557              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'get_loc_setting'
     558              : 
     559              :       INTEGER                                            :: handle, is, ispin, istate, max_states, &
     560              :                                                             nmo, nmoloc, nstate, nstate_list(2)
     561            4 :       INTEGER, ALLOCATABLE, DIMENSION(:, :)              :: state_list
     562            4 :       REAL(dp), DIMENSION(:, :), POINTER                 :: center_array
     563              :       TYPE(linres_control_type), POINTER                 :: linres_control
     564              :       TYPE(localized_wfn_control_type), POINTER          :: localized_wfn_control
     565            4 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     566              :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
     567              :       TYPE(section_vals_type), POINTER                   :: dcdr_section
     568              : 
     569            4 :       CALL timeset(routineN, handle)
     570              : 
     571              :       CALL get_qs_env(qs_env=qs_env, &
     572              :                       linres_control=linres_control, &
     573            4 :                       mos=mos)
     574              : 
     575              :       ! Some checks
     576            4 :       max_states = 0
     577            4 :       CALL get_mo_set(mo_set=mos(1), nmo=nmo)
     578            4 :       max_states = MAX(max_states, nmo)
     579              : 
     580              :       ! check that the number of localized states is equal to the number of states
     581            4 :       nmoloc = SIZE(linres_control%qs_loc_env%localized_wfn_control%centers_set(1)%array, 2)
     582            4 :       IF (nmoloc .NE. nmo) THEN
     583            0 :          CPABORT("The number of localized functions is not equal to the number of states.")
     584              :       END IF
     585              : 
     586              :       ! which center for the orbitals shall we use
     587            4 :       dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
     588            4 :       CALL section_vals_val_get(dcdr_section, "ORBITAL_CENTER", i_val=dcdr_env%orb_center)
     589            8 :       SELECT CASE (dcdr_env%orb_center)
     590              :       CASE (current_orb_center_wannier)
     591            4 :          dcdr_env%orb_center_name = "WANNIER"
     592              :       CASE DEFAULT
     593            4 :          CPABORT(" ")
     594              :       END SELECT
     595              : 
     596            4 :       qs_loc_env => linres_control%qs_loc_env
     597            4 :       CALL get_qs_loc_env(qs_loc_env, localized_wfn_control=localized_wfn_control)
     598              : 
     599           16 :       ALLOCATE (dcdr_env%centers_set(dcdr_env%nspins))
     600           12 :       ALLOCATE (dcdr_env%center_list(dcdr_env%nspins))
     601           16 :       ALLOCATE (state_list(max_states, dcdr_env%nspins))
     602           24 :       state_list(:, :) = HUGE(0)
     603           12 :       nstate_list(:) = HUGE(0)
     604              : 
     605              :       ! Build the state_list
     606            8 :       DO ispin = 1, dcdr_env%nspins
     607            4 :          center_array => localized_wfn_control%centers_set(ispin)%array
     608            4 :          nstate = 0
     609           20 :          DO istate = 1, SIZE(center_array, 2)
     610           16 :             nstate = nstate + 1
     611           20 :             state_list(nstate, ispin) = istate
     612              :          END DO
     613              :          nstate_list(ispin) = nstate
     614              : 
     615              :          ! clustering the states
     616            4 :          nstate = nstate_list(ispin)
     617            4 :          dcdr_env%nstates(ispin) = nstate
     618              : 
     619           12 :          ALLOCATE (dcdr_env%center_list(ispin)%array(2, nstate + 1))
     620           12 :          ALLOCATE (dcdr_env%centers_set(ispin)%array(3, nstate))
     621           64 :          dcdr_env%center_list(ispin)%array(:, :) = HUGE(0)
     622           68 :          dcdr_env%centers_set(ispin)%array(:, :) = HUGE(0.0_dp)
     623              : 
     624            4 :          center_array => localized_wfn_control%centers_set(ispin)%array
     625              : 
     626              :          ! point to the psi0 centers
     627            4 :          SELECT CASE (dcdr_env%orb_center)
     628              :          CASE (current_orb_center_wannier)
     629              :             ! use the wannier center as -center-
     630            4 :             dcdr_env%nbr_center(ispin) = nstate
     631           20 :             DO is = 1, nstate
     632           16 :                istate = state_list(is, 1)
     633           64 :                dcdr_env%centers_set(ispin)%array(1:3, is) = center_array(1:3, istate)
     634           16 :                dcdr_env%center_list(ispin)%array(1, is) = is
     635           20 :                dcdr_env%center_list(ispin)%array(2, is) = istate
     636              :             END DO
     637            4 :             dcdr_env%center_list(ispin)%array(1, nstate + 1) = nstate + 1
     638              : 
     639              :          CASE DEFAULT
     640            4 :             CPABORT("Unknown orbital center...")
     641              :          END SELECT
     642              :       END DO
     643              : 
     644            4 :       CALL timestop(handle)
     645           12 :    END SUBROUTINE get_loc_setting
     646              : 
     647              : ! **************************************************************************************************
     648              : !> \brief Initialize the dcdr environment
     649              : !> \param dcdr_env ...
     650              : !> \param qs_env ...
     651              : ! **************************************************************************************************
     652           24 :    SUBROUTINE dcdr_env_init(dcdr_env, qs_env)
     653              :       TYPE(dcdr_env_type)                                :: dcdr_env
     654              :       TYPE(qs_environment_type), POINTER                 :: qs_env
     655              : 
     656              :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'dcdr_env_init'
     657              : 
     658              :       INTEGER                                            :: handle, homo, i, isize, ispin, j, jg, &
     659              :                                                             n_rep, nao, natom, nmo, nspins, &
     660              :                                                             nsubset, output_unit, reference, &
     661              :                                                             unit_number
     662           24 :       INTEGER, DIMENSION(:), POINTER                     :: tmplist
     663              :       LOGICAL                                            :: explicit
     664           24 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
     665              :       TYPE(cp_fm_type)                                   :: buf
     666              :       TYPE(cp_fm_type), POINTER                          :: mo_coeff
     667              :       TYPE(cp_logger_type), POINTER                      :: logger
     668           24 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_ks, matrix_s
     669              :       TYPE(dft_control_type), POINTER                    :: dft_control
     670           24 :       TYPE(mo_set_type), DIMENSION(:), POINTER           :: mos
     671           24 :       TYPE(molecule_type), DIMENSION(:), POINTER         :: molecule_set
     672              :       TYPE(mp_para_env_type), POINTER                    :: para_env
     673              :       TYPE(neighbor_list_set_p_type), DIMENSION(:), &
     674           24 :          POINTER                                         :: sab_all, sab_orb
     675           24 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     676              :       TYPE(qs_ks_env_type), POINTER                      :: ks_env
     677              :       TYPE(section_vals_type), POINTER                   :: dcdr_section, loc_section, lr_section
     678              : 
     679           24 :       CALL timeset(routineN, handle)
     680              : 
     681              :       ! Set up the logger
     682           24 :       NULLIFY (logger, loc_section, dcdr_section, lr_section)
     683           24 :       logger => cp_get_default_logger()
     684           24 :       loc_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%LOCALIZE")
     685           24 :       dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
     686              :       dcdr_env%output_unit = cp_print_key_unit_nr(logger, dcdr_section, "PRINT%APT", &
     687              :                                                   extension=".data", middle_name="dcdr", log_filename=.FALSE., &
     688           24 :                                                   file_position="REWIND", file_status="REPLACE")
     689              : 
     690           24 :       lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
     691              :       output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
     692           24 :                                          extension=".linresLog")
     693           24 :       unit_number = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", extension=".linresLog")
     694              : 
     695           24 :       IF (output_unit > 0) THEN
     696           12 :          WRITE (output_unit, "(/,T20,A,/)") "*** Start DCDR calculation ***"
     697              :       END IF
     698              : 
     699           24 :       NULLIFY (ks_env, dft_control, sab_orb, sab_all, particle_set, molecule_set, matrix_s, matrix_ks, mos, para_env)
     700              :       CALL get_qs_env(qs_env=qs_env, &
     701              :                       ks_env=ks_env, &
     702              :                       dft_control=dft_control, &
     703              :                       sab_orb=sab_orb, &
     704              :                       sab_all=sab_all, &
     705              :                       particle_set=particle_set, &
     706              :                       molecule_set=molecule_set, &
     707              :                       matrix_s=matrix_s, &
     708              :                       matrix_ks=matrix_ks, &
     709              :                       mos=mos, &
     710           24 :                       para_env=para_env)
     711              : 
     712           24 :       natom = SIZE(particle_set)
     713           24 :       nsubset = SIZE(molecule_set)
     714           24 :       nspins = dft_control%nspins
     715           24 :       dcdr_env%nspins = dft_control%nspins
     716              : 
     717           24 :       NULLIFY (dcdr_env%matrix_s)
     718              :       CALL build_overlap_matrix(ks_env, matrix_s=dcdr_env%matrix_s, &
     719              :                                 matrix_name="OVERLAP MATRIX", &
     720              :                                 nderivative=1, &
     721              :                                 basis_type_a="ORB", &
     722              :                                 basis_type_b="ORB", &
     723           24 :                                 sab_nl=sab_orb)
     724              : 
     725           24 :       NULLIFY (dcdr_env%matrix_t)
     726              :       CALL build_kinetic_matrix(ks_env, matrix_t=dcdr_env%matrix_t, &
     727              :                                 matrix_name="KINETIC ENERGY MATRIX", &
     728              :                                 basis_type="ORB", &
     729              :                                 sab_nl=sab_orb, nderivative=1, &
     730           24 :                                 eps_filter=dft_control%qs_control%eps_filter_matrix)
     731              : 
     732              :       ! Get inputs
     733           24 :       CALL section_vals_val_get(dcdr_section, "DISTRIBUTED_ORIGIN", l_val=dcdr_env%distributed_origin)
     734           24 :       CALL section_vals_val_get(loc_section, "_SECTION_PARAMETERS_", l_val=dcdr_env%localized_psi0)
     735           24 :       CALL section_vals_val_get(lr_section, "RESTART", l_val=qs_env%linres_control%linres_restart)
     736           24 :       CALL section_vals_val_get(dcdr_section, "Z_MATRIX_METHOD", l_val=dcdr_env%z_matrix_method)
     737              : 
     738           96 :       dcdr_env%ref_point = 0._dp
     739              : 
     740              :       ! List of atoms
     741           24 :       NULLIFY (tmplist)
     742           24 :       isize = 0
     743           24 :       CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", n_rep_val=n_rep)
     744           24 :       IF (n_rep == 0) THEN
     745           72 :          ALLOCATE (dcdr_env%list_of_atoms(natom))
     746           96 :          DO jg = 1, natom
     747           96 :             dcdr_env%list_of_atoms(jg) = jg
     748              :          END DO
     749              :       ELSE
     750            0 :          DO jg = 1, n_rep
     751            0 :             ALLOCATE (dcdr_env%list_of_atoms(isize))
     752            0 :             CALL section_vals_val_get(dcdr_section, "LIST_OF_ATOMS", i_rep_val=jg, i_vals=tmplist)
     753            0 :             CALL reallocate(dcdr_env%list_of_atoms, 1, isize + SIZE(tmplist))
     754            0 :             dcdr_env%list_of_atoms(isize + 1:isize + SIZE(tmplist)) = tmplist
     755            0 :             isize = SIZE(dcdr_env%list_of_atoms)
     756              :          END DO
     757              :       END IF
     758              : 
     759              :       ! Reference point
     760           24 :       IF (dcdr_env%localized_psi0) THEN
     761              :          ! Get the Wannier localized wave functions and centers
     762            4 :          CALL get_loc_setting(dcdr_env, qs_env)
     763              :       ELSE
     764              :          ! Get the reference point from the input
     765           20 :          CALL section_vals_val_get(dcdr_section, "REFERENCE", i_val=reference)
     766           20 :          CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", explicit=explicit)
     767           20 :          IF (explicit) THEN
     768            0 :             CALL section_vals_val_get(dcdr_section, "REFERENCE_POINT", r_vals=ref_point)
     769              :          ELSE
     770           20 :             IF (reference == use_mom_ref_user) &
     771            0 :                CPABORT("User-defined reference point should be given explicitly")
     772              :          END IF
     773              : 
     774              :          CALL get_reference_point(rpoint=dcdr_env%ref_point, qs_env=qs_env, &
     775              :                                   reference=reference, &
     776           20 :                                   ref_point=ref_point)
     777              :       END IF
     778              : 
     779              :       ! Helper matrix structs
     780              :       NULLIFY (dcdr_env%aoao_fm_struct, &
     781           24 :                dcdr_env%momo_fm_struct, &
     782           24 :                dcdr_env%likemos_fm_struct, &
     783           24 :                dcdr_env%homohomo_fm_struct)
     784              :       CALL get_mo_set(mo_set=mos(1), mo_coeff=mo_coeff, &
     785           24 :                       nao=nao, nmo=nmo, homo=homo)
     786              :       CALL cp_fm_struct_create(dcdr_env%aoao_fm_struct, nrow_global=nao, &
     787              :                                ncol_global=nao, para_env=para_env, &
     788           24 :                                context=mo_coeff%matrix_struct%context)
     789              :       CALL cp_fm_struct_create(dcdr_env%homohomo_fm_struct, nrow_global=homo, &
     790              :                                ncol_global=homo, para_env=para_env, &
     791           24 :                                context=mo_coeff%matrix_struct%context)
     792           24 :       dcdr_env%nao = nao
     793              : 
     794           72 :       ALLOCATE (dcdr_env%nmo(nspins))
     795          100 :       ALLOCATE (dcdr_env%momo_fm_struct(nspins))
     796           76 :       ALLOCATE (dcdr_env%likemos_fm_struct(nspins))
     797           52 :       DO ispin = 1, nspins
     798              :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, &
     799           28 :                          nao=nao, nmo=nmo, homo=homo)
     800              :          CALL cp_fm_struct_create(dcdr_env%momo_fm_struct(ispin)%struct, nrow_global=nmo, &
     801              :                                   ncol_global=nmo, para_env=para_env, &
     802           28 :                                   context=mo_coeff%matrix_struct%context)
     803              :          CALL cp_fm_struct_create(dcdr_env%likemos_fm_struct(ispin)%struct, &
     804           28 :                                   template_fmstruct=mo_coeff%matrix_struct)
     805           80 :          dcdr_env%nmo(ispin) = nmo
     806              :       END DO
     807              : 
     808              :       ! Fields of reals
     809           72 :       ALLOCATE (dcdr_env%deltaR(3, natom))
     810           48 :       ALLOCATE (dcdr_env%delta_basis_function(3, natom))
     811           72 :       ALLOCATE (dcdr_env%apt_el_dcdr(3, 3, natom))
     812           48 :       ALLOCATE (dcdr_env%apt_nuc_dcdr(3, 3, natom))
     813           48 :       ALLOCATE (dcdr_env%apt_total_dcdr(3, 3, natom))
     814              : 
     815          960 :       dcdr_env%apt_el_dcdr = 0._dp
     816          960 :       dcdr_env%apt_nuc_dcdr = 0._dp
     817          960 :       dcdr_env%apt_total_dcdr = 0._dp
     818              : 
     819          312 :       dcdr_env%deltaR = 0.0_dp
     820          312 :       dcdr_env%delta_basis_function = 0._dp
     821              : 
     822              :       ! Localization
     823           24 :       IF (dcdr_env%localized_psi0) THEN
     824           16 :          ALLOCATE (dcdr_env%apt_el_dcdr_per_center(3, 3, natom, dcdr_env%nbr_center(1)))
     825           16 :          ALLOCATE (dcdr_env%apt_el_dcdr_per_subset(3, 3, natom, nsubset))
     826           12 :          ALLOCATE (dcdr_env%apt_subset(3, 3, natom, nsubset))
     827          644 :          dcdr_env%apt_el_dcdr_per_center = 0._dp
     828          484 :          dcdr_env%apt_el_dcdr_per_subset = 0._dp
     829          484 :          dcdr_env%apt_subset = 0.0_dp
     830              :       END IF
     831              : 
     832              :       ! Full matrices
     833          100 :       ALLOCATE (dcdr_env%mo_coeff(nspins))
     834           76 :       ALLOCATE (dcdr_env%dCR(nspins))
     835           76 :       ALLOCATE (dcdr_env%dCR_prime(nspins))
     836           76 :       ALLOCATE (dcdr_env%chc(nspins))
     837           76 :       ALLOCATE (dcdr_env%op_dR(nspins))
     838              : 
     839           52 :       DO ispin = 1, nspins
     840           28 :          CALL cp_fm_create(dcdr_env%dCR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
     841           28 :          CALL cp_fm_create(dcdr_env%dCR_prime(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
     842           28 :          CALL cp_fm_create(dcdr_env%mo_coeff(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
     843           28 :          CALL cp_fm_create(dcdr_env%chc(ispin), dcdr_env%momo_fm_struct(ispin)%struct)
     844           28 :          CALL cp_fm_create(dcdr_env%op_dR(ispin), dcdr_env%likemos_fm_struct(ispin)%struct)
     845              : 
     846           28 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff)
     847           52 :          CALL cp_fm_to_fm(mo_coeff, dcdr_env%mo_coeff(ispin))
     848              :       END DO
     849              : 
     850           24 :       IF (dcdr_env%z_matrix_method) THEN
     851           36 :          ALLOCATE (dcdr_env%matrix_m_alpha(3, nspins))
     852           16 :          DO i = 1, 3
     853           34 :             DO ispin = 1, nspins
     854           18 :                CALL cp_fm_create(dcdr_env%matrix_m_alpha(i, ispin), dcdr_env%likemos_fm_struct(1)%struct)
     855           30 :                CALL cp_fm_set_all(dcdr_env%matrix_m_alpha(i, ispin), 0.0_dp)
     856              :             END DO
     857              :          END DO
     858              :       END IF
     859              : 
     860              :       ! DBCSR matrices
     861           24 :       NULLIFY (dcdr_env%hamiltonian1)
     862           24 :       NULLIFY (dcdr_env%moments)
     863           24 :       NULLIFY (dcdr_env%matrix_difdip)
     864           24 :       NULLIFY (dcdr_env%matrix_core_charge_1)
     865           24 :       NULLIFY (dcdr_env%matrix_s1)
     866           24 :       NULLIFY (dcdr_env%matrix_t1)
     867           24 :       NULLIFY (dcdr_env%matrix_apply_op_constant)
     868           24 :       NULLIFY (dcdr_env%matrix_d_vhxc_dR)
     869           24 :       NULLIFY (dcdr_env%matrix_vhxc_perturbed_basis)
     870           24 :       NULLIFY (dcdr_env%matrix_hc)
     871           24 :       NULLIFY (dcdr_env%matrix_ppnl_1)
     872           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%perturbed_dm_correction, nspins)
     873           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%hamiltonian1, nspins)
     874           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%moments, 3)
     875           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_difdip, 3, 3)
     876           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_core_charge_1, 3)
     877           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp, 3)
     878           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_nosym_temp2, 3)
     879           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_s1, 4)
     880           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_t1, 4)
     881           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_apply_op_constant, nspins)
     882           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_d_vhxc_dR, 3, nspins)
     883           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis, nspins, 6)
     884           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_hc, 3)
     885           24 :       CALL dbcsr_allocate_matrix_set(dcdr_env%matrix_ppnl_1, 3)
     886              : 
     887              :       ! temporary no_symmetry matrix:
     888           96 :       DO i = 1, 3
     889           72 :          CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp(i)%matrix)
     890              :          CALL dbcsr_create(dcdr_env%matrix_nosym_temp(i)%matrix, template=matrix_ks(1)%matrix, &
     891           72 :                            matrix_type=dbcsr_type_no_symmetry)
     892           72 :          CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp(i)%matrix, sab_all)
     893           72 :          CALL dbcsr_set(dcdr_env%matrix_nosym_temp(i)%matrix, 0._dp)
     894              : 
     895           72 :          CALL dbcsr_init_p(dcdr_env%matrix_nosym_temp2(i)%matrix)
     896              :          CALL dbcsr_create(dcdr_env%matrix_nosym_temp2(i)%matrix, template=matrix_ks(1)%matrix, &
     897           72 :                            matrix_type=dbcsr_type_no_symmetry)
     898           72 :          CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_nosym_temp2(i)%matrix, sab_all)
     899           96 :          CALL dbcsr_set(dcdr_env%matrix_nosym_temp2(i)%matrix, 0._dp)
     900              :       END DO
     901              : 
     902              :       ! moments carry the result of build_local_moment_matrix
     903           96 :       DO i = 1, 3
     904           72 :          CALL dbcsr_init_p(dcdr_env%moments(i)%matrix)
     905           72 :          CALL dbcsr_copy(dcdr_env%moments(i)%matrix, matrix_ks(1)%matrix, "dcdr_env%moments")
     906           96 :          CALL dbcsr_set(dcdr_env%moments(i)%matrix, 0.0_dp)
     907              :       END DO
     908           24 :       CALL build_local_moment_matrix(qs_env, dcdr_env%moments, 1, ref_point=[0._dp, 0._dp, 0._dp])
     909              : 
     910           96 :       DO i = 1, 3
     911          312 :          DO j = 1, 3
     912          216 :             CALL dbcsr_init_p(dcdr_env%matrix_difdip(i, j)%matrix)
     913          216 :             CALL dbcsr_copy(dcdr_env%matrix_difdip(i, j)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
     914          288 :             CALL dbcsr_set(dcdr_env%matrix_difdip(i, j)%matrix, 0.0_dp)
     915              :          END DO
     916              :       END DO
     917              : 
     918           52 :       DO ispin = 1, nspins
     919           28 :          CALL dbcsr_init_p(dcdr_env%hamiltonian1(ispin)%matrix)
     920           28 :          CALL dbcsr_init_p(dcdr_env%perturbed_dm_correction(ispin)%matrix)
     921           28 :          CALL dbcsr_init_p(dcdr_env%matrix_apply_op_constant(ispin)%matrix)
     922           28 :          CALL dbcsr_copy(dcdr_env%matrix_apply_op_constant(ispin)%matrix, matrix_ks(1)%matrix)
     923           52 :          CALL dbcsr_copy(dcdr_env%perturbed_dm_correction(ispin)%matrix, matrix_ks(1)%matrix)
     924              :       END DO
     925              : 
     926              :       ! overlap/kinetic matrix: s(1)    normal overlap matrix;
     927              :       !                         s(2:4)  derivatives wrt. nuclear coordinates
     928           24 :       CALL dbcsr_init_p(dcdr_env%matrix_s1(1)%matrix)
     929           24 :       CALL dbcsr_init_p(dcdr_env%matrix_t1(1)%matrix)
     930              : 
     931           24 :       CALL dbcsr_copy(dcdr_env%matrix_s1(1)%matrix, matrix_s(1)%matrix)
     932           24 :       CALL dbcsr_copy(dcdr_env%matrix_t1(1)%matrix, dcdr_env%matrix_t(1)%matrix)
     933              : 
     934           96 :       DO i = 2, 4
     935           72 :          CALL dbcsr_init_p(dcdr_env%matrix_s1(i)%matrix)
     936           72 :          CALL dbcsr_copy(dcdr_env%matrix_s1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
     937              : 
     938           72 :          CALL dbcsr_init_p(dcdr_env%matrix_t1(i)%matrix)
     939           96 :          CALL dbcsr_copy(dcdr_env%matrix_t1(i)%matrix, dcdr_env%matrix_nosym_temp(1)%matrix)
     940              :       END DO
     941              : 
     942              :       ! j=1...3: derivative wrt nucleus A, 4...6: wrt nucleus B
     943           52 :       DO ispin = 1, nspins
     944          220 :          DO j = 1, 6
     945          168 :             CALL dbcsr_init_p(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix)
     946          196 :             CALL dbcsr_copy(dcdr_env%matrix_vhxc_perturbed_basis(ispin, j)%matrix, dcdr_env%matrix_s1(1)%matrix)
     947              :          END DO
     948              :       END DO
     949              : 
     950           96 :       DO i = 1, 3
     951           72 :          CALL dbcsr_init_p(dcdr_env%matrix_hc(i)%matrix)
     952              :          CALL dbcsr_create(dcdr_env%matrix_hc(i)%matrix, template=matrix_ks(1)%matrix, &
     953           72 :                            matrix_type=dbcsr_type_symmetric)
     954           72 :          CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_hc(i)%matrix, sab_orb)
     955           96 :          CALL dbcsr_set(dcdr_env%matrix_hc(i)%matrix, 0.0_dp)
     956              :       END DO
     957              : 
     958           96 :       DO i = 1, 3
     959           72 :          CALL dbcsr_init_p(dcdr_env%matrix_ppnl_1(i)%matrix)
     960              :          CALL dbcsr_create(dcdr_env%matrix_ppnl_1(i)%matrix, template=matrix_ks(1)%matrix, &
     961           72 :                            matrix_type=dbcsr_type_symmetric)
     962           72 :          CALL cp_dbcsr_alloc_block_from_nbl(dcdr_env%matrix_ppnl_1(i)%matrix, sab_orb)
     963           96 :          CALL dbcsr_set(dcdr_env%matrix_ppnl_1(i)%matrix, 0.0_dp)
     964              :       END DO
     965              : 
     966           96 :       DO i = 1, 3
     967          156 :          DO ispin = 1, nspins
     968           84 :             CALL dbcsr_init_p(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix)
     969          156 :             CALL dbcsr_copy(dcdr_env%matrix_d_vhxc_dR(i, ispin)%matrix, dcdr_env%matrix_s1(1)%matrix)
     970              :          END DO
     971              : 
     972           72 :          CALL dbcsr_init_p(dcdr_env%matrix_core_charge_1(i)%matrix)
     973           72 :          CALL dbcsr_copy(dcdr_env%matrix_core_charge_1(i)%matrix, dcdr_env%matrix_s1(1)%matrix)
     974           96 :          CALL dbcsr_set(dcdr_env%matrix_core_charge_1(i)%matrix, 0.0_dp)
     975              :       END DO
     976              : 
     977              :       ! CHC
     978           52 :       DO ispin = 1, nspins
     979           28 :          CALL get_mo_set(mo_set=mos(ispin), mo_coeff=mo_coeff, nao=nao, nmo=nmo)
     980           28 :          CALL cp_fm_create(buf, dcdr_env%likemos_fm_struct(ispin)%struct)
     981              : 
     982           28 :          CALL cp_dbcsr_sm_fm_multiply(matrix_ks(ispin)%matrix, mo_coeff, buf, nmo)
     983              :          ! chc = mo * matrix_ks * mo
     984              :          CALL parallel_gemm('T', 'N', nmo, nmo, nao, &
     985              :                             1.0_dp, mo_coeff, buf, &
     986           28 :                             0.0_dp, dcdr_env%chc(ispin))
     987              : 
     988           80 :          CALL cp_fm_release(buf)
     989              :       END DO
     990              : 
     991              :       CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
     992           24 :                                         "PRINT%PROGRAM_RUN_INFO")
     993              : 
     994           24 :       CALL timestop(handle)
     995           72 :    END SUBROUTINE dcdr_env_init
     996              : 
     997              : ! **************************************************************************************************
     998              : !> \brief Deallocate the dcdr environment
     999              : !> \param qs_env ...
    1000              : !> \param dcdr_env ...
    1001              : ! **************************************************************************************************
    1002           24 :    SUBROUTINE dcdr_env_cleanup(qs_env, dcdr_env)
    1003              : 
    1004              :       TYPE(qs_environment_type), POINTER                 :: qs_env
    1005              :       TYPE(dcdr_env_type)                                :: dcdr_env
    1006              : 
    1007              :       INTEGER                                            :: ispin
    1008              :       TYPE(cp_logger_type), POINTER                      :: logger
    1009              :       TYPE(section_vals_type), POINTER                   :: dcdr_section
    1010              : 
    1011              :       ! Destroy the logger
    1012           24 :       logger => cp_get_default_logger()
    1013           24 :       dcdr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES%DCDR")
    1014           24 :       CALL cp_print_key_finished_output(dcdr_env%output_unit, logger, dcdr_section, "PRINT%APT")
    1015              : 
    1016           24 :       DEALLOCATE (dcdr_env%list_of_atoms)
    1017              : 
    1018           24 :       CALL cp_fm_struct_release(dcdr_env%aoao_fm_struct)
    1019           24 :       CALL cp_fm_struct_release(dcdr_env%homohomo_fm_struct)
    1020           52 :       DO ispin = 1, dcdr_env%nspins
    1021           28 :          CALL cp_fm_struct_release(dcdr_env%momo_fm_struct(ispin)%struct)
    1022           52 :          CALL cp_fm_struct_release(dcdr_env%likemos_fm_struct(ispin)%struct)
    1023              :       END DO
    1024           24 :       DEALLOCATE (dcdr_env%momo_fm_struct)
    1025           24 :       DEALLOCATE (dcdr_env%likemos_fm_struct)
    1026              : 
    1027           24 :       DEALLOCATE (dcdr_env%deltar)
    1028           24 :       DEALLOCATE (dcdr_env%delta_basis_function)
    1029              : 
    1030           24 :       IF (dcdr_env%localized_psi0) THEN
    1031              :          ! DEALLOCATE (dcdr_env%psi0_order)
    1032            4 :          DEALLOCATE (dcdr_env%centers_set(1)%array)
    1033            4 :          DEALLOCATE (dcdr_env%center_list(1)%array)
    1034            4 :          DEALLOCATE (dcdr_env%centers_set)
    1035            4 :          DEALLOCATE (dcdr_env%center_list)
    1036            4 :          DEALLOCATE (dcdr_env%apt_subset)
    1037              :       END IF
    1038              : 
    1039           24 :       DEALLOCATE (dcdr_env%apt_el_dcdr)
    1040           24 :       DEALLOCATE (dcdr_env%apt_nuc_dcdr)
    1041           24 :       DEALLOCATE (dcdr_env%apt_total_dcdr)
    1042           24 :       IF (dcdr_env%localized_psi0) THEN
    1043            4 :          DEALLOCATE (dcdr_env%apt_el_dcdr_per_center)
    1044            4 :          DEALLOCATE (dcdr_env%apt_el_dcdr_per_subset)
    1045              :       END IF
    1046              : 
    1047              :       ! Full matrices
    1048           24 :       CALL cp_fm_release(dcdr_env%dCR)
    1049           24 :       CALL cp_fm_release(dcdr_env%dCR_prime)
    1050           24 :       CALL cp_fm_release(dcdr_env%mo_coeff)
    1051           24 :       CALL cp_fm_release(dcdr_env%chc)
    1052           24 :       CALL cp_fm_release(dcdr_env%op_dR)
    1053              : 
    1054           24 :       IF (dcdr_env%z_matrix_method) THEN
    1055            4 :          CALL cp_fm_release(dcdr_env%matrix_m_alpha)
    1056              :       END IF
    1057              : 
    1058              :       ! DBCSR matrices
    1059           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%perturbed_dm_correction)
    1060           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%hamiltonian1)
    1061           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%moments)
    1062           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_difdip)
    1063           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_core_charge_1)
    1064           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp)
    1065           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_nosym_temp2)
    1066           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s)
    1067           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t)
    1068           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_s1)
    1069           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_t1)
    1070           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_apply_op_constant)
    1071           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_d_vhxc_dR)
    1072           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_vhxc_perturbed_basis)
    1073           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_hc)
    1074           24 :       CALL dbcsr_deallocate_matrix_set(dcdr_env%matrix_ppnl_1)
    1075              : 
    1076           24 :    END SUBROUTINE dcdr_env_cleanup
    1077              : 
    1078              : END MODULE qs_dcdr_utils
        

Generated by: LCOV version 2.0-1