LCOV - code coverage report
Current view: top level - src - ex_property_calculation.F (source / functions) Hit Total Coverage
Test: CP2K Regtests (git:0de0cc2) Lines: 104 108 96.3 %
Date: 2024-03-28 07:31:50 Functions: 1 1 100.0 %

          Line data    Source code
       1             : !--------------------------------------------------------------------------------------------------!
       2             : !   CP2K: A general program to perform molecular dynamics simulations                              !
       3             : !   Copyright 2000-2024 CP2K developers group <https://cp2k.org>                                   !
       4             : !                                                                                                  !
       5             : !   SPDX-License-Identifier: GPL-2.0-or-later                                                      !
       6             : !--------------------------------------------------------------------------------------------------!
       7             : 
       8             : ! **************************************************************************************************
       9             : !> \brief Routines for property calculations of excited states
      10             : !> \par History
      11             : !>       02.2020 Adapted from ec_properties
      12             : !> \author JGH
      13             : ! **************************************************************************************************
      14             : MODULE ex_property_calculation
      15             :    USE atomic_kind_types,               ONLY: atomic_kind_type,&
      16             :                                               get_atomic_kind
      17             :    USE cell_types,                      ONLY: cell_type,&
      18             :                                               pbc
      19             :    USE cp_control_types,                ONLY: dft_control_type
      20             :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      21             :                                               dbcsr_deallocate_matrix_set
      22             :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      23             :                                               cp_logger_get_default_unit_nr,&
      24             :                                               cp_logger_type
      25             :    USE cp_output_handling,              ONLY: cp_p_file,&
      26             :                                               cp_print_key_finished_output,&
      27             :                                               cp_print_key_should_output,&
      28             :                                               cp_print_key_unit_nr
      29             :    USE cp_result_methods,               ONLY: cp_results_erase,&
      30             :                                               put_results
      31             :    USE cp_result_types,                 ONLY: cp_result_type
      32             :    USE dbcsr_api,                       ONLY: dbcsr_add,&
      33             :                                               dbcsr_copy,&
      34             :                                               dbcsr_create,&
      35             :                                               dbcsr_dot,&
      36             :                                               dbcsr_p_type,&
      37             :                                               dbcsr_release,&
      38             :                                               dbcsr_set,&
      39             :                                               dbcsr_type
      40             :    USE distribution_1d_types,           ONLY: distribution_1d_type
      41             :    USE input_section_types,             ONLY: section_get_ival,&
      42             :                                               section_get_lval,&
      43             :                                               section_vals_get_subs_vals,&
      44             :                                               section_vals_type,&
      45             :                                               section_vals_val_get
      46             :    USE kinds,                           ONLY: default_path_length,&
      47             :                                               default_string_length,&
      48             :                                               dp
      49             :    USE message_passing,                 ONLY: mp_para_env_type
      50             :    USE moments_utils,                   ONLY: get_reference_point
      51             :    USE mulliken,                        ONLY: mulliken_charges
      52             :    USE particle_types,                  ONLY: particle_type
      53             :    USE physcon,                         ONLY: debye
      54             :    USE qs_environment_types,            ONLY: get_qs_env,&
      55             :                                               qs_environment_type
      56             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      57             :                                               qs_kind_type
      58             :    USE qs_moments,                      ONLY: build_local_moment_matrix
      59             :    USE qs_p_env_types,                  ONLY: qs_p_env_type
      60             :    USE qs_rho_types,                    ONLY: qs_rho_get,&
      61             :                                               qs_rho_type
      62             : #include "./base/base_uses.f90"
      63             : 
      64             :    IMPLICIT NONE
      65             : 
      66             :    PRIVATE
      67             : 
      68             :    ! Global parameters
      69             : 
      70             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'ex_property_calculation'
      71             : 
      72             :    PUBLIC :: ex_properties
      73             : 
      74             : CONTAINS
      75             : 
      76             : ! **************************************************************************************************
      77             : !> \brief ...
      78             : !> \param qs_env ...
      79             : !> \param matrix_pe ...
      80             : !> \param p_env ...
      81             : ! **************************************************************************************************
      82         550 :    SUBROUTINE ex_properties(qs_env, matrix_pe, p_env)
      83             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      84             :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_pe
      85             :       TYPE(qs_p_env_type)                                :: p_env
      86             : 
      87             :       CHARACTER(LEN=*), PARAMETER                        :: routineN = 'ex_properties'
      88             : 
      89             :       CHARACTER(LEN=8), DIMENSION(3)                     :: rlab
      90             :       CHARACTER(LEN=default_path_length)                 :: filename
      91             :       CHARACTER(LEN=default_string_length)               :: description
      92             :       INTEGER                                            :: akind, handle, i, ia, iatom, idir, &
      93             :                                                             ikind, iounit, ispin, maxmom, natom, &
      94             :                                                             nspins, reference, unit_nr
      95             :       LOGICAL                                            :: magnetic, periodic, tb
      96             :       REAL(KIND=dp)                                      :: charge, dd, q, tmp
      97         550 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: mcharge
      98             :       REAL(KIND=dp), DIMENSION(3)                        :: cdip, pdip, pedip, rcc, rdip, ria, tdip
      99         550 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
     100             :       TYPE(atomic_kind_type), POINTER                    :: atomic_kind
     101             :       TYPE(cell_type), POINTER                           :: cell
     102             :       TYPE(cp_logger_type), POINTER                      :: logger
     103             :       TYPE(cp_result_type), POINTER                      :: results
     104         550 :       TYPE(dbcsr_p_type), DIMENSION(:), POINTER          :: matrix_p, matrix_s, moments
     105             :       TYPE(dbcsr_type), POINTER                          :: matrix_pall
     106             :       TYPE(dft_control_type), POINTER                    :: dft_control
     107             :       TYPE(distribution_1d_type), POINTER                :: local_particles
     108             :       TYPE(mp_para_env_type), POINTER                    :: para_env
     109         550 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     110         550 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
     111             :       TYPE(qs_rho_type), POINTER                         :: rho
     112             :       TYPE(section_vals_type), POINTER                   :: print_key
     113             : 
     114         550 :       CALL timeset(routineN, handle)
     115             : 
     116         550 :       rlab(1) = "X"
     117         550 :       rlab(2) = "Y"
     118         550 :       rlab(3) = "Z"
     119             : 
     120         550 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     121         550 :       tb = (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)
     122             : 
     123         550 :       logger => cp_get_default_logger()
     124         550 :       IF (logger%para_env%is_source()) THEN
     125         275 :          iounit = cp_logger_get_default_unit_nr(logger, local=.TRUE.)
     126             :       ELSE
     127             :          iounit = -1
     128             :       END IF
     129             : 
     130             :       print_key => section_vals_get_subs_vals(section_vals=qs_env%input, &
     131         550 :                                               subsection_name="DFT%PRINT%MOMENTS")
     132             : 
     133         550 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key), cp_p_file)) THEN
     134             : 
     135             :          maxmom = section_get_ival(section_vals=qs_env%input, &
     136         154 :                                    keyword_name="DFT%PRINT%MOMENTS%MAX_MOMENT")
     137             :          periodic = section_get_lval(section_vals=qs_env%input, &
     138         154 :                                      keyword_name="DFT%PRINT%MOMENTS%PERIODIC")
     139             :          reference = section_get_ival(section_vals=qs_env%input, &
     140         154 :                                       keyword_name="DFT%PRINT%MOMENTS%REFERENCE")
     141             :          magnetic = section_get_lval(section_vals=qs_env%input, &
     142         154 :                                      keyword_name="DFT%PRINT%MOMENTS%MAGNETIC")
     143         154 :          NULLIFY (ref_point)
     144         154 :          CALL section_vals_val_get(qs_env%input, "DFT%PRINT%MOMENTS%REF_POINT", r_vals=ref_point)
     145             :          unit_nr = cp_print_key_unit_nr(logger=logger, basis_section=qs_env%input, &
     146             :                                         print_key_path="DFT%PRINT%MOMENTS", extension=".dat", &
     147         154 :                                         middle_name="moments", log_filename=.FALSE.)
     148             : 
     149         154 :          IF (iounit > 0) THEN
     150          77 :             IF (unit_nr /= iounit .AND. unit_nr > 0) THEN
     151           0 :                INQUIRE (UNIT=unit_nr, NAME=filename)
     152             :                WRITE (UNIT=iounit, FMT="(/,T2,A,2(/,T3,A),/)") &
     153           0 :                   "MOMENTS", "The electric/magnetic moments are written to file:", &
     154           0 :                   TRIM(filename)
     155             :             ELSE
     156          77 :                WRITE (UNIT=iounit, FMT="(/,T2,A)") "ELECTRIC/MAGNETIC MOMENTS"
     157             :             END IF
     158             :          END IF
     159             : 
     160         154 :          IF (periodic) THEN
     161           0 :             CPABORT("Periodic moments not implemented with TDDFT")
     162             :          ELSE
     163         154 :             CPASSERT(maxmom < 2)
     164         154 :             CPASSERT(.NOT. magnetic)
     165         154 :             IF (maxmom == 1) THEN
     166         154 :                CALL get_qs_env(qs_env=qs_env, cell=cell, para_env=para_env)
     167             :                ! reference point
     168         154 :                CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
     169             :                ! nuclear contribution
     170         154 :                cdip = 0.0_dp
     171             :                CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
     172         154 :                                qs_kind_set=qs_kind_set, local_particles=local_particles)
     173         484 :                DO ikind = 1, SIZE(local_particles%n_el)
     174         725 :                   DO ia = 1, local_particles%n_el(ikind)
     175         241 :                      iatom = local_particles%list(ikind)%array(ia)
     176             :                      ! fold atomic positions back into unit cell
     177        1928 :                      ria = pbc(particle_set(iatom)%r - rcc, cell) + rcc
     178         964 :                      ria = ria - rcc
     179         241 :                      atomic_kind => particle_set(iatom)%atomic_kind
     180         241 :                      CALL get_atomic_kind(atomic_kind, kind_number=akind)
     181         241 :                      CALL get_qs_kind(qs_kind_set(akind), core_charge=charge)
     182        1294 :                      cdip(1:3) = cdip(1:3) - charge*ria(1:3)
     183             :                   END DO
     184             :                END DO
     185         154 :                CALL para_env%sum(cdip)
     186             :                !
     187             :                ! electronic contribution
     188         154 :                CALL get_qs_env(qs_env=qs_env, rho=rho, matrix_s=matrix_s)
     189         154 :                CALL qs_rho_get(rho, rho_ao=matrix_p)
     190         154 :                nspins = SIZE(matrix_p, 1)
     191         154 :                IF (tb) THEN
     192           4 :                   ALLOCATE (matrix_pall)
     193           4 :                   CALL dbcsr_create(matrix_pall, template=matrix_s(1)%matrix)
     194           4 :                   CALL dbcsr_copy(matrix_pall, matrix_s(1)%matrix, "Moments")
     195           4 :                   CALL dbcsr_set(matrix_pall, 0.0_dp)
     196           8 :                   DO ispin = 1, nspins
     197           4 :                      CALL dbcsr_add(matrix_pall, matrix_p(ispin)%matrix, 1.0_dp, 1.0_dp)
     198           4 :                      CALL dbcsr_add(matrix_pall, matrix_pe(ispin)%matrix, 1.0_dp, 1.0_dp)
     199           8 :                      CALL dbcsr_add(matrix_pall, p_env%p1(ispin)%matrix, 1.0_dp, 1.0_dp)
     200             :                   END DO
     201           4 :                   CALL get_qs_env(qs_env=qs_env, natom=natom)
     202             :                   ! Mulliken charges
     203          12 :                   ALLOCATE (mcharge(natom))
     204             :                   !
     205           4 :                   CALL mulliken_charges(matrix_pall, matrix_s(1)%matrix, para_env, mcharge)
     206             :                   !
     207           4 :                   rdip = 0.0_dp
     208           4 :                   pdip = 0.0_dp
     209           4 :                   pedip = 0.0_dp
     210          20 :                   DO i = 1, SIZE(particle_set)
     211         128 :                      ria = pbc(particle_set(i)%r - rcc, cell) + rcc
     212          64 :                      ria = ria - rcc
     213          16 :                      q = mcharge(i)
     214          68 :                      rdip = rdip + q*ria
     215             :                   END DO
     216           4 :                   CALL dbcsr_release(matrix_pall)
     217           4 :                   DEALLOCATE (matrix_pall)
     218           4 :                   DEALLOCATE (mcharge)
     219             :                ELSE
     220             :                   ! KS-DFT
     221         150 :                   NULLIFY (moments)
     222         150 :                   CALL dbcsr_allocate_matrix_set(moments, 4)
     223         750 :                   DO i = 1, 4
     224         600 :                      ALLOCATE (moments(i)%matrix)
     225         600 :                      CALL dbcsr_copy(moments(i)%matrix, matrix_s(1)%matrix, "Moments")
     226         750 :                      CALL dbcsr_set(moments(i)%matrix, 0.0_dp)
     227             :                   END DO
     228         150 :                   CALL build_local_moment_matrix(qs_env, moments, 1, ref_point=rcc)
     229             :                   !
     230         150 :                   rdip = 0.0_dp
     231         150 :                   pdip = 0.0_dp
     232         150 :                   pedip = 0.0_dp
     233         316 :                   DO ispin = 1, nspins
     234         814 :                      DO idir = 1, 3
     235         498 :                         CALL dbcsr_dot(matrix_pe(ispin)%matrix, moments(idir)%matrix, tmp)
     236         498 :                         pedip(idir) = pedip(idir) + tmp
     237         498 :                         CALL dbcsr_dot(matrix_p(ispin)%matrix, moments(idir)%matrix, tmp)
     238         498 :                         pdip(idir) = pdip(idir) + tmp
     239         498 :                         CALL dbcsr_dot(p_env%p1(ispin)%matrix, moments(idir)%matrix, tmp)
     240         664 :                         rdip(idir) = rdip(idir) + tmp
     241             :                      END DO
     242             :                   END DO
     243         150 :                   CALL dbcsr_deallocate_matrix_set(moments)
     244             :                END IF
     245             : 
     246         616 :                tdip = -(rdip + pedip + pdip + cdip)
     247             : 
     248         154 :                IF (unit_nr > 0) THEN
     249          77 :                   WRITE (unit_nr, "(T3,A)") "Dipoles are based on the traditional operator."
     250         308 :                   dd = SQRT(SUM(tdip(1:3)**2))*debye
     251          77 :                   WRITE (unit_nr, "(T3,A)") "Dipole moment [Debye]"
     252             :                   WRITE (unit_nr, "(T5,3(A,A,F14.8,1X),T60,A,T67,F14.8)") &
     253         308 :                      (TRIM(rlab(i)), "=", tdip(i)*debye, i=1, 3), "Total=", dd
     254         308 :                   WRITE (unit_nr, FMT="(T2,A,T61,E20.12)") ' DIPOLE : CheckSum  =', SUM(ABS(tdip))
     255             :                END IF
     256             :             END IF
     257             :          END IF
     258             : 
     259         154 :          CALL get_qs_env(qs_env=qs_env, results=results)
     260         154 :          description = "[DIPOLE]"
     261         154 :          CALL cp_results_erase(results=results, description=description)
     262         154 :          CALL put_results(results=results, description=description, values=tdip(1:3))
     263             : 
     264             :          CALL cp_print_key_finished_output(unit_nr=unit_nr, logger=logger, &
     265         154 :                                            basis_section=qs_env%input, print_key_path="DFT%PRINT%MOMENTS")
     266             :       END IF
     267             : 
     268         550 :       CALL timestop(handle)
     269             : 
     270        1100 :    END SUBROUTINE ex_properties
     271             : 
     272             : END MODULE ex_property_calculation

Generated by: LCOV version 1.15