LCOV - code coverage report
Current view: top level - src - ex_property_calculation.F (source / functions) Coverage Total Hit
Test: CP2K Regtests (git:42dac4a) Lines: 96.3 % 108 104
Test Date: 2025-07-25 12:55:17 Functions: 100.0 % 1 1

            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 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_api,                    ONLY: dbcsr_add,&
      21              :                                               dbcsr_copy,&
      22              :                                               dbcsr_create,&
      23              :                                               dbcsr_p_type,&
      24              :                                               dbcsr_release,&
      25              :                                               dbcsr_set,&
      26              :                                               dbcsr_type
      27              :    USE cp_dbcsr_contrib,                ONLY: dbcsr_dot
      28              :    USE cp_dbcsr_operations,             ONLY: dbcsr_allocate_matrix_set,&
      29              :                                               dbcsr_deallocate_matrix_set
      30              :    USE cp_log_handling,                 ONLY: cp_get_default_logger,&
      31              :                                               cp_logger_get_default_unit_nr,&
      32              :                                               cp_logger_type
      33              :    USE cp_output_handling,              ONLY: cp_p_file,&
      34              :                                               cp_print_key_finished_output,&
      35              :                                               cp_print_key_should_output,&
      36              :                                               cp_print_key_unit_nr
      37              :    USE cp_result_methods,               ONLY: cp_results_erase,&
      38              :                                               put_results
      39              :    USE cp_result_types,                 ONLY: cp_result_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          552 :    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          552 :       REAL(KIND=dp), ALLOCATABLE, DIMENSION(:)           :: mcharge
      98              :       REAL(KIND=dp), DIMENSION(3)                        :: cdip, pdip, pedip, rcc, rdip, ria, tdip
      99          552 :       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          552 :       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          552 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
     110          552 :       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          552 :       CALL timeset(routineN, handle)
     115              : 
     116          552 :       rlab(1) = "X"
     117          552 :       rlab(2) = "Y"
     118          552 :       rlab(3) = "Z"
     119              : 
     120          552 :       CALL get_qs_env(qs_env=qs_env, dft_control=dft_control)
     121          552 :       tb = (dft_control%qs_control%dftb .OR. dft_control%qs_control%xtb)
     122              : 
     123          552 :       logger => cp_get_default_logger()
     124          552 :       IF (logger%para_env%is_source()) THEN
     125          276 :          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          552 :                                               subsection_name="DFT%PRINT%MOMENTS")
     132              : 
     133          552 :       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          552 :       CALL timestop(handle)
     269              : 
     270         1104 :    END SUBROUTINE ex_properties
     271              : 
     272              : END MODULE ex_property_calculation
        

Generated by: LCOV version 2.0-1