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
      10             : ! **************************************************************************************************
      11             : MODULE qs_loc_dipole
      12             :    USE atomic_kind_types,               ONLY: get_atomic_kind
      13             :    USE cell_types,                      ONLY: cell_type,&
      14             :                                               pbc
      15             :    USE cp_control_types,                ONLY: dft_control_type
      16             :    USE cp_log_handling,                 ONLY: cp_logger_type
      17             :    USE cp_output_handling,              ONLY: cp_iter_string,&
      18             :                                               cp_p_file,&
      19             :                                               cp_print_key_finished_output,&
      20             :                                               cp_print_key_should_output,&
      21             :                                               cp_print_key_unit_nr
      22             :    USE cp_result_methods,               ONLY: cp_results_erase,&
      23             :                                               get_results,&
      24             :                                               put_results
      25             :    USE cp_result_types,                 ONLY: cp_result_type
      26             :    USE input_section_types,             ONLY: section_get_ival,&
      27             :                                               section_vals_get_subs_vals,&
      28             :                                               section_vals_type,&
      29             :                                               section_vals_val_get
      30             :    USE kinds,                           ONLY: default_string_length,&
      31             :                                               dp
      32             :    USE mathconstants,                   ONLY: twopi
      33             :    USE moments_utils,                   ONLY: get_reference_point
      34             :    USE particle_types,                  ONLY: particle_type
      35             :    USE physcon,                         ONLY: debye
      36             :    USE qs_environment_types,            ONLY: get_qs_env,&
      37             :                                               qs_environment_type
      38             :    USE qs_kind_types,                   ONLY: get_qs_kind,&
      39             :                                               qs_kind_type
      40             :    USE qs_loc_types,                    ONLY: qs_loc_env_type
      41             : #include "./base/base_uses.f90"
      42             : 
      43             :    IMPLICIT NONE
      44             :    PRIVATE
      45             : 
      46             :    ! Global parameters
      47             :    CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_loc_dipole'
      48             :    PUBLIC :: loc_dipole
      49             : 
      50             : ! **************************************************************************************************
      51             : 
      52             : CONTAINS
      53             : 
      54             : ! **************************************************************************************************
      55             : !> \brief Computes and prints the Dipole (using localized charges)
      56             : !> \param input ...
      57             : !> \param dft_control ...
      58             : !> \param qs_loc_env ...
      59             : !> \param logger ...
      60             : !> \param qs_env the qs_env in which the qs_env lives
      61             : ! **************************************************************************************************
      62          80 :    SUBROUTINE loc_dipole(input, dft_control, qs_loc_env, logger, qs_env)
      63             :       TYPE(section_vals_type), POINTER                   :: input
      64             :       TYPE(dft_control_type), POINTER                    :: dft_control
      65             :       TYPE(qs_loc_env_type), POINTER                     :: qs_loc_env
      66             :       TYPE(cp_logger_type), POINTER                      :: logger
      67             :       TYPE(qs_environment_type), POINTER                 :: qs_env
      68             : 
      69             :       CHARACTER(len=*), PARAMETER                        :: routineN = 'loc_dipole'
      70             : 
      71             :       CHARACTER(LEN=default_string_length)               :: description, descriptionThisDip, iter
      72             :       COMPLEX(KIND=dp)                                   :: zeta
      73             :       COMPLEX(KIND=dp), DIMENSION(3)                     :: ggamma, zphase
      74             :       INTEGER                                            :: handle, i, ikind, ispins, j, n_rep, &
      75             :                                                             reference, unit_nr
      76             :       LOGICAL                                            :: do_berry, first_time, floating, ghost
      77             :       REAL(KIND=dp)                                      :: charge_tot, theta, zeff, zwfc
      78             :       REAL(KIND=dp), DIMENSION(3)                        :: ci, dipole, dipole_old, gvec, rcc, ria
      79          80 :       REAL(KIND=dp), DIMENSION(:), POINTER               :: ref_point
      80             :       TYPE(cell_type), POINTER                           :: cell
      81             :       TYPE(cp_result_type), POINTER                      :: results
      82          80 :       TYPE(particle_type), DIMENSION(:), POINTER         :: particle_set
      83          80 :       TYPE(qs_kind_type), DIMENSION(:), POINTER          :: qs_kind_set
      84             :       TYPE(section_vals_type), POINTER                   :: print_key
      85             : 
      86          80 :       CALL timeset(routineN, handle)
      87             : 
      88          80 :       print_key => section_vals_get_subs_vals(input, "DFT%LOCALIZE%PRINT%TOTAL_DIPOLE")
      89          80 :       IF (BTEST(cp_print_key_should_output(logger%iter_info, print_key, first_time=first_time) &
      90             :                 , cp_p_file)) THEN
      91          16 :          NULLIFY (cell, particle_set, qs_kind_set, ref_point, results)
      92             :          CALL get_qs_env(qs_env=qs_env, &
      93             :                          cell=cell, &
      94             :                          particle_set=particle_set, &
      95             :                          qs_kind_set=qs_kind_set, &
      96          16 :                          results=results)
      97             : 
      98          16 :          reference = section_get_ival(print_key, keyword_name="REFERENCE")
      99          16 :          CALL section_vals_val_get(print_key, "REF_POINT", r_vals=ref_point)
     100          16 :          CALL section_vals_val_get(print_key, "PERIODIC", l_val=do_berry)
     101          16 :          description = '[DIPOLE]'
     102          16 :          descriptionThisDip = '[TOTAL_DIPOLE]'
     103          16 :          CALL get_reference_point(rcc, qs_env=qs_env, reference=reference, ref_point=ref_point)
     104             : 
     105          16 :          dipole = 0.0_dp
     106          16 :          IF (do_berry) THEN
     107          64 :             rcc = pbc(rcc, cell)
     108          16 :             charge_tot = REAL(dft_control%charge, KIND=dp)
     109         256 :             ria = twopi*MATMUL(cell%h_inv, rcc)
     110          64 :             zphase = CMPLX(COS(ria), SIN(ria), KIND=dp)**charge_tot
     111          64 :             ggamma = CMPLX(1.0_dp, 0.0_dp, KIND=dp)
     112             : 
     113             :             ! Nuclear charges
     114         102 :             DO i = 1, SIZE(particle_set)
     115          86 :                CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
     116          86 :                CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost, floating=floating)
     117         188 :                IF (.NOT. ghost .AND. .NOT. floating) THEN
     118          86 :                   CALL get_qs_kind(qs_kind_set(ikind), core_charge=zeff)
     119          86 :                   ria = pbc(particle_set(i)%r, cell)
     120         344 :                   DO j = 1, 3
     121        1032 :                      gvec = twopi*cell%h_inv(j, :)
     122        1032 :                      theta = SUM(ria(:)*gvec(:))
     123         258 :                      zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(zeff)
     124         344 :                      ggamma(j) = ggamma(j)*zeta
     125             :                   END DO
     126             :                END IF
     127             :             END DO
     128             : 
     129             :             ! Charges of the wfc involved
     130             :             ! Warning, this assumes the same occupation for all states
     131          16 :             zwfc = 3.0_dp - REAL(dft_control%nspins, dp)
     132             : 
     133          32 :             DO ispins = 1, dft_control%nspins
     134         148 :                DO i = 1, SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array, 2)
     135         116 :                   ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3, i), cell)
     136         480 :                   DO j = 1, 3
     137        1392 :                      gvec = twopi*cell%h_inv(j, :)
     138        1392 :                      theta = SUM(ria(:)*gvec(:))
     139         348 :                      zeta = CMPLX(COS(theta), SIN(theta), KIND=dp)**(-zwfc)
     140         464 :                      ggamma(j) = ggamma(j)*zeta
     141             :                   END DO
     142             :                END DO
     143             :             END DO
     144          64 :             ggamma = ggamma*zphase
     145          64 :             ci = AIMAG(LOG(ggamma))/twopi
     146         208 :             dipole = MATMUL(cell%hmat, ci)
     147             :          ELSE
     148             :             ! Charges of the atoms involved
     149           0 :             DO i = 1, SIZE(particle_set)
     150           0 :                CALL get_atomic_kind(particle_set(i)%atomic_kind, kind_number=ikind)
     151           0 :                CALL get_qs_kind(qs_kind_set(ikind), ghost=ghost, floating=floating)
     152           0 :                IF (.NOT. ghost .AND. .NOT. floating) THEN
     153           0 :                   CALL get_qs_kind(qs_kind_set(ikind), core_charge=zeff)
     154           0 :                   ria = pbc(particle_set(i)%r, cell)
     155           0 :                   dipole = dipole + zeff*(ria - rcc)
     156             :                END IF
     157             :             END DO
     158             : 
     159             :             ! Charges of the wfc involved
     160             :             ! Warning, this assumes the same occupation for all states
     161           0 :             zwfc = 3.0_dp - REAL(dft_control%nspins, dp)
     162             : 
     163           0 :             DO ispins = 1, dft_control%nspins
     164           0 :                DO i = 1, SIZE(qs_loc_env%localized_wfn_control%centers_set(ispins)%array, 2)
     165           0 :                   ria = pbc(qs_loc_env%localized_wfn_control%centers_set(ispins)%array(1:3, i), cell)
     166           0 :                   dipole = dipole - zwfc*(ria - rcc)
     167             :                END DO
     168             :             END DO
     169             :          END IF
     170             : 
     171             :          ! Print and possibly store results
     172             :          unit_nr = cp_print_key_unit_nr(logger, print_key, extension=".Dipole", &
     173          16 :                                         middle_name="TOTAL_DIPOLE")
     174          16 :          IF (unit_nr > 0) THEN
     175           8 :             IF (first_time) THEN
     176             :                WRITE (unit=unit_nr, fmt="(A,T31,A,T88,A,T136,A)") &
     177           8 :                   "# iter_level", "dipole(x,y,z)[atomic units]", &
     178           8 :                   "dipole(x,y,z)[debye]", &
     179          16 :                   "delta_dipole(x,y,z)[atomic units]"
     180             :             END IF
     181           8 :             iter = cp_iter_string(logger%iter_info)
     182           8 :             CALL get_results(results, descriptionThisDip, n_rep=n_rep)
     183           8 :             IF (n_rep == 0) THEN
     184           3 :                dipole_old = 0._dp
     185             :             ELSE
     186           5 :                CALL get_results(results, descriptionThisDip, dipole_old, nval=n_rep)
     187             :             END IF
     188           8 :             IF (do_berry) THEN
     189             :                WRITE (unit=unit_nr, fmt="(a,9(es18.8))") &
     190          88 :                   iter(1:15), dipole, dipole*debye, pbc(dipole - dipole_old, cell)
     191             :             ELSE
     192             :                WRITE (unit=unit_nr, fmt="(a,9(es18.8))") &
     193           0 :                   iter(1:15), dipole, dipole*debye, (dipole - dipole_old)
     194             :             END IF
     195             :          END IF
     196          16 :          CALL cp_print_key_finished_output(unit_nr, logger, print_key)
     197          16 :          CALL cp_results_erase(results, description)
     198          16 :          CALL put_results(results, description, dipole)
     199          16 :          CALL cp_results_erase(results, descriptionThisDip)
     200          16 :          CALL put_results(results, descriptionThisDip, dipole)
     201             :       END IF
     202             : 
     203          80 :       CALL timestop(handle)
     204             : 
     205          80 :    END SUBROUTINE loc_dipole
     206             : 
     207             : END MODULE qs_loc_dipole